Height–diameter relationship of plantation-grown juvenile black locust trees is differentiated according to their growth rate, which is positively affected by spacing

Total height and breast-height diameter are among the most important tree characteristics and their relationship has been widely investigated primarily because it allows an indirect estimation of the vertical stand structure (i.e., tree height distribution), since direct measurement would be difficult and time consuming. Simple linear or non-linear relationship usually satisfies the requirement for adequate tree height estimation at stand level (‘local height–diameter model’) and in most local models used to study the height–diameter relationships in Bulgarian forests, linear functions of the variables have been tested (Dimitrov 2003). However, the non-linear equation form is generally preferred in this type of study worldwide (Calama and Montero 2004; Lei et al. 2009). Curtis (1967) generalised that the height–diameter curves should have a positive slope that approaches zero as diameter becomes large and pass through the origin, but not necessarily, if there is no interest in small trees, or if they are not present in the stand. On the contrary, constraining the height–diameter curve to pass through (0, 1.3) is especially important when measurements include very young trees (Curtis 1967). Yuancai and Parresol (2001), on the other hand, advised that the functions used to model height–diameter relationships must increase monotonically, have an upper asymptote and an inflection point, but Paulo et al. (2011) questioned the third requirement, because although diameter and height growth curves have an inflection point, this may not be necessarily so for the relationship between height and diameter. In addition, Annighöfer et al. (2016) pointed out that the upper-asymptote requirement should not be regarded as obligatory in height–diameter modelling of juvenile-age trees.

Although height–diameter relationship has been extensively investigated for many tree species within mature stands (e.g., Huang et al. 1992; López Sánchez et al. 2003; Navroodi et al. 2016; Lebedev and Kuzmichev 2020; Han et al. 2021), only a few studies have been conducted to develop height–diameter models for juvenile trees (Lei et al. 2009; Böhm et al. 2011; Annighöfer et al. 2016). The exploration of the relationship between tree height and diameter of young fast-growing deciduous species, used in plantations for wood and biomass production at temperate climate (poplar hybrids, willows, black locust, paulownia etc.), contains some specific challenges as well. In black poplar hybrids, for example, it would be of interest to determine whether the height–diameter relationship is clone-specific and to what extent this allometric relationship is affected by branch removal during the early stages of plant development. Paulownia, on the other hand, is known to be characterised by a strong tapering of the stems and for other broadleaves at juvenile age (e.g., Robinia pseudoacacia L., Populus alba L.) the stem is often difficult to identify because of active branch growth and lack of dominance (Blujdea et al. 2012). These issues raise the question of whether the diameter at the stem base or that at breast height can be good enough as predictors of height and which one is more appropriate.

While simple height–diameter models can be sufficiently accurate at stand level, expanding the predictions to a wider region would probably lead to biased predictions, as the relationship is highly dependent on the growth conditions and stand characteristics, such as plant age, site quality, stand density and ecoregional specificities (Huang et al. 2000; López Sánchez et al. 2003; Crecente-Campo et al. 2010). Generalised model forms and mixed-effects modelling are usually applied to localise the height–diameter relationship to specific stands (Weiskittel et al. 2011). The generalised model forms consider expansion of the simple height–diameter function to a multiple-predictor relationship by inclusion of relevant stand- and/or site-specific independent variables, such as stand density and/or basal area (Staudhammer and LeMay 2000; Shröder and Álvarez-González 2001; Temesgen and Gadow 2004), stand age (Lenhart 1968) and dominant or average stand height and diameter (Schnute 1981; Pienaar et al. 1990; Cimini and Salvati 2011). Mixed-effects models, on the other hand, specify the relationship by differentiating the model parameters through inclusion of specific parameter components and can be applied to both simple (Trincado et al. 2007) and generalised model forms (Calama and Montero 2004; Castedo-Dorado et al. 2006; Crecente-Campo et al. 2010).

