Skip to main content

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


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.


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.


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.


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.


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.


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 (TrT1/3) where the smallest trees are removed while ensuring equal spacing between the remaining trees, the tree release (TrTElite) where either 50 or 100 dominant trees per ha were released from competition on all sides and the balsam fir priority selection (TrTBF) 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, TrT1/3, TrTElite and TrTBF). 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 m2 (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 (TrT0, 8 natural regenerated plots, 10 plantation plots), thinning from below (TrT1/3, 8 natural regenerated plots, 10 plantation plots) and tree release (TrTElite, 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 m2 (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 (TrT0, 5 plots), thinning from below (TrT1/3, 9 plots) and tree release (TrTElite, 3 plots). An additional fourth treatment was attributed to this experiment: balsam fir priority selection (TrTBF, 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 (NSap, stems/ha) and diameter distributions of each plot as well as merchantable stand-level characteristics: density of trees (NMer, stems/ha), and mean diameter (Dmean, in cm), diameter variance (Dvar, in cm2) and diameter standard deviation (Dsd, 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

NSap 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 (NSap, 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. $$

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)} $$
Fig. 1
figure 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 } $$

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
figure 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 (PDiaGsp) 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)} $$

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) $$

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 NSap, 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 (R2) 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 NSapGsp. The predictive accuracy of the final model was evaluated by comparing the predicted data to the observed data. RMSE and R2 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}} $$

where a0a9 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 (a1 = 2.794, Table 4), the absence of saplings is more likely to occur in plantations than in naturally regenerated stands. TrTElite also reduced the probability of sapling occurrence (a5 = 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 (a9 = − 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}} $$

where a0a18 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 (a4 = 5.035, Table 4) and similar for HW (a4 = − 0.096, Table 4). Also, we observed that the density of WS saplings was higher in the plantations (a1 = 1.827, Table 4) what is also the case for HW (a1 = 1.827 and a10 = − 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 (a1 = 1.827, a4 = 5.035 and a10 = − 3.659, Table 4).

Commercial thinning reduces sapling density of WS when compared to the control. This effect was slightly negative for TrTBF (a5 = − 0.177, Table 4). It was stronger for both TrTElite (a5 = − 0.459, Table 4) and TrT1/3 (a5 = − 0.415, Table 4). However, the effect varied between Gsp, with the most important reduction in BF (a5 = − 0.018 and a17 = − 1.057, Table 4).

Model evaluation was carried out using the combination of Eqs. 6 and 7. The RMSE was 262 and R2 0.63. When validated by Gsp, the RMSE was, respectively, 98, 314 and 275 and the R2, respectively, 0.63, 0.68 and 0.50 for WS, BF and HW. Finally, the RMSE was found to be 173 and 353 and theR2 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
figure 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}} $$

