Describing the Structure and Relationship of Height and Diameter in an Old Unmanaged Eucalyptus spp. Plantation

The aim of this study was to model the height and diameter distributions and the height-diameter relationship in an old unmanaged Eucalyptus stand in southeastern Brazil. The diameter at 1.30 m above soil (DBH) and total height (H) of all 469 trees in the stand were measured. For the DBH-H modeling, besides testing various traditional linear and nonlinear models, the inclusion of diameter class as random effect was also testes. For the DBH and H distributions modeling, seven probability density function were tested. The DBH-H relationship was best represented by the Henriksen’s linear model. Including diameter classes as a random effect did not improve estimates. The three parameters Weibull distribution was the best model to describe the diameter and height in old and unmanaged planting of Eucalyptus. The DBH distribution showed the highest frequency of individuals in lower diametrical classes, influenced by competition, while H had distribution similar to normal.


INTRODUCTION
The Brazilian commercial plantations with Eucalyptus genus occupy over 7 million hectares (Ferraz Filho et al., 2018) and stimulates the national economy with a representative part of gross domestic product (IBÁ, 2019). In addition to its role in the national economy, these forests are of significant importance for small and medium-sized producers, who depend on this raw material as a source of income.
Besides the intensively managed plantations of large companies, Eucalyptus stands are also found in National Forests (Azevedo et al., 2011), small and medium-sized rural properties, roadside and experiments. Due to the tax incentives and reforestation policies carried out in the 1970s (Antonangelo & Bacha, 1998), when there was the establishment of the first seminal plantings, Eucalyptus have been used as wood stock for eventual use . These stands are generally of advanced age, unmanaged, from seminal plantings, with diversity of species and variability between individuals, presence of large individuals (diameter and height) and with a large individual volume of wood .
Forest inventory of Eucalyptus plantations, which includes records of diameter at breast height (DBH) and total height of the trees, is essential for planning and management of forest resources (Ferraz Filho et al., 2018), since this information is used in forecasting the volume of wood and production estimate. In addition to the practical/commercial aspect, these measures are important in analysing the distribution and structure of the forest . However, in unmanaged stands, measurements become more difficult and estimates imply greater error . Thus, the forest inventory becomes even more expensive and time consuming, especially in relation to obtaining the total height (Ferraz Filho et al., 2018).
Mathematical models can be used to estimate height as a function of diameter, known as height-diameter (DBH-H) models, and thus reduce the time and cost of sampling. The height-diameter are generally in the form of an ascending curve, with the slope and concavity characteristic of each stand; the curve is dynamic, as over the years it reduces its slope, moves the maximum response (top) to the right (Araújo et al., 2012a). The curve is influenced by the plantation age, quality of site, spacing, sociological position, canopy size, species and genetic material (Machado et al., 2008) and, therefore, tends to have a weak biological relationship. To reduce variability, hypsometric models are usually fitted by species and for a given age of the stand.
Eucalyptus plantations in Brazil normally have a cycle of eight years or less, depending on the forest purpose. Studies of these models are rare in older plantations and without any type of management (Ferraz Filho et al., 2018). Despite the efficiency and versatility of some models, the quality of fit in different species and population conditions must be tested by appropriate statistics and thus, identify precise and exact options for this situation (Azevedo et al., 2011). Furthermore, the use of random effect variables may increase the performance of the hypsometric models, as species , stands or locations, ages, basal area, site index (Melo et al., 2017) and diameter classes (Mendonça et al., 2015). However, few studies used diameter classes as a random factor. In General, they use the diameter classes as a stratification factor. The inclusion of Diameter classes may improve estimates, especially in old-growth stands, which have a wide amplitude of DBH and height . Therefore, allowing the model parameters (one or more) to vary according to the diameter classes, especially in wide classes such as 5 and 10 cm. For example, it may be that the trees in the 10-20 class develop a different DBH-H relationship from the trees in the 40-50 class and then from the 80-90 class because of the dominance-suppression process. As diameter classes' effect is unknown, it is essential to be tested as a random effect variable.
In addition to the height-diameter relationship, the description the structure of these variables is a valuable piece of information for the management of such stands. The diameter and height structure can be analysed using probability density functions (PDF) (Debastiani et al., 2019) and allows for the estimation of the probability of trees occurring within previously defined diametric intervals .
In addition to above mentioned functions, other PDF's such as 3-parameter Lognormal and Generalized Extreme Value (GEV) have been used in hydrological studies with good performance. These PDF's represent data with asymmetric distribution and have the potential to represent data of diameter and height in old and unmanaged stands, which tend to have higher data variability.
Although there are several studies on height-diameter curves in Eucalyptus stands Melo et al., 2017;Ferraz Filho et al., 2018), on the distribution and structure of these stands (Castro et al., 2016;Santos et al., 2019), most are directed to commercial plantations, managed and aged up to 8 years. In contrast, the most suitable PDF's to describe the structure of old unmanaged Eucalyptus stands may be considerably distinct.
In this context, the aim of this study was to model the height-diameter relationship and the structure of these variables in an old unmanaged Eucalyptus stand.

