Stand density biases the estimation of the site index especially on dry sites

The site index (SI) has been widely used in forest management and silviculture. It relies on the assumption that the height of dominant trees in a stand is independent from the local density. However, research on climate change suggests that under certain moisture stress conditions, this may not hold. Here, based on 29 plots from five long-term research experiments, we tested the effect of local stand density on the SI of Norway spruce (Picea abies (L.) H. Karst). With generalized additive models (GAMM), we analyzed the effect of stand structure and climate predictors on SI. The two evaluated models revealed that local stand density and age had a significant effect on SI (p ≤ 0.001), showing a clear negative trend especially significant on sites with poor and dry soils, which may reduce the SI by a maximum of approximately 4 m for an increase in density of between 400 and 600 trees/ha. We stress that the physiological characteristics of Norway spruce, flat-rooting system and xeromorphism, especially when growing in pure stands, may explain these effects. Thus, density control and growth in mixtures may help to reduce the water stress and losses in height growth under future climate conditions.


Introduction
Measuring the quality of a site has been a major concern of forestry science since its beginning.
The height of dominant trees at a defined reference age, in even-aged stands, known as the site index, has been widely used as a measure of site quality (MacFarlane et al., 2000). First experiences in forest science have assumed that the site-index equations, the dominant height versus age relationship, remain unchanged by stand density (Huang and Titus, 1993). In the 60s, Vincent (1961) recognized in his comprehensive review that the use of total height over age as an index of site quality is not as simple as some textbooks on forest mensuration might lead one to believe. Besides, many things may happen to a tree which can alter the height-age relationship so as to obscure any differences ascribed to site quality (Puhlick et al., 2013;Vincent, 1961).
While it is recognized that stand density affects average stand height, the evidence on the dominant height is yet unclear (Burkhart and Tomé, 2016). Many authors have discussed the influence of stand density on the height of the dominant trees, finding both no effects and negative effects depending on the experiment and species carried out in the respective studies. Studies like MacFarlane et al. (2000) or Antón-Fernández et al. (2011) maintained that assessments of site index are confounded by the influence of population density on height growth, even when only the tallest trees in the stand are used to obtain the index.
Due to the general lack of long-term studies with a broad density variation in most regions, Nelder plots (Nelder, 1962) have been commonly used to study tree growth under different densities.
However, studies are frequently based on the very early growth stages. For example, in Kuehne et al. (2013), it was observed that the total tree height in Nelder plots from oak is lower at higher densities, especially due to the presence of smaller suppressed trees. In Woodruff et al. (2002), for a plantation density experiment of Douglas-fir (Pseudotsuga menziesii (Mirbel) Franco), it was observed how high initial planting density has a positive effect on height growth. However, Newton (2015) observed a negative relationship between density and height, in Nelder plots of Jack pine and black spruce, which tends to be asymptotic in nature. This means that the differences in height growth declined with age, and initial density has a different impact on growth compared to later development phases. In an experiment on a Nelder plot in Scots pine (Pinus sylvestris L.), which was first reported in Dippel (1982), and then repeated in Spellmann and Nagel (1992), it was observed that in advanced stages of development, the heights in the middle of the plot, where the density is the highest, were lower than in the less dense parts. This indicated a possible effect of higher density, and showed the effects on a later stage of development, compared with other studies.
In Nelder studies, during the earliest measurement periods, the relationship of height increment to tree density may be a result of the extremely high densities in the plots (Hummel, 2000). Thus, the impact of high densities on height growth in young trees is demonstrated through such experiments (Woodruff et al., 2002). Also studies focused on the dominant height, like MacFarlane et al. (2000), which was extended by Antón-Fernández (2011), questioned the density-independence assumption of the dominant height in even aged stands of loblolly pine (Pinus taeda L.), stating that high densities may have a negative impact on dominant height. Antón-Fernández (2011) even stated that forest managers should adapt the site index for different planting densities, especially if high densities are involved. The general consensus states that these effects are compensated during stand development (micro-site occupation, best individuals grow faster, the others die).
However, the question whether stress factors on poor sites that affect the rooting system, such as drought, may have an impact on growth at high densities remains unclear.
Density can be controlled through initial escapement and thinning operations. In studies of stress by drought, it can be observed how thinning mitigates growth reduction, i.e. increases resistance, during a drought year, and enhances recovery of stem growth during following years (Sohn et al., 2016). It has also been proved by Sohn et al. (2013), D'Amato et al. (2013 and Gebhardt et al. (2014) that a reduction in density can enhance the capacity of spruce stands to cope with lack of water availability. However, this mitigating effect depends on the time span between thinning and drought: the extra water availability decreases as trees grow and increase their leaf area and fine root biomass. Another indicator that points out towards impacts of density on the height-age relationship is the hydraulic limits to the height growth (Ryan and Yoder, 1997), which may play an even more relevant role under conditions of water stress.
An indication that the reduction of height growth with higher competition for water in denser stands may be especially relevant on poor sites is explained as follows. As trees compete mainly for belowground resources, and density especially increases competition for them, height growth may be reduced. This means, a prioritization of the height growth does not bring any competitive advantage und thus, both, diameter and height growth may decrease monotonically (Pretzsch, 2019). This can be especially relevant in monocultures of Norway spruce (Picea abies (L.) H. Karst), whose physiological characteristics (flat rooting system and isohydric) make this species particularly sensitive to water stress conditions (Lyr et al., 1992;Puhe, 2003).
With increasing drought stress such negative effects of stand density and water limitation on height growth (Briffa et al., 2009;Spiecker, 2003) and bias in the estimation of the site index are prone to increase. Moreover, in addition to stand density and ageing, climatic factors are being reported to affect growth partitioning, with species-specific differences (Trouvé et al., 2017). Thus, it can be expected that, especially on sites under water stress, the effects of stand density on height growth may introduce substantial variances on the estimation of the site index, also motivated from different allocation strategies (Uhl et al., 2015).
In order to assess these effects, in this study, we used statistical modelling based on the long-term Bavarian experimental plots. Inventory data from over 40 years in Norway spruce stands, distributed along the southern region of Bavaria, will offer new insights to address the following specific questions: Q1-Is the estimation of the site index biased by a reduction of height growth due to local density? Q2-Does water availability or drought modulate the effect of the local density on the estimation of the site index?

