Temporal Dynamics of Root Reinforcement in European Spruce Forests

: The quantiﬁcation of post-disturbance root reinforcement (RR) recovery dynamics is of paramount importance for the optimisation of forest ecosystem services and natural hazards risk management in mountain regions. In this work we analyse the long-term root reinforcement dynamic of spruce forests combining data of the Swiss National Forest Inventory with data on root distribution and root mechanical properties. The results show that root reinforcement recovery depends primarily on stand altitude and slope inclination. The maximum root reinforcement recovery rate is reached at circa 100 years. RR increases continuously with different rates for stand ages over 200 years. These results shows that RR in spruce stands varies considerably depending on the local conditions and that its recovery after disturbances requires decades. The new method applied in this study allowed for the ﬁrst time to quantify the long term dynamics of RR in spruce stands supporting new quantitative approaches for the analysis of shallow landslides disposition in different disturbance regimes of forests.


Introduction
In mountain regions, many forests protect human settlements and activities against gravitational natural hazards [1]. At the same time, forest ecosystems worldwide are increasingly perturbed by climate-driven natural disturbances [2], which may determine a temporary reduction in the provision of forest ecosystem services, including protection from natural hazards [3]. One of the priorities of forest management is therefore to avoid or shorten temporary gaps in forest's protective effectiveness following forest fires or stand-replacing silvicultural interventions (e.g., coppicing, clearcutting) [4,5]. Given the limited resources for forest management, prioritizing silvicultural measures based on the long-term dynamics of forest ecosystems and their protective effect may be of paramount importance in the context of an integrated risk management of natural hazards [1].
The relationship between post-disturbance root reinforcement recovery dynamics in a forest stand and the hazards of shallow landslides is a prime and documented example in this context. The long-term dynamics of root reinforcement in mountain forest management has important implications for the decision makers, not only for post-disturbance management measures but also for optimising forest ecosystem services of forests in general. Post-event time windows of increased shallow landslide disposition may range from a few years [3,6] to several decades [5,7] depending on the disturbance severity, tree species concerned, and site characteristics. While the recovery of vegetation cover in the short term has in general an important effect in reducing the runoff in the contribution area of a potential landslide and in the related effect on building of pore water pressure [8], in the long term, root reinforcement becomes the dominant factor for slope stability in areas of potential shallow landslides [9,10].
A large number of studies quantify the short-term dynamics of root reinforcement after a disturbance in forest stands [3,[11][12][13][14] considering both changes in the spatiotemporal distribution of roots [13,15] and the change in their mechanical properties due to decay [16,17]. Vergani et al. [18] showed that after a stand-replacing fire in scots pine (Pinus Sylvestris) stands, the reduction in root reinforcement efficiency was mainly due to the loss of the mechanical properties (such as tensile strength and stiffness) of coarse roots, implying the total loss of the reinforcement within 5-10 years. Similar rates of decay have been reported following clearcut coppicing [19]. Studies dealing with dynamics of root reinforcement recovery are very rare [4]. Ziemer [20] reported that a 50% recovery of root reinforcement was reached within 15-25 years after clearcutting of lodgepole pine (Pinus contorta). Gehring et al. [14] found that the gap in protection against shallow landslides due to the loss of root reinforcement in beech (Fagus sylvatica) stands following wildfires may extend over 35 years in case of medium to high severity fires. Liu et al. [10] showed that the recovery of root reinforcement in the first four post-disturbance years in closed forest stands at high altitudes (i.e., 1800 m a.s.l.) was slow or almost absent, suggesting that spruce stands at this altitude need particular attention for the study of root reinforcement dynamics in the long term.
In this paper, we analyse the long-term dynamics of root reinforcement in spruce forests (Picea abies). This species is known to have a relatively low specific root reinforcement compared to other tree species such as beech or chestnut (Castanea sativa) [14,19,21]. However, spruce is the most common tree species of the reported Alpine protective forest area [22], being dominant in Swiss and Austrian protective forests [23,24]. Moreover, spruce is considered one of the most sensitive species to the ongoing climate change and is expected to be increasingly vulnerable to disturbances in future decades [25], particularly due to drought stress and bark beetle outbreaks [26,27]. In addition, in its upper natural mountain range, spruce forests usually display slow stand dynamics and regeneration rates [28,29]. In the context of an integrated forest management in the Alpine region, analytical models for the quantification of post-disturbance root reinforcement dynamic of this species are particularly needed in assessing future protective capacity, as well as for planning silvicultural measures.
The aim of this paper is thus to characterise the long-term root reinforcement dynamics of spruce forests along an altitude gradient by assessing the post-disturbance recovery of lateral root reinforcement at the stand scale. To this purpose, we combined data of the Swiss National Forest Inventory (Swiss NFI) and other European forest inventories with published data on root distribution and root mechanical properties of spruce to calibrate the Root Bundle Model [30] and to upscale the related root reinforcement effect at stand scale [14,21].

Workflow and Methodological Approaches
The flowchart in Figure 1 shows the overall methodological steps needed for the calculation of the long-term root reinforcement dynamics at the stand scale. In the first step, we test different single tree growth models to predict the mean annual diameter at breast height (DBH) increment of spruce trees using as reference Swiss NFI trees. In a second step, we calibrate a stand-scale model for the estimation of mean stand age based on observed forest structure and mean tree size (stem diameter). In the third step, stand characteristics and root reinforcement data are combined to calculate the lateral root reinforcement at stand scale. The calculated root reinforcement at the stand scale is then paired to the corresponding stand ages. Finally, in a fifth step, the temporal dynamics of root reinforcement are calculated using an empirical model.

Data Sources
Long-term spruce stands data are retrieved from the Swiss NFI and other European forest inventories. For about 40 years, the Swiss NFI has collected data at tree as well as at stand level using systematic and well-described methods [31,32]. The first Swiss NFI (NFI1) was conducted between 1983 and 1985, followed by three additional surveys in 1993-1995 (NFI2), 2004-2006 (NFI3), and 2009-2017 (NFI4). Since 2009, data have been collected over a nine-year period. Currently, the fifth inventory is underway from 2018 to 2026 (NFI5). The terrestrial data of NFI1 were collected on some 11,000 plots located on a systematic 1 km × 1 km grid covering the entire country [33]. Since NFI2 the number of sample plots is reduced, and field measurements are conducted on a 1.4 km × 1.4 km subgrid of the NFI1 grid resulting on about 6500 forest plots in NFI4. On each plot, two concentric circles (i.e., 200 m 2 and 500 m 2 in size, respectively) are defined to measure all trees and shrubs with a DBH between 12 cm and 36 cm in the inner circle, respectively, and all trees with a DBH larger than 36 cm in the outer one [33].
Site characteristics such as altitude, aspect, and slope as well as the stand density index (SDI) were also retrieved from Swiss NFI data to be used as explanatory variables for modeling the DBH growth of individual trees. SDI is defined by Reineke et al. [34] using the allometric coefficient described in Schütz and Zingg [35]. Tree density (number of stems per hectare) and mean size (quadratic mean DBH) were also collected from the Swiss NFI for the upscaling of the root reinforcement at stand level. Furthermore, 2640 plots from Italy, Germany, Czech Republic, and France were included in order to compare the results of the Swiss stands with values obtained for the rest of the continent [36]. This dataset is named "European" in this manuscript.
Data on spruce root distribution and mechanical properties (maximum tensile strength and apparent secant stiffness) were retrieved from published studies on root reinforcement of spruce stands [13,37,38]. The spatial distribution of roots in spruce stands was calculated following Schwarz et al. [21] using data from 95 soil profiles in six different locations in Switzerland [37,38]. The mechanical properties of roots were quantified based on 93 field pull-out tests described in Vergani et al. [37] and Moos et al. [38].

Stem Diameter Growth Modelling
To quantify the average DBH growth of spruce trees under the existing Swiss environmental conditions, we built a simplified model using diameter increment data of more than 50,600 individuals belonging to 8631 spruce-dominated plots and measured in at least two consecutive Swiss NFI inventories (i.e., NFI1-NFI2, NFI2-NFI3, NFI3-NFI4). We defined spruce-dominated stands as Swiss NFI-plots with at least 80% of the basal area represented by spruce trees.
In our modeling approach, we first modelled tree growth as a function of individual tree DBH and then refined the model using stand and site parameters such as SDI, altitude, aspect, and slope, implementing either a nonlinear [29] or a linear modeling approach.
The mean annual tree diameter increment (in cm) at the stand scale can be considered as the approximation of a continuous growth rate, i.e., where τ represents the beginning time of the measurements of the DBH. If we denote by t the age (years) of a spruce tree, we have that: where c 1.3 is the mean time needed by seedlings to reach 1.3 m height. Data for the estimation of c 1.3 were sourced from an unpublished dataset from the Swiss Alps [39], containing totally 231 spruce. The spruce at <1200 m.a.s.l. reaches 1.3 m height in 9 years, while this time increases to 12 and 17 years respectively for high montane (1200-1600 m.a.s.l.) and subalpine forests (1600-2100 m).
We present hereafter two suitable models to model DBH growth. The nonlinear model (NLM) is based on the Gompertz function: where β 1 and β 2 are model coefficients. We consider a small error term ε error ≈ 0 and take as initial measurement: i.e., the measured initial DBH of the spruce at the age t = τ + c 1.3 years is small but different from 0. NLM will not allow to use a more simplified initial condition and then reduce all to ε DBH = 0. Imposing ε DBH = 0 the differential equation will lead to the trivial solution where r indicates the growth rate and K the carrying capacity (or upper asymptote for growth). Relations among coefficients of Equations (3) and (5) are as follows: According to Schelhaas et al. [29] the parameters β 1 and β 2 are dependent upon a number of environmental and site variables: where θ i,2 is a weighted array of model coefficients. This approach allows the weighting of individual variables (X i ) to be discussed in a multivariable model. The analytical solution of this model is given by: where C is the constant related to the initial condition DBH(0) = ε BHD : The linear model (LM) proposed as a simplification of NLM: Since the measurement of DBH begins at τ = 0 (i.e., at age c 1.3 years, when the spruce reaches the height of 1.3 m) we have again: One of the advantages of LM is the possibility to use a suitable simplified approximated initial condition without losing the sense of the growth process: we can here impose DBH(τ = 0) = 0. This initial condition makes sense even if inserted in Equation (12): The "initial" growth rate at τ = 0 of the spruce tree is given by the strictly positive coefficient β 1 .
The general solution of the classical ordinary linear differential Equation (12) is: where From (15) and (16), we obtain where parameters β 1 and β 2 were again determined using Formulas (8) and (9) for NLM.

Estimation of Mean Stand Age
Considering Equation (14) and its solution (17) to the mean annual DBH growth rate of the stand, it is possible to determine a formula for τ. In particular, it is possible to represent the mean age of a spruce stand as a function of its quadratic mean diameter (QMD): The age of the stand, t (years), is therefore:

Upscaling of Root Reinforcement
The Root Bundle Model-Weibull (RBMw) was used as a basis for the calculation of the root reinforcement at tree and upscaled to stand level following Schwarz et al. [21], as described in Vergani et al. [13] and most recently in Gehring et al. [14]. The upscaling consists of two parts. On one hand, (1), roots of several trees are counted and measured for their thickness in soil trenches at different distance from tree stems. On the other hand, (2), the mechanical properties, tensile strength, and apparent secant stiffness of the roots are needed to determine the maximum pulling force of the roots [30,40]. The root mechanical properties for different root diameter (up to 40 mm) are determined by field pull-out tests [13,19,38,41].
From these parameters (i.e., number of roots and root mechanical properties) a model of the maximum tensile force at the level of a single tree root system can be generated.
The maximum root reinforcement RRmax is calculated with the following formula implemented in the SlideforNET model [14,15] : where Γ stands for the gamma distribution, and d To quantify root reinforcement, stand data from the Swiss NFI surveys were used. The lateral root reinforcement was calculated from the 7191 (filtered from the original 8631, considering only stand with RRmax > 0) stands using the corresponding input parameters and SlideforNET.

Model for the Temporal Dynamics of Root Reinforcement after Disturbances
The data on calculated stand age (Sections 2.3 and 2.4) and the data on upscaled stand root reinforcement (Section 2.5) were used to model the recovery of root reinforcement as function of time since disturbance. The model is based on the cumulative distribution function of the two parameter Weibull distribution: where s_1, q1, q2, k, and λ are coefficients fitted after Schelhaas et al. [29].
The coefficients of the model were fitted using the nls() function of the R software [42]. The upper and lower 95% prediction interval were calculated using standardized residuals [43], in order to consider heteroscedasticity. Age 0 is considered as the year of first establishment of spruce regeneration after a disturbance, which may depend on the type of successional dynamics that takes place at any specific site.

Modeling of Growth Rate
Based on the analysis of the Swiss NFI data, the yearly diameter growth of single spruce trees ranges between 0.1 and 1 cm/year. The overall median value is 0.305 cm/year, whereas the mean value is 0.377 cm/year. Most of the data refer to a range of DBH of single trees between 25 and 50 cm. The results of the model fitting show that the performance of the two models (NLM and LM) are similar as well as the significance of the model variables (Tables 1 and 2). For both models the SDI is the most important parameter, followed by slope and altitude. Aspect has no statistically significant influence on growth rate. Overall, assuming a median value for each predictor, LM predicts an almost constant growth rate of about 0.4 cm/year (Figure 2). The variability of model results decreases with increasing DBH.
In view of these results, we chose to drop the "aspect" predictor (see Table 3 for coefficient) and to use LM for further analyses.

Calculation of Root Reinforcement
The information on mean DBH and stem density for each stand were used to calculate the minimum values of lateral root reinforcement at the stand scale for the Swiss NFI and European datasets (Figure 3). The two ranges of root reinforcement overlap up to a DBH class of about 50 cm, while for larger DBH the returned root reinforcement was higher in the European dataset. In general, lateral root reinforcement of spruce stands in Europe can be expected to be linearly proportional to the mean DBH of the stand and to have values up to ca. 30 kN/m.

Modelling Root Reinforcement Dynamic
The calculated lateral root reinforcement (Equation (19)) increases continuously with increasing stand age up to median values of about 10 kN/m at 200 years ( Figure 4). 50% of the root reinforcement values (i.e., from the 25% to the 75% quantiles) for stand ages higher than 150 years range between 6 and 13 kN/m. All the parameters of the RR˜age model have a p-value lower than 0.001 (Table 4), indicating a statistically significant influence on the results. The negative values of the coefficient q1 and q2 indicate that lateral root reinforcement decreases with increasing altitude and slope. The scaling factor λ of the Weibull distribution has a value of 150, and the shape parameter k is 2.017. The variability of the predicted values increases with increasing root reinforcement, whereby up to values of 12 kN/m, the range of 50% predicted values is in the range of approximately 2 kN/m ( Figure 5). A sensitivity analysis in accord to Schelhaas et al. [29] of the RR˜age model shows that altitude influenced the results up to approximately 30% of their median value, whereas the slope influences up to approximately 15% ( Figure 6). The inflection point of the curve is at 106 years (corresponding to the mode value of the Weibull cumulative distribution).

Comparison of the Two Growth Rate Models
The two models applied to describe the DBH increment in this paper are related to the classical topic in the mathematics of differential equation for population growth models. The choice of a model to fit a growth process can be done stepwise considering which type of complexity can be afforded and at the same time, which prediction accuracy is requested for the model applications. The use of a linear model is considered in general to give a good approximation of an initial growing phase or for a very slow growth process. Vice versa, "logistic" or "Gompertz" models are expected to fit better for more complex growing behaviours [29,44]. The very large dataset considered in the analysis supports, within the aim of this study, the linear growth model approximation to perform best, allowing us to reduce the model complexity by however achieving a suited prediction accuracy. In our specific case, incertitudes due to the 1 cm measurements accuracy of the DBH are compensated by the high number of trees measured and the years between two inventories. The linear behaviour of the growth rate data in this work is influenced by two additional characteristics: firstly, measurements are taken for trees with DBH > 12 cm, meaning that the usually higher growth rate of the first years is not considered, and secondly, the growth rate is calculated as the mean considering the DBH increment of single trees over 4-6 years, leading to a smoothing of the annual growth rate, especially in the first years.
In selecting the model variables we paid attention to inserting only easily measurable parameters in order to make the proposed approach useful for practical applications such as the assessment of forest recovery after disturbances. As a result, we ended up with a reduced number of variables in comparison to Schelhaas et al. [29] and we used for instance the easy retrievable altitude instead of the correlated temperature data. However, our approach makes it more difficult to project current dynamics into the climate change scenarios. Nevertheless, present studies did not reveal significant changes in the growth rate of spruce along altitude gradients in the European Alps, indicating no climate change effects so far [45].
In the model, the SDI for the estimation of the stand age is considered as a static variable, although it may change over time in reality. In the scope of this work, however, the SDI calculated from the Swiss NFI data is considered representative for predicting the root reinforcement dynamics after a disturbance. Future studies on this topic could implement more variables such as soil properties and type of forest management to improve the performance of the growth rate model. The competition between trees and climate extremes considerably influences the growth rate and eventually also the survival chance of single individuals [46]. However, the effects of self-thinning mechanisms and/or silvicultural interventions have in general small and limited-in-time impacts on the overall growth rate of spruce stands [47,48]. Considering the results by Uhl [49], this behaviour seems limited to less productive sites, whereas in favourable conditions for spruce, single trees show significant growth increment after thinning.
Based on the NFI data used for calibration, the proposed model assumes site specific optimisation of the overall stand growth rate at the stand scale, which results from the balancing competition for resources and the related population densities [34]. Pretzsch et al. [50] discussed the reaction pattern of forest stands or cohorts as more than just the growth response of individual trees depending on their competition for resources. The results of our study for spruce show a complete symmetry of growth rate versus single tree size (DBH). This result could be interpreted as a stronger competition of trees for limited below-ground resources (independently of the tree size) [51] and less than competition for light (which is expected to lead to a stronger correlation between growth rate and single tree dimension). Pretzsch et al. [50] state that "simulations for archetypical stands yield a transition from size-asymmetric to size-symmetric competition along the gradient from fertile to poor sites". In combination with our results, this statement leads to the conclusion that spruce in the NFI dataset tend to compete for soil resources, more than for light. A detailed analysis of dendrochronological data could be used in future works to further investigate these aspects.

Correlation between Mean DBH and Stand Age
The general trend of a constant growth rate is confirmed by the linear correlation between age and mean DBH of dominant trees in a stand, as reported in yield tables for spruce [52]. Figure 7 compares model results of this work with the different growing classes ("Bonitäten", growing potential of the tree species on a given site) of the Swiss yield tables. Similar trends are shown in Stöcker [53], Ammann [54], Schelhaas et al. [55], Simmler [56], and Kindermann et al. [57]. Castagneri et al. [58] obtained DBH˜age curves of single spruce trees in Norway from dendrochronological analysis. Their results show that fast growing young trees generally tend to die, whereas older trees have a more regular linear growth. This supports the conclusion that mean growing rate of forest stand is dominated in the long term by trees with a constant linear growth. The results of this study show that site altitude and slope can explain most of the growing capacity of the concerned spruce stands. Figure 7. Comparison of the increments of the basal mean stem (dg) of spruce according to the Swiss yield table [52] with integration of the BHD increments resulting from the Swiss NFI data. The red dots stand for the subalpine altitudinal stage, blue stand for high montane, light green upper/lower montane, and dark green stands for the submontane/colline altitudinal stage.

Upscaling of Root Reinforcement
We use the minimum lateral root reinforcement as a proxy variable for the characterisation of the mechanical effect of root on slope stability at the stand scale. In reality, lateral root reinforcement is highly variable within a stand, but potential tension cracks of a landslide can be assumed to follow the weakest paths within a stand. To improve and extend this calculation to the stabilization effect of roots, basal root reinforcement should be considered too. In contrast to the lateral reinforcement which acts on the vertical edge of a potential landslide, the basal root reinforcement interacts along the shear plane of a potential landslide (parallel to slope) [15]. The basal root reinforcement is thus strongly influenced by characteristics of the site [59,60] and the depth of the potential shear plane of a landslide [14,61]. In particular, in the case of very shallow landslides (with depth < 0.5 m) low values of lateral root reinforcement may result in high values of basal root reinforcement. High values of lateral root reinforcement may on the contrary result in a null basal root reinforcement if the landslide is deep enough. Gehring et al. [14] present an approach that shows the effect of different basal root reinforcement on slope stability.
In consideration of the variables used to calculate the lateral root reinforcement (mean stem DBH and mean distance between stems), it is important to mention that the dimension of the Swiss NFI plots can be considered big enough to be representative of a forest stand and small enough to avoid extreme heterogeneous pattern of vegetation cover and structures usually present in larger areas of forest.

Root Reinforcement Dynamics
The post-disturbance window of susceptibility for shallow landslides mainly depends on root decay of the dying trees, recovery rate of the new stand in terms of roots and canopy, and the frequency and magnitude of triggering rainfall events. Considering that total root decay is usually achieved in a temporal range of 5 to 15 years, the window of susceptibility can in general be expected to persist over 50 to 100 years. In reality, this temporal window tends to be shorter in stable areas (where landslides are triggered by extreme events only) and longer where high values of root reinforcement are needed to stabilise the slope. A detailed and spatial-explicit quantification of this effect is fundamental for an effective post-disturbance management that optimises investments of resources and mitigates the risks due to shallow landslides.
The evaluation of RR-recovery depends on the original conditions of a stand too. For example, Liu et al. [10] state that recovery in forest gaps is faster than in closed forest, whereby the term of comparison is completely different. For instance, in forest opening root reinforcement is dominated by small roots that regrow faster after disturbances due to grass regrowth and tree regeneration. Whereas, in old and close forest stands, decayed coarse roots need much more time to be replaced. Moreover, several publications confirm that the contribution of fine roots to the overall root reinforcement is negligible if coarse roots are present [10,41].

Conclusions
The results of this work show the possible range of root reinforcement dynamic calculated based on existing forest stands (Swiss NFI) in Switzerland. We do not know the "story" of each stand, but in general the upper boundary of the calculated root reinforcement can be considered as the maximum potential condition that an "ideal" forest can reach in function of time. In this case, the definition of the "ideal" forest considering the upper boundary of the estimated root reinforcement is considerably different from the one applied in the Swiss guidelines for protection forest management. While the latter is defined on the basis of ecological and structural characteristics of a forest stand (such as tree stability and the presence of renovation), in this work we propose a more quantitative approach based on the effect of the vegetation on the hazard process. This allows a better characterisation of the time required for a forest stand to recover its protection function and to evaluate the effectiveness of the protection function with respect to the potential of a stand. In reality, the recovery of the protective function of the forest is often influenced by many other factors (almost independently of the characteristics of the stand) such as the concurrence with other species (in some cases neophytes), the damages due to ungulates, the presence of seed trees, and the mobility of seeds. Funding: This research was funded by the "Wald und Holz Forschungsfond" within the scope of the Project "Rutschungen nach Waldbrand" and the "Waldbrandmanagement auf der Alpennordseite AWN-1" founded by the Wyss Academy Fundation for Nature at the University of Bern.  [36] Juvenil growth data = [39].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: