Christian Kuehne, Aaron R. Weiskittel, Aksel Granhus. Examining approaches for modeling individual tree growth response to thinning in Norway spruce[J]. Forest Ecosystems, 2022, 9(1): 100060. DOI: 10.1016/j.fecs.2022.100060
Citation:
Christian Kuehne, Aaron R. Weiskittel, Aksel Granhus. Examining approaches for modeling individual tree growth response to thinning in Norway spruce[J]. Forest Ecosystems, 2022, 9(1): 100060. DOI: 10.1016/j.fecs.2022.100060
Christian Kuehne, Aaron R. Weiskittel, Aksel Granhus. Examining approaches for modeling individual tree growth response to thinning in Norway spruce[J]. Forest Ecosystems, 2022, 9(1): 100060. DOI: 10.1016/j.fecs.2022.100060
Citation:
Christian Kuehne, Aaron R. Weiskittel, Aksel Granhus. Examining approaches for modeling individual tree growth response to thinning in Norway spruce[J]. Forest Ecosystems, 2022, 9(1): 100060. DOI: 10.1016/j.fecs.2022.100060
Using periodic measurements from permanent plots in non-thinned and thinned Norway spruce (Picea abies (L.) H. Karst.) stands in Norway, individual-tree growth models were developed to predict annual diameter increment, height increment, and height to crown base increment. Based on long-term data across a range of thinning regimes and stand conditions, alternative approaches for modeling response to treatment were assessed. Dynamic thinning response functions in the form of multiplicative modifiers that predict no effect at the time of thinning, a rapid increase followed by an early maximum before the effect gradually declines to zero could not be fitted to initially derived baseline models without thinning related predictors. However, alternative approaches were used and found to perform well. Specifically, indicator variables representing varying time periods after thinning were statistically significant and behaved in a robust manner as well as consistent with general expectations. In addition, they improved overall prediction accuracy when incorporated as fixed effects into the baseline models for diameter and height to crown base increment. Further, more simply, including exponentially decreasing multiplicative thinning response functions improved prediction accuracy for height increment and height to crown base increment. Irrespective of studied attribute and modelling approach, improvement in performance of these extended models was relatively limited when compared to the corresponding baseline models and more pronounced in trees from thinned stands. We conclude that the largely varying and often multi-year measurement intervals of the periodic data used in this study likely prevented the development of more sophisticated thinning response functions. However, based on the evaluation of the final models' overall performance such complex response functions may not to be necessary to reliably predict individual tree growth after thinning for certain conditions or species, which should be further considered in future analyses of similar nature.
Regulating stand density through thinning is one of the most common silvicultural approaches to achieve forest management objectives (Zeide, 2001). The proper removal of trees through thinning enhances the growth and vigor of the remaining trees in the stand. Moreover, regulation of stand density through initial spacing or thinning interventions can also affect the trees' resistance towards windthrow (Gardiner et al., 1997; Achim et al., 2005). The risk of windthrow is of particular importance to shallow-rooted tree species like Norway spruce (Picea abies (L.) Karst), which in Norway is by far the commercially most important tree species, making up near one half of the growing stock by volume (Svensson et al., 2021). With a warming climate, less soil frost in winter and hence poorer root anchorage might render trees more susceptible, even if the frequency and intensity of storms do not change from the current-day situation. Yet another climate related factor to consider that comes along with warmer winters in boreal and/or mountainous regions, is the risk posed by wet snow events, which may cause extensive stem breakage and/or overturning of trees, occasionally destroying entire stands (Nykänen et al., 1997). In view of the ongoing climate change and the above-mentioned associated risks of abiotic damage, there is a need to develop density management guidelines for Norway spruce stands that can aid managers to balance production targets and stability goals.
To assess the risks of windthrow and/or snow damage in response to different spacing or thinning regimes using current mechanistic models such as e.g. HWIND (Peltola et al., 1999) and ForestGALES (Gardiner et al., 2000), information is needed on how certain stability-related tree-level variables are affected by silvicultural treatments (Gardiner et al., 2008). Important variables are amongst other tree height and a suite of other metrics that are strongly correlated with the resistance to uprooting or breakage, such as tree diameter and the height to diameter ratio. Metrics describing the size of the tree crown, such as crown volume or the ratio of crown length to total tree height, are needed as these characteristics correspond well with the wind- and/or snow-loading experienced by the trees. In turn, these metrics are largely shaped by the initial spacing and subsequent thinning(s) in the stand. While investigations of thinning in Norway spruce stands by use of long-term permanent trials have been conducted in Norway since the late 19th century (Andreassen et al., 2018), individual tree growth models which explicitly incorporate the effects of thinning, have yet to be developed for Norwegian forests (Andreassen and Tomter, 2003). The development of such models is crucial for understanding post-thinning individual tree and forest growth as well as for understanding how different stability related metrics are affected under different stand density and thinning scenarios.
Alternative approaches for modeling tree-level response to thinning have been explored (e.g. Weiskittel et al., 2011). For example, Kuehne et al. (2016) developed detailed thinning response functions for red spruce (Picea rubra Sarg.) and balsam fir (Abies balsamea [L.] Mill.) to accurately predict and quantify change in tree growth after thinning when compared to growth in trees from similar but non-thinned stands, i.e. stands of similar age, species composition, density, and site productivity (Weiskittel et al., 2011). The derived dynamic response functions in Kuehne et al. (2016) are continuous and able to exactly predict growth at any point in time after thinning mainly because of the annual resolution of the data used to build the functions. Such exhaustive data, however, is often not available given the large amount of resources necessary to collect them. Instead, permanent forest inventory or experimental plots are often monitored periodically with measurement intervals usually stretching five to ten years. Sophisticated approaches have been developed in recent years to build reliable stand- and tree-level growth and yield models based on such periodic data (McDill and Amateis, 1993; Amaro et al., 1998; Cieszewski and Bailey, 2000; Flewelling and Monserud, 2002). Among them, annualizing periodic observations has been proven to be a reliable method for modelling individual tree growth (e.g. Weiskittel et al., 2007; Kuehne et al., 2020). However, it is questionable whether such an annualization approach is capable to detect and quantify the subtle differences in growth between trees from thinned and non-thinned stands when using periodic forest growth data with varying multi-year measurement intervals.
Consequently, this work aimed to develop individual-tree models for diameter and height increment as well as crown recession using data from periodically monitored permanent plots in thinned and non-thinned Norway spruce stands in Norway. We aimed to address three research questions, namely whether i) complex thinning response functions can be fitted to data from periodic measurements with widely varying growth intervals, ii) more simple thinning modifiers including thinning indicator variables can be used in case the approach outlined in i) does not work, and iii) complex response functions are generally necessary to reliably predict individual tree growth after thinning.
2.
Material and methods
2.1
Study area
Majority of the experimental stands used in this study were distributed within and covering the natural distribution range of Norway spruce in Norway, which encompasses most of southeastern and central Norway as well as northern Norway up to approximately 66°N latitude (Fig. 1). While most trials were in stands previously dominated by Norway spruce, the trials in northern Norway north of 66°N were established in first-generation spruce stands. Climatic conditions vary in the study area with mean annual temperature ranging from approximately 0.5 ℃–7.2 ℃ and mean annual precipitation varying from slightly below 400 mm and up to 1, 800 mm (1981–2010, Meteorological Institute of Norway). The current study excluded data from the southwestern and western counties with their maritime climate (Fig. 1), as growing conditions in oceanic western Norway are considered to be very different to the rest of the country (Vestjordet, 1967; Bauger, 1995).
Figure
1.
Locations of Norway spruce forest experimental sites used in this study.
Data used in this study were from long-term forest management trials maintained by the Norwegian Institute of Bioeconomy Research (Andreassen et al., 2018). Majority of these silvicultural experiments was established in pure, even-aged Norway spruce stands in the early 1970s as part of a larger thinning trial network (Braastad and Tveite, 2001). These trails were initially established to examine the effects of thinning from below with focus on thinning intensity and timing of thinning on stand-level volume production in 18- to 47-year-old stands. Treatments included one to three thinnings from below with different combinations of target residual trees per hectare (TPH) of 2070, 1600, 1100, and 800 at stand dominant heights of 8, 12, 16, and 20 m. However, deviations from the original experimental design resulted in a diverse dataset covering a wide range of thinning regimes varying in number of thinning interventions, time of thinning, and thinning intensity (Supplementary Material Fig. S1). We further added measurements from experimental trials from spruce stands that were never thinned. Although the trials in non-thinned stands often were established prior 1970, we only included measurements from that year onwards. At that time, these non-thinned stands were between 12 and 111 years old. Some of the plots used here contained individual trees of species other than Norway spruce, including Scots pine (Pinus sylvestris L.), birch (Betula pendula Roth, B. pubescens Ehrh.), and goat willow (Salix caprea L.). However, Norway spruce accounted for an average of 99% of total plot basal area over all plots. Plots with excessive mortality as a result of major windstorms were only included up until the occurrence of the storm event.
The final dataset comprised a total of 116, 773 diameter increment (with 41, 870 from non-thinned and 74, 903 from thinned stands), 43, 850 height increment (13, 901 from non-thinned and 29, 949 from thinned stands), and 33, 499 (8, 141 from non-thinned and 25, 358 from thinned stands) height to crown base increment (crown recession) observations (growth intervals) from 204 permanent plots in 60 trials. The length of the studied growth intervals varied between 1 and 42 years with an average of 6 years. Measurement records for all trees comprised species, vital status (live or dead) and diameter at breast height (DBH, cm). Total tree height (HT, m) and height to crown base (HCB, m) were measured for a randomly selected subsample of all live trees during every measurement campaign with about half of all trees measured for total height on average across all plots.
2.3
Data preparation
HT measurements missing in the final database were imputed using a total tree height equation similar to model 6 in Temesgen et al. (2007) and a mixed-modeling approach with trial as random effect. Similarly, HCB measurements missing in the final database were imputed using an equation similar to the one in Kuehne et al. (2016). Individual tree crown ratio (CR) was then calculated as the ratio of (live) crown length (HT-HCB) and HT. Basal area in larger trees (BAL, m2·ha−1) was also computed for each individual tree (Kuehne et al., 2019). Individual tree measurements were then summarized for each plot to calculate stand-level metrics including dominant height (mean height of the 100 thickest trees·ha−1, HTDOM, m) and total basal area (BA, m2·ha−1) (Table 1). Site index (SI, m) at base age 40 years was determined using the dominant height model from Allen et al. (2020). We further calculated treatment related variables including time since thinning (TST), BA (BABEFORE, BAAFTER) and TPH before as well as after each thinning intervention (TPHBEFORE, TPHAFTER), respectively, and also quantified dominant height (HTTHIN) and stand age (AGE) at the time of thinning (AGETHIN) when plot records indicated a thinning took place. An overview of the stand- and tree-level variables for each dataset used in this work is provided in Table 1.
Table
1.
Average, standard deviation (SD), and minimum-maximum of individual tree- and stand-level attributes at the beginning of the measurement period of the analyzed Norway spruce diameter increment (ΔDBH, n = 116, 773), height increment (ΔHT, n = 43, 850), and height to crown base increment (ΔHCB, n = 33, 499) datasets of this study.
We started our analysis following the approach outlined in Kuehne et al. (2016). First, a baseline model (ƒBASE) not specifically accounting for thinning, i.e. not including treatment related variables was developed for all three studied response variables, i.e. annual diameter increment (ΔDBH), annual height increment (ΔHT), and annual height to crown base increment (ΔHCB). Using all observations, i.e. measurements from thinned and non-thinned stands, the following general model form was used to derive the sought-after ƒBASE equations:
Y=exp(Xβ)
(1)
where Y is the response variable (ΔDBH, ΔHT, or ΔHCB), Xβ is the model-specific explanatory variable design matrix (linear predictor; Zuur et al., 2009) with the associated estimated fixed (βij) and trial-specific random parameters (bij) for equation i and explanatory variable j. Random effects and residuals of the derived models were assumed to be normally distributed. Explanatory variables of Xβ comprised the previously described tree- and stand-level attributes which were added to each model in a stepwise manner depending on significance and biological plausibility.
Next, we aimed to fit multiplicative thinning modifiers (ƒTHIN) following a Type 1 response as described in Snowdon (2002). Also called treatment response functions, ƒTHIN modify the outcome of ƒBASE (i.e. ƒBASE × ƒTHIN or ƒBASE + ƒTHIN for additive modifiers; Weiskittel et al., (2011) based on treatment variables with time since thinning (TST) as the most influential one. A Type 1 thinning modifier should predict a relatively rapid increase and often early maximum of the thinning effect shortly after the thinning followed by a more gradual lessening of the response over time after thinning until the thinning effect eventually diminishes and becomes zero again. The thinning effect is supposed to be positive for ΔDBH and potentially negative for ΔHT and ΔHCB (Kuehne et al., 2016). To correspond with the general response pattern for thinning interventions, we used a Type I combined exponential - power function of the following form:
fTHIN=1±(γ21×γTST22×TSTγ23)
(2)
where γi are fixed effect parameters and TST = 0 for non-thinned plots. To account for thinning intensity, Eq. 2 can be extended by incorporating e.g. the ratio of BAREM and BABEFORE (Kuehne et al., 2016):
The response predicted in Eq. 4 roughly follows the one outlined for Eq. 2 but depicts a more gradual initial increase and also needs to be restrained according to the duration parameter γ41 derived during the fitting process (Liu et al., 1995).
In contrast to Eqs. 2-4, the multiplicative thinning modifier by Hann et al. (2003) predicts a continuous, exponential treatment effect approaching 1 with increasing TST after a maximum response right after the thinning intervention:
fTHIN=1±[γ50×(BAREMBABEFORE)γ51×exp(−γ52×TST)]
(5)
Similarly, the modifier by Short and Burkhart (1992) follows the same exponential response trajectory:
fTHIN=(BAAFTERBABEFORE)(γ61×(AGETHINAGE))
(6)
If complex modifiers as in Eqs. 2-6 could not be fitted in a biologically plausible manner and/or estimated parameters γi signaling insignificance, we tested whether TST could be fitted into ƒBASE by iteratively adding indicator variables representing varying time periods after treatment with e.g. TST2 indicating the second year after thinning and TST3‒5 representing years three to five after thinning (Bianchi et al., 2022). Final models including a thinning modifier or TST indicator variables are referred to full models (ƒFULL) hereafter.
All models were derived using nonlinear mixed effects modelling and the package 'nlme' (Pinheiro et al., 2021) of the statistical computing software R (version 4.1.2, R Development Core Team, 2021). A variance structure was incorporated to account for the variability in the most influential explanatory variable of each equation (DBH, HT, and CR for ΔDBH, ΔHT, and ΔHCB, respectively). To overcome problems of the varying measurement intervals observed in our datasets, parameters were annualized using an iterative mixed-effects technique of Weiskittel et al. (2007). Based on Cao (2000), the right side of the equation was a function that summed the predicted annualized ΔDBH, ΔHT or ΔHCB estimates, respectively, over the number of growing seasons during the observed growth period using the updated parameter estimates from the optimization algorithms. For each growing season during the growth period, DBH, HT or HCB, respectively, was subsequently updated using the annual increment estimates, while all other explanatory variables were linearly interpolated between their beginning values and ending values. SI, however, was assumed to be constant over time.
2.5
Model evaluation
We calculated mean error (ME), mean absolute error (MAE), and relative MAE (MAE%) to evaluate and compare model prediction accuracy:
ME=∑ni=1(yi−ˆyi)n
(7a)
MAE=∑ni=1(|yi−ˆyi|)n
(7b)
MAE%=∑ni=1((|yi−ˆyi|)yi×100)n
(7c)
where yi is the observed DBH, HT, or HCB at the end of the measurement period, respectively, ˆyi is the predicted DBH, HT, or HCB, respectively, and n is the number of observations.
3.
Results
3.1
Diameter increment (ΔDBH)
Besides DBH, the final ΔDBH ƒBASE also included CR, BA, BAL, and SI as explanatory variables. Fitting complex thinning modifier as in Eqs. 2-4 did not result in biologically plausible outcomes, while a modifier following Eq. 5 produced reasonable results (data not shown). However, incorporating indicator variables representing the third year, years four and five as well as six and seven after thinning, respectively, were significant and improved model performance. Extending these thinning indicator variables to account for thinning intensity by multiplying the ratio of BAREM and BABEFORE further improved prediction accuracy and resulted in the following best performing ƒFULL (Table 2, Fig. 2):
Table
2.
Fixed effects parameter estimates and statistics of the annual tree breast height diameter increment (ΔDBH, cm·yr−1) baseline (ƒBASE) and full (ƒFULL) mixed effects models for Norway spruce in Norway without and with indicator variables for time since thinning, respectively.
Figure
2.
Individual tree-level annual DBH increment (ΔDBH), annual total height increment (ΔHT), and annual height to crown base increment (ΔHCB) thinning response functions for Norway spruce growing in Norway. Curves were derived for an exemplary thinning removing 20% of total stand-level basal area.
While adding the thinning indicator variables slightly improved prediction accuracy overall and for observations from thinned stands (Table 3), we found no obvious bias when plotting the residuals of ΔDBH ƒBASE (Supplementary Material Fig. S2) and ΔDBH ƒFULL over the explanatory variables including TST or other potentially influential factors such as age, elevation, and latitude (Fig. 3).
Table
3.
Prediction accuracy statistics including mean error (ME), mean absolute error (MAE) and relative MAE (MAE%) for DBH at the end if the measurement period derived from the Norway spruce individual tree annual diameter increment models without (ƒBASE) and with indicator variables for time since thinning (ƒFULL) for observations from non-thinned, thinned and all stands.
Figure
3.
Standardized residuals for individual tree diameter (cm) at the end of the measurement period over predicted diameter, diameter at the beginning of the measurement period, time since thinning (TST) at the end of the measurement period, and latitude for thinned and non-thinned (TST = 0) even-aged stands of Norway spruce growing in Norway. Residuals were derived from the diameter increment model with indicator variables for TST.
Besides HT, the final ΔHT ƒBASE also included CR, BAL, and SI as explanatory variables. Adding a thinning effect to ƒBASE resulted in mixed and contradicting results with indicator variables representing varying time periods after thinning often suggesting a slight increase in height increment after thinning (data not shown). Fitting a response function following Eq. 5 resulted in a predicted limited negative thinning effect lasting over 50 years (data not shown). The best performing plausible ƒFULL included a multiplicative modifier following Eq. 4 predicting a marginal decrease in height increment for a time period of approximately 30 years following thinning (Table 4 and 5, Fig. 2):
Table
4.
Fixed effects parameter estimates and statistics of the annual total tree height increment (ΔHT, m·yr−1) baseline (ƒBASE) and full (ƒFULL) mixed effects models for Norway spruce in Norway without and with a thinning modifier (response function) following Gyawali and Burkhart (2015), respectively.
Table
5.
Prediction accuracy statistics including mean error (ME), mean absolute error (MAE) and relative MAE (MAE%) for total tree height at the end if the measurement period derived from the Norway spruce individual tree annual height increment models without (ƒBASE) and with a thinning modifier (response function) following Gyawali and Burkhart (2015) (ƒFULL) for observations from non-thinned, thinned and all stands.
There was no obvious bias when plotting the residuals of ΔHT ƒBASE (Supplementary Material Fig. S3) and ƒFULL over TST and any other explanatory variable or environmental factor (Fig. 4).
Figure
4.
Standardized residuals for individual total tree HT (m) at the end of the measurement period over predicted total tree height, total tree height at the beginning of the measurement period, time since thinning (TST) at the end of the measurement period, and latitude for thinned and non-thinned (TST = 0) even-aged stands of Norway spruce growing in Norway. Residuals were derived from the height increment model with a thinning modifier (response function) following Gyawali and Burkhart (2015).
The final ΔHCB ƒBASE included CR, HT, BA, and SI as explanatory variables. Fitting complex thinning modifier as in Eqs. 2-4 did not result in biologically plausible outcomes. However, an indicator variable representing the first eight years after thinning (TST1‒8 = −0.113047) was significant and slightly improved prediction accuracy for trees from non-thinned and thinned stands (data not shown). Adding a modifier following Eq. 5 produced the best performing ƒFULL of the following form (Table 6 and 7, Fig. 2):
Table
6.
Fixed effects parameter estimates and statistics of the annual height to crown base increment (ΔHCB, m·yr−1) baseline (ƒBASE) and full (ƒFULL) mixed effects models for Norway spruce in Norway without and with a thinning modifier (response function) following Hann et al. (2003), respectively.
Table
7.
Prediction accuracy statistics including mean error (ME), mean absolute error (MAE) and relative MAE (MAE%) for total tree height at the end if the measurement period derived from the Norway spruce individual tree height to crown base increment models without (ƒBASE) and with a thinning modifier (response function) following Hann et al. (2003) (ƒFULL) for observations from non-thinned, thinned and all stands.
There was no obvious major bias when plotting the residuals of ΔHCB ƒBASE (Supplementary Material Fig. S4) and ƒFULL over the explanatory variables and other environmental factors (Fig. 5).
Figure
5.
Standardized residuals for individual tree height to crown base (m) at the end of the measurement period over predicted height to crown base, crown ratio at the beginning of the measurement period, time since thinning (TST) at the end of the measurement period, and latitude for thinned and non-thinned (TST = 0) even-aged stands of Norway spruce growing in Norway. Residuals were derived from the height to crown base increment model with a thinning modifier (response function) following Hann et al. (2003).
This study developed dynamic growth models to predict various individual tree attributes from periodic measurements in non-thinned and thinned Norway spruce stands in Norway. Irrespective of the attribute examined, baseline models without thinning treatment related predictors performed well showing no bias over time since thinning (TST). However, adding thinning related variables either as simple indicator variables for diameter increment (ΔDBH) or multiplicative thinning modifiers (response functions) for height increment (ΔHT) and height to crown base increment (ΔHCB) slightly improved overall prediction accuracy when compared to the corresponding baseline models. Sophisticated modifiers that predict a rapid increase, an early peak, and gradual decline of the thinning effect could not be fitted here.
Multiplicative treatment modifiers are the preferred approach to predict treatment responses because they are potentially the most accurate and detailed way of modelling silvicultural treatment effects including thinning. Multiplicative thinning modifiers as in Eqs. 2-6 do not alter the corresponding baseline model but rather function as modifiers that adjust the outcome of the baseline model in a proportional manner. The resulting relative effect of the modelled treatment response is independent of the baseline equation and thus constant but depends on the outcome of the baseline function in absolute terms. If deemed appropriate, thinning response functions thus can be used in other, separate applications individually, i.e. independent from the corresponding co-developed baseline model.
Continuous and biologically plausible modifiers that predict no effect at the time of thinning, a swift increase followed by an early maximum before the effect gradually declines and becomes zero again, could not be fitted in our study (Kuehne et al., 2016). The modified thinning response function after Liu et al. (1995) fitted to ΔHT somewhat followed the outline response curve but depicted a more symmetrical trajectory over TST. Further, the modelled response effect further needs to be retained after 32 years as it does not approach zero with increasing TST. We believe the experienced difficulties in fitting dynamic modifiers in our study were a result of the periodic data used with its varying multi-year growth intervals that likely concealed the subtle growth differences between trees of non-thinned and thinned stands after accounting for stand density. Previous similar studies that were able to fit such sophisticated modifiers as outlined above either used annual measurements (Kuehne et al., 2016) or short growth intervals of constant, fixed time periods (Liu et al., 1995; Gyawali and Burkhart. 2015).
The final, best performing multiplicative modifiers fitted to our ΔHT and ΔHCB baseline models including an exponentially decreasing thinning effect after a maximum right after thinning for ΔHCB (Hann et al., 2003) seem to be reasonable for ΔHCB and not unlikely for ΔHT. The ideal, i.e. most realistic modifier for ΔHT, however, would potentially depict a delay in the thinning effect of one year because of the partial predefinition of height growth in the growing season prior to thinning. Based on our findings, a basal area removal of 20% would reduce ΔHT and ΔHCB by approximately 5% and 29% approximately nine years after thinning and in the first year following thinning, respectively. Comparable (maximum) values have been reported in similar studies (Liu et al., 1995; Hann et al., 2003; Mäkinen et al., 2006; Kuehne et al., 2016). According to Weiskittel et al. (2011), thinning can either reduce, maintain, or increase tree height growth and the varying behavior is driven by species, crown class, age, and stand density prior to thinning. Our mixed results for ΔHT with indicator variables for various time periods after thinning often signaling a minor increase in ΔHT seem to corroborate this. Insensitivity of ΔHT to thinning has been reported in several previous works (Hynynen, 1995; Westfall and Burkhart, 2001; Hann et al., 2003; Mäkinen and Isomäki, 2004). Given the rather small magnitude and the very short duration, Sharma et al. (2006) and Gyawali and Burkhart (2015) concluded that the thinning induced changes in (dominant) ΔHT observed in their studies were neglectable. Such a conclusion also seems reasonable for our study. A potential decrease in ΔHT and ΔHCB is likely triggered by a reallocation of resources with lateral crown expansion favored over leader growth. In addition, parts if not the entire crown need to adjust to a significantly altered light regime causing a thinning shock to formerly, i.e. pre-thinning shaded foliage.
Although multiplicative response functions could not be fitted for ΔDBH in this study, the indicator variables for TST that were successfully incorporated into the final ΔDBH full model still are proportional and multiplicative in nature. This is a result of the exponential model form used to analyze and predict ΔDBH. TST3, TST4‒5, and TST6‒7 of the full ΔDBH model predict an annual increase in ΔDBH by 40%, 12%, and 8%, respectively, for a thinning intervention removing 20% of pre-treatment basal area. The general magnitude of the modelled ΔDBH thinning effect of our study thus appears to be approximately in line with results from previous similar works (Hynynen, 1995; Jonsson, 1995; Hann et al., 2003; Kuehne et al., 2016). However, the thinning treatment effect usually lasted longer, i.e. up to 15–25 years, in these other studies while a similar delay of the thinning response as found here was reported previously (Thorpe et al., 2007; Mehtätalo et al., 2014). Such a delay is likely to be caused by physiological acclimation to the altered post-thinning growing conditions and the associated resource allocation to root and shoot growth.
Irrespective of the attribute studied here, the residuals of all baseline models without a thinning modifier or thinning indicator variables, respectively, showed no obvious bias when plotted over TST. This was in part to be expected as the main effect of a thinning intervention, i.e. the reduction in stand density is accounted for in each of the baseline models developed here (Weiskittel et al., 2011). All baseline models of this study include at least one metric of competition, i.e. basal area (BA) and/or basal area in larger trees (BAL), representing site occupation and/or stand density. Consequently, a reduction in stand density as a result of a thinning is automatically reflected in the baseline model input variables through altered competition values. Various studies at the stand- as well as the individual tree-level thus found no need to incorporate thinning modifiers into dynamic growth models (Hynynen, 1995; Barrio Anta et al., 2006). However, whether baseline models accurately predict growth of non-thinned as well as thinned stands and trees thereof appears to depend on the kind of forest or species and attribute examined (Weiskittel et al., 2011).
Thinning modifiers or indicator variables representing certain time periods after thinning as in this study still have been examined and successfully added to ΔDBH and stand-level basal area increment baseline equations (Bailey and Ware, 1983; Amateis et al., 1989; de Souza et al., 2022) because of additional, minor effects other than stand density reduction associated with thinning activities (Weiskittel et al., 2011). Among them, fertilization and selection effects appear to be most influential (Hynynen, 1995). Decomposition of organic material including logging slash is often accelerated as a more open canopy following thinning results in a more favorable microclimate on the ground leading to a greater nutrient supply for a limited time. In soil moisture limited forest ecosystems, thinnings also have the potential to improve water availability for residual trees. Further, conventional thinning interventions aim to remove unhealthy, damaged, and thus less vigorous trees. Residual trees not only grew at higher rates before thinning but are also able to better respond to the modified post-thinning growing conditions.
Likely because the baseline models derived in our study are able to capture the main thinning effect, i.e. the reduction in stand density, improvement in model prediction accuracy was mostly rather marginal when extending the models by adding treatment related predictors irrespective of the attribute and approach, which is in agreement with previous studies (Kuehne et al., 2016). While marginal on an annual basis, improvements in prediction accuracy are likely to sum and scale up for longer-term simulation runs of several decades – particularly at the stand level (Kuehne et al., 2016). Still, predictions of our baseline models appear to be reliable because we found no obvious major bias when plotting residuals over TST or any other potentially influential factor. This suggests that complex treatment response functions might not be needed for certain conditions or species, which highlights the often complex and highly varied nature of tree growth following thinning (e.g. Bose et al., 2018).
We thus conclude that i) analyzing periodic growth data of varying multi-year measurement intervals using an annualization method proved to be successful in the detection and quantification of thinning effects not captured in the modification of competition variables and ii) sophisticated, continuous thinning modifiers as initially aimed for here may not be necessary to predict these other individual tree growth-altering thinning effects and thus to build reliable models.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
The authors are indebted to the many people who have been involved in the installment and monitoring of the permanent sample plots analyzed in this work. We also like to thank the reviewer for constructive comments on the previous version of this paper and the associate editor for handling the manuscript.
Achim, A., Ruel, J. -C., Gardiner, B., Laflamme, G., Meunier, S., 2005. Modelling the vulnerability of balsam fir forest to wind damage. For. Ecol. Manag. 204, 35-50.
[2]
Allen, Ⅱ, M.G., Antón-Fernández, C., Astrup, R., 2020. A stand-level growth and yield model for thinned and unthinned managed Norway spruce forests in Norway. Scand. J. For. Res. 35(5-6), 238-251.
[3]
Amaro, A., Reed, D., Tomé, M., Themido, I., 1998. Modeling dominant height growth: Eucalyptus plantations in Portugal. For. Sci. 44, 37-46.
[4]
Amateis, R.L., Burkhart, H.E., Walsh, T.A., 1989. Diameter increment and survival equations for loblolly pine trees growing in thinned and unthinned plantations on cutover, site-prepared lands. South. J. Appl. For. 13(4), 170-174.
[5]
Andreassen, K., Tomter, S.M., 2003. Basal area growth models for individual trees of Norway spruce, Scots pine, birch and other broadleaves in Norway. For. Ecol. Manag. 180(1-3), 11-24.
[6]
Andreassen, K., Allen, M., Holt Hanssen, K., Kvaalen, H., Støtvig, S., 2018. Langsiktig feltforsøk i skog ved NIBIO. [NIBIO's long-term forest field trials]. NIBIO Rapport 4/61/2018. (In Norwegian).
[7]
Bailey, R.L., Ware, K.D., 1983. Compatible basal-area growth and yield model for thinned and unthinned stands. Can. J. For. Res. 13(4), 563-571.
[8]
Barrio Anta, M., Castedo, F., Diéguez-Aranda, U., Álvarez González, J.G., Parresol, B.R., Rodríguez, R., 2006. Development of a basal area growth system for maritime pine in northwestern Spain using the generalized algebraic difference approach. Can. J. For. Res. 36, 1461-1474.
[9]
Bauger, E. 1995. Funksjoner og tabeller for kubering av stående trær: furu, gran og sitkagran på Vestlandet (Tree volume functions and tables: Scots pine, Norway spruce and Sitka spruce in western Norway). Rapport fra Skogforsk 16/95. (In Norwegian).
[10]
Bianchi, S., Huuskonen, S., Hynynen, J., Siipilehto, J., Niemistö, P., 2022. Tree-level differences in Norway spruce and Scots pine growth after extreme thinning treatments. Scand. J. For. Res. .
[11]
Bose, A.K., Weiskittel, A., Kuehne, C., Wagner, R.G., Turnblom, E., Burkhart, H.E., 2018 Tree-level growth and survival following commercial thinning of four major softwood species in North America. For. Ecol. Manag. 427, 355-364.
[12]
Braastad, H., Tveite, B., 2001. Tynning i gran- og furubestand: effekt pa tynning pa volumproduksjon, middeldiameter og diameter av de 800 groevste traer per ha. [Thinning in spruce and pine stands: effect of thinning on volume production, average diameter, and the diameter of the largest 800 trees per hectare]. Rapport fra Norsk institutt for skogforskning 10/01. (In Norwegian).
[13]
Cao, Q.V., 2000. Prediction of annual diameter growth and survival for individual trees from periodic measurements. For. Sci. 46(1), 127-131.
[14]
Cieszewski, C.J., Bailey, R.L., 2000. Generalized algebraic difference approach: theory based derivation of dynamic site equations with polymorphism and variable asymptotes. For. Sci. 46, 116-126.
[15]
de Souza, H.J., Miguel, E.P., Nascimento, R.G.M., Cabacinha, C.D., Rezende, A.V., dos Santos, M.L., 2022. Thinning-response modifier term in growth models: An application on clonal Tectona grandis Linn F. stands at the Amazonian region. For. Ecol. Manag. 511, 120109.
[16]
Flewelling, J.W., Monserud, R.A., 2002. Comparing methods for modeling tree mortality, in NL Crookston and RN Havis (eds). Second Forest Vegetation Simulator Conference. Proc. RMRS-P-25, USDA Forest Service, Rocky Mountain Res. Sta., Ogden, UT, pp. 168-175.
[17]
Gardiner, B.A., Stacey, G. R, Belcher, R.E. Wood, C.J., 1997. Field and wind tunnel assessments of the implications of respacing and thinning for tree stability. Forestry 20(1), 233-252.
[18]
Gardiner, B., Peltola, H. Kellomäki, S., 2000. Comparison of two models for predicting the critical wind speeds required to damage coniferous trees. Ecol. Model. 129, 1-23.
[19]
Gardiner, B., Byrne, K., Hale, S., Kamimura, K., Mitchell, S.J., Peltola, H. Ruel, J. -C. 2008. A review of mechanistic modelling of wind damage risk to forests. Forestry 81(3), 447-463.
[20]
Gyawali, N., Burkhart, H.E., 2015. General response functions to silvicultural treatments in loblolly pine plantations. Can. J. For. Res. 45(3), 252-265.
[21]
Hann, D.W., Marshall, D.D., Hanus, M.L., 2003. Equations for predicting height-to-crown-base, 5-year diameter-growth rate, 5-year height-growth rate, 5-year mortality rate, and maximum size-density trajectory for Douglas-fir and western hemlock in the coastal region of the Pacific Northwest. Oregon State University, College of Forestry, Forest Research Laboratory, Corvallis, OR.
[22]
Hynynen, J., 1995. Predicting the growth response to thinning for Scots pine stands using individual-tree growth models. Silv. Fenn. 29, 225-246.
[23]
Jonsson, B., 1995. Thinning response functions for single trees of Pinus sylvestris L. and Picea abies (L. ) Karst. Scand. J. For. Res. 10(4), 353-369.
[24]
Kuehne, C., Weiskittel, A.R., Wagner, R.G., Roth, B.E., 2016. Development and evaluation of individual tree-and stand-level approaches for predicting spruce-fir response to commercial thinning in Maine, USA. For. Ecol. Manag. 376, 84-95.
[25]
Kuehne, C., Weiskittel, A.R., Waskiewicz, J., 2019. Comparing performance of contrasting distance-independent and distance-dependent competition metrics in predicting individual tree diameter increment and survival within structurally-heterogeneous, mixed-species forests of Northeastern United States. For. Ecol. Manag. 433, 205-216.
[26]
Kuehne, C., Russell, M.B., Weiskittel, A.R., Kershaw Jr, J.A., 2020. Comparing strategies for representing individual-tree secondary growth in mixed-species stands in the Acadian Forest region. For. Ecol. Manag. 459, 117823.
[27]
Liu, J., Burkhart, H.E., Amateis, R.L., 1995. Projecting crown measures for loblolly pine trees using a generalized thinning response function. For. Sci. 41 (1), 43-53.
[28]
Mäkinen, H., Isomäki, A., 2004. Thinning intensity and long-term changes in increment and stem form of Norway spruce trees. For. Ecol. Manag. 201(2-3), 295-309.
[29]
Mäkinen, H., Isomäki, A., Hongisto, T., 2006. Effect of half-systematic and systematic thinning on the increment of Scots pine and Norway spruce in Finland. Forestry 79(1), 103-121.
[30]
McDill, M.E., Amateis, R.L., 1993. Fitting discrete-time dynamic models having any time interval. For. Sci. 39(3), 499-519.
[31]
Mehtätalo, L., Peltola, H., Kilpeläinen, A., Ikonen, V.P., 2014. The response of basal area growth of Scots pine to thinning: A longitudinal analysis of tree-specific series using a nonlinear mixed-effects model. For. Sci. 60(4), 636-644.
[32]
Nykänen, M. -L., Peltola, H., Quine, C.P., Kellomäki, S., Broadgate, M., 1997. Factors affecting snow damage of trees with particular reference to European conditions. Silv. Fenn. 31(2), 193-213.
[33]
Peltola, H., Kellomäki, S., Väisänen, H., Ikonen, V.P., 1999. A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce and birch. Can. J. For. Res. 29, 647-661.
[34]
Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., 2021. nlme: Linear and nonlinear mixed effects models. R package version 3.1-153.
.
[35]
R Development Core Team, 2021. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.
.
[36]
Sharma, M., Smith, M., Burkhart, H.E., Amateis, R.L., 2006. Modeling the impact of thinning on height development of dominant and codominant loblloly pine trees. Ann. For. Sci. 63, 349-354.
[37]
Short Ⅲ, E.A., Burkhart, H.E., 1992. Predicting crown-height increment for thinned and unthinned loblolly pine plantations. For. Sci., 38(3), 594-610.
[38]
Snowdon, P. 2002. Modeling Type 1 and Type 2 growth responses in plantations after application of fertilizer or other silvicultural treatments. For. Ecol. Manag. 163(1): 229-244.
[39]
Svensson, A., Eriksen, R., Hylen, G., Granhus, A., 2021. Skogen i Norge. Statistikk over skogforhold og skogressurser i Norge for perioden 2015-2019 (Norway's forests. Statistics on forest conditions and forest resources in Norway for the period 2015-2019). NIBIO Rapport 7(142), 1-53. [In Norwegian]. .
[40]
Temesgen, H., Hann, D.W., Monleon, V.J., 2007. Regional height-diameter equations for major tree species of southwest Oregon. West. J. Appl. For. 22(3), 213-219.
[41]
Thorpe, H.C., Thomas, S.C., Caspersen, J.P., 2007. Residual-tree growth responses to partial stand harvest in the black spruce (Picea mariana) boreal forest. Can. J. For. Res. 37(9), 1563-1571.
[42]
Vestjordet, E., 1967. Funksjoner og tabeller for kubering av staende gran (Functions and tables for volume of standing Norway spruce trees) In: Alf Langsaeter: festskrift utgitt til aere for skogdirektoer, professor, dr. agr. Alf Langsaeter pa hans 70-ars dag. p. 539-570. (In Norwegian).
[43]
Weiskittel, A.R., Garber, S.M., Johnson, G.P., Maguire, D.A., Monserud, R.A., 2007. Annualized diameter and height growth equations for pacific northwest plantation-grown Douglas-fir, western hemlock, and red alder. For. Ecol. Manag. 250(3), 266-278.
[44]
Weiskittel A.R., Hann D.W., Kerhsaw Jr. J.A., Vanclay J.K., 2011. Forest growth and yield modeling. Wiley.
[45]
Westfall, J.A., Burkhart, H.E., 2001. Incorporating thinning response into a loblolly pine stand simulator. South. J. Appl. For. 25 (4), 159-164.
[46]
Zeide, B., 2001. Thinning and growth: a full turnaround. J. For. 99(1), 20-25.
[47]
Zuur, A.F., Leno, E.N., Walker, N.J., Saveliev, A.A., Smith, G.M., 2009. Mixed Effects Models and Extensions in Ecology with R. Springer, New York, NY, p. 574.
Table
1.
Average, standard deviation (SD), and minimum-maximum of individual tree- and stand-level attributes at the beginning of the measurement period of the analyzed Norway spruce diameter increment (ΔDBH, n = 116, 773), height increment (ΔHT, n = 43, 850), and height to crown base increment (ΔHCB, n = 33, 499) datasets of this study.
Table
2.
Fixed effects parameter estimates and statistics of the annual tree breast height diameter increment (ΔDBH, cm·yr−1) baseline (ƒBASE) and full (ƒFULL) mixed effects models for Norway spruce in Norway without and with indicator variables for time since thinning, respectively.
Table
3.
Prediction accuracy statistics including mean error (ME), mean absolute error (MAE) and relative MAE (MAE%) for DBH at the end if the measurement period derived from the Norway spruce individual tree annual diameter increment models without (ƒBASE) and with indicator variables for time since thinning (ƒFULL) for observations from non-thinned, thinned and all stands.
Table
4.
Fixed effects parameter estimates and statistics of the annual total tree height increment (ΔHT, m·yr−1) baseline (ƒBASE) and full (ƒFULL) mixed effects models for Norway spruce in Norway without and with a thinning modifier (response function) following Gyawali and Burkhart (2015), respectively.
Table
5.
Prediction accuracy statistics including mean error (ME), mean absolute error (MAE) and relative MAE (MAE%) for total tree height at the end if the measurement period derived from the Norway spruce individual tree annual height increment models without (ƒBASE) and with a thinning modifier (response function) following Gyawali and Burkhart (2015) (ƒFULL) for observations from non-thinned, thinned and all stands.
Table
6.
Fixed effects parameter estimates and statistics of the annual height to crown base increment (ΔHCB, m·yr−1) baseline (ƒBASE) and full (ƒFULL) mixed effects models for Norway spruce in Norway without and with a thinning modifier (response function) following Hann et al. (2003), respectively.
Table
7.
Prediction accuracy statistics including mean error (ME), mean absolute error (MAE) and relative MAE (MAE%) for total tree height at the end if the measurement period derived from the Norway spruce individual tree height to crown base increment models without (ƒBASE) and with a thinning modifier (response function) following Hann et al. (2003) (ƒFULL) for observations from non-thinned, thinned and all stands.
DownLoad:
Email Alerts
RSS Feeds