Black locust (Robinia pseudoacacia L.) is a multipurpose forest tree, which can be considered a completely naturalised species for Bulgaria due to more than 200 years of adaptation to the local conditions (Kalmukov 2006). Rédei et al. (2015) qualify black locust as the best tree for short-rotation wood energy production plantations, because of its excellent features, such as vigorous growing potential in juvenile phase, very good coppicing ability, high wood density and dry matter production, favourable combustibility, relatively fast drying, easy harvesting and wood processing. The mass propagation of black locust seedlings in Bulgaria began by the end of 19th century (Naydenov and Dimitrova 2018) and initiated the process of expanding afforestation, especially in the lower forest belt (Panayotov et al. 2006). Black locust currently occupies an area of more than 150 000 hectares, which represents around 4% of the total forest area in the country (Petkova et al. 2017), over an altitudinal range of 0 to 1500 m a.s.l. (Kalmukov 2006; Panayotov et al. 2006). The main objective of this study was to explore the relationship between the height and the breast-height diameter of plantation-grown juvenile black locust trees and to answer the following questions: Can the tree height be adequately predicted at stand level from the breast-height diameter and which is the most appropriate functional form? Can the predictions be expanded to a wider region by mixed-effects modelling and which is the most relevant level for model localisation? Can the resultant most adequate mixed-effects model formulation be adapted for practical application by calibration of the random parameter components with stand-level variables and which of them would be appropriate as a predictor?

Collection of data took place in juvenile even-aged black locust plantations established through planting of one-year-old seedlings (Tab. 1). Four rectangular sample plots were installed along the banks of river Maritsa in South-Eastern Bulgaria and three plots of the same shape, each comprising an open-pollinated progeny of a black locust clone, were established in a 5-year-old experimental plantation in Northern Bulgaria. Two of the parental genotypes of the latest data source were locally selected, while the third genotype was the Hungarian clone ‘Roszin Varga’. Growing space per plant was assessed according to the planting scheme and tree age was determined from two to six destructively sampled trees in situ. Breast-height diameters and total heights of 50% of the trees, randomly selected within the plots, were measured with 0.1 cm and 0.1 m precision, respectively. The total number of height–diameter measurements per plot varied from 48 to 70 depending on the plot size and density (Tab. 1).

Based on visual examination of the scatter plot of the height–diameter data (Figs. 1 and 2) and considering the modelling rationale, defined for this allometric relationship, we selected seven one-predictor functions (Tab. 2) of one- or two-parameters to avoid problems with model convergence, and we tested these functions in modelling the height–diameter relationship of the juvenile back locust trees. The first two models (M1 and M2) are polynomials of first and second degree respectively, of constant intercept equal to 1.3 that assures predicted tree height of 1.3 m when breast-height diameter equals zero. Model M3 is widely known as ‘the allometric equation’ or ‘relative growth equation’ (Huxley 1972), while model M6 can be regarded as its generalisation constrained, like the rest of the tested model forms (M4, M5 and M7), in the same manner and for the same reason as models M1 and M2. Equation M6 has been regarded also as the model by Stage (1975) (Lei et al. 2009). Model M4 was suggested by Wykoff et al. (1982), M5 originated from a simple relationship, known as Schumacher function or Michailoff function (Gadow and Hui 1999) and has been referred also as the model by Loetsch et al. (1973) (Lei et al. 2009) and M7 is based on Michaelis–Menten rational function (Bolker 2008). Two of the models (M4 and M5) have an inflection point and three of the selected functions (M4, M5 and M7) have an upper asymptote. All, but the parabolic relationship (M2), are monotonic and all, but the allometric equation (M3), assure height prediction of 1.3 m for breast-height diameter value of 0. The use of the second-degree polynomial (M2) was justified by the juvenile growth stage of the plants, as suggested by Annighöfer et al. (2016).

Figure 1

Figure 2

