Post-cutting Mortality Following Experimental Silvicultural Treatments in Unmanaged Boreal Forest Stands

Montoro Girona, Miguel; Morin, Hubert; Lussier, Jean-Martin; Ruel, Jean-Claude

doi:10.3389/ffgc.2019.00004

ORIGINAL RESEARCH article

Front. For. Glob. Change, 15 March 2019

Sec. Forest Disturbance

Volume 2 - 2019 | https://doi.org/10.3389/ffgc.2019.00004

Post-cutting Mortality Following Experimental Silvicultural Treatments in Unmanaged Boreal Forest Stands

$\r\nMiguel Montoro Girona,*$ Miguel Montoro Girona^1,2^*

Hubert Morin²

Jean-Martin Lussier³

Jean-Claude Ruel⁴

¹Restoration Ecology Group, Department of Wildlife, Fish and Environmental Studies, Swedish University of Agricultural Sciences, Umeå, Sweden
²Département des Sciences Fondamentales, Université du Québec à Chicoutimi, Chicoutimi, QC, Canada
³Canadian Forest Service, Canadian Wood Fibre Centre, Natural Resources Canada, Quebec, QC, Canada
⁴Faculté de Foresterie et Géomatique, Université Laval, Quebec, QC, Canada

Partial cutting has been recommended as an alternative harvesting method to ensure the sustainable management of boreal forests. The success of this approach is closely linked to the survival of residual trees as additional losses through mortality could affect post-cutting timber production at harvest. To better quantify post-cutting mortality in previously unmanaged boreal forests, we addressed two main questions: (1) what is the level of mortality 10 years after cutting? and (2) what ecological factors are involved in this phenomenon? Even-aged black spruce [Picea mariana (Mill.) B.S.P.] stands in the Canadian boreal forest were subjected to three experimental shelterwood treatments, a seed-tree treatment and an untreated control. Tree status (live/dead) was recorded prior to cutting and 10 years after cutting. Dead trees were classified as standing dead, overturned or broken. Ten years after experimental seed-tree treatment, 60% of residual trees were dead, compared to 30% for the shelterwood cuttings. Windthrow (overturned and broken trees) represented 80% of residual tree mortality; only the amount of overturning was influenced by treatment. Broken trees were associated with small-diameter trunks, stands having high growth prior to cutting, younger stands or forest plots located near to adjacent cuts (<200 m). Overturning was associated with a high harvesting intensity and large-diameter trees. Standing dead mortality was the most difficult to explain: it was related to untreated plots having suppressed and small-diameter trees. Based on these results, applying intermediate levels of harvest intensity could reduce post-cutting damage. Understanding tree mortality after cutting is essential to reduce economic losses, improve silvicultural planning and stand selection and ensure ultimately the sustainable harvest of North American boreal forests.

Introduction

Currently, clear-cutting is used as the harvesting method in 93% of all harvested areas in the North American boreal forest (National Forest Database, 2016). However, there has been a recent paradigm change in forest management aiming to reduce the impact of harvesting practices on forest ecosystems (Messier et al., 2013; Puettmann et al., 2015). Ecosystem-based management proposes an increased use of partial cutting to reduce differences between natural and managed ecosystems (Gauthier et al., 2009). Sustainable management of the boreal forest and viable yields of timber from partial cutting require that most retained trees survive and grow rather than being lost to windthrow shortly after harvest.

In the boreal spruce-fir forest of Eastern Canada, harvesting often takes place in previously unmanaged stands, where the post-cutting mortality after partial cutting is unknown. Partial cutting can affect the risk of tree mortality as canopy openings increase wind penetration (Gardiner et al., 1997; Ruel, 2000). Studies in Europe and North America have quantified mortality in residual forests subjected to different partial cutting treatments and have reported that losses range anywhere from 15 to 74% of stems (Scott and Mitchell, 2005; Hautala and Vanha-Majamaa, 2006; Rosenvald et al., 2008; Solarik et al., 2012; Wallentin and Nilsson, 2013; Anyomi and Ruel, 2015). One to six years after cutting, mortality was assessed to be 2.5 to 4 times higher in partially cut stands than in untreated stands (Bladon et al., 2008). However, few studies have quantified the impact of pre-cutting stand characteristics on losses following partial cutting in previously unmanaged forests (Riopel et al., 2010). This lack of knowledge is particularly true for black spruce forests where partial cutting has not been commonly practiced (Solarik et al., 2012; Urgenson et al., 2013; Prévost and Dumais, 2014).

The shelterwood system is a promising partial-cutting treatment with intermediate level of harvest intensity that could allow some harvesting of timber and favor the development of advanced regeneration prior to the final felling in black spruce stands, while maintaining or developing key ecological attributes of mature stands (Vanha-Majamaa et al., 2007; Bouchard, 2009; Gauthier et al., 2009; Montoro Girona et al., 2017, 2018b). Seed-tree harvesting, on the other hand, remove most of trees from a stand, maintaining them only in small groups or as dispersed individuals to provide an even-aged forest (Nyland, 2016). The remaining trees are chosen to provide sufficient seed sources after harvesting; the remaining cover is low, thus enabling light to reach the soil surface. The shelterwood and seed-tree systems could offer potentially sustainable silvicultural practices. However, their impact on regeneration, stand growth and mortality needs to be assessed to confirm this silvicultural option for boreal forests.