b0b16 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 (b1 = 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
figure 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 TrT1/3 and TrTBF (b2 = 0.646 and 2.818, respectively, Table 5). The sapling density of WS slightly reduced the effect for TrTBF (b14 = − 0.007, Table 5). The effect on βBF was significantly negative for two treatments TrTElite and TrTBF (b2 = − 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
figure 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 (b7 = − 0.0004, Table 5).

Equation 8 predicts the diameter distribution parameters for each Gsp separately. For α predictions, we obtained R2 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 R2 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 R2 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, TrT1/3 and TrTBF 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 TrT1/3 and TrTBF, 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 (TrT1/3 and TrTBF) 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. TrTElite 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 (TrT1/3 and TrTBF) 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 TrTElite and TrT1/3. It was expected because these treatments removed the smallest trees or the competition around elite trees without targeting a specific species, unlike TrTBF which specifically targeted BF. This explains why the strongest negative effect affecting BF was with TrTBF. In an objective of restoring WS and limiting the presence of BF in forests (Grondin et al. 2003), TrTBF 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 TrT1/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 TrT1/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 TrTElite, 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.

Data availability

The data that support the findings of this study are available from Robert Schneider but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are, however, available from the authors upon reasonable request and with permission of Robert Schneider.


  • Bailey RL, Dell TR (1973) Quantifying diameter distributions with the Weibull function. For Sci 19(2):97–104

    Google Scholar 

  • Bergeron Y, Harvey B, Leduc A, Gauthier S (1999) Forest management guidelines based on natural disturbance dynamics: stand-and forest-level considerations. For Chron 75(1):49–54

    Google Scholar 

  • Bliss CI, Reinker KA (1964) A lognormal approach to diameter distributions in even-aged stands. For Sci 10(3):350–360

    Google Scholar 

  • Boucher Y, Arseneault D, Sirois L (2006) Logging-induced change (1930-2002) of a preindustrial landscape at the northern range limit of northern hardwoods, eastern Canada. Can J For Res 36(2):505–517

    Google Scholar 

  • Boucher Y, Arseneault D, Sirois L, Blais L (2009) Logging pattern and landscape changes over the last century at the boreal and deciduous forest transition in eastern Canada. Landsc Ecol 24(2):171–184

    Google Scholar 

  • Burnham KP, and Anderson DR (2002) Model selection and multimodel inference: a practical information-theoretic approach. Edited bySpringer-Verlag

    Google Scholar 

  • Cao QV (2004) Predicting parameters of a weibull function for modeling diameter distribution. For Sci 50(5):682–685

    Google Scholar 

  • Clutter JL, and Bennett FA (1965) Diameter distributions in old-field slash pine plantations. Georgia Forest Research Council

  • Dupuis S, Arseneault D, Sirois L (2011) Change from pre-settlement to present-day forest composition reconstructed from early land survey records in eastern Québec, Canada. J Veg Sci 22(3):564–575

    Google Scholar 

  • Eriksson S, Hammer M (2006) The challenge of combining timber production and biodiversity conservation for long-term ecosystem functioning—a case study of Swedish boreal forestry. For Ecol Manag 237(1–3):208–217

    Google Scholar 

  • Faraway JJ (2016) Extending the linear model with R: generalized linear, mixed effects and nonparametric regression models. CRC press

  • Fortin M, DeBlois J (2007) Modeling tree recruitment with zero-inflated models: the example of hardwood stands in southern Québec, Canada. For Sci 53(4):529–539

    Google Scholar 

  • FSC, Forest Stewardship Council (2006) FSC principles and criteria for forest stewardship

  • Gauthier S, Vaillancourt M-A, Leduc A, De Grandpré L, Kneeshaw D, Morin H, Drapeau P, Bergeron Y (2009) Ecosystem management in the boreal forest. Presses de l’Université du Québec, Québec, Can

    Google Scholar 

  • Gauthier M-M, Barrette M, Tremblay S (2015) Commercial thinning to meet wood production objectives and develop structural heterogeneity: a case study in the spruce-fir forest, Quebec, Canada. Forests 6(2):510–532 Multidisciplinary Digital Publishing Institute

    Google Scholar 

  • Getzin S, Wiegand T, Wiegand K, He F (2008) Heterogeneity influences spatial patterns and demographics in forest stands. J Ecol 96(4):807–820

    Google Scholar 

  • Gouvernement du Québec (2015) Loi sur l’aménagement durable du territoire forestier. Chapitre A-18:1

    Google Scholar 

  • Grondin P, Cimon A (2003) Les enjeux de biodiversité relatifs à la composition forestière. In: Gouvernement du Québec, Ministère des Ressources naturelles. Parcs, de la Faune et des

    Google Scholar 

  • Grondin F, Drouin N (1998) Optitek sawmill simulator-user’s guide. Forintek Canada Corporation, Québec

    Google Scholar 

  • Grondin P, Noël J, and Hotte D (2003) Raréfaction de l’épinette blanche dans la sapinière de la forêt boréale. P. Grondin et A. Cimon (éd.), Les enjeux de biodiversité relatifs à la composition forestière, ministère des Ressources naturelles, de la Faune et des Parcs, Direction de la recherche forestière et Direction de l’environnement forestier

  • Harvey BD, Leduc A, Gauthier S, Bergeron Y (2002) Stand-landscape integration in natural disturbance-based management of the southern boreal forest. For Ecol Manag 155(1–3):369–385

    Google Scholar 

  • Hjeljord O, Hövik N, Pedersen HB (1990) Choice of feeding sites by moose during summer, the influence of forest structure and plant phenology. Ecography 13(4):281–292

    Google Scholar 

  • Lambert D (1992) Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics 34(1):1–14

    Google Scholar 

  • Mäkelä A, Mäkinen H (2003) Generating 3D sawlogs with a process-based growth model. For Ecol Manag 184(1–3):337–354

    Google Scholar 

  • Makela A, Landsberg J, Ek AR, Burk TE, Ter-Mikaelian M, Agren GI, Oliver CD, Puttonen P (2000) Process-based models for forest ecosystem management: current state of the art and challenges for practical implementation. Tree Physiol 20(5–6):289

    PubMed  Google Scholar 

  • Man R, Lieffers VJ (1999) Effects of shelterwood and site preparation on microclimate and establishment of white spruce seedlings in a boreal mixedwood forest. For Chron 75(5):837–844

    Google Scholar 

  • Nelder JA, Mead R (1965) A simplex method for function minimization. Comput J 7(4):308–313

    Google Scholar 

  • Nelson TC (1964) Diameter distribution and growth of loblolly pine. For Sci 10(1):105–114

    Google Scholar 

  • Newton PF, Lei Y, Zhang SY (2005) Stand-level diameter distribution yield model for black spruce plantations. For Ecol Manag 209(3):181–192

    Google Scholar 

  • O’brien RM (2007) A caution regarding rules of thumb for variance inflation factors. Qual Quant 41(5):673–690

    Google Scholar 

  • Pommerening A (2006) Evaluating structural indices by reversing forest structural analysis. For Ecol Manag 224(3):266–277

    Google Scholar 

  • Pretzsch H (1997) Analysis and modeling of spatial stand structures. Methodological considerations based on mixed beech-larch stands in lower saxony. For Ecol Manag 97(3):237–253

    Google Scholar 

  • Pretzsch H (2005) Diversity and productivity in forests: evidence from long-term experimental plots. In Forest diversity and function Springer pp:41–64

  • Pretzsch H, and Forrester DI (2017) Stand dynamics of mixed-species stands compared with monocultures. In Mixed-species forests. Springer. pp. 117–209

  • R Core Team (2017) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. Available from

  • Rennolls K, Wang M (2005) A new parameterization of Johnson’s SB distribution with application to fitting forest tree diameter data. Can J For Res 35(3):575–579

    Google Scholar 

  • Robitaille A, and Saucier J.-P. (1998) Paysages régionaux du Québec méridional. Forestiers, Québec (Province). Direction de la gestion des stocks; Gouvernement du Québec, Ministère des ressources naturelles

  • Rodriguez JD, Perez A, Lozano JA (2009) Sensitivity analysis of k-fold cross validation in prediction error estimation. IEEE Trans Pattern Anal Mach Intell 32(3):569–575

    Google Scholar 

  • Ruel J-C, Roy V, Lussier J-M, Pothier D, Meek P, Fortin D (2007) Mise au point d’une sylviculture adaptée à la forêt boréale irrégulière. For Chron 83(3):367–374

    Google Scholar 

  • Schneider R, Bégin J, Danet A, Doucet R (2013) Effect of selective precommercial thinning on balsam fir stand yield and structure. For Chron 89(6):759–768

    Google Scholar 

  • Schütz J (2002) Silvicultural tools to develop irregular and diverse forest structures. Forestry 75(4):329–337

    Google Scholar 

  • Stephens MA (1974) EDF statistics for goodness of fit and some comparisons. J Am Stat Assoc 69(347):730–737

    Google Scholar 

  • Wang M, Rennolls K (2005) Tree diameter distribution modelling: introducing the logit logistic distribution. Can J For Res 35(6):1305–1313

    Google Scholar 

  • Weibull W (1951) A statistical distribution function of wide applicability. J Appl Mech 18(3):293–297

    Google Scholar 

  • Weiskittel AR, Hann DW, Kershaw Ja, Vanclay JK (2011) Forest growth and yield modeling. John Wiley & Sons, Ltd, Chichester

    Book  Google Scholar 

  • West GB, Enquist BJ, Brown JH (2009) A general quantitative theory of forest structure and dynamics. Proc Natl Acad Sci U S A 106(17):7040–7045

    CAS  PubMed  PubMed Central  Google Scholar 

  • Yamamoto S-I (2000) Forest gap dynamics and tree regeneration. J For Res 5(4):223–229

    Google Scholar 

Download references


We thank MFFP (ministère des Forêts, de la Faune et des Parc) for providing some of the data and for the time allowed to work on this manuscript.


This study was funded by Fonds québécois de la recherche sur la nature et les technologies (FQRNT).

Author information

Authors and Affiliations


Corresponding author

Correspondence to Emmanuel Duchateau.

Ethics declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Additional information

Handling Editor: Céline Méredieu

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This article is part of the Topical Collection on Mensuration and modelling for forestry in a changing environment

Contributions of the co-authors Emmanuel Duchateau: designed figures and tables in the manuscript.

Stéphane Tremblay, Robert Schneider: conceived and designed the experiments.

Emmanuel Duchateau, Laurie Dupont-Leduc: analyzed the data.

Emmanuel Duchateau, Robert Schneider, Stéphane Tremblay, Laurie Dupont-Leduc: wrote and revised the paper.

Robert Schneider: funding acquisition.

Rights and permissions

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Duchateau, E., Schneider, R., Tremblay, S. et al. Density and diameter distributions of saplings in naturally regenerated and planted coniferous stands in Québec after various approaches of commercial thinning. Annals of Forest Science 77, 38 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Stand structure
  • Diameter distribution
  • Weibull distribution
  • White spruce plantation
  • Naturally regenerated softwood stand
  • Thinning treatment