Data
In this study, we have used data from the long-term forest experimental plots network from Bavaria. We have chosen for this purpose, monocultures of Norway spruce spread throughout the state. These correspond to a diverse range of thinning and planting experiments (see Table B1).
For each research experiment, different kinds of thinning (from above, below and selective) and strength of thinning (unthinned, slight, moderate to strong density reduction) have been tested in independent plots, which lie in close vicinity. For this reason, a wide range of densities was available, from sparse densities, where trees are growing under low competition, to dense unthinned plots. In the following, research experiment will denote the specific experiment per site, and research plot the individual thinning or plantation density plots within the experiment.
The inventories were repeated every 10 years. Some of the plots have been measured since 1880.
However, as the climatic data is only available since 1975 we have used inventory data only since this year. The spatial position of all the trees within the plot is known, together with the diameter at breast height (dbh), which is measured for every tree. Stand age was calculated based on the time passed since the establishment of the experiment. Since all experiments are based on evenaged stands, stand age can be assumed equal for all tress in each plot.
During the standard inventory of the research experiments, tree height was only measured for a subset of trees. Such measurements are designed to cover a representative spectrum of heights and diameters. In order to have height estimates for all trees, we fitted a log-function = + ln ( ℎ) (where denotes the tree height, ℎ the diameter at breast height and a and b are empirical parameters) to the measured height-diameter relationship for each year, and applied it to the rest of the individuals. This was performed at the research (thinning) plot level, in order to have enough measured individuals to ensure a representative curve for each experiment and silviculture regime. A summary of the main characteristics of the stands and sites is provided in Tables 1-3, and detailed information regarding the number of plots and subplots (Table B2) and the specific structure characteristics depending on the thinning experiment (Table B1) are included in Appendix B.
We have then grouped the experiments in two classes of water availability, we used the information provided by the forest administration data, according to the type of substratum, depth of the soil and overall water availably. We named these groups "moist" and "dry".  Table 2). The dry group is characterized by low water availability, low precipitation and, especially by low quality soils. Vohenstrauß 622 presents only 375 mm of precipitation during the vegetation period, and it is additionally located on sandy soils.
Weißenburg 613 has a better substratum but also not optimum soil, being very shallow and with a low water availability. The rest of the sites present optimum soil conditions, over moraines and with high precipitation, indicating optimum growth conditions for Norway spruce. A special case was Eurach 606. This site is a desiccated swamp over moraines. The exact soil site is unknown, but records of the forest enterprise suggest good water availability.
In order to track differences in local density, and thanks to the known location of each tree, we have virtually extracted sub-plots with 10 m radii, overlaid on the actual research (thinning) plot ( Figure 1). First, we decided the radius to be constant instead of height dependent, in order to have the same probability of having a dominant tree in each circle and across the research plots. Second, we considered that a 10 m radius, results in enough replications per research plot, which typically have dimensions between 500-1000 m 2 . Thus, as a kind of moving window, we have moved the circle in 5 m steps in the x and the y directions to cover the entire research plot. Within the mentioned circle, we have calculated at every step the following (local) stand structural variables: age (years), number of trees per hectare ( ) in tree/ha, quadratic mean diameter ( ) in cm, dominant height ( ) in m, calculated as the mean height of 100 thickest trees per ha, and the stand density index ( ), which was calculated according to Reineke (1933), but using the allometric exponent (-1.664) estimated for Norway spruce by Pretzsch and Biber (2005): (1)

