Density and diameter distributions of saplings in naturally regenerated and planted coniferous stands in Québec after various approaches of commercial thinning

Abstract

Key Message

A model describing species composition, density and diameter distribution of saplings was developed from operational inventory data. It could be used as an input into growth models calibrated exclusively with merchantable trees to correct some recruitment bias. Important differences in distributions were found between plantations and naturally regenerated stands. Longer-term monitoring would be required to observe the effects of thinning treatments on saplings.

Context

Saplings play important ecological and structural roles in forest stands. They also constitute the pool of candidate trees that are responsible for recruitment of merchantable sized trees. However, sapling information is often absent from regular inventory measurements (e.g. where no saplings are measured) even though they could be used as inputs in predicting recruitment in merchantable trees.

Aims

The objectives were to develop models describing density and diameter distribution of saplings from operational inventories, e.g. having only merchantable tree inventory, and to evaluate how stand type (naturally regenerated stands and plantations) and various thinning treatments influence these distributions.

Methods

Using data from both white spruce (Picea glauca [Moench] Voss) plantations and naturally regenerated balsam fir (Abies balsamea (L.) Mill.) stands having been commercially thinned, a zero-inflated poisson regression was used to model the stand density and a gamma regression to predict the two parameters of the Weibull used to characterize the diameter distribution.

Results

Despite the fact that the operational inventory data is often limited (e.g. species, dbh, height), the accuracy of the models was good and unbiased. It could be integrated into growth models to complete missing sapling distributions and help to correct some recruitment bias. Important differences in species composition, density and diameter distribution were found between plantations and naturally regenerated stands, but only a moderate response in diameter distribution was observed with thinning treatments.

Conclusion

These models will enable managers to estimate saplings in intermediate aged softwood forests of eastern Quebec using harvesting inventories or National Forest Inventory. Characterization of differences between plantations and naturally regenerated stands will be useful for integrating intensive plantation silviculture with ecosystem-based management. Longer-term follow-up would be needed to better evaluate the effects of thinning treatments.

1 Introduction

Forests are complex ecosystems that require accurate information in order for forest managers to make the best decisions possible (Pretzsch and Forrester 2017) and evaluate the efficiency of various management approaches. Growth and yield models have become the corner stone of most of the decision support systems (Makela et al. 2000). Furthermore, these models also help improve understanding of complex stands characterized by multi-aged and mixed species (Weiskittel et al. 2011; Pretzsch and Forrester 2017). They are also essential to predict whether evolutionary trajectories following the application of silvicultural treatments will make it possible to achieve specific management objectives. To have accurate computer-based systems that are capable of simulating dynamics at the stand level, at least three models are needed: growth of existing trees, mortality and recruitment of new trees (Fortin and DeBlois 2007). However, prior to simulating stand development, some models require initialization processes in order to feed all the input variables into the model (Mäkelä and Mäkinen 2003). Model initiation can be used to fill in information either at the tree or stand level. At the stand level, forest structure generators are used to either spatialize inventory data (i.e. attribute XY coordinates to trees) (Pretzsch 1997) or complete the diameter distribution (i.e. information on the size of all the trees is needed whereas only stand level information is available or those of trees of a certain size are measured).

Diameter distributions are important to study both forest dynamics (West et al. 2009) and to evaluate stand level wood properties or value (Weiskittel et al. 2011). Diameter distribution models have also been used to scale down the stand-level information to the tree level (Newton et al. 2005). Various probability density functions have been used to quantify tree diameter distributions (log-normal (Bliss and Reinker 1964), Gamma (Nelson 1964), Beta (Clutter and Bennett 1965), Weibull (Bailey and Dell 1973), logit-logistic (Wang and Rennolls 2005), Johnson SB (Rennolls and Wang 2005)). But they are often built to simulate merchantable sized trees (i.e. trees above a certain diameter). This stems from the fact that non-merchantable trees are generally not measured, or if so, only in small areas in forest surveys, despite the fact that they play an important ecological role (Hjeljord et al. 1990). Recruitment in most of growth and yield models predicts the number and size of the new merchantable sized trees that reach the merchantable size threshold during a simulation step. However, these models only rely on the tally of merchantable sized trees, and thus are blind to either stands with either an abundance or absence of non-merchantable sized trees (Fortin and DeBlois 2007). This lack of information can lead to a recruitment bias, whereby the simulator is not able to simulate changes in the stand notably of composition after silvicultural treatments or natural disturbances.