Natural disturbance regimes determine the dynamics, structure and composition of forests and also control ecosystem functioning (Bergeron et al., 1999; Montoro Girona et al., 2018a; Portier et al., 2018). In forests characterized by long fire cycles, insect outbreaks and windthrow play important roles (Waldron et al., 2013; De Grandpré et al., 2018; Navarro et al., 2018). Windthrow can be difficult to predict (Klaus et al., 2011; Díaz-Yáñez et al., 2016) as it is a complex phenomenon influenced by multiple factors including climate, topography, forestry practices as well as stand and soil characteristics. Windthrow could contribute to post-cutting mortality, as many of these factors change in a partially harvested forest. Tree vulnerability to windthrow differs in terms of social status (dominant vs. suppressed trees) and spatial position, e.g., the edge effect (Rich et al., 2007; Thorpe et al., 2008; Lavoie et al., 2012; Solarik et al., 2012; Seidl et al., 2014). Tall, large-diameter dominant trees located along the edges of skid trails are more susceptible to overturning (Coates, 1997; Rich et al., 2007; Mitchell, 2012). Differences also exist between conifer species in boreal forests in terms of their susceptibility to windthrow (Ruel, 2000). At the stand level, forest structure can affect post-cutting mortality as open and older stands are more susceptible to wind (Scott and Mitchell, 2005; Prévost and Dumais, 2014). Geographical and topographical characteristics, such as forest fragmentation, climate and exposure also determine the frequency and scale of windthrow events (Achim et al., 2005; Waldron et al., 2013). The fragmentation of the forest canopy from adjacent cuttings increases the edge area, in turn influencing wind exposure (Scott and Mitchell, 2005; Harper et al., 2016). High wind speeds and topographic exposure have been associated with the “wind tunnel effect” and catastrophic windthrow after partial cutting (Waldron et al., 2013; Seidl et al., 2014; Anyomi et al., 2016).

Significant progress has been made in mechanical modeling of windthrow (Gardiner et al., 2008). The developed models consider stem breakage and tree uprooting as two different processes that must be modeled separately. Mortality in old, unmanaged stands can also occur as trees die standing (Cimon-Morin et al., 2010). According to Aakala et al. (2006), standing mortality is a major process in the development of old-growth balsam fir [Abies balsamea (L.) Mill.]–black spruce forests. Not considering this type of mortality could lead to substantial underestimations of mortality after partial cutting (Thorpe et al., 2008); thus, models trying to predict mortality must take this into account. Thorpe et al. (2008) suggested that standing dead mortality and windthrow could be compensatory mechanisms, but the importance and the processes behind this type of mortality remain to be determined.

This study examines the amount and nature of post-cutting mortality, 10 years after various experimental shelterwood and seed-tree cuttings, in unmanaged even-aged black spruce stands. Our aim is to (i) quantify post-cutting mortality across a gradient of silvicultural treatments ranging from 0 to 75% harvest intensity and (ii) determine the factors influencing the different types of mortality in the residual stands.

Materials and Methods

Study Area

The study was conducted in the boreal forest of Quebec, Canada (Figure 1). The stands were located in the Monts-Valin and North Shore regions that include two bioclimatic zones: the balsam fir–white birch (Betula papyrifera Marsh.) and the eastern spruce–feathermoss zones (Saucier et al., 1998). Climate is subhumid subpolar and mean annual temperatures range between −2 and 1.5°C with an annual precipitation of 950–1,350 mm (Robitaille and Saucier, 1998). Annual mean wind speed ranges from 3.0 to 4.3 m·s⁻¹ (Environment and Climate Change Canada, 2017). The topography is undulating and heterogeneous with altitudes ranging from 460 to 680 m. Surface deposits consist primarily of thick glacial till. Rocky outcrops are common along the top of steep slopes (Robitaille and Saucier, 1998). As fire cycles are long—between 270 to >500 years (Bouchard et al., 2008)—windthrow and spruce budworm (Choristoneura fumiferana Clem.) outbreaks are the main natural disturbances (De Grandpré et al., 2000).

FIGURE 1

Figure 1. Study area location (a) in the boreal zone of North America and (b) the location of study blocks by stand type in the Monts-Valin and North Shore regions of Quebec, Canada.

Experimental Design

To evaluate post-cutting mortality 10 years after experimental silvicultural treatments, a factorial split-plot experimental design with complete randomized blocks was set up in previously unmanaged, even-aged black spruce-dominated stands in 2003 (Figure S1 and Table 1). This forest type is typical of the regional landscape and represents 14–31% of all forests in Eastern Canada (Cyr et al., 2007; Bouchard et al., 2008). We selected two stand types based on developmental stages (age and density criterion): (i) younger and dense stands (80–100 years, 2600 trees ha⁻¹) and (ii) older and more open stands (110–160 years, 1,500 trees ha⁻¹) (Figure S2). According to Bergeron et al. (2007), the first category would consist of the first cohort of trees that establishes after fire whose structure could be maintained through clear-cutting. Partial cutting could be used to develop structures corresponding to the second-cohort phase. The second category of stands in our study relates more to the understory re-initiation stage of Oliver (1981) and the beginning of the second-cohort group of Bergeron et al. (2007). Partial cutting in this second category could be especially useful for maintaining the forest's current structure. In all cases, black spruce accounted for at least 75% of the stand basal area (Table 2). Six blocks were studied following a stratified random sampling (three of each stand type), each including five experimental units with a replicate of each silvicultural treatment and an untreated control (a total of 30 experimental units) (Figure 3a). Silvicultural treatments were randomly assigned to each experimental unit. The experimental units consisted of 3 ha permanent square plots, and were relatively homogeneous and comparable within the same block in terms of tree composition and density (Tables 1, 2). Sites with steep slopes (≥30%) were excluded due to the constraints on mechanized harvesting. The experiment was set up in operational harvesting areas, some experimental units were thus bordered by clearcuts. In such case, a 30 m buffer of untreated stand was left between clearcuts and the experimental units (Figure 3a). More information about our experimental design here: https://visitesvirtuelles.partenariat.qc.ca/monts-valin/.