9
Finally, the site index ( ) by Assmann and Franz was calculated based on the dominant height at age 100 according to yield tables for spruce in Bavaria (Assmann and Franz, 1972). The was calculated dynamically for each year, so it reflects changes in site quality with time. The suitability of this can be observed in Figure 2 where we have plotted the development of the dominant height ( 100) for the research experiments under study, compared with the curves by Assmann and Franz (1972). These SI curves were fitted with data that goes back to the 1880s, being the most appropriate curves to track forest development before the climate change effects became substantial.
Due to the size of the circular sub-plots only three trees per subplot are selected for the estimation of the dominant height, which can lead to underestimations of the 100 and as consequence in the (García, 1998). Even if for the rather homogenous vertical structure of the research experiments studied here such biases could be negligible ( 100 ≅ ), the effect must be taken into account when discussing the results.
Two climatic variables were additionally estimated: mean precipitation of the vegetation period (TotPrepVP) in mm and the mean temperature of the vegetation period (MeanTempVP) in °C.
The vegetation period was defined as the months with a mean temperature above 10°C. These data were estimated using the Gridded Agro-Meteorological Data from the European Commission, Joint Research Center -JRC with 25x25 km resolution grid (CGMS, 2014), and assigned to each experimental plot based on their geographical coordinates. The data were available from 1975 until the present time.

Statistical analysis
We have chosen generalized additive mixed models (GAMM) (mgcv package) (Wood, 2017; R Core Team, 2018) as our main statistical tool to examine the impact of stand density on the site index. GAMs offer a convenient way to incorporate non-linear effects (Biber et al., 2013). Due to spatial and temporal correlation effects between the observation units, we adopted a nonlinear mixed-effects model with random effects potentially at each of the hierarchical plot units (moving sub plot, plot and research experiment) and/or the year. To scrutinize which levels are optimum we carried out an Akaike's Information Criterion (AIC) (Akaike, 1974) analysis for all levels of nesting showing that including the research (thinning) plots and the "moving window circle" subplots would not improve the estimation. Therefore, we accounted for the autocorrelation effect by introducing a random effect only at the research experiment level (RE) and the year. The approach is logical given that the data from the research experiments, included in this study, have been collected for over 40 years, so main random effects are assumed to be reflective of the different measuring teams and climate effects. The random effect at the research experiment can additionally account for potential site effects and different population clusters.
The predictor function for a GAMM has the following form: , ⋯ , , , ⋯ , where ⋯ , is a set of explanatory variables, ⋯ is a set of explanatory variables, ⋯ is a set of nonparametric smoother functions, are regression parameters, the research experiment and year specific random effect ( = , = ), and the error term. The indices and denote the -th observation year on the -th research experiment. In this work, the dependent variable will be SI and the explanatory variables the structural and climatic variables we have introduced above, and their potential interactions. These are: , , , , TotPrepVP, MeanTempVP and the type of soil as a categorical variable. As , and are variables explaining local stand density they may be co-dependent. Thus, we analyzed possible collinearity effects calculating the correlation between each pair. For pairs with a correlation coefficient | | > 0.7 one of the predictors was excluded from the model (Dormann et al., 2013).
Only significant variables ( < 0.05) were included in the final model specification as predictors; non-significant variables were excluded during the model runs. Yet, GAMs are best interpreted by visual examination instead of the statistical significances, as these are less expressive in GAMs than in other model types (Biber et al., 2013). Variable importance was assessed based on the relative proportion of variance explained by each predictor. In order to quantify the explained variance for some relevant variables, we run AIC tests between models, which were identical, except for including or not the target variable.

