Even-aged or uneven-aged modelling approach? A case for Pinus brutia

Abstract

Context

The past management of Pinus brutia forests in Lebanon has led to diverse stand structures that cannot be easily classified as even-aged (EA) or uneven-aged (UA). Most stands are between these stand types, and they may be called as “semi-even-aged”. This is a very common characteristic throughout the Mediterranean conifer forests and makes the choice between the EA and UA approaches problematic, in both management and modelling. However, previous research has devoted little attention to the performance of growth and yield models when applied to transitional stand structures.

Aims

The aim of this study was to find the best modelling approach and to recommend equations for simulating the dynamics of the semi-even-aged P. brutia stands of Lebanon on an individual-tree basis.

Methods

Fifty sample plots were measured in Lebanon. Individual-tree growth models were fitted to the whole dataset using either UA or EA modelling approach. Models were also fitted using two sub-samples containing the most EA and the most UA plots. The performance and accuracy of the two modelling approaches were evaluated in all three datasets.

Results

The article provides the first complete growth model for uneven-aged P. brutia stands. The EA sub-models presented better statistical fitting. However, the UA sub-models enabled more accurate predictions of wood production and were almost as good as the EA sub-models when predicting stand dynamics of the EA plots. The EA approach provided poor predictions, and the errors were high when it was applied to UA stands.

Conclusions

In structurally complex stands, the UA modelling approach is to be preferred since it predicts the whole stand dynamics more accurately and enables simulations of a broader range of silvicultural treatments.

1 Introduction

Stand structure can be described by the way the trees within a stand are distributed into age and diameter classes (Oliver and Larson 1996) and has traditionally been a key issue when aiming at predicting forest stand dynamics. It is also tightly related to forest management as the current and future state of the forest resources is a basic matter in decision-making in forestry. According to their structure, forest stands are conventionally classified either as even-aged (EA) or as uneven-aged (UA). In even-aged stands, most trees belong to a single age class, and, although their diameters at breast height (dbh) present certain variation, they cluster around the average dbh and their frequencies diminish at larger and smaller diameter classes. On the other hand, uneven-aged stands are represented by trees of different ages and sizes that often follow the reverse J-shaped diameter distribution (Peng 2000). A specific case of uneven-aged stands that is often regarded as the theoretical uneven-aged stand structure are the so-called all-aged stands, in which all age classes are represented. Although uneven-aged stands may comprise only few ages, they may look and behave much like the theoretical all-aged ones. In even-aged forestry, low thinnings (thinnings from below) are the most common silvicultural treatments, and the forest cover is usually removed in the final cutting, whereas uneven-aged forest management is characterised by repeated high thinning (thinnings from above), continuous regeneration and continuous forest cover.

When aiming at predicting forest stand dynamics, the individual-tree growth modelling techniques enable adaptable and thorough simulations of the stand development and, therefore, are more suitable to deal with transitional stand structures (Sterba 2004). Different models are required for the simulation of even-aged or uneven-aged stand dynamics. A frequent set of distance-independent equations for predicting the dynamics of even-aged stands consists of the following sub-models (even-aged modelling approach): site index curves, height-diameter, diameter increment and mortality/survival. Similarly, the dynamics of uneven-aged stands can be predicted and simulated by using the following set of sub-models (uneven-aged modelling approach): ingrowth, height-diameter, diameter increment and mortality/survival (Vanclay 1994). Stand age and dominant height are not suitable variables to be included in the model set for uneven-aged stands, and tree age distribution is rarely available in forest management. Soil and topographical variables can be used instead to characterize site quality (Clutter et al. 1983). On the other hand, ingrowth is not considered in even-aged stand dynamics since all trees are supposed to rise more or less simultaneously.

However, under natural conditions, forest stands present gradations between the two theoretical stand structures (e.g. Smith et al. 1996). It is, in fact, a very common characteristic in the Mediterranean region that forest stands cannot be easily classified as even-aged or as uneven-aged. In many areas dominated by light-demanding species (mainly Pinus sp.), the pure even-aged silvicultural scheme has seldom been followed since large dominant trees have been cut to provide better incomes than low thinnings. As a consequence of the opening of the canopy layer, regeneration was boosted. During the past decades, these periodical selection cuttings have shaped the forests into two-aged or uneven-aged stands (González 2005).

