Modeling dominant height growth using permanent plot data for Pinus brutia stands in the Eastern Mediterranean region

Aim of study: At current, forest management in the Eastern Mediterranean region is largely based on experience rather than on management plans. To support the development of such plans, this study develops and compares site index equations for pure even-aged Pinus brutia stands in Syria using base-age invariant techniques that realistically describe dominant height growth. Materials and methods: Data on top height and stand age were obtained in 2008 and 2016 from 80 permanent plots capturing the whole range of variation in site conditions, stand age and stand density. Both the Algebraic Difference Approach (ADA) and the Generalized Algebraic Difference Approach (GADA) were used to fit eight generalized algebraic difference equations in order to identify the one which describes the data best. For this, 61 permanent plots were used for model calibration and 19 plots for validation. Main results: According to both biological plausibility and model accuracy, the so-called Sloboda equation based on the GADA approach showed the best performance. Research highlights: The study provides a solid classification and comparison of Pinus brutia stands growing in the Eastern Mediterranean region and can thus be used to support sustainable forest management planning.


Introduction
The classification of forest stands in terms of their productivity is highly important for foresters as management practices depend on the site quality, which is also known as site index (Álvarez-González et al., 2005). Due to the difficulty of estimating site productivity via direct measures depending on abiotic factors such as light conditions, heat, soil moisture, nutrient availability, soil depth and aeration (Stearns-Smith, 2001), indirect methods are more frequently used. Site index equations based on the height growth of dominant trees are the classical way of indirectly estimating site quality in forest management. In this sense, a site index is defined as the average height of the co-/dominant trees per unit area at a given stand age (Álvarez-González et al., 2005). The dominant height at this age is relatively stable and the least affected by stand density or thinning treatments (Assmann, 1961;Elfving & Kiviste, 1997). Three main methods can be used to fit site index equations to data, being: (1) a guide curve, (2) parameter prediction, and (3) differential equation method (Clutter et al., 1983).
For developing site index equations, the differential equation approach has been widely used (e.g., Álvarez-González et al., 2005). It can be applied to any height-age equation to produce families of anamorphic or polymorphic curves (Burkhart & Tome, 2012). The method requires data from permanent plots or stem analyses, which are considered highly dynamic though expensive (Clutter et al., 1983). It is important to note that site index equations can be static or dynamic. While static models assume a fixed base-age, i.e. a specific age at which site index is quantified, dynamic models are base-age invariant (Cieszewski, 2001). The latter category of models provides a set of desirable properties such as biological plausibility (i.e. height should be zero at age zero and equal to site index at reference age), polymorphism, multiple horizontal asymptotes, logical behavior, a single inflection point (Elfving & Kiviste, 1997) and an increasing curve behavior over time (Diéguez-Aranda et al., 2005). To get one or more of these desirable characteristics, two approaches are used: the Algebraic Difference Approach (ADA) (Bailey & Clutter, 1974) and the Generalized Algebraic Difference Approach (GADA) (Cieszewski & Bailey, 2000). The main limitation of ADA is that it produces a polymorphic site index equation with a single asymptote for all sites (Cieszewski, 2001), whereas GADA expands base equations with various theories about growth characteristics (e.g., multiple asymptotes, growth rate), and is thereby able to reflect local, site-specific height growth trajectories (Cieszewski & Bailey, 2000).
Pinus brutia is an important natural tree species of the eastern Mediterranean countries that covers about 55,000 ha in Syria and 17,000 ha in Lebanon (MFWA, 2012). In addition to a wide-ranging utilization of its wood in various wood-based industries and its' environmental importance, e.g., for soil protection, P. brutia produces resin, and is used in reforestation and restoration programs because of its good performance on poor sites (Nahal & Zahoueh, 2005).
Sustainable forest management of P. brutia stands is not widely practiced in the region; instead of being based on plans, experience and simple calculations often guide forest management. Effective management of P. brutia forests requires stratifying forest areas into productivity classes and building forest growth and yield models, as they allow determining annual growth, predicting future stand development and comparing different silvicultural treatments. Given that site index is among the main factors that influence the growth and yield of trees, one of the first steps to develop growth and yield models and determining the economic and environmental effects of forest management is the construction of reliable site index equations (Pretzsch, 2009;Skovsgaard & Vanclay, 2013). A set of site index equations has been developed for P. brutia, for example for Turkey (Misir & Misir, 2007) as well as for Greece and Cyprus (Kitikidou et al., 2011(Kitikidou et al., , 2012. In the Eastern Mediterranean region, the only published site index equation for P. brutia used the guide curve method that produces anamorphic curves that represent dominant height-age relationships only to a limited extent (Suliman, 2013). Hence, this study adopts a (G)ADA approach as a more realistic prediction method for site index of pure and even-aged P. brutia stands in the Eastern Mediterranean region. More specifically, a comparison of eight potential site index equations was performed. It is expected that equations based on GADA outperform those based on ADA because of its more dynamic and flexible characteristics.

Study area
This study focuses on pure even-aged Pinus brutia plantations in the coastal region of Syria (Fig. 1), where 38.7% of the Syrian forests is located and P. brutia is the dominating tree species. The elevation in the region ranges between 0 and 1600 m and the climate can be characterized as Mediterranean. Over the past ten years, the annual precipitation varied between 600 and 1200 mm; the mean annual temperature between 14 and 16°C (Drought and Natural Disasters Fund Directorate, 2016).

Data
Inventory data from 80 plots that were established in P. brutia stands was used. All sites, except for four, were unmanaged. At these four sites, management interventions took place before the year 2000 and were based on experience rather than on a management plan. The plots (Table S1 [suppl.]), which were mainly of circular shape (except for five rectangular ones), were measured in 2008 and 2016. Plots varied in size from 70 to 1963 m 2 depending upon stand density, so that about 50-75 trees were measured per plot. Stand density ranged between 137 and 3144 and the average was 1036 trees ha -1 . Stand age varied between 16 years and 129 years; and the average was 56 years in the first inventory.
In each plot, all trees were numbered, stand age was determined and diameter at breast height (d 1.3 ) measured. For a selection of trees (i.e. 10-11 per plot), tree height was measured as well. Using these data, the quadratic mean diameter (Dq) was calculated and used to derive the top stand height (H 100 ) following Michailow (1943). Top stand height estimations were done separately for groups of 61 and 19 plots that were later on used for model calibration and validation (Table S1 [suppl.]). Calibration and validation plots were selected in such a way that both groups are comparable in terms of stand characteristics.

Fitting the site index equation
To fit a site index equation, a base-age-invariant method (Bailey & Clutter, 1974) was used. The approach considers that all height trajectories follow a similar functional form with parameters shared by all individuals (global parameters) and parameters specific for each individual (local parameters).
A difference algebraic form of the height-age relationship (i.e. differential equation) was developed, where H 100 (t 2 ) is expressed as an equation of the re-measurement age (t 2 ), the initial age (t 1 ) and the height at the initial measurement H 100 (t 1 ). We selected eight base-invariant equations (M1-M8, see Table S2 [suppl.] for model equations) that are commonly used in forest research (Sloboda, 1971;McDill & Amateis, 1992;Cieszewski & Strub, 2008;Kitikidou et al., 2011;Sharma, 2013;Rojo Alboreca et al., 2017), half of which based on ADA and the other half based on GADA, and compared their performance. Parameters were estimated with non-linear least square regression analysis, applying the Levenberg-Marquardt algorithm due to correlation among parameter estimates (More, 1977).

Model selection criteria
Model performance was assessed while focusing on: (1) the biological meaning of the model and its parameters (e.g., polymorphism, sigmoid growth pattern, and multiple asymptotes), (2) the behavior of model residuals, and (3) the evaluation of statistical indices that describe the goodness-of-fit (Amaro et al., 1998). For the latter, model bias, relative bias, model precision, relative model precision, model accuracy and relative model accuracy were used (cf. Pretzsch, 2009), while the coefficient of determination (R 2 ) was used as an index for model efficiency (Table S3 [ suppl.]). Data analysis was conducted with R 3.4.0 (R Core Team, 2017).

Results and Discussion
In this study, data on top height and stand age were obtained in 2008 and 2016 from 80 permanent plots from pure even-aged Pinus brutia stands in the coastal region of Syria to develop site index equations that can be used in yield modeling. Eight dynamic site index equations were developed (M1-M8, Table S2 [suppl.]). The analysis of the fit statistics revealed that there is no serious breakdown of homoscedasticity (Fig. S1 [suppl.]). Though effects of stand density on dominant height growth have been reported frequently (Nilsson et al., 2010), this could not be confirmed in our study given that relative bias values were particularly low (Table S4 [suppl.]). A likely explanation is that height growth of dominant trees is unaffected by competition for a large range of stand densities, confirming the validity of the recommendation of Assmann (1961) to use dominant height growth as indicator of stand productivity.
Further, it was found that all equations performed similarly for calibration and validation data with a slight superiority of the base-age invariant Sloboda equation (M5, Fig. 2, Table S4 [suppl.]). As the other GADA equations, M5 is characterized by multiple asymptotes for different site index values (i.e. polymorphism). Given that 31% of our plots had a stand age over 70 years, the inclusion of asymptotes was needed for a realistic description of the relationship between age and tree height at older stand age. As described by the hydraulic limitation hypothesis (Ryan & Yoder 1997), photosynthetic rates drop when trees get older or grow taller, which is reflected by an asymptote in the height curve. ADA equations are conceptually able to model such a single asymptote. In reality, however, there are also factors other than hydraulic limitation that affect height growth, e.g., nutrient availability (Manso et al., 2021). These factors make that asymptotes may differ between forest stands. GADA is able to model multiple asymptotes, and can therefore be considered as an ecologically more meaningful approach. The culmination of height increment was found to differ with site class, and for M5 occurred between 9 and 14 years in the site classes 30 to 10 ( Fig. S3 [suppl.]). Competition is a likely cause of the earlier culmination of height increment on sites with higher site class, as increased competition at such sites may promote height increment at young age.
A prominent feature of all developed equations is that they contain the initial age and initial height as independent variables in addition to the standard age. This feature allows computing dominant heights for different standard ages and site indices from the same equation. Nevertheless, a standard age of 50 years is recommended in this study, because the reference age should be close to the rotation age and should be as young as possible in order to help in earlier decision making of the silvicultural treatments (Diéguez-Aranda et al., 2005). Besides, a standard age of 50 years is commonly used for coniferous species (Suliman, 2013;Gezer, 1986), and the relative error according to a test proposed by Diéguez-Aranda et al. (2005) was lowest at this value ( Fig. S4 [suppl.]).
A comparison of M5 with previously developed site index equations for P. brutia in Syria (Suliman, 2013) and Cyprus (Kitikidou et al., 2012) revealed that at an age of 100 years, the site class curves of M5 reached slightly higher heights than those developed for the Rabiaa region in Syria using the Richard-Chapman function (Fig. S5a [suppl.]). Curves developed for central Cyprus using the Korf function showed generally lower growth rates at younger and older ages (Fig. S5b [suppl.]), which might be related to higher precipitation amounts in Cyprus as compared to Syria.
Irrespective of the site index equation tested, younger generation trees showed a stronger height increment as trees from older generations (Fig. S2 [suppl.]), which may be attributed to the climatic component of site index (Sharma, 2013), to factors such as eutrophication (cf. Pretzsch et al., 2014), or could be related to the natural and anthropogenic selection of the best adapted genotypes (Matisons et al., 2017). More thorough studies on possible causes are, however, needed.
In conclusion, this study presented a dynamic base-age invariant site index equation (Sloboda equation) for P. brutia in Syria that can serve as a baseline for classifying and comparing stands in the wider Mediterranean basin. In comparison to other equations, it provides more reliable estimates and can be used in developing individual-tree growth and mortality models. The relatively simple form and convenient use makes it furthermore suitable as a scientific basis for decision support in forestry and land-use planning.