Bertalanffy-Pütter Models for the Growth of Tropical Trees and Stands

The 
Bertalanffy-Putter (BP) five-parameter growth model provides a versatile 
framework for the modeling of growth. Using data from a growth experiment in 
literature about the average size-at-age of 24 species of tropical trees over 
ten years in the same area, we identified their best-fit BP-model parameters. 
While different species had different best-fit exponent-pairs, there was a 
model with a good fit to 21 (87.5%) of the data (“Good fit” means a normalized 
root-mean-squared-error NRMSE 
below 2.5%. This threshold was the 95% quantile of the lognormal distribution 
that was fitted to the NRMSE 
values for the best-fit models for the data). In view of the sigmoidal character of this model 
despite the early stand we discuss whether the setting of the growth experiment may have 
impeded growth.


Introduction
There are multiple factors that influence the growth of trees and forest stands [1] and therefore also a variety of growth equations used in forest science [2] [3] [4] and [5].
Practitioners often use simple three-parameter models, e.g. of Brody [6], Gompertz [7], or Verhulst [8], as these are numerically tractable. The four-parameter Chapman-Richards growth function [9] is popular, too (110,000 hits in Google Scholar). However, the use of many different growth models makes the comparison of the outcomes difficult. Therefore, here we consider the five-parameter Bertalanffy-Pütter (BP) model, as it generalizes these models and provides a more unified approach using growth curves with better fits to the data.  [11] and [12], whence a single sigmoidal growth curve is capable of modeling only one of these phases.
For each species we identify the best-fit parameters for the BP-model, from which we compute additional parameters with a silvicultural relevance, such as the location of the inflection point (maximal yearly growth) and asymptotic height (final expected size of the first growth phase). As to the above-mentioned research question, we speculate that the experimental setup of the growth study, from which we have taken the data, may have impeded the growth of most trees.
For, the analysis of the inflection points indicates an early slowing down of growth (unexpected at this early phase) and we identified a single three-parameter model that had a "good fit" (defined later) to most data.

Data
We used the literature data of Table 1 (see also Figure 5) from Devaranavadgi [10]: An unspecified number of trees from 24 species was planted in 1990 and tree heights were measured annually till 2000 (11 data per species). The table informs about the species T01 -T24 and reports their average-height-at-age data (The source paper provides additional information, such as soil composition. There is no information about the standard deviation of the heights). The data are of additional interest, as often tree size is measured by other parameters [15]. Table 1 also informs about the annual mean temperatures (for 1990, …, 1999) for the larger region around Vijaya Pura. We used Mathematica to retrieve them from the Wolfram Alpha database.
The growth data were obtained during a study at the Regional Agricultural Research Station of Vijaya Pura (district Bijapur, Karnataka, India), located in the Deccan plateau. Figure 1 pinpoints the study site. It has a semi-arid climate with temperatures ranging from 15˚C to 42˚C and average annual rainfall of 594 mm with 39 rainy days. The region suffers from deforestation owing to poor management practices, low fertility of soil, and a harsh climate [16]. The study therefore searched for species with added economic or ecological value that were suitable for re-forestation. Note that some of the species that the study considered may become invasive in more humid areas (e.g. mesquite).

Growth Model
The Bertalanffy-Pütter (BP) model describes tree growth by means of the differential Equation (1) of Pütter [17] for the tree height h(t) at time t. The differential equation can be solved analytically, though in general not by means of elementary functions [18].
Tree growth models for temperate climates often use difference equations, i.e.
In temperate climates this assumption is warranted, as the sharp seasonal changes are apparent from the growth rings. We use the differential equation, as we consider tropical trees. Further, the yearly growth data showed random fluctuations, while the height-at-age data used for the fitting of model (1) were comparably smooth.
The model parameters of Equation (1) have no meaning a priori. They are to be determined from fitting the model to height-at-age data: Four parameters are displayed in the equation, namely the exponent-pair a < b and the scaling constants p and q. An additional parameter is the initial value at age 1, meaning h(1) = c > 0. While in forestry literature also negative exponents were considered (a < 0, b = 1), as thereby a growth model of Schnute [19] would fit into the Open Journal of Modelling and Simulation BP-framework [20], this paper assumes non-negative exponents.
In this form, the BP-model was popularized by Bertalanffy [21] [22] [23] as a model of ontogenetic growth. According to Bertalanffy, the growth of animals, plants and biomass would be governed by certain biophysical principles, which in their most general form would be embodied by Equation (1). Specifically, the growth of trees would be based on an allometric relation between living biomass and photosynthetic area [24], whereby the exponent-pair (a, b) would be related to plant metabolism.
Each exponent-pair (a, b) defines a unique three-parameter model BP(a, b), using the parameters p, q, c. For comparison, Figure 2 plots the exponent-pairs of well-known models and compares them with the exponent-pairs that this paper scanned in an initial search for the optimal model parameters.
These parameters are computed from the parameters of Equation (1) as follows; for t infl numerical equation solving is used. Growth is unbounded, if q = 0 and it is not sigmoidal, if a = 0.
Some authors [28] were concerned that asymptotic length would be unreliable if it exceeded the maximal observed length substantially. We therefore compared h max with typical tree heights.

