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

: 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

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 3 ) 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 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.

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 2 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). 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.

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]: 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 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: 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: 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.

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]: 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 , β 0 is the intercept term, and g −1 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 ): 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: 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: 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: 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 µ ij , ϕ 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].

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.

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: 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).

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.

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. 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.

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 nutrientrich 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.

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 nutrientrich 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.

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.  (FVZs 2-7). Shadow areas depict 95% confidence intervals of means. Red and blue lines represent acid and nutrient-rich sites, respectively.

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)(3)(4)(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).

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

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 lowerlocated 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 nottoo-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].

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.

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

Acknowledgments:
We sincerely thank the Forestry Management Institute (FMI) for sharing raw data from the National Forest Inventory (NFI) analysed in the presented study. Authors are indebted to anonymous reviewers for their valuable comments, significantly improving the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the study's design, data collection, analysis and interpretation, manuscript writing, or publication.