From Acid Rain to Low Precipitation: The Role Reversal of Norway Spruce, Silver Fir, and European Beech in a Selection Mountain Forest and Its Implications for Forest Management

: Research Highlights: We make use of long term observation data from a selection forest in Bavaria. Despite the changing environmental conditions, stand level productivity remains constant over time. Maintaining species and structural diversity by forest management can contribute to resilient forest ecosystems. Background and Objectives: Forests in mountains are similarly affected by environmental changes like those in northern latitudes as species are closer to the edge of their ecological niche. There are recent studies that report species-speciﬁc responses to climate change in unmanaged, mono-layered mountain mixed forests. We analyze how environmental changes modify the growth of multi-layered, managed selection forest, which are often targeted for stabilization and risk prevention. We pose the central hypothesis that different species-speciﬁc susceptibility to disturbances and structural diversity contribute to ecosystem stability. Materials and Methods: Based on the long-term experiment Freyung 129 in the montane zone of the Bavarian Forest, Germany we analyze long term chronologies of periodic single tree and stand growth of Norway Spruce, silver ﬁr, and European beech in dependence of environmental factors and forest management. Results: First, we show that despite environmental changes in terms of air pollution and drought stress, productivity at stand level persists constantly because of structural diversity and species traits. Second, we show that the species-speciﬁc contribution to total stand growth and growth distribution among stem diameter classes may change over time; the species interactions balance total growth. Third, we reveal a role reversal of tree species growth pattern. N. spruce was superior in growth in the ﬁrst half and was replaced by s. ﬁr in the second half of the survey period. Fourth, we identify the interplay of different stress factors on species-speciﬁc growth as the main cause for species-speciﬁc asynchronous but growth stabilizing reaction pattern. Finally, we show that density regulation was limited in its impacts to mitigate prevailing stress factors. Conclusions: We discuss the reasons for the observed stability of productivity. We interpret results, where especially the diversity of species and structure typical for selection forests result in stable productivity and wider plateau of the density-productivity relationship, and the suitability of the selection forest concept for risk prevention and stress resilience. We conclude that species composition and stand structure of selection forestry in mixed mountain contribute to climate smart forestry.


Introduction
Structurally heterogeneous mixed-species stands are widely favored as they fulfil many ecosystem services better than mono-cultures and are supposed to be more stable under changing environmental conditions. Therefore, homogenous and mono-specific forest stands currently are frequently transformed into more heterogeneous stands [1][2][3]. Stand structure in (a) even-aged normal forest and (b) selection forest. The age-class forest (a) consists of a spatial mosaic of stands with varying age, tree composition, and structure. The selection forest (b) has a wider range of species and structural diversity at the stand level, but retains a similar structure across space and time (after [17]).
The species diversity in selection forest composed of N. spruce, s. fir, and E. beech may have a stabilizing effect as N. spruce is fast growing under normal conditions, less sensitive to acid rain than s. fir but highly susceptible to windthrow, drought, and bark beetle attacks [18]. S. fir shows higher mechanical stability, higher drought resistant due to deeper rooting, but is particularly sensitive to SO 2 -imissions [19]. E. beech is comparatively slow growing, less prone to windthrow but reactive to increased ozone concentrations [20]. In contrast to N. spruce, E. beech behaves anisohydric under drought; i.e., continuous assimilation and growth even at water deficit, whereas isohydric N. spruce and s. fir close stomata and reduce growth in the early phase of drought [21]. The mixture may help to protect beech against frost [22,23], may shade trees from strong radiation and sunburn [24,25], may improve humus and nutrient cycling [26], and augment water supply [27,28].
The structural diversity in selection forest composed of N. spruce, s. fir, and E. beech is based on the medium to high shade tolerance [29,30] of all three species and morphological plasticity especially of fir and beech [31]. All three species can remain for decades or even centuries in the understory [32], accelerate growth after canopy openings even after long periods of suppression [33], and endure shading and lateral crown limitation in the middle layer [8]. This enables a multi-layered canopy structure with all three species present from predominant to understory positions [34]. The inter-specific diversity of stress susceptibility together with their ability to form multiple-layered stands creates the potential of high resilience and stability against abiotic as well as biotic disturbances [35,36].
In changing environments with shifts in disturbance regimes e.g., triggered by pollution or climate change these traits are getting more prominent. In a fluctuating environment with frequent stress events, they may fulfil the insurance hypothesis [37] in a twofold way. On the one hand, the species mixture means a risk distribution; species show different reactions to stress and the growth decrease or drop out of one species may be mitigated or compensated by other less affected neighboring tree species. On the other hand, the multiple layers promote stability in growth; damages or tree losses in one layer may be compensated by the remaining trees in another layer. For instance, canopy opening caused by wind throws or bark beetle may be immediately compensated by the previously subdominant or suppressed trees of lower layers. So, the growth resilience is higher than in age class forests where similar disturbance may cause open spaces, productivity loss, until a new stand generation is established.
Nevertheless, the overwhelming majority of research into the stress resistance and resilience of tree species [38,39] and the silvicultural prescriptions for adaptation [40] refer to mono-specific stands or even-aged, mono-layered mixtures [41]. For mono-layered mixed mountain forests of N. spruce, s. fir, and E. beech throughout Europe [42] recently showed that the growth remained relatively constant over the last 30 years on a level of 9.8 m 3 ha −1 year −1 . This was about 20% above the productivity expected by common yield tables and amazingly constant despite various stress impacts such as acid rain [43], drought [44], and late frost events [45] during this period. However, the stable productivity at stand level was accompanied by a species-specific shift in growth with N. spruce declining and s. fir increasing in growth. The database included experiments established and re-measured between 1906 and 2017. Another study [46], based on tree cores, revealed a general increase of growth for N. spruce, E. beech, and s. fir in mountain regions whereas the increase was more pronounced for E. beech and s. fir than for N. spruce. Since these findings refer mainly to mono-layered stands, it is still an open question how uneven-aged mixtures, and especially selection forests that are widely pursued by climate adaptive forest management [1,2,17,47], perform under environmental changes.
For further insight into the response of uneven-aged and mixed species forests to environmental stress, we here put to test the growth stability of selection forests under environmental disturbances. We pose the central hypothesis that the different species-specific responses to natural and human-induced disturbances contribute to a stability of growth and flexibility of harvest options. Based on the long-term selection forest experiment Freyung 129 with species N. spruce, s. fir, and E. beech, we tested the following hypotheses.