TABLE 1

Table 1. Tree characteristics and species composition by silvicultural treatment for each stand type before cutting (mean ± standard error).

TABLE 2

Table 2. Stand characteristics by silvicultural treatment for each stand type before and 1 year after cutting (mean ± standard error).

Silvicultural Treatments

Four mechanized cutting treatments were applied in 2003: mini-strip shelterwood, distant selection, close selection and seed-trees (Montoro Girona et al., 2018b). The first three treatments are variants of the uniform shelterwood system (Smith et al., 1997). The studied treatments differed in terms of harvesting intensity, the spatial distribution of the skidding trails and the width of residual strips (Figure 2 and Table 3). Prescribed harvest intensity was 50 and 75% of the basal area for shelterwoods and seed-trees, respectively (Tables 2, 3). Mini-strip consisted of a succession of 5 m wide cut strips interspersed by 5 m wide residual strips. In the case of close selection and distant selection, 5 m wide trails were set at 15 and 25 m intervals, respectively. Distant selection had secondary trails transverse to the main harvesting trails, each separated by 10 m. In the case of distant selection and close selection, trees located at the residual strip were partially harvested on either side of the trails at a maximum distance of 5 m from the trail edge (Montoro Girona et al., 2017). Tree selection was performed by the harvesting operators, sequentially counting trees in groups of three and removing the largest one. Seed-trees had wider 15 m cut strips with 5 m wide intact residual strips. Harvest operations were performed with a cut-to-length system using a harvester and a forwarder. In all instances, skid trails were set in the direction of the slope.

FIGURE 2

Figure 2. Spatial patterns of studied treatments in hemispheric photographs and diagrams at the ground level. Open bands represent total harvested surface or skidding trails, and dark bands indicate the residual strips. Double arrows indicate the width of the skidding trails and residual strips.

TABLE 3

Table 3. Main characteristics of experimental treatments.

Study Variables: Measurements and Compilation

Tree and Stand Measurements

Rectangular permanent sampling plots (10 × 60 m) were established in the center of each experimental unit, perpendicular to the skidding trails. Post-cutting mortality was measured 1 year before cutting and then 1 and 10 years after cutting on trees with a diameter at breast height (DBH) ≥9 cm (minimal DBH for merchantable stems in boreal forest). Each dead merchantable stem located in the sampling plot was noted and classified into three possible types of mortality: broken, standing dead and overturned, where broken and overturned trees were considered as windthrow. The current study focuses only on black spruce, which composed 93% of the sampled trees of the first inventory.

To take tree measurements, two inventories were conducted in the 30 experimental units that made up this experiment (2 stand types × 3 blocks × 5 treatments). The first was general for all trees (n = 3,739 and 2,243, before and after cutting, respectively) and included species identification, DBH, wound severity (five levels from intact to major damage) and tree status (live/dead). The second inventory involved a subsample of randomly selected trees and included tree height and crown length (n = 201 and 143, before and after cutting, respectively).

The second inventory allowed the calibration of a prediction model of height with DBH, to make height estimations of all the tallied trees in the sample plots. The Naslund model with random slope effects was retained among six candidate models (Model S1). Based on the crown length measurements before cutting, a categorical variable was established to classify the trees in crown classes (full, reduced and leaning). The stand variables (volume, density, harvest intensity, mortality, basal area and dominant height) were estimated from the inventory data (Table 4). Both inventories were conducted 1 year before cutting (2002) and 10 years after silvicultural intervention (2013).

TABLE 4

Table 4. Grouping and description of the variables examined.

Growth before cutting (GBC) and stand age were used to provide information about the status of trees and forest structure (Thorpe et al., 2008; Montoro Girona et al., 2016). Based on dendroecological approaches, we measured GBC and tree age using tree-ring series (Baral et al., 2016; Khakimulina et al., 2016). In the summer of 2014, increment cores were extracted at breast height using an increment borer from 34 to 38 randomly selected trees per plot, for a total of 1,039 samples. Samples were air-dried, mounted on wood boards and sanded. Tree-ring width was measured using WinDendro™(version 2009, Regent Instruments, Quebec) or a manual Henson micrometer (LINTAB™, Rinntech, Heidelberg, Germany) having an accuracy of 0.01 mm. The individual tree-ring series were cross-dated using TSAP-Win™ (Rinntech, Heidelberg, Germany). GBC was established as the mean growth over 20 years for each tree. An average of GBC was estimated at the stand level for each experimental unit (n = 30). For the study of stand age, wood discs (n = 349) were collected at the root collar of trees located in square 400 m² (20 × 20 m) plots in a cutting area within each block (n = 6). Tree age was determined using a binocular microscope to count the tree rings.