Goodness of Fit and Method of Calibration
We sought for parameters that minimized SSE, the sum of squared errors for fitting the BP-growth function to height-at-age-data. If h(t) is a solution of Equation (1), using certain exponents a < b and parameters p, q, c, and if h i are the n = 11 average weights, then SSE is defined by Equation (4): The data-fitting exercise for the BP-model is more challenging than for the Richards model, where standard optimization routines may run into difficulties [29]. In recent papers a method of data-fitting was developed for the BP-model [30] [31] [32] and [33]. This method was based on a grid-search, whereby we searched the best-fitting exponent-pairs (a, b) on a grid with step size 0.01 in both directions (Figure 1). For each grid point we identified the best fitting model parameters (p, q, c) that minimized SSE using a custom-made variant of the method of simulated annealing [34]. Simulated annealing alone could be used to optimize for all five parameters (a, b, p, q, c) at once, but often the so computed parameters achieved a suboptimal fit. Note that we optimized also for the initial values c, whence in Table 2 the best-fit values for c slightly differed from the observed initial heights h 1 in Table 1.
For each time series, the best-fit parameters (a min , b min , p min , q min , c min ) achieved the least value of SSE, namely SSE min . Thereby, for the exponents we aimed at an accuracy of 0.01 (defined from the grid), while the other parameters were identified with a higher accuracy.
Our method of obtaining the five best-fit parameters of the BP-model requires for each dataset the consideration of many three-parameter models whereby for all these models the best-fit parameters need to be computed. We utilize the surplus information from this approach and ask, if one three-parameter model would fit for all (or for most) species, if the notion of "fit" was somewhat relaxed; i.e.  "good fit" rather than best fit. We thereby identify, for each species, those exponent-pairs, where the corresponding three-parameter growth model has a "good fit", and then we form the intersection of these 24 (or of fewer) sets. Our definition of a "good fit" is based on the normalized root-mean-squared-error NRMSE, as it allows to compare the fit across different species of different height. NRMSE expresses the root-mean-squared-error (RMSE; 11 is the number of data-points) as a percentage of the maximal observed height (h 11 of the last data-point): Further, to define a "good fit" we use a threshold for NRMSE, based on a statistical analysis of the distribution of the observed NRMSE-values of the best-fit models. We found that a lognormal distribution would fit to these values (Section 3.2) and we therefore use the threshold corresponding to the 95% quantile of this distribution. As explained below, for the present data this threshold is 2.5%.

Materials
We used Mathematica 12.1 (Wolfram Research) for computer algebra, including optimization and statistical analysis, and to access meteorological data (Wolfram repository). We used nonparametric methods, as not all data were normally distributed: Spearman rho for rank correlation and Spearman rank test for nonzero correlation (using the permutation method with 1000 Monte-Carlo simulation steps). Where we used a parametric distribution (lognormal distribution), we first tested the distribution assumption using the Anderson-Darling statistic [35] and the threshold of p = 0.05 (95% confidence).

Best Fit Parameters
The size-at-age data of Table 1 (c.f. Figure 5) inform about 11 annual height measurements. It started in 1990 with average heights h 1 = 11 cm for "bastard teak" T08 to 96 cm for mesquite T21 and it ended in 2000 with h 11 = 2.56 m for black plum T23 to 9.60 m for river tamarind T04. Table 2 informs about the growth data, the parameters of the best-fit models, and the goodness of fit.
The best-fit exponent-pairs (a, b) in Table 2 identify those three-parameter models BP(a, b) that achieved the best fit to the data for T01 to T24. Figure 3 Figure 3. Line a = b (blue) in the region a < 1.05, b < 3 and best-fit exponent-pairs (red), exponent-pairs with NRMSE below 5% for all data (green), and with NRMSE below 2.5% for at least 21 data (blue); plot using Mathematica 12.1.

N. Brunner, M. Kühleitner Open Journal of Modelling and Simulation
plots those 15 exponent-pairs that were in the region a ≤ 1.05, b ≤ 3 (Not in the plot were three exponent-pairs with 1.05 < a ≤ 1.7 and six exponent-pairs with 3 < b ≤ 16.74).
As is evident from this picture, none of these exponent-pairs defined a best-fit model for all data, nor did the exponent-pairs concentrate anywhere. Further, the best-fit exponent-pairs differed clearly from those of the named three-parameter models (i.e. Brody, Gompertz etc.). The distances of these exponent-pairs were smaller to the lines of exponent-pairs defining the four-parameter models, but none of these models was best-fit, either. Table 2 informs also about derived parameters, asymptotic height h max and age t infl and size h infl of the inflection point. Three-parameter models have the disadvantage that by Equation (3) there is a fixed ratio between the height at the inflection point and asymptotic height, e.g. 50% ratio t infl /h max for logistic growth.
By contrast, for the best-fit models either there was no inflection point for one data (lemon-scented eucalyptus T14) or this ratio varied widely between 27.6% for black plum T23 and 76.7% for anjan tree T16.
When compared with the sample data, asymptotic height of the best-fit model was in the range of 94% -130% of the maximal observed height h 11 of Table 1.
In view of the Anderson-Darling test, for the present data the quotient h max /h 11 was lognormally distributed (no refutation of this hypothesis: p = 0.26). Using the maximum-likelihood estimates for its location and shape parameters (0.082 and 0.079, respectively) we concluded that with 99.5% probability the maximal projected height of any tree was below 133% of h 11 . However, except for three species (gum-arabica tree T04, orchid tree T07 and sweet inga T18) the usual heights for trees of the considered species (by Table 1 in most cases 20 m or more) exceeded the asymptotic heights considerably (e.g. by the factor 3 -5 for rain tree T22).

