A Parsimonious Generalised Height-Diameter Model for Scots Pine Plantations in bulgaria: a Pragmatic Approach

Considering the state-of-the-art of forest inventory in Bulgaria, our investigation pursued development of a parsimonious generalised height-diameter model for the Scots pine plantations in the country. A number of 2-, 3- and 4-predictor candidate models were examined and compared based on their goodness-of-fit statistics. Data records obtained in variable-sized sample plots, established throughout the distribution range of the plantations and covering the variety of sites, densities and growth stages were used to fit the models. Two hundred twenty-four plot-level measurements and 3056 tree height-diameter pairs were utilised for parameterization. An independent data set of tree-level measurements and two sets of dominant height-diameter pairs, estimated for differently defined top height tree collectives, were used for model validation. Statistical analyses were carried out using packages nlstools, moments, equivalence, car, nlme, stats and the results were illustrated with ggplot2 and graphics packages of R software environment. A modified form of Gaffrey’s model was selected, which estimates the height of a tree through the breast-height tree diameter, mean stand height and diameter, and accounts for the tree social status. It was fitted by generalised non-linear least squares method, with residual variance weighted by a product of tree diameter and mean stand height exponential functions. An adjusted coefficient of determination of 0.917 and residual standard error of 0.794 m indicated the high predictive potential of the derived model. Validation tests showed that the estimated regression line is very well fitted to the independent data and is appropriate to forecast dominant stand heights. The range of errors, relative to the predicted dominant height values, was narrow, ±25-30%, with low magnitude of the average of their absolute values (4-5%). The equivalence tests rejected the null hypothesis of dissimilarity regarding model bias (observations-predictions line intercept) for all validation data sets, for a region of equivalence as narrow as ±5%. The 3-predictor generalised height-diameter model developed in our study needs information readily available from the inventories and therefore can be broadly used. Its application in dominant stand height prediction is recommended.