Study location and data analysed
The study area is located in Pedralva city, Minas Gerais State, Brazil, latitude: 22.60º S; longitude: 45.42º W, with an area of 1.1 ha, an average altitude of 1,278 meters, in a hilly area with a slope greater than 20% (strong and undulating relief). The climate of the region is classified as Cwb type, according to Köppen (Alvares et al., 2013); the average precipitation is approximately 1,600 mm; and the average annual temperature is 19º C.
Data were collected in a stand of Eucalyptus spp., which was planted in the 1960s. There were no periodic silvicultural operations throughout stand development and only very light interventions were carried out in few occasions such as the felling of some trees with posterior planting of new seedlings. Therefore, this stand can be categorized as unmanaged.
In 2010, when the trees were between 50 and 60 years old, a census was performed, in which 469 trees of different species of Eucalyptus were measured. The circumference was measured at 1.30 m above the ground, using a tape measure, later converted to DBH (diameter at height 1.30 m above the ground), and the total height (h), using a HAGA altimeter. The mean, standard deviation, kurtosis and skewness for the DBH data were 32.64 cm, 15.41 cm, 1.86 and 1.34 and, for total height, 30.55 m, 9.82 m, -0.53 and 0.24, respectively.

Height-diameter modelling
The data were divided into two subsets: the first subset (236 trees) to fit and the second subset to validate (233 trees) the height-diameter curves. The split into fit and validation data set was performed randomly within diametric classes. For model fitting, the representation of subsample of trees in all diameter and height classes was guaranteed. The Chi-square test of independence was performed between the fit and the validation datasets to test the null hypothesis (H 0 ) that the datasets are independent. In this way, there is no relationship between the fit and validation samples. Knowing the value of one variable does not help to predict the value of the other variable. Linear and nonlinear models (Table 1) were tested to relate height to diameter, whose coefficients were estimated via ordinary least square and Gauss-Newton methods, respectively. Where: DBH = diameter at 1.30 m above the ground (cm); H = height (m); ln = natural logarithm; e = Euler constant; β i = coefficients and ɛ = random error.
The homoscedasticity of the linear and nonlinear model residuals was verified by graphic analysis and the Breusch-Pagan test (Breusch & Pagan, 1979) was performed to the linear models. Residuals normality was verified by the Shapiro-Wilk test for linear and nonlinear models. The level of significance adopted in these tests was 5%. The goodness of fit was assessed based on the following indices: correlation coefficient (r) between the estimates and observed values, Akaike's information criterion (AIC), mean absolute error (MAE), standard error of the estimate (Syx), root of mean square error (RMSE) and PBIAS. Based on the results of these indexes, a ranking was prepared, in which the model with the best performance, in each index, received a score of 1, the second-best performance model received a score of 2, and so on. Thus, the best model was the one with the lowest sum in the ranking of statistical indexes. Finally, the fitted models were validated with the second subset for which the statistics r, AIC, MAE and Syx were calculated.
In addition, to allow for greater flexibility of the heightdiameter models, which could possibly result in better estimates, the inclusion of diameter class as a random effect associated to all parameters was evaluated. The amplitude of the diameter classes tested were 3, 5, and 10 cm. The mixed effect models were fit via restricted maximum likelihood. The significance of the random effects was tested by the likelihood ratio test, with a 5% significance level.
Model fitting was run in R, in which the Stats package was used for linear models (R Core Team, 2020) and the nlme package was used for nonlinear models and the mixed modelling (Pinheiro et al., 2020).

Modelling diameter and total height distribution
Seven probability density functions were fitted to the data of the 469 measured trees distributed classes of amplitude of 5 cm for diameter and 5 m for total height ( Table 2). The observed frequency of individuals per class is shown in Table 3.
The PDF's parameters were estimated via maximum likelihood using the EnvStats package (Millard, 2013) of the R software (R Team Core, 2020). The adequacy of the probabilistic models to the DBH and H series was verified by the Kolmogorov-Smirnov (KS) (Téo et al., 2011;Bueno-Coelho et al., 2017) and chi-square (χ 2 ) (Silva, 2012) tests, at a level of 5% significance. The null hypothesis (H 0 ) of these tests is written as: the DBH and H data comes from the specified distribution.

