Mixing Effects in Norway Spruce — European Beech Stands Are Modulated by Site Quality , Stand Age and Moisture Availability

Although mixing tree species is considered an efficient risk-reduction strategy in the face of climate change, the conditions where mixtures are more productive than monocultures are under ongoing debate. Generalizations have been difficult because of the variety of methods used and due to contradictory findings regarding the effects of the species investigated, mixing proportions, and many site and stand conditions. Using data from 960 plots of the Swiss National Forest Inventory data, we assessed whether Picea abies (L.) Karst–Fagus sylvatica L. mixed stands are more productive than pure stands, and whether the mixing effect depends on siteor stand-characteristics. The species proportions were estimated using species proportion by area, which depends on the maximum stand basal area of an unmanaged stand (BAmax). Four different alternatives were used to estimate BAmax and to investigate the effect of these differing alternatives on the estimated mixture effect. On average, the mixture had a negative effect on the growth of Picea abies. However, this effect decreased as moisture availability increased. Fagus sylvatica grew better in mixtures and this effect increased with site quality. A significant interaction between species proportions and quadratic mean diameter, a proxy for stand age, was found for both species: the older the stand, the better the growth of Fagus sylvatica and the lower the growth of Picea abies. Overyielding was predicted for 80% of the investigated sites. The alternative to estimate BAmax weakly modulated the estimated mixture effect, but it did not affect the way mixing effects changed with site characteristics.


Introduction
In recent years, functions and ecosystem services provided by forests have been subject to increasing research efforts [1].Tree species diversity is correlated with many forest functions [2].Mixed forests are often considered more resistant and resilient to biotic and abiotic threats [3][4][5][6].In the face of changing climate and increasing natural disturbances, promoting mixed-species stands is considered an efficient strategy to ensure adaptation and to mitigate against the expected risks [7,8].The growth and yield of mixtures is therefore of great interest to forest managers and scientists considering the increasing needs for wood supply and carbon storage [4].
While mixtures are often more productive than monocultures [9], there are also many examples where they are not [10,11].This makes generalizations difficult because mixing effects depend on the mixture investigated [12], the site quality [13,14], the silvicultural treatment [15], the age of the stand [16], the climate [17], and the proportions of the different species [18].The productivity of mixed stands is often described in terms of "overyielding" or "underyielding".Overyielding (or underyielding) occurs when the productivity of the mixed stand is higher (or lower) than the mean productivity of the pure stands.In some cases, the mixed stand is even more productive than a pure stand of the most productive species within the mixture; this situation is called "transgressive overyielding" [15].
In Central Europe, the most investigated mixture in terms of growth is that of Norway spruce (Picea abies (L.) Karst.) and European beech (Fagus sylvatica L.) [13,19].The different studies on this mixture have shown a variety of mixing effects, depending on the study and the site conditions.In references [12,13], the authors found cases ranging from transgressive overyielding to underyielding depending on the site quality, with less frequent overyielding at high quality sites.A temporal variation in the mixing effect was also found [20].In Switzerland, Norway spruce is the most common species, such that 18% of the Swiss National Forest Inventory (NFI) plots are covered by pure Norway spruce stands and the mixture of Norway spruce and European beech is the second most common two-species mixture type [17].The climatic and edaphic conditions of these sites cover wide ranges due to the variable topography; however, stand level studies on the mixing effect on these species have never been carried out in this country.
In the literature, many different approaches have been used to calculate species proportions.These include simple approaches such as proportions in terms of stem number, basal area, stem volume, and aboveground biomass, to more elaborate indices such as the relative density index [21], basal area weighted by the species-specific wood density [22] or proportions by area occupied by the tree [20].However, the mixture effect has been shown to depend on the method employed to define species proportion [11,23] and inappropriate approaches can lead to biased predictions of over-or under-yielding [16].In [23], the authors suggested to use proportions in terms of leaf area, canopy surface area, leaf biomass or root biomass.Since these variables are often unavailable, especially when considering inventory data, many authors suggested the use of species proportion by area as a proxy for the area and resources available for each species [23,24].Nevertheless, this method requires the calculation of the maximum basal area BA max that would occur in a fully stocked pure stand on the given site.In the literature, BA max has been estimated using Reineke's self-thinning rule [25] and the extended competition density rule [26].In [27], the authors showed how the resulting BA max differs between these two methods.However, it is unknown how these differences influence the estimated mixing effect, i.e., how sensitive the results in terms of estimated mixing effects are to the method used to estimate BA max .
Given the uncertainties remaining on the mixture effect in Norway spruce and European beech mixtures, and the methodological aspects of the calculation of species proportion by area, the main aims of the current study were: 1.
To compare the productivity of mixed and pure stands of Norway spruce and European beech along Switzerland's large environmental gradients and to identify possible interactions between mixing effects and site, climate, age, and stand density.For this purpose, we focused on two hypotheses: (a) On average, mixed stands are more productive than pure stands, and (b) the positive effect decreases as the site conditions become more favorable for growth.2.
To investigate whether mixing effects differ depending on the method employed to estimate BA max (required for the calculation of species proportion by area).We hypothesized that (c) the methods employed lead to a different estimation of the mixing effect and the resulting overor under-yielding.
For this study, we used data from the Swiss NFI, which provide an unbiased database derived from a systematic sampling across Switzerland.They are well suited to study large scale mixing effects because they include repeated measures over time, and cover a large area including a wide range of environmental and climatic gradients [16,28,29].

Materials and Methods
As this study has both an ecological as well as a methodological focus, the materials and methods section is rather long.To enhance readability, we first briefly explain the content of the following subsections and how they are connected to the investigations outlined in the introduction.Section 2.1 contains details about the Swiss NFI data and which plots were selected for the study.The methods used to derive species proportions by area are explained in Section 2.2.For the purpose of investigating whether mixing effects differ depending on the method to estimate BA max , four different alternatives were used to estimate the species-specific BA max occurring in pure stands.The four alternatives were derived from two main rules: the self-thinning rule (Section 2.2.1) and the competition density rule (Section 2.2.2).A modelling approach was used to investigate the productivity of mixed and pure stands: Section 2.3 contains the statistical approach to investigate mixture effects at the species level, and Section 2.4 describes how the statistical models at the species level are then applied to investigate if the mixture effects lead to over-or under-yielding at the total stand level.Finally, Section 2.5 explains how it was investigated if the estimated over-or under-yielding at the stand level differs depending on the four alternatives to estimate BA max .

