The Mechanical Stability of Pure Norway Spruce Stands along an Altitudinal Gradient in the Czech Republic

Abstract

Norway spruce stands are established and managed along various site conditions in central Europe. Currently, spruce often grows at locations outside of its ecological optimum, resulting in extensive damage elicited by harmful abiotic and biotic factors, which relatively shortens the time to change this adverse status in the adaptation frame by foresters. Except for the rapid change in species composition through clear-cuts, another way is possible, i.e., stabilising current (especially young) spruce stands to extend the time required to implement adaptation measures. The assumption that different site conditions will have to be respected as part of this adaptation was confirmed by our study based on NFI data of the Czech Republic. A semiparametric generalized linear model (GAM) was used to model the relationship between the height-to-diameter ratio and forest stand age, differentially considering particular forest vegetation zones. Spruce stands with lower elevations attain a lower stability (expressed by their height-to-diameter ratio; HDR) than those in the mountains. The HDR culminated in lower and middle altitudes in the first half of the rotation period, representing the most critical timing and effectivity of silvicultural measures. Contrary to previous findings, we found higher HDR values at nutrient-rich sites than those at acid ones, especially up to 50–60 years old. Therefore, more research should be devoted to the issue concerning the same thinning regime under different site conditions.

1. Introduction

European forest ecosystems are threatened by abiotic (storm, snow, frost, drought, and fire, e.g., [1,2,3,4]) and biotic factors (pest insects, diseases, and ungulates, e.g., [3,5,6,7]) under ongoing global climate change conditions, characterised by increasing mean annual air temperatures and changes in the distribution of precipitation within the growing season [3,8], resulting in a decrease in the provision of multiple ecosystem services [9]. Therefore, the mechanical stability of forest ecosystems, as one of the most productive components in the planetary biosphere, plays a pivotal role [10,11,12].

Natural disorders significantly affect forest structure and dynamics [13,14]. Forest vulnerability to natural disturbances has recently increased, with around 58% of all European forests being at risk of biomass loss [15]. This represents 39 million hectares of anthropogenically or naturally disturbed forests (in the 1986–2016 period), corresponding to 17% of the total European forest area [16], with an evident increasing trend in recent decades [17]. Specifically, in the Czech Republic, salvage cuttings removing damaged timber accounted for 86.9% of the total felling volume (from 28.72 million m³) in 2020 [18]. Windstorms elicit significant disturbances in all European forest ecosystems [11,19], resulting in severe economic and ecological consequences [20,21] and representing more than 50% of the volume of all the damaged forests in Europe [22]. These disturbances significantly affect forest productivity, carbon storage, and biodiversity [23,24,25,26]. On the contrary, snow damage routinely occurs in northern European spruce stands, mainly in mountainous areas [27]; however, it is complicated to differentiate snow damage from disturbances caused by wind, because both usually coincide [2].

European conifers (gymnosperms), mainly Norway spruce [28], are less stable than angiosperms [29]. Nonetheless, comparison is tricky, because conifers occupy sites that are usually inappropriate for many broad-leaved tree species (wind exposure, high altitudes, shallow soils, and high water tables; [12]). In addition, forest management can significantly enhance mechanical stability and persistence against forest stand disturbances [11] through silvicultural treatments [17,30,31,32,33]. The tending of young stands positively affects their resistance to wind and snow at the tree and stand levels [30,34,35,36,37,38,39,40,41]. This has also been confirmed for other conifers (e.g., [42,43,44]) or broad-leaved tree species [45,46]. Generally, the mechanical stability of forest stands against abiotic factors and stem radial growth is related to crown development [47,48,49]. In Norway spruce stands, radial increment and resistance to abiotic and biotic factors substantially increase after reductions in stand density by thinning (e.g., [50]), especially in pole stage stands [30]. However, heavy silvicultural measures applied in pole or older stands are connected to the risk of stand disintegration, especially if the treatments are the first in this growing phase [51,52], which also applies to pine stands [31]. Since the main risk factors affecting mechanical stability are tree species, tree height, and the diameter at the breast height (i.e., at 1.3 m above ground level, DBH [51,53]), the slenderness ratio, defined as the height-to-DBH ratio (HDR), is a crucial parameter for evaluating the mechanical stability at the tree and stand levels [5,54,55,56]. Generally, intensive stand tending can amplify a favourable HDR value. However, the damage risk may be greater several years after the interventions (before the HDR value decreases [57,58]). Therefore, in-time stand tending (no later than when a 7 m stand dominant height is reached [59]) is necessary to attenuate HDR culmination and lower the tree’s centre of gravity in a young growing phase [57]. Trees with lower HDR values withstand natural hazards, such as wind and snow damage, because of their more extended live crown and lower located centre of gravity [60]. Heavy thinning is more effective than mild thinning, but the timing of this intervention is crucial. The earlier the thinning application, the more efficiently it affects the HDR. Once the dominant heights exceed about 20 m, the thinning effect is comparatively slight [61].