Competition between black spruce was estimated 10 years after cutting for a random sub-sample of 240 trees in each experimental unit. Hegyi's competition index (CI_i) was selected as it is strongly correlated with basal area growth in black spruce stands (Mailly et al., 2003). The distance (Dist_ij) and DBH of each neighboring tree (_j) within a 4 m radius of the subject tree (_i) were measured (n = 240) to calculate the CI_i:

\begin{array}{rcl} C I_{i} = \sum_{j = 1}^{n} (\frac{D B H_{i}}{D B H_{J}} \times \frac{1}{D i s t_{i j}}) & (1) \end{array}

Topographic and Geographic Factors

The variability of site characteristics was accounted for by measuring topographical elements (slope and the topographical exposure index), climatic conditions (regional wind speed) and cutblock configuration (forest gap and adjacent cuts). Slope measurements were determined for each experimental unit using a clinometer before cutting. The topographical exposure index (Topex; Quine and White, 1998) was estimated using the spatial database of the Ministère des Ressources Naturelles et des forêts (MRNF) following the method of Ruel et al. (2002). Data were calculated for a 100 × 100 m grid having a limiting distance of 500 m before cutting.

The Canadian Wind Energy Atlas (Environment and Climate Change Canada, 2017) was used to characterize regional mean annual wind speeds at 10 m above ground at a resolution of 5 km (Table 4). To consider any fragmentation caused by anthropic or natural disturbances and the configuration of cutblocks, three variables were created: the radius of forest patches, the distance to adjacent cuts and the simple fetch index (Figure 3b). These variables were inspired from Scott and Mitchell (2005). The width radius of forest patches was quantified using the mean distance extending from the center of the experimental unit to the edge of open areas (harvested areas, trails, lakes, burned areas, avalanche chutes, rock falls and dams). The distance to adjacent cuts was estimated to evaluate the effect of cutover size, and it corresponds to the distance between the edge of the study forest patch (experimental unit) and the nearest forest area for a distance up to 600 m. Finally, the simple fetch index is equal to the ratio between the forest patch radius and the mean distance to adjacent cuts. These variables were estimated in each experimental unit for the 8 main cardinal directions using measurements derived from aerial photographs 10 years after cutting. The adjacent cuts were performed the same year as the experimental silvicultural treatments (2003).

FIGURE 3

Figure 3. (a) Representation of 3-ha experimental units (red squares) and location of permanent plots (yellow rectangle) for a study block, where number refer to: (1) control, (2) mini-strip shelterwood, (3) seed-trees, (4) distant selection, and (5) close selection. (b) Method based on aerial photographs for estimating the mean radius of a forest patch and the mean distance to adjacent cuts in an experimental unit following the eight main cardinal directions.

Data Analysis

All analyses were performed at the tree level, with the exception of a polynomial regression analysis done at the stand level. To achieve the first objective, stand type and silvicultural treatment effects on the probability of tree mortality (MORT), 10 years after cutting, were evaluated using Generalized Linear Mixed Models (GLMM) that accounted for binomial distributions (dead vs. alive). The blocks (BLOCK) and the experimental units (EU) nested in blocks were considered as random effects on the intercept. Fixed effects included stand types (TYPE), treatments (TREAT) and their pairwise interactions (Script S1). The R package lme4 was used, considering a binomial distribution and a logit link function (Bates et al., 2014).

Another set of GLMM was performed for each mortality type (broken, standing dead and overturned) to determine the effects of treatment and stand type on the probability of each mortality type. A binomial distribution of a binary variable identifying trees belonging to the mortality type was considered. In these models, stand type, silviculture treatment and their interactions were considered as fixed effects, while blocks and experimental units were considered as random effects. The GLMMs were performed using the software R 3.4.4 (R Core Team, 2013) and the lme4 package (Bates et al., 2014). In a first step, the significance of each random effect was evaluated using a likelihood ratio test (Bolker et al., 2009). When none of the candidate random effects were significant, a general linear model was fitted using the glm function. Hypothesis testing for each fixed effect was performed using a bootstrapped likelihood ratio test (SSCC, 2016) via the pbnm package (Banghart, 2015) with 500 iterations and applying the model retained after selection of the significant random effects. Prior to this analysis, overdispersion was tested using the function proposed by Bolker (2018). Differences between the levels of significant fixed effects were verified using post-hoc Tukey tests with the package multcompView and the function cld (Graves et al., 2012) on the marginal means estimated by the model using the emmeans package (Lenth, 2018). A second-degree polynomial regression was fitted between the percentage of removed trees and post-cutting frequency of dead trees in each experimental unit at the stand level.

The second objective of the study was addressed with two analyses. First, an exploratory analysis was conducted to identify the most influential factors on the relative probability of trees for the mortality types (broken, standing dead, and overturned) and also living trees (n = 2,051), using the RandomForest (RF) method (Liaw and Wiener, 2002) within the caret package in R (Kuhn and Johnson, 2013). The variables used were DBH, estimated height, slenderness, crown classes, wound severity, basal area, realized harvest intensity, age, average GBC, slope, topex, wind speed, forest patch, adjacent cuts, simple fetch index.