Hypothesis 1 (H I).
Stand growth in surveyed stands remained stable over the last 40 years.

Hypothesis 2 (H II).
Species specific growth distribution among diameter did not change over time.

Hypothesis 3 (H III).
The long-term growth trend at the individual tree level shows no change over time.

Hypothesis 4 (H IV).
If observed, any growth trends can be assigned to environmental factors.

Hypothesis 5 (H V).
Diversity and stand density regulation can reduce environmental stress effects.
We discuss the potential of selection forest management to cope with various stress regimes. We further argue that species and structural diversity contribute to stabilize managed forest ecosystems.

Material
The selection forest experiment Freyung was established in 1980 in the Bavarian Forest 5 km northeast of the city of Freyung (longitude 13.58 • E, latitude 48.85 • N). It is located in the Kreuzberg Forest at the southwest border of the ecoregion "11.3 Innerer Bayerischer Wald". This region belongs to the montane to high-montane zone of the German-Czech border mountains and geologically to the Bohemian mass [7] (pp. [26][27]. The total 1.5 ha sized experiment is located in a 160 ha area which has been managed as a selection forest for the last centuries [9,48,49]. The forest owner presently transforms additional parts of its in total 560 ha forest to selection forest due to the advantageous productivity and stability of this type of silviculture. The selection forest experiment Freyung 129 is positioned at a slight southeast exposed slope at 720-730 m elevation a.s.l. The stands are stocking on podsol brown soil of sandy loam material originating from weathered gneiss and granite parent material. The annual precipitation amounts to 1200 mm year −1 (700 mm during the vegetation period from May-October) and has a mean temperature of 6.5 • C. The area lies in the transition zone between Atlantic and continental climates. Between 700-1200 m a.s.l. N. spruce, s. fir, and E. beech are naturally associated with each other (vegetation unit Luzulo-Fagion according to [50]). Vitality and growth of E. beech is occasionally impaired by frost, growth of N. spruce by bark beetle. In the 1950-1980s, vitality of s. fir was strongly affected by SO 2 immissions [51,52].
The site has experienced environmental changes in the last 40 years ( Figure 2). Mean temperature increased by approx. 1 • C, whereas mean precipitation slightly decreased from 1980 to 2020 [53,54]. Both trends result in a negative progression of soil water availability as indicated by the soil moisture index [55]. Meanwhile, national emissions of SO 2 were drastically reduced by 96% (2018) of the value in 1980 [56,57]. Emissions caused a large-scale mean annual atmospheric concentration of 25 to 50 µg m −3 (regional: 50 to 75 µg m −3 ) in the 1980s. Now, the concentrations fall largely below the threshold for protecting ecosystems of 20 µg m −3 .
The experiment is divided into three main plots; of which each is split into two subplots. Subplots are numbered first according to main plot (1-3) and then according to subplot (1)(2). Each subplot has a size of 0.25 ha. Due to the trial's aim to assess the interplay between stand density and stand structure, all subplots were kept under different but more or less constant stand densities levels. Desired densities were established by selective cuttings on each plot. In this way, plots at lower densities (FRY129/11, FRY129/12, FRY129/21) faced more frequent and more intensive selective thinning activities than the others (FRY129/22, FRY129/31, FRY12932). Table 1 summarizes the main measurement variables and metrics used in this study. From each tree, we recorded the species identity, measured the x and y stem coordinates at the first survey and all stem diameters during each of the 7 surveys. Tree height, h, height to the crown base, hcb, were sampled (>30 trees per species) at each survey.  To provide an overview of stand development, we first evaluated the consecutive inventories using standard procedures [58][59][60]. Single tree volume was calculated based on basal area, tree height, and species-specific form factors. Tree height was retrieved from height-diameter relationships of sampled tree heights distributed over the diameter distribution. Both information was calculated for remaining trees and removal trees per observation period. Single tree values were summed up to stand level. The results encompass mean tree dimensions as well as stand volume, and volume growth.

Species Level Evaluation
At the species level, we calculated the shares of the three tree species N. spruce, s. fir, and E. beech of the total periodic annual volume increment (PAIV). We calculated the respective shares by dividing the PAIV of a species by the total PAIV of all tree species taken together (PAIV species /PAIV total ). We made these calculations for all six subplots and each observation period.

Diameter Distribution Level
For the exploration and modeling of the size-dependent distribution of the periodic annual stem volume growth, we used 10-cm diameter classes. We calculated the distribution for each subplot and survey period (see exemplarily Figures S4-S6). For further statistical analysis of the stem volume growth-diameter class distributions, we used the respective mean diameter of each DBH-class (e.g., 5 cm DBH in DBH-class 0-10 cm). To model the distribution and its temporal development, we fitted the distributions with generalized additive models (GAM) [61] separately for the three tree species N. spruce, s. fir, and E. beech and also the distribution for all species together. By this analysis, we revealed any change in the growth distribution, deviations from the steady state and compensation of one tree species by the other in term of their contribution to the stand structure and growth. Here, we focused on the subplots FRY 129/22, 129/31, and 129/32 which were kept at higher stand densities from 1980 to 2018 and faced only moderate density fluctuations in course of the periodic selection cuts to avoid cofounding effects of thinning on single tree growth. Since beech only played a minor role in the plots, it was of particular interest to investigate whether there occurred shifts in the contribution of the two species N. spruce and s. fir to the volume growth in the respective diameter classes. We investigated possible shifts using the relative periodic volume increment (RPAIV) per diameter class. This is the modeled periodic annual volume increment (PAIV) of the stronger species divided by the PAIV of the weaker species in each diameter class minus one; RPAIV = PAIV stronger species /PAIV weaker species − 1). In the case that N. spruce was the stronger of the two species in a diameter class, we also multiplied the corresponding relative PAIV by minus one.

