- Research Paper
- Published:

# 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*
**volume 77**, Article number: 38 (2020)

## Abstract

###
*Key Message*

*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^{2} (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_{0}, 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).

The second experiment included 26 plots in two operational white
spruce plantations. The plot area was 900 m^{2}
(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_{0}, 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^{2}) 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.

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

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

And the density component as a Poisson regression (Eq. 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.

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

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

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*^{2}) 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*^{2} 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):

where *a*_{0}–*a*_{9} are the fixed effect parameter
estimates found in Table 4.

As the stand type parameter was positive (*a*_{1} = 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*_{5} = 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*_{9} = − 0.002, Table 4). However, no differences were observed between
Gsp.

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

where *a*_{0}–*a*_{18} 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*_{4} = 5.035, Table 4) and similar for HW (*a*_{4} = − 0.096, Table 4). Also, we observed that the density of WS
saplings was higher in the plantations (*a*_{1} = 1.827, Table 4) what is also the case for HW (*a*_{1} = 1.827 and *a*_{10} = − 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} = 1.827, *a*_{4} = 5.035 and *a*_{10} = − 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*_{5} = − 0.177,
Table 4). It was stronger for both
TrT_{Elite} (*a*_{5} = − 0.459, Table 4) and TrT_{1/3} (*a*_{5} = − 0.415, Table
4). However, the effect varied
between Gsp, with the most important reduction in BF (*a*_{5} = − 0.018 and *a*_{17} = − 1.057, Table 4).

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

α_{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}.

*b*_{0}–*b*_{16} 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).

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} = 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.

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*_{2} = 0.646 and 2.818,
respectively, Table 5). The sapling
density of WS slightly reduced the effect for TrT_{BF}
(*b*_{14} = − 0.007,
Table 5). The effect on
β_{BF} was significantly negative for two treatments
TrT_{Elite} and TrT_{BF} (*b*_{2} = − 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.

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*_{7} = − 0.0004,
Table 5).

Equation 8 predicts the
diameter distribution parameters for each Gsp separately. For α predictions, we
obtained *R*^{2} 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*^{2} 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*^{2} 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.

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

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

## References

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

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

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

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

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

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

*Edited by*Springer-VerlagCao QV (2004) Predicting parameters of a weibull function for modeling diameter distribution. For Sci 50(5):682–685

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Pretzsch H (2005) Diversity and productivity in forests: evidence from long-term experimental plots.

*In*Forest diversity and function Springer pp:41–64Pretzsch H, and Forrester DI (2017) Stand dynamics of mixed-species stands compared with monocultures.

*In*Mixed-species forests. Springer. pp. 117–209R Core Team (2017) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. Available from https://www.r-project.org/

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

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

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

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

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

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

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

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

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

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

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

## Acknowledgments

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.

## Funding

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

## 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

## About this article

### 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). https://doi.org/10.1007/s13595-020-0929-5

Received:

Accepted:

Published:

DOI: https://doi.org/10.1007/s13595-020-0929-5

### Keywords

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