Defining the Threshold for a Good Fit
For the best-fit models NRMSE varied between 0.6% for sweet inga T18 and 3.2% for yellow flame-tree T19 ( Table 2). The good fit was also illustrated by Using the Anderson-Darling test, the NRMSE-values of Table 2 followed a lognormal distribution (no refutation owing to p = 0.93). We used the maximum-likelihood method to estimate its location and shape parameters (−4.38 and 0.42, respectively). The 95% quantile of this lognormal distribution led to the threshold NRMSE = 2.5%. We used it to define "good fit". In this sense, the fits of anjan tree T16 and yellow flame-tree T19 were not good. For instance, based on the lognormal distribution, yellow flame-tree T19 was an outlier with a probability of only 1.3% for NRMSE > 3.2%. Indeed (Figure 4), the growth data for T19 were insofar atypical, as height did not always increase (for average data this can occur if e.g. the largest trees are removed) and as growth finally seemed Open Journal of Modelling and Simulation to accelerate (e.g. more light for the smaller trees, when the largest ones are removed).
There was no exponent-pair (a, b), where a single model BP(a, b) could have a good fit to all 24, to 23, or to 22 of our data. The maximal number of data with a common good-fitting model was 21: Figure 3 (blue region) plots the exponent-pairs corresponding to models with good fits to 21 of our data (whereby for each data different best-fit parameters c, p, q were used). Most of the best-fit exponent-pairs were outside this region (but for their data NRMSE in general was below 2.5%). Further, none of the exponent-pairs of the named three-parameter models (e.g. Brody, Gompertz, etc.) was in the blue region of Figure 3. However, the West-model had a good fit to 19 models and amongst the named three-parameter models this was the best outcome.
For comparison, we also considered the threshold 5%. Using it there were exponent-pairs whose BP-models could be fitted to all data with NRMS < 5%; the green region in Figure

Growth and Temperature
In order to explore, to what extent the growth was affected by environmental factors (there was no information to this end in the source paper), we used the annual mean temperatures from

Discussion and Speculation
The results of this paper invite two opposite threads for further speculation, namely that perhaps the West model is a universal tree growth model or, alternatively, that perhaps something is wrong with the data.
The former speculation has been forwarded repeatedly in different contexts. This question leads to the second speculation mentioned above. One reason for the common growth pattern may be the setting of the growth experiment, as it could have impeded the growth of the trees in about the same way. For, the growth of all trees was sensitive to environmental factors, such as ambient temperature (see above: negative correlation). That the growth may have been impeded is suggested from the fact that for all species, except one. Table 2 displays an inflection point. This means that from this moment on, for most species within t infl = 2 to 6 years, growth began to stagnate. This is atypical for the mean heights that are normally observed for young stands with all trees still in the early growth phase. Further, the typical heights for most species (Table 1) clearly indicate much higher trees than reported in Table 2 as asymptotic heights h max Open Journal of Modelling and Simulation We therefore scrutinized the experimental setting as described in the original data sources [10] [16]. The trees have been planted in 2 m distance from each other on very dry red soil with 30 -35 cm depth. A main concern was the preservation of soil moisture, which was achieved by the method of planting. Thus, there are two obvious reasons, why the growth of the trees may have beendecayed: Moisture preservation may have failed, or the selected species were not sufficiently adapted to drought. We refute these hypotheses, as such observations would have been reported in the source paper as major findings. Another conceivable explanation is competition between trees for sunlight, which might reduce height growth of suppressed trees at later stages, affecting the mean height. However, 2 m distance between trees below 10 m height usually does not affect height growth much. Thus, in an arid region there remains the explanation of competition belowground. Indeed, soil was shallow and 2 m distance at 30 cm depth result in 4 m 2 per tree and about 1 m 3 soil for its roots, which may not be enough for trees that have the potential to grow to heights and crown diameters of 10 -30 m. Of course, this proposed "bonsai effect" needs to be tested, but the source paper does not provide sufficient information.
We conclude that BP-models are a useful tool to analyze the height growth of trees and stands. However, the growth curves depend not only on the species but also on the environmental situation. Thereby, as for the present data, the experimental setting may curb the growth to an extent that the same model may fit to many different species.
Finally, we have added Figure 5 as a graphical summary of the paper, plotting the data in different colors and the best-fitting growth curves.