Swiss NFI Sampling Design and Plot Selection
We used data from four Swiss NFI surveys collected between 1983 and 2016.The sample plots were established during the first NFI (NFI1;1983-1985) according to a systematic 1 km square grid across Switzerland.For the second (1993-1995), the third (2004)(2005)(2006) and the fourth NFI (2009-2017; in progress), the surveys were carried out on a subsample of the NFI1 grid.Half of the plots were resampled under a √ 2 km square grid diagonal to the one used for the NFI1.The trees' characteristics (e.g., diameter at breast height DBH, tree species) were recorded within two concentric circles of 200 m 2 and 500 m 2 where the sample trees included all the trees with 12 ≤ DBH < 36 cm and DBH ≥ 36 cm, respectively.The height was measured for a subsample of the sample trees.Missing tree heights were estimated using site-specific models [30].Stand characteristics (e.g., stand structure, stand damages) were assessed by the field crew within a 2500 m 2 square around the plot center.The detailed methods of survey and analysis of the Swiss NFI data can be found in [31][32][33] and on the NFI-website [34].
To define species proportion, we estimated the species-specific maximum density N max [25] and the potential maximum basal area BA max [26], which is the highest basal area achievable by pure stands at a given site.For this purpose, we used sample plots from NFI1 fulfilling the following criteria: (1) the basal area (BA) of the investigated species exceeded 90% of the total stand BA; (2) the stand had an even-aged structure with high or normal crown closure; (3) the stand was not or only slightly damaged (grade based on the weighted mean of the damages assessed on the individual trees); (4) the last silvicultural treatment occurred more than 10 years ago; (5) the plot contained at least 4 trees; and (6) more than 50% of the sample circle was located inside the forest.The selection resulted in 550 plots for Norway spruce and 158 for European beech (see Table S1).
For investigating tree growth, we used data from all four NFI periods.A stand was considered as pure when the BA of the investigated species accounted for more than 90% of the total BA.For mixed stands, the following criteria were chosen: the BA of Norway spruce and European beech together exceeded 90% of the total BA and each investigated species alone represented at least 10% of the BA.Furthermore, we restricted the plots to those showing an even-aged stand structure, high or normal crown closure, no or only slight damage (see paragraph above), and more than 75% of the sample circle located inside the forest.In addition, we excluded plots where the stand BA decreased by more than 5% between two surveys to avoid stands where severe silvicultural treatments or natural disturbances had occurred.Table 1 shows the main stand and site variables of the selected plots.

Species Proportions by Area
To calculate species proportions in mixed stands, we accounted for species-specific density and growing space requirements [15].When spatial information at the tree level is not available, species proportion by area is often used to assess the area occupied by each species at the stand level [10,23]: BA i BA maxi (1) where P i is the proportion by area of species i, BA i is the observed basal area of species i (m 2 •ha −1 ), BA max i is the maximum basal area of species i (m 2 •ha −1 ) and n is the number of species in the stand.
The species-specific maximum basal area BA max is the basal area of a fully stocked pure stand; a species with a high BA max has low growing space requirements.P i can be considered as the proportion of the area and resources available for the species i [35].We estimated BA max for Norway spruce and European beech using four different alternatives based on the two methods described below.

Self-Thinning Rule
Reineke's self-thinning rule (STR) relates the stem number to the quadratic mean diameter of the stand ( [25]; Equation (2)).In logarithmic scales, the relationship between quadratic mean diameter and stem density becomes linear (Equation ( 3)).The concept of STR was developed for fully stocked pure stands.For stands that are not at their stocking capacity, the self-thinning boundary line can be estimated using quantile regression [36].The coefficients were estimated for the 95th percentile using the function rq from the package {quantreg} [37] in the statistical software R [38].BA max was then calculated using Equation (4).
where N max is the species-specific maximum stem number (ha −1 ), Dg is the quadratic mean diameter (cm), BA max is the species-specific maximum basal area (m 2 •ha −1 ), and C and E are the coefficients of the STR.Without the need to apply allometric relationships derived from other models, the use of Dg has been shown to be more appropriate than mean stem volume or mean stand height [39].
It is generally assumed that the parameters of the STR vary between species [27,40], but in the literature different parameters for the same species also exist [21,27,41,42].In order to assess how sensitive mixing effects are to the parameters of the STR, we did not use only our own parameters but alternatively recalculated BA max (Equation (4)) using the value for C and E described in Charru et al. [43].

Extended Competition Density Rule
The competition density rule (CDR) is another principle widely used in plant ecology ( [44]; Equation ( 5)).For a given dominant height, Dg decreases as the density (N) increases.Sterba [26] showed that BA max can be obtained from the CDR and introduced the CDR into the expression of BA showing that when we derive this "extended CDR" ∂BA/∂N and set the first derivative to 0, we obtain an expression of BA max dependent only on H dom (Equation ( 6)).However, this 4-parameter CDR can only be fitted using large datasets [27].Sterba [26] also found that the coefficients of Reineke's STR and the CDR are correlated (Equation ( 7)), which reduces the number of parameters to be estimated.Thus, alternatively both or only one of the two coefficients of the STR can be introduced into the CDR resulting in a 3-parameter (Equation ( 8)) or a 2-parameter extended CDR (Equation ( 9)).
where Dg is the quadratic mean diameter (cm), N is the stem number per hectare (ha −1 ), H dom is the dominant height (m), BA max is the species-specific maximum basal area (m 2 •ha −1 ), a 0 , a 1 , b 0 , b 1 are the coefficients of the CDR, and C and E are the coefficients of the STR (Equation ( 2)).The dominant height was calculated according to Assmann [45].The 3-parameter and 2-parameter CDR were fitted for each species using nonlinear regression and the coefficients of the STR that were obtained with the NFI dataset.We used the function nls from the R package {stats} [38].BA max was then calculated for each species and each alternative using Equation (6).
To enhance readability, the different alternatives are henceforth called "3P" for the 3-parameter extended CDR, "2P" for the 2-parameter extended CDR, "Direct-SNFI" for the alternative where BA max is calculated directly from the STR using coefficients derived with NFI data (SNFI stems for Swiss National Forest Inventory), and "Direct-Charru" for the alternative where BA max is calculated directly from the STR but using coefficients by Charru et al. in France [43] who reported the highest slopes for both species compared to other studies.Table 2 summarizes the four different alternatives and their BA max formulae.The 3P alternative will be the reference alternative for the selection of the growth models.

