Mathematical modelling to determine the greatest height of trees

This study aimed to analyse the critical height of a column whose weight varies vertically in order to obtain a simple scaling law for a tree where the weight distribution considered. We modelled trees as cantilevers that were fixed to the ground and formulated a self-buckling problem for various weight distributions. A formula for calculating the critical height was derived in a simple form that did not include special functions. We obtained a theoretical clarification of the effect of the weight distribution of heavy columns on the buckling behaviour. A widely applicable scaling law for trees was obtained. We found that an actual tree manages to distribute the weight of its trunk and branches along its vertical extent in a manner that adequately secures its critical height. The method and findings of this study are applicable to a wide range of fields, such as the simplification of complicated buckling problems and the study of tree shape quantification.

( ρ(x) = const. ). However, in this study, the density was treated as a function of height to clarify the relationship between the weight distribution and the critical height of a tree. We considered the following density functions: where ρ 0 is the density at the reference point [kg/m 3 ] and n is a parameter that controls the shape of the density distribution. An example of the distribution of each model is shown in Fig. 1. The coordinate axis was aligned with the neutral axis, with x = 0 at the free end and x = l c at the fixed end.
Model A is a model in which the density at the free end is constant and independent of parameter n ( ρ(0) = ρ 0 ) whereas the density at the fixed end is n times the density at the free end ( ρ(l c ) = nρ 0 ) . Model B in Eq. (3) is the inverted form of Model A, i.e., the density at the fixed end is constant (ρ(l c ) = ρ 0 ) whereas the density at the free end is n times the density at the fixed end ( ρ(0) = nρ 0 ). For the same value of n in Models A and B, both models have the same total weight. Based on the observation that the total weight of the branches rarely exceeds the total weight of the trunk 40,41 , we assumed that the range of n in Models A and B was 0 ≤ n ≤ 3 . Models C and D in Eqs. (4) and (5) are models in which the total weight of the system is independent of the parameter n if the height is invariant. The difference between Models C and D is whether the distribution is linear (2) Model A: ρ(x) = n − 1 l c x + 1 ρ 0 , www.nature.com/scientificreports/ or curved. In these models, the density distribution is constant in the vertical direction, as in Greenhill's equation with n = 0 . Moreover, the weight is concentrated at the top when n = −1 and at the bottom when n = 1 . A summary of these results appears in Table 1. Not only to clarify the effect of weight distribution on the critical height, but also to ensure future expandability, we used the above four density distribution functions as they are geometrically simple and can represent various distribution shapes. The calculation model was a cantilevered elastic column with a radius r l [m] and a length l [m] (see Fig. 2). Greenhill's formula 5 that was verified by McMahon 6 and has been broadly used in various research, it had been derived by modelling trees as non-tapered cylinders. Furthermore, the geometrically simple models are desirable to obtain the wide applicability and extensibility. For the above reasons, we modelled trees as non-taper cylinders. The weight W(x) from the upper end to any point x is given by:  where θ is the deflection angle. Note that the deflection angle θ is not included in the integration of Eq. (6) with respect to x because this is only a volume integration to determine the weight W(x) . If we assume that θ is very small, sinθ can be approximated as θ . Therefore, Eq. (7) can be written as From the elastic curve equation of the beam, the bending moment M(x) at any point x can be obtained as where y is the deflection. Because the deflection angle θ is very small, Eq. (9) can be approximated as: Governing differential equation. From the relationship between the shearing force S(x) of Eq. (7) and the end bending moment M(x) of Eq. (10), the governing differential equation is obtained as Equation (11) is a second-order ordinary differential equation that has an independent variable x and dependent variable θ(x) . We integrated the equations that are obtained by substituting the density functions in Eqs.
(2)-(5) for ρ(x) in Eq. (11) and then applied the following transformation:  where ω is a constant. The governing differential equations were obtained in the following form: If we consider the governing equations in their simplest form, the constant ω for all the models can be defined as Critical height formula. The general solution of the governing equations was obtained by the series solution method using Mathematica. For example, in Model C, the series solution of Eq. (15) is obtained as follows: where c 1 and c 2 are arbitrary constants.
We applied the following boundary conditions of the cantilever to the general solution of the governing differential equations, deriving the following algebraic equation for the critical height: By applying the boundary condition at the free end, we found that c 2 = 0 in Eq. (18). Moreover, by applying the boundary condition in the fixed end, Eq. (18) can be rewritten as By considering the condition for finding a non-trivial solution ( c 1 = 0 ), we obtained the following algebraic equation with respect to ξ: When n = 0 , Eq. (21) corresponds to the equation of von Karman and Biot 11 . The critical height l c can be determined for any value of n by calculating the smallest positive real number ξ c that satisfies Eq. (21) and substituting that value of ξ c into the following equation (which follows from Eqs. (11), (15), and (17)): where C is Greenhill's constant ( C ≈ 1.959 ). Equation (22) is expressed in the form of the product of the coefficient varying with the density distribution and the critical height in the constant-density column introduce by Greenhill. Equation (21) has only two variables, ξ c and n . In other words, ξ c varies only with n . This property is the same for all density models, and the critical height can be obtained using the same method.
Buckling mode. In this section, we describe the calculation of the buckling modes. First, we integrated Eq. (18) with respect to ξ , in order to obtain the deflection y(ξ ): dθ dξ = 0(at x = 0) θ = 0(at x = l c ). www.nature.com/scientificreports/ where c 3 is an arbitrary constant. By considering the boundary condition at the fixed end ( y(ωl c ) = 0 ), c 3 is given by: Next, we substituted Eq. (24) into Eq. (23) and applied the condition ξ = ωl c , The deflection y(x) was obtained as follows: Because Eq. (25) includes the unknown coefficient c 1 , we used the deflection ratio y(x)/y max to eliminate it. The maximum deflection y max appears at x = 0 , and the deflection ratio y/y max to draw the buckling modes is given by Therefore, by substituting l c into Eq. (26), the buckling modes under self-weight could be obtained. We obtained the buckling modes for other density models with the same method.
Example: application to real trees. In this section, we consider the effect of weight distribution on critical height in real trees by using the critical height formula derived by the aforementioned methods. Based on the investigation on the weight distribution of real trees 40,41 , we used the following functions to add the trunk weight to Models C and D, respectively: where ρ B is the density of the branches [kg/m 3 ], and ρ T is the density of the trunk [kg/m 3 ]. An example of the distributions in Eqs. (27) and (28) is shown in Fig. 3. In this study, we expressed the density ratio of the branches and the trunk, ρ B /ρ T , as a weight ratio W R because the density ratio of the branches and trunk is equal to their weight ratio. By using the same method as for the other models, the governing differential equations were obtained as follows: Calculation. Method for solving the critical height equation. Algebraic solving of high-order equations such as Eq. (21) is impossible. Therefore, we used numerical computation to find the smallest positive real-valued solution that satisfies the high-order equation presented in this study. For the numerical method, we adopted the secant method, which does not require differentiation to allow for future extensions.
We rewrote Eq. (18) as follows: By increasing the value of ξ , we explored the interval [ξ 0 , Using these values of ξ 0 and ξ 1 as initial values, we iterated using the following equation: We used the following equation as a convergence criterion: where N is the increment of the expansion order and ξ (N) is the minimum positive real solution in the highorder equation (Eq. (18)) when the expansion order is N . In many studies, the height of a tree is measured in centimetres 42,43 . Therefore, we used the following judgement function to determine the expansion order: We increased N by N until the above condition was satisfied. The smallest expansion order N that satisfied Eq. (35) was the expansion order used in this study.
Relationship between critical height and parameter n. In this section, we clarify the relationship between the critical height l c and the density distribution parameter n . The critical height ratio f (n) was defined as follows: where l c(n) was the critical height in the variable density model, and l cs was the critical height in the model with constant density. If the specification of each model is the same, the critical height is larger than the critical height of the density constant state when f (n) is larger than 1. In contrast, the critical height is smaller than the critical height of the state with constant density when f (n) is smaller than 1. Because l c(n) and l cs are given by Eqs. (22) and (1), respectively, the critical height ratio f (n) was obtained as follows: In order to clarify the relationship between the critical height ratio f (n) and the density distribution parameter n , we obtained the regression curve (discretely calculated) from f (n) . By using this method, the critical height ratio f (n) was expressed as a simple function of n . By applying it to Eq. (22), the relationship between the critical height l c and the weight distribution was clarified. In Sect. 2.5, the critical height ratio f (n, W R ) was defined as follows:  1). The quantity l c(n,W R ) is the critical height when the density distribution parameter is n , and the weight ratio is W R . In other words, Eq. (38) replaces ξ c (n) with ξ c (n, W R ) in Eq. (37). We assumed the following models as regression models; each parameter was determined by regression analysis using the non-linear least-squares method. Furthermore, the significance level α was 0.05: The exponential model and the power-law model in Eqs. (39) and (42) are simple models expressed in monomial form. The linear and polynomial models in Eqs. (40) and (41) are more complicated than the above models, as they are both expressed in a polynomial form. In particular, the polynomial model in Eq. (41) is a quadratic polynomial and is the most complicated among the four models. The power-law model in Eq. (42) cannot be applied to data with n ≤ 0 . Therefore, we use only the data with n > 0 to perform regression analysis in the case of the power-law model.

