Climate Sensitive Tree Growth Functions and the Role of Transformations

The aim of this study is to develop climate-sensitive single-tree growth models, to be used in stand based prediction systems of managed forest in Switzerland. Long-term observations from experimental forest management trials were used, together with retrospective climate information from 1904 up to 2012. A special focus is given to the role of transformation of modelling basal area increment, helping to normalize the random error distribution. A nonlinear model formulation was used to describe the basic relation between basal area increment and diameter at breast height. This formulation was widely expanded by groups of explanatory variables, describing competition, stand development, site, stand density, thinning, mixture, and climate. The models are species-specific and contain different explanatory variables per group, being able to explain a high amount of variance (on the original scale, up to 80% in the case of Quercus spec.). Different transformations of the nonlinear relation where tested and based on the mean squared error, the square root transformation performed best. Although the residuals were homoscedastic, they were still long-tailed and not normal distributed, making robust statistics the preferred method for statistical inference. Climate is included as a nonlinear and interacting effect of temperature, precipitation and moisture, with a biological meaningful interpretation per tree species, e.g., showing better growth for Abies alba in warm and wet climates and good growing conditions for Picea abies in colder and dryer climates, being less sensitive on temperature. Furthermore, a linear increase in growth was found to be present since the 1940s. Potentially this is an effect of the increased atmospheric CO2 concentration or changed management in terms of reduced nutrient subtractions from forest ground, since industrialization lowered the demand of residue and slash uptake.


