Germany ranges from the North Sea to the Bavarian Alps in the south encompassing lowlands in the north, uplands in central Germany and low mountain ranges mainly in southern Germany. 30 % of the area is covered by temperate forests. In the northwest and the north the climate is oceanic, whereas in the east there is a strong continental influence. In central and southern Germany the climate varies from moderately oceanic to continental. The Alps and some low mountain ranges have a mountain climate with lower temperatures and higher precipitation.

National Forest Inventory Data

To estimate biomass increment data of the second (2002) and third (2012) NFI were used. NFI in Germany is based on a permanent nationwide 4 km × 4 km grid. Each grid point in forest area is the center of an angle-count sampling ( BMELV, 2011 ). There are trees that were included in the angle-count sampling (basal area factor 4) in NFI 3 but had not been thick enough to be included in NFI 2. Other trees were measured for the NFI 2 but were missing in the NFI 3. Diameter at breast height (dbh) and height of these trees were predicted for the middle of the period between NFI 2 and NFI 3 ( Jenkins et al. , 2001 ; Dahm, 2006 ) using the function of Sloboda ( Riedel et al. , 2017 ). Thus, plots where thinning occurred between the inventories are included in our dataset. However, plots where all trees that had been surveyed in NFI 2 were missing in NFI 3 due to harvest or mortality were excluded. For the study plots with a basal area proportion of spruces ≥ 70 % and stand age (calculated as mean of the age estimations of the sample trees weighted by the stem numbers per ha that they represent) between 30 and 150 years were selected. Plots where the climate signal is likely to be confounded by extreme soil characteristics (gley soils, pseudogley soils and moor soils) were discarded. Finally, 3830 plots remained for analysis.

Above-ground wood biomass was estimated using species-specific functions of dbh and height. We chose the functions of Zell (2008) , as they were developed based on German NFI data. The functions estimate total above-ground wood biomass, i.e. comprise both stem wood biomass as well as branch biomass. The increment of above-ground wood biomass per year was determined for each tree as the difference between NFI 3 and NFI 2 divided by the period length. These values were extrapolated to 1 ha and summed up at plot-level resulting in one ΔB assigned to each plot ( Jenkins et al. , 2001 ; Dahm, 2006 ). In summary, ΔB represents total above-ground wood biomass net growth of the stand, i.e. turnover of plant organs is not considered. A detailed description of how ΔB is derived based on the angle-count sample is presented in Appendix 1.

The stand density index of Reineke (SDI) ( Reineke, 1933 ; Zeide, 2005 ) with an exponent of -1.605 was used as a measure of stand density. As tree species differ in their requirements of growing space, SDI values are species specific. In order to allow comparisons between different species or to use the SDI for mixed stands, it is necessary to weight the SDI. For each species the 95-percentile of the SDI distribution of pure stands was determined. Weighting factors were calculated dividing the 95-percentile value of spruce (used as reference species) by the 95-percentile value of the respective tree species. For each NFI plot species specific SDI values were multiplied by the weighting factors and then summed up to the overall SDI of the respective plot. A detailed description of the calculation of the SDI is presented in Appendix 2. Statistical values of the NFI data are summarized in Table 1.

Environmental Data

Regionalized daily climate ( Böhner et al. , 2018 ) and soil data ( von Wilpert et al. , 2017 ) are available at the NFI plots. Based on the daily climate data, for each NFI plot annual values of the climate variables presented in Table 2 were calculated and then averaged over the measurement period between NFI 2 and NFI 3.

In order to address the first research question (What stand variables explain the variability in ΔB for a given site index?), we explored the variation in ΔB not already explained by SI. We aimed at identifying the stand and tree characteristics that differ between plots of greater and lesser ΔB but of the same SI: Is greater productivity mainly due to greater stem numbers or do stem num­bers not differ that much, but trees are thicker and radial growth of single trees is faster? First, for each plot SI was determined by estimating the top height and extrapolating it to age 100 applying the Chapman-Richards function ( Brandl et al. , 2018 ). Second, a generalized additive model (GAM) was fitted explaining ΔB in dependence on SI (package mgcv ( Wood, 2011 ) in R 3.3.2 ( R Core Team, 2016 )). Stand age was included as additional covariate in order to account for the influence of age on ΔB. The residuals of this model correspond to the variation in ΔB not explained by SI and age. Third, we divided the residuals in quartiles and tested if stand and tree parameters differed significantly between the quartiles using Kruskal Wallis and post-hoc Nemenyi-Test (significance level p = 0.01), as the data were not normally distributed. On plot level we considered SDI, stem number (N), standing above-ground wood biomass and quadratic mean diameter (dg), on single tree level we considered height, dbh and relative dbh increment, i.e. dbh increment between NFI 2 and NFI 3 divided by the dbh measured at NFI 2. Relative dbh increment was only assessed for trees measured at both inventories. In order to be able to compare height and dbh of trees of varying ages height and dbh had to be rescaled: A 95%-quantile regression was applied describing height or dbh respectively as a fourth order polynomial of age. Then, each tree's height or dbh respectively was divided by the predicted 95%-quantile of height or dbh respectively at the tree's age resulting in a relative measure independent of age.

