Can forest structural diversity be a response to anthropogenic stress? A case study in old-growth fir Abies alba Mill. stands

Abstract

Key message

From 1973 to 1991, Polish SO ₂ emissions above 3250 Gg/year resulted in a decline of fir Abies alba Mill. After stresses connected with these emissions, five main diameter at breast height (DBH) structural types were described. This type of heterogeneous forest structure is supposed to increase forest resistance and resilience to abiotic, biotic and anthropogenic disturbances.

Context

The analyses of forest structure are important under the current scenario of global change, since heterogeneous structures allow for better responses to disturbances. Forests with more complex structures should present greater vitality.

Aims

The main hypotheses were as follows: (1) the temporal changes of atmospheric SO₂ emissions caused (a) the abrupt changes in the tree DBH radial increment and (b) the death of fir trees; and (2) atmospheric SO₂ emissions and related fir decline and recovery processes ultimately result in the development of stands characterised by diverse DBH structures.

Methods

Radial growth trends of 49 and 215 fir trees in the older and younger generations, respectively, and 85 dead fir trees were evaluated. Using data collected in 32 stands, the DBH structural types were identified, and the shapes of these types were illustrated.

Results

The structural diversification of forest patches may influence forest resistance and resilience to disturbances; five main structural types were identified: ML1 and ML2 represent DBH distributions of multi-layered stands, and OS, TS1 and TS2 represent DBH distributions of one- and two-storied stands.

Conclusion

Structural diversity of forests was a response to SO₂ emissions; fir trees had the ability to increase their radial growth, although there were still high SO₂ emissions.

1 Introduction

Considering only the tree components of a forest, structural diversity can be observed in the heterogeneity of forest patches, primarily as a variation in tree diameter at breast height (DBH) distributions. Structural diversity plays an important role in the analysis of biological processes and functions, e.g. in explaining forest dynamics (Korpeľ 1995). Diameter structures are often used to characterise developmental stages and phases (sensu Korpeľ 1995); this knowledge is important to understanding the changes taking place in the forest developmental cycles (Coomes and Allen 2007; Ashton and Kelty 2017).

The multifactorial nature of tree mortality means that simple associations among variables such as various anthropogenic factors, or climate changes and forest dynamics, may be often unrealistic (Lloret et al. 2012). The number of studies closely documenting forest shifts is still too low to provide a solid framework for studying changes in forest dynamics. The response of trees to various factors is highly variable and specific, especially if these factors do not exceed critical values (Lines et al. 2010). Forest ecosystems use ecophysiological and demographic-stabilising mechanisms, which can compensate for dying. However, such mechanisms should not underestimate the potential of these events to induce dramatic changes, especially if they are more frequent (Lloret et al. 2012). In addition to looking for thresholds above which impact factors cause significant changes in forests, research on the global change drivers should be directed to the study of forest structures ensuring the greatest forest stability. Heterogeneous forests are expected to be more resistant and resilient, owing to, for example, forest patches with uneven-aged structures allowing for the permanence of high regeneration cover in the stand and to complement in response to disturbance (e.g. O’Hara 2006; Lafond et al. 2013). The impact of SO₂ pollution was usually greater in vertically and horizontally homogeneous forests with a greater tree density (Vacek and Lepš 1996). This observation confirms the beneficial effect of thinnings shaping the forest structure on the vitality of stands exposed to SO₂ emissions (Rydval and Wilson 2012).

The fir Abies alba Mill. is an ecologically and economically important forest tree species, especially in lower-mountain forests. In Central Europe, the death of firs has occurred since approximately 1950; forest stands dominated by firs have often disappeared. The mortality of fir trees and DBH radial increment depression reached their maxima between 1970 and 1990. These phenomena were usually not mono-causal, but atmospheric SO₂ emissions were probably the most important factor (e.g. Pitelka and Raynal 1989; Becker et al. 1990; Ellenberg 1996; Wilson and Elling 2004). Studies have focused on the impact of SO₂, nitrogen deposition, and tropospheric O₃, as well as meteorological and hydrological data (e.g. Elling et al. 2009; Paoletti et al. 2010; Bošeľa et al. 2014). In Central Europe, between 1970 and 1990, the level of anthropogenic SO₂ emissions was among the highest on the continent (e.g. Cienciala et al. 2016; Altman et al. 2017). In Poland, anthropogenic SO₂ emissions peaked in the 1980s and decreased after 1990 (Smith et al. 2011a, b).