On the contrary, trees with higher HDR values are more slender and less stable, resulting in a higher vulnerability to abiotic damage [62]. Generally, HDR values below 80 indicate a sufficient tree mechanical stability against external abiotic factors. In contrast, HDR values exceeding 100 signify specimens that are significantly vulnerable to natural hazards [63].

Although Norway spruce (Picea abies L. Karst.) is immensely affected by abiotic damage [64], followed by bark beetle (Ips typographus L.) outbreaks in the Czech Republic nowadays [65], its representation will be relatively high with significant economic value in the ongoing century in central Europe [20]. It still represents an ecologically and economically important tree species in central Europe (e.g., [30,66,67]). However, only a few papers have been devoted to the mechanical stability of Norway spruce stands along an altitudinal (ecological) gradient with various site conditions (i.e., acid vs. nutrient-rich sites). Thus, the presented paper aims to find the most suitable sites for spruce growing related to its mechanical stability (HDR; Crown ratio, CR), based on the assessment of an extensive dataset from the National Forest Inventory (NFI) performed in the Czech Republic. The results will also help to formulate silviculture recommendations that could contribute to better mechanical stability parameters in existing spruce stands.

We hypothesise:

H1. 

Height-to-diameter ratio (HDR) values culminate in 40-year-old and younger pure Norway spruce stands.

H2. 

HDR increases and CR decreases with lower altitudinal vegetation zones.

H3. 

A higher HDR will be observed at nutrient-rich sites, promoting intensive height growth.

2. Materials and Methods

2.1. Data Source

All the biometric data were obtained within the first NFI in the Czech Republic and performed by the Forest Management Institute (government organisation) in 2001–2004 [68]. Circular inventory plots with an area of 500 m² were selected randomly within inventory squares with a spacing of 2 km × 2 km (two plots per square). For this study, 14,220 inventory plots were established and measured across broad forest stand conditions countrywide—the field measurements involved all trees with a DBH above 12 cm occurring in the inventory plots.

We selected 3969 inventory plots with dominating Norway spruce and a minimum of 15 spruces per plot, using information about stand age, Forest Vegetation Zones (FVZ), and Ecological Series (ES) at the plot level. In the Czech typological system, FVZs are defined by the relationship between the climate (elevation is a proxy for climate conditions) and the tree layer composition of forest stands along an altitudinal gradient [69]. For modelling purposes, six FVZs (from second to seventh) with a sufficient data pool from the NFI were selected (

Table 1. Basic characteristics of forest vegetation zones (FVZs) in the Czech Republic.

Forest Vegetation Zone (FVZ)	Mean Annual Air Temperature (°C)	Mean Sum of Annual Precipitation (mm)	Growing Season Length (Days)
2nd Beech—Oak zone	7.5–8.0	600–650	160–165
3rd Oak—Beech zone	6.5–7.5	650–700	150–160
4th—Beech zone	6.0–6.5	700–800	140–150
5th Fir—Beech zone	5.5–6.0	800–900	130–140
6th Spruce—Beech zone	4.5–5.5	900–1050	115–130
7th Beech—Spruce zone	4.0–4.5	1050–1200	100–115

).