The Pinus brutia forests of Lebanon constitute a good example of such diverse stand structures. This species is the most widespread conifer in the Eastern Mediterranean, where it is ecologically important and produces wood for several purposes (Gezer 1985; Fischer et al. 2008). Despite this, there is a lack of knowledge concerning its stand dynamics. Regardless of being a light-demanding conifer that often regenerates after wild fire and is supposed to be characterised by even-aged dynamics (Boydak 2004), the Lebanese P. brutia forests are neither strictly even- nor uneven-aged due to the past forest management (Assaf 2010). In contrast, most stands are between the two typical stand types, and they may be called as “semi-even-aged”. This particularity makes a straightforward application of one of the modelling approaches questionable. Nevertheless, most models fitted for Pinus species have assumed the stands to be either even-aged using even-aged modelling approach, or uneven-aged using uneven-aged modelling approach. The growth and yield modelling studies on P. brutia have mostly applied the even-aged approach (de-Miguel et al. 2010; Shater et al. 2011) or have provided a partial set of models for uneven-aged stands (Palahí et al. 2008). Notwithstanding, de-Miguel et al. (2010) already highlighted the existence of a certain degree of heterogeneity within the diameter distributions of Lebanese P. brutia stands when comparing them to stands inventoried in Syria. Palahí et al. (2008) also discussed the possibility that the stands in north-eastern Greece were not completely even-aged. Not many authors (e.g. Groot 2002; Sterba 2004) have explicitly dealt with transitional stand structures in general and few studies (if any) have analysed the suitability of each modelling approach for Mediterranean semi-even-aged stands. Hence, to date, there is not enough information to conclude which approach enables more accurate predictions of the stand dynamics of P. brutia or other similar Mediterranean forest systems. As stated by Vanclay (2003), “it is important to test models, to establish their strengths and weaknesses, and to demonstrate to users the range of conditions over which reliable projections can be expected”.

The aim of this study was to find the best modelling approach and to recommend equations for simulating the dynamics of the semi-even-aged P. brutia stands of Lebanon on an individual-tree basis. Furthermore, this study raises new research questions of interest for other forest species and ecosystems.

2 Materials and methods

2.1 Modelling data

Fifty temporary sample plots placed throughout the country were measured so as to capture a wide range of variation in site quality, stand age and stand density (Tables 1 and 2). The sample plots were circular with different radius depending on the stand density. Around 75 trees were measured in the plots. In every plot, the current dbh and diameter growth of the past 20 years divided into two consecutive 10-year growth periods were recorded for each tree using increment coring; tree height and bark thickness were measured for 10–11 sample trees; and age was measured for five dominant trees. In order to assess site quality and productivity, the following variables were also measured in every plot: altitude, slope, aspect, soil type and average soil depth. During the data preparation process, additional stand and tree level variables (i.e. stand basal area, basal area of trees larger than the subject tree) and several transformations of these variables were calculated for every plot or tree. Finally, tree and stand variables were backdated to the beginning of the past 10-year growth periods. As expected, the diameter distributions of the inventoried plots showed a wide diversity of stand structures. Although in some cases the diameters closely followed the expected theoretical distribution of either an even-aged or an uneven-aged stand, the diameter distribution in most plots was semi-even-aged (Fig. 1).

Table 1 Summary of the data used to model diameter increment (7568 observations) and height-diameter relationship (581 observations)

Table 2 Summary of the data used to model ingrowth (100 observations)

Fig. 1

Selection of the 25 most even-aged and 25 most uneven-aged stands based on the criterion SD + 2 SK. The black dots represent the stands classified as uneven-aged and the empty dots represent those stands classified as even-aged. The displayed diameter distributions represent the contour of the dataset and different degrees of semi-even-agedness, that is, different states of the gradation between the theoretical even-aged and uneven-aged stand structure. The dots in the bottom-left “corner” of the cloud are the most even-aged, whereas the dots in the upper-right side of the figure are the most uneven-aged

2.2 Growth models fitting