Hypsometric modelling
All of the estimated coefficients of the height-diameter equations were significant (p < 0.05) by the t test (Table 4). The residuals of the fitted curves followed the assumption of normality (p > 0.05), according to the Shapiro-Wilk test (Table 4) and homoscedasticity (p > 0.05), by the graphical analysis ( Figure 1) and/or Breusch-Pagan test (Table 4). For all models, a higher frequency of residuals was found between two standard deviations, indicating that the equations met the presumption of homoscedasticity. The estimated curves with the observed fitting and validation data partitions are shown in Figure 2. The Breusch-Pagan test is not applicable to nonlinear models (Breusch & Pagan, 1979). β 0 , β 1 and β 2 are parameters of the hypsometric models described in Table 1, p-value is the probability of obtaining a test statistic equal to or more extreme than that observed in a sample, under the null hypothesis (H 0 ), i.e., probability of rejecting H 0 .  In general, the ranking of the goodness-of-fit statistics showed the same trend among the models in the fitting and validation steps, with few changes of position (Table 5). Model 1 (Henriksen, 1950) performed best in both fitting and validation steps, although no large differences in goodnessof-fit was detected among the tested models. The inclusion of diametric classes (amplitudes of 3, 5, and 10 cm) as a random effect did not result in great improvement in the goodness-of-fit statistics and was not significant for neither model. As an example, Table 6 shows the goodness-of-fit statistics and the result of the likelihood ration test for the Henriksen model.

Diameter and height distribution
The diameter and height distribution showed positive asymmetry (g = 1.34 and g = 0.24, for DBH and H, respectively). With regards to kurtosis, the DBH presented a platikurtic distribution (k = 1.86), while the height, a leptokurtic distribution (k = -0.53).
The diameter showed a higher frequency of individuals in the lower classes, while the height showed a more homogeneous distribution around the mean (Figure 3). This difference in the distribution of DBH and H corroborates with the dispersion of the H/DBH ratio (Figure 2).
Because of the difference in conservativeness, there were some divergences between the results of the KS and χ 2 tests. The 3-parameter Weibull function was the only non-significant for both tests, which reinforces its better performance to represent the diameter structure in the studied stand.
For the DBH distribution, all functions with 3 parameters (3-parameter Weibull and GEV) were suitable to by the Kolmogorov-Smirnov (KS) test, while only the 3-parameter Weibull function was suitable by the χ 2 . For the height, all functions were suitable in representing their distribution by the KS test, while by the χ 2 test, only the 3-parameter Weibull was suitable (Table 7).