The RF procedure provides an average model from the fit of multiple classification trees obtained from subsamples of observations and variables. The average model can then be used to weigh the importance of the predictors, considering their contribution to each single classification tree to form the final model. The RF method provides robust results for correlated predictor variables and small sample sizes (Cutler et al., 2007). Prior to the analyses, categorical predictors were re-coded using binary dummy variables; variables having near-zero variance and those having high correlation (r > 0.75) were removed (Kuhn and Johnson, 2013). The model training function had a 10-fold cross-validation using Kappa as a performance metric for optimizing the number of variables randomly sampled as candidates at each split. Predictor importance, a measure of the contribution of each predictor to the accuracy of the classification model, was estimated using the package caret (Kuhn, 2012). To estimate the confusion matrix and the accuracy and Kappa indices, 20% of the observations were set aside for validation. The remaining observations (80%) were used for model calibration. Second, to analyse the effects of the most important variables after the RF classification, single nominal logistic regressions were conducted at the tree scale with each factor (realized harvest intensity, basal area, age, adjacent cuts, simple fetch index, slope, average GBC, slenderness, DBH, wound severity, estimated height, topex, density, and competition), using the relative probability of dead trees for each mortality type as a predicted variable (living trees not included). Analyses were conducted using JMP Pro 12 software (SAS Institute Inc., Cary, NC, USA).

Results

Post-cutting Mortality and Windthrow

Prior to any treatment, the abundance of dead trees in the study stands was uniform and low, ranging between 2 and 8% (Table 5). GLMM analysis, based on 2,242 tree observations, showed a significant fixed effect of the cutting treatments on the probability of dead trees 10 years after cutting with no significant effects of stand type or from the interaction of these two factors (P > 0.05, Table 6). Post-hoc comparisons showed that the seed-tree treatment had a significantly higher mortality rate (60%), whereas the control and the shelterwood variants were comparable with mortality ranging from 9 to 43% (Figure 4 and Table 5).

TABLE 5

Table 5. Living and dead tree densities by silvicultural treatment for each stand type before and after cutting (mean ± standard error).

TABLE 6

Table 6. Bootstrapped likelihood ratio test of fixed effects from the Generalized Linear Mixed Model (GLMM) for the probability of post-cutting mortality and mortality types (overturned, broken, and standing dead) in black spruce residual stands, 10 years after treatment.

FIGURE 4

Figure 4. Marginal mean probabilities for post-cutting mortality and overturned trees for each treatment level with their corresponding 95% confidence interval, 10 years after cutting, classified by silvicultural treatment. Vertical bars show the standard error. Letters represent significant differences (P ≤ 0.05) for the treatment effect estimated by post-hoc Tukey analyses, where b > a. Means and standard errors were obtained by the fitted model, using the emmeans package in R.

The post-cutting mortality showed a polynomial relationship with the percentage of removed trees (R² = 39.9%, P = 0.0003). Three groupings of plots were identified, separating into control, shelterwood and seed-trees plots (Figure 5). Post-cutting mortality was lowest in the control plots although two plots were separated slightly from the main group. Shelterwood treatments, especially mini-strips, were highly variable relative to seed-trees in terms of percentage of removed trees. However, the lowest variation for the frequency of dead trees was identified in control plots (8–27%), while shelterwood and seed-trees had a similar variability.

FIGURE 5

Figure 5. Quadratic relationship between tree frequency of post-cutting mortality and percentage of removed trees, 10 years after silvicultural intervention. Filled circles correspond to the study treatment plots. Post-cutting mortality represents the percentage of dead trees within each experimental unit.

In terms of the type of mortality, windthrow accounted for 75% of dead trees, whereas only 25% of trees were standing dead (Figure 6). Breakage was the dominant form of mortality for almost half of dead trees (44%), while overturned trees accounted for a third of post-cutting mortality after silvicultural treatments (31%). The probability of mortality was also more variable in older stands than in younger ones. Additional GLMM analyses evaluated the stand type and silvicultural treatment effects on each of the three mortality types 10 years after cutting (Table 6). Treatments had significant effects on the probability of being overturned, whereas stand type had no effect, either as a main effect or as an interaction with the stand type (P > 0.05). Overturned trees were more prevalent in the shelterwood variants and the seed-tree treatments than in the control sites; however, there was no difference between the silvicultural treatments (Figure 4 and Table 6). Stand type and cutting treatment, either as single effects or as interactions, had no significant effect on the proportion of broken and standing dead trees (Table 6). In terms of the proportion of standing dead trees, neither of the two random effects was significant. Therefore, a generalized linear model analysis was performed for this mortality type. For all the other cases, the random effect of experimental units nested into the blocks was significant; however, the random effect of blocks was not significant (α = 0.05).

FIGURE 6

Figure 6. Probability of post-cutting mortality as represented by stand type, treatment and typology.

Factors Influencing the Probability of Mortality Types

The RF model classified trees by mortality type (living, standing dead, broken, overturned). Based on the classification of trees in the validation dataset, the confusion matrix showed good concordance between the model's prediction and the actual observations, most observations lining up on the diagonal, leading to an accuracy of 0.76 and a Kappa of 0.51 (Table S1). The classification model was able to distinguish between living and dead trees (concordance of 89% between the prediction and the reference). However, model performance was only 60% for broken and overturned trees. Standing dead was the most difficult mortality type to classify correctly, having the lowest level of concordance (10%). The most important variables in the RF analyses were slenderness, estimated height, DBH and average GBC (Figure 7). Most of these tree variables were strongly related: a tree having a short crown will likely have a smaller DBH, a slenderer stem and a lower growth rate. A second set of important variables related to openings within and around the treated plots (harvest intensity, basal area, patch, adjacent cuts and fetch index). However, other variables such as age and wind speed had the lowest importance in the RF analyse.