Tree Level Evaluation
Tree height estimation: For estimating the species-specific height, h, of each tree depending on the stem diameter, d, we parameterized the allometric relationship: The assessment was based on n = 2919, 2360, and 1216 combined hand dmeasurements of N. spruce, s. fir, and E. beech between 1980 and 2018. The measurements were distributed over the whole diameter range. The model fitting resulted in the parameters a 0 = 0.31 and a 1 = 0.80 for N. spruce, a 0 = 0.16 and a 1 = 0.86 for s. fir, and a 0 = 0.96 and a 1 = 0.62 for E. beech with R 2 values of 0.95, 0.93, and 0.85, respectively (Supplementary Table S1). The individual tree heights predicted by this model were used as basis for the delimitation of the following neighborhood analysis.
Edge correction: Before calculating single tree-specific local competition values and mixing proportions for neighborhood analysis, we established a toroidal shift of the plot to all eight directions of the plot periphery for edge bias compensation [59,62,63]. Using the toroidal shift, we extended the same plant position pattern and distances in all eight directions and avoided any overestimation of density, as it could result from other techniques [63].
Neighborhood analyses within sample circles: To characterize the individual trees' competitive status, we quantified the local stand basal area (BAL). For this purpose, we defined an influence zone by a circle with search radius sr 1 = 0.25 × h 1 around the position of each tree. All trees within the circle except the center tree were used to calculate the local basal area (ba) on the circle area a. BAL = ba/a was the respective basal area upscaled to one hectare. In the constructed circles, there were, on average, 20-30 trees and at least 10-15 most impactful neighbors [64]. The BAL values were calculated with and without the removed trees (BAL pre , BAL post ); based on both values, we quantified the competition release by the tree removal: Mixing proportions: The trees sampled in the circle were also used to calculate local mixing proportions. The mixing proportions mportion 1 . . . mportion n should reflect the proportions of two or more species in the observed mixed stands [65,66]. Tree number, basal area, or volume proportions are only appropriate for this purpose if the mixed species have similar growing area requirements [67]. The considered tree species vary per se in the growing area requirement and maximum stand density in fully stocked stands. For example, a E. beech with a stem diameter of 25 cm may require approximately double the growing space as a N. spruce of the same diameter. This means the density in terms of trees per hectare is only half of that of N. spruce. To standardize the density and to calculate the unbiased area-related mixing proportions, we applied the equivalence factors according to [68]. These factors for the main tree species assemblages are shown in Supplementary Table S2.
Soil moisture Index (SMI): The Soil Moisture Index (SMI) is an indicator to quantify agricultural droughts. The SMI is simulated by the hydrologic model mHM using current weather data. The model mHM consists of a digital elevation model, geological map, soil map, and land cover and leaf area information derived from satellite data (see [69] for further explanation of the model). The model provides 30-day soil moisture indices SMI [70,71]. The range of the values is SMI = 0-1. A value SMI < 0.2 is considered as drought. The index is primarily used to provide information on the current moisture status of the soil. Applied mainly in climate research, the index is used for national drought monitoring in Germany [69]. Furthermore, the SMI has been successfully related to tree-ring series in dendroecological studies to identify drought events [72].