In the typological system of the Czech Republic, forest site units are organised by edaphic categories, which are further divided into ecological series (ES) based on the similarity of the soil conditions [69]. Within the statistical software environment R 4.1.0 [70], the two most frequent ESs were analysed in the study:

(1)

nutrient-rich—this represents mesotrophic soils without a significant soil water influence.

(2)

acidic—this represents oligotrophic soils with no significant soil water effect.

2.2. Aggregation of Tree-Level Data

We used the DBH, total height, and crown base height data at the individual tree level. The diameter–height relationship was modelled individually per plot according to the Levakovic function [71]:

$$h=a\times {\left(\frac{d}{d+1}\right)}^b+1.3$$

(1)

where h is the tree height (m), d is the DBH (cm), and a and b are regression coefficients.

The dominant stem diameter was calculated as the quadratic mean diameter for the 100 largest trees per hectare, i.e., the five thickest trees per plot (500 m²).

$$d_{dom}=\sqrt{\frac{\sum d_i^2}{5}}$$

(2)

where d_dom is the dominant stem diameter, d_i is the individual tree diameter of the dominant trees, and 5 is the number of dominant trees per plot. Subsequently, the dominant stem height (h_dom) was calculated using the Levakovic function’s parameters [71].

The height-to-diameter ratio (HDR) of the dominant stem at the plot level was calculated as follows:

$$HDR=\frac{h_{dom}}{d_{dom}}\times 100$$

(3)

where HDR is the height-to-diameter ratio, h_dom is the dominant tree height (m), and d_dom is the dominant stem diameter (cm).

The mean crown ratio (CR) at the plot level was calculated as follows:

$$CR=\frac{\sum _{i=1}^n\frac{l_i}{h_i}}{n}$$

(4)

where CR is the crown ratio, l_i is the crown length of a tree (m), h_i is a tree’s total height (m), and n is the number of trees per plot.

2.3. HDR-Age Relationship Modelling

A semiparametric generalised additive model (GAM) from the R mgcv library [72] was used to model the HDR–age relationships. A GAM can be viewed as a GLM (generalised linear model), in which the linear predictor partly and linearly depends on some unknown smooth function f. GAMs are estimated by a penalized version of the method used to fit GLMs. In general, a GAM has the form [73]:

$$E\left(Y\right)=g^{-1}({\beta }_0+{\sum }_{j=1}^J{f}_j(x_j)$$

(5)

where E(Y) is the mean value of the response Y (with an appropriately defined distribution and link function g), f_j is a smooth function of the covariate x_j, β₀ is the intercept term, and g⁻¹ is the inverse link function. Each smoother f_j is represented by a sum of K basis functions (b_j,k) multiplied by coefficients (β_j,k):

$$f_j\left(x_j\right)={\sum }_{k=1}^K{\beta }_{j,k}{b}_{j,k}\left(x_j\right)$$

(6)

where K determines the maximum complexity of each smoother and the large basis size K could lead to overfitting, which is compensated for by a smoothing penalty [73]. In the mgcv library, thin plate regression spline (TPRS) [72,74] uses a penalty matrix based on the integral of the squared derivates of the basis functions [72,73,75] in the form:

$$\lambda \int {\left[f_j^m\left(x\right)\right]}^{2}dx$$

(7)

where f_j^m is an m-th order derivation (usually second-order) and λ is a penalty constant (also called a smoothing parameter). The integral (7) is a roughness penalty [76]. The integral (7) can be rewritten into linear algebra form:

$$\lambda {\beta }^{T}S\beta$$

(8)

where S is the penalty matrix multiplied by the parameter vector β (T denotes matrix transposition), and the penalty term (8) is then subtracted from the model log-likelihood [73]. The restricted maximum likelihood method (REML) was chosen for searching the optimal λ value, because the default method in the mgcv library (GCV—generalized cross-validation) tends to overfit smaller or noisy datasets. Overfitting tends to look like spline fits that are too wiggly [73].

The GAM allows us to model the nonlinear relationship between the HDR and stand age (covariate) for each level of factor variable (FVZ). As a smoother, we specifically used thin-plate regression spline [72,74] in the mgcv R library, which is appropriate for models with many independent variables, both continuous and categorical. The dependent variable (HDR) was modelled as it comes from a gamma distribution (with a logarithmic link transform function), because the model with a normal distribution showed skewness and heteroscedasticity in its residuals. The model can be informally written as follows:

$$\log \left({HDR}_{ij}\right)=\alpha +f_j\left({Age}_i\right)+{FVZ}_j+f\left({Age}_i\right):{FVZ}_j$$

(9)

$${HDR}_{ij}~Gamma\left({\mu }_{ij},\phi \right)$$

(10)

where α is the intercept term, f is the j-th smoothing term for the i-th age, Age is a smoothed covariate, FVZ is the factor variable with the j-th level of forest vegetation zones, and Age: FVZ is the interaction. Including an interaction term in the equation guarantees the computation of a smooth term for each factor variable (FVZ) level with its penalties. The term log(HDR) indicates that the dependent variable was calculated on a logarithmic scale. The notation: ${HDR}_{ij}~Gamma\left({\mu }_{ij},\phi \right)$ indicates that the dependent variable came from a gamma distribution with a mean of µ for the i-th age and j-th FVZ, and φ is the dispersion parameter. The gamma distribution is a continuous, two-parameter distribution belonging to the exponential family and is used for positively skewed and strictly positive real data [77].

In R code, the model is specified as:

model_9_10 <- gam (HDR ~ s(Age, by = FVZ, bs = ‘tp’) + FVZ, method = ‘REML’, family = ‘Gamma’(link = log)).

2.4. Modelling HDR–Age Relationship by Ecological Series (ES)

Inventory plots in two ES series were selected for modelling the HDR–age relationships in terms of ES: nutrient-rich (1492 plots) and acidic (1598 plots). As mentioned above, we added the factor ES to the HDR–age relationship model.

$$\log \left({HDR}_{ijk}\right)=\alpha +f_j\left({Age}_i\right)+{FVZ}_j+{ES}_k+f\left({Age}_i\right):{FVZ}_j$$

(11)

$${HDR}_{ijk}~Gamma\left({\mu }_{ijk},\phi \right)$$

(12)

In R code, the model is specified as:

model_11_12 <- gam (HDR ~ s(Age, by = FVZ, bs = ‘tp’) + FVZ + ES, method = ‘REML’, family = ‘Gamma’(link = log)).

2.5. CR–Age Relationship Modelling

The CR–age relationship model was built analogically, the same as the HDR–age model, but the distribution function was switched from Gamma to Beta with a logit-link function. The beta distribution is a two-parameter distribution beyond an exponential family distribution and is appropriate for response variables ranging between 0 and 1 [77]. The logistic function is a canonical link function for beta distribution models. The model can be written as follows:

$$\log \left(\frac{{\pi }_{ij}}{{1-\pi }_{ij}}\right)=\alpha +f_j\left({Age}_i\right)+{FVZ}_j+f\left({Age}_i\right):{FVZ}_j$$

(13)

$${CR}_{ij}~Beta\left({\mu }_{ij},\phi \right)$$

(14)

where log(π_ij/(1 − π_ij) is the logistic transformation of the response variable CR, which is assumed to come from a beta distribution with the parameters µ (mean) and φ (dispersion parameter).

In R code, the model is specified as:

model_13_14 <- gam (CR ~ s(Age, by = FVZ, bs = ‘tp’) + FVZ, method = ‘REML’, family = ‘betar’(link = logit)).

3. Results

3.1. Development of Dominant Trees’ HDR by Forest Vegetation Zones (FVZ)

The HDR development of the dominant spruce trees differed in terms of the FVZ (Figure 1), although the trend of an initial increase and a later HDR decrease was noted for almost all the analysed FVZs. Generally, the highest FVZ (i.e., highest elevation) corresponded to the lowest HDR of the dominant trees. The lowest HDR values (below 65 throughout the observed period, i.e., for 20–100-year-old stands) were found for “montane” FVZ 7. In the downwardly following FVZ 6, the HDR continually increased from 62 to 70 in the analysed age period.

An initial steeper HDR value increase was found in FVZs 2–5. The HDR in FVZ 5 culminated before the age of 40 years at the level below 75 and slowly decreased in older stands. A similar pattern was shown in the results from FVZs 4 and 3. However, the HDR culminated at approximately 45 years and at higher values of 77 and 79 in FVZs 4 and 3, respectively. The highest HDR values, indicating the lowest mechanical stability, were found in the lowest-located FVZ 2. In this “lowland” zone, the HDR of the spruce dominant trees increased to 83 in 55-year-old stands. Subsequently, it decreased to level 77.

The above-mentioned order of the HDR values of the dominant trees in the studied FVZs is also evident from Figure 2. Empirical cumulative density curves (ECDF) showed a similar pattern, while the most favourable HDR distribution was evident in FVZ 7, and the most adverse was observed in FVZ 2.

3.2. Development of Dominant Trees’ HDR in Acidic and Nutrient-Rich Ecological Series (ES)

The differences between the HDR values of the dominant trees at acidic and nutrient-rich sites were found mainly for the first half of the rotation period, i.e., up to the age of 50–60 years (Figure 3). The HDR culminated in lower values at the nutrient-rich sites compared to the acidic ones for all the analysed FVZs. Therefore, more propitious conditions at the nutrient-rich sites significantly contributed to higher radial and lower height increments, leading to lesser HDR values.

3.3. Development of Mean Crown Ratio (CR) by Forest Vegetation Zones (FVZs)

The results revealed that the CR naturally decreased with an increasing age in all the analysed FVZs (Figure 4). The highest CR values (0.7–0.8) were detected in 20-year-old stands and the lowest (below 0.5) in 100-year-old ones. This age-related decrease trend was more gradual in the lowest (2–5) compared to the highest (6–7) FVZs. Thus, the most favourable CR values (i.e., the highest) were detected in “montane” (FVZ 7) throughout the studied stand age range.

The empirical cumulative density curves (ECDF) showed a similar pattern (Figure 5) as well, while the CR distribution in FVZ 7 was the most favourable (approximately 75% of all the analysed inventory plots showed a mean CR exceeding 0.5). The percentage decreased with lower FVZs (only around 30% of all the inventory plots showed a CR of > 0.5 in FVZ 2).

4. Discussion

NFI data are routinely used to monitor forest stands’ productivity (e.g., [33,50,78]) or evaluate ecosystem services and biodiversity (e.g., [79,80,81]). The presented study used NFI data from the Czech Republic to assess the mechanical stability of Norway spruce. Recently, a nonlinear mixed-effects HDR model was also developed on the basis of these Czech NFI data by Sharma et al. [55] for several tree species. Additionally, the predictive performance of this model was evaluated using data from research sample plots. The dominant height and diameter, relative spacing index, and DBH-to-quadratic mean DBH ratio were identified as the most important predictors of HDR variations. Our approach to the NFI data processing differed from the above, because we used the HDR–age and CR–age relationships for the modelling.

Our results confirmed that the HDR changed naturally with age, culminating in the first half of the standard rotation period for spruce in lower and middle elevations, supporting hypothesis H1 at these altitudes. This was also consistent with the findings of Cremer et al. [57], where increases in the HDR values were most pronounced up to dominant heights of 20 m or 25 m, and for greater heights, the HDR values levelled off and even declined with further growth.

Dominant trees are naturally more stable, i.e., the HDR values are the lowest for dominant and highest for suppressed trees [82], but especially at lower elevations (FVZ 2), where we found the HDR culmination of the dominant spruce trees to be over the value of 80, where the risk of abiotic damage (firstly from snow and later from wind) is increased [63].

In this study, higher altitudes (represented by higher FVZs) generally represented lower HDR values throughout all the analysed age periods (see hypothesis H2). Although spruce stands in the Czech Republic do not always correspond to their original populations (due to the historical use of spruce outside of its natural distribution area), the results confirmed that these stands are more stable in higher locations [83,84], primarily due to their lower height growth, which is the result of adaptation to the given environment [62,85,86,87,88], and the higher radiation use efficiency of mountainous spruce stands [89], which probably allocate carbon to the stem (i.e., to a radial increment resulting in a lower HDR value and higher stand stability). In contrast, silvicultural measures promoting crown morphology and the growth of young spruce stands are less effective in the mountains than lower-located sites [90].

A lower elevation also generally represents more favourable site conditions (mainly through deeper soils with higher available nutrients). These conditions lead to larger diameters, but especially to height growth (expressed, for instance, by site index), i.e., a higher site index corresponds with a lower mechanical stability—higher HDR values [91,92]. On the contrary, we found in our analyses that the spruce trees attained more favourable HDRs at nutrient-rich sites compared to acid ones at the same age (hypothesis H3 was rebutted). The differences were not high, but they were confirmed for all the analysed FVZs in the first half of the rotation period. More favourable growth conditions promoted that height growth did not substantially overwhelm radial growth, resulting in a not-too-intense increase in the HDR value in youth (thickets and small pole stands). The increased growth of spruce and beech recorded in Europe [93] could also play a specific role, especially at nutrient-rich sites. On the other hand, Dobbertin [94] in Switzerland and Bouchard et al. [95] in Canada observed that windthrow occurred more frequently at fertile sites than poor sites. Therefore, our finding leads to the need for further research, focusing on the characteristics of stands in different conditions (site indexes), but under identical thinning management.

The analysed NFI plots are managed “in a standard way”, and the detailed thinning regimes for individual locations are unknown. Stand density and management applied through thinning treatments are essential from a stability point of view, expressed in terms of HDR or CR. The competition was recorded in ponderosa pine [87], loblolly pine [96], or Engelmann spruce [88] dense stands as a significant factor affecting the stem and crown parameters in connection with the risk of damage from snow or wind. Thus, windstorms often suffer well-stocked mature Norway-spruce-dominated stands [97] or dense coniferous stands [94,95].

The HDR or CR values are also influenced by a tree’s distance from the stand edge [98]. Trees growing at the stand’s edge with longer live crowns are naturally more stable (a lower HDR and higher CR) than trees inside the stand.

These results should be, together with the knowledge of the thinning effects on stability [51,99], the basis of practical recommendations for forest practices to reduce the risk of forest damage. Of course, besides thinning applications as a mitigation measure, it is necessary to focus on forest regeneration, which is a principal measure enabling species composition change [100,101] in order to reduce the proportion of spruce (especially in lowlands and middle elevations) and increase the share of mixtures, characterised by a higher annual volume increment [102] and having a more positive effect on productivity [103] and increased temporal stability than monocultures [104].

5. Conclusions

Norway spruce stands are damaged by abiotic and biotic harmful factors in many European localities, often due to their historical extension to sites outside of the ecological optimum for this tree species. In addition to a gradual change in species composition towards mixtures, it is essential to manage spruce stands to minimise the risk of their large-scale disintegration and allow enough time for the abovementioned conversion. In forestry management, it is advisable to adapt procedures to the different site conditions where spruce stands grow. The presented study confirmed the facts based on the NFI data from the Czech Republic. Spruce stands in lower elevations show a lower stability (expressed in terms of HDR) than those in the mountains. The HDR culminated in lower and middle elevations in the first half of the rotation period, representing the time with potentially the highest silviculture measures (thinning) effectiveness. The study did not confirm predicted higher HDR values at nutrient-rich sites than at acid ones, especially for the first half of the rotation period. Therefore, this topic should be the subject of further research, mainly concerning the same thinning regime at different site conditions.