DISCUSSIONS
The height-diameter ratio (DBH-H) can be very variable, even in homogeneous stands, such as clonal stands. But especially in poorly formed or poorly conducted old stands, there is a weaker relationship between height and diameter resulting in trees with the same diameter presenting contrastingly different heights and vice versa. In the present study, this variation in the DBH-H ratio was found, as these are trees with advanced age, without the proper management to favour their growth.
The linear growth phase has already been passed and, in this case, the simple linear model (model 2) does not represent this relationship mathematically, justifying its inferior results (Table 5, Figure 2). However, it is worth highlighting that linear models have been used to express the relationship between height and DBH, showing satisfactory performance (Azevedo et al., 2011;Liu et al., 2017), even in stands with advanced age, e.g., in Eucalyptus (Azevedo et al., 2011) and Pinus stands (Araújo et al., 2012b).
Despite the parabolic trend, with downward concavity (β 2 < 0), the second-degree model (model 3) presented an intermediate performance in fit and validation steps. This might have occurred due to the weak correlation between height and diameter in unmanaged forests (Scolforo, 2005).
The DBH-H ratio tends to be stabilized in old even-aged plantations, thus, it is best represented by logarithmic or nonlinear models. The Henriksen model (model 1) performed best in both fitting and validation. Although this model has a linear structure, the DBH is in the logarithmic form, which smooths out extreme values and, consequently, reduces the variability of the diameter data. This model also stood out in the height-diameter curves in Pinus taeda L. (Nicoletti et al., 2016) and Eremanthus erythropappus (Araújo et al., 2012a).
The great fit of nonlinear models should also be highlighted. These models represented well the upward relationship, with parameters depicting the slope and concavity characteristic of the height diameter trend, as observed in some studies (Mendonça et al., 2015;Sena et al., 2015;Melo et al., 2017). Nonlinear models are usually pointed out as preferable for representing the height-diameter relationship, as they are parsimonious and have the advantage of interpreting their parameters (Melo et al., 2017).
The insertion on the diameter classes in different amplitudes (3, 5 and 10 cm) as random effects associated with the parameters of the Henriksen model did not result in significant improvements in the accuracy and precision on the estimates. It was expected that, due to the large amplitude and variation in diameter and height, the inclusion of diameter classes, as a random effect, would promote a significant improvement, as found by Mendonça et al. (2015), because it would allow for flexibility in the DBH-H relationship within the diameter classes.
The stand hereby studied developed a very distinguishable canopy vertical structure, in contrast to the single layered industrial Eucalyptus stands commonly managed under 7-year long rotation or less. Since trees at lower phytosociological positions (i.e. lower canopy layers) may grow under considerable competition stress conditions which can affect significantly the DBH-H relationship (Scolforo, 2005), we suspected that the allowing the equation parameters to vary according to diameter class would allow the model to account for this competition effect between trees developing in different phytosociological positions.
Alternatively, variables such as field, age or species, when used as variables with random effect, provided performance gains from the tested hypsometric models (Mendonça et al., 2015;Oliveira et al., 2015;Melo et al., 2017). However, in old unmanaged Eucalyptus stands, normally present in small and medium-sized properties, the use of these variables can be rendered unfeasible, either due to the absence of stands, or lack of resources for identifying species and ages, for example.
Regarding the distribution of DBH and H, the DBH showed a positive skew distribution, i.e., higher frequency of trees in the lower classes, while H presented an approximately normal distribution, with a greater number of trees close to the average (30.55 m). This demonstrates that the lack of proper management negatively affects DBH growth, since it is known that the diameter has a higher growth rate in systems managed with thinning, as this procedure provides an increase in spacing between individuals (Assmann, 1970;Schneider et al., 2015;Pezzutti et al., 2016).
This difference in the DBH and H distribution may have influenced the results of the height-diameter curve fits, as previously mentioned. The behaviour of the distribution of DBH and height in old unmanaged stands is greatly influenced by competition and tends to reflect on the height-diameter relationship.
The best PDF choice was made based on the chi-square test (χ 2 ) to represent the distribution of DBH and H, which is more rigorous than the Kolmogorov-Smirnov (KS) test; KS test is permissive in accepting the null hypothesis of adherence, which has been portrayed in studies in the area of engineering and the environment (Beskow et al., 2015;Abreu et al., 2018). Some studies sought to increase the significance level of KS test to 20%, aiming to increase its rigor in accepting the null hypothesis (Silva et al., 2002). The 3-parameter Weibull function was the only PDF that showed no significance for χ 2 test, for the distribution in DBH and H.
The 3-parameter Weibull location parameter (k) is associated with the minimum diameter and, in several studies, it is disregarded by the premise that the function passes through the origin (Binoti et al., 2010). This parameter showed a value of 9.757 due to occurrence of trees in the diameter class of 5 to 10 cm ( Figure 3a). This shows that this parameter was important for a better fit of 3-parameter Weibull in the diametric distribution of the studied old unmanaged Eucalyptus stand, which showed asymmetry on the left. The 3-parameter Weibull parameter for H was 3.189, indicating individuals in the 5 to 10-meter class, suppressed due to lack of management (Figure 3b).
The Weibull shape parameter (β) is also known as the Weibull slope, because the value of β is equal to the slope of the 3-parameter Weibull probability density function curve. Different values of the β parameter may have marked effects on the behaviour of the 3-parameter Weibull distribution. When β = 1, for example, the 3-parameter Weibull reduces to the exponential distribution, which has reverse "J" shape. If β < 1 or β > 1, the function has a density rate that decreases with the decrease/increase of the DBH and H classes. The β was greater than 1, both for DBH (β = 1.568) and for H (β = 3.036), being greater in magnitude for the distribution of H, justifying a behavior closer to the normal distribution of H.
The scale parameter (α) showed a lower value for DBH compared to H. This difference causes an effect on the abscissa axis, in which a skewed behavior was observed on the left for DBH. In old and unmanaged stands, trees tend not to grow in diameter due to competition, and the growth in height is more evident, before stagnation (Assmann, 1970).

CONCLUSION
The linear model of Henriksen was the most suitable in predicting height as a function of the diameter at breast height in the studied old unmanaged Eucalyptus stand. The inclusion of diameter classes as random effect did not result in significant improvement of the model.
The 3-parameter Weibull distribution was the best model to describe the diameter and height distributions.
The diameter distribution was characterized as leptokurtic and positive skewed, while the total height distribution was mostly platykurtic and positive skewed.