Constructing site quality curves and productivity assessment for uneven-aged and mixed stands of oriental beech (Fagus oriental Lipsky) in Hyrcanian forest, Iran

ABSTRACT In the domain of sustainable forest management, site productivity assessment is a major forestry topic. The reliable estimates of site quality are crucial for improved predictions of timber yields and for meaningful simulation studies. Site quality curves were constructed and assessed for beech forests in a Hyrcanian forest, located in the north of Iran. A measure of the site quality concept known as site form (i.e. the expected mean height at 75 cm diameter at breast height) was used. One hundred and twenty seven 0.1 ha circular sample plots were established covering soil and topographic variations in Tarbiat Modares University Experimental Forest. Site productivity curves were constructed using a difference equation approach in the context of the Chapman-Richards model and the Edminster el al. (1991) method. After estimating site form for each plot, they were divided into three categories (poor, medium, and good). Ordinal logistic regression model of site form on soil and topographic variables were fitted using the free statistical programming language ‘R’. The results showed that only the percentage of clay, altitude, and phosphorous were significant at 0.05 significance level (Nagelkerke pseudo-R-squared = 51.58%). The results of this study suggest that site form appears to be a potentially useful indicator of site quality in mixed and uneven-aged forests.


Introduction
Measuring site productivity is an important issue in the domain of sustainable forest management (Bontemps and Bouriaud 2014). It is obvious that some sites support luxuriant forests whilst others are capable only of supporting poor forests (Vanclay 1992); hence, one of the first tasks of the forest growth modeler may be to find a method to quantify the site differences. Variation in site productivity may be assessed within the conceptual framework of site quality, defined as the capacity of a ground to produce vegetation of given characteristics, as determined by existing environmental conditions (Carmean 1975;Hagglund 1981). Forest management needs a reliable measure of site productivity in order to forecast the growth and yield (Vanclay 1994), because any bias in its estimation may transmit through the growth, mortality, and recruitment functions and affect all the modeling results. Any single estimate of site productivity must be accurate, because it summarizes several multi-dimensional factors of the environment as a single index. Therefore, the challenge is to find the most useful index.
There are many ways to estimate the productivity of a site including forest dimensions (Clutter et al. 1983;Vanclay 1994) and abiotic factors (Clutter et al. 1983;Schonau and Aldworth 1991;Vanclay 1994), as well as the visual appearance of a stand, indicator species, and composition of ground vegetation (Daubenmire 1976;Schafer 1989), and a comparison of their merits requires some definitions and a classification of the major alternatives. Most of the site quality studies have been done in even-aged stands (Herrera et al. 1999); however, the method which is most frequently used in these stands is the estimation of site index (Carmean 1975;Ortega and Montero 1988;Herrera et al. 1999), being simple to measure and a good predictor of the utilizable production from the site (Herrera et al. 1999). But several problems may arise in some forests. Sometimes some species such as Fagus orientalis form false annual growth rings and attempts to determine the age are quite subjective. Other problems include the disparity response of this index to climatic factors and the dependence of site index on forest density (Bontemps and Bouriaud 2014). For these reasons, the site index cannot be used to quantify the site productivity and some other measures are needed. On the other hand, in unevenaged and mixed stands, tree height in relation to age cannot be used to express the site quality. The height growth of a species in this type of stand is not closely related to age but more to varying stand conditions by which it has been affected during its life. McLintock and Bickford (1957) considered several alternatives for evaluating the site quality in uneven-aged stands in their study of site quality for red spruce in the northeastern United States. Their results showed that the relationship between tree height and diameter at breast height (DBH) in a stand was the most sensitive and reliable measure of site quality. Site index according to this concept is then defined as the height attained by dominant trees at a standard DBH (Husch et al. 1963). The height at a defined index diameter was used as a measure of site quality and was named site form, to avoid confusion with site index resulting from height-age relationships. Vanclay and Henry (1988) concluded that site form has the following desirable features which can be considered as good measures of site productivity: (1) reproducible and consistent over long periods; (2) indicative of site and not influenced by stand condition or management history; (3) correlated with the site's productive potential; and (4) simply determined by normal inventory measurements which mean that site form estimations do not require data from permanent plots.
Reliable estimates of site productivity are essential for the sustainable management of forest resources (Skovsgaard and Vanclay 2013). Site productivity depends both on natural factors inherent to the site and on management related factors (Skovsgaard and Vanclay 2013). Two primary factors affecting forest site productivity are soil and topography characteristics (Hererra et al. 1999). Many studies investigated the correlation between soil and topographic factors with site quality (e.g. Brown and Loewenstein 1978;Monserud et al. 1990;Hererra et al. 1999, Herrera-Fern andez et al. 2004Aertsen et al. 2010Aertsen et al. , 2011. One of the most abundant and economic hardwood genera in Hyrcanian forests is the beech genus. Because of the importance of oriental beech (Fagus orientalis Lipsky) as one of the main timber species of Hyrcanian forest, evaluating its site productivity and identifying the most important factors influencing it, may help foresters to manage these forests sustainably. This study aims at developing site quality curves for beech forests in a mixed and uneven-aged Hyrcanian forest located in the north of Iran. This study takes advantage of the modified guide curve methodology developed by Edminster et al. (1991) and height-diameter relationship suggested by Vanclay and Henry (1988).

Study site
This research was conducted in Tarbiat Modares University (TMU) experimental forest station located in a temperate forest of Mazandaran province in the north of Iran (36 31 0 56 00 -36 32 0 11 00 N, 51 47 0 49 00 -51 47 0 56 00 E). TMU forest extends from 100 to 1700 m above sea level (a.s.l.). Mean annual rainfall is 1500 mm with maximum and minimum falls in October and June, respectively. Temperatures vary from a mean monthly minimum of 6.6 C in December to a maximum of 25 C in June, with an annual mean of 15.9 C. The study area's forests are mixed and uneven-aged, dominated by Fagus orientalis associated with Carpinus betulus, Acer velutinum, Parottia persica (DC.) C. A. Mey., Sorbus torminalis L., Quercus castaneifolia, Alnus subcordata C.A.Mey., Acer leatum C.A.Mey., Prunus avium L., Ulmus glabra L., and Tilia begonifolia Stev. These forests are managed as close-to-nature with single tree harvesting methods. Beech forests accounts for c. 17.6% of the total forest area, 30% of the standing volume and 23.6% of the stem number in Hyrcanian forests in Iran. The average beech volume per ha ranges between 480 m 3 and 740 m 3 in pure stands and 600 m 3 and 700 m 3 in mixed stands (Sagheb-Talebi et al. 2004).

Sampling
Data was collected from a systematic network of temporary sample plots established in 2014. As part of the systematic sampling approach, two grids were drawn and the interval between grid lines was set to 100 m in both longitudinal and latitudinal directions. One hundred and twenty seven 0.1 ha circular sample plots were laid out for this study. The following criteria were taken into account for plot establishment: it is tried sampling of all soil and topographic gradients in the research area: homogeneous soil and topographic conditions within each sample plot with no evidence of disturbance in the plot. The plots extended over the altitudinal belt from 1000 to 1500 m a.s.l.
All trees with DBH > 7 cm within the plots were measured using calipers to the nearest millimeter and total height of all beech trees was measured by a Vertex IV height meter with accuracy of 0.1 m. Forked trees or those with damaged tops were not included in the analysis. For quantifying nutrient availability, five soil samples (0-10 cm) were randomly taken within each plot. Soil samples were mixed and analyzed in the laboratory. Soil variables were determined by the following methods: bulk density (by clod method ;Plaster 1985); texture (by Bouyoucos hydrometer method; Bouyoucos 1962); pH (soil:water ratio 1:2.5); total organic C (by Walkley and Black method; Allison 1975); total N (by Kjeldahl method; Bremner and Mulvaney 1982); available phosphorous by using the Olsen method (Homer and Pratt 1961); and available K by a flame atomic absorption spectrophotometer (AA500F, PG Instruments Ltd, China). Physiographic variables such as slope (in %), aspect (u), and elevation were recorded in each plot. Aspect, as the azimuth measured from true north, was transformed to a topographic radiation index using the equation TRASP = [1cos((p/180)(u-30))]/2 (Aertsen et al. 2010).
Methods for constructing a site index system Ahmadi et al. (2013) assessed six non-linear growth functions (i.e. Chapman-Richards, Schnute, Lundqvist/Korf, Weibull, Modified Logistic and Exponential) to tree height-diameter data of oriental beech in the study area. The predictive performance of these models showed that the Chapman-Richards function had better predictive performance than others for predicting the height of oriental beech trees. According to Clutter et al. (1983) most of the approaches used to fit the site index curves can be viewed as special cases of three general development techniques: (1) the guide-curve method; (2) the parameter-prediction method; and (3) the difference equation method.
The guide curve method assumes proportionality (anamorphic) among curves for different site qualities and is used to generate a set of anamorphic site productivity index curves. The parameter prediction method is based on fitting a growth function tree-by-tree or plot-by-plot and relating the parameters of the fitted curves to the site productivity index. The difference equation method is based on the fact that observations of the same plot or dominant tree should belong to the same site productivity index curve. A difference algebraic form of a height-DBH or differential equation is developed where height at re-measurement (H 2 ) is expressed as a function of the re-measurement DBH (DBH 2 ), the initial DBH (DBH 1 ), and the height at the initial measurement (H 1 ). The algebraic difference form is obtained through substitution of one parameter in the differential equation. A modified guide curve methodology developed by Edminster et al. (1991) was used to develop site productivity index curves for oriental beech in this study.
After selecting the best nonlinear regression model, the next step in developing the family of site index curves is to examine the assumption that height-DBH curves for different site indexes are proportional to the shape of the guide curve, using the method employing the coefficient of variation of height suggested by Chapman and Meyer (1949).
The assumption of proportional shape was not supported since there was a definite curvilinear downward trend in a plot of the coefficient of variation of height with increasing DBH. For developing the final height growth guide curve, the lower 95% prediction interval (HL) was constructed using the oriental beech height-DBH data and the chosen model. Applying the method described by Edminster et al. (1991), site productivity index curves for oriental beech were then constructed by determining the distance from the guide curve to both the lower and upper 95% confidence bounds for the reference DBH of 75 cm. The distance proportion (P) was calculated as the difference between the guide curve and decadal site productivity index height values, divided by the difference between the guide curve and the lower 95% confidence bound for each height-DBH sample as represented by the equation: where HG = height (m) of the guide curve at the sample DBH (cm); HT = total height (m) of the sample tree; and HL = height (m) of the lower 95% confidence bound at the sample DBH (cm). Finally, given P, breast height diameter, and height for any measured site tree, the equation for estimation of site productivity index is represented by: where SF = site form; HG 75 = height of the guide curve at DBH 75 cm; and HL 75 = height of the lower 95% confidence bound at DBH 75 cm. After calculating the site form in each plot based on mean DBH and height, classInt package (Bivand 2015) with Fisher style for classifying the site form values into three classes (poor, medium, and good) in R software (R Core Team 2015) was used. The classInt package is commonly used for choosing univariate class intervals for mapping or other analysis purposes. The Fisher style uses the algorithm proposed by Fisher (1958) and discussed by Slocum et al. (2005) as the Fisher-Jenks algorithm.

Ordinal regression
When the dependent variables are ordinal there are many options to build a model. Since the response variable is a set of ordered site quality categories (poor, medium, and good), ordinal logistic regression (often just called 'ordinal regression') could be used to predict the dependent variable given one or more independent variables. Ordinal logistic regression was fit using polr function in MASS package (Venables and Ripley 2002). As in linear regression, simple R-squared cannot be used in ordinal regression. R-squared gives the information about how much variance is explained by the independent variables. However, variance is split into categories, hence Cox and Snell's, Nagelkerke's, and McFadden's pseudo-R 2 statistics are used in ordinal regression for estimating the variance explained by the predictors. For this study Nagelkerke's pseudo-R 2 was reported.

Results
After fitting data to the height-DBH data using the Chapman-Richards model, the coefficients for the final guide curve equation were calculated and obtained as bl = 47.773, b2 = 0.011 and b3 = 0.611. Regression results were significant (P < 0.05) with an approximate R 2 value of 75.5 and an approximate RMSE of 3.30. Residuals from the guide curve regression were normally distributed.
As mentioned, the assumption of proportional shape was not supported as there was a definite curvilinear downward trend in a plot of the coefficient of variation of height with increasing DBH. Hence, site productivity index curves for oriental beech were constructed using lower and upper 95% confidence bounds for the reference DBH. The height growth guide curve, bounds of the lower and upper 95% confidence intervals for individual predictions, and height-DBH data used to develop the oriental beech site productivity index curves are illustrated in Figure 1.
At the reference DBH (75 cm), the height of the guide curve was 35.90 m. The height of the lower 95% confidence bound for an individual prediction at the reference DBH was 29.18 m. The difference between these two heights was 6.72 m. Therefore, using the guide curve method, and Equation (1) with an index DBH 75 cm, the site form equation was: SF ¼ 35:90 À ð35:90 À 29:18ÞP ( 3 ) The source data used to construct the site form curves is presented in Table 1. The resulting site form curves (constructed for DBH 10 to 120 cm and site form from 10 to 50 m) Figure 1. Observed total heights over DBH (a) and guide curve and bounds for 95% confidence intervals for the mean and individual prediction (b) related to oriental beech height-DBH.
for oriental beech are polymorphic in which the heights vary between different site form classes (Figure 2).
At the next step, site form for each plot was estimated referring to Table 1 and based on mean DBH and height and then classified into three categories by classInt package in R statistical software. Three site form categories were classified as poor ((30.45-30.55), medium (33.55-35.83), and good (35.83-39.30).
Environmental predictors for modelling site form are summarized in Table 2. Before running the ordinal regression, the variance inflation factor (VIF) was used to determine the collinearity between the independent variables (Hererra et al. 1999), where VIF < 10 indicates no collinearity problems. All the independent variables showing VIF > 10 were discarded for further analysis.
The results of ordinal regression are presented in Table 3. Of the 11 non-collinear variables entered into the model, phosphorous, clay, and altitude were the most important variables to predict the oriental beech site productivity at 0.05 significance level. Residual deviance, AIC and BIC for ordinal regression were 162.75, 188.75, and 223.25 respectively. The predictive performance of ordinal regression in terms of coefficient of determination (R 2 ) was 51.58%.
The means of all potential predictors were calculated for each ordinal categories of the response variable and plotted against Figure 3.

Discussion
This study presents a site form model for F. orientalis stands in Hyrcanian forest in the north of Iran. The difference equation method was used to develop dynamic site form curves. Although the Chapman-Richards model performed the best in the study area, but due to changes in site conditions and populations, the relationship between the diameter and height of a tree in different forest stands can also vary. Thus, relying only on the height-diameter model is not suitable for prediction of all relationships in a forest. Developing an individual height-diameter model for each stand, and/or using the generalized height-diameter model, are two alternative approaches to this problem. Based on the results, site form curves in this study created the polymorphic form. Durkaya and Durkaya (2007) developed the polymorphic site form curves for F. orientalis in Turkey. Duan et al. (2013) applied a two-parameter Korf equation for building a polymorphic site index equation for Chinese fir plantations.
Site form is commonly used as an indicator of forest productivity in which the mean height of dominant trees, defined as the height of the 100 largest diameter trees per hectare, is used. Site form for each sample plots was estimated using polymorphic height growth curves for oriental beech and then related to edaphic and physiographic variables using ordinal regression. Ecological variables used in many forest site productivity studies are edaphic and  topographic variables (Herrera et al. 1999;Herrera-Fern andez et al. 2004). Several authors have found correlations between these attributes and site productivity using multiple regression (Wang 1995), principal components analysis (Sanchez-Rodriguez et al. 2002), tree classification models (Verbyla and Fisher 1989), and generalized additive models (Aertsen et al. 2010). This study also showed that altitude, phosphorous, and clay are among the most important factors affecting the beech forest site productivity. Evaluation of site productivity based on environmental variables may help in forest management, with consideration of local growth conditions (Socha 2008). According to the results, clay has a negative effect on site productivity. Clay is one of the most important soil physical properties because it influences the ability of a soil to supply water, nutrients, and aeration necessary for plant growth. Means plots of site form against clay indicated that the high productive stands of oriental beech trees have between 30% and 35% clay in the soil. Herrera et al. (1999) reported a negative correlation between site productivity of Vochysia ferruginea in a secondary forest and soil clay content. Several studies mention the significant effect of clay on site productivity (Monserud et al. 1990;Corona et al. 1998). Because clay-textured soils have a higher percentage of small capillary pores, they do not release water as compared to loam-textured soils. Most plants are best adapted to extract available water from loam-textured soils as opposed to clay-or sand-textured soils. This may additionally hinder plants' ability to access oxygen in the soil. Thus, roots become stressed due to oxygen deprivation, reducing plant growth (Nepstad et al. 1994;Laurance et al. 2010). This study also showed that local site productivity is partly limited by phosphorus availability, which tends to be critically limiting to plant growth (Vitousek 2004;Lambers et al. 2008;Turner 2008). Means plots of site form against phosphorous showed that forest sites having good productivity occur in soils with high phosphorous levels and vice versa. This result contrasts with those reported by Herrera et al. (1999) who found a negative correlation (r = -0.59, P = 0.001) between phosphorous and the dominant height (as an indicator of site productivity) of Vochysia ferruginea in a secondary forest. Future forest growth studies focusing on the role of phosphorus on site productivity could thus play an important role in clarifying a crucial issue for sustainable forest production over wide forest territories.
Altitude is another important environmental variable affecting many aspects of forest ecosystems, manifesting in different ways such as changes in species composition, growth rates, productivity, form, and reproductive capacity. McGrath (1975) states that altitude, along with soil texture and other ecological parameters, should explain site quality. The most productive beech forests are located at altitudes between 1210 m a.s.l. and 1240 m a.s.l. The effects of altitude on forest site productivity were also found in previous studies conducted in boreal forests (Holmgren 1994;Klinka and Chen 2003;Fontaine et al. 2007) and in Mediterranean areas (Bravo-Oviedo and Montero 2005). All of these studies show that forest productivity decreases with increasing altitude. Klinka et al. (1996) found that elevation has a strongly negative correlation with forest productivity and temperature, and is considered as a limiting factor to forest growth in high altitude. Dahl (1986) hypothesized that respiration, or some factors closely associated with respiration, limits plant developments and growths in high altitudes.

Conclusion
There are several methods for estimating forest site productivity and constructing their productivity curves. The difference equation approach in the context of the Chapman-Richards model and the Edminster et al. (1991) method were applied in building site form curves for the current study. Since the Chapman-Richards model performed well in this study it can therefore be used to quantify the height-DBH relationship in various sites in Hyrcanian forests. Although the data were collected from a specific region, the site form curves constructed can be extended to various beech forest sites in the Hyrcanian forest to give a satisfactory estimation of potential productivity. Second, our results highlight the fact that the relationships between beech forest site productivity and the easily measurable topographic and edaphic factors may be applied in the estimation of potential productivity of beech forests in Hyrcanian forests of Iran. Three environmental factors (altitude, phosphorus, and clay) are important and affect the forest site productivity. Finally, this study provides more information about relationships between forest productivity and environmental factors, to help forest managers achieve sustainable management of beech forests in the north of Iran.