A tree diameter distribution is a structural feature employed to describe the forest heterogeneity (e.g. Westphal et al. 2006; Wang et al. 2009). The biological reason behind the shape of the DBH distributions is linked to the reduction rate with an increasing DBH bin (Diaci et al. 2011). In multi-layered forests, the reduction rate can be constant (negative exponential distribution), decreasing (negative power function) or variable (rotated sigmoid shape) (overview Westphal et al. 2006). The shape of DBH distributions in mixed, multi-layered forests can also be multimodal and irregularly descending (Podlaski 2010). In this type of forest, DBH structures are usually too complex to be modelled by a single function (Zasada and Cieszewski 2005). Among the parametric methods, the most precise and flexible are the mixtures of probability density functions (e.g. Podlaski 2011a, b). Parameter estimation can be especially problematic with mixed distributions. One of the ways to solve this problem is to use the gamma shape mixture (GSM) model and a general Bayesian approach for estimating the unknown parameters (Venturini et al. 2008).

The main hypotheses were as follows: (1) the temporal changes of atmospheric SO₂ emissions caused (a) the abrupt changes in the tree DBH radial increment and (b) the death of fir trees and (2) atmospheric SO₂ emissions and related fir decline and recovery processes ultimately result in the development of stands characterised by diverse DBH structures. This structural differentiation plays a key role in the functioning of forest ecosystems; therefore, recommendations for management strategies are also discussed.

2 Materials and methods

2.1 Study area

The sample plots were located in the Świętokrzyski National Park (Central Poland; geographical coordinates: 50° 50′–50° 58′ N, 20° 48′–21° 08′ E) (Fig. 1). In the study area, there are three strictly protected forest reserves, created in 1920s. The purpose of their creation was to protect old-growth forests that include fir. Here, the term ‘old-growth’ describes ‘stands composed entirely of trees which have developed in the absence of allogenic processes’ (Oliver and Larson 1996). In the above definition, the ‘allogenic processes’ refers to large-scale disturbances (Wirth et al. 2009). These stands are also relatively old, with their structural and compositional features witnessing self-replacement (see also Mosseler et al. 2003).

Fig. 1

A grid of the System of Information on Natural Environment (SINUS) covering the whole of Poland; P₀ tessellation divides the country using 10′ × 10′ longitude and latitude blocks (P₀ coded 5041, 5042 and 5043 include the area of the Świętokrzyski National Park). Study area: a Święta Katarzyna forest section, b Święty Krzyż forest section, c Chełmowa Góra forest section

The main soil types are Distric Cambisol and Haplic Luvisol (sub-types according to FAO, ISRIC, and ISSS 2006). The mean annual temperature is 5.9 °C, the growing season is about 182 days (Olszewski et al. 2000). Mean annual precipitation ranges from 700 to 850 mm (Olszewski et al. 2000). The study area is dominated by the following plant associations: Dentario glandulosae-Fagetum Klika 1927. em. W. Mat. 1964, and Abietetum polonicum (Dziub. 1928). Br.-Bl. & Vlieg 1939 (in the Święta Katarzyna and Święty Krzyż forest sections), as well as Tilio-Carpinetum abietetosum with Larix decidua Mill. subsp. polonica (Racib.) Domin (in the Chełmowa Góra forest section); nomenclature after Matuszkiewicz (2008).

2.2 Data collection

2.2.1 Sulphur dioxide data

Emissions data for Poland were obtained from Smith et al. (2011a, b). Annual estimates of anthropogenic SO₂ emissions have been constructed spanning 1850–2005 using a bottom-up mass balance method, calibrated to country-level inventory data (Smith et al. 2011a, b). The western, south-western and southern winds prevail in the investigated area. The direction of the dominant winds suggest that large industrial regions situated to the west and south of the Świętokrzyski National Park played a significant role in shaping radial growth of fir and the structure of stands. Their combined share in total SO₂ emissions in Poland exceeded 30% during the period considered; in addition, transborder emissions were deposited (Wertz 2012).

2.2.2 Sample points sampling