Alternative
BA max Formula STR-Coefficients Used C and E from Swiss NFI data Direct-SNFI C and E from Charru et al. [43] BA max : maximum stand basal area of an unmanaged stand; SNFI: Swiss National Forest Inventory; STR: self-thinning rule; C and E are the coefficients of the STR.3P: 3-parameter extended competition density rule; 2P: 2-parameter extended competition density rule; Direct-SNFI: BA max is calculated directly using STR coefficients derived with SNFI data; Direct-Charru: BA max is calculated directly using STR coefficients by Charru et al. [43].

Modelling the Stand Basal Area Increment
Growth was quantified as the annual BA increment (BAI) of the survivor trees, i.e., the trees that were alive at two consecutive inventories.The BA increment between two surveys was divided by the number of growing seasons and the species proportion: where iG i is the annual increment per unit species proportion of species i (m 2 •ha −1 •year −1 ), aBAI i is the average annual basal area increment of species i (m 2 •ha −1 •year −1 ), and P i is the proportion of species i in the stand (according to Equation ( 1)).
The growth was assumed to depend on climate, site, stand characteristics, and species proportions.We used the following model as a basis for the growth investigations: where iG i is the annual growth increment per unit species proportion of species i, P i is the proportion of species i (Equation ( 1)), ∑ m j=1 f ij × v j × P i represents the interactions between P i and the m potential explanatory variables v j included in the model, and a i , b i , c i , d i , e i and f ij are model coefficients to be estimated.These interactions were added in the growth models in order to assess whether the mixing effect depends on climate-, site-or stand-conditions.

Climate Variables
Daily data for precipitation sum, moisture index (ratio between actual and potential evapotranspiration ETa/ETp; the lower the dryer) and minimum, mean and maximum temperature were obtained for each NFI plot [46].The growing season was defined as the period when the monthly mean of the daily mean temperature exceeded 5 • C and the monthly mean of the daily minimum exceeded 0 • C [11].The effect of climate was investigated through the following expression: where t is the sum of the daily mean temperature above 3 • C per growing season (degree days•year −1 ), p is the precipitation sum per growing season (mm•year −1 ), and mi is the mean moisture index per growing season (%).The pairwise Pearson correlations among the explanatory variables did not exceed |0.42|.The t 2 and p 2 terms were introduced into the model in order to test for a possible maximum or minimum growth below the upper range values of temperature and precipitation.

Site Variables
The NFI does not assess any site index or edaphic variables (except pH) that could be used as a proxy for site quality.In Switzerland, a stand growth capacity model was used as a proxy for site quality while studying the mixing effect on the growth of Norway spruce and Abies alba Mill. in a previous study using NFI data [11].This model takes into account, amongst other variables, elevation, azimuth, and slope [47].As we were interested in possible interactions between mixing effects and those site variables, we refrained from using this model.
For this reason, we developed an alternative proxy for site quality based on the method by Trasobares and Pukkala [48] in a slightly modified form.This method is based on predicting the growth of a tree on an average site and setting this in relation to its observed growth [48].The same kind of modelling has been used as a productivity index at the stand level [49] or by using modified residuals of the individual tree growth model [50].To predict tree growth on an average site, we modelled the annual past radial increment (APRI) of a tree of a certain species as a function of species-specific competition, tree size, and climate but omitting site variables (Equation ( 13)).BAL sp was added to the model to quantify intraspecific competition.Combined with BAL already included in the model, we explicitly accounted for the possibility of considering both inter-and intra-specific competition.The main variables used for the individual-tree growth modelling are shown in Table S2.Due to possible inaccuracy in DBH measurement, observed values of APRI could take on negative values, which is problematic when APRI should be used to derive a growth index (see below).We thus added a constant to APRI in the individual-tree growth models to avoid negative values.The minimum APRI was −12.24 for Norway spruce and −8.28 for European beech (Table S2); we therefore added constants equal to 13 and 9 for Norway spruce and European beech, respectively.Since the trees were repeatedly measured, we used linear mixed effect models with a tree random intercept.The models were fitted in R [38] using the function lme from the package {nlme} [51].
where APRI i is the tree annual past radial increment of species i (mm•year −1 ), K i is a constant (13 for Norway spruce and 9 for European beech), DBH is the diameter at breast height of the tree (cm), BA is the stand basal area (m 2 •ha −1 ), BAL is the basal area of the trees larger than the subject tree (m 2 •ha −1 ), BAL sp is the basal area of the trees of the same species and larger than the subject tree (m 2 •ha −1 ), t is the sum of the daily mean temperature above 3 • C during the growing season (degree days•year −1 ), p is the sum of the precipitation during the growing season (mm•year −1 ), mi is the mean daily moisture index during the growing season (%), β i0 − β i8 are the model coefficients, z t is the random intercept and ε is the residual error term.
Finally, the stand growth index was calculated for each investigated species as the mean ratio per plot of the observed growth of an individual tree and its model-predicted growth: where GI i is the stand growth index of species i, n i is the number of trees of species i in the plot, APRI ij is the observed annual past radial increment of the tree j of species i (mm•year −1 ), and AP " RI ij is the predicted past radial increment of the tree j of species i using the individual-tree model (Equation ( 13)) (mm•year −1 ).GI, according to Equation (14), is based on the relative residuals of the growth model, thus it accounts for all the factors influencing growth that are not included in the model.Since we used tree size, climate and competition measures as explanatory variables in the model, GI can be seen as a proxy for site properties such as soil and topography, which were not included in the model.A GI value of 1 stands for average growth while large GI values (>1) indicate a better than average growth (i.e., a good site quality), and small GI values (close to 0) indicate a poor site quality.
The effect of site was included in the initial model (Equation ( 11)) as follows: where GI i is the growth index of species i.