Statistical Models
For testing the overall growth trend of the stands (H I) we fitted the models: in order to reveal any changes of IV by the volume of the removal trees, V remove , or the calendar year. Whether the long-term growth trend at the individual tree level has changed over time (H III) was tested by fitting the model: For the revelation of the causes of the long-term effects (H IV) we fitted the model: Another model was fitted to relate growth and driving variables of growth for the periods 1980-1993 (Model 6a) as well as 1993-2018 (Model 6b) in order to reveal any changes in the species-specific reaction patterns in the first compared to the second half of the survey period.
Model 6 reveals the effects for the first and second half of the survey period how the stress effect modified the growth of the three species and how selection cutting were able to modify the environmental impact (H V).
For testing H I, H III, H IV, and H V, we applied linear mixed models. The dependent variables were either the stand stem volume growth, IV (Models 2 and 3), or the mean annual stem diameter growth, id (Models 4-6). The independent variables were in the case of Models 2 and 3 the standing stem volume, V, the volume of the removal trees, V remove , the soil moisture index, SMI, and the calendar year, year. In the case of Models 4-6, the independent variables were the individual tree diameter, d, local basal area, BAL, release by tree removal, ∆BAL, mixing proportion, mportion, soil moisture index, SMI, and the calendar year, year. In all equations, the indexes i and k represent the k th observation of the i th plot (Models 2 and 3) or tree (Models 4-6). The fixed effects were covered by the parameters a 0 -a n . With the random effect b i ∼ N 0, τ 2 , we cover the correlation between the single observations on plot (Models 2 and 3) and tree (Models 4-6) level. In preliminary model formulations of the Models 4-6, we also worked with random effects on plot level, i.e., one additional nesting level. As this caused confounding effects with the fixed effect, we constrained ourselves to the simpler random effect structure of Equations (3)- (6). With ε ik , we denote the independently and identically distributed errors.
Hypothesis H II, whether the overall growth distribution on the trees in different diameter classes was modified over the entire survey period 1980-2018 was tested by fitting a generalized additive model (GAM) to the data of the subplots FRY 129/22, 31, and 32.
PAIV DBHclass = period + s(DBH class mean ) + s(DBH class mean , by = period) Here, we analyzed any growth change in the respective DBH-classes of the distribution over time (period = single observation period); no significant differences between the observation periods would indicate a steady state. Any deviation between observation periods in each DBH-class would indicate temporal changes.
For assessing whether there were changes in the PAIV within the DBH-classes between the species, we first fitted a generalized additive model for each observation period individually.
PAIV DBHclass = Species + s(DBH class mean ) + s(DBH class mean , by = Species) Using this model, we were able to derive the corresponding relative annual periodic volume increments (RPAIV) in the 10 cm diameter classes (see Section 2.3). To test whether there were temporal changes in the RPAIV values in the respective diameter classes, we used another generalized additive model for the plots Freyung 129/22, 31, and 32.
RPAIV DBHclass = s(DBH class mean ) + s(period) + s(DBH class mean , period) We restricted evaluation for testing HII and HIII on data of the subplots FRY129/22, 31, and 32 where the standing volume was kept continuously at a high level since the beginning of the survey. In this way, we could eliminate silvicultural effects on growth dynamic and elucidate environmentally triggered changes.