The site index sub-model to be used in the even-aged modelling approach was obtained from de-Miguel et al. (2010) and includes an indicator variable for the specific country effects of Lebanon. In addition, new equations for 10-year individual-tree diameter increment (with variables describing site quality, tree size and competition) and for height-diameter relationship were fitted to the entire sample using either uneven-aged or even-aged modelling approach. To complete the set of sub-models for the uneven-aged modelling approach, a two-equation ingrowth sub-model (see Table 3) was also fitted to the whole dataset in order to predict the number and mean dbh of recruited trees at the end of a 10-year time step. The diametric threshold for ingrowth was 10 cm.

Table 3 Sets of models available for simulating P. brutia stand dynamics under each modelling approach

Nonlinear regression analysis was applied using R software (R development core team 2011). The criteria considered when evaluating the suitability of each equation were: (a) agreement with current biological knowledge, (b) parsimony and robustness, (c) statistical significance (p value < 0.05), (d) absence of bias, (e) logical behaviour in extrapolations out of the range of the modelling data and in long-term simulations, (f) homoscedasticity and normal distribution of residuals and (g) absence of multicollinearity (except in those cases in which multicollinearity was voluntarily sought due to the specific meaning of some highly correlated variables).

2.3 Evaluation of uneven-aged and even-aged modelling approaches

2.3.1 Simulation procedure

The selected sets of growth sub-models using both uneven-aged and even-aged approaches were used separately to simulate a 20-year growth period in every sample stand using the known backdated stand conditions 20 years ago as the starting point for the simulation process and running the simulation until the current stand conditions. A taper equation (de-Miguel et al. 2011) was used to compute the stem volumes of trees and, by aggregation, the total stand volume. The predictions obtained by the simulation process were then compared with the measured value of current stand volume. A detailed step-by-step procedure for simulating P. brutia stand dynamics by combining the models under an even-aged assumption can be found in Shater et al. (2011). The procedure was slightly different in the UA approach, which included the ingrowth process as an additional element. One 10-year time step is simulated as follows: (1) predict the 10-year diameter increment for each tree and add it to the current diameter, (2) calculate the number and initial diameter of ingrowth trees and add these trees to the stand, and (3) compute the new tree heights based on the height–diameter relationship. Survival was not simulated since backdated characteristics of current survivors were used as input data; there was no mortality in the data.

2.3.2 Preparation of evaluation datasets

The whole 50-plot sample was split into two sub-samples of 25 plots containing, respectively, the most even-aged and the most uneven-aged stands according to the standard deviation (SD) and skewness (SK). A high standard deviation of dbh indicates “uneven-agedness” even if the diameter distribution is bell-shaped. SK describes the degree of asymmetry of a distribution, the typical uneven-aged J-shaped diameter distributions having positive SK values. Plots with high SD and SK are therefore characterised by higher degrees of uneven-agedness. Several classification criteria based on SD and SK were tested, and the standard deviation plus two times the skewness (SD + 2 SK) was finally selected to bisect the plots as even-aged and uneven-aged (Fig. 1). As a result, a 50-plot sample containing all the stands and two 25-plot sub-samples containing the most even-aged and the most uneven-aged stands were obtained to evaluate the performance of the two modelling approaches.

2.3.3 Evaluation methods

The performance of each modelling approach was evaluated by comparing the model-based predictions and the measured values in three different ways: (a) assuming that all the stands were either uneven-aged or even-aged, that is, testing the predictions of each modelling approach against the observed values in all the 50 stands (“overall-performance”); (b) testing the predictions of each approach against the measured values of the 25 stands corresponding to the same stand structure as the approach (“self-performance”), that is, examining the accuracy of both uneven-aged and even-aged approaches on the most uneven-aged and even-aged stands, respectively; and (c) testing the predictions of each approach against the measured values of the 25 stands corresponding to the opposite stand structure (“cross-performance”), that is, analysing the performance of both uneven-aged and even-aged approaches on the most even-aged and uneven-aged stands, respectively.

2.3.4 Statistical analysis

The statistical and graphical analyses were mainly based on the partitioning of the mean square deviation (MSD). Kobayashi and Salam (2000) suggested dividing the MSD into the following three components: squared bias (SB), squared difference between standard deviations (SDSD) and lack of correlation weighted by the standard deviation (LCS). Gauch et al. (2003) proposed a slightly different partitioning, namely SB, nonunity slope (NU) and lack of correlation (LC), arguing that these new components are distinct and additive, with clearer statistical interpretation and transparency related to the regression parameters. These statistics are calculated as follows:

$$ {\text{MSD}} = \frac{{\sum\limits_{{i = 1}}^n {{{\left( {{X_n} - {Y_n}} \right)}^2}} }}{N} $$

(1)

$$ {\text{SB}} = {\left( {\overline X - \overline Y } \right)^2} $$

(2)

$$ {\text{NU}} = {\left( {1 - b} \right)^2}\left( {\sum {\frac{{x_n^2}}{N}} } \right) $$

(3)

$$ {\text{LC}} = \left( {1 - {r^2}} \right)\left( {\sum {\frac{{y_n^2}}{N}} } \right) $$

(4)

where $\overline X$ and $\overline Y$ are the means of the simulated (X) and measured values (Y); b is the slope of the least-squares regression of Y on X; r ² is the square of the correlation coefficient; x _n is equal to ${X_n} - \overline X$, y _n is equal to ${Y_n} - \overline Y$, and N is the number of observations.

This splitting of the MSD enables a proper assessment of the deviations from the perfect equality by analysing separately the sources of discrepancy between the observed and the predicted values, namely translation, rotation and scatter. Furthermore, the more widespread straightforward evaluation of the regression coefficients (intercept and slope) and the Pearson product–moment correlation coefficient of the scatter plot of measured vs. simulated values (e.g. Piñeiro et al. 2008) were also used.

Finally, a multi-criteria decision analysis (MCDA) was done on the basis of quantitative criteria to produce a ranking to facilitate the decision-making concerning the best modelling approach (Render and Stair 1992). The different components and parameters utilized in the above-mentioned model evaluation methods are interrelated. Therefore, to avoid redundancy in the assessment, namely the fact that the same discrepancy from the perfect equality is reflected in more than one criterion, solely the following criteria were used: root of the mean square deviation (RMSD), SB, NU and LC. The same weight was given to all four criteria, and a relative scale was used to rate separately the performance of each criterion. The final scores for each modelling approach were first calculated by aggregation of criteria for each performance level, and, finally, a “global performance” rank was computed by adding the rates of the overall-, cross- and self-performances. The results produced by the MCDA were subsequently contrasted with the information yielded by the other evaluation methods.

3 Results

3.1 Growth models

Table 3 summarizes the two sets of equations representing the two modelling approaches.

3.1.1 Diameter increment equations

In the uneven-aged approach, the residual standard error of the diameter increment model was 1.75 cm. Tree size was represented by diameter (d), and site quality is represented by altitude, slope, soil depth and soil type, whereas competition was represented by stand basal area and basal area in larger trees (BAL). The residual standard error was 1.64 cm for the even-aged approach. In this model, tree size was represented by diameter and age, site quality by site index, slope, soil depth and soil type, and competition by stand basal area and basal area in larger trees.

Variables representing competition or tree size as well as those site predictors which describe challenging conditions for a tree had negative signs in the equations (provoking a decrease in the 10-year growth rate), whereas variables positively related to site quality had a positive sign. The altitude and its logarithmic transformation, considered within the equation of the uneven-aged approach, illustrate an increasing–decreasing trend with a maximum in some point along the altitudinal range. In both models, the ratio BAL/d characterizes the competition status of each individual tree. The dominant trees within a stand are characterised by a low BAL (BAL is equal to zero for the largest tree), and a high d and, therefore, the negative effect on diameter growth in the equation is reduced. Non-dominant and suppressed trees have higher BAL and lower d, which entails a reduction in diameter growth. All the predictors included in the models were highly significant (p value < 0.001).

3.1.2 Height–diameter equations

The height–diameter equation fitted for the uneven-aged approach was an adaptation of the “Hossfeld I modified” function (Peschel 1938). The residual standard error was 2.688 m. The equation fitted for modelling the height–diameter relationship under the even-aged management assumption was a power equation based on Stoffels and Van Soest (1953) modified by Tomé (1989) that forces the model to pass through the point determined by dominant diameter and dominant height. The residual standard error was 2.164 m.

3.1.3 Ingrowth sub-model