Details on the methodology are given in Appendix 3.

We modelled ΔB in dependence on site and stand characteristics using generalized additive models with a gamma error distribution and log-link function. A variety of climate and soil variables was offered to variable selection (Table 2). Climate variables comprise annual precipitation sum (P_yr), summer precipitation (P_wq), precipitation during growing season (P_5to9), mean annual temperature (T_yr), summer temperature (T_wq, Tmax_wm), temperature during growing season (T_5to9), winter temperature (Tmin_cm, T_cq) as well as temperature variability (T_sd, T_range). Soil parameters include base sa­turation (BS), soil texture variables (clay, silt, sand), C/N-ratio (CN) and available water capacity (AWC) of the first 60 cm. We selected the best model of all possible combinations of explanatory variables using AIC as criterion. Combinations including highly correlated variables ( Dormann et al. , 2013 ) had been discarded beforehand.

ΔB strongly depends on stand density. Stand density itself depends both on thinning regime and environmental conditions, since favorable sites allow a greater stand density than unfavorable sites ( Pretzsch, 2002 ). We wanted to find a measure of site productivity that is independent of forest management and solely reflects differences in site quality. Stand density could be included as covariate into the model and set to a fixed value for predictions. But as stand density is not independent of site quality, it weakens the explanatory power of the environmental variables. Therefore, we tried to separate the effect of environmental conditions from the effect of forest management on stand density. Our approach follows the methodology applied to experimental plots when characterizing density on plots of varying thinning grades on the same site. Density of a given stand is expressed by the ratio of the basal area of the stand and the maximum basal area observed on the same site ( Pretzsch, 2002 ). Regarding NFI plots as a huge experimental design we identified plots of similar site conditions using k-means clustering. The k-means method partitions the observations into a specified number of groups (i.e. clusters) so that the sum of squares from the observations to the assigned cluster centers is minimized. Based on a comprehensive set of climatic (P_yr, P_wq, P_cv, T_yr, T_wq, Tmin_cm, T_sd, ETpm_5to9) and soil variables (BS, AWC, silt, sand, CN) observations were assigned to 21 clusters using the algorithm of Hartigan & Wong (1979) and trying 1000 initial random sets of cluster centers. The optimal number of clusters had been determined according to the Bayesian information criterion for expectation-maximization, initialized by hierarchical clustering for parameterized Gaussian mixture models using the R package mclust ( Scrucca et al. , 2016 ). For each cluster, interpreted as one experimental plot with a set of comparable site conditions but different thinning grades, the 95-percentile of the SDI distribution (SDI ₉₅ ) was determined. Again we chose the 95-percentile instead of the maximum in order not to give potential outliers too much influence. Still, SDI ₉₅ is interpreted as the maximum SDI that can be reached under the corresponding site conditions. Then, for each NFI plot the ratio of its SDI and the SDI ₉₅ of the corresponding cluster was calculated resulting in a relative density (RD) that reflects the effect of thinning on density. ΔB can then be explained by RD, age and the environmental variables.

The model's predictive performance was evaluated by calculating root mean squared error (RMSE) based on a 10-fold cross validation (data splitting train data : test data = 9 : 1) ( e.g. Mellert et al. , 2016 ). Besides, we checked for systematic errors by determining the slope of a least squares regression without intercept of observed ΔB against predicted ΔB both at the scale of the linear predictor, i.e. the log scale ( e.g. Dolos et al. , 2015 ).

For external validation independent data of 78 long-term experimental plots on 14 locations in Bavaria were available. From these data the increment periods which were close to the inventory periods of the NFI were used. Table 3 comprises the stand characteristics of the experimental plots.