FIGURE 7

Figure 7. Importance index of the 12 most important predictors according to the RandomForest multinomial model of tree distribution by mortality type and living trees.

Complementing the previous analyses, simple nominal logistic regressions described how the relative probability of the mortality types varied according to the stand, geographical and tree variables (Figure 8 and Figure S3). Trees with a high slenderness (height/DBH ratio) were more vulnerable to breakage, while more tapered stems were more likely to be overturned. Similar trends were observed for DBH and estimated height. In all cases, higher predictor values were associated with a lower proportion of broken trees and a higher proportion of overturned trees. DBH was important in explaining standing dead mortality; standing dead trees made up 30% of the mortality for trees having a diameter <10 cm, a percentage 10× greater than for trees having a DBH >26 cm. Standing dead trees made up 30% of slow growing trees (GBC = 0.2 cm year⁻¹) but represented <5% of trees having a GBC = 0.7 cm year⁻¹. Finally, 50% of trees having low GBC values were overturned, and almost 90% of trees having high GBC values were classified as broken trees. Wound severity also influenced mortality type, although this had a lower importance in the RF analyses. Sixty percent of post-cutting mortality for undamaged trees with a DBH <12 cm were broken trees, while 80% of the mortality for severely damaged trees with larger DBH values corresponded to overturned trees.

FIGURE 8

Figure 8. Variation in the probability of mortality types at the tree level (n = 672) according to different stand, geographical, and tree variables as obtained through nominal logistic regression.

Half of the topographical and climatic variables were not significant (P > 0.05), e.g., wind speed (P = 0.1608). However, geographical variables, such as the width of adjacent cuts or the simple fetch index, were determinant in explaining the changes in the probability of mortality types (Figure 8). Among dead trees, overturning was more frequent in stands having large areas cut around the residual forest (>300 m), this pattern being the opposite for broken trees. The logistic regressions using Topex, initial stand density and Hegyi competition index as predictors were significant (P < 0.0001), although they had low importance values in the RF classification model (Figure S3).

Discussion

Post-cutting Mortality and Windthrow

Partial cutting is likely to play a more important role in the future management of boreal forests (Gauthier et al., 2009; Kuuluvainen, 2009; Bose et al., 2014; Bergeron et al., 2017). Therefore, quantifying losses due to post-cutting mortality becomes essential for the planning of the long-term timber supply. Post-cutting mortality, 10 years after experimental shelterwood treatments in unmanaged black spruce stands, ranged from 21 to 43% of residual trees, similar to previous studies in boreal and cool temperate zones (Ruel et al., 2003; Scott and Mitchell, 2005; Rosenvald and Lohmus, 2008; Thorpe et al., 2008; Lavoie et al., 2012; Urgenson et al., 2013). Our results showed less variability in post-cutting mortality than other studies, e.g., 15 to 46% for Norway spruce [Picea abies (L.) Karst] stands (Hautala and Vanha-Majamaa, 2006). This reduced variability could be explained by our study area having more uniform climatic, geographical and stand characteristics.

Mortality after silvicultural treatments was proportional to harvest intensity, consistent with previous research (Rosenvald and Lohmus, 2008; Steventon, 2011; Lavoie et al., 2012; Urgenson et al., 2013; Anyomi and Ruel, 2015). Post-cutting mortality was greatest for the treatment having the highest harvesting level compared to the control and other treatments. In the boreal mixed wood forests in Alberta, Bladon et al. (2008) reported post-cutting mortality levels 2.5–4× higher in stands having a 10% retention level relative to untreated plots, matching our results obtained for seed-tree cuttings. Our experimental shelterwood cuttings had a post-cutting mortality around 30% at 10 years after cutting, confirming that 50% harvest intensity successfully reduced the proportion of tree loss compared to the seed-tree system. The lack of difference between shelterwood variants is linked to the within-treatment variation associated with some minor differences in harvest intensity during operational harvesting. This reflects the random variability expected from “real-life” mechanized operations with no tree or trail marking prior to the harvest. The causes of such variations were related to site topography, which does not allow regularly spaced trails, and the different choices of operators in the selection of trees when applying the silvicultural treatments. These elements are part of the experimental error, and are assumed as such. The absence of significant difference in post-cutting mortality between the control and shelterwood plots may be related to snow damages that occurred in two sites in 2004. It is likely that this event increased mortality in the control plots; snow damage occurs often as broken trees (Wallentin and Nilsson, 2013). Similar results were obtained after partial cutting with small retention (20–30%) in a natural, structurally diverse, multispecies forest in British Columbia (Coates et al., 2018).