Results and discussion
In this study, based on the physical characteristic values found in 43  In Model A, the critical height l c slowly decreased with increasing n , whereas the critical height l c showed a sharp decrease with increasing n in Model B. Comparing the two models in more detail, the critical in Model B was higher than that in model A in the interval 0 ≤ n < 1 , whereas the two models showed agreement at n = 1 . In the interval 1 < n ≤ 3 , the critical height in Model A was higher than that in Model B. The reason for these results was that Model A concentrated more weight at the top of the tree in the interval 0 ≤ n < 1 , whereas Model B concentrated more weight at the top in the interval 1 < n ≤ 3 . This result indicated that increasing the weight of the upper portion of the tree significantly reduced the critical height, as compared to the result for an increase in the weight of the lower portion.
In Models C and D, which have the same total weight regardless of n , the critical height l c displayed a curvilinear increase n . Retaining the equal total weights, but with the weight distributed such that the upper and lower parts were lighter and heavier, respectively, the critical height showed a dramatic increase. The difference between Model C (straight line) and Model D (curved line) was negligible when n < 0 , the value for Model C being slightly higher. The difference between the two models became apparent at higher values of n . The values for Model D were higher than those for Model C when n > 0.
For all models, when the density was constant (in agreement with Greenhill's model), we obtained a value of l c = 60.455[m] , which is almost the same as Greenhill's solution. This indicates that when the governing differential equation is solved by the series solution method, the solutions are as accurate as those obtained by solving exactly if the number of terms N ≥ 25 is used. The solutions of von Karman and Biot 11 were derived using N = 6.
Derivation a simple scaling law by regression analysis. The relationship between the critical height ratio f (n) and the density distribution parameter n , and the regression curves obtained by non-linear regression analysis are displayed in Fig. 6. The values of the critical height ratio f (n) are indicated on the vertical axes in the figure, and the values of n are indicated on the horizontal axes. Detailed results of the regression analysis are listed in Table 2. All the parameters of all the regression models ( P 1 ∼ P 9 ) are valid because the p-values are consistently smaller than the significance level α = 0.05. www.nature.com/scientificreports/ In Model A, the critical height ratio (the ratio of the actual critical height l c to the value corresponding to the constant-density case) varies from 1.09 to 0.88 (depending on the parameter n ). Moreover, all regression models, except for the power-law model, are consistent with the theoretical solution. Based on the Akaike Information Criterion (AIC), the polynomial model is the most accurate model. However, for the purpose of obtaining a simple scaling law for trees, the simplicity of the regression model is very important. Therefore, the exponential model is a good regression model because it is superior in terms of both accuracy and simplicity.
In Model B, the critical height ratio varies from 1.60 to 0.73. The polynomial and power-law models are consistent with the theoretical solution. In contrast, the exponential and linear models are not accurate in the neighbourhood of n = 0 and n = 3 . Notably, the power-law model is not valid at n = 0. From the standpoint of AIC, the most accurate model is the power-law model. This model has good accuracy and simplicity.
In Model C, the critical height ratio varies from 1.25 to 0.87. All regression models, except for the power-law model, are consistent with the theoretical solution. The power-law model breaks down when n < 0 , and its error with respect to the theoretical solution is larger than that of the other models. Based on the AIC, the most accurate model is the polynomial model. For the purpose of obtaining a simple scaling law for trees, the exponential model is a good regression model because it is superior in both accuracy and simplicity.
In Model D, the critical height ratio varies from 1.38 to 0.85. Regarding the regression models, all models except for the power-law model are consistent with the theoretical solution. The power-law model cannot account for the critical height when n < 0 , and its error with respect to the theoretical solution is larger than that of the other models. The exponential model shows good accuracy and simplicity.
To derive a simple scaling law such as Greenhill's, we considered the use of simple fractions to express the critical heights. Simple expressions of the regression models considered optimal for each density model in terms of appropriate fractions are shown in Table 3. The limits of application for each formula are also listed in the table.
Application to real trees. The results of the application to a real tree, in which the balance of branch and trunk weights and their distribution on the critical height were investigated, are displayed in Fig. 7. The results for Models E and F are shown in panels (a) and (b), respectively.
In both models, the critical height increased with increasing n . Focusing on W R , it was possible to increase the critical height by making the top lighter, even when the same number of branches and leaves were used. www.nature.com/scientificreports/ Comparing Models E and F, the quadratic distribution was more favourable when the top was lighter, and the linear distribution was more favourable when the top was heavier. Based on King et al. 's study of the weight ratio 40 and Hirata's study of the weight distribution of trees 41 , the value ranges we used for W R and n for real trees were: 0.1 ≤ W R ≤ 0.6 and 0 ≤ n ≤ 1 . We found the following maximum and minimum values of f (n, W R ) in these ranges: The physical meaning of the quantities f max and f min in Eqs. (43) and (44) may be expressed as "the maximum and minimum values of the critical height ratio with respect to the case when all branches and leaves are cut off ". Even if a real tree has the most unfavourable configuration ( n = 0, W R = 0.6 ), the mechanical critical height is reduced by only about 15%.
Niklas reported that the critical height l c obtained by the Greenhill equation has a safety factor of approximately S = 4.0 ( S = l real /l c ) with respect to the actual height of the tree, l real 45 . Combining the results of Eqs. (43) and (44) with the results of Niklas et al., we can estimate that the safety factor for the critical height equation, considering the weight distribution of trees and the weight ratio of the trunk to the branches, lies approximately in the following range: Actual trees have highly diverse weight distributions. However, by performing calculations using a density distribution function that can represent such a wide range of weight distributions and then refining the results using actual measurement data, it is possible to estimate the critical height more accurately than using the critical height equation which ignores existing branch and leaf weights.
Buckling modes under self-weight. The buckling modes under self-weight for all the density models are represented in Fig. 8. The vertical axis in each panel represents the dimensionless coordinate x/l c , and    www.nature.com/scientificreports/ the horizontal axis in each panel represents the dimensionless deflection y/y max . The buckling modes for n = 0, 1, 2 and 3 are shown for models A and B. The buckling modes for n = −1, −0.5, 0, 0.5 and 1 are shown for models C, D, E and F. There were no significant differences in the buckling modes for any of the density models. The results are similar to those of the first-order mode in Euler buckling, where the deflection uniformly increases from the fixed end to the free end. No significant difference in the mode shapes was observed when the density distribution parameters n and the weight ratios W R were changed.

Conclusions
In this study, to clarify the effect of the weight distribution of trees on their critical height, the critical height for self-weight buckling was formulated for cylindrical models with various weight distributions and the critical height equations that include the influence of branch weight was derived for the first time. Regression analysis was performed on the discrete theoretical solutions, and a simple relationship between the weight distribution and the limit height was derived for each density model. Furthermore, using the obtained critical heights, the mode shapes under buckling under self-weight for various weight distributions were obtained. As a result, the following new findings were obtained: Table 3. Simplified formulas for critical height.

Model Formula Applicable limit
A l c = 11 10 e −3n/40 C E γ r 2 Figure 7. Relationship between the critical height and the branch distribution (contained in the weight ratio). www.nature.com/scientificreports/ (1) When a certain number of branches and leaves are distributed in the height direction, the risk of buckling due to self-weight can be reduced if the upper part of the trees is lighter and the lower part is heavier. This method of weight distribution can increase the critical height by a maximum of approximately 1.25 times for a linear distribution and by a maximum of approximately 1.38 times for a curvilinear distribution, compared with a constant weight distribution in the vertical direction. In addition, the total allowable weight is approximately 2 times larger in the case of a linear distribution and approximately 2.6 times larger in the case of a curved distribution. (2) The buckling mode when self-weight buckling occurs is similar to the first buckling mode in the case of general long columns, regardless of the weight distribution in the height direction. In addition, the mode shape hardly changes even if the distribution form changes in the same density model. (3) Based on the weight distributions of real trees and the measurements of branch weights, it can be said that trees distribute their branches and leaves in a very rational way that has little effect on the critical height. Even if we consider the most unfavourable condition of a real tree ( n = 0, W R = 0.6 ), the critical height is only reduced by approximately 15% from the condition with no branches or leaves at all. Based on the findings of Niklas 45 , it can be estimated that the safety factor of the actual trees is approximately 3.4 to 3.9 in relation to the theoretical critical height of the present study, considering the branches. The method in this study can obtain more accurate critical height than previous methods that ignore branch weight because a safety factor of approximately 4 was reported in previous research 45 .
In future work, based on the methods and findings of this study, the shape law and rationality of plant morphology will be further explored from a theoretical point of view. Currently, we plan to investigate the effects of crown weight and position using more complex density function models and to extend the method to hollow cylindrical models such as bamboo and to solid and hollow models with a taper. In the future, we intend to summarise these findings and propose the next generation of material-saving and high-performance structural designs and rational new design methods.