In the data there was a clear trend to larger quadratic mean diameter (dg), dbh, standing biomass and ΔB with increasing SI. Thus, in general greater ΔB coincided with higher SI. However, there was considerable variation in ΔB that was not explained by SI and stand age. This residual variation could be related to stand variables (Table 4): Differences in ΔB were largely due to differences in stand density (Fig. 1a). Sites with greater ΔB generally had a higher stem number per ha, whereas there was no clear trend for quadratic mean diameter. Standing above-ground wood biomass significantly differed between the quartiles of the distribution of the residuals and showed an increasing trend. Trees on sites with greater ΔB but same SI did not have greater single tree diameters on average, but relative dbh increments were higher (Fig. 1b). There was no clear trend in single tree heights. Thus, in general, at a given SI greater ΔB was mainly due to greater stand density: Production was higher, because stem number and standing biomass was higher. In addition, faster dbh-growth contributed to the greater ΔB.

The final model can be described with:

where f denotes a regression spline (Table 5).

ΔB strongly depends on stand density and stand age, but plausible effects of site conditions can be fitted as well (Fig. 2). Relative density (RD) has a strong, approximately linear positive effect. ΔB decreases with increasing stand age from 30 years onwards (for model fitting plots with stand ages between 30 and 150 years were used). ΔB increases with rising mean annual temperatures (T_yr). The increase is stronger in the low and medium temperature range, whereas the slope flattens at higher temperatures. The effect of precipitation is smaller. Low summer precipitation (P_wq) clearly limits ΔB. As above a value of about 800 mm confidence intervals become very wide, no conclusions should be drawn from the subsequent curve progression. ΔB is reduced at both extremes of temperature seasonality (T_sd). Optimum ΔB is reached at medium base saturation (BS), whereas high base saturation has a negative effect on ΔB. To a lesser extent low base saturation reduces ΔB as well. Low sand content (sand) has a positive effect on ΔB, whereas the effect of very high sand content is negative. ΔB decreases nearly linearly with rising C/N ratio.

Cross-validation resulted in a RMSE of 1.996 t ha ^-1 yr ^-1 . The slope of the regression of observed against predicted ΔB was nearly 1 (0.988). External validation of the model with an independent data set revealed that differences in ΔB can be predicted quite reliably (Fig. 3). The R² of the linear relationship is 0.753. RMSE was 1.652 t ha ^-1 yr ^-1 .

As the trend to structurally diverse mixed stands and thinning from above reduces the informative value of SI ( Pretzsch, 2009 ), it makes sense to look for complementary measures of site productivity (Bon­temps & Bouriaud, 2014). We chose above-ground wood biomass increment (ΔB): On the one hand, ΔB encompasses height and dbh increment as well as stand density, and on the other hand it is feasible to estimate ΔB based on NFI data. We preferred ΔB to volume increment for two reasons: First, it constitutes a physiological measure. Second, wood density is taken into account which facilitates the comparison between different species. The downside of the use of ΔB is that its calculation draws on assumptions about allometric relationships between different tree compartments. In comparison to stem volume increment, this leads to higher uncertainty in the estimated productivity measure.

ΔB serves as an indicator or proxy of site produc­tivity. Therefore, when interpreting our results, we relate them to productivity. However, it has to be kept in mind that there are more aspects to net primary productivity like below-ground biomass growth and turnover of plant organs that are not taken into account. Both the allocation of NPP to different tree compo­nents as well as the turnover depend on stand density and stand age. For instance, declining woody biomass increment towards older ages does not necessarily mean that NPP is declining in the same way, but it likely reflects a change of allocation between stem biomass and the rest of the tree. This has to be kept in mind when interpreting the results.