The effect of stand type was non-significant, in contrast to the findings of Anyomi and Ruel (2015). We argue that stand types in our study were more similar in terms of age and density because the difference between the oldest block of younger stands and the youngest block of older stand were of <20 years and 350 trees·ha⁻¹ (Figure S2), this could be considered as a limitation of this experimental design (Montoro Girona et al., 2016). In fact, our old stands would correspond only to the beginning of the second-cohort phase of Bergeron et al. (2007). Many previous studies dealing with windthrow after partial cutting were conducted in managed forests, reflecting the fact that partial cutting has not been used to a large extent in old unmanaged boreal forests. As such, post-cutting mortality in these forests, especially in black spruce stands, remains unquantified for a wide range of silvicultural scenarios (Coates, 1997; Ruel et al., 2003; Bladon et al., 2008; Thorpe et al., 2008; Lavoie et al., 2012; Hämäläinen et al., 2016). Our study provides a rigorous assessment of post-cutting mortality within previously unmanaged black spruce forests for several partial cutting treatments and coupled with replicated experiments over comparable conditions.

After natural or anthropic disturbances (e.g., harvest operations), trees are expected to have an increased risk of mortality (with wind as the main cause of death) when stand density is reduced or new edges are created (Thorpe et al., 2008; Urgenson et al., 2013). In our study, most of this post-cutting mortality (around 80%) was due to windthrow, as has been previously documented (Ruel et al., 2003; Lavoie et al., 2012). Thus, the applied shelterwood and seed-tree treatments contributed to an increased vulnerability to windthrow through an increase in overturning, without any significant effect on the other types of mortality. As only overturning was significantly affected, treatment effects on overturning parallel those for total mortality. Windthrow had a major impact on post-cutting mortality according to other studies that showed how 20–30% of the trees were affected by this natural disturbance in our study area between 2006 and 2013 (De Grandpré et al., 2009; Waldron et al., 2013). However, a wet snow event also affected two study sites and the added weight of snow increases the turning moment when the tree starts leaning due to the wind action. Black spruce has a narrow crown with drooping branches that generally sheds the snow, except in certain episodes of wet snow occurring in late fall or early winter. Thus, even though snow added to the wind action, this damage could not be separated from the wind impact alone and both were included as windthrow. Considering the future, we recommend improving the understanding the impact of wet and heavy snowfall events on tree mortality since it is expected that this kind of events will increase in a climate change context (Lute et al., 2015; Woods et al., 2017).

Factors Influencing the Probability of Mortality Types

Forecasting and minimizing post-treatment mortality require an understanding of the relative importance of the various mortality types and of the factors influencing this relative importance. Mortality is one of the most complex ecological processes to model in post-cutting stand dynamics. Despite this complexity, our RF analysis was able to account for 89% of the mortality at the tree level. Previous studies have only explained <60% of the variation (Solarik et al., 2012; Anyomi and Ruel, 2015). This performance may be due to our inclusion of a broad set of potential predictors and the range of harvest intensities in our experimental design. Unexplained variation can be attributed to the randomness of wind factors and other unaccounted tree-level factors, such as decay. The segregation by mortality type suggests differences in the complexity and explanation of these phenomena; the RF analysis was quite efficient in predicting overturning and stem breakage but performed poorly for standing dead trees. Based on our results, overturning, standing dead and breakage correspond to three distinct processes, thereby justifying the relevance of analyzing them separately. Mechanistic models are available to estimate the increased risk of overturning and stem breakage after partial harvesting (Gardiner et al., 2008), but similar models are lacking for standing dead trees.

In terms of the relative importance of mortality types, overturning and breakage showed an opposite response to the different variables. Overturning was the only type of mortality directly associated with the treatment effect. Overturning increased with harvest rate, as observed by Jönsson et al. (2007). Losses through overturning also increased with the size and proximity to adjacent clear-cuttings, also matching previous observations (Ruel et al., 2003; Scott and Mitchell, 2005). This trend might be explained by the increased wind speed at new edges around large gaps (Gardiner et al., 2005). However, the higher wind speed would normally also lead to increased stem breakage, but the trend for this type of damage was the opposite. In many previous studies, an assessment of wind damage was presented without specifying the type of damage. Since no solid explanation about these opposite responses can be suggested, additional studies separately considering these two types of mortality remain relevant.

Based on the RF analysis, tree-level factors describing tree size and vigor had a much greater effect on the relative importance of mortality types in comparison with treatment. In our study, tall, tapered, faster-growing dominant trees were more vulnerable to overturning in agreement with previous studies (Bladon et al., 2008; Thorpe et al., 2008; Steventon, 2011; Lavoie et al., 2012; Solarik et al., 2012). Wound proportions influenced overturning probability, in contrast to Urgenson et al. (2013). Trees damaged by machinery would be more abundant in proximity to skid trails and have wounds such as abrasion, crown damage and root compaction (Gea-Izquierdo et al., 2004; Thorpe et al., 2008). With increased wind speeds in the harvest trails due to a funneling effect (Ruel et al., 2001), trees (especially damaged trees) located along the edges of residual strips are more exposed and vulnerable to overturning than trees located toward the interior of these strips (Gardiner et al., 1997; Ruel et al., 2001; Thorpe et al., 2008).

In this study, breakage was the dominant type of mortality, although it was not significantly affected by treatment. As mentioned previously, a wet snow event in 2004, likely increased the occurrence of broken trees across most of the study sites. In black spruce stands, the smaller diameter trees were more vulnerable to breakage, while trees having low slenderness ratios were more resistant to stem breakage, as suggested by Wood (1995) (see Valinger and Fridman (2011) for an opposite result). However, Dunham and Cameron (2000) found similar slenderness ratios for all trees among mortality types. Thus, there is no solid evidence able to explain the relationship between slenderness ratios and breakage, probably because although theoretically more tapered trees should be more resistant to breakage it is necessary to also account for changes in wind loading with changes in crown size (Gardiner et al., 1997). The risk of breakage is subject to many factors, including the stature of the trees, the aerodynamic properties of the stand, the general and local incidence of winds, and soil conditions (Cremer et al., 1982), thus, more work is needed to better understand this phenomenon after silvicultural treatments.