Introduction
Projections of tree growth support management decisions and are the basis for forest policy.Pro-active and adaptive managers [1] are particularly interested in effects on future productivity under a changing climate [2].Single tree stand simulators (e.g., Prognosis model [3], Prognaus [4], Silva [5], BWinPro [6], SwissStandSim [7]) are useful tools and the basis of forestry decision support systems [8] in the development of future management strategies.One of the core processes of tree simulators is tree growth [3].
Tree growth is a complex process, depending on multiple and interacting factors.It was shown to be influenced by site conditions [9,10], climate [11], competition and management [12] and also nitrogen deposition [13].Some single tree growth models share analogies to the model developed from Wykoff [3], where groups of explanatory variables are used to describe biological meaningful reference to tree growth [4,14,15].Site is then usually modeled as being constant over time, explained by location specific variables like height above sea level, latitude or longitude.Just recently the importance of climate change on tree growth helped to intensify the modeling of retrospective climates, enabling more flexibility in modeling tree growth under changed climate [16,17].
In many empirical studies [17][18][19] the climate effect depends on a substitution of space for time (i.e., climate observed over a large spatial gradient within a short time-frame).Predicting with those models is assuming that a change in climate is the same as moving in space.Rohner et al. [17] for example show that site, climate, competition, management and nitrogen deposition contribute significantly to tree growth over large gradients in central Europe.Another approach to quantifying effects on productivity are physiological process-based models that rely on biological principles [20].Such a model was e.g., used by Reyer et al. [21], and indicates that productivity can increase or decrease under a climate change in central Europe, depending on the tree species.For both approaches, long-term observations are essential in providing validation evidence to support or reject findings.
Over the long rotation periods of forests, such as in central Europe, changes in climate potentially occur more rapidly than forests are able to react with their natural adaptation potential [22].The effects of climate change in forests are therefore difficult to assess via typical experimental design.Long-term permanent observations of tree growth have the potential to quantify climate and many other effects simultaneously.Empirical models, fitted to such data [5,23,24] help our understanding of the ongoing shifts in growth and enable the prediction of possible future forest states.
Diameter at breast height (d) is a precise and easy measure of tree dimension that contains information on individual tree growth, over long terms and large regions.Growth can be expressed by diameter or basal area increment (bai).In terms of predictive precision, it does not matter which is being modeled [25].Two conceptually different approaches are common when modeling bai for single-tree stand simulators.The semi-empirical approach, in which maximum growth is reduced by means of generally accepted principles, often makes use of previously identified parameters and models [5,26].This is in contrast to the direct empirical approach, where mean growth is modeled in one equation while simultaneously taking into account all potential explanatory variables.This is referred to as the "composite-model" approach in Wykoff [3] and is also used in [4,17].The latter was favored in the present study, since it requires fewer assumptions, and a best model of growth can be identified in a straightforward manner through statistically valid model selection.
Bai is nonlinearly related to the size of a tree, measured by diameter at breast height (d) or age.Consequently, a log-transformation (base e) is often used [3,4,[27][28][29][30][31] to linearize parameter estimation by means of linear regression.This comes with many advantages such as the estimation problem now being identifiable, a global optimum exists, residuals are homogenized and model selection can be automatized.In addition, ridge regression can be used to analyze multi-collinearity [3] and lasso regression may be used for model selection [32].However, the log-function has a steep gradient at small numbers, resulting in very strong transformations of smallest bai values.Wykoff [3]-using log-transformation-reports on systematic overestimation of small trees.Potentially this could be an effect of the transformation.
The data used in the present study showed that a nonlinear model formulation is relevant, and a model formulation according to Quicke [33] was used.Different transformations were used in order to normalize the residuals distribution and hence stabilize the parameter estimation.This helped to construct a wealth of different models with increasing model complexity, meaning that the original formulation of [33] was widely expanded by groups of explanatory variables, resulting in more complex models.
The question of which transformation to choose for tree growth data is not often discussed or analyzed in tree-growth in the literature.It has importance since the different transformations change the error distribution of the models and hence influence statistical inference.Technically it is possible to use the most flexible transformation by the formulation given by Box and Cox [34], where the transformation is described by an estimable parameter.The flexibility of this transformation results in functions, which are nearly impossible for a direct understanding and interpretation, like applied in Fischer [35], with examples of tree growth.Therefore, the focus was on two "families" of transformations: one is based on the log-transformation and the other one on root-transformation.The log-transformation needs strictly positive numbers, and hence the negative observed bai needs to be shifted into the positive range, e.g., by adding a small number (α) e.g., [30,35,36].However, the potential influence of α on the mean squared error was not investigated so far.
The overall aim of this study was to find climate-sensitive single-tree growth models, with good predictive ability.Simultaneously, all information stored in the long-term observations of trees should be used to construct models.All potential and interpretable drivers of tree growth (i.e., between-tree competition, stand development, site, thinning, mixture and climate) should be incorporated for use in a single-tree stand simulator.Further transformations were investigated and used to predict bai unbiased and with minimum mean squared error over the full spectrum of tree sizes.

Experimental Forest Research Data
The tree-and stand-specific data were extracted and calculated from the Experimental Forest Management (EFM) trials in Switzerland, including the first systematic growth and yield research plots that were established, containing re-measured trees since 1880.The plots had a standard area of 2500 m 2 and the standard interval length for re-measurement was 5 years.Immediately before inventories, managers mark trees for harvesting, so that the intensity and timing of the intervention is precisely recorded.The entire data set included 374 plots, at altitudes between 385 and 1930 m a.s.l., covering a wide range of climate types.They also cover the Swiss ecophysiological regions of Jura, Plateau, Pre-Alps and Alps.Their locations are spread all over Switzerland and are shown in Figure 1.Lat The single trees were marked individually for re-measurement.Breast height diameters of all trees in the plots were measured crosswise with a caliper at a height of 1.3 m.The bai was then calculated for the standard 5-year interval using: In order to use as much of the information inherent in the data as possible, negative values were not excluded, since they are not, per se, outliers.Only a small data trimming of 0.1% was applied (i.e., discarding a fixed proportion symmetrically on both sides of the distribution of bai).That resulted in a total of of 574,000 observations, and an overview on the range and quantiles of the data can be found in Table 1.However, still after trimming, a large number of potentially influencing data points were still present.The data contain errors, which are not only due to the measurement error of the caliper (e.g., tree confusions, transposed digits, . . .).The large amount of data and their long history made it impossible to reconstruct the supposed true measured numbers.To handle the contaminated data, robust statistics was applied [37].This method uses most of the information inherent in the data, by an appropriate and objective construction of weights.

Explanatory Variables
The stand-and tree-specific calculations resulted in 47 potential explanatory variables, which are listed with abbreviations, calculations and units in the chapter Abbreviations at the end of the text.They are grouped for interpretation and modeling, with an overview on the grouping in Table 2 and summary statistics in Table 1.

Competition
The competition variables, based on diameter, were basal area of larger trees (bal), relative basal area of larger trees (rbal) and the cumulative basal area (cba).A crown competition index (cc) was also calculated.This was evaluated at 66% percent of the subject tree's height (0.66h), based on Biging's demonstration [12] that this relative height is best in terms of predictive ability to describe competition in growth.It relates the sum of cross sectional areas (csa) at 66% height of a particular tree j to total area of the stand (A).
The derivation of crown sectional area (csa) is based on assumptions regarding crown form, total tree height, crown base and crown width.The crown forms from Pretzsch ([38], p. 307) were used.Tree height, crown base and crown width were measured only on a sub-sample and were imputed with tree-specific models, based on the observed values.

Stand Development
The stand development was described by stand age (a), the quadratic mean diameter of the thickest trees ddom, and total crown surface area of a stand cs.

Site
Historical N-deposition (in kg/ha/a) (ndep) was estimated for each plot using a modeling approach.It combines emission inventories, dispersion models, spatially interpolated monitoring data and deposition models for nitrogen compounds [39].In addition, a set of factor variables describing soil characteristics were derived by a spatial join with geological maps (https: //www.bfs.admin.ch/bfs/de/home/dienstleistungen/geostat/geodaten-bundesstatistik/bodennutzung-bedeckung-eignung/abgeleitete-und-andere-daten/bodeneignungskarte-schweiz.html)(1:200,000), resulting in factors for water holding capacity, water permeability and nutrient availability with three levels each.

Stand Density
A number of explanatory variables were calculated as potential measures of stand density.They included stand density sdi as defined by Reineke [40], total basal area per hectare (bat) and total stem number per hectare (nt).Minimum diameter at breast height was zero cm.

Thinning
Potential thinning effects were expressed by a dummy variable treat.It defined whether thinning was performed in a previous inventory.Thinning intensity was described by the basal area of harvested trees (bah) as a percentage of its total.Similarly, nt contained the relative number of trees being cut.The change in the crown competition index ccc was also a measure of thinning intensity.Therefore, cc calculation was performed before and after harvesting took place.

Mixture
Mixture effects were calculated by two variables.mixB contains the basal area share of broadleaved trees in a stand, whilst lead.con is a dummy variable, with a value of 1 if the leading tree species, in terms of maximum basal area, is coniferous.

Climate
For each EFM plot, monthly temperature and precipitation for the years 1900-2012 were calculated by Remdund et al. [41].Yearly precipitation sum (p, in mm) and yearly temperature (t, in • C) were then calculated over period length, as given by the re-measurement intervals on the plots.Instead of calender years, we used physiological year spans, meaning that the sum was taken from October until September, which is understood to have a better signal on bai [42].Furthermore, an aridity index (as defined by deMartonne ( [38], p. 520) was calculated, where A = p t+10 .Since this index means smaller values have higher aridity, its inverse interpretation is more intuitive.Therefore, it was labeled "moisture-index" (m).Using yearly values in this case is not independently measuring moisture; in the model it would just be an interaction between yearly temperature and precipitation, therefore moisture was calculated over spring months (March-June).

Years Past 1940
After the first round of final models, the residuals showed a clearly accelerating growth trend after around 1940, which was present for all species, except pine.Therefore, a variable was added containing the years past 1940 (calculated as current year−1940, otherwise 0).

Base Model
A decreasing trend in bai over d was not found in the data (as was assumed in [3] or observed in [43]), see Figure 2. As a result, a sigmoidal curve, approaching an asymptote over d was the base model.Its dependence on d is the same as in Quicke [33], but largely expanded by other explanatory variables.This model formulation is intrinsically nonlinear.Nonlinearity in parameter estimation has some drawbacks, however.Heteroscedasticity requires the quantification of variance for the appropriate weighting of residuals.In the present investigation, because of the contaminated data, weights have already been used.The complex models, including variance functions, did not always converge (the problem may be unidentifiable).Further, starting values need to be supplied and there is no guarantee that a global maximum can be found.Both of these concerns become more problematic with increasing model complexity.These considerations call for the transformation of the dependent variable.
Log-transformations needs strictly positive numbers, hence the observed values are shifted into the positive range by adding a small number α.The analysis of the growth model thus depends on the appropriate choice of α in connection with the best transformation of the dependent variable.The focus was on a subset of all possible transformations, with reasonable and interpretable functions.
The basic relationship between bai and d is given by: where g(.) and f (.) are functions.g(.) is especially important since it influences the random variable and hence the random error structure of the model.On the right-hand side, β 0 is the intercept, which was sequentially enlarged by linear combinations of potential explanatory variables.In its original model formulation f (.) = exp() and g(.) = Identity, meaning that all further explanatory variables to replace β 0 affect bai in a multiplicative way.In that formulation, the parameter estimation is unstable (β 1 an β 2 are highly correlated), and convergence problems occur, especially when the variance function is estimated simultaneously.As a result, transformation of the random variable bai is a plausible option.Here we tested a set of different transformations g(.), although any transformation f (.) applied to the right-hand side will never result in a linear model as it is intrinsically non-linear in terms of parameter estimation.

Transformation
Since the log-function requires positive values, bai was shifted into the positive range by adding a small constant α.Adding the smallest possible values (min(bai)) will result in a very strong transformation (in this case −∞), resulting in biased predictions.However, very large values of α will also introduce bias.Thus, an optimum value for α can be expected.Box-Cox [34] transformations can be extremely flexible in that they use a parameter λ to describe the transformation.Here we focused on a subset of all possible transformations, with reasonable and interpretable formulae.The following transformations g(.) were tested:

Model Selection
The nonlinear model structure and the large number of potential explanatory variables is problematic since the models can not directly be derived using an automated process.Building all models with 47 potential explanatory variables is infeasible.Therefore, a forward model selection process was started, where groups of explanatory variables were added stepwise.The groups are listed in Table 2. Preferably, only one explanatory variable, which held the most information of the interpretable group, was added at each step.In the case of nonlinear behavior of explanatory variables or interacting explanatory variables, or in the case of site-variables (which are assumed to be independent and interpretable), more variables were allowed to enter the model.In the last round, final models were reduced.Due to the large number of influential data points (the residuals were long-tailed and not normal), robust non-linear methods [44] were applied, within the R environment [45].Since they iteratively minimize the sum of squares, model comparisons were made using F-tests.
At each step, the current full model was further analyzed using a generalized additive model ("gam") [46] to enable a visual inspection of the linearity of the effects.If nonlinear effects were present, quadratic terms were added.Interactions were assumed to be present between total basal area (bat) and stems per hectare (nt) and within the climate variables.

Model Validation
Since the goal of the study is the prediction of bai, the models ability to predict independent data was evaluated by a 10-fold cross-validation.Since it can be assumed that the locations of the plots are independent, the random splitting of the data was based on the locations.In each round approx.90% of the locations were used to build the model, predicting bai for the remaining 10%.Based on the full set of predicted values, the R 2 and root mean squared error (rmse) of the predicted values was computed and compared to the simpler version of the in-sample predictions, based on the reported full models, see Table 2.The 80 resulting models can be found the Supplement "Supplementary2.xlsx".

Role of Transformation
The final models were refitted with the different transformations.Their back-transformed predictions were then compared with the observed values and the residuals sum of squares were calculated (which is the same result as mean squared error, since the number of observations are fixed).Figure 3 shows exemplary for spruce and beech the effect of the transformations.Weaker transformations (square root is weaker than cube root) generally showed better results.Root transformation do not require a shift since their domain includes 0. However, if the shift is included, the residual sum of squares increases with α (not shown in the figure).In the case of log-transformation, there is an optimum value for α.Back-transformation was applied directly and with a correction based on [47] to avoid transformation bias; resulting in the same conclusion.Based on that result, final models are all based on the square root transformation.

Overview
The expected mean growth curves for all tree species are presented in Figure 4, which correspond to the parameters β 1 and β 2 in Table 2, whilst the other explanatory variables where held constant at their mean values, which influence the asymptotes of the curves.The progression in the curves is plausible and all β 2 are negative, thus approaching meaningful asymptotes.The predictions for each model can be calculated by: bai where the . . . is a species-specific sum of the estimated parameter (β 3 − β 45 ) multiplied by the value of the explanatory variable, see Table 2.
The resulting summary statistics for the models on the original scale (back-transformed) can be found in Table 3.The amount of explained variance (in-sample, R 2 in−sample ) was relatively high (74%, 62% and 68% for beech, spruce and fir, respectively) and did not drop very much compared to predicting independent data (R 2 predict ).The same can be observed for the difference between the in-sample and prediction based rmse, showing stability in predicting new data.
In Table 2, the eight resulting final models are presented with their estimated coefficients (standard errors of the models can be found in the supplement, "supp1.docx").

Effects Competition
Competition measured by bal or rbal always had a negative sign.Established trees in a stand show a clear increase in growth, due to better access to resources.The same interpretation was found for cba, where higher values mean that the tree is in the upper quartile of the basal area distribution of a stand, which is related to low between-tree competition.The estimated coefficients were therefore all positive.

Stand Development
To describe the influence of age in the stand development on growth, non-linearity was found to be present.Figure 5 illustrates the effect of age.The blue background is the two-dimensional density of diameter and age as given in the observations, overlaid with the expected growth curves, whilst holding all other explanatory variables constant on their mean values.The effect of age was always negative, while age 2 was positive.For small diameters, the estimated growth is large for small ages, but heavily decreases with increasing age (the shown gradient is steep).The same trend is found with large diameters, but with a much smaller gradient.In the case of fir, ddom was superior compared to age-related effects.This is reasonable, since fir is highly shade tolerant and thus age is potentially a bad measure of growth vigor.For small diameters, an increase in ddom means slightly faster expected growth, while in larger trees this effect is stronger.

Site
Nitrogen deposition (ndep) had an linear positive effect in maple and a very small negative effect in douglas fir-model, see Table 2. Since the effect was also assumed to have a non-linear behavior, the resulting effect sizes are graphically shown in Figure 6.Here an accelerating effect of nitrogen deposition was observed, being most present in beech, oak and larch.The effects were getting stronger around 30-40 kg/ha/a.A point of saturation was not found.
Slope (in %) showed a clear negative effect in the douglas fir-model, while in the spruce-model, it had a small positive effect.Eastness was found to have positive effects, which is in accordance with [17], whilst northness was a negative effect for beech, which is in contrast to the very small positive effect found in [17].The reference level for relief is flat terrain, and is included in the intercept (β 0 ).All other effects are contrasts.Depression had a negative effect for beech, slopes and summits had a negative effect on maple and beech, whilst oak shows the positive effect of summits.Since oak is drought tolerant [48], this effect is plausible.

Stand Density
The stand density index was favored in the case of the maple and larch models, with plausible negative signs.More complex relationships between total basal area (bat) and stems per hectare (nt) were found to be present in the other tree species.Spruce, as shown in Figure 7, revealed a quadratic effect of bat on bai, independent of nt (horizontal lines).For beech and spruce, the increase in bai intensifies with smaller bat.In contrast, this is much less pronounced in fir.Larger bai in beech and fir can be found with small bat and large nt.In the case of beech, this decreases more rapidly with small nt.

Thinning
Treatment (measured by treat and intensity bah) both had reasonable positive signs, except for oak.The negative sign in the change of crown competition in fir makes sense because it describes a lower crown density after thinning.

Mixture
Maple, beech, spruce and oak showed the clear negative effect of mixtures with broadleaved trees (mixB), which is in accordance with the findings in Houpert et al. [49].Positive effects were found when the leading species was coniferous for larch and fir.This is also supported by [49], although only the presence of coniferous trees was used rather than the proportion.

Climate
Interacting and nonlinear climate effects are presented in Figures 8 and 9.The first figure shows temperature and precipitation, overlaid by a contour line of predicted bai.This is indicating that beech and fir growth principally accelerates with temperature.Both higher temperature and higher precipitation had a strong impact on fir growth.In contrast, good growing conditions for beech were found with higher temperatures, but only medium-level precipitations.In the case of spruce, the dependence on temperature was clearly less pronounced as growth is nearly constant across temperatures.
The interacting effects of temperature and moisture as presented in Figure 9 show that higher moisture (in spring months) had a clearly positive effect on spruce and fir, and a medium effect on beech.The spruce optimum was found in cold and moist areas, fir optimum in warm moist areas and beech optimum in warm and medium-high moist areas.Non Significant Effects Some explanatory variables were never included in the final models and are not included in Table 2.They were cc and csa, which are variables that describe growth processes based on crown forms.Only the change of crown competition, ccc, had a strong effect in the fir model.Obviously, the models that impute crown forms were less informative than those with measured values.
The soil related variables wat.hold.cap,wat.perm, nut.av were not present in the final models.This could be due to the coarse resolution of the data (1:200,000) or simply because soil information is poor given that it comes from geological land classification.

Transformation
Transformations of bai helped to stabilize parameter estimation of the investigated complex non-linear models.Convergence problems in model fitting occurred less frequently.In addition, residuals were homoscedastic, making further variance functions unnecessary.Fischer [35] used Box-Cox [34] transformation, which are difficult for direct interpretation.Therefore, logarithm and root functions were used in the present study.The result showed that square root transformation is best in terms of mean squared error, even better than log transformation.This also holds for different added constants (α) to bai, which is necessary to shift the data into the positive range.The gradient of the log-function is very steep in small numbers, potentially resulting in an overestimation of small trees, as e.g., described in Wykoff [3].A back-transformed logarithmic model results in multiplicative effects.Back-transforming a square-root-transformed model results in both additive and multiplicative effects.It may be that this type of transformation comes closer to the nature of the data in describing the growth process.
An analysis of the constant α was not found in tree-growth literature, even though it is commonly applied [30,35,36].The application of α comes with an optimum in mean squared error, which depends on the data constellation.Therefore, it is recommended to test different values of α if it seems necessary to be used.This may be the case when the model assumes that only positive growth is possible but the data contains negative observations, or when transformations appear adequate, but the applied function requires positive numbers.

Estimation Method
Bai data used for modeling contain measurement error from the caliper and potentially confused trees and transposed digits and other sources of errors, originating from the long history of these data.Hence the residuals of the models were long-tailed and not normally distributed.Since it is the goal of the present study to use as much of the information present in the data, a stronger trimming does not solve the problem.For this reason, robust nonlinear statistics [44] was used.In this approach, a weighting is applied such that most sensitive observations of the parameter estimates have less weight.The number of potential outliers was large and the share of down-weighted observations ranged between 19.5% and 22.6% depending on the tree species (see Table 3).For that reason, the repeated observations of trees in plots [50] were assumed to be independent.If the models were able to fully describe the effects at the plot level, this would result in independent residuals.However, plotting the the residuals across the different plots showed systematic deviations between the plots, making a non-linear mixed model [51,52] more appropriate.To date, a combination of non-linear mixed models with robust statistics is unknown.Only robust mixed linear models have been undertaken thus far [53].Incorporating the estimated weights from the robust approach into the mixed-effects modeling approach did not lead to model convergence.The non-linear mixed-effects models would be influenced by the outliers.Thus, in both approaches, we can assume that the estimated standard errors of the coefficients are potentially too small, i.e., overoptimistic.At this point, it is unclear which of the two approaches would yield better predictive accuracy.It is beyond the scope of the present article to develop new methods or to analyze the predictive behavior of the two mentioned approaches under different data constellations.

Explained Variance
The explained variance in predicting independent data is relatively large (74%, 62% and 68% for beech, spruce, fir), compared to Rohner [17] (43%, 37% and 53%).A part of that may be due to the fact that the variance EFM-data is smaller compared to National Forest Inventory (NFI) data since the trials are under full control and the NFI contains all possible forest states.Furthermore, the model selection approach used in the present paper results in species-specific models, and the explanatory variable(s) per interpretable group can be very different between the species.This was avoided in the model selection approach used in [17], where tree models were described with the same set of potential explanatory variables (for comparison between the trees).Additionally, the models in the present paper contain nonlinear explanatory variables and potential interaction terms.Another difference is with respect to the transformation.That which is used in this paper results in more stable parameter estimation, allowing for more complex growth models.

Model Selection
The pre-selection of models in [17] is based on the multi-collinearity of explanatory variables.If the interpretation of the sign of the effect is of central concern, this approach may be feasible.However, it does not necessarily lead to the best model for the purpose of predicting data.Potentially excluded variables still contain information that may enhance prediction accuracy.For example, the effects of stand density on growth may be incorporated into one explanatory variable, in our case sdi, according to [40], condensing the effects of basal area and stem number into one index.However, our analysis reveals that the effect of stand density is often better reflected by the interaction of basal area with stem numbers.It should be noted that the benefits of approach may come at a price.For some parameters to be estimated, a more complex interpretation is required.This holds true for climate effects.Here we used interaction and the nonlinearity of climate variables.Although a direct interpretation of parameters is not straightforward, the Figures 8 and 9 show that interaction can be interpreted and is biologically meaningful.Spruce is known to have good growing conditions in colder and moist climates and fir is known to have good growing conditions in warm and moist climates [48].

Nitrogen
Unfortunately the nitrogen values are estimated based on models and are hence random variables, which in turn means that the reported effects are potentially underestimating the true effect.A correction was not possible, since the variance is not known.

Years Past 1940
At first glance, it is surprising that an accelerating growth effect was found for the years following 1940, even though climate change was incorporated into the models.This was a linear effect that was present in all species except for pine (see Table 2).Since we were not able to explicitly find an explanatory variable for this effect, it can be hypothesized that it is the effect of increased CO 2 in the atmosphere, which can stimulate growth.Alternatively, it could be the result of changes in management.After industrialization, the pressure on forests to collect residues and slash clearly lessened, resulting in soils regenerating from the extraction of nutrients.

Conclusions
In the present study, a variety of effects on bai were quantified for eight species.Most identified effects were plausible and well reflected species-specific physiological meaningful biological interpretation.The models are complex and contain many explanatory variables on tree growth.This results in a high R 2 , which remains on high level for predicting new, independent data.Hence, it can be recommended to use them in predicting, e.g., in single tree growth simulators, especially for the analysis of future changes in environmental conditions with respect to climate and N-deposition.

Abbreviations
The following abbreviations are used in this manuscript: .Since this index means smaller values have higher aridity, its inverse interpretation is more intuitive.Consequently it was labeled "moisture-index" (m).In contrast to precipitation and temperature, moisture shows higher explanatory power in the spring months [42]; hence only the spring months (March-June) over periods were used to calculate the moisture index.

Figure 1 .
Figure 1.Locations of the Experimental Forest Management (EFM) trials in Switzerland, used in this study.Reproduced by permission of swisstopo (JA100118).

Figure 2 .
Figure 2. Observed basal area increment (bai) over diameter at breast height (d) for all species.Fifty-three percent of the variance is explained by a simple model (Equation (3), with g = √ .and f = Identity), using only d as explanatory variable.The blue background is the two-dimensional density of the observed d and bai.Dark blue contains more data, the 100 points indicate areas with the lowest density.

3 Figure 3 .
Figure 3.Residual sum of squares (in millions) for spruce and beech final models under different transformations of growth.Back-transformation was applied directly and with a correction factor[47].The log-transformation depends on α.Horizontal lines show the root-transformations, independent of α.The fourth root transformations for the beech model are outside the plotting range.Their values are 88.5 (direct) and 89.8 (Snowdon-correction).

Figure 4 .
Figure 4. Expected mean growth for tree species.

Figure 5 .
Figure 5. Effects of the stand development measured by age or ddom (quadratic mean diameter of dominant trees) on mean expected growth.The blue background is the two-dimensional density of the observed age and diameter values.Dark blue contains more data.

Figure 6 .
Figure 6.Effects of nitrogen deposition on basal area increment, per tree species.Only diameter and nitrogen deposition effects are shown, all other explanatory variables are held constant (mean).

Figure 7 .
Figure 7. Interacting effects of basal area total (bat) and total stem number per hectare (nt).The blue in the background is the two-dimensional density of the observed bat and nt.Dark blue contains more data.

Figure 8 .Figure 9 .
Figure 8. Predicted basal area increment for interacting effects of temperature (t) and precipitation (p).t and p are mean values over physiological years (see Abbreviations).The blue in the background is the two-dimensional density of the observed t and p. Dark blue contains more data.

Table 2 .
Selected models and estimated coefficients ("Co.").A blank space means that the variable was not part of the final model.The standard errors of the coefficients can be found in the Supplement "Supp1.docx".
single tree-specfic) tree diameter (at 1.3 m height, breast height) in cm ddom dominant diameter in cm: quadratic mean of the 100 thickest trees per hectare A total plot area in hectares (10 4 m 2 )) t i −t i−1 rbairelative basal area increment: bai/ba cba cumulative basal area: cba j = 1 bat ∑ n i=1 (1 d i ≤d j ba i ) bal basal area of larger trees: bal j = 1 A ∑ n i=1 (1 d i >d j ba i ) rbal relative basal area of larger trees: bal/bat bah percentage of basal area previously harvested in relation to bat nt total number of trees per hectare in a plot nh percentage of stems harvested since last inventory cs crown surface area in m 2 : π r 6 l2 (4 l2 + r2 ) 3/2 − r3 with r the imputed crown radius, l the imputed crown length cb crown base (height where the living crown starts), in m csa crown cross sectional area in m 2 sdi stand density index: nt(dq/25) 1.605 cc crown competition index Equation (2) ccc change in crown competition index after harvesting has taken place Equation (2) lead.conleading tree species is coniferous (holding maximum basal area in a stand) year.past1940if the observation year is before 1940, it is zero, otherwise the current year of observation−1940.p yearly mean precipitation sum over period length in mm.Calculated for physiological year, meaning that the sum was taken from October until September, in contrast to calender years [42].t mean physiological year temperature, over intervals such as p. m Aridity index (defined by deMartonne ([38], [p.520])).A =