Sites with similar SI and stand age showed no­ticeable variation in biomass increment: Greater ΔB was mainly due to higher stem numbers, reinforced by larger relative dbh increments. If more productive sites at similar SI and age differ more in stem number and only to a lesser degree in diameters from less productive ones, site productivity is better captured looking at ΔB of the stand than at the increment of single trees alone or mere stand height. It has to be kept in mind that this effect was found for sites of similar SI and is not a general principle. When looking at the entire data set i.e. the whole range of site indices and ages there is a clear trend to larger dg with increasing SI. The differences in productivity at same SI cannot immediately be traced back to differences in site con­ditions and thus be interpreted as subdivided specific yield levels, as most forests in Germany are managed and therefore differences in stand density leading to differences in productivity are mainly due to thinning. Still, maximum stand density i.e. carrying capacity on a given site depends on site conditions ( Pretzsch, 2002 ). Favorable sites would show greater dbh increment than unfavorable sites given the same stand density. But as forest owners might tend to keep higher stem numbers at favorable sites, better site conditions are sometimes not expressed as much in greater dbh increment but in higher stand density. Thus, exploring the relationship between the variability in ΔB and stand variables at a given SI and stand age emphasized the importance of adequately dealing with stand density when modelling ΔB. Therefore, we differentiated between management effects and environmental effects on stand density by calculating a relative density in the modelling approach. This allowed us to develop a model that separates the effects of thinning from the effects of site conditions on productivity. Of course, this is an idealization as the effects of thinning regime and site quality can never be separated completely and there are many influences on stand density not encompassed by the explanatory variables used in this study. We modelled ΔB in dependence on age, relative stand density and environmental variables in one step and can predict potential ΔB by setting the relative density and age to reference values, just as height can be modelled in dependence on age and environmental variables in one step and SI can be predicted by setting age to a reference age ( e.g. Brandl et al. , 2014 ; Vallet & Perot, 2016 ). Actual ΔB can be predicted by setting the relative density and age to the current values of the given plot. An alternative approach would be to first estimate a productivity measure detrended from age and density effects and in a second step model this detrended productivity index in dependence on environmental variables ( e.g. Watt et al. , 2010 ; Charru et al. , 2014 ), which is in analogy with the approach of first deriving the SI of a stand, i.e. detrending height of the age effect, and then modelling SI in dependence on environ­men­tal variables ( e.g. Albert & Schmidt, 2010 ).

Overall the model shows a high goodness of fit and validation on an independent data set showed that it reliably predicts differences in ΔB. The effects of age and relative density on ΔB in the model are clear and ecologically plausible: Since stand density is directly connected to leaf biomass ( Pretzsch et al. , 2014b ), dense stands reach maximum leaf area and thus maximum light interception ( Zeide, 2001 ). Therefore, it makes sense that productivity increases with increasing stand density. This result is in contrast to Pretzsch (2006) who found a unimodal optimum relationship between stand density and growth. This contradiction might be due to our use of NFI data instead of data of experimental plots. As most German forests are managed the pro­portion of unthinned NFI plots with such high stand densities as to cause reductions in growth is too small to influence the model effect.

One would expect net primary productivity for a given stand to increase until an age of about 50 years and then decline again due to the changing balance between gross primary productivity and respiration during stand development ( Barnes et al. , 1998 ). But in this study ΔB declines monotonously with stand age within the age range considered (30 until 150 years). On the one hand, this might indicate that the age dependence of above-ground wood biomass increment differs from the age dependence of NPP due to changes in allocation with age. On the other hand, it can be explained by our use of cross-sectional data instead of time series, i.e. we did not follow the trajectory of one stand through time. Plots of the same age can differ in their developmental stage ( Mehtätalo, 2004 ). Replacing age by dominant height as an indicator for developmental stage reveals the expected pattern with an increase in ΔB at low dominant heights followed by a slow decline at greater heights (not shown). Still, in order to compare and predict site productivity it is preferable to use stand age in the model ( Mehtätalo, 2004 ).

Adding climate and soil parameters as explanatory variables renders plausible effects on ΔB. On a global scale aboveground NPP is relatively low in cold and dry climates and rapidly rises as both temperatures and water availability increase ( Barnes et al. , 1998 ). This global-scale pattern can also be observed on a German scale, although we are looking at ΔB here. Temperature regime and water supply clearly are limiting factors, and as both increase, productivity rises. The most influential environmental factor in our study is mean annual temperature. This is consistent with other current studies. For instance, Pretzsch et al. (2014a) concluded that mainly rising temperatures and extended growing seasons increase growth. Kauppi et al. (2014) identified the spatial and temporal variation of growing degree days as the main causal factor affecting variations in forest growth. At the high end of the temperature range the increase of ΔB with rising temperatures slows down and approximates a more or less constant level. One might expect a decline at very high temperatures due to drought stress ( Dolos et al. , 2015 ). However, due to high risks and adapted forest management spruce dominated stands in Germany simply do not occur in sufficient numbers at very high temperatures in order to clearly detect such an effect. On a global scale water supply is a crucial variable constraining biomass ( Stegen et al. , 2011 ). The effect of precipitation in our study is rather weak, as water is, except in extreme drought years, not generally the growth limiting factor in our data set (spruce dominated NFI plots in Germany): Annual precipitation of 92 % of the plots exceeds the threshold value of 800 mm given by Mayer (1992) for the optimum growth range. Still, ΔB decreases when summer precipitation is low. Within the same climate ΔB differs, since it is influenced by soil properties, species composition and the stage of ecosystem development ( Barnes et al. , 1998 ). The effect of base saturation on ΔB follows an optimum relationship. On acidic soils the supply of basic cations reduces growth, whereas on calcareous sites Ca-K-antagonism ( Rehfuess, 1990 ) and immobilization of phosphor ( Mellert & Ewald, 2014 ) can occur. Low sand content has a positive effect on ΔB, whereas high sand contents affect ΔB negatively. The effect of sand content might both reflect effects of nutrient and water supply. Soils with high sand content often have a low available water capacity and are poor in nutrients.