To select the function that most adequately describes the height–diameter relationship of the juvenile black locust trees, we followed the methodological procedure established in the study by Stankova and Diéguez-Aranda (2013). The seven simple regression models were fitted by plots employing (linear or non-linear) least squares method. Three test statistics (formulae below Tab. 2) were used to evaluate and compare the model performance: the mean bias estimated as an average of the different plot biases, mean squared error (MSE) calculated per plot and aggregated to a single value for each model and the value of the 90th percentile of the absolute values of the relative biases evaluated for the plots. These three characteristics were chosen for assessment of the compared models since they account for both the magnitude (mean bias, aggregated MSE) and the range (aggregated MSE, 90th percentile of the relative bias) of model errors.

The regression model that showed the lowest estimates of the three statistical criteria, was fitted over the entire data set in nine different mixed-effects model forms. Each of the model parameters separately or both of them together were presented as composed of fixed and random components. In addition, the random components were specified at three different levels: plot within geographic region, plot regardless geographic region and geographic region (Tab. 3). The statistical significance of the fixed coefficients and the variance components of the mixed-effects model formulations was assessed and Likelihood Ratio Test was performed to elect the most adequate mixed-effects model of high predictive power. The best mixed-effects model served also as a basis for development of a generalised deterministic model with two predictor variables by expressing the mixed-effects parameter as a function of a plantation-level variable. The best mixed-effects model and the generalised deterministic relationship were compared through their estimates of model bias, model efficiency (ME) and the average of the absolute values of the relative errors (MARE%) (formulae below Tab. 4). Their goodness-of-fit was also considered by performing tests for model bias and diagnostics of residuals. Unbiasedness of the models was judged according to the t-test for mean error equals zero and simultaneous F-test for slope equal to 1 and zero intercept of the linear regression relating observed and predicted values. Normality of errors was examined by Anderson–Darling analytical test and by inspection of the Quantile–Quantile plot. Homoscedasticity of errors was evaluated by exploring the plots of residuals against the predictors and the predicted values.

In addition, calibration of the random parameter component at each occasion with a height–diameter measurement of one tree, following the procedure described by Trincado et al. (2007), was also attempted. The range (0th, 25th, 50th, 75th and 100th percentiles) and the magnitude (Bias, MARE%) of the errors estimated at each occasion after calibration with one-tree measurement were assessed.

Statistical analyses were carried out using msm, nlstools, nortest, car, nlme and stats packages of R software environment (Jackson 2011; Baty et al. 2015; Gross and Ligges 2015; Fox and Weisberg 2019; Pinheiro et al. 2021; R Core Team 2021).

Models M2 and M5 showed the lowest magnitude of errors, while M4 and M5 revealed the narrowest error range (Tab. 2). The function by Loetsch et al. (1973) (M5) outperformed all other tested relationships in all three analysed statistical criteria (Tab. 2); therefore, it was chosen to represent the height–diameter relationship of the juvenile black locust trees and was fitted in mixed-effects model forms.

Three of the sample plots were situated along Danube River in Northern Bulgaria, while the other four were near the banks of river Maritsa in Southern Bulgaria and for this reason we attempted localisation of the random parameter components according to geographic region (Northern vs. Southern Bulgaria). Although all three regionally localised mixed-effects models (GR.01, GR.0, GR.1 in Tab. 3) were superior to the simple relationship applied to all plots, they were inferior to all other mixed-effects model forms (Tab. 3). When the random parameter components were specified according to both plot and region (GR.PP.01, GR.PP.0, GR.PP.1 in Tab. 3), these model forms did not outperform the respective mixed-effects model forms localised only at plot level (Tab. 3). Consequently, height–diameter relationship localisation at plot level, regardless the geographic region, was most suitable for the investigated juvenile black locust data. Expansion of the rate (slope) parameter b₀ alone with a random component (PP.0 in Tab. 3) produced better goodness-of-fit than the shape parameter b₁ expansion alone (PP.1 in Tab. 3), but both models were surpassed by the mixed-effects model with two expanded parameters PP.01 (Tab. 3). Close examination of the parameter components of the best performing mixed-effects formulation PP.01 showed that the fixed coefficient estimates were significantly different from 0 (b₀ = 8.09, SE(b₀) = 0.66, b₁ = −2.55, SE(b₁) = 0.25, Sign. level P < 0.001). The standard errors of variances of random effects, however, showed that the specific component of parameter b₀ was marginally significant (σ²_u = 2.67, SE(σ²_u) = 1.58, P = 0.09) and that of parameter b₁ was not statistically significant (σ²_v = 0.28, SE(σ²_v) = 0.20, P = 0.17), which indicated that the shape parameter random component does not contribute substantially to model predictability improvement. The second best according to the Likelihood Ratio Test mixed-effects model PP.0 showed statistically significant fixed coefficient estimates and marginally significant random component of b₀ (Tab. 4) and therefore it was preferred as a final mixed-effects model choice.