A grid of the System of Information on Natural Environment (SINUS) consists of ellipsoidal trapezoids 10′ × 10′ (blocks, marked with symbol P₀) (Fig. 1). The P₀ blocks are divided into 54 P₁ sub-blocks (100.00″ × 66.67″; 6 × 9 = 54); each P₁ sub-block is divided into 4 P₂ sub-blocks (50.00″ × 33.33″; 2 × 2 = 4); each P₂ sub-block is divided into 4 P₃ sub-blocks (25.00″ × 16.67″; 2 × 2 = 4); etc. (Fig. 2). The blocks P₀, comprising the area of the Świętokrzyski National Park, coded 5041, 5042 and 5043, have the dimensions 11,697 × 18,537 m; P₁ 1950 × 2060 m; P₂ 975 × 1030 m and P₃ 487.5 × 515.0 m. These blocks and sub-blocks were drawn on the management maps of the Święta Katarzyna, Święty Krzyż and Chełmowa Góra forest sections (for details, see Podlaski 2005). The sample points were selected in the P₂ (Święta Katarzyna and Święty Krzyż) or P₃ (Chełmowa Góra) sub-blocks, and a simple random sampling with replacement (SRSWR) was used (Cochran 1977). On the map, the position of each point, using 2.5-m stretches, was marked on the x-axis (P₂: 975 m = 390 stretches, P₃: 487.5 m = 195 stretches) and on the y-axis (P₂: 1030 m = 412 stretches, P₃: 515.0 m = 206 stretches) (Fig. 2). Using a random number generator, two numbers were drawn. The first number determined the position of a 2.5-m stretch on the x-axis, and the second number determined the position on the y-axis. If either number was greater than its respective maximum border value, the pair was rejected, and the next pair was considered. The intersection of the middle points of the two sampled stretches (on the x- and y-axes) determined the coordinates of a sample point (Fig. 2). In each P₂ or P₃ sub-block, 10 sample points were selected at random according to the sampling scheme described. In marginal sub-blocks, a proportionally smaller number of sample points was drawn. A total of 251 permanent sample points were selected, traced out in the field and marked in the stands (Fig. 3).

Fig. 2

A division of the block P₀ into sub-blocks P₁, P₂ and P₃, and the method of random selection of sample points and fir trees from the older and younger generation, dead fir trees and plots. An exemplary sampling in the sub-block P₃, when, in the surrounding area, the sample points are a fir trees from the older and younger generation, and lying dead fir trees; b standing dead fir trees; c fir trees from the older and younger generation, and the conditions required to create the sample plots are met

Fig. 3

Study area and sampled sample points. a Święta Katarzyna forest section. b Święty Krzyż forest section. c Chełmowa Góra forest section

2.2.3 Sample fir tree and sample plot selection; tree and stand level data

In the area surrounding each sample point, sample fir trees were selected, one tree from the older generation (trees aged from 136 to 300 years) and one tree from the younger generation (trees aged from 45 to 135 years) (Fig. 2). These trees represented Kraft’s second class in one-story stands or the upper canopy layer (> 2/3 h_max; 100 according to the International Union of Forest Research Organisations, IUFRO) in stands of complex structure. Furthermore, one dead fir tree (standing or lying on the ground) was chosen. Finally, 49 fir trees of the older generation, 215 fir trees of the younger generation, and 85 dead fir trees were selected according to this pattern. The number of selected sample fir trees was lower than the number of sample points because there were no suitable fir trees near some of the sample points.

Two increment cores per tree bored to the pith were taken from each sample fir tree (in 1995—from all trees, in 2005, 2015, and 2017—only from live trees). For most trees, especially for dead trees, due to wood decay, it was not possible to obtain complete, undamaged increment cores.

When all the sample fir trees had been selected, sample plots were chosen. Based on forest management maps, sample points located in stands with a share of fir trees larger than 80% were selected. Next, during the field inspection, sample points situated too close to (1) open areas (forest glades, large gaps); (2) borders between forest patches of the various vertical stand structures and (3) paths and forest roads were rejected. Finally, 15 sample points in the Święta Katarzyna forest section and 17 in the Święty Krzyż forest section were selected (altogether 32 sample points). The sample points indicated the central points of the circular sample plots from 0.20 to 0.35 ha. While the plots were created, the stages and phases of forest development were determined. The main criteria used to determine the stages and phases were the following (Korpeľ 1995): (1) vertical stand structure, (2) tree age distribution and (3) tendency of volume increment (increasing or decreasing). Tree height had been approximately determined, and vertical stand structure was identified based on the number of trees in three main stories: (1) upper (> 2/3 h_max), (2) middle (between 1/3 h_max and 2/3 h_max) and (3) lower (< 1/3 h_max). Tree age was measured in the individual layers in the increment cores from several trees in each main story. The tendency of the volume increment was derived from vertical stand structure and age distribution. The size of the sample plots was chosen in such a way that each of the plots represented a homogenous forest patch of the same vertical stand structure.

In each sample plot, DBH data were collected for every live tree above a minimum DBH of 6.9 cm. Measurements of trees below this threshold add field time and provide relatively minimal significant ecological data at the level of the stand.

2.3 Data analysis