trees growing at higher densities will probably have smaller diameters than those in less dense stands due to competition (i.e. density-dependence), trees having the same diameter at different times will belong to sociologically different classes (i.e. age-dependence) and the height curves for good quality sites are expected to have steeper slopes than those for poor quality sites (i.e. site quality-dependence) . Analysis of the ecoregion-based height-diameter models for white spruce in Alberta's boreal forests has suggested that the height-diameter relationships are different in different ecoregions, probably because of the very different biogeoclimatic conditions (Huang et al. 2000). The authors found that applying a height-diameter model fitted from one ecoregion to different ecoregions resulted in overestimations between 1.10% and 29.05%, or underestimations between 1.92% and 21.92%. Generalized model forms and mixedeffects modelling are usually applied to localize the heightdiameter relationship to specific stands (Weiskittel et al. 2011). Referring to earlier investigations, Crecente-Campo et al. (2010) generalised that the incorporation of stand variables in a local height-diameter model reduces bias and increases precision. Ahmadi and Alavi (2016) concluded that the inclusion of stand characteristics improved the prediction accuracy of tree height estimation for Fagus orientalis Lipsky trees, while Staudhammer and LeMay (2000) found that the introduction of stand density variables resulted in increased accuracy for predicting heights of alder. Temesgen and Gadow (2004) reported that the expansion of the simple heightdiameter relationships with stand-level attributes reduced the root mean square errors (RMSE) by 30.0 cm. Temesgen et al. (2007), who developed regional height-diameter equations for major tree species of southwest Oregon, estimated decrease in the RMSE of the expanded models by a minimum of 6.1% to a maximum of 22%.
By definition, the mean stand height corresponds to the quadratic mean diameter (or mean basal area diameter), due to its practical importance in volume calculation, and therefore its estimation is straightforward. The definition and, consequently, the estimation of the dominant stand height, however, may vary widely. While Pretzsch (2009) limits the definition for dominant height to 3 versions, according to the top height tree collective defined, Van Laar and Akça (2007) mention 6 different measures of stand dominant height. Poryazov (2009) registered 19 different ways to define and calculate dominant stand height. In Bulgaria, similar to other countries such as Estonia (Tarmu et al. 2020), mean stand height is the measure commonly used in the forest management practice and there is no standard or officially accepted protocol for dominant height estimation. Duhovnikov (1972) determined dominant stand height as the height corresponding to the average diameter of the 20% thickest trees in the stand, while Shikov (1974) and Ferezliev and Tsakov (2010) calculated dominant height as an arithmetic average of the heights of the 100 thickest trees per hectare. Petrin (1987) and Tonchev (2022), on the other hand, applied the definition of dominant height as the one corresponding to the quadratic mean diameter of the 100 thickest trees per hectare. Stankova et al. (2006) derived a relationship of the dominant on mean stand height of Scots pine plantations, using sub-samples of total height measurements, taken on 20% of the trees in the plots that represent the average and dominant diameter trees in equal proportions (i.e. dominant height corresponding to the average height of the 10% thickest trees in the plot). It becomes clear that, for situations where a standard procedure for dominant height estimation is not established, such as in Bulgaria, a height-diameter relationship based on dominant stand height as a predictor would be of limited applicability. A generalised height-diameter model, based on tree diameter, dominant stand diameter and height, stand age and density has been developed for Scots pine plantations (Stankova and Diéguez-Aranda 2013). Its confident application, however, assumes sufficient precision of the stand age estimate and requires additional field-collected data on stand density, as advised by Stankova (2012). In addition, the aforementioned lack of coherence in dominant height estimation in Bulgaria casts doubt on both the consistency of the data used for parameterization and the reliability of the model predictions for independent data sets.
The national forest inventories in Bulgaria are characterised as stand-wise inventories for local forest management planning (FMP) (Groen et al. 2013), which are implemented by management units at approximately 10year intervals. The inventory description of each forest stand includes several principal quantitative parameters: stand age (years), basal area (m 2 ·ha -1 ), standing stock (m 3 ·ha -1 ), quadratic mean diameter (cm), mean stand height (m), stocking rate, relative to "normal" (fully stocked) stand (0-1) and site index expressed as a categorical variable in the range I-V and based on functional relationship of mean height to age. Stand density, expressed by the number of trees per hectare, is rarely recorded and the highest precision for stand age is usually 5 years. According to the Ordinance for inventory and planning of the forest territories of Bulgaria, in the vast majority of cases (stand types by age classes) estimates of the stand basal area and volume are obtained from the growth and yield tables according to the in situ evaluations of mean height, age and stocking rate. The specifics of forest inventory in Bulgaria suggest that the stand-level parameters that are of sufficient precision and are readily available for implementation into generalized height-diameter functions are the average stand height and diameter, and perhaps basal area.
In forest inventories, the mean stand height is needed to estimate the stand volume, while dominant stand height is considered a more suitable measure for predicting site quality, because it is less sensitive to thinnings (Van Laar andAkça 2007, Tarmu et al. 2020). Consequently, the accuracy of the dominant stand height predictions by a height-diameter model is important. Sharma et al. (2002), who investigated top height as defined in 7 different ways, found out that for both thinned and unthinned loblolly pine plantations, there were significant differences in the respective top height estimates. Moreover, with the exception of a few cases at certain measurements, the site index predictions based on them also differed significantly. Therefore, it is important to assess also the height-diameter model accuracy in predicting stand dominant height defined in various ways, i.e. from differently assembled largest trees collectives.
The main objective of our study was to derive a practically-oriented, parsimonious generalized heightdiameter model for Scots pine plantations in Bulgaria and to examine and validate its applicability for prediction of tree height, particularly of the largest trees that form the collective for the estimation of dominant stand height.