Standing dead mortality was the least frequent yet the most complex type of mortality, in contrast with Anyomi and Ruel (2015) who found this type to be the most important. Trees can die standing through a stand self-thinning process, senescence, insect or pathogen attacks. In our study, this mortality type was relatively more important for trees having a small DBH and GBC in uncut stands. This suggests that this type of mortality was more important in suppressed trees and that it was related to self-thinning through intraspecific competition, as described by Lussier et al. (2002) in unmanaged black spruce stands. Partial harvesting could reduce suppressed tree mortality by favoring an increase in resource availability (light, nutrients, etc.) and decreased competition (Prévost and Dumais, 2014). This pattern has been previously described in studies demonstrating the strong growth response of most of the previously suppressed trees (Montoro Girona et al., 2016, 2017). However, the probability of standing dead trees was similar between the treated and untreated plots (P > 0.05). This finding contrasts with that of Anyomi and Ruel (2015) who observed that standing mortality increased with harvesting intensity. Their study was conducted in old irregular stands where self-thinning and senescence could act at the same time in different parts of the stand.

Standing dead mortality in our stands was lower than that found by Thorpe et al. (2008) and Anyomi and Ruel (2015). These latter studies were conducted in stands having trees of variable age (not even-aged stands) likely marked by a greater tree senescence than within our study sites. Standing dead mortality can also be related to physiological stress resulting from increased light and wind exposure as the stand becomes more open after harvesting or due to other types of mortality, e.g., overturned trees (Thorpe et al., 2008). In old-growth black spruce stands, Anyomi and Ruel (2015) noted that windthrow and standing dead dynamics had common facilitating factors. However, given the poor performance of the RF analysis in explaining this type of mortality, additional studies are required to better understand the underlying processes as well as to consider the interaction between silvicultural treatments with the factors involved in the different types of mortality (Elith et al., 2008).

Conclusion

This study quantified black spruce post-cutting mortality along a silvicultural gradient, from untreated to seed-tree cuttings, in a previously unmanaged boreal forest. Our study demonstrated that post-cutting mortality increased with harvest intensity and that shelterwood treatments of 50% tree removal reduced windthrow damage in comparison with seed-tree cuttings. This study highlighted the complexity of mortality in these stands as total mortality was composed of three types that each respond differently to harvesting. The increase of mortality with cutting intensity was directly linked to an increase in overturning. Mortality is often considered to be a stochastic event (Taylor and MacLean, 2007). However, our results indicate that post-harvest losses were largely dependent on preharvest stand conditions, tree characteristics and harvest intensity. Most topographical and geographical characteristics had a minimal effect with only the nature of adjacent stands being important.

Given the observed differences between breakage, standing dead and overturning mortality, we suggest that future research must account for these differences in mortality to improve our understanding of post-partial cutting mortality. As mechanistic models already exist for stem breakage and overturning, emphasis should be placed on understanding standing mortality.

This research has contributed to improve the evaluation of the studied treatments as a silvicultural alternative for sustainable management in Canadian boreal forests. So far, the stands submitted to the experimental shelterwoods showed an excellent growth response, high levels of regeneration (specially mini-strip shelterwood), and intermediate levels of post-cutting mortality according to our previous studies (Montoro Girona et al., 2016, 2017, 2018b). However, more studies should be conducted to evaluate the economic cost and the effect on biodiversity and wildlife and to define the conditions of applications of these silvicultural treatments in a global change context.

Author Contributions

MMG and J-CR: conceived and designed the study; MMG: carried out lab and fieldwork; MMG and J-ML: analyzed the data; MMG, J-CR, J-ML, and HM: provided interpretation; MMG: project administration; HM: supervision; J-ML and HM: established the long-term experiment and provided substantial input to the design; HM, J-ML, and MMG: contributed to funding; MMG: wrote the draft version of the manuscript. All authors contributed to the writing and revision of the manuscript.

Funding

This project was supported by the Fonds de Recherche du Québec–Nature et Technologies (FQRNT), the Programme de mise en valeur des ressources forestières (MFFPQ), the Canadian Wood Fibre Centre of the Canadian Forest Service (Natural Resources Canada), and the Forest Complexity Modeling program of the Centre for Forest Research.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank to R. Gagnon and G. Grosbois for fieldwork assistance. We extend special thanks to L. De Grandpré, A. Leduc, S. Rossi, C. Krause, A. Hof, and J. Hjältén for their suggestions on the manuscript, S. Gosselin and F. Marchand for technical support, M. Hay for reviewing the English of the manuscript as well as the reviewers and the editorial board for the help in improving the quality of this manuscript. This manuscript is part of the Ph.D. thesis of MMG (Montoro Girona, 2017).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2019.00004/full#supplementary-material

References

Aakala, T., Kuuluvainen, T., Grandpré, L. D., and Gauthier, S. (2006). Trees dying standing in the northeastern boreal old-growth forests of Quebec: spatial patterns, rates, and temporal variation. Can. J. For. Res. 37, 50–61. doi: 10.1139/x06-201