2.3.1 Radial growth data

Tree ring widths were measured with a precision of 0.01 mm (two instruments were used—GP-3 and CODIMA). The individual tree ring width series were cross-dated. As quality criteria, the t value (Baillie and Pilcher 1973) and the collinearity of increments (Gleichlaüfigkeit; Buras and Wilmking 2015) were considered. The matching of chronologies was accepted as reliable when it reached a minimum t value of 3.5 and a minimum collinearity of 70% for a 50-year overlap. In the case of the dead sample fir trees, due to wood decay, only 28 chronologies were reliable. The radial increment trends of the investigated sample fir trees from the older and younger generation were analysed using the generalised additive model (GAM):

$$ T{R}_{\bullet t}={\beta}_0+s(t)+{\varepsilon}_t $$

(1)

where TR_•t ≡ TR_OGt andTR_•t ≡ TR_YGt are the mean tree-ring width (in mm) of the sample fir trees from the older and younger generations, respectively, at the year t and TR_•t~normal, s(•) is the smooth function, β₀ is an intercept and ε_t is the residual error. The models were estimated using the restricted maximum likelihood (REML). The number of knots was 20. The thin plate regression splines tp were employed. The tp splines have knots that differ significantly from conventional knots, and thus, a truncated eigen-decomposition is used to achieve the rank reduction (for detailed information, see Wood 2017).The radial growth periods of significant environmental change were identified using the first derivative of the fitted trend. Derivatives of the fitted spline are not easily available analytically, but they can be estimated using the method of finite differences. The GAM model detects periods of significant change as those time points where the Bayesian confidence interval on the first derivative does not include zero (Simpson 2018). These intervals can be obtained by simulation from the posterior distribution of the first derivative. A (1 − α) 100% confidence interval contains in its entirety 1 − α of all random draws from the posterior distribution (Simpson 2018). Such intervals are known as simultaneous intervals; here, α = 0.05 (for detailed information, see Wood 2017).

2.3.2 Sulphur dioxide emissions—radial growth relationships

In Poland, from 1965 to 1995, SO₂ emissions reached high values, above 2000 Gg/year, and from 1973 to 1991 peaked above 3250 Gg/year (Smith et al. 2011a, b; see also arrows and vertical, dashed lines in Figs. 4, 5 and 6). The relationships between tree-ring widths and anthropogenic SO₂ emissions were analysed for 1964–1996 using the GAM model:

$$ \log \left(T{R}_{\bullet it}\right)={\beta}_0+s\left({\mathrm{SDE}}_t\right)+{\varepsilon}_{it} $$

(2)

where TR_•it ≡ TR_OGit and TR_•it ≡ TR_YGit are the tree-ring width (in mm) of the sample fir tree i from the older and younger generations, respectively, at the year t and TR_•it~gamma, SDE_t is the SO₂ emissions in Poland (in Gg) at the year t, s(•) is the smooth function, β₀ is an intercept and ε_it is the residual error. The smoothing parameter estimation problem and the choice of the number of knots were solved using the generalised cross-validation (GCV) criterion. The thin plate regression splines tp were employed.

Fig. 4

a Polish anthropogenic SO₂ emissions, based on the dataset from Smith et al. (2011a, b). The distance between arrows = high values of emissions, the distance between two vertical, dashed lines = maximum values of emissions, above 3250 Gg/year. b The mean radial increment and the GAM radial growth trend of the older (small dots and a thin, dashed curve) and younger generation (large dots and a thick, solid curve) of fir. c The mean radial increment of dead sample fir trees selected in 1995, a region between two vertical, dashed lines = the period from 1975 to 1980, when all these trees died

Fig. 5

First derivatives of the GAM radial increment trends calculated using finite differences (thick, solid lines), 95% simultaneous confidence intervals (thin, dashed lines), and 500 random draws from the posterior distribution of the first derivatives (sets of thin, solid lines) for the older (a, b) and younger (c, d) generation of fir. The distance between arrows = the simultaneous confidence interval does not include 0, the GAM radial increment trends detect periods of significant temporal change (see also Fig. 4b)

Fig. 6

The GAM model for the older (a) and younger (b) generations of fir. The response is the tree ring width (TR_•it ≡ TR_OGit and TR_•it ≡ TR_YGit), which is expressed in mean deviation form (the value of zero on the y-axis is the mean of the response) and is centred; thus, each panel represents how the response changes relative to its mean with changes in emissions. Panels show the estimated effects (solid curves), with 95% confidence limits (dotted lines). The number in each y-axis caption is the effective degrees of freedom. A vertical, dashed line = the border between high and maximum values of emissions