Data Sets
The data set used to derive the height-diameter relationship of the trees in the Scots pine plantations in Bulgaria was generated from both personally collected and published data records. These records were obtained in variable-sized sample plots (Table 1) of circular or rectangular form, which were established randomly throughout the distribution range of the plantations, with a primary criterion to cover the variety of sites, densities and growth stages of the Scots pine plantations in Bulgaria. One hundred and twenty-one of the plot-level measurements were collected in once-measured plots, 111 of which were personally established and 10 were inventoried by other researchers (Efremov 2006). The remaining 103 plot-level data records were provided from the Appendices to the Forest Inventory Plans and other published data sources (Marinov 2008) and comprised data from permanent sample plots measured 1 to 4 times. Beside the total height (h, m) -breast-height diameter (d, cm) measurements of trees, sampled in random or systematic manner in each plot, the stand-level parameters: basal area (G, m 2 ·ha -1 ), stand density (N, trees·ha -1 ), quadratic mean diameter (D m , cm) and mean stand height (H m , m) were used to fit the regression models.  In addition to the parameterization data set, we employed 3 data sets for model validation (Table 1). The first of them, Validation Data Set 1, was composed of completely independent data, personally collected in 46 plots measured once, following the same principle of plot distribution as of the parameterization data set. The subsamples of fewer tree heights measured per plot (Table  1) were not appropriate to estimate dominant heights for differently defined largest tree collectives and therefore the applicability of the model was only examined with the 329 pairs of tree-level height-diameter measurements recorded (abbreviated as "Valid.Data" in the results). The second and the third sets for validation were created from subsets of the parameterization data set, which allowed dominant height estimation from differently assembled largest tree collectives. Two types of largest tree groups were compiled. The tree-level data of Validation Data Set 2 was used to calculate the arithmetic averages of the tree diameters and heights of the 10% thickest trees in the plot (as estimated in Stankova et al. 2006), referred hereafter as "top diameter" and "top height" and abbreviated as D d and H d . The tree-level data of Validation Data Set 3 were used to estimate the quadratic mean diameter of the 100 thickest trees per hectare, recalculated according to the plot size, and its respective basal area-weighted (Lorey's) height, referred hereafter as "dominant diameter" and "dominant height" and abbreviated as D 0 and H 0 . The ability of the height-diameter relationship to predict the top/dominant stand height from the top/dominant stand diameter was tested with 111 top height-top diameter pairs (abbreviated as "Dom.10perc" in the results) and 100 dominant heightdominant diameter pairs (abbreviated as "Dom.D0H0" in the results), which were able to be extracted from, respectively, Validation Data Set 2 and Validation Data Set 3 (Table 1).

Model Development, Estimation and Validation
Considering the state-of-the-art of forest inventory in Bulgaria, we concluded that the stand-level parameters which are the most appropriate to be used as predictors in the generalized height-diameter functions are the average stand height and diameter, because they are measured during the inventories and are available in the inventory descriptions of the stands. Generalised height-diameter models that include mean stand diameter and height (or dominant stand diameter and height) are usually derived from simple height-diameter relationships by formulating a second equation of the same functional form, but constrained to pass through the mean (or dominant) stand height and mean (or dominant) stand diameter. The system of the two equations is then solved for tree height, converting the simple one-predictor relationship into a 3-predictor regression model (Mønness 1982, Krumland and Wensel 1988, Cimini and Salvati 2011. Alternatively, a generalized height-diameter model can be derived by expressing the parameters of a simple relationship as functional forms of the selected stand-level variables (Harrison et al. 1986, Castedo-Dorado et al. 2001. A number of such 3-predictor height-diameter models, which are based on tree diameter, mean (or dominant) stand height and mean (or dominant) stand diameter and their modifications, developed from simple predecessor functions that have shown appropriate for modelling the height-diameter relationship, can be found in the literature (e.g., López , Lei et al. 2009). We classified the functions that we extracted into 3 groups: 1) based on D m and H m (Table A1); 2) based on D 0 and H 0 (Table A2) and 3) based on D m and H 0 (Table A3). Two of the third-group relationships, by Harrison et al. (1986) and by Gadow and Hui (1999), did not include a stand-level diameter measure as an independent variable, i.e. they described 2-predictor functions. We examined the models from the first group (7 models) in their original form. The second group formulations (18 models) were tested keeping their functional forms, but replacing D 0 and H 0 by D m and H m , considering that the second group models were derived in an analogous manner to the first group relationships. We fitted the relationships from the third group (7 models) in their base form after substituting H m for H 0, but we also examined their expanded forms where H 0 was substituted by a function of H m . We assumed 2 different forms of the relationship: linear function of H m (H 0 =a 0 H m +a 1 ), as it has been often approximated for practical application (M.L.W. F. 1980, Mihov 1986) and nonlinear function of H m and N (ln(H 0 /H m )=a 0 [ln(N/150)] a 1 ), as suggested by García (2017). Given the importance of basal area as a measure of stand stocking, we considered the possibility for its incorporation into the generalised heightdiameter model as well. Therefore, we examined eight more relationships (Table A4) that included in their originally published form tree diameter, measures of stand stocking (basal area and/or stand density) and, in some cases, standlevel measures of height and diameter (mean or dominant).
We first fitted the formulated generalised heightdiameter relationships by non-linear least squares (NLS) method, we checked if the assumptions to the residuals were met and we examined the goodness-of-fit of the regression models. Homoscedasticity of errors was evaluated by exploring the plots of residuals against the independent variable values and against the predicted values. In case of diagnosed heteroscedasticity of errors, that is often the case for biological data, the model was not rejected, but the option to resolve the issue through refitting by generalised non-linear least squares (GNLS) method was assumed. The data from the Scots pine plantations in Bulgaria have generally revealed a tendency to leptokurtic residual distributions for most of the elaborated biometric models (e.g., Stankova and Diéguez-Aranda 2013, 2017. The relative robustness of the regression analyses against even considerable departures from normality has been considered, provided that the sample size is sufficiently large (Frost, 2014) or the residual distribution is not severely asymmetric (Sokal and Rohlf 1995). Therefore, we examined the normality of errors by separate tests for skewness and kurtosis and we discarded those tested models that showed significantly skewed residual distributions as assessed by D'Agostino's skewness test (D'Agostino 1970). In addition, we visually inspected the histograms of the residual distributions as well as the quantile-quantile plots.
Thereafter, the models were tested for bias by a t-test for mean error equal to zero and simultaneous F-test for slope equal to 1 and zero intercept of the linear regression relating the observed and predicted values. Stability of parameter estimates was assessed by the Percent Relative Standard Error (PRSE%), which is the ratio (in percent) between the standard error and the absolute value of the fitted our Scots pine data and did not include a stand-level diameter measure. It was parameterized successfully after substitution of dominant stand height by mean height, while its expanded model forms failed to converge. The model by Gaffrey (1988) as well as its modification by Diéguez-Aranda et al. (2006) (Gafr2 and Gafr1 in Table 2, respectively) were fitted successfully when the dominant stand height was substituted by a linear function of mean height. Two expanded forms of the model by Pienaar et al. (1990) also showed adequate. The first of them was a function with three independent variables, where dominant height was presented by a linear relationship of the mean height (Pien1 in Table 2) and the second function had four predictors, where a product of mean height with density-based term was considered instead: H m e [ln(N/150)]a1 (Pien2 in Table 2).
The model by Shröder-Álvarez-González (2001) and its predecessor, the model by Mirkovich (1958), supplied four more adequate regressions: derived by a simple substitution of dominant with mean stand height (Mir1, ShrA1 in Table 2) or by incorporation of a product of mean stand height with power term of density -H m (N/150) a 0 -on its place (Mir2, ShrA2 in Table 2). Finally, a model by Sharma and Parton (2007) including the ratio N/G (Table A4, ShaPar2 in Table  2) as well as its modified form where this ratio was replaced by the quadratic mean diameter (Table A2, ShaPar1 in Table  2) were also elected as possible candidates to describe the height-diameter relationship of the trees in the Scots pine plantations in Bulgaria. All adequate generalized models showed high coefficients of determination (above 0.9 in most of the cases) and root mean squared errors in a narrow range: 1.29 -1.47m (Table 2). However, the model Gafr2, which was based on the function by Gaffrey (1988), proved superior to the other compared models, as indicated by Akaike weights (Table 2).
Unbiased, symmetric residual distribution was revealed by both graphical and analytical tests for this model, but heteroscedasticity of residuals was diagnosed ( Figure 1). To assure higher precision of the parameter estimates as well as their standard errors, generalized non-linear least squares method was applied. A variance function that was the product of exponential relationships to tree diameter and mean stand height was selected and applied (Table 3). The likelihood ratio test indicated an increased predictive power of the model fitted via generalized non-linear least squares (Table 3) and the plots showed an improvement in the residual variance and model adequacy (Figures 1 and 2). Parameter Relative Standard Errors (PRSE%) attained values below 20%, assuring the stability of the parameter estimates (Table 3).
High rates of model efficiency were calculated for all examined validation data sets (Table 4). They showed that the estimated regression line is very well fitted to the independent data and is appropriate to forecast top and dominant stand heights. The range of errors, relative to the predicted dominant and top height values, was narrow, ± 25 -30% (Table 4), with low magnitude of the average of their absolute values (4 -5%). A tendency to underestimate the tree heights at the upper size range and overestimate those at the lower size range for the independent validation data can be seen in Table 4 and Figure 3A. It was asserted also by the equivalence test, which showed that the hypothesis for regression parameter and must attain values below 25%. The PRSE% served indirectly also to control the effect of possible outliers and influential observations (Sileshi 2014). Collinearity was handled through the condition number that must obtain estimates below 30.
All models that proved adequate were compared by their regression statistics adjusted coefficient of determination (Adj. R 2 ), residual standard error (RMSE, m) and Akaike Information Criterion (AIC). To select the relationship best suited to the data of our study, we employed Akaike weights (AIC weights, %) that represent the relative likelihoods of the models, providing strength of evidence in favour of one model over the other (Wagenmakers andFarrell 2004, Sileshi 2014). In case of diagnosed heteroscedasticity of errors, the selected model was refitted by the generalised nonlinear least squares method. Variance-correcting functions, based on combinations of different functional forms of the predictors, were tested and compared in consideration of both residual homoscedasticity and overall goodness-of-fit of the regression model. The improvement of the residual variance and the model adequacy by GNLS as compared to NLS was assessed by the residual plots, and the model predictive abilities, as fitted by the alternative methods (NLS vs. GNLS), were contrasted by a likelihood ratio test.
To validate the newly developed model, we fitted it to each validation data set separately and estimated several test statistics (formulae shown under Table 4). Model efficiency (ME), which is a relative measure of model performance similar to the coefficient of determination, was assessed. The quartiles of the relative error distributions as well as the averages of the absolutes values of these errors were evaluated, which indicate the range and the magnitude of prediction errors relative to the predicted heights. Finally, to compare the observed with the predicted values, we performed equivalence tests (Robinson et al. 2005, Weiskittel et al. 2011) which allow goodness-of-fit judgment by combining test of model bias, assessing equality of means (test for the intercept) and test of model proportionality, evaluating similarity of individual observations (test for the slope). The regions of equivalence were set as narrow as ±5% to be compared with the approximate joint two onesided 95% confidence intervals for the slope and intercept for 100 bootstrap replicates.
Statistical analyses were carried out using packages nlstools, moments, equivalence, car, nlme, stats and the results were illustrated with ggplot2 and graphics packages of R software environment (Baty et al. 2015, Komsta and Novomestky 2015, Robinson 2016, Wickham 2016, Fox and Weisberg 2019, Pinheiro et al. 2021, R Core Team 2021).

rESULTS
One model with two independent variables, eight threepredictor models and two models with four independent variables met the goodness-of-fit criteria ( Table 2). Nine of these 11 models came from the third group of tested relationships that originally included dominant height and quadratic mean diameter as stand-level predictor variables (Table A3). The model by Harrison et al. (1986) (abbreviated "Harr" in Table 2) was the only 2-predictor model that   Abbreviations: GNLS -generalized non-linear least squares, NLS -non-linear least squares, AIC -Akaike Information Criterion, BIC -Bayesian Information Criterion, logLik -log-likelihood, Adj. R 2 -adjusted coefficient of determination, RMSE -Residual standard error (m), d -diameter at breast height (cm), H m -mean stand height (m), a 0 , a 1 , b 0 , θ, η -model parameters, SE -standard error, PRSE % -Parameter Relative Standard Error (%). Table 3. Regression estimates and goodness-of-fit statistics of the selected height-diameter model, fitted by generalized non-linear squared method.
dissimilarity of the slope was not rejected for this data set for ±5% region of similarity. However, as seen in Figure 3A, the lower-size range of this validation data set is underrepresented in the data set used for model parameterization. In addition, the range of the stocking rates of the sampled plantations in Validation Data Set 1 exceeds that of the parameterization data (see minimum basal area and maximum stand density in Table 1); both disparities could have caused the validation outcome. Further investigation showed that, given the ±5% region of similarity for the intercept, the smallest region of indifference that would reject the null hypothesis of dissimilarity of the slope is ±12%, which suggests acceptable accuracy. The residual errors of the predicted dominant and top stand heights, on the other hand, are symmetrically distributed across the range of the predicted heights ( Figures  3B, 3C) and equivalence tests confirmed rejection of the hypotheses of dissimilarity for both model bias and model proportionality (Table 4).

DISCUSSION
Eleven 2-, 3-and 4-predictor models proved adequate to fit the examined height-diameter relationship, with coefficients of determination and residual standard errors varying within narrow ranges. Comparisons of various equations, regarding height-diameter modelling, across multiple species have shown that most of them predict with a similar degree of precision when extensive data are available (Huang et al. 1992, Sonmez 2009, Weiskittel et al. 2011. Stand density (treesˑha -1 ) is considered the most obvious factor affecting the height-diameter relationship (Crecente-Campo et al. 2010) and stand basal area (m 2 ·ha -1 ) is regarded as another measure of stand stocking (Clutter et al. 1983). A study by Nguyen et al. (2019) indicated that the height-diameter relationship of Pinus koraiensis Sieb. et Zucc. differed between stand density levels and therefore could be conditioned on stand density. An option to estimate stand density from the distance to the 3rd neighbour (Priesol 1970, Anuchin 1977 is suggested by the regulations in Bulgaria for the inventory of relatively homogenous even-aged stands, such as the Scots pine plantations. Both stand stocking parameters were considered as predictors in the generalised model formulations that we tested and were included in four of the selected relationships. Studies by other investigators proposed modelling of dominant stand height as a function of mean stand height and density (García 2017, Tarmu et al. 2020, although other results revealed that the mean-dominant height relationship can often be handled only through the functional form of a simple model (e.g., M.L. W.F. 1980, Van Laar andAkça 2007). Staudhammer and LeMay (2000), who developed a generalised height-diameter model of improved accuracy for alder after inclusion of density as an independent variable, explained the observed result with the shadeintolerance of the species coupled with rapid early growth and restricted longevity.
The final comparison of the 11 candidate models of our study elected as the most adequate a modification of Gaffrey's model (Gaffrey 1988) that included two stand attributes: mean stand height and diameter. Adamec (2015) pointed out the advantageous application of such model as compared to the elaboration of a local height-diameter relationship: the number of measured heights necessary to determine the mean height is lower than the total number of measured heights needed to fit the stand-level height curve. Mean stand height can also be regarded as an indicator of the growth stage of an even-aged stand. Unlike age, larger mean height reflects not only the time since the plantation establishment, but also better site quality on the one hand and stronger growth potential due to genetic factors (intrinsic growth rate, resistance and tolerance to adverse conditions) on the other. Consequently, it can as well be viewed as a composite quantitative variable, a product of the interaction between the time since establishment and the growth conditions (Stankova et al. 2016). Stand mean diameter, on the other hand, is a parameter negatively correlated with stand density and therefore its value reflects the stand stocking rate. It can also be viewed as a composite quantitative variable, a product of the interaction between the growth potential of the trees, speeding up the self-thinning process when high, and stand stocking. A study by Zhang et al. (2021) revealed that the inclusion of the interaction effects of stand density and site index could significantly improve the prediction accuracy of the height-diameter model for Larix olgensis Henry. In addition, the modified Gaffrey's model derived here specifies the estimated tree height value according to the social status of the tree as well. Indeed, both addends of the exponent (1-D m /d) and (1/D m -1/d) additionally amplify the height estimates of trees of diameters above the stand average and reduce those of the smaller-sized trees. The model was obtained from Gaffrey's equation (1988) after replacing dominant stand height by a linear function of mean height and constraining the parameter multiplied by (1/D m -1/d) to be equal to 1. This modification removed the bias revealed by the model predecessor and produced significant pa-rameter estimates.
Gaffrey's equation and other models such us those by Sloboda et al. (1993) and Šmelko et al. (1987)   height, ME -model efficiency, , ARE% -absolute value of the relative error, , where y i, ŷ i are experimental and predicted height values of the i-th measurement and yrepresents the mean observed height value; MARE% is the average of the absolute values of the relative errors; Perc0, Perc25, Perc50, Perc75, Perc100 denote the 0 th , 25 th , 50 th , 75 th and 100 th percentile of the relative errors , respectively.
* The approximate joint two one-sided 95% confidence intervals for the slope and intercept are shown and the ±5% intervals of equivalence beneath them in parentheses.
from a simple relationship, very popular in forest modelling and known as Schumacher function or Michailoff function (Gadow and Hui 1999). It possesses all desired features recommended by Yuancai and Parresol (2001) for functions used to model height-diameter relationships: to increase monotonically, to have an upper asymptote, and to have an inflection point. Van Laar and Akça (2007) examined and compared the model with a similar, but monotone increasing relationship without upper asymptote. The authors used both estimated functions to calculate the regression height of the tree with the arithmetic and quadratic mean diameter, the median, the diameter of the tree with the mean volume, the central area basal area tree, the mean derived from the Weise rule, for the quadratic mean of the 100 thickest trees per hectare, and those of the 10th and 90th percentile of the diameter distribution. They found that the differences between the height estimates based on both equations were almost negligible for the central values, but substantial for top height. The estimates of model efficiency and prediction errors as well as the equivalence tests conducted to validate the generalised height-diameter model derived in our study asserted its adequacy in forecasting top and dominant stand heights of Scots pine plantations in Bulgaria.

CONCLUSIONS
Considering the state-of-the-art of forest inventory in Bulgaria, our investigation pursued development of a parsimonious generalised height-diameter model for the Scots pine plantations in the country. After comparison of a number of 2-, 3-and 4-predictor candidate models, a modified form of Gaffrey's model (Gaffrey 1988) was selected. It predicts the height of a tree from the breastheight tree diameter, mean stand height and diameter, and accounts for the tree social status. The model was derived using large and representative data set and its adequacy was substantiated by verification and validation tests. The model can be broadly used, because it requires information readily available from the inventories. Its application in dominant stand height estimation is recommended. Abbreviations: h -tree height (m), d -diameter at breast height (cm), H m -mean stand height (m), D m -quadratic mean stand diameter (cm), b 0 , b 1 , b 2 , b 3 , b 4 , b 5 -model parameters.         Abbreviations: h -tree height (m), d -diameter at breast height (cm), H m -mean stand height (m), D m -quadratic mean stand diameter (cm), H 0dominant stand height (m), G -stand basal area (m 2 ·ha -1 ), N -stand density (trees·ha -1 ), b 0 , b 1 , b 2 , b 3 , b 4 , b 5 -model parameters.