The ingrowth sub-model consisted of two equations that predict the number of trees that pass the 10-cm dbh limit during the next 10-year period and the mean diameter of those trees at the end of the 10-year period. Added to the previous models, it enables a complete simulation of uneven-aged stand dynamics. The stand basal area was the main predictor for both the mean diameter of ingrowth trees (D _in) and the number of ingrowth trees per hectare (F _in). The higher the stand basal area, the lower is the number and diameter of ingrowth trees. The residual standard errors were 3.231 cm for the equation to predict D _in and 109.9 trees ha^-1 for the equation to predict F _in.

3.2 Evaluation of the modelling approaches

The results of evaluating the behaviour of each modelling approach when predicting wood production are summarized in Tables 4 and 5. The RMSD is presented instead MSD due to more straightforward interpretation of its units. It also has a geometric interpretation since it equals with the standard deviation of the deviations around the 1:1 line (perfect equality line) in a plot depicting the predicted against the observed values. The discrepancies from perfect fit between simulated and measured values for each type of performance are illustrated using regression lines in Fig. 2.

Table 4 Root mean square deviation (RMSD) in cubic metres per hectare and its partitioning into the three components SB, NU and LC for each modelling approach (uneven-aged and even-aged) and for each performance type

Table 5 Bias (cubic metres per hectare) and parameters obtained from the ordinary least-squares regressions of measured vs. simulated values, namely intercept (a), slope (b) and correlation coefficient (r), for each modelling approach and for each performance type

Fig. 2

Discrepancy between the regressions of measured vs. predicted values (dashed line) and the 1:1 perfect equality line (solid line) for each modelling approach (uneven-aged or even-aged) and separately for the overall-, cross- and self-performance

3.2.1 Overall-performance

When the measured and the model-based values were compared assuming all 50 stands either as uneven-aged or as even-aged (overall-performance), the RMSD was practically the same for both modelling approaches (only 0.21 m³ ha⁻¹ difference). However, the predictions based on the EA approach were much more biased (higher SB). The EA approach led to considerable underestimation (positive bias) of wood production mainly in stands with rather low stocking, whereas the UA approach tended to overestimate production (negative bias) in stands with intermediate and high stocking. The UA models met much better the nonunity slope principle (smaller NU). In contrast, the LC was lower when the stand dynamics were simulated using the even-aged assumption. The regression line of the UA approach (Fig. 2, top left) almost crossed the origin, and its slope was closer to 1. On the contrary, the correlation coefficient was higher for the EA approach.

3.2.2 Cross-performance

When the predictions based on one of the modelling approaches were compared with the measured values of plots representing the opposite stand structure (e.g., EA models in UA plots), the RMSD was considerably lower (better) for the UA approach. Moreover, this approach was less biased and also performed better according to the NU criterion. Similar to the overall performance, the EA approach tended to underestimate wood production (mainly in intermediate and low stocking stands), whereas the UA one tended to overestimate it. The LC was higher for the UA modelling approach. When using the UA approach, the regression line was closer to the origin and its slope was also closer to 1. In contrast, the correlation coefficient was higher for the EA approach.

3.2.3 Self-performance

When the predictions based on one of the modelling approaches were compared with the measured values of those stands representing the same stand structure as the modelling approach (e.g., EA models in EA plots), the RMSD was considerably lower (better) for the EA approach. However, this approach was much more biased (underestimation) and performed worse also with respect to the NU criterion. In fact, the simulation based on the UA set of models was almost non-biased. The LC was again higher for the UA approach. Consequently, the regression line was closer to the origin, and the slope was closer to 1 when the UA modelling approach was used to predict the most UA stands. Although the correlation coefficient was still higher when modelling EA stands with the EA modelling approach, the r of the UA approach was much closer than in the previous performance categories.

3.2.4 Ranking of modelling approaches on the basis of the MCDA

The UA modelling approach was ranked better when analysing the overall performance due to its smaller SB and NU. The RMSD got the same score in both modelling approaches, but the EA one had lower LC. The analysis of the cross-performance also ranked the UA approach better due to its lower RMSD, SB and NU. Nevertheless, when analysing the self-performance, both modelling approaches were ranked equal, as the EA approach presented lower RMSD and LC, but the UA one was less biased and better met the NU criterion. The global performance (aggregation of the scores obtained for each performance type) was better for the UA modelling approach, i.e., it turned out to be the best way to simulate and predict the stand dynamics of P. brutia semi-even-aged stands (Fig. 3).

Fig. 3

Performance ranking of each modelling approach when predicting wood production of semi-even-aged P. brutia stands in Lebanon

4 Discussion

4.1 Growth models

This article, together with some earlier models (de-Miguel et al. 2010), presents a complete set of models enabling simulation and prediction of uneven-aged stand dynamics of P. brutia. The weakest part of the model set is survival since individual-tree survival models are missing. The existing model (de-Miguel et al. 2010) is a self-thinning equation which is the most suitable in strictly EA stands. However, since there is mortality already before the self-thinning limit is reached, the best model type for both approaches would be an individual-tree survival model. Since this model would most probably be very similar in both approaches, the lack of mortality in the data and simulations of the current study should have no impact on the conclusions.

To our knowledge, this study is also the first systematic evaluation of the predictive capability of the EA vs. UA modelling approaches in forests characterized by semi-even-aged stand structures. The statistical fitting of the growth equations was always better for the EA approach. This is understandable as the EA models use more measured information such as stand age and dominant height. The combination of both variables results in the site index concept which has been used as a smart predictor of site quality in EA forest management (Skovsgaard and Vanclay 2008). Stand age, dominant height and site index were used as predictors in all the equations of the EA modelling approach. As UA stands, by definition, do not have stand age, tree size was represented just by tree diameter, and site quality was assessed using several variables within a geo-centric approach, which showed less predictive capacity.

The interaction between dolomitic sand and sand soil types with the stand basal area was highly significant and always negative in all diameter increment equations. A plausible interpretation is that sandy soil constitutes a limiting factor for P. brutia development, and, when the competition within a stand increases (represented by higher G), the growing conditions become still more challenging. In Syrian P. brutia forests, diameter increment correlated significantly and positively with two parent rock types, namely dolomitic limestone and radiolarite (Shater et al. 2011). This agrees with the Lebanese results since most plots where the soil type was not sand were on calcareous soil. The most favourable sites for P. brutia growth are those characterised by gentle slopes, deeper soils and, in accordance with previous research (e.g. Fontaine et al. 2007), an altitude around 500–600 m above sea level. Eastern aspect seemed to indicate better growing conditions for this species. However, aspect was not among the most significant variables, and the sample was unbalanced due to small representation of eastern aspects. Therefore, aspect was not used as a predictor in the final version of the equations.

4.2 Performance of the modelling approaches

Despite the better statistical fitting of the sub-models based on the EA approach, the set of sub-models fitted under the UA assumptions performed better when simulating the whole stand dynamics. Basically, the results revealed that, if one single modelling approach has to be chosen, the models for UA stands enable more accurate predictions of wood production. The EA modelling approach provided poorer predictions in simulation, and the errors were very high when the approach was applied to uneven-aged stands. It could be concluded that the trade-off between using models with slightly higher statistical accuracy (e.g. use of site index as predictor) and using less accurate models but describing with higher fidelity the inherent stand dynamics (e.g. use of an ingrowth model) is clearly positive for the latter option. This might be caused by the fact that the assumptions that underlie the EA modelling approach are more constraining than those underlying the UA one. As a result, the UA approach is more flexible to slight deviations from the theoretical foundations of the method and, therefore, is more suitable to predict the stand dynamics of semi-even-aged stands. Applying the principle of caution, the use of the UA modelling approach might be safer to ensure sound simulations and predictions of stand dynamics.