Stand Variables
Quadratic mean diameter was used as a proxy for stand age [11].The level of competition was estimated using the stocking degree [42]: where sd is the stocking degree, BA i is the basal area of species i (m 2 •ha −1 ), BA maxi is the species-specific maximum basal area (m 2 •ha −1 ) of species i, and n is the number of species in the stand.Based on our plot selection criteria, our stands can carry a maximum of 10% of other species.No BA max values were available for these minor species, so we used BA max of Norway spruce for the other coniferous species and BA max of European beech for the other broadleaved species.Stand characteristics were added in the model as follows: where Dg is the quadratic mean diameter (cm), and sd is the stocking degree (Equation ( 16)).

Model Fitting
We used linear mixed-effect models with a random intercept and the sample plot as the grouping factor [11,52] to account for the repeated observations on the same plots.In order to exclude the potential outliers and obtain a more robust estimation, we excluded 1% of the plots that were located at the large and small end of the annual BAI distribution, respectively.The initial regression model was the following: where iG i is the annual basal area increment per species proportion of species i (m 2 •ha −1 •year −1 ), t is the sum of the daily mean temperature above 3 • C per growing season (degree days•year −1 ), p is the precipitation sum per growing season (mm•year −1 ), mi is the mean moisture index per growing season (%), GI i is the growth index of species i, Dg is the quadratic mean diameter (cm), sd is the stocking degree, P i is the proportion of investigated species i (Equation ( 1)), ∑ m j=1 f ij × v j × P i represents the interactions between P i and the m potential explanatory variables v j included in the model, z p is the random intercept at plot level, and ε is the residual error term.
Variable selection was performed separately for Norway spruce and European beech following a 3-step procedure: (1) we fitted Equation (18) excluding species proportions and interactions, with ordinary least squares regression (OLS), and used a stepwise-backward approach to exclude non-significant variables (p-value > 0.05); (2) the pre-selected variables from step 1, the species proportions and the interactions between species proportions and the pre-selected variables from step 1 were added in the model, and the non-significant interactions were excluded stepwise-backward using OLS; and (3) the resulting model from step 2 was fitted using the lme function from the package {nlme} [51] in R [38] and the non-significant variables were excluded step by step.

Comparison of the Productivity of Pure and Mixed Stands
Two stand basal area increment models (Figure 1) were constructed using the fixed effects of the best growth models described in Section 2.3: G mixed , the predicted growth of a mixed stand and G 0 , the predicted growth of a mixed stand under the null hypothesis (i.e., considering that no mixing effects occur).G 0 was the sum of the predicted growth of pure stands (in the selected growth models, the proportions were set to 1) which were weighted by their share in the hypothetical mixed stand (P Ns and 1 − P Ns for Norway spruce and European beech, respectively).G mixed was obtained by adding the predicted growth of both species in mixed stands (in the growth models, the proportions were set to P Ns for Norway spruce and 1 − P Ns for European beech) which were weighted by their share in the mixture P Ns and 1 − P Ns .For given species proportions, G mixed − G 0 > 0 indicate that the growth of the mixed stand exceeds the growth of pure stands (i.e., overyielding).The reverse situation implies underyielding.Transgressive overyielding can occur when G mixed > GP Ns and G mixed > GP Eb , i.e., the mixed stand is even more productive than a pure stand of the most productive species [15].
These stand basal area increment models were applied using hypothetical 50%-50% proportions for all the 192 NFI plots where Norway spruce and European beech occur, to avoid extrapolations.We selected 50%-50% mixture proportions because on average, the patterns we highlighted were the strongest.

Comparison of the Productivity of Pure and Mixed Stands
Two stand basal area increment models (Figure 1) were constructed using the fixed effects of the best growth models described in Section 2.3: Gmixed, the predicted growth of a mixed stand and G0, the predicted growth of a mixed stand under the null hypothesis (i.e., considering that no mixing effects occur).G0 was the sum of the predicted growth of pure stands (in the selected growth models, the proportions were set to 1) which were weighted by their share in the hypothetical mixed stand (PNs and 1 − PNs for Norway spruce and European beech, respectively).Gmixed was obtained by adding the predicted growth of both species in mixed stands (in the growth models, the proportions were set to PNs for Norway spruce and 1 − PNs for European beech) which were weighted by their share in the mixture PNs and 1 − PNs.For given species proportions, Gmixed − G0 > 0 indicate that the growth of the mixed stand exceeds the growth of pure stands (i.e., overyielding).The reverse situation implies underyielding.Transgressive overyielding can occur when Gmixed > GPNs and Gmixed > GPEb, i.e., the mixed stand is even more productive than a pure stand of the most productive species [15].
These stand basal area increment models were applied using hypothetical 50%-50% proportions for all the 192 NFI plots where Norway spruce and European beech occur, to avoid extrapolations.We selected 50%-50% mixture proportions because on average, the patterns we highlighted were the strongest.Comparing the productivity of pure and mixed stands.Gmixed, the predicted growth of a mixed stand and G0, the predicted growth of a mixed stand under the null hypothesis (i.e., considering that no mixing effects occur).

Comparison of Mixing Effects Depending on the BAmax Alternative
The growth increment (Equation ( 10)) was recalculated using 2P, Direct-SNFI, and Direct-Charru proportions (Table 2).The final models for Norway spruce and European beech, obtained as described in Section 2.3 using the reference alternative 3P, were refitted for each of these growth increments with the corresponding proportions.The resulting under-or overyielding with 50%-50% proportions were then estimated as described in Section 2.4 for the same sites and the results of the four alternatives were compared using pairwise t-tests.