2.3.3 Identification of DBH structural types

Choosing the appropriate clustering variables is extremely important when defining homogeneous DBH structures based on empirical DBH distributions. There are different statistics characterising empirical DBH distributions. A new approach is to use suitable sums of mixture weights and the parameter θ of the GSM model (1/θ is a scale parameter) as clustering variables in the hierarchical cluster analysis (HCA) (Venturini et al. 2008). The sums of mixture weights and the parameter θ are related to the proportion of trees in appropriate DBH classes (Fig. 7).

Fig. 7

Illustration of the relationship between the sums of mixture weights, the parameter θ (1/θ is a scale parameter) and the proportion of trees in appropriate DBH classes (see Eq. 3). a The GSM model fits the empirical DBH data; this model consists of 250 gamma distributions j = 1, …, J, J = 250. b The share of each gamma distribution j in the GSM model is determined by the corresponding mixture weight (proportion) π_j, and therefore, the share of trees in the DBH class is closely related to the specified sums of mixture weights and the parameter θ

The GSM model is defined as (Venturini et al. 2008):

$$ f\left(x|{\pi}_1,...,{\pi}_J,\theta \right)=\sum \limits_{j=1}^J{\pi}_j{f}_j\left(x|\theta \right) $$

(3)

where J is the number of mixture components (known and fixed); π₁, ..., π_J are mixture weights (proportions) (unknown) and 1/θ is a scale parameter for the whole GSM model (unknown). The gamma distribution f_j(x| θ) has a probability density function (PDF) given by:

$$ {f}_j\left(x|\theta \right)=\frac{\theta^j}{\varGamma (j)}{x}^{j-1}{e}^{-\theta\;x} $$

(4)

Each gamma distribution in the GSM model is indexed by a component-specific shape parameter (j) and has a single-scale parameter (1/θ).

A general Bayesian approach for estimating the unknown parameters of the GSM model is used. The posterior distribution is estimated using a Gibbs sampler, and the parameter θ is derived analytically through integration (Venturini et al. 2008). Fitting DBH data with the GSM model requires three hyperparameters: the number of components, J, and the α and β from the conjugate prior on θ. It was assumed that the value of J = 250 and that the weight of the prior information ω = 0.35 (ω values between 0.2 and 0.5 are usually chosen). With these assumptions, the α and β values were calculated (for detailed information, see Venturini et al. 2008). Generally, the GSM model is very useful for modelling differentiated data sets (Venturini et al. 2008; Podlaski 2017).

The HCA was used with the Manhattan measure and Ward’s minimum variance agglomeration method (Murtagh and Legendre 2014). To group the empirical DBH distributions representing stand structural diversity, the following 11 clustering variables were employed: sums of mixture weights of the GSM model (π₁, ..., π_J; J = 250) where var1 was from j = 1 to j = 25, var2 from j = 26 to j = 50, var3 from j = 51 to j = 75, var4 from j = 76 to j = 100, var5 from j = 101 to j = 125, var6 from j = 126 to j = 150, var7 from j = 151 to j = 175, var8 from j = 176 to j = 200, var9 from j = 201 to j = 225, var10 from j = 226 to j = 250 and var11 ≡ the parameter θ. To avoid the effect of considerable differences between the scales of the variables from var1 to var10 and var11, the values were standardised to z scores, with a mean of 0 and a standard deviation of 1. The homogeneity among stand clusters was evaluated using two non-parametric procedures: analysis of group similarities (ANOSIM; Clarke 1993) and multiple response permutation procedures (MRPP; Mielke 1991). The ANOSIM R statistic is scaled into the range from − 1 to + 1 (R = 0 indicating independence). The MRPP A statistic is based on the observed average of within-group dissimilarity d and its expected value E (d) assessed from permutations. The statistic A = 1 when all items are identical within groups; A > 0.3 is relatively high in ecology (Warton et al. 2012).

2.3.4 Modelling DBH structural types by finite mixture models

To illustrate the shapes and to emphasise the differences between the DBH distributions of distinguished DBH structural types, finite mixture models consisting of two or three gamma distributions were employed. The appropriate PDF can be written as:

$$ {f}_X\left(x|\boldsymbol{\Psi} \right)={\pi}_1{f}_1\left(x|{\boldsymbol{\uptheta}}_1\right)+\left(1-{\pi}_1\right){f}_2\left(x|{\boldsymbol{\uptheta}}_2\right) $$

(5)

and