In order to dispel any doubts about the soundness of these results, it is worth mentioning that the MCDA settings were, in the end, quite conservative to rank the UA approach as best. As pointed out by Gauch et al. (2003), when the purpose of a model is prediction, “this imposes special interest in the 1:1 equality line”. Those criteria in which the UA approach performed better, namely SB and NU (Table 4), are the ones most related to the perfect equality line. The LC, that is, the third component of the MSD, simply represents the dispersion around the regression line, but it does not describe critical deviations from the 1:1 line. In fact, the lower values of LC systematically found in the EA approach are just reflecting the fact of the better statistical fit of those models due to, as explained above, the use of more information, namely stand age and total height growth. Therefore, if higher weights had been assigned to criteria closely related to the 1:1 line (SB and NU), the balance would have shifted even more in favour of the UA approach. This can be perceived also from the standard output of the regression of observed vs. predicted values (Table 5): The intercept is closer to 0; the slope is closer to 1, and the correlation coefficient is lower in the UA approach for any class of performance. The SDSD and LCS of Kobayashi and Salam (2000) were also calculated but not reported since they had no clear interpretation as already pointed out by Gauch et al. (2003). To further prove the soundness of these results, additional growth models were fitted under both EA and UA assumptions using separately the two 25-plot sub-samples as modelling data and as independent testing dataset. This additional verification serves as a validation test since the performance of each modelling approach in an independent dataset was similar as when using the whole 50-plot dataset, namely lower bias and slope closer to 1 for the uneven-aged approach.

One of the main concerns during the analyses was the method for splitting the whole dataset into two 25-plot sub-samples for evaluation purposes. Some authors (e.g. Buongiorno et al. 2004) have used the Shannon–Wiener entropy index to assess the structural diversity. However, this index depends on the width of diameter class. Moreover, such a method focuses only on the number of diameter classes and the amount of trees per diameter class which may not be sufficient to describe whether a diameter distribution is bell-shaped or J-shaped. In contrast, the criterion finally applied (SD + 2 SK) after testing several combinations of criteria including different weights of SD and SK, also the Kurtosis coefficient, successfully classified the stands either as the most even-aged and most uneven-aged. This was graphically confirmed from the diameter distributions of the classified stands as presented in Fig. 1. In addition, another indicator of the efficiency of the classification was the fact that both modelling approaches were equally ranked in the self-performance, whereas they were unequal in the other performance types.

4.3 Additional considerations

Forest modelling could be defined as the art of understanding, emulating and predicting forest stand dynamics. The systematic analysis carried out in this study enabled us to propose the best modelling approach and to recommend equations for simulating the dynamics of the typical semi-even-aged P. brutia stands of Lebanon on an individual-tree basis. As mentioned above, the UA modelling approach was regarded better in the prediction of wood production, which suggests that the underlying natural stand dynamics are better emulated by this approach. However, other aspects may be considered when choosing the set of models. For instance, it may be appropriate that the models are able to simulate the effect of diverse silvicultural practices on the stand dynamics. In this respect, the UA approach is also to be preferred as it can be used to simulate any kind of thinning. On the contrary, the use of the EA set of models assumes that the dominant trees are not removed in thinnings. Consequently, the models may not simulate properly the development of a stand when high thinnings are used.

Furthermore, the used modelling approach should reflect the real forest conditions and forestry practices to be useful as a decision support tool when comparing forest management alternatives. Therefore, changes in forest management policies may also affect the suitability of a modelling approach. Three major alternative forest management scenarios may be envisaged from the current semi-even-aged stand structures: (1) to move towards an even-aged forestry, (2) to maintain the status quo of semi-even-agedness and (3) to move towards an uneven-aged forestry. If even-aged forestry is adopted, the models developed by de-Miguel et al. (2010) for Middle East can be recommended to be applied in Lebanon. However, this would be an unlikely development in Lebanon and in the other countries in Middle East. If the current stand structures were expected to be maintained or conventional UA forestry is practiced, then the set of models presented in this article based on the UA modelling approach, would be the most suitable tool for predicting stand dynamics.

Further improvements of this set of models might include new equations in order to, for instance, take into account different kinds of risks (see Jactel et al. 2009) or to properly simulate the regeneration dynamics of P. brutia, especially in the absence of fire (Fyllas et al. 2008). Models to predict mortality or survival on an individual-tree basis are also needed. Finally, another interesting question emerges concerning explicitly the suitability of each modelling approach to be applied on transitional stand structures. What degree of semi-even-agedness constitutes a realistic threshold determining the suitability of a modelling approach? Or, in other words, what is the case-by-case optimal modelling method providing more accurate predictions of stand dynamics? This questioning can be expanded to the broader range of those forest species and ecosystems which often show structural diversity or semi-even-agedness.