Test of collinearity and variable selection
After the collinearity test among the density predictors, we observed potential non-linear effects between − and − so we decided to estimate the correlation between the logtransformed pairs. The correlation coefficients between , and and their transformations were lower than | | = 0.7, < 0.05, except for ln( ) and ln( ) (| | = 0.75, p = 0.00).
Moreover, as we decided to include the interaction between and , we also tested the correlation between the density variables and tree age. In this case, a correlation coefficient of | | = 0.90 (p = 0.00) between and age was observed. For this reason, we decided to exclude during model selection. Correlations plots with the correlation coefficients and p-values for the mentioned pairs can been seen in Figure 3.
Preliminary model selection was based on biological plausibility including the structure variables with the climatic variables TotPrepVP and MeanTempVP, and differentiating between the effects that may be a consequence of the soil water availability. These were models that included the density indicators, linearly or as smoother terms, and the interactions between pairs. Specifically, we compared the relative importance of and the interaction ( , In model 1, all predictors included as smoother terms were grouped by type of soil (dry or moist).
The coefficient of determination for this model resulted in = 0.31 and the abovementioned = 7521.481. The parameter estimates (Eq. 3) are summarized in Table 4. With a closer look into Figure 4, we can observe with a higher detail the non-linear predictors.
In Figure 3, we observed that the correlation coefficient between ln( ) and ln( ) was low (| | = 0.43) but reflected a non-linear effect due to collinearity, and it showed some clustering effect that can be caused by the site effects. A GAM model can account for these non-linear effects; however, we decided to test a second model that includes a smoother term with the interaction between SDI and N resulting in model 2 (Eq. 4). Model 2 had higher coefficient of determination ( = 0.46) and also a lower AIC (6901.675). Is the estimation of the site index biased by a reduction of height growth due to local density? (Q1) The reduction of height growth was evaluated looking especially at the effect that the interaction between and had on the site index. In model 1, we first observed a negative effect in the interaction between and , indicating that density decreases the value of for a constant age.
Model 2, showed a higher biological plausibility that is worth discussing. Most importantly, the slope in the interaction between and , which also showed a negative effect, decreased (see Even if the effects of both TotPrepVP and MeanTempVP showed a linear effect, we decided to keep them in the models as smoother terms in order to assess the differences observed between the soil moisture types (the distinction by type with the R mgcv package is only possible when the predictors are included as smoother terms (Wood S, 2007)). Tested in model 1, if both terms were considered linear and independent from the soil type. TotPrepVP and MeanTempVP, with a negative and positive slope, appeared to be non-significant ( = 0.1946) and highly significant ( = 0.0101), respectively. However, when we included both terms depending on the soil moisture type, we observed that then MeanTempVP appeared significant ( = 0.0053) but only for the moist soil type, while TotPrepVP was significant ( = 0.0015), but only for the dry soil type, indicating that only when water is limiting in the soil, lower precipitation reduces . For the dry sites, the effect has been negative below 560 mm, and positive above this level. Accordingly, the opposite effect was observed for MeanTempVP on the moist sites: above 15 degrees, higher temperatures have a negative effect on the SI while the effect is linearly positive below this level.
In model 2 the effects of the weather variables (MeanTempVP and TotPrepVP) remained unchanged with similar levels of statistical significance (Table 5).
With the selected models we performed model predictions,

Discussion
Explanation of the findings Site index and growth can be misjudged when using common systems when not acknowledging the potential effect of stand density. In this study, the effect of local stand density in the estimation of was found significant for spruce monocultures in Bavaria, southern Germany (Q1 densities. This is explained in the theory of yield levels and can be expressed using the (Franz 1967, Bergel 1989). Thus, here, we hypothesized that sites with higher availability of water are able to support a higher number of individuals. For this reason, such effect can be potentially explained by the , included in the models. Moreover, thanks to mixed models, the codependence and thus non-linear effects between the interaction terms ( , ) and ( , ) were accounted for, and thus, being able to answer our first research question (Q1). This manifested the different effects that has on (and as a conclusion on height growth), but also, depending on the climate and water availability in the soil (Q2).
Specifically, the models we tested provided consistent results that supported the assumptions made prior to this study. First, the interaction between and have remained robust independently from the predictors chosen in the model. Particularly interesting was the difference in the effects of this interaction between groups of sites with a lower and higher water availability in the soil between the two models. When looking at model 2 ( Figure 5) we could observe that SI decreases almost linearly when increasing . Such attenuation effect is much stronger for the moist soils than for the dry ones. In this case, in moist soils after a density ~ 1000 , is almost not affected for an increase in N for the same age. In model 2 this attenuation is even stronger, indicating a higher biological plausibility (supported by a low AIC) and more clearly distinguishing the hypothesized differences in the SI response to density between moist and dry soils (Q2).
With the model predictions in Figure 6 and Figure 7 the dominant influence of the interaction between and , could be clearly observed, as expected from the variance importance assessment. Such predictions manifested the differences between the potential effects of local density depending on the quality of the site, in terms of water availability, and this depending on the precipitation. On sites with poorer water supply the effect increased up to a maximum reduction in the order of ~4 m dominant height for an increase of ~300 trees/ha . As abovementioned, we could observe, in both models, how such maximum potential density biases in the estimation of SI would occur at the lower densities, that is, stands with a greater number of large trees. The observed linear and negative relationship between TotPrepVP and (in sites with lower water availability), may also indicate that, under future trends of reductions in precipitation (especially during the vegetation period) and drought events, reductions of height growth are expected to be stronger. At the same time, under climate change conditions, the reduction of precipitation can be strong enough to limit water availability even on sites where currently the water supply in the soil is not limiting tree growth.
As pointed out by García (1998) and later by Ritchie et. al (2012), the smaller plot size (respect to the plot size used for development of the SI curves) here selected, may incur an underestimation of the top height ( 100) and therefore the . However, this effect would not affect the conclusions we derived here. Especially in the case of few large trees, an overall underestimation of the dominant height would mean that the impact that density has on the underestimation of SI, could even be larger than estimated. On the one hand, and as abovementioned, this bias is expected to be low in homogenous pure stands, and as such, it is systematic over all plots. On the other hand, the differentiation in height growth we found between sites with high and low soil water availability, must therefore stand. Antón-Fernández et. al (2011), also pointed out, that the definition of the dominant height and the estimation of SI are sensitive to initial planting density.
Here we decided to use the dominant height ( 100) and site index curves by Assmann and Franz (Assmann and Davis, 2014) as they are the most common ones applied in Bavaria. Those are especially well suited for spruce monocultures (see Figure 2), as they have been originally derived from the research experiments used in this study.
Our findings, that dominant trees reduce height growth especially on poor sites, is also substantiated by findings of Carl et al. (2018), Ding et al. (2017) or Pretzsch and Biber (2010).
Under water stress, the growth partitioning of the trees in forest stands becomes more sizesymmetric. The dominant trees get drought stressed, make less use of their privileged access to light, and may reduce growth (Pretzsch et al. 2018).

Ecological implications
These results also manifest the specific physiological characteristics of Norway spruce. Norway spruce is very sensitive to competition effects (Klimo et al., 2000). The shallow rooting system (Puhe, 2003) enhances the competition for water and mineral nutrients in the soil, which is particularly relevant in monocultures (Schmid and Kazda, 2002). Especially when repeatedly cultivated on the same sites, as it is common in many European forests, the shallow root system can contribute to soil compaction and acidification that may further reduce the access to the storage of water (Wiedemann, 1923). This results in a strong height zonation and root concentration in the soil, therefore reducing the portion of soil for water absorption. Thus, the described effect of stand density on growth on dry sites, may be less than for other deeper rooting species such as European beech or Scots pine (Gale and Grigal, 1987).
Xeromorphism of needles may be a feature to preserve water in the tree, once the stomata have closed. However, Spruce has proven to be more drought-susceptible than for example beech (Kölling et al., 2007), despite xeromorphic foliage. Under non-limiting water supply, the lower leaf-level transpiration rate of spruce is counteracted by higher leaf biomass and leaf area index at the stand level relative to most other tree species in Europe (Ellenberg et al. 1986;Lyr et al. 1992).
Moreover, Norway spruce is an isohydric species, which reduces stomatal conductance at early stages of soil drought in order to prevent cavitation and temporary loss of the water conducting system. Leaf stomata in trees with higher resistance (taller trees) close earlier in the day or earlier in a drought cycle (Ryan and Yoder, 1997). Early stomata closure under drought means reduction of carbon uptake and growth. In addition, xeromorph needles contribute to better preserve water in the tree, under water stress. Both may increase the revealed height growth reduction under water limitation, especially under high competition for water due to high stand density. This means that periods with a reduction in height growth, due to stomata closure, will be longer for higher densities with increasing water stress, and thus stronger competition for it. In summary, we could observe how especially sensitive Norway spruce is to stand density, decreasing height growth when the limiting resources are in the soil, but not so much when water availability is sufficiently high so it still can allocate more resources to compete for the light.

Consequences for forest management
An underestimation of the site index of ~4 m for density changes in the range of 200-300 trees/ha on moist limiting sites, may underestimate the prognosis of future productivity (> 2 SI-classes), and increase impression in annual cut estimates, and associated forest management decision making. Norway spruce is still the main productive tree species in Southern Germany, so such consequences may translate in important deviations, which will potentially increase as climate becomes dryer.
The effects observed on higher densities, above ~1000 trees/ha, are particularly affected by the soil conditions and precipitation (especially according to model 2). This range of densities, which correspond mostly with a lower age range ( Our results may also explain the sit-and-wait phenomenon already reported for lodgepole pine on dry sites and high density (Cole 1975, Cole andKoch, 1995)). Due to extreme water competition under such circumstances height growth nearly stops. A density reduction may finally improve the water supply of the remaining individuals and cancel the standoff condition (Sohn et al., 2013). In our stands, the water stress of Norway spruce was not so extreme to completely stop height growth, but it substantially reduced it. For this reason, the reduction of local density in spruce stands under drought periods could then reduce the stress caused by the lack of water. This, at the same time, would improve height growth, especially on sites where the availability of water in the soil is generally low. Moreover, these results suggest that maintaining stands at medium density level (reducing rotation length and minimize risk of exposure to extreme climatic events), may serve as a potential effective measure to reduce water stress (Trouvé et al., 2017).
In the natural range, mountain regions with high precipitation, the previously discussed processes of water stress, height growth reduction, reduced competitiveness due to drought stress and its shallow root system may be irrelevant. However, on moisture limiting sites and particularly during dry years,, those ecological traits and limitations become critical in the lowlands. Moreover, intraspecific neighbourhood can mean niche similarity, niche restriction and limited resource uptake and reduction of overall competitive ability.

Main conclusions
Especially under dramatic climatic change, well-established forestry rules like the independence of the site index estimation from the local stand density must be questioned and revised. In this work, we carried on a study to tackle this issue, using the long-term experimental plots from Bavaria, and we chose thinning experiments of Norway spruce monocultures, one of the economically most important but also sensitive tree species. Here, we saw that especially on poor sites with shallow soils, local density must be considered in the estimation of . Such effects may result in a misleading , when it is used as criterion for reference plot selection as, especially on such sites, site effects may be mixed with density effects.
The special physiological characteristics of Norway spruce makes it especially sensitive to drought conditions when growing in monocultures. For example, unmanaged unthinned plots of thinning experiments may develop lower site index with progressing stand development, which may suggest site differences between the plots, even though the sites are similar. Moreover, especially on dry sites, instead of using traditional models to predict growth, spatial explicit growth models may be more suitable given that they account for the effects of local density, and therefore biases in the estimation of the dominant height can be corrected as indicated in this study.            predicted results are displayed for the soil type "dry" and on the right side for the type "moist" for intervals of total precipitation during the vegetation period (TotPrepVP) of 100 mm. predicted results are displayed for the soil type "dry" and on the right side for the type "moist" for intervals of total precipitation during the vegetation period (TotPrepVP) of 100 mm.