Examining approaches for modeling individual tree growth response to thinning in Norway spruce

Christin Kuehne,Aron R.Weiskittel,Aksel Grnhus

a Norwegian Institute of Bioeconomy Research,Division of Forestry and Forest Resources,P.O.Box 115,NO-1431Ås,Norway

b University of Maine,Center for Research on Sustainable Forests,5755 Nutting Hall,Orono,ME,04469,USA

Keywords:Growth modeling Treatment response functions Multiplicative modifiers Picea abies Norway

A B S T R A C T 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.

1.Introduction

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¨anen 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 rubraSarg.)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 firstgeneration spruce stands.Climatic conditions vary in the study area with mean annual temperature ranging from approximately 0.5–°7C.–27 °.2C a°Cn da nmde mane aann naunanlu palr epcriepciitpaittiaotnio vna vrayrinyign fgr ofrmom sl isglihgthlytl yb ebleolwow 4 0400 m0 mmm 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).

Fig.1.Locations of Norway spruce forest experimental sites used in this study.

2.2.Data

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 sylvestrisL.),birch(Betula pendulaRoth,B.pubescensEhrh.),and goat willow(Salix capreaL.).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 1Average,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.

2.4.Data analysis

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ƒBASEequations:

whereYis 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 equationiand explanatory variablej.Random effects and residuals of the derived models were assumed to be normally distributed.Explanatory variables ofXβ 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,ƒTHINmodify the outcome ofƒBASE(i.e.ƒBASE׃THINorƒBASE+ƒTHINfor 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:

where γiare 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 BAREMand BABEFORE(Kuehne et al.,2016):

We further examined other multiplicative modifiers including a slightmodification of Eq.5 in Gyawali and Burkhart(2015):

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 γ41derived 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:

Similarly,the modifier by Short and Burkhart(1992)follows the same exponential response trajectory:

If complex modifiers as in Eqs.2–6 could not be fitted in a biologically plausible manner and/or estimated parameters γisignaling insignificance,we tested whether TST could be fitted intoƒBASEby iteratively adding indicator variables representing varying time periods after treatment with e.g.TST2indicating the second year after thinning and TST3‒5representing 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

2.5.Model evaluation

We calculated mean error(ME),mean absolute error(MAE),and relative MAE(MAE%)to evaluate and compare model prediction accuracy:

whereyiis the observed DBH,HT,or HCB at the end of the measurement period,respectively,^yiis the predicted DBH,HT,or HCB,respectively,andnis the number of observations.

3.Results

3.1.Diameter increment(ΔDBH)

Besides DBH,the finalΔDBHƒBASEalso 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 BAREMand BABEFOREfurther improved prediction accuracy and resulted in the following best performingƒFULL(Table 2,Fig.2):

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.

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ƒFULLover the explanatory variables including TST or other potentially influential factors such as age,elevation,and latitude(Fig.3).

Table 2Fixed 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 3Prediction 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.

3.2.Height increment(ΔHT)

Besides HT,the finalΔHTƒBASEalso included CR,BAL,and SI as explanatory variables.Adding a thinning effect toƒBASEresulted in mixed and contradicting results with indicator variables representing varying time periods after thinning often suggesting a slight increase in height

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

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 ƒFULLincluded a multiplicative modifier following Eq.4 predicting a marginal decrease in height increment for a time period of approximately 30 years following thinning(Tables 4 and 5,Fig.2):

There was no obvious bias when plotting the residuals ofΔHTƒBASE(Supplementary Material Fig.S3)andƒFULLover TST and any other explanatory variable or environmental factor(Fig.4).

3.3.Height to crown base increment(ΔHCB)

The finalΔHCBƒBASEincluded 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 nonthinned and thinned stands(data not shown).Adding a modifier following Eq.5 produced the best performingƒFULLof the following form(Tables 6 and 7,Fig.2):

Table 4Fixed 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 5Prediction 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 6Fixed 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 7Prediction 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ƒFULLover the explanatory variables and other environmental factors(Fig.5).

4.Discussion

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.

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

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¨akinen 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¨akinen and Isom¨aki,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.

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

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‒7of 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¨atalo 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).

Fig.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)evenaged 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).

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 growthaltering thinning effects and thus to build reliable models.

Funding

This research was financially supported by The Research Council of Norway(Norges Forskningsrådet,Project#301745).

Author contributions

CK:Data curation;Methodology;Formal analysis;Writing-original draft.ARW:Formal analysis;Writing-review & editing.AG:Project administration;Funding acquisition;Writing-review&editing.

Declaration of competing interest

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.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.i.org/10.1016/j.fecs.2022.100060.