Environmental Drivers and Age Trends in Site Productivity for Oak in Southern Poland

: Site productivity provides critical information for forest management practices and is a fundamental measure in forestry. It is determined using site index (SI) models, which are developed using two primary groups of methods, namely, phytocentric (plant-based) or geocentric (earth-based). Geocentric methods allow for direct site growth modelling, in which the SI is predicted using multiple environmental indicators. However, changes in non-static site factors—particularly nitrogen deposition and rising CO 2 concentration—lead to an increase in site productivity, which may be visible as an age trend in the SI. In this study, we developed a geocentric SI model for oak. For the development of the SI model, we used data from 150 sample plots, representing a wide range of local topographic and site conditions. A generalized additive model was used to model site productivity. We found that the oak SI depended predominantly on physicochemical soil properties—mainly nitrogen, carbon, sand, and clay content. Additionally, the oak SI value was found to be slightly shaped by the topography, especially by altitude above sea level, and topographic position. We also detected a signiﬁcant relationship between the SI and the age of oak stands, indicating the long-term increasing site productivity for oak, most likely caused by nitrogen deposition and changes in climatic conditions. The developed geocentric site productivity model for oak explained 77.2% of the SI variation.


Introduction
Site productivity arises from the complex interactions between biotic and abiotic drivers in a particular location [1]. It provides key information for forest management and is the primary determinant of decisions concerning species composition selection, the rotation period and allowable cut estimation, the planning of silvicultural prescriptions, and forecasting the timber yield [2][3][4][5]. Under various management alternatives, site productivity provides a decision support tool and knowledge about species variability resulting from sensitivity to environmental factors. Thus, research on the site productivity is crucial for operational purposes and contributes to our understanding of site factors affecting the site productivity for given species.
The most commonly used indirect proxy for forest site productivity is the site index (SI) [1,[6][7][8], which can be estimated using two methods, namely, phytocentric (plantbased) or geocentric (earth-based) [1]. In the case of even-aged stands, the SI is most commonly estimated using phytocentric methods based on stand age and height [6]. However, for non-forest areas or stands of young age or uneven-aged with complex species composition, phytocentric methods are subject to significant uncertainty or even become impossible. In the case of such limitations, geocentric methods can provide a much better

Sample Plot Data
The research material collected in the years 2013-2014 consisted of data from 150 cir cular sample plots established in oak-dominated (Quercus sessilis and Quercus robur stands. Plot size varied between 0.02 and 0.12 ha, depending on the age of the stand, and was adjusted so that the number of trees on the plot was at least 20 ( Table 1). The numbe of trees per hectare ranged from 75 to 2350 (Table 1). Diameter at breast height (>7 cm) and height of all trees were measured in each sam ple plot. The stand age was estimated on the basis of increment core samples collected in each sample plot from the 5 trees with the largest diameter.

Sample Plot Data
The research material collected in the years 2013-2014 consisted of data from 150 circular sample plots established in oak-dominated (Quercus sessilis and Quercus robur) stands. Plot size varied between 0.02 and 0.12 ha, depending on the age of the stand, and was adjusted so that the number of trees on the plot was at least 20 ( Table 1). The number of trees per hectare ranged from 75 to 2350 (Table 1). Diameter at breast height (>7 cm) and height of all trees were measured in each sample plot. The stand age was estimated on the basis of increment core samples collected in each sample plot from the 5 trees with the largest diameter. The top height (TH) was defined as the mean height of 100 thickest trees per hectare. In order to eliminate the dependence of TH on the plot size, we used the adjusted largest trees (ALT) method, developed by García and Batho [39], in the calculations. According to the ALT method, TH for a sample plot of size A (ha) is calculated on the basis of 160 × A − 0.6 of the thickest trees. The SI was determined for the sample plots using the SI model developed for oak in Poland [40], on the basis of the dynamic formulation of the Korf equation [41], which describes changes in the top height with age Equation (1): where b 0 and b 1 are parameters of the model; H 0 is the height (in m) at age t 0 ; and t 1 is the base age, equal to 100 years.

Soil Data and Physiography
Soil samples were taken as a collective sample mixed of 5 sub-samples sampled within each plot from the mineral layer (0-40 cm). Prior to analysis, samples were sieved through a 2 mm mesh. The 2 mm sieved samples were used to measure the particle size composition (the percentage content of sand, silt, and clay), as determined using laser diffraction (Fritsch Analysette 22), and pH in water (pH H 2 O) and in 1M KCl solution (pH KCl), as determined potentiometrically. The sub-samples were ground to fine fraction using a ball mill (Fritsch) to improve their homogeneity. We measured the content of soil carbon (C) and nitrogen (N) with a LECO CNS True Mac analyzer in the fine sub-samples. The soil type was determined using a digital soil map for Poland with the scale of 1:500,000 [42].
The topographies of particular sample plots were described on the basis of indices derived from the European Digital Elevation Model (EU-DEM) [43]. On the basis of the EU-DEM, we calculated a set of derived topography variables: the hillshade, determining the degree of shading of the terrain [44]; landform, which provides information about categories of surface features [45]; and terrain convexity metrics, the Topographic Position Index (TPI), specifying the position of any point on the earth surface in relation to the neighboring points located at a given distance, and slope position. Hillshade was calculated according to the methodology proposed by Siwek and Wacławik [44]. The landform index was calculated on the basis of TPI values calculated for the 2 neighbor scales of 900 and 1500 m [45]. TPI was computed on the basis of a digital terrain model using the ArcGIS 10.3 software. The slope position was calculated using the TPI and its standard deviation [45].
The geological substratum type was determined using a geological map of Poland with a scale of 1:500,000 [46]. The basic characteristics of the environmental variables are summarized in Table 2.

Model Development
Site productivity may be modelled through a wide range of modelling approaches; however, we proceed with different types of variables (quantitative and qualitative). Therefore, collinearity of predictor variables may occur and problems with fulfilling criteria for the normal distribution of residual values of analyzed variables. Aertsen et al. [47] showed that generalized additive models (GAM) are one of the best methods for geocentric site productivity modelling, which offers good model performance and the best ecological interpretability. Type of geology substrate 1-Quaternary aeolian sand; 2-Quaternary sand and residual gravel with boulders; 3-Tertiary Neogene Krakowieckie clay; 4-Quaternary South Polish Pleistocene sand, silt, and gravel; 5-Quaternary Pleistocene glacial till; 6-Quaternary Holocene fluvial deposits; 7-Quaternary dune aeolian sand; 8-Quaternary Pleistocene sand and gravel with boulders; 9-Quaternary Holocene preglacial sand and gravel; 10-Quaternary Pleistocene glacial till residuum; 11-Quaternary South Poland Pleistocene glacial till and its residuum; 12-Quaternary Pleistocene loess; 13-Cambrian mudstone, claystone, and sandstone; 14-Devonian limestone; 15-Quaternary Pleistocene sand, glaciofluvial, and fluvial gravel; 16-Quaternary Pleistocene loess with fossil soil and Baltic glaciation; 17-Cambrian claystone, clay-slate, mudstone, sandstone, and quartzitic sandstone; 18-Quaternary Pleistocene periglacial deposits; 19-Quaternary Holocene fluvial deposits in general.
GAM models [48] enable the estimation of multivariate variables using the additive approximation of the regression function by substituting the explanatory variable linear function with a non-parametric function; this can be estimated using a smoothing spline function (spline). Application of a GAM method allows us to explain the effects of independent variables by analyzing the model graphs' partial residual spline curves, depending on the independent variable values. This approach allows for establishing a direct relationship of the SI with individual variables while holding all other variables fixed [49]. This property of GAM allows for analyzing the specific patterns of site productivity in response to environmental gradients. The GAM method allows for interpretability, flexibility and automation, and regularization [49]. Mainly, when the model contains nonlinear effects, GAM provides a regularized and interpretable solution, whereas other methods generally lack at least one of these three features. GAM strikes a nice balance between the interpretable, yet biased, linear model and the extremely flexible "black box" learning algorithms [49].
Considering that GAM models are not fully resistant to predictors collinearity, we applied the correlation matrix to assess the presence of collinearity for quantitative variables with a cut-off value of 0.9. The variables for which collinearity was found were eliminated from the set of independent variables. The procedure was performed using the R package caret [50].
An established set of variables was used for the developed final GAM models according to Equation (2) in order to determine the relationship between the SI and the multiple variables characterizing tree growth conditions for oak stands.
where G(E(Y)) is a function that links the expected value of the variable Y E(Y)) with additive functions of the explanatory variables For examination, the model performance, and possible overfitting in calculating adjusted R 2 , we also used 10-fold cross-validation. The procedure was performed using the R language package gam and package caret [50].
Next, we examined the effect of age on SI in order to detect an age trend. For this purpose, besides soil properties, geology, and topography, the age of oak stands was accounted for in the model. Using the GAM method, we were able to examine the effect of age of stands on the SI whilst holding all site variables fixed. Therefore, this approach provides the ability to remove the effects of other variables that affect site productivity, which may lead to over-representation of the oldest stands on infertile sites and more productive young stands [51,52].
In the last step, we assessed the GAM model performance; for this purpose, the following measures were used: where y i is the observed value andŷ i is the predicted value.
• adjusted coefficient of determination: where y i are the observed values,ŷ is the model value, p denotes the number of parameters used in the model, and R 2 is the coefficient of determination. We also assessed the adequacy of the SI model fitted to the whole data set by analyzing the plot of the model residuals and the histogram of the error distribution.

Results
For the analyzed sample plots, the SI ranged between 21.09 and 37.77 ( Table 2). As a result of modelling, it was found that most of the variables selected for the analysis had a significant impact on the SI. Only TPI900, pH H 2 O, silt [%], and C/N were the variables that were not significant or showed collinearity; therefore, they were not included in the model. Analyzing individual environmental variables, we found that the SI of oak stands was highly determined by the percentage of nitrogen and carbon content in the soil. Nitrogen content was positively correlated with the site productivity for oak. Increasing the nitrogen content in the range from about 0.01 to 0.08 caused a rapid increase in the SI; however, after reaching 0.10, further increasing the nitrogen content did not increase the SI. A high level of carbon content in soil indicated lower site productivity for oak ( Figure 2). The relatively fast decrease in SI was observable, with increasing carbon content up to 1%; however, after reaching this level, the SI decrease slowed down significantly. Analysis of the spline curves of the model partial residuals for nitrogen and carbon content in soil demonstrated that the variation of SI was in the range of 5 m; however, in the case of soils with higher nitrogen and carbon contents, considerable uncertainty was visible ( Figure 2). Our analysis also indicates the influence of soil physical properties on the SI. We found an alteration in site productivity observable with sand and clay content shifts in the soil ( Figure 3). Sand content in soil was negatively correlated with the SI. These results indicate a reduction of site productivity, especially for soils with sand content above 50%. The percentage content of the clay fraction in the soils had an opposite effect on the site productivity for oak. The positive effect of the percentage of clay content in soil was pronounced until the level of 40% was reached; after that, a further increase of the clay fraction in the soil did not significantly improve the site productivity for oak. Besides the soil properties, the obtained results showed a significant influence of soil type on the site productivity variability. The highest site indices were found in the stands growing on gley-podzol mucky soils, while the lowest valuations were observed on gleypodzol soils appropriate for coniferous sites (Figure 4). Other soil types showed less variation of SI, with differences not exceeding 1 m, compared with the brown rusty soil benchmark in the GAM model. Our analysis also indicates the influence of soil physical properties on the SI. We found an alteration in site productivity observable with sand and clay content shifts in the soil ( Figure 3). Sand content in soil was negatively correlated with the SI. These results indicate a reduction of site productivity, especially for soils with sand content above 50%. The percentage content of the clay fraction in the soils had an opposite effect on the site productivity for oak. The positive effect of the percentage of clay content in soil was pronounced until the level of 40% was reached; after that, a further increase of the clay fraction in the soil did not significantly improve the site productivity for oak. Our analysis also indicates the influence of soil physical properties on the SI. We found an alteration in site productivity observable with sand and clay content shifts in the soil (Figure 3). Sand content in soil was negatively correlated with the SI. These results indicate a reduction of site productivity, especially for soils with sand content above 50%. The percentage content of the clay fraction in the soils had an opposite effect on the site productivity for oak. The positive effect of the percentage of clay content in soil was pronounced until the level of 40% was reached; after that, a further increase of the clay fraction in the soil did not significantly improve the site productivity for oak. Besides the soil properties, the obtained results showed a significant influence of soil type on the site productivity variability. The highest site indices were found in the stands growing on gley-podzol mucky soils, while the lowest valuations were observed on gleypodzol soils appropriate for coniferous sites (Figure 4). Other soil types showed less variation of SI, with differences not exceeding 1 m, compared with the brown rusty soil benchmark in the GAM model. Besides the soil properties, the obtained results showed a significant influence of soil type on the site productivity variability. The highest site indices were found in the stands growing on gley-podzol mucky soils, while the lowest valuations were observed on gley-podzol soils appropriate for coniferous sites (Figure 4). Other soil types showed less variation of SI, with differences not exceeding 1 m, compared with the brown rusty soil benchmark in the GAM model.  We found high differentiation in site productivity, depending on the geological substratum type ( Figure 5). The highest mean values of site productivity were obtained for Quaternary Holocene fluvial deposits (6), forming fertile soil formations. The lowest potential productivity can be expected in the Quaternary South Poland Pleistocene sands, silt, and gravel (4); Quaternary Pleistocene sand and gravel with boulders (8); as well as Quaternary southern Poland Pleistocene tills and their residues (11) and Quaternary Pleistocene loess (12), which are associated with sandy soil formations. Other types of geology substrata showed a smaller variability of site productivity. We found high differentiation in site productivity, depending on the geological substratum type ( Figure 5). The highest mean values of site productivity were obtained for Quaternary Holocene fluvial deposits (6), forming fertile soil formations. The lowest potential productivity can be expected in the Quaternary South Poland Pleistocene sands, silt, and gravel (4); Quaternary Pleistocene sand and gravel with boulders (8); as well as Quaternary southern Poland Pleistocene tills and their residues (11) and Quaternary Pleistocene loess (12), which are associated with sandy soil formations. Other types of geology substrata showed a smaller variability of site productivity.
The site productivity for oak was also slightly affected by elevation a.s.l. and the TPI1500 index. The optimum range of site productivity was found to be between 225 and 275 m a.s.l., with further elevation increases associated with a gradual decrease in the SI ( Figure 6). The relationship between the residuals of the geocentric GAM model and TPI was relatively low. However, some regularity was identified. Higher values of the TPI were related to a slightly higher average SI, as determined on the basis of the spline curve ( Figure 6). stratum type ( Figure 5). The highest mean values of site productivity were obtained for Quaternary Holocene fluvial deposits (6), forming fertile soil formations. The lowest potential productivity can be expected in the Quaternary South Poland Pleistocene sands, silt, and gravel (4); Quaternary Pleistocene sand and gravel with boulders (8); as well as Quaternary southern Poland Pleistocene tills and their residues (11) and Quaternary Pleistocene loess (12), which are associated with sandy soil formations. Other types of geology substrata showed a smaller variability of site productivity.  The site productivity for oak was also slightly affected by elevation a.s.l. and the TPI1500 index. The optimum range of site productivity was found to be between 225 and 275 m a.s.l., with further elevation increases associated with a gradual decrease in the SI ( Figure 6). The relationship between the residuals of the geocentric GAM model and TPI was relatively low. However, some regularity was identified. Higher values of the TPI were related to a slightly higher average SI, as determined on the basis of the spline curve ( Figure 6). The results of GAM modelling showed a strong age trend in the SI for oak (Figure 7). From the age of 60, we found that SI demonstrated a sharp and successive decrease, reaching 2 m in every subsequent age class. Therefore, a 140-year-old stand growing in the same site conditions obtained a SI reduced by 6 m, compared with 60-year-old (youngest age) stands. The results of GAM modelling showed a strong age trend in the SI for oak (Figure 7). From the age of 60, we found that SI demonstrated a sharp and successive decrease, reaching 2 m in every subsequent age class. Therefore, a 140-year-old stand growing in the same site conditions obtained a SI reduced by 6 m, compared with 60-year-old (youngest age) stands.

Discussion
Our results demonstrate that the developed geocentric model explained a significant part of the variability of the site productivity for oak. Among all investigated site variables, a significant impact on SI was observed for predictors concerning soil properties. The nitrogen content had the most substantial effect on the site productivity for oak. In investigated stands, an increase in the nitrogen content in soils was associated with a considerable rise in the SI to a certain level. Similar results have been stated in studies involving pine stands, in which it was found that, after reaching a certain threshold, a further increase in nitrogen concentration does not result in a proportional increase of the SI [53]. The effect of nitrogen content in soil on the site productivity is biologically plausible and compliant with oak ecological requirements. Due to the high site requirements, this species prefers soils rich in nutrients and clay; therefore, the highest site productivity was observed with soils that represented such conditions, especially those with the high nitrogen content. The relationship between SI and N content indicates the possibility of improving the site productivity for oak resulting from the increasing nitrogen deposition. The opposite relationship between SI and carbon content in soil ( Figure 2) turned out to be a good indicator of the site productivity for oak. A high level of carbon content is characteristic of low fertility soils with low biological activity and a slow rate of soil processes; therefore, a high carbon content indicates that the soil is less suitable for oak growth. Sites with high carbon content were unfavorable for oak, which was reflected in its correlation with SI.
Our results showed shifts in site productivity along with the physical properties of soil and soil type. The types of soils characterized by high clay content, moderate moisture, and pH close to neutral demonstrated the highest site productivity ( Figure 5). The lowest productivity of oak sites was observed with gley-podzol soils-typical of coniferous forest sites-with high sand content. The optimal substrate for oak is fertile, rich, and moderately moist with a neutral or close to alkaline pH soil. Therefore, sites with such properties showed the highest SI values. Besides soil properties, we found that oak stand productivity was strongly affected by the type of geological substrate. In comparison with Quaternary aeolian sands-used as a reference level for modelling-the highest SI was obtained for the stands growing on Quaternary Holocene fluvial deposits, which are the most fertile. Fine-grained and loamy material contains a large amount of the smallest soil fractions rich in macro-and micronutrients [54] and, therefore, is associated with higher productivity. Additionally, research involving indigenous oak stands [55] indicated that the penetration of roots of both oak species through gravelly and loamy soils was hampered while riverside meadows were hampered where clay soils prevail-oak roots could easily reach the groundwater table. The type of geological substrate accounts for a large proportion of the SI variation and, along with soil sub-type, it should be taken into account when evaluating the potential productivity of oak stands.
Our results also suggest that the topography affects the site productivity of oak. Topography, expressed as altitude and TPI, can be used as an indirect measure of climatic condition variability [56]. The gradual decrease of SI from 260 m a.s.l. is likely caused by interactions between factors affecting moisture and temperature that can be observed at different elevations above the sea level. The research conducted by Lyr [57] indicated that the optimal temperature for the growth of aerial parts of common oak is 25 • C, while, in the period from April to September, this tree species requires average temperatures reaching at least 12 • C. The increase of altitude up to 250 m is connected with the limited occurrence of frost hollows and positively correlates with the total annual precipitation, which is the most likely reason for exerting a positive effect on the oak growth [52]. Additionally, higher SI values were obtained for higher TPI values, corresponding to the location above the surrounding area. It can be assumed that this resulted from the effect of topography on moisture and thermal conditions. Higher locations, however, are associated with worse thermal conditions, which is the likely cause of the decline in site productivity for stands located above 250 m above sea level.
Interestingly, the oak SI displayed a strong negative correlation with stand age. The results presented in our studies are almost tantamount to observations from central Europe. Changes in the SI reported for oak stands in southwestern Germany [52] indicated the average increase in the SI amounted to 2.1 m per each 20-year period. Similar age trends have also been previously detected and usually interpreted as an effect of the changes in site conditions [24,27,28,30,52] associated with the effects of non-static factors such as nitrogen deposition and rising CO 2 concentration. However, highly productive forests are usually managed under short rotation. Therefore, old stands are often found on poor sites [30,51,52], which may also result in the negative correlation of SI with stand age. In the present research, this problem was solved by the use of GAM models. The advantage of GAM is that the influence of each explanatory variable can be examined individually, while all other variables are constant [49]. Using this approach, we examined the effect of germination age on SI, keeping all site variables fixed. This allowed us to remove the effect of site variables and individually assess age influence on the SI.
The occurrence of an age trend means that the stand density increases in a shorter time and trees more quickly reaching specific dimensions [58]. Therefore, its detection can be a powerful tool to develop forest management guidelines concerned with planning the rotation age and silvicultural treatments. Furthermore, understanding this phenomenon effect may be crucial in predicting climate and environmental changes and the impact of non-static factors on forest productivity. We found the age trend in SI, indicating an increase of site productivity for oak. SI is strongly affected by nitrogen content in the soil. However, after reaching a specific level, a further increase in nitrogen content did not increase SI (Figure 3). Therefore, it is probable that the observed increase in the productivity of oak stands will stabilize at a certain level. This phenomenon may be of great importance for future carbon sequestration in oak stands, and thus it is a call for further research.
The results of the geocentric modelling of site productivity demonstrated in this article are crucial to precisely define the conditions under which oak stands will achieve the highest productivity in Central Europe. Information concerning trends in site productivity may provide essential information for management, as well as adaptation and mitigation strategies against changing climatic conditions [23]. There have been many site productivity simulations under the changing climate conditions of European forests [5,59,60]. The observed trends of increasing the average temperature and reducing rainfall during the growing season adversely affect North American species (e.g., lodgepole and white spruce) [61]. The range of distribution of other species, such as oak (Quercus L.)-which are more thermophilic and relatively resistant to drought or strong winds-may increase their economic importance [36,52]. Our study points out the drivers of site productivity for oak, as related to terrain and soil properties. The relationships between site productivity and features demonstrated in this study confirm this possibility of directly applying selected environmental characteristics for assessing the potential productivity of oak at the local scale. The limitation of this study is related to the relatively small scale of the research area. Data of this type do not allow for extrapolation of the results to large areas. However, research limited to a smaller area also has some advantages, as it allows for the exclusion of the influence of many factors that differ at a large scale, such as, for example, climatic conditions [37]. Because we excluded the influence of variation of climate conditions, it was possible to describe the influence of local drivers of site productivity, mainly nitrogen content, and microclimate represented by the topography. The results enrich the knowledge about the autecology of oak; however, they cannot be uncritically related to the entire area of its occurrence. Our results also highlight the role of non-static environmental factors in affecting forest site productivity, as expressed by an age trend. The strong age trend in SI confirms increasing site productivity for oak for several decades in the research area. This calls for further research in order to assess the relationships between SI and environmental factors and age on a larger spatial scale and for other tree species.

Conclusions
We developed a geocentric model of site productivity, providing an appropriate explanation of site variable impact on the site productivity for oak. The oak SI value depends predominantly on the physicochemical soil properties-mainly nitrogen, carbon, sand, and clay content. The oak SI is also slightly shaped by the topography, especially altitude (a.s.l.). We also detected a significant relationship between SI and age, indicating the long-term increasing site productivity for oak, which is most likely caused by nitrogen deposition and climate conditions changes.