Calculation of the Maximum Basal Area and the Growth Index GI
The estimates of the coefficients of the STR and their standard error are shown in Table 3.The coefficients obtained using the NFI dataset considerably differed from the ones found by Charru et al. [43] resulting in a higher range of BAmax when calculated using the Direct-SNFI alternative (Figure 2).The alternatives 3P and 2P do not show such large differences between each other and to the Direct-SNFI alternative (see Table S3 and Figure S1).Comparing the productivity of pure and mixed stands.G mixed , the predicted growth of a mixed stand and G 0 , the predicted growth of a mixed stand under the null hypothesis (i.e., considering that no mixing effects occur).

Comparison of Mixing Effects Depending on the BA max Alternative
The growth increment (Equation ( 10)) was recalculated using 2P, Direct-SNFI, and Direct-Charru proportions (Table 2).The final models for Norway spruce and European beech, obtained as described in Section 2.3 using the reference alternative 3P, were refitted for each of these growth increments with the corresponding proportions.The resulting under-or overyielding with 50%-50% proportions were then estimated as described in Section 2.4 for the same sites and the results of the four alternatives were compared using pairwise t-tests.

Calculation of the Maximum Basal Area and the Growth Index GI
The estimates of the coefficients of the STR and their standard error are shown in Table 3.The coefficients obtained using the NFI dataset considerably differed from the ones found by Charru et al. [43] resulting in a higher range of BA max when calculated using the Direct-SNFI alternative (Figure 2).The alternatives 3P and 2P do not show such large differences between each other and to the Direct-SNFI alternative (see Table S3 and Figure S1).Table 3. Coefficients of the self-thinning rule estimated using quantile regression (τ = 0.95) with the NFI data and coefficients found by Charru et al. [43] using stochastic frontier analysis.(a) (b) With regard to the growth index GI, all variables of the annual past radial growth increment models were significant (p-value < 0.05) except ln(DBH) for European beech that has been removed.The estimates of the parameters as well as the sample size and goodness of fit statistics can be found in Tables S4 and S5.For Norway spruce and European beech, the mean GI values were 1.000 and 1.004 with standard deviation of 0.030 and 0.052, respectively (see Table S6).

Stand Basal Area Increment Modelling
Table 4 shows the estimated fixed-effect coefficients of the best-fitting growth models for the investigated species.The corresponding sample sizes and goodness of fit statistics can be found in Table S7).The distributions of standardized residuals are shown in Figure S2 for Norway spruce and Figure S3 for European beech; no heteroscedasticity was detected.
Both species' growth was influenced by species proportions.The growth of Norway spruce increased the higher its own proportions in the stand, i.e., Norway spruce grew better in pure stands than in mixed stands.A significant interaction between the quadratic mean diameter (Dg) and species proportion was found.The slope of the regression line (Figure 3a) was negative (i.e., positive mixing effect) for stands with small Dg (i.e., young stands) whereas the mixing effect was found to be positive and increasing with higher Dg.We also found a significant interaction between species proportions and the moisture index (mi).The negative mixing effect (positive slopes) decreased slightly as the moisture index increased (Figure 3b), meaning that the negative mixing effect was lower at wet sites than at dry sites.
The growth of European beech decreased when its own proportions increased, denoting With regard to the growth index GI, all variables of the annual past radial growth increment models were significant (p-value < 0.05) except ln(DBH) for European beech that has been removed.The estimates of the parameters as well as the sample size and goodness of fit statistics can be found in Tables S4 and S5.For Norway spruce and European beech, the mean GI values were 1.000 and 1.004 with standard deviation of 0.030 and 0.052, respectively (see Table S6).

Stand Basal Area Increment Modelling
Table 4 shows the estimated fixed-effect coefficients of the best-fitting growth models for the investigated species.The corresponding sample sizes and goodness of fit statistics can be found in Table S7).The distributions of standardized residuals are shown in Figure S2 for Norway spruce and Figure S3 for European beech; no heteroscedasticity was detected.
Both species' growth was influenced by species proportions.The growth of Norway spruce increased the higher its own proportions in the stand, i.e., Norway spruce grew better in pure stands than in mixed stands.A significant interaction between the quadratic mean diameter (Dg) and species proportion was found.The slope of the regression line (Figure 3a) was negative (i.e., positive mixing effect) for stands with small Dg (i.e., young stands) whereas the mixing effect was found to be positive and increasing with higher Dg.We also found a significant interaction between species proportions and the moisture index (mi).The negative mixing effect (positive slopes) decreased slightly as the moisture index increased (Figure 3b), meaning that the negative mixing effect was lower at wet sites than at dry sites.effect was positive and increased as Dg increased (Figure 4a).Additionally, we found a significant interaction between species proportions and the growth index GI.The slope of the regression line (Figure 4b) was negative and slightly increased with increasing GI; this means that positive mixing effects were higher for stands with higher site quality.The growth of European beech decreased when its own proportions increased, denoting positive mixing effects on the growth of this species.As for Norway spruce, we found a significant interaction between the proportions and Dg, but in opposite directions.The mixture had a slight negative effect (positive slope) on the growth of European beech in young stands while the mixing effect was positive and increased as Dg increased (Figure 4a).Additionally, we found a significant interaction between species proportions and the growth index GI.The slope of the regression line (Figure 4b) was negative and slightly increased with increasing GI; this means that positive mixing effects were higher for stands with higher site quality.

Productivity in Pure and Mixed Stands
Due to the significant interactions with other factors, the mixture effects at the species level can lead to over-or under-yielding.The diagrams in Figure 5 show two examples of basal area increment per species and for the total stand at increasing proportions of beech (0: pure Norway spruce stand; 0.5: 50%-50% mixed Norway spruce-European beech stand; 1: pure European beech stand).The annual basal area increments were calculated for 3P proportions and with fixed climate and site conditions.Figure 5a shows the site with the lowest underyielding.This site presents a poor site quality for European beech (i.e., low growth index) and dry conditions (i.e., low moisture index).Growth of Norway spruce was negatively affected by the mixture while European beech was affected slightly positively.This led to an underyielding of −0.157 m 2 •ha −1 •year −1 for 50%-50% proportions.Figure 5b illustrates a case of transgressive overyielding (overyielding of 0.061 m 2 •ha −1 •year −1 for 50%-50% proportions) at medium site quality.The mixing effect was slightly negative for Norway spruce but highly positive for European beech.The total stand growth is even higher than the productivity of pure stands of both species; this is a situation of transgressive overyielding.