However, the consequences of these disturbances on future forest dynamics have to be understood (Pommerening 2006; Getzin et al. 2008). Forest managers also need this information to properly assess the differences between natural and managed stands in order to establish ecosystem-based management guidelines. Indeed, several jurisdictions have now implemented close-to-nature forestry. For example, the province of Quebec, in Eastern Canada, requires that all management must follow ecosystem-based principles since 2013 (Gouvernement du Québec 2015). The goal is to strike a balance between wood harvesting and maintenance of ecosystem structure and functions (Bergeron et al. 1999; Gauthier et al. 2009). This approach is based on minimizing the differences between natural and managed landscapes to ensure long-term maintenance of ecosystem functions. This requires the characterization of key composition (abundance or absence of a species, relative proportion), structural attributes (age or diameter structure, trees spatial distribution) and functional attributes of the forested ecosystems in order to choose the correct way to manage forest stands (Eriksson and Hammer 2006). This assessment of naturalness and respecting sustainable management is also a major concern for forest certifications such as the one proposed by the Forest Stewardship Council (FSC 2006). Gap dynamics is one of the main processes by which natural disturbances (e.g. insect epidemics, fires, windfalls, etc.) modify the canopy structure (Harvey et al. 2002; Boucher et al. 2009). Forests under the influence of gap dynamics have heterogeneous and irregular structures (Schütz 2002; Ruel et al. 2007). Interventions imitating natural forest dynamics must thus be used in order to maintain or reintroduce the heterogeneity and irregularity of the forest structure (Schütz 2002).

In the Bas-Saint-Laurent region of the province of Quebec, Canada, balsam fir (Abies balsamea (L.) Mill.) and hardwoods dominate the current extensively managed forests. The forest landscape is also covered by 15% of plantations. However, studies have shown that regular harvests and logging have led to important changes in the composition of the forests of Eastern Canada (Boucher et al. 2009) and current unmanaged forests cannot always be used as the reference of a natural state of the forests due to this long history of management (Boucher et al. 2006). Dupuis et al. (2011) showed that there has been a rarefaction of certain species such as northern white-cedar (Thuja occidentalis L.), red spruce (Picea rubens Sarg.), white spruce (Picea glauca [Moench] Voss), white pine (Pinus strobus L.) and red pine (Pinus resinosa Aiton), which is a biodiversity issue (Grondin and Cimon 2003).

In order to meet the requirements of the ecosystem-based forest management, a portion of the actual even-aged stands will be shifted towards an irregular/uneven-aged forest management path. A framework of thinning approaches with various intensities was tested, and in addition to a wood production objective, treatments can also help in restoring or maintaining some of the declining species such as white spruce (Grondin et al. 2003). Three thinning treatments were used in this study: the thinning from below (TrT_1/3) where the smallest trees are removed while ensuring equal spacing between the remaining trees, the tree release (TrT_Elite) where either 50 or 100 dominant trees per ha were released from competition on all sides and the balsam fir priority selection (TrT_BF) where all the balsam fir are harvested.

The aim of this study was to develop a unique flexible model describing the diameter distribution of saplings for intermediate aged softwood stands, for white spruce, balsam fir and hardwood species grouped together. A first objective was to generate a diameter distribution of saplings for operational inventories that only contain merchantable sized trees in order to complete the diameter distribution. The second objectives were to evaluate and to compare stand structural complexity of naturally regenerated stands and of plantations. And a third objective was to assess if the structure of both stand types varied with four various commercial thinnings (control, TrT_1/3, TrT_Elite and TrT_BF). We hypothesized that commercial thinning should modify the diameter distribution of the stands, and therefore increase their heterogeneity.

2 Materials and methods

2.1 Study site

The naturally regenerated softwood stands and plantations (stand types) used in this study are located in the eastern balsam fir-yellow birch bioclimatic subdomain of the boreal mixed-wood zone (Robitaille and Saucier 1998) in eastern Quebec, Canada. In this area, forests are characterized by mixed stands of yellow birch (Betula alleghaniensis Britton), balsam fir, white spruce and northern white-cedar (Grondin and Drouin 1998). Some other hardwoods can be observed locally (trembling aspen (Populus tremuloides Michx.), paper birch (Betula papyrifera Marsh.), red maple (Acer rubrum L.) and sugar maple (Acer saccharum Marsh.)). A total of 98 sample plots were used to calibrate the diameter distribution model. They are part of two commercial thinning experiments.

The first experiment included 72 plots of which 40 were in two operational white spruce plantations and 32 in two naturally regenerated balsam fir stands. The plot area was 450 m² (15 × 30 m). The plantations were mainly composed of planted and naturally regenerated white spruce (78% of the stand basal area, Table 1) with some balsam fir (21%, Table 1) and hardwood (1%, Table 1) natural regeneration. Balsam fir was predominant (91% of the stand basal area, Table 2) in naturally regenerated stands with some white spruce (7%, Table 2) and hardwood (2%, Table 2). Each plot had one of the three commercial thinning treatments (TrT) that were randomly attributed: control (TrT₀, 8 natural regenerated plots, 10 plantation plots), thinning from below (TrT_1/3, 8 natural regenerated plots, 10 plantation plots) and tree release (TrT_Elite, 16 natural regenerated plots and 20 plantation plots).

Table 1 Mean characteristics of plantations after treatment and for all commercial species, white spruce, balsam fir and hardwoods (standard error)

Table 2 Mean characteristics of natural stands after treatment and for all commercial species, white spruce, balsam fir and hardwoods (standard error)

The second experiment included 26 plots in two operational white spruce plantations. The plot area was 900 m² (15 × 60 m). We observed disparities between plots; however, the plots were mainly composed of planted and naturally regenerated white spruce (89% of the stand basal area, Table 1) with some balsam fir (6%, Table 1) and hardwood (5%, Table 1) natural regeneration. Each plot had one of the three commercial thinning treatments (TrT) that was randomly attributed: control (TrT₀, 5 plots), thinning from below (TrT_1/3, 9 plots) and tree release (TrT_Elite, 3 plots). An additional fourth treatment was attributed to this experiment: balsam fir priority selection (TrT_BF, 9 plots).

2.2 Data acquisition

In the province of Quebec, Canada, the merchantable diameter at breast height (dbh) limits is set to 9 cm for conifers. Trees with dbh of at least 9.1 cm are considered merchantable, whereas the others are non-merchantable in size (saplings). In most inventories, sapling characteristics are not measured. In these experiments, within each permanent sample plot, all trees with dbh greater than 1 cm were tallied by species and their dbh were measured in millimetres with a measuring tape immediately after treatment. Trees were grouped by species into three groups (Gsp): white spruce (WS), balsam fir (BF) and hardwoods (HW, white birch, yellow birch, sugar maple, red maple, ash and trembling aspen). Some rare species (northern white-cedar) and non-commercial hardwoods (i.e. elderberry (Sambucus canadensis L.), rowan (Sorbus americana Marsh.), speckled alder (Alnus rugosa (L.) Moench)) were not considered in the analysis.

The inventory data was used to obtain saplings stand-level characteristics: density of saplings (N_Sap, stems/ha) and diameter distributions of each plot as well as merchantable stand-level characteristics: density of trees (N_Mer, stems/ha), and mean diameter (D_mean, in cm), diameter variance (D_var, in cm²) and diameter standard deviation (D_sd, in cm). The mean plantation characteristics are presented in Table 1 (natural stand characteristics, Table 2). Table 3 presents the abbreviations of the variables used in this study.

Table 3 Definition and abbreviation of the variables

2.3 Model development

A two-step approach was used to predict the number and size of the saplings. The first step consists of predicting the density of sapling trees and the second in describing the diameter distribution by predicting the parameters of the distribution. In both cases, a list of potential stand-level explanatory variables including silvicultural treatment, stand type (natural or plantation), merchantable tree information and interaction between the selected variables was tested. Variables that were highly correlated were excluded (variance inflation factor (VIF) > 10, O’brien 2007). These exclusions varied by Gsp.

The GLM methods in the stats library of the R statistical programming environment (R Core Team 2017) were used to fit the models. A backwards elimination based on the corrected form of Akaike’s information criterion (AICc, Burnham and Anderson 2002) was used to select variables.

2.3.1 Density of saplings

N_Sap was modelled using a zero-inflated poison regression (Lambert 1992; Fortin and DeBlois 2007) which models together excess of zeros and count data (Eq. 1). This approach accounts for the absence of saplings observed for some Gsp in several plots, resulting in an excess of zeros, and for the density of saplings per hectare (N_Sap, Faraway 2016) (Fig. 1). It is considered that saplings are present as soon as we observed a minimum of one tree.