Mixed-effects model application in practice requires calibration of the random parameter component and therefore, we attempted to express the rate parameter of model PP.0 that varied among the plots, as a function of plantation-level variables, thus converting the model into purely deterministic. The estimated values of the rate parameter were plotted against plantation age and against growing space. While no age-related dependence was observed, increase in the parameter with spacing, with a linear pattern of the relationship, was suggested by the charts (data not shown). This observation was used to formulate and test a two-predictor deterministic function, where the slope coefficient was expressed as a linear function of the growing space. Owing to observed heteroscedasticity of residuals, the regression was refitted by generalised non-linear least squares method, with a variance function that was a power function of tree diameter (Tab. 4). The mixed-effects and the generalised deterministic functions showed adequate regarding model bias and residual diagnostics. However, the mixed-effects model revealed a higher predictive power as evidenced by the lower magnitude of the mean prediction error and the higher model efficiency value (Tab. 4), and illustrated by the prediction trajectories (Figs. 1 and 2).

Another approach to calibrate the random parameter component is by additional tree data and one-tree height–diameter measurement for calibration of the plot-level rate parameter was attempted (Tab. 5). All measured trees in the plots were used for calibration one by one and the final estimates were based on all trials. The results showed that MARE% ranged between 5.423 and 14.612% and bias was between −0.363 and 0.506 m. Tendency to tree height overestimation was noticed for PP1.

Huang et al. (1992) observed that depending on the sample sizes and the species groups, many functions perform well in describing the height–diameter relationship and Lei et al. (2009) concluded that little differences between models existed when breast-height diameter was the only independent variable. We compared seven simple linear and non-linear functions that revealed similar goodness-of-fit statistics, but the model by Loetsch et al. (1973), referred elsewhere as the model by Strub (1974), by Burk and Burkhart (1984) or by Buford (1986), performed best in modelling the height–diameter relationship of juvenile black locust trees. The same function was selected to express the dependence for adult Pinus sylvestris var. mongolica trees 13- to 60-years old (Han et al. 2021) and was shown as an adequate two-parameter function in modelling the height–diameter relationship of birch (Lebedev and Kuzmichev 2020). In a study on yield prediction of young black locust plantations, grown on post-mining area, Böhm et al. (2011) established linear relationship between shoot length and diameter, while Redei et al. (2006) fitted age-specific logarithmic height–diameter functions for 16-and 21-year-old Robinia pseudoacacia L. trees, grown in mixed with Populus alba L. plantations.

Our study showed that the height–diameter relationship for the investigated black locust data can be best specified at plot level, distinguishing the plantations according to their growth rate. The black locust trees from the locally selected clones ‘Karaisen’ and ‘Tsarevets’ (PP6 and PP7) showed substantially faster increase in height than the unimproved plant stock as well as the progeny of the introduced Hungarian clone, as indicated by the positive values of the specific parameter components (Tab. 5) and the regression trajectories (Fig. 1). Superiority in growth and productivity of the progenies of ‘Karaisen’ and ‘Tsarevets’ over the introduced ‘Roszin Varga’ at 5, 7 and 11 years of age was reported also in the studies by Dimitrova (2017) and Dimitrova and Kalmukov (2014, 2019).