Productivity in Pure and Mixed Stands
Due to the significant interactions with other factors, the mixture effects at the species level can lead to over-or under-yielding.The diagrams in Figure 5 show two examples of basal area increment per species and for the total stand at increasing proportions of beech (0: pure Norway spruce stand; 0.5: 50%-50% mixed Norway spruce-European beech stand; 1: pure European beech stand).The annual basal area increments were calculated for 3P proportions and with fixed climate and site conditions.Figure 5a shows the site with the lowest underyielding.This site presents a poor site quality for European beech (i.e., low growth index) and dry conditions (i.e., low moisture index).Growth of Norway spruce was negatively affected by the mixture while European beech was affected slightly positively.This led to an underyielding of −0.157 m 2 •ha −1 •year −1 for 50%-50% proportions.Figure 5b illustrates a case of transgressive overyielding (overyielding of 0.061 m 2 •ha −1 •year −1 for 50%-50% proportions) at medium site quality.The mixing effect was slightly negative for Norway spruce but highly positive for European beech.The total stand growth is even higher than the productivity of pure stands of both species; this is a situation of transgressive overyielding.Underyielding (or overyielding) is the negative (or positive) difference between the predicted basal area increment in mixture (continuous line) and the expected increment under the assumption of no mixing effects (dashed line).GINs: growth index for Norway spruce.
Overyielding was found on 153 plots out of the 192 NFI mixed plots (79.7%); the mean overyielding was 0.012 m 2 •ha −1 •year −1 (Figure 6).In terms of relative gain or loss compared to pure stands, the productivity of mixed stands ranged between −15% and 22% with a mean positive effect of 2%.Transgressive overyielding occurred on 13 plots of the 192 NFI mixed plots.Underyielding (or overyielding) is the negative (or positive) difference between the predicted basal area increment in mixture (continuous line) and the expected increment under the assumption of no mixing effects (dashed line).GI Ns : growth index for Norway spruce.
Overyielding was found on 153 plots out of the 192 NFI mixed plots (79.7%); the mean overyielding was 0.012 m 2 •ha −1 •year −1 (Figure 6).In terms of relative gain or loss compared to pure stands, the productivity of mixed stands ranged between −15% and 22% with a mean positive effect of 2%.Transgressive overyielding occurred on 13 plots of the 192 NFI mixed plots.

Comparison of Different Alternatives to Estimate BAmax
The best growth models obtained using the reference alternative 3P (Table 4) were refitted for the other alternatives.The estimated fixed-effect parameters, and the goodness of fit statistics are shown in Table S8.The predicted growth increments were recalculated for each of the 192 NFI mixed plots using 50%-50% proportions and the distributions of the resulting under-or over-yielding are shown in Figure 7.The different alternatives showed similar distributions except Direct-SNFI that gave significantly different results (p-value < 0.05).Diagrams of site productivity for the site showing the lowest underyielding, the mean overyielding and a case of transgressive overyielding are shown in Figures S4-S6, respectively.The mixing effects obtained with the Direct-SNFI alternative differed from the others: the negative mixing effect was stronger for Norway spruce and the positive effect was slightly stronger for European beech, resulting in a lower overall productivity when we aggregated the productivity of both species at the total stand level.

Comparison of Different Alternatives to Estimate BA max
The best growth models obtained using the reference alternative 3P (Table 4) were refitted for the other alternatives.The estimated fixed-effect parameters, and the goodness of fit statistics are shown in Table S8.The predicted growth increments were recalculated for each of the 192 NFI mixed plots using 50%-50% proportions and the distributions of the resulting under-or over-yielding are shown in Figure 7.The different alternatives showed similar distributions except Direct-SNFI that gave significantly different results (p-value < 0.05).Diagrams of site productivity for the site showing the lowest underyielding, the mean overyielding and a case of transgressive overyielding are shown in Figures S4-S6, respectively.The mixing effects obtained with the Direct-SNFI alternative differed from the others: the negative mixing effect was stronger for Norway spruce and the positive effect was slightly stronger for European beech, resulting in a lower overall productivity when we aggregated the productivity of both species at the total stand level.

Comparison of Different Alternatives to Estimate BAmax
The best growth models obtained using the reference alternative 3P (Table 4) were refitted for the other alternatives.The estimated fixed-effect parameters, and the goodness of fit statistics are shown in Table S8.The predicted growth increments were recalculated for each of the 192 NFI mixed plots using 50%-50% proportions and the distributions of the resulting under-or over-yielding are shown in Figure 7.The different alternatives showed similar distributions except Direct-SNFI that gave significantly different results (p-value < 0.05).Diagrams of site productivity for the site showing the lowest underyielding, the mean overyielding and a case of transgressive overyielding are shown in Figures S4-S6, respectively.The mixing effects obtained with the Direct-SNFI alternative differed from the others: the negative mixing effect was stronger for Norway spruce and the positive effect was slightly stronger for European beech, resulting in a lower overall productivity when we aggregated the productivity of both species at the total stand level.

Mixing Effects in Norway Spruce-European Beech Stands
Using NFI data, we found that European beech benefitted from the mixture.This finding is comparable to the results of previous studies [12,13,42].European beech exhibits a low self-tolerance and it generally benefits from the admixture of Norway spruce due to a release from severe intraspecific competition [53,54].A previous study showed that European beech can develop larger crowns in mixture than in pure stands [55]; this can facilitate light interception and can partly explain the positive mixing effect that we found for this species.In our study, the mixing effect on the productivity of European beech increased with increasing site quality (Figure 3b).Another study on the same species mixture in central Europe [13] found a similar trend suggesting that the intraspecific competition is higher at stands with high site quality and thus European beech benefits more from the admixture of Norway spruce on these sites.The mixing effect on European beech growth also depended on the quadratic mean diameter, which is a proxy for the age of the stand.At very young development stages, the growth of European beech suffered from the mixture but at medium old development stages, European beech benefitted from the mixture.Similarly, the mixing effect for European beech with Pinus sylvestris (L.) was predicted to increase with stand age as the early dominance of Pinus sylvestris, in terms of height, was lost [56].The mixture therefore not only affects the productivity of the mixed stands but also the temporal dynamics of the tree species interactions [19].
The admixture of European beech led to a decrease in the growth of Norway spruce.However, this negative mixing effect diminished with increasing moisture index, i.e., with increasing water availability.This pattern may result because European beech tends to have a high above-and below-ground competitiveness [13,57,58].In this mixture type, Norway spruce was shown to exhibit a shallower root system [59] whereas European beech appeared to develop longer and more fine roots and to extend the space exploited to soil layers less occupied by competitors [60].This interaction is likely to be less harmful for Norway spruce at sites where water resources are not strongly limiting growth [61], and this is consistent with our results.Our hypothesis (b) is therefore not supported.As for European beech, the mixing effects on Norway spruce growth depended on the proxy for stand age, but in the opposite direction: a positive mixture effect occurred on young stands whereas the mixture had a negative effect on medium and old stands.
Regarding the estimated gain or loss of productivity, compared to pure stands, a 50%-50% mixed stand productivity ranged between −15% and 22% and led to overyielding on 80% of the NFI mixed plots (mean effect of 2%).On average, the high performance of European beech in mixture [62] was greater than any loss in productivity of Norway spruce.These results are in line with the stated hypothesis (a).However, our predicted overyielding is considerably smaller than previous studies on this mixture type [12,13].From a forester point of view, the results of the present study might not be relevant in terms of wood production especially on highly productive sites for Norway spruce.Nevertheless, it should not veil the interactions at the species level and the benefits in terms of biodiversity, resistance and resilience that reduce financial risks [4].Concerning carbon storage, below-ground carbon stocks have been shown to be influenced by species composition [63,64].A previous study demonstrated a complementary effect occurring on soil organic carbon stocks in Norway spruce-European beech mixtures [64].Norway spruce litter increased the stock in the forest floor (low degradation rates) whereas the root turnover of European beech increased soil organic carbon stock in the upper mineral soil layer.

Influence of Methodology on the Estimation of Mixing Effects
The different alternatives to define the species proportion resulted in similar predictions for the growth of each species and for the overall productivity (Figure 7) except Direct-SNFI.The estimated productivity of the mixed stand depended on the alternative used to estimate BA max and this result partly confirms our hypothesis (c).However, even when the mean of the overall predictions was significantly different between Direct-SNFI and the other alternatives, the choice of the BA max approach appeared to be of minor relevance due to the similarity of the main finding: the mixture led to overyielding more often than underyielding.
Whether one alternative should be preferred to another is a question that cannot be clearly answered.A previous study suggested that both CDR and STR can be used to estimate BA max [27].Concerning the Direct alternatives, only the quadratic mean diameter is required the dominant height H dom is necessary for 3P-and 2P-CDR.The Direct alternatives might be a reliable approach when H dom is not available in the dataset or if the H dom estimation is not highly robust.Moreover, H dom itself can be influenced by species mixture [65].Given the poor robustness of H dom estimation in the present study (see below), the Direct alternatives should be the most reliable.However, in this study, the Direct-SNFI alternative gave significantly different results.This could be explained by the fact that we found considerably low values for the STR slope and intercept leading to biased BA max estimation (see below).When using 3P-and 2P-alternatives, the underestimated STR coefficients (from the NFI dataset) might have been balanced by the CDR coefficients leading to similar predictions than Direct-Charru.As a conclusion, it seems relevant to avoid the Direct alternative when the STR coefficients are underestimated or to use values from the literature.When it is not the case, the Direct or the 3P-and 2P-CDR alternatives should be preferred depending on the availability and robustness of H dom .

Other Methodological Aspects
Using NFI data, we obtained relatively low slopes for the STR compared to the slopes found in the literature, in particular for European beech (Table 5).This might be due to the statistical design because most forests in Switzerland are managed, so we could not restrict our analysis to fully stocked pure stands as suggested in [25].Furthermore, even when considering the managed forests with the highest stocking, there were relatively few stands, so we were only able to fit a quantile regression to the 95th percentile and not the 99th percentile.In addition, the Swiss NFI uses a DBH-threshold of 12 cm.Stands with comparably small quadratic mean diameter and high stem density are therefore not represented in our data set, likely resulting in underestimated slopes of the STR.
Table 5.Values of the self-thinning rule slopes found using NFI data and found in the literature [11,21,[25][26][27][40][41][42][43]66].As one alternative, we used the coefficients of the study showing the highest slopes for both species [43] in order to conduct an analysis of sensitivity of the estimated BA max and mixing effects depending on the value of the slope.The STR parameters employed had a strong impact on the resulting BA max .The use of the NFI slopes led to higher BA max especially at high Dg (Figure 2).BA max might be overestimated as a consequence of using underestimated slopes of the STR.

Species
Additionally, tree height measurements were limited to a subsample of trees and the missing tree heights were estimated using site-specific equations [30].H dom was then calculated based on the heights of the 1 to 5 biggest trees on the plot (mean weighted by the representative stem number of the sample trees).Consequently, the estimation of H dom might not be robust because of the small sample size and the use of models to estimate the missing tree heights.Since H dom is a key variable for the modelling of CDR, the robustness of the 3P and 2P BA max models might be limited by the low robustness of the H dom estimation.
In order to use GI as a proxy for site quality, we included climate variables in the tree growth model so that the majority of the remaining variability would likely be associated with soil properties and topography.Furthermore, the basal area of the larger trees of the same species was included as a measure of intraspecific competition.The resulting GI was a significant predictor in our stand growth models and increased the goodness of fit.However, this approach assumes that site quality was evenly distributed over tree size, climate, and competition.As this is unlikely, some of the site effect might be removed when rating the observed growth by the predicted growth [48].Moreover, this index might include missing variables other than soil properties and topography.We thus refrained from calling it a "fertility index".