$$ {y}_i\sim \left\{\begin{array}{ll}0& \kern1em \mathrm{with}\ \mathrm{probability}\ \mathrm{p}\ \\ {} Poisson\left(\lambda \right)& \kern1em \mathrm{with}\ \mathrm{probability}\ 1-\mathrm{p}\ \end{array}\right. $$

(1)

where the parameters p and λ can be modelled as functions of the stand-level explanatory variables.The probability of p was modelled as a binomial regression with a logit-link function (Eq. 2).

$$ p=\frac{e^{\mathbf{Z}\gamma }}{\left(1+{e}^{\mathbf{Z}\gamma}\right)} $$

(2)

Fig. 1

Link between density of merchantable trees and density of saplings. Black dots are the measured data. White dots are the prediction of Eqs. 2 and 3. The lines represent a locally weighted regression with the predicted values. The shaded areas represent the standard errors

And the density component as a Poisson regression (Eq. 3).

$$ {N}_{Sap}={e}^{\mathbf{X}\omega } $$

(3)

where Z and X are the design matrices of covariates and ω, γ the fixed-effects parameter vectors to be estimated.

2.3.2 Diameter distribution of saplings

The complete distribution of diameter classes varies greatly among plots. However, it is observed that these distributions seemed to be very often composed of at least two modes (i.e. of different cohorts: saplings and merchantable trees). In some cases, segmentation between modes is hardly visible (e.g. WS and HW in Fig. 2a, b) but can be very important in other cases (e.g. BF in Fig. 2c, d). In addition, in Fig. 2b, only a few HW merchantable trees are present whereas the sapling distribution is more important. These examples illustrate that the sapling diameter distribution cannot be easily predicted only by extending the distribution of the merchantable trees and it should be analyzed separately.

Fig. 2

Distribution of trees per diameter classes and per Gsp from four stands illustrating various diameter structure and Gsp density. a,b Plantations, where the presence of hardwoods was important in the sapling section in a and with few saplings in b. c, d Natural stands where BF was dominant

Foremost, the density of saplings in each 1-cm dbh classes was measured and divided by the total density of saplings per Gsp and per plot. The obtained proportion (PDia_Gsp) was then characterized by a two-parameter Weibull distribution (Weibull 1951; Cao 2004).

$$ f\left(x;\alpha, \beta \right)=\frac{\alpha }{\beta }.{\left(\frac{x}{\beta}\right)}^{\left(\alpha -1\right)}.\mathit{\exp}{\left(\frac{-x}{\beta}\right)}^{\left(\alpha \right)} $$

(4)

where χ is the 1-cm dbh class, α the shape parameter and β the scale parameter.

The optimization function Optim in the R statistical programming environment (R Core Team 2017) was used to search the best combination of α_Gsp and β_Gsp for each plot by minimizing the residual sum of squares (Nelder and Mead 1965). In some plots, saplings were absent for a given Gsp or had most of their saplings concentrated in only very few dbh classes (fewer than three classes). Equation 4 converged on 79, 64 and 90 plots for WS, HW and BF, respectively.

Then, α_Gsp and β_Gsp were both modelled as a function of the available inventory variables and variables predicted with the model for the density of saplings. As they are continuous and positive variables, we used a gamma regression model for the two variables (Faraway 2016) (Eq. 5).

$$ y=\mathit{\exp}\left(\mathbf{A}\theta \right) $$

(5)

where y is either α_Gsp or β_Gsp, A is the design matrix of covariates and θ the fixed-effects parameter vectors to be estimated.

2.3.3 Model validation

The predictive accuracy of the models (i.e. the predictions of the variables N_Sap, of α_Gsp and β_Gsp) was calculated by using a repeated 10-fold cross validation (Rodriguez et al. 2009). The dataset was randomly split into 10 subsets of equal size. The first nine subsets were used to calibrate the model, and the tenth to validate it. The process was repeated till all subsets were used to validate the model. The root mean square error (RMSE) and the R-square (R²) were then calculated using the validation subset. The random segregation of the plots into the subsets was repeated 50 times. The validation result is the average of all the repetitions.

The final model was calculated by predicting the density of saplings per diameter classes. The predicted values of α_Gsp and β_Gsp were inserted in the Weibull equation and the result was multiplied by the predicted value of NSap_Gsp. The predictive accuracy of the final model was evaluated by comparing the predicted data to the observed data. RMSE and R² were calculated. In addition, a Kolmogorov-Smirnov test (Stephens 1974) was done to evaluate the accuracy of the obtained distributions. The test will return two values (D-statistic and p value). The closer the D-statistic is to 0, the more likely it is that the two samples were drawn from the same distribution. When the p value is inferior to 0.05, the null hypothesis that the two samples were drawn from the same distribution is rejected.

3 Results

3.1 Density of saplings

The density of saplings varied between stands (Fig. 1). Considering only the control plots, 36% of the trees were saplings in natural stands (Table 2), and they constituted 30% of the total density in plantations (Table 1). Furthermore, the proportion of trees between 5 and 9.1 cm in dbh (i.e. trees that should contribute fairly quickly to recruitment in the merchantable sized trees) accounted for 19% of the trees in the natural stands and 15% in the plantations. Thinning treatments tended to increase the proportion of saplings in natural stands whereas this proportion was quite constant in plantations. We noticed that the proportion per diameter class varied with Gsp. The proportion of HW does not vary between the natural stands and plantations. However, WS saplings were quite rare in natural stands and represent more than 1/3 of trees in the plantations. BF saplings were present in both stand types, but their proportion is larger in natural stands, especially in the smaller diameter classes.

Considering control plots, when merchantable trees were present in natural stands, saplings were absent for 60% of the plots for WS, 0% for BF and 12.5% for HW (respectively, 0% for WS, 6.7% for BF and 13.3% for HW in plantations). On the opposite, when saplings were present in natural stands, merchantable trees were absent for 25% of the plots for WS, 0% for BF and 25% for HW (respectively, 0% for WS, 0% for BF and 20% for HW in plantations). Those results indicated that the presence/absence of saplings cannot be only directly deducted from the presence/absence of the merchantable trees of the same species.

Among all the possible covariates tested, the binomial part of Eq. 2 with the lowest AICc was found to be (Eq. 6):

$$ {\displaystyle \begin{array}{rr}{\mathbf{Z}\boldsymbol{\upgamma}}_{Gsp}=& {a}_0+{a}_1\times StandType+{a}_3\times Nme{r}_{Gsp}+{a}_5\times TrT+\\ {}& {a}_6\times Dmea{n}_{Gsp}+{a}_9\times Nme{r}_{Gsp}\times Dmea{n}_{Gsp}\end{array}} $$

(6)

where a₀–a₉ are the fixed effect parameter estimates found in Table 4.

Table 4 Parameter estimation (standard errors) of the lowest AIC selected models (Eqs. 2 and 3) (.P < 0.1; *P < 0.05; **P < 0.001; ***P < 0.0001)

As the stand type parameter was positive (a₁ = 2.794, Table 4), the absence of saplings is more likely to occur in plantations than in naturally regenerated stands. TrT_Elite also reduced the probability of sapling occurrence (a₅ = 2.657, Table 4), when compared to control plots. The probability of sapling occurrence also decreased with the density of merchantable trees in interaction with the mean diameter stand for each Gsp (a₉ = − 0.002, Table 4). However, no differences were observed between Gsp.

The count component (Eq. 3) final form is given in Eq. 7.

$$ \kern1.25em {\displaystyle \begin{array}{ll} NSa{p}_{Gsp}=\mathit{\exp}\Big(& {a}_0+{a}_1\times StandType+{a}_2\times Nme{\mathrm{r}}_{tot}+{a}_3\times Nme{r}_{Gsp}+{a}_4\times Gsp+\\ {}& {a}_5\times TrT+{a}_6\times Dmea{n}_{Gsp}+{a}_7\times Ds{d}_{tot}+{a}_8\times Dva{r}_{tot}+\\ {}& {a}_{10}\times Gsp\times StandType+{a}_{11}\times Nme{r}_{Gsp}\times StandType+\\ {}& {a}_{12}\times Dmea{n}_{Gsp}\times StandType+{a}_{13}\times Nme{r}_{tot}\times StandType+\\ {}& {a}_{14}\times Nme{r}_{tot}\times Nme{r}_{Gsp}+{a}_{15}\times Nme{r}_{Gsp}\times Gsp+\\ {}& {a}_{16}\times Nme{r}_{tot}\times Gsp+{a}_{17}\times Gsp\times TrT+{a}_{18}\times Dmea{n}_{Gsp}\times Gsp\Big)\end{array}} $$

(7)

where a₀–a₁₈ are the fixed effect parameter estimates found in Table 4.

With WS in control stands as a reference, we observed that the density of saplings was more important for BF (a₄ = 5.035, Table 4) and similar for HW (a₄ = − 0.096, Table 4). Also, we observed that the density of WS saplings was higher in the plantations (a₁ = 1.827, Table 4) what is also the case for HW (a₁ = 1.827 and a₁₀ = − 0.477, Table 4). The sapling density was also affected by the merchantable density of all species and of the same Gsp as well as by with the mean diameter of commercially sized trees of the same Gsp (Table 4) in plantations. On the other hand, the density of BF saplings was lower in the plantations (a₁ = 1.827, a₄ = 5.035 and a₁₀ = − 3.659, Table 4).

Commercial thinning reduces sapling density of WS when compared to the control. This effect was slightly negative for TrT_BF (a₅ = − 0.177, Table 4). It was stronger for both TrT_Elite (a₅ = − 0.459, Table 4) and TrT_1/3 (a₅ = − 0.415, Table 4). However, the effect varied between Gsp, with the most important reduction in BF (a₅ = − 0.018 and a₁₇ = − 1.057, Table 4).

Model evaluation was carried out using the combination of Eqs. 6 and 7. The RMSE was 262 and R² 0.63. When validated by Gsp, the RMSE was, respectively, 98, 314 and 275 and the R², respectively, 0.63, 0.68 and 0.50 for WS, BF and HW. Finally, the RMSE was found to be 173 and 353 and theR² 0.52 and 0.66 for the plantations and naturally regenerated stands, respectively.

3.2 Diameter distribution of saplings

The observed mean value of α_Gsp and β_Gsp by stand type and by silvicultural treatment are presented in Tables 1 and 2. The diameter distribution of HW saplings was similar between plantations and natural stands and between treatments. For WS, the mean diameter of saplings was bigger, and the distribution was more extended in plantations whereas the mean density and the mean diameter were smaller in natural stands. For BF, the mean density of saplings was much more important in natural stands than in plantations (Fig. 3).

Fig. 3

Mean characteristics of the sapling diameter distribution by stand type and per silvicultural treatment. The mean values observed on the sapling density and on the two Weibull parameters (α and β) are detailed in Tables 1 and 2

α_Gsp and β_Gsp are both modelled by a gamma distribution. As the final model form was similar for the two parameters, results are presented in Eq. 8 where y is the general notation for α_Gsp and β_Gsp.

$$ {\displaystyle \begin{array}{ll}y=\mathit{\exp}\Big(& {b}_0+{b}_1\times StandType+{b}_2\times TrT+{b}_3\times {\upalpha}_{Gsp}+{b}_4\times Nme{r}_{tot}+\\ {}& {b}_5\times Nme{r}_{Gsp}+{b}_6\times NSa{p}_{WS}+{b}_7\times NSa{p}_{BF}+{b}_8\times NSa{p}_{HW}+\\ {}& {b}_9\times Dva{r}_{HW}+{b}_{10}\times Ds{d}_{WS}+{b}_{11}\times Dmea{n}_{HW}+{b}_{12}\times StandType\times {\upalpha}_{WS}+\\ {}& {b}_{13}\times StandType\times NSa{p}_{WS}+{b}_{14}\times TrT\times NSa{p}_{WS}+{b}_{15}\times NSa{p}_{WS}\times NSa{p}_{BF}+\\ {}& {b}_{16}\times NSa{p}_{HW}\times {\upalpha}_{HW}\Big)\end{array}} $$

(8)

b₀–b₁₆ are the fixed effect parameter estimates. Their values are presented in Table 5 for both α_Gsp and β_Gsp. β_Gsp used α_Gsp as a predictive variable; therefore, α_Gsp was modelled first. We observed large variations in the significant variables affecting the three Gsp (Table 5).

Table 5 Parameter estimation (standard errors) of the lowest AIC selected models for diameter distribution of sapling (Eq. 4) (.P < 0.1; *P < 0.05; **P < 0.001; ***P < 0.0001)

Stand type had a strong influence on α_Gsp and β_Gsp for WS and BF but not for HW. A positive effect on β_Gsp was found in the two Gsp (b₁ = 1.991 for WS, 0.527 for BF, Table 5). This leads to an increase in the larger diameter classes. As the value was higher for WS, there were slightly more trees in the larger diameter classes on plantations. However, the stand type can affect directly the density of saplings which is a variable of the diameter distribution model. To illustrate the global stand type effect on diameter distribution, the two stand types were tested in density and diameter distribution models while keeping constant all the other stand characteristics (Fig. 4). The distribution shape was slightly left-shifted for WS. For HW, the distribution shape was quite similar in both stand types, only the density of saplings varied. The most important effect was on BF where the distribution was left-shifted producing many more small-diameter recruits in natural stands.

Fig. 4

Simulation illustrating the effect of the stand type (natural stand or plantation) on the sapling distribution predictions for the three Gsp. All of the other stand characteristics remain identical, only the stand type changes. The density of merchantable trees per hectare and their mean dbh in cm are set to 1000 and 16 for WS, 1000 and 16 for BF, 150 and 14 for HW

No significant effect of commercial thinning was found on WS and HW. On BF, we observed a positive effect on α for the treatments TrT_1/3 and TrT_BF (b₂ = 0.646 and 2.818, respectively, Table 5). The sapling density of WS slightly reduced the effect for TrT_BF (b₁₄ = − 0.007, Table 5). The effect on β_BF was significantly negative for two treatments TrT_Elite and TrT_BF (b₂ = − 0.253 and − 0.345, respectively, Table 5). The combination of those effects results in a flattening of the distribution and a shift to smaller diameter. However, even in cases where the treatment effect was not significant, the distribution was influenced by the sapling density. As this variable was itself influenced by the treatment effect, it was difficult to interpret each variable individually. This effect of the commercial thinning on the sapling diameter distribution is illustrated in Fig. 5 while keeping constant all the other stand characteristics. The main effect of thinning treatments was on the predicted density of saplings (Eqs. 6 and 7) and was very low on the diameter structure itself. Despite the absence of a significant effect of the treatment variable, an important effect was observed in natural stands where the diameter distribution of WS was strongly left-shifted in thinned stands. The effect of the sapling density variable that was itself affected by the type of treatment appeared to be the main effect.

Fig. 5

Simulation illustrating the effect of the thinning treatment into the sapling distribution predictions for the three Gsp on the two stand types (natural stand and plantation). For each stand type, all of the stand characteristics remain identical, only the treatment changes. On the plantation (natural stand), the density of merchantable trees per hectare and their mean dbh in cm are set to 1460 (62) and 16 (16) for WS, 244 (1661) and 16 (17) for BF and 36 (38) and 14 (12) for WH

The value of α_Gsp and β_Gsp was mostly influenced by the characteristics of the studied Gsp. However, we noticed that the different Gsp can interact. The density of BF saplings was negatively related to α_WS (b₇ = − 0.0004, Table 5).

Equation 8 predicts the diameter distribution parameters for each Gsp separately. For α predictions, we obtained R² values of 0.53, 0.60 and 0.52 and RMSE of 3.5, 1.1 and 2.5, respectively, for WS, BF and HW. For β predictions, we obtained R² values of 0.42, 0.34 and 0.28 and RMSE of 1.4, 3.8 and 1.4, respectively, for WS, BF and HW.

Finally, the completely predicted diameter distribution was evaluated. The model predictions (presence/absence of saplings, density and diameter distribution) were used to generate the distribution. Firstly, the total density of saplings (Eqs. 6 and 7) was predicted, then used as an independent variable to obtain the predicted values of α_Gsp then β_Gsp (Eq. 8). α_Gsp and β_Gsp were then integrated into the Weibull function to obtain the relative proportion of saplings by diameter class. This relative proportion was multiplied by the total density of saplings. We compared the predictions to the observed density of saplings per diameter classes. On average, the results were unbiased and with R² of 0.72, 0.81 and 0.54, respectively, for WS, BF and HW and RMSE of 13.7, 42.1 and 27.3. If we analyzed the results per stand type, we observed that the model was good in plantations for the three Gsp. However, in natural stands, the results were poorer for WS and HW, considered as companion species. The Kolmogorov-Smirnov test was also performed given a D-value (and a p value) equal to 0.051 (0.332), 0.0375 (0.628) and 0.095 (0.013), respectively, for WS, BF and HW. This indicated that for WS and BF, the predicted and the observed values were drawn from the same distribution. For HW, the low D-value indicates that the distribution was quite similar but cannot be considered as drawn from the same distribution. As shown in Table 6, we observed result disparities by stand type and by treatment, but the distribution was good for WS and BF. For HW, we observed that in plantations and with the treatments, TrT_1/3 and TrT_BF predictions were less effective.

Table 6 Evaluation of the goodness-of-fit of the complete distribution for each Gsp, by stand type and per treatment

4 Discussion

Saplings contribute to the structure as well as for the renewal of the stand. To assess the impacts of management decisions, especially in close-to-nature systems, a better understanding of saplings is needed. In this study, our objective was to model the complete tree diameter distribution of a stand using operational inventory data containing only information from merchantable sized trees. We thus developed a two-step approach that first predicts the density of saplings and their diameter distribution. Despite the fact that the operational inventory data is often limited (species, dbh, height), the accuracy of the model was good, and the final model was unbiased. This approach will enable managers to estimate saplings in intermediate aged softwood forests of eastern Quebec, as it can be applied to either pre- or post-harvesting inventories or to the National Forest Inventory.

The choice to group all the hardwood species together was motivated by the fact that many species were poorly represented or not present in all the sampled plots. This grouping allowed us to understand the general structure of hardwood saplings, where balsam fir and white spruce are dominant. However, these hardwood species (yellow birch, trembling aspen, paper birch, red maple, sugar maple) can present a diversity of traits (shade tolerance, regeneration, growth, competitiveness) and commercial interest. In order to refine the specific structure of each of these species, in-depth analysis would be needed as well as larger sampling.

The two-parameter Weibull distribution was used to model the complete diameter distribution for each species group. The model fits the data well and was able to describe the various diameter distributions that were observed, especially for WS and BF (Cao 2004). For HW, the prediction of the sapling distribution was good in natural stands; however, the Kolmogorov-Smirnov test showed that in plantations and for TrT_1/3 and TrT_BF, predictions were less accurate. This case seemed to occur mainly in plots where the thinning treatment targeted saplings big enough to be in competition for resources with merchantable trees (TrT_1/3 and TrT_BF) limiting the impact on the smaller trees. The lowest amount of data available for these treatments may also explain these results.

As the models were built from data sampled in naturally regenerated stands and in plantations, the comparison of how species composition and stand diameter structure are affected by anthropic disturbances can be undertaken. The observed species composition and density were quite different between the two stand types (i.e. plantation vs natural regeneration). Models showed that the probability of having no saplings was higher in plantations and affected the three species groups in the same way. These results were expected. Indeed, microclimate, site variations (Man and Lieffers 1999) and canopy opening (Yamamoto 2000) have positive effects on the establishment and survival rate of seedlings. These conditions are more susceptible to occur in natural stands where the stand is more heterogeneous than in plantations.

When saplings were present, the density of merchantable trees was the most important parameter to describe their density. The effect of the stand type was different on the three species groups (Fig. 1). The density of merchantable HW trees was quite similar in both stand types, and if the density of the saplings slightly increased for HW in plantations, the difference was low. On the contrary, the effect was strong and opposite on BF and WS. We observed for BF a lower density of merchantable trees in plantations than in natural stands. The effect is reversed for WS where the density of merchantable trees is lower in natural stands. These differences affected saplings in the same way. In plantations, WS sapling density was higher and BF lower, whereas the opposite was observed in natural stands (Fig. 3).

Stand types also affected the diameter distribution that was characterized by the two parameters of the Weibull distribution (i.e. α and β). The two parameters were simultaneously influenced by stand type for WS and BF. In plantations, α is negatively influenced and β positively, which results from a more left-skewed and platykurtic. As a result, non-merchantable trees are more dispersed throughout the diameter classes. On the other hand, only β is influenced by the stand type for BF and no distribution parameters changed with the stand type for HW. In plantations, the distribution is shifted to the right for BF. These differences could be explained by the fact that in natural stands, the density of BF regeneration and lesser number of WS seed trees might have limited the possibilities for WS to regenerate, whereas in plantations, with site preparation, the regeneration of all Gsp start at the same time allowing the WS regeneration to be more important.

These plantations are to a certain degree an artificialization of natural forests. They contribute to a wood production objective (Grondin et al. 2003). However, in the Bas-Saint-Laurent region (Quebec, Canada), plantations are not monocultures and should not be considered completely different from current unmanaged stands. Indeed, up to 25% of the trees are not WS, and if the proportion of HW is similar in both stand types, we have a reversal of the proportion of WS and BF. Dupuis et al. (2011) showed a rarefaction of WS in current unmanaged forests which are now dominated by BF and HW. These results could therefore be used to integrate plantation silviculture into ecosystem-based management.

In both stand types, the sampled plots were also part of thinning treatment experiments established 5 years ago. These treatments aimed at increasing tree growth. TrT_Elite also tried to increase stand structural diversity. No interaction between treatment and stand type was found in the model; the thinning seemed to affect similarly plantations and naturally regenerated stands. The diameter distribution model showed that the thinning treatments (TrT_1/3 and TrT_BF) increased the probability of observing saplings. Our results suggested slight impact on the tree diameter distribution too. No direct impact of the thinning on diameter distribution was observed; the proportion of saplings was slightly modified. Thinning significantly affected the stand density and had different impacts on the three species groups. For WS, thinning reduced the density of saplings. The strongest effect was with TrT_Elite and TrT_1/3. It was expected because these treatments removed the smallest trees or the competition around elite trees without targeting a specific species, unlike TrT_BF which specifically targeted BF. This explains why the strongest negative effect affecting BF was with TrT_BF. In an objective of restoring WS and limiting the presence of BF in forests (Grondin et al. 2003), TrT_BF could be the most appropriate. More unexpectedly, we observed that the other thinning treatments slightly increase BF and HW density. As the selection of the elite trees or the selection of the competition to be removed in TrT_1/3 was distributed across all the stand, the spatial structure of these species could be a possible explanation. However, more data will be necessary to evaluate this spatial effect.

The limited amount of data for some treatments may also explain the weak response. However, as the data were only taken a short period after thinning (i.e. 5 years), it is also possible that the effects are not visible yet. The modification of stand density could affect the competition for resources by the saplings, but it is difficult to differentiate the true response to thinning from the effect of thinning (i.e. the chainsaw effect). It is therefore difficult to determine if the small differences observed in stand structure compared to the control stands will be maintained, or if the differences will increase or decrease in time. However, Schneider et al. (2013) analyzed the effects of selective pre-commercial thinning on growth yield and diameter distribution over a 28-year period and showed no significant effects. Gauthier et al. (2015) showed after a TrT_1/3, the structural heterogeneity increased 10 years after the harvest, but this was mainly due to the growth of larger trees and no effect on saplings was described. These stands should be monitored during many years to evaluate the real impact of these anthropic disturbances over the long run (Pretzsch 2005). Indeed, according to Schütz (2002), to convert regular/even-aged stands to irregular/uneven-aged stands with TrT_Elite, many entries are necessary.

5 Conclusion

In conclusion, our model succeeded in characterizing the density and the diameter distribution of saplings. Using operational inventory data, it could be used as an input into growth models calibrated exclusively with merchantable trees to correct some recruitment bias. Important differences were found between plantations and naturally regenerated stands. These results help to improve our understanding of the dynamic of these two stand types and are important to adapt our management decisions. After 5 years, the effects of the thinning treatments are quite low and it is not yet possible to evaluate these treatments with respect to ecosystem-based management objectives. However, the small responses of the diameter distribution indicate trends. Longer-term monitoring would be required to observe the effects of these treatments on saplings.