The proportion of explained variance by environ­mental variables is small, but therein comparable with other studies ( e.g. Condés & García-Robredo, 2012 ; Charru et al. , 2014 ). If we could look at total volume production the effect of site conditions on productivity would be accumulated over the whole life of the stand. The same applies to stand height. Differences in site conditions cannot be reflected as distinctly in ΔB between a time span of 10 years. For instance, when looking at a rather short time span, it is more likely that weather variability between the years does not reflect average climate conditions and thus blurs the effect of climate on growth. However, this time span in combination with corresponding climate data allows to assess short-term growth response, which can also be perceived as an advantage of this approach. Environmental data are regionalized and thus introduce uncertainty into the analysis. Environmental influences on forest growth must be summarized into a few quantifiable factors. It is no wonder that their effect is small considering the complexity of tree growth. Complex interactions between site conditions and forest management, extreme events as well as genetic variability may greatly affect productivity.

Comparisons with studies about site factors in­fluencing biomass ( e.g. Chave et al. , 2003 ; Keith et al. , 2009 ; Stegen et al. , 2011 ) are only possible to a certain degree, as biomass and biomass growth may react differently to environmental influences. For instance, an extension of the growing season will increase biomass growth as long as water supply is not limiting. Forest biomass may stay the same, since trees only move faster along their life's trajectory and die at a younger age, but self-thinning lines remain constant ( Pretzsch et al. , 2014a ).

Productivity is often estimated based on height information alone ( e.g. Seynave et al. , 2005 ; Albert & Schmidt, 2010 ; Nothdurft et al. , 2012 ), thus taking only the vertical aspect of productivity, i.e. height growth, into account. The results of this study illustrate the importance of the horizontal aspect of productivity, i.e. diameter increment and the strongly correlated branch increment as well as stand density, as sites that do not differ significantly in SI and age can still differ in productivity. Recent analyses of Norway spruce stands in Bavaria (Southern Germany) based on NFI data could be interpreted in the light of these findings: Based on NFI data similar site indices are estimated for the two Bavarian forest eco-regions Swabia and Spessart. A SI-model based on Bavarian NFI data also predicts similar site indices for these two regions ( Brandl et al. , 2014 ). However, in forestry practice Swabia is generally considered the better site for spruce. Looking at our data we found that sites of similar stand age (Spessart 80 years, Swabia 76 years) and SI (Spessart 36.7 m, Swabia 36.4 m) in Swabia indeed have greater above-ground wood biomass (Spessart 296 t ha ^-1 , Swabia 403 t ha ^-1 ) and show greater above-ground wood biomass increment (ΔB) (Spessart 7.8 t ha ^-1 yr ^-1 , Swabia 9.9 t ha ^-1 yr ^-1 ). In contrast to the mentioned SI-model our model predicts significant differences in productivity. Actual ΔB can be predicted using actual stand density and age (Spessart 8.4 t ha ^-1 yr ^-1 , Swabia 10.4 t ha ^-1 yr ^-1 ). Setting a fixed age ( e.g. 80 years) and a fixed relative stand density ( e.g. 0.7) potential ΔB can be predicted (Spessart 8.7 t ha ^-1 yr ^-1 , Swabia 9.5 t ha ^-1 yr ^-1 ) resulting in a difference of 9.2 % due to climate and soil. Potential ΔB cannot be compared to the measured values, but reveals differences in site potential. This example illustrates the benefit of not only looking at SI but also at ΔB when assessing site productivity.