Conclusions
Using NFI data, we studied the annual basal area increment of mixed stands of European beech and Norway spruce at the stand level and compared these growth rates to those of pure stands.From a methodological perspective, we compared how the estimated mixing effect and over-or under-yielding depended on the method used to derive BA max .
We found that the growth of European beech was positively influenced by the admixture of Norway spruce whereas the growth of Norway spruce was lower in mixed stands compared to pure stands.Overall, the growth in mixed stands increased with site quality and moisture index and was modulated by stand age.A 50%-50% mixture was more productive than the respective pure stands on 80% of the sites where both European beech and Norway spruce occurred.The stand overyielding was due to the high performance of European beech in mixtures, which was greater than the loss of productivity for Norway spruce in mixtures but the overyielding was comparably weak (on average, 2%).
With the exception of Direct-SNFI, the different alternatives to define the species proportions resulted in similar predictions for the productivity of both species.All the alternatives appeared to be consistent while estimating BA max .The Direct alternative or 3-and 2-parameter CDR alternatives should be preferred depending on the dominant height availability and the value of the self-thinning rule coefficients.

Supplementary Materials:
The following are available online at www.mdpi.com/1999-4907/9/2/83/s1, Figure S1: Maximum basal area over quadratic mean diameter fitted using two alternatives (3P, 2P), Figure S2 S1: Main stand variables of the sample plots used to estimate Reineke's maximum stand density and maximum basal area, Table S2: Main variables for the modelling of the annual past radial growth increment, Table S3: Coefficients of the 3-and 2-parameter extended competition density rule estimated using nonlinear regression, Table S4: Fixed effect parameter estimates for the Norway spruce and European beech annual past radial increment models, Table S5: Sample size and goodness of fit statistics for the annual past radial increment models, Table S6: Statistics of the resulting GI, Table S7: Sample size and goodness of fit statistics for the growth models using 3P-SNFI proportions, and Table S8: Fixed effect parameter estimates for the growth models using different alternatives.

Figure 1 .
Figure 1.Comparing the productivity of pure and mixed stands.Gmixed, the predicted growth of a mixed stand and G0, the predicted growth of a mixed stand under the null hypothesis (i.e., considering that no mixing effects occur).

Figure 1 .
Figure 1.Comparing the productivity of pure and mixed stands.G mixed , the predicted growth of a mixed stand and G 0 , the predicted growth of a mixed stand under the null hypothesis (i.e., considering that no mixing effects occur).

Figure 2 .
Figure 2. Maximum basal area over quadratic mean diameter derived using two Direct alternatives: (a) for Norway spruce; (b) for European beech.The "SNFI" parameters were estimated using data of the Swiss National Forest Inventory (SNFI).The "Charru" parameters stem from Charru et al.[43].Dg: quadratic mean diameter (cm).

Figure 2 .
Figure 2. Maximum basal area over quadratic mean diameter derived using two Direct alternatives: (a) for Norway spruce; (b) for European beech.The "SNFI" parameters were estimated using data of the Swiss National Forest Inventory (SNFI).The "Charru" parameters stem from Charru et al.[43].Dg: quadratic mean diameter (cm).

Figure 3 .Figure 3 .
Figure 3.Effect of species proportion on Norway spruce basal area growth per unit species proportion (a) depending on the quadratic mean diameter and (b) depending on the moisture index.All the other explanatory variables are set to their mean values.mi: moisture index.

Figure 3 .
Figure 3.Effect of species proportion on Norway spruce basal area growth per unit species proportion (a) depending on the quadratic mean diameter and (b) depending on the moisture index.All the other explanatory variables are set to their mean values.mi: moisture index.

Figure 4 .
Figure 4. Effect of species proportion on European beech basal area growth per unit species proportion (a) depending on the quadratic mean diameter and (b) depending on the growth index.All the other explanatory variables are set to their mean values.GI Eb : growth index for European beech.

Figure 4 .
Figure 4. Effect of species proportion on European beech basal area growth per unit species proportion (a) depending on the quadratic mean diameter and (b) depending on the growth index.All the other explanatory variables are set to their mean values.GIEb: growth index for European beech.

Figure 5 .
Figure 5. Basal area increment of two NFI sites using 3P proportions: (a) the site showing the lowest underyielding for 50%-50% proportions; and (b) a site showing a transgressive overyielding.The basal area increment was predicted for Dg = 35 cm, sd = 0.65 and the ecological conditions indicated on the respective graphs for European beech (blue), Norway spruce (red) and the total stand (black).Underyielding (or overyielding) is the negative (or positive) difference between the predicted basal area increment in mixture (continuous line) and the expected increment under the assumption of no mixing effects (dashed line).GINs: growth index for Norway spruce.

Figure 5 .
Figure 5. Basal area increment of two NFI sites using 3P proportions: (a) the site showing the lowest underyielding for 50%-50% proportions; and (b) a site showing a transgressive overyielding.The basal area increment was predicted for Dg = 35 cm, sd = 0.65 and the ecological conditions indicated on the respective graphs for European beech (blue), Norway spruce (red) and the total stand (black).Underyielding (or overyielding) is the negative (or positive) difference between the predicted basal area increment in mixture (continuous line) and the expected increment under the assumption of no mixing effects (dashed line).GI Ns : growth index for Norway spruce.
: Predicted versus observed data plots and standardized residuals versus explanatory variables plots for the Norway spruce growth model, Figure S3: Predicted versus observed data plots and standardized residuals versus explanatory variables plots for the European beech growth model, Figure S4: Replacement diagrams of the NFI site showing the lowest underyielding (using the 3P alternative), Figure S5: Replacement diagrams of the NFI site showing the mean overyielding (using the 3P alternative), Figure S6: Replacement diagrams of the NFI site showing a transgressive overyielding, Table

Table 1 .
Main stand variables of the sample plots used for the growth modelling.

Table 2 .
Summary of the four alternatives employed to estimate BA max : corresponding formulas and origins of the coefficients introduced.

Table 3 .
[43]ficients of the self-thinning rule estimated using quantile regression (τ = 0.95) with the NFI data and coefficients found by Charru et al.[43]using stochastic frontier analysis.
n: sample size.The standard errors for the coefficients obtained using the SNFI data are given in parentheses.

Table 4 .
Fixed effect parameter estimates for the Norway spruce and European beech growth models using 3P proportions.

Table 4 .
Fixed effect parameter estimates for the Norway spruce and European beech growth models using 3P proportions.
SE: standard error.