$$ {f}_X\left(x|\boldsymbol{\Psi} \right)={\pi}_1{f}_1\left(x|{\boldsymbol{\uptheta}}_1\right)+{\pi}_2{f}_2\left(x|{\boldsymbol{\uptheta}}_2\right)+\left(1-{\pi}_1-{\pi}_2\right){f}_3\left(x|{\boldsymbol{\uptheta}}_3\right) $$

(6)

where θ_i = (α_i, β_i, γ) are the parameters of the gamma distribution f_i(•); the parameter γ sets as fixed, γ = min  − 0.1, where min is the minimum value of the DBH of all trees in the plot investigated; π_i is the weight (fraction); Ψ is a complete parameter set for the overall distribution; and i = 1, 2, or i = 1, 2, 3.

To maximise the likelihood function, the multistart method as well as the EM algorithm with the NT method were employed. The computational procedures were implemented in R (R Core Team 2017); the GSM (Venturini 2014), vegan (Oksanen et al. 2013), and mixdist (Macdonald and Du 2004) packages of R were also used.

3 Results

The following stand and DBH distribution variables were calculated for each sample plot: share of fir assessed on the basis of a basal area (s_f), number of trees per hectare (N), stand basal area (G), minimum DBH (d_min), average DBH ($\overline{d}$), maximum DBH (d_max), 1 quartile DBH (d_0.25), 2 quartile DBH (≡median) (d_0.50) and 3 quartile DBH (d_0.75). Summary statistics, including the minimum, mean and maximum values and SD of the variables are shown in Table 1.

Table 1 Stand variables of the data set used for analysing the degree of DBH structural diversity

Anthropogenic SO₂ emissions in Poland have steadily increased, with a peak in the 1980s and a decrease after 1990 (Fig. 4). The emissions increased in rough proportion to the activity levels. The change after 1990 is connected to the removal of SO₂ from various waste gas streams including those from incineration, coal burning and metal smelting operations. The radial growth trends and the relationships between tree ring widths and SO₂ emissions for fir trees were similar, but the greater and longer decrease occurred in fir trees from the older generation compared to that of the fir trees from the younger generation (Figs. 4, 5 and 6). The GAM models used to analyse the radial increment trends highlight that a growth decrease during the highest SO₂ emissions was significant and exceptional in relation to the entire examined period (Figs. 4, 5). The greatest decrease in the DBH increment occurred from 1965 to 1982 for fir from the older generation and from 1971 to 1978 for fir from the younger generation (Fig. 5). Polish anthropogenic SO₂ emissions below 3250 Gg/year (more specifically, from 2000 to 3250 Gg/year of SO₂) caused a slight reduction in fir radial increments (Figs. 4 and 6). Increasing these emissions above 3250 Gg/year resulted in fir mortality and in a strong reduction in tree radial growth (Figs. 4 and 6). The GAM models used to analyse the relationships between tree ring widths and anthropogenic SO₂ emissions show that the shape of these relationships is nonlinear; at first, the form of the curves is close to horizontal and then the curves drop slightly. After exceeding the SO₂ value of 3250 Gg/year, the curves sharply decrease (Fig. 6). The last segments of the curves are nearly linear (Fig. 6).

In the studied fir stands, after stresses connected with SO₂ emissions, five main DBH structural types were identified (A–E; Fig. 8). The ANOSIM analysis showed that each DBH structural type was different from the other types (P = 0.001) for five identified clusters R = 0.763. The MRPP results were significant (P = 0.001), and the A statistic values were 0.39, 0.41, 0.53 and 0.57 for two to five identified clusters, which indicated an increasing similarity in the DBH structures within clusters. These tests note substantial differences between the identified clusters and high levels of homogeneity within these clusters.

Fig. 8

Hierarchical cluster analysis (HCA) of investigated stands. To group the stands, suitable sums of mixture weights and the parameter θ of the gamma shape mixture (GSM) model (1/θ is a scale parameter) were used. Cluster A = the DBH structural type ML1. Cluster B = the DBH structural type ML2. Cluster C = the DBH structural type OS. Cluster D = the DBH structural type TS1. Cluster E = the DBH structural type TS2