A positive relationship between the rate parameter of the model by Loetsch et al. (1973), fitted to our black locust data, and growing space was observed in our study that allowed model expansion to two-predictor relationship. Black locust is a fast growing, light-demanding species, whose height growth peaks around 5 years of age (Redei et al. 2008) and therefore the positive correlation with the available growing area could have been expected. Kalmukov (2006) derived the conclusion that black locust responds positively to increase in the growing space although stand stocking and tree diameter seem to be more affected than tree height. Redei et al. (2013) observed small differences in height and diameter growth of 5-year-old black locust trees at spacing ranging from 0.6 to 1.6 m², while general tendency of height increase with growing area is distinguished in the data reported by Kalmukov (2013) for 6-year-old plants grown on Fluvisols at 9 spacings in the range from 0.5 to 4 m².

Since the practical implementation of the mixed-effects models requires calibration of the random parameter components usually with a subsample of supplementary measurements, we examined this option assuming one additional observation for calibration. Trincado et al. (2007) who studied the height–diameter relationship of loblolly pine using mixed-effects modelling concluded that the most significant improvements in accuracy were observed when one sample tree was used for predicting random parameters. In addition, only a marginal gain in accuracy was observed when two or three sample trees were used for calibration that was mainly because of a reduction in bias. On the contrary, in a study on mixed-effects modelling of radiata pine height–diameter relationship, Castedo-Dorado et al. (2006) concluded that obviously, the greater the number of measurements included in the subsample for calibration, the greater the decrease in residual mean squared error (RMSE) and bias and the increase in the coefficient of determination. They found out that the random selection of three trees per plot reduced the bias and the RMSE of the predictions by almost 50% and more than 10%, respectively. Calama and Montero (2004), who investigated the height–diameter relationship of stone pine also confirmed that when subsample size is increased, prediction error decreases. They found that the largest rate of decrease in sum of squares error is obtained when comparing calibration from two random trees with calibration from four random trees and therefore prediction when height measurements for calibration are taken from four randomly selected trees can be considered a good option.

Beside the sample size for random-component calibration, the size of the trees selected for calibration is also important. Subedi and Sharma (2011) found that calibration with small trees underestimated the predictions of the diameter growth of all trees, while calibration with dominant and co-dominant trees overestimated the diameter growth of small trees. On the other hand, the biases in predictions due to calibration with random and average-sized subsamples were smaller. In studying the height–diameter relationship of Norway spruce Adamec (2015) came to the conclusion that well-functioning calibration was the one that selects trees for calibration in three diameter intervals. Three trees should be measured in each interval (nine trees in total) to achieve the practically acceptable goodness-of-fit of the mixed-effects model. According to the author, calibration based on the measurement of trees only in the interval of mean diameter proved inoperable.

Another way to calibrate the random parameter components of the mixed-effects models is by expressing them as functions of stand-level variables. Stankova and Diéguez-Aranda (2010) attempted calibration of the random components of the rate parameter in a Weibull diameter distribution model as a linear function of quadratic mean diameter, dominant height, plantation age and density. In the present study, we preferred to express the mixed-effects model parameter instead of its random component alone as a function of the growing space. Similarly, Paulo et al. (2011) developed non-linear fixed-effects model for cork oak expanding the regression parameters of a simple one-predictor relationship by linear functions of the reciprocal of stand density and of dominant height. However, when the authors fitted the expanded height–diameter model in a mixed-effects form, its predictive power showed superior to that of the fixed-effects model. The superiority of the mixed-effects model of our study to the purely deterministic two-predictor relationship also derives the importance of the subject-specific variation.

Our study derived a mixed-effects and a two-predictor deterministic models, based on the simple non-linear function by Loetsch et al. (1973), to describe the height–diameter relationship of juvenile plantation-grown black locust trees. The rate parameter of the models served to localise the variation among the stands. A positive correlation between the height increase and the growing space was distinguished through the two-predictor deterministic function. The specific component of the rate parameter in the mixed-effects model form was best localised at stand level differentiating the plantations according to their growth potential.