Tree and Stand Characteristics
The dendrometric stand characteristics (Table 2) show that the tree sizes are quite similar on the different plots (e.g., mean stem diameter 21.46-28.87 cm) but the stand density differs strongly. As a result of the silvicultural density regulation, the tree number, basal area, or standing volume of plots with high density show doubled values compared to the most sparsely stocked plot (e.g., 373.57 vs. 873.11 m 3 ha −1 ). However, the volume growth differs less strongly ( Due to the lack of printing space, we restrict the introduction of the experimental stands to the information essential for understanding this study. For details of the stand structure, overstorey growth, and regeneration, see [9,17,48,49]. Typically, for selection forests, the trees on the plots belong to very different age phases and social positions. Table 3 reflects the wide range of tree state variables and stem diameter increment represented on the six plots. The tallest trees achieve 89.5 cm stem diameter, 42.8 m tree height, 31.40 m crown length, and 16.18 m crown diameter. E. beech has on average lower tree heights but wider lateral crown expansion, longer crowns, and higher crown ratios. Altogether, more than 6000 stem records from the repeated surveys and their variation in size and tree development and competitive state provide a solid database for the subsequent evaluation of the tree growth depending on tree status, competition, and environment.
The mean local densities before and after thinning, BAL pre and BAL post , respectively, are higher for the highly shade tolerant s. fir and E. beech compared with N. spruce, whereas the minimum and maximum values are similar for all species. The release by selection cutting, ∆BAL, was the highest for s. fir. The mean proportion of other species in the neighborhood, mportion was 0.50, 0.56, 0.75 for N. spruce, s. fir, and E. beech; this reflects that the stands were dominated by the two conifers.
The mean stem diameter growth plotted over the calendar year revealed a clear species-specific trend (Figure 3a-c). Obviously, N. spruce decreased in growth since the start of the survey in 1980; in the middle of the survey period, the decrease was nearly linearly, during the last 20 years, the decrease lessened. S. fir, in contrast, continuously increased in growth since 1980; the mean growth rate more than doubled. E. beech showed a unimodal course of growth with a peak in the 1990s when growth of N. spruce was already reduced and s. fir not yet reached the maximum. Additionally, to the course of the mean stem diameter growth depending on calendar year (Figure 3), the Supplement Figures S1-S3 show the growth of all individual trees versus calendar year and the growth plotted over stem diameter. The shown speciesspecific growth trends may be co-determined by the stand structure and silvicultural treatment of the plots. To avoid respective biases, we considered those effects in the following analyses.  Figure 4a shows that the stem volume growth on the plots FRY 129/11-32 ran on different levels and varied, but remained on a rather constant level of about 10-15 m 3 ha −1 yr −1 since 1980. There was a tendency that the plots with the higher standing volume at the beginning of the respective survey periods had higher productivity than the plots with lower stand volume (Figure 4b).  Table 4); we inserted mean values for the stand stem volume, V, in case of Model 3. Note, that in (d) lines for different SMI are on the same level as it had no significant effect (see Table 4). Model 2: ln(IV) = a 0 + a 1 × ln(V ik ) + a 2 × ln(V remove ik ) + b i + ε ik ; Model 3: ln(IV) = a 0 + a 1 × ln(V ik ) + a 2 × ln(year i ) + a 3 × ln(SMI i ) + b i + ε ik ; Model 4: ln(id) = a 0 + a 1 × ln(d ik ) + a 2 × ln(year i ) + a 3 × ln(d ik ) × ln(year i ) + b i + ε ik ; Model 5: ln(id) = a 0 + a 1 × ln(d ik ) + a 2 × ln(BAL ik ) + a 3 × ln(∆BAL ik ) + a 4 × ln(year i ) + a 5 × ln(SMI i ) + a 6 × ln(year i ) × ln(SMI i ) + b i + ε ik ; Model 6: ln(id) = a 0 + a 1 × ln(d) + a 2 × ln(BAL) + a 3 × ln(∆BAL) + a 4 × ln(mportion) + a 5 × ln(SMI) + b i + ε ik .

Development of Stand Level Growth (H I)
We further analyzed the dependency of stand stem volume growth on volume and removed volume by selection cuttings (Figure 4c, Table 4 Model 2). Stem volume growth slightly decreased with lower stand densities; the model coefficient a 1 indicates that a density reduction of the volume by 50% reduced the productivity only by~30%. This indicates a subproportional effect of density reduction on stand growth; in case of a proportional effect, a 1 would equal 1. The removal stand volume at the beginning of the survey period increased the volume growth. Especially moderate selection cutting by removal of 50-100 m 3 ha −1 increased the productivity, whereas stronger density reductions had only minor additional effects (saturation curve).
Neither the calendar year nor the soil moisture index had a significant effect on the stand volume production ( Figure 4d); indicated by slopes of zero in the regression of stem volume versus calendar year and SMI ( Our results about the shares of the three tree species N. spruce, s. fir, and E. beech in the total volume growth showed that they have changed substantially since the beginning of the experiment in 1980 ( Figure 5). The share of s. fir concerning volume growth has steadily increased in all subplots. The proportion of N. spruce has decreased to almost the same extent as the proportion of s. fir has increased. E. beech, on the other hand, has remained almost constant in its share of volume growth. Although the shares in volume growth on the subplots FRY 129/11, 12, and 21 were more strongly co-determined by selection cuttings, it is of special interest that we found this pattern also on the moderately thinned and highly stocked subplots FRY 129/22, 31, and 32.

Size Class Contribution to Stand Growth over Time (H II)
By fitting Model 7 to the total volume growth distributions (m 3 ha −1 year −1 for DBHclasses) of all tree species together on the subplots FRY 129/22, 31, and 32, we found a significant deviation of the distribution only for the first observation period 1980-1986 (Table 5, Figure 6) in the 50-60 cm, 60-70 cm, and 80-90 cm diameter classes (the confidence interval excludes zero). Since 1986, the total volume growth distribution has not shifted significantly on these plots ( Figure 6). From 1986 to 2018, we could not detect any significant changes in any of the 11 diameter classes.  If all three tree species are considered together, the volume increment distributions on the fully stocked subplots FRY 129/22, 31, and 32 are almost constant since 1986 ( Figure 6). Figure 5 showed, however, that the shares of the two tree species N. spruce and s. fir in the total PAIV have changed substantially even after 1986.
It was of particular interest to see in which diameter classes this shift has taken place. Our results from Model 9 ( Table 6) showed on the one hand that there was a significant temporal change in the relative PAIV between N. spruce and s. fir, but they also provide information on the diameter classes in which these shifts took place (Figure 7). S. fir had at any time higher RPAIV in smaller diameter classes (10-50 cm DBH). In the diameter classes 50-60, 60-70, and 70-80 cm DBH, we observed a change in growth dominance between the two species N. spruce and s. fir. In the first surveys, N. spruce still had a higher share of growth in these diameter classes. However, during the last surveys, it was s. fir that showed a higher share of growth in the respective diameter classes. N. spruce had at any time higher RPAIV in higher diameter classes (80-110 cm DBH). However, it is remarkable that the growth dominance of spruce in these diameter classes has declined substantially since 1980. In the 110-120 cm diameter class, the situation was almost balanced between the two species.  7. A comparative plot of the relative annual periodic volume increment (RPAIV) between the tree species Norway spruce (red) and silver fir (grey) modeled according to Equation (9). Bars within a DBH-class represent, from left to right, the observation periods 1980-1986, 1986-1993, 1993-1999, 1999-2005, 2005-

Change of Growth and Growth Partitioning from 1980ies to Present (H III)
The visual comparison between the growth-size relationship of the three species (Figure 8a,b) shows a fundament change of the species ranking from the first survey in 1980 to the last in 2018. Under constant environmental conditions and steady state, structure and volume on the plots the species-specific level and slope of the id-d relationship should be similar in 2018 compared to 1980. However, s. fir was inferior to N. spruce in the past and shows superior growth at present. E. beech also improved in growth, but to a lower extent than s. fir. For N. spruce, applies the opposite; it was superior in the past and decreased in the id-d relationship.
In addition to the visual comparison, Model 4 revealed that the growth changed significantly depending on the calendar year within the survey period from 1980-2018. For N. spruce and E. beech, the stem diameter growth decreased since 1980 (see Model 4 in Table 4, a 2 = −178.24 and a 2 = −158.50, respectively). For s. fir, the regression coefficient in Table 4 reveals a significant increase of the level (a 2 = 140.24) of the stem diameter growth relationship. Both N. spruce and E. beech show a decrease of the level and increase of the slope of the id-d relationship on the high density plots from 1980-2018. As we used only data from subplots less affected by selective cuttings, we identify environmental changes responsible for the role reversal in favor of s. fir within the last 40 years.  Table 6).

Impact of Environmental Factors on Growth Development (H IV)
S. fir continuously improved in stem diameter growth since 1980 according to the evaluation by Model 5 (Figure 9a). Note that we modified the calendar year, but set the other model variables constant to show this significant trend. N. spruce decreased significantly, and E. beech showed no trend between 1980-2018. Thus, the ranking between the species in terms of growth changed from initial N. spruce > E. beech > s. fir to the opposite at the end of the survey period. At parity of diameter, local stand density, and SMI stem growth of silver was more than threefold in 2018 compared with 1980. All three species react positively to SMI, however, s. fir much stronger than N. spruce and E. beech (Figure 9b). In dry years, s. fir reduced growth stronger than N. spruce and E. beech, but benefitted most from improving water availability.
The growth of all three species depended similarly on the local stand density (Figure 9c). Stem diameter growth decreased continuously with increasing local stand density BAL. This reaction pattern was shown for the year 2000, SMI conditions of 0.5, and trees with mean diameter of d = 25 cm, however, the relationship between the species changed with the SMI conditions.

Tree Growth Depending on Environemnatl Conditions in the Past and at Present (H V)
In the period 1980-1993 when the atmospheric SO 2 -load was high, we found a high sensitivity of silver to drought, indicated by a strong reduction of stem diameter growth with decreasing SMI (Figure 10a). For N. spruce and E. beech, we found no dependency of growth on SMI (broken lines). This caused a strong decrease of the relationship between id and d of s. fir even by a moderate reduction from SMI = 0.50 to SMI = 0.45 (Figure 10b,c). In the period 1980-1993, different tree species in the neighborhood were beneficial for N. spruce and E. beech, but not for s. fir (Figure 10d). Density reduction slightly increased the stem diameter growth of N. spruce and s. fir but not the growth of E. beech (Figure 10e). Selection interventions had a slightly positive effect on the growth of s. fir and E. beech, whereas N. spruce reacted independently from the intervention intensity (Figure 10f). Figure 10. Stem diameter growth, id, of Norway spruce, silver fir, and European beech on the selection forest experiment Freyung 129 in the period with high atmospheric SO 2 -load (1980-1993). Visualized is the modification of id by (a) Soil Moisture Index, (b) stem diameter growth in periods with SMI = 0.50, (c) stem diameter growth in periods with SMI = 0.45 (d) mixing portion in the neighborhood, (e) local stand density, and (f) density reduction by selection cutting (according to Model 6a, see Table 4). Broken lines indicate non-significant, invariant relationships.
In the period with low atmospheric SO 2 -load (1994-2018), s. fir grew rather independent from SMI, whereas N. spruce and E. beech reduced their growth significantly under drought. The growth stability of s. fir under drought may partly result from the reduced water consumption by the neighboring species (Figure 11a). In wet years, N. spruce is superior in growth (Figure 11b), however, under drought, s. fir turns the best and N. spruce and beech strongly reduce their stem growth. The slight increase of fir in dry years may result from a competition reduction due to the more drought affected neighboring N. spruces and beeches (Figure 11c). This is in line with the finding that silver benefits from species admixture in the neighborhood, whereas N. spruce and E. beech hardly react on increasing mixing proportion in the vicinity (Figure 11d). S. fir show the highest growth rates over the whole range of local stand densities (Figure 11e) and also reacts most strongly on a given density release by selection cutting (Figure 11f). Figure 11. Stem diameter growth, id, of Norway spruce, silver fir, and European beech on the selection forest experiment Freyung 129 in the period with low atmospheric SO 2 -load (1993-2018). Visualized is the modification of id by (a) Soil Moisture Index, (b) stem diameter growth in wet periods, (c) stem diameter in dry periods, (d) mixing portion in the neighborhood, (e) local stand density, and (f) density reduction by selection cutting (according to Model 6b, see Table 4). Broken lines indicate non-significant, invariant relationships.

Individual Tree Resilience Versus Stand Level Resilience
In contrast to the species-specific reaction patterns of individual tree growth, stand productivity remained rather constant since 1980. In essence, stand productivity depended on density and structure regulation but hardly on environmental conditions. Obviously, the species-specific stress and growth reduction of one species was buffered by the increase of another resulting in a neutral net effect on stand productivity. Note, that an analysis at the individual tree level indicated a growth decrease for s. fir in 1980-1990 and for N. spruce in the 2010-2018. A deduction of stand growth based on the finding at the individual tree level would be misleading without considering the counteracting and counterbalancing effect between the tree species. Upscaling of growth reactions from tree level to stand level should be made cautiously; a diversity of species traits, stress susceptibility, growth compensation between different layers and asynchrony of stress reactions can buffer and neutralize growth patterns at the stand level that are alarming at the tree level. On the other hand, our analysis showed that findings of a constant productivity at the stand level may be misguiding; thus, we were able to show with our results that, if all tree species are considered together, the distribution of volume growth across the diameter classes has not changed since 1986. However, we were also able to show that the shares of the three tree species in volume growth in the respective diameter classes have shifted significantly over the last 40 years. Especially in the smaller and medium diameter classes, N. spruce showed a declining share. This declining share of N. spruce was compensated by the restrengthening of s. fir, so that the volume growth distribution at the stand level remained constant. However, it also led to the fact that fir has become dominant in terms of standing stock and volume growth. These findings point out that despite constant accumulative growth, the contribution of different species and layers may change and thus the system structure may also change.

Diversity Promotes Stability
The recovery of s. fir after reduction of SO 2 -emissions is a model example for the society's capacity to advocate, decide, and act in favor of nature including forest ecosystems. Starting with the environmental political measures of industry emission reduction in the 1970s and 1980s focusing on SO 2 -emissions of the coal electric power plants air pollution has been significantly reduced and simultaneously the far-distant transport of pollution to the low and high mountains among other to the Bavarian Forest [51,52]. This enabled especially the recovery of s. fir which is highly susceptible to damage caused by sulfur emissions [19].
In selection forests, particularly dominant s. firs had suffered from the emissions; subdominant und understory s. firs survived, could recover after emission reduction, and recaptured their key role for the sustainable productivity and functioning beside N. spruce and E. beech. The ability of s. fir to "sit and wait" and survive in subdominant position enabled its comeback after the strong decline due to acid rain.
In the 1980s, s. fir was on a low growth level in all diameter classes. In addition, it reacted very negatively on dry years. This means that the vitality reduction by SO 2 was linked with an increased drought sensitivity [51]. Evaluations of the crown vitality showed a high crown transparency due to needle losses in the 1980s and 1990s; this needle shedding probably caused a growth reduction, deeper irradiation into lower canopy layers, and maybe promoted N. spruce and E. beech which in return may have increased the competition for water disadvantageous for s. fir.
Maintaining tree species, even when vitality is reduced, may be advantageous in future times. Keeping and promoting s. fir in phases of strong emissions despite damages creating a potential for strong recovery soon afterwards. A further pro for the s. fir was its lower susceptibility to drought stress when other stress factors are absent. That became relevant in the second half of the survey period. Since the 1990s, s. fir not only stabilized stand growth by its recovery from emission damages but also by compensating the growth losses of N. spruce and E. beech caused by drought events, among others in the years 2003 and 2015.

Management of Diversity
The stable productivity within the study plots were so far based on their specific tree species diversity. In the 1980s, N. spruce had a share of 56% of the standing stock and 63% of the stand productivity, so that the growth losses of s. fir caused by SO 2 could be compensated by N. spruce. The other way round, the increasing drought stress of N. spruce since 2000 was compensated by the steadily recovering and more drought resistant s. fir. Due to its reviving vitality, s. fir became dominant in terms of standing stock and growth of volume, whereas N. spruce decreased. E. beech was kept on a rather low level, due to its inferior wood quality. The overrepresentation of s. fir and its essential contribution to the stand productivity is advantageous now, however, suppose it would be impaired by new biotic or abiotic disturbances the stands might strongly suffer due to this dominance. Anticipatory and preventive silvicultural management might balance the proportions of the three species or even increase the species diversity by adding stress resistant provenances of present species or even additional species such as Douglas-fir (Pseudotsuga menziesii) or yew tree (Taxus baccata), to ensure genetic or species diversity and to stabilize productivity according to the insurance theory [37]. The presently high productivity may be a misleading indicator for general productivity and system stability without considering how it is backed by tree species diversity. Maybe stand productivity divided by the tree species evenness would be a better indicator.
When stand density is reduced in selection forests, productivity eventually also decreases; however, the decrease starts off at much lower density and proceeds much slower than in even-aged monocultures. The reason for this broader saddle and lower steepness of the stand density-productivity relationship lies in the trees of the medium and lower layer, which can buffer and compensate the growth losses of the removal stand immediately [77].
The primarily tree-wise interventions and removals only cause small gaps, which can be closed quickly by trees of lower layers and by neighbors. The halving of the standing stock from 600 to 300 m 3 ha −1 on one subplot of the Freyung experimental resulted in a growth reduction of only 20-30% compared to the fully stocked stand (Figure 4c). It indicates a high growth resilience of the selection forest to silvicultural or natural disturbances compared with the age class forest. In the latter, any gaps in the canopy may be closed much slower due to the mono-layering and lack of natural regeneration.
The potential to mitigate impacts of stress situations on stand level growth by density management needs to take species-specific traits concerning stress susceptibility into account; our results show that the species-specific effect of density regulation very much depended on the prevailing stress factor.

Methodological Consideration
Compared to other competition indices, the local stand basal area BAL is easy to interpret and to link with practical decision making; stand basal area is a well known measure which is easy to assess e.g., by relascope sampling [78] and frequently used in silvicultural models and prescription [59,79]. The search radius of half of the tree height caused a size dependent delineation of neighborhood analyses as proposed by [80] and resulted in BAL and ∆BAL values that correlated closely (R 2 > 0.6) with the stem diameter growth. By application of density equivalence factors [65], we considered the speciesspecific growing space requirements when merging/summarizing the effect of different competitive strength of neighboring species in one index. Alternative indices, e.g., based on crown size or leaf area, would have required more detail survey data; BAL and ∆BAL used only stem diameter and height which were measured most frequently since the start of the survey.

Relevance for Climate Smart Forest Management
The concept of Climate Smart Forestry [81] aims at enhancing the adaptive and mitigative potential of forest in the frame of climate change, concurrently providing multiple benefits for society. The stable growth of the here investigated selection forest despite changing environments and associated disturbance patterns indicates a high level of resilience of the system. This level is achieved by mixing species of different traits and susceptibility to stress factors as well as by sustaining the specific forest structure. Applying the indicators of Climate Smart Forestry [81] to the specific plots resulted in generally high and stable over time smartness values [82] for selected indicators. Those selected indicators refer to the underlying criteria of sustainable forest management, global carbon cycle, health and vitality, productive function, and biological diversity [83]. In that sense, selection forestry seems to permit a high potential of adaptability and sustainable delivery of forest functions.

Conclusions
The successful reduction of SO 2 -emissions and recovery of tree growth and vitality demonstrate that human engagement in environmental policy can revive forest ecosystem health. Selection forests of N. spruce, s. fir, and E. beech that are common in the submontane to subalpine zone proved to be well adapted to environmental changes. Their combination of species diversity and structural heterogeneity can stabilize growth and promote climate smartness. The study object corroborates that uneven-aged, multi-layered mixed stands, which are the target stands of many present transitioning programs can result in rather stress resistant ecosystems and in stable productivity on the given sites. Next, steps should encompass analogous analyses regarding other kinds of stress scenarios, other species combinations, and site conditions. The study also reveals that for long-lived managed ecosystems, long-term observations are essential to understand underlying system processes extractable to improve sustainable forest management.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/f12070894/s1, Figure S1: Single tree periodical mean annual stem diameter growth of Norway Spruce, silver fir, and European beech over calendar year on the experimental plot FRY 129. Figure  S2: Species-specific mean of periodical mean annual stem diameter growth of Norway Spruce, silver fir, and European beech over diameter on the experimental plot FRY 129. Figure S3: Single tree periodical mean annual stem diameter growth of Norway Spruce, silver fir, and European beech over diameter on the experimental plot FRY 129. Figure S4: Species-specific periodic annual volume increment distribution over DBH-classes of the plot Freyung 129/22. Figure S5: Species-specific periodic annual volume increment distribution over DBH-classes of the plot Freyung 129/31. Figure  S6: Species-specific periodic annual volume increment distribution over DBH-classes of the plot Freyung 129/32. Table S1: Species-specific estimates for modelling tree height in dependence of stem diameter.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.