The A and B clusters (multi-layered stands—DBH structural types ML1 and ML2, respectively) represent DBH distributions that consist of three sub-populations, and therefore, they were adequately approximated by the three-component gamma mixtures (Figs. 8 and 9). These DBH structural types are characterised by a significantly higher proportion of thin and medium thickness trees compared to that of thick ones (Fig. 9). Furthermore, the DBH distributions representing the thin and medium thickness trees are partially overlapping (for the type ML1; Fig. 9a), or they are clearly separated (in the case of the type ML2; Fig. 9b). The range of DBHs is also different, at 7–111 cm for ML1 and 7–67 cm for ML2 (Table 1, Fig. 9). These structures are characteristic of stands composed of trees of all ages with multi-layered canopies and for complex and stratified stands with fir trees aged 45 to 300 years. Natural processes of fir regeneration are less intensive in stands of ML1 compared to stands of ML2 (see the share of the thin fir trees from the smallest DBH classes; Fig. 9).

Fig. 9

The DBH structural type ML1 (a) and ML2 (b) of the multi-layered fir stands, see also the A and B clusters in Fig. 8. The approximations of the summarised empirical DBH data using the gamma model (a thick curve) with three components (three thin curves). The empirical DBH distributions of the sample plots belonging to the appropriate cluster and the least (1) and most (2) deviating from the summarised distribution

The C–E clusters (one- and two-storied stands—DBH structural types OS, TS1 and TS2, respectively) represent DBH distributions that consist of two sub-populations, and therefore, they were adequately approximated by the two-component gamma mixtures (Figs. 8 and 10). The DBH structural type OS does not have two maxima (the DBH distributions representing two main sub-populations are partially overlapping). The DBH structural types TS1 and TS2 have two maxima; the first maximum is global for TS1, and the second one is global for TS2 (Fig. 10). These structures are typical of one-storied, relatively uneven-aged stands (type OS) and two-storied stands with various shares of trees in the two main layers (types TS1 and TS2). Generally, the curves of mixture models reach a global maximum (connected with the youngest tree cohorts), and then they show a decrease (Figs. 9 and 10a, b).

Fig. 10

The DBH structural type OS (a), TS1 (b) and TS2 (c) of the one- and two-storied fir stands; see also the C–E clusters in Fig. 8. The approximations of the summarised empirical DBH data using the gamma model (a thick curve) with two components (two thin curves). The empirical DBH distributions of the sample plots belonging to the appropriate cluster and the least (1) and most (2) deviating from the summarised distribution

4 Discussion and conclusions

The SO₂ emissions explained a larger part of variability in the radial growth of fir and changes in the structure of stands, but various abiotic (e.g. temperature, water deficiency) and biotic (e.g. insect outbreaks) factors were significant as well (Jaworski et al. 1995). The decline of fir was also varied depending on altitude and orography (Jaworski and Skrzyszewski 1986). In Poland, the investigated changes were affected more by SO₂ emissions than by any other factors; SO₂ emissions initiated the sequence of events leading to the decline of fir (Jaworski 1991). In Europe, there have been several observations suggesting that SO₂ emission played a key role in the decline of fir and other tree species in the 1960s and 1970s of the twentieth century (e.g. Pitelka and Raynal 1989; Wilson and Elling 2004; Elling et al. 2009).

The patches of different DBH structures, lying side by side, were created as a result of disturbances. Fir trees growing in these new conditions had the ability to increase their radial increment. Relatively small, even-aged patches surrounded by uneven-aged structures, compared to even-aged stands, also showed greater resistance and resilience to still high SO₂ emissions. Transformation from a mono-structural form to a complex forest, consisting of even- and uneven-aged patches of various sizes, can significantly affect the stability of the forest.

The phenomenon of fir decline caused the creation of untypical periods, significantly different from the stages and phases of undisturbed developmental cycles. In investigated forests, there are stages and phases, which have not been described or explained by typical developmental cycles (Podlaski 2008). In temperate forests with fir, a natural disturbance regime is usually characterised by small-scale disturbances, while intermediate-scale disturbances are episodic. This is why authors presenting variations on the mosaic cycle concept of forest dynamics (e.g. Watt 1947; Leibundgut 1982, 1993; Mueller-Dombois 1987; Koop 1989; Oldeman 1990; Korpeľ 1995) did not consider these disturbances. The intermediate-scale disturbances are mainly caused by wind, snow or ice storms, but fir trees are rather resistant to these. The reconstruction of the real dynamics of investigated forests is possible only when typical developmental cycles are modified (Podlaski 2008). The identified DBH structures are specific to intermediate-scale disturbances. The main consequence on forest dynamics of the decrease in SO₂ emissions at the end of the twentieth and the beginning of the twenty-first century is the return to typical developmental cycles.

An abrupt shift in forest dynamics would occur when the relative abundance of fir trees significantly changes over the range of natural temporal variability (Lloret et al. 2012). The SO₂ emissions resulted in higher mortality rates of fir than natural variability, which was later compensated by a higher survival of the remaining population and by higher regeneration in gaps. The heterogeneous forest structure supports the mortality compensation after the perturbation.

To appoint homogeneous DBH structures of forests containing tree generations formed by firs of various ages, new types of clustering variables are required. A new approach is to use the sums of mixture weights and the parameter θ of the GSM model (1/θ is a scale parameter) fitting empirical DBH distributions as clustering variables. The GSM model employs a mixture of gamma density functions (Venturini et al. 2008). A general Bayesian approach allows the creation of a flexible model characterised by a single-scale parameter for all the gamma components and the ordinary set of mixture weights (Venturini et al. 2008). This method significantly improves the predictive performance in estimating local extremes, connected with tree generations, compared to that of standard approaches employing mixture distributions with a few components. A particularly important advantage of the GSM model is the possibility to use a great number of mixture components; e.g. 250 or more (Venturini et al. 2008). This methodology can be used to better identify or imitate complex forest structures.

The investigated forests were characterised by the heterogeneous DBH distributions. It is interesting to compare DBH structures of fir-only stands and mixed forests with fir and beech, Fagus sylvatica L. (Podlaski 2010). Fir and beech are classified as shade-tolerant canopy species, but fir is more tolerant than beech is. It is mainly for this reason that there is a greater number of age generations in fir stands than that in mixed stands with fir and beech. Fir trees can respond with increased growth when cyclic changes occur in their competitive environment. This feature is one of the main factors enabling the development of multi-layered stands with fir, which are characterised by very diverse DBH structures.

Forest patches possessing variable DBHs are likely to also have their own set of habitat requirements for water and soil nutrients. A multi-layered forest structure specifically allows a more efficient utilisation of light, water and soil nutrients at the stand level; as a result, a multi-layered forest structure increases aboveground C storage (Zhang and Chen 2015; Ali et al. 2016). In temperate forests, the structural diversity, rather than species diversity, increases aboveground biomass and C storage over time (Szwagrzyk and Gazda 2007; Ali et al. 2016).

The presence of trees of different sizes creates structural diversity. The number of trees, basal area and biomass are usually disproportionately distributed within DBH classes. If the large-diameter trees have been removed or are dead, then the smaller-diameter trees often take a shorter time to advance to the upper story. The abundance of small-diameter trees suggests that the forest could have the ability to respond to disturbances. The demography of small-diameter trees is important for predicting the long-term changes in the structural diversity and in the overall ecosystem functions; however, forest structural properties very strongly depend on the locations of the largest 1 to 2% of the trees (Lutz et al. 2013). Large trees play a key role in the functioning of the ecosystem (Bošeľa et al. 2013). Increasing structural diversity makes it unlikely that disturbances will affect all trees in a stand, making forests more resilient (Franklin et al. 2000; Gustafsson et al. 2010; O’Hara and Ramage 2013).

In Central Europe, simple silvicultural systems (clear-cuttings, shelterwood systems with short regeneration periods) were actively implemented until the 1990s when deliberate forest management shifted towards an ecosystem management approach. This transition enabled the introduction of close-to-nature silviculture. The main idea of this silvicultural concept is to mimic the structure and functions occurring in unmanaged, protected forests; its fundamentals were presented by Karl Gayer in the nineteenth century (Gayer 1886; Schütz 1994).

Structural diversity can mainly be achieved using systems with long regeneration periods and transformation cuts (Schütz 2001). Transformation from a regular forest to an irregular form, which is stable, long-lasting and diverse, follows a set sequence of stages (Schütz 2001): (1) the stage of differentiation, in which the main aim is to promote each existing valuable element, which ensures structural development; (2) the stage of regeneration promotion, in which the principal focus is on favouring new decentralised regeneration groups; (3) the stage of structural development, in which the focus is to achieve good horizontal and vertical distribution of structural elements and (4) the stage of structure achievement, in which the focus is to achieve vertical individualisation of the remaining groups. Maintaining structural diversity is connected to close-to-nature silviculture, including three main systems: (1) single-tree selection, (2) group selection and (3) a shelterwood system with long regeneration period. In forests with fir, close-to-nature silviculture should implement appropriate cutting methods that imitate disturbances, especially at the fine scale, which, in the long term, shape complex DBH structures. Out of the noted three main systems, the shelterwood method is the least useful, because in its classic form, it tends to produce even-aged mono- or two-layered stands. In Central Europe, most of the stands with fir are still managed using the shelterwood system; therefore, it is necessary to introduce an irregular form. The presented recommendations for management strategies may positively influence the structural diversity of the fir stands.