Success Factors for Experimental Partial Harvesting in Unmanaged Boreal Forest: 10-Year Stand Yield Results

Over the past two decades, partial harvesting has been increasingly used in boreal forests as an alternative to clearcutting to promote irregular stand structures and maintain a balance between biodiversity preservation and continued timber production. However, relatively little is still known about the silvicultural potential of partial harvesting in Canada’s boreal forest, especially in areas prone to organic matter accumulation (paludification), and most prior research has focused on biodiversity responses. In this study, we assess the effects of partial harvesting on stand development (recruitment, growth, and mortality) ten years after harvesting in previously unmanaged black spruce stands and quantify its effectiveness in reducing the impacts on ecosystem structures. Our analyses revealed that pre-harvest stand structure and site characteristics, especially initial basal area, sapling density, tree diameter, and organic layer thickness (OLT) were major factors involved in stand development ten years following these partial harvesting treatments. Depending on pre-harvest structure and site characteristics, partial harvesting can result in either an increase in post-harvest tree recruitment and growth or a loss of stand volume because of standing tree mortality. To increase the chances of partial harvesting success in ensuring an increase in decennial stand yield after harvest in black spruce forest stands, sites prone to paludification (i.e., where OLT >17 cm) should be left unharvested. This study illustrates the importance of taking into account pre-existing structure and site characteristics in the selection of management strategies to maximize the potential of partial harvesting to achieve sustainable forest management in black spruce stands.


Introduction
In recent decades, the use of forestry practices other than clearcutting has been explored as a means to mitigate the major issues facing sustainable management of boreal forests, i.e., the impacts on ecosystem structures, functions and dynamics, and the preservation of biodiversity [1][2][3][4][5]. The predominant use of clearcutting in the boreal forest during the last half-century, and the lack of consideration of unmanaged black spruce stands [57]. This study aimed to assess the effects of partial harvesting on stand development ten years after harvesting and quantify its effectiveness for maintenance of stand structural attributes. Our specific objectives were to: (i) Determine whether net stand yield is primarily a product of residual tree growth, recruitment into the merchantable tree class, or tree mortality ten years after partial harvesting. (ii) Identify the main factors influencing recruitment of trees into the merchantable tree class, mortality of residual trees, and tree growth. Factors examined included harvesting intensity, site characteristics, and pre-harvest stand structural characteristics. (iii) Determine how site and stand structure characteristics influence recruitment of merchantable trees, residual tree mortality, and tree growth after harvesting individually to determine thresholds to help managers determine whether partial harvesting can fulfill stand yield objectives.

Study Area
The study took place in the northern portion of the Clay Belt of northeastern Canada ( Figure 1, Table 1) where forest stands are dominated by black spruce in association with companion species such as jack pine (Pinus banksiana Lamb), trembling aspen (Populus tremuloides Michx.), balsam fir (Abies balsamea [L.] Mill), and white birch (Betula papyrifera Marsh) [58]. This region is prone to paludification as a result of poorly drained clay-dominated soil. The Clay Belt is a major physiographic region created by the clay deposits left by Lakes Barlow and Ojibway, and by the Cochrane till in the northern part, i.e., a compact till made up of a mixture of clay and gravel, after their maximum extension during the Wisconsin glaciation, approximately 8000 years before the present [59]. The topography is generally flat, and the climate is moderately humid and cold (between 1981 and 2010 mean annual precipitation and mean annual temperature were 928 mm and 1 • C, respectively) [60]. Large wildfires are an integral part of the natural dynamics in this region. The fire cycle was about 135 years between 1850 and 1920, and has since extended to around 400 years, and average stand age is currently around 150 years [61]. The forest mosaic has been shaped by both fire and harvest over the last 50 years [35,62].  Table 1. Site level variables based on the data contained in eco-forestry maps [63] and characteristics of partial harvesting treatments and stand characteristics for each site before and after harvesting given as mean (standard error). * There is no pre-harvest inventory at this site. HARP: Harvesting with Advance Regeneration Present; CCCC: Cut with Conservation of Canopy Cover.

Experimental Design
The study is part of the Réseau d'expérimentation des coupes partielles en Abitibi (RECPA, http: //recpa.uqat.ca/), i.e., the Abitibi Partial Harvesting Experimental Network, which comprises 11 sites, established between 1998 and 2007 with the objective to better understand the short-and mid-term effects of partial harvesting on ecosystem functioning, biodiversity, and stand yield [14,64]. Five sites (i.e., Muskushii, Fenelon, Maïcasagi, Gaudet, and Puiseaux) for which ten years post-harvest data was available are included in this study ( Figure 1 and Table 1). These five sites are characterized by stands >120 years old, with an irregular diameter structure, that were established after stand-replacing fires [53,57]. These stands had never been harvested before and their irregular structure made them suitable for testing partial harvesting treatments [65].
One partial harvest 25 ha block was established on each site. In this study, the data set included 93 permanent sampling plots (11.28 m radius, 400 m 2 ) where harvested basal area varied from 45 to 80%, i.e., retained basal area varied from 20 to 55% (Table 1). Prior to harvesting, sampling plots were randomly established across the five selected sites, i.e., 51 plots in Muskushii, 18 in Fenelon, 17 in Maïcasagi, 17 in Gaudet, and 17 in Puiseaux. One of the five sites (Muskushii) was harvested with Coupe avec Protection des Petites Tiges Marchandes (CPPTM treatment) (known as Harvesting with Advance Regeneration Present or HARP in Ontario), which is a partial harvest that involves the removal of stems with diameter at breast height (DBH) greater than a minimum diameter limit (in our case, >14 cm DBH). HARP aims at promoting the recruitment of advance regeneration and merchantable tree following canopy removal by retaining small-diameter trees (in our case <14 cm DBH) on the cutover. HARP is operationally used in irregular boreal stands with saplings and small merchantable stems [15,24,25,66]. At the other four sites (Fenelon, Maïcasagi, Gaudet, and Puiseaux), stems were harvested with Cut with Conservation of Canopy Cover (CCCC treatment or variable retention). CCCC is a type of partial harvesting where stems are harvested from all diameter classes to maintain a similar pattern as before harvest. This second technique generally results in a more even cover and the retention of larger trees than HARP, and is more dependent on the accessibility of the land and operators experience [57,67].

Forest Inventory Data
Forest inventories were carried out in the autumn before harvest, during the first spring following harvest, and then again ten years after harvest (Table 1). Prior to harvesting, we measured in each plot organic layer thickness (i.e., from the soil surface to the mineral soil or bedrock) and diameter at breast height (DBH) of all live commercial trees (i.e., DBH ≥9 cm). We then computed initial mean DBH prior to harvesting, in addition to initial stand density and initial stand basal area of live commercial trees ( Table 2). At remeasurement in the spring of the year following harvesting, we measured mean DBH of live residual trees, and also determined residual stand density and basal area ( Table 2). This allowed us to evaluate harvesting intensity (i.e., percentage of harvested basal area). To assess the abundance of saplings (1 cm ≤ DBH < 9 cm) corresponding to the measured forest canopy, one 40 m 2 plot was established at the center of each 400 m 2 plot in which all saplings were identified to species, and their density computed (Table 1). Trees and saplings were marked with paint or tags and their location recorded with GPS to ensure that the same trees were sampled in subsequent sampling. Ten years post-harvest, tree mortality, residual tree growth, and recruitment into the merchantable tree class were evaluated. Decennial net stand yield was defined as the net gain in basal area ten years after harvesting compared to residual basal area at remeasurement in the spring of the year following harvest ( Table 2).
Recruitment into the merchantable tree class corresponded to the density of marked-live saplings that reached DBH >9 cm during the ten years after harvest ( Table 2). Recruitment was measured at the plot scale and then expressed in density of stems per hectare.
Tree mortality (TM) was considered on the basal area of dead trees ten years after harvest ( Table 2). Mortality was calculated only on tallied, live stems that had DBH >9 cm before the first growing season of the year following harvesting.
Residual tree growth was considered as the average annual radial growth (ARG) of living merchantable trees ten years after harvesting ( Table 2). Annual radial growth is the tree growth per year as measured by the increase in tree radius. The annual radial growth was measured using permanently numbered post-harvest trees that were still alive ten years later. Annual radial growth was calculated at tree scale as: DBHt and DBHr corresponded respectively to mean diameter at breast height ten years after harvesting, and mean diameter at breast height at remeasurement before the first growing season in the year following harvesting.

Data Analyses
First, the influence of residual tree growth, recruitment to the merchantable tree class, and tree mortality on net stand yield, ten years after harvesting, was evaluated using multiple regression models (nlmer) used to fit linear mixed-effect models ( Table 3). The R package lme4 was used, with a Gaussian distribution and Satterthwaite's method (lmerMod LmerTest) [68]. The potential effect of site location was considered a random factor. Fixed effects included residual tree growth, recruitment, and tree mortality (Appendix A). We first verified the effect of each explanatory variable individually (Appendix A). Tree recruitment was log-transformed to homogenize variances. Based on univariate models with the smallest Akaike's information criterion (AIC) values, we subsequently built models with two variables. The best two variable models were then selected and nested in an improvedmodel with three variables (Appendix A). Second, to identify structural and site factors influencing tree growth, recruitment of merchantable trees, and residual tree mortality, three other sets of lmer models were performed. In these three model sets, silvicultural treatment intensity (i.e., percentage of harvested basal area), site characteristics (e.g., organic layer thickness), and residual stand structural characteristics (tree basal area, sapling density, mean DBH of trees) before the first growing season of the year following harvest were considered fixed effects, and site was considered a random effect (Tables 4-6). We verified the multi-collinearity of the explanatory variables using the variance inflation factor (VIF) [69]. Third, to determine site and stand structure characteristic thresholds, a recursive partitioning approach was used [70], with four non-parametric regression trees [71,72] that were constructed using the "package rpart, function rpart, method anova" in CRAN R 3.6.0 [73]. Based on the results of the first two objectives (Tables 3-6; Appendix A), four regression trees were derived to partition the data and determine how the factors identified in the first two sections of analyses, and that significantly affecting stand yield and its components, control the recruitment of merchantable trees, tree mortality, and tree growth. To partition the data and identify the interactions among tree mortality, residual tree growth, and recruitment into the merchantable tree class leading to an increase in stand yield ten years after harvest, the first regression tree was derived using these three yield components. Then, we developed the last three regression trees using the factors (identified in the first two sections of analyses) that significantly affected the recruitment of merchantable trees, the residual tree mortality, and the tree growth. A regression tree analysis is a type of repeated regression analysis in which the data are recursively partitioned to maximize the homogeneity within the resulting groups using different levels of the explanatory factors. In the single node at the top of the tree, the complete data set is split according to the most influential factor. Then, repeated binary splitting of the data set of the preceding node grows the tree. Each split is defined by a simple rule, usually based on a single explanatory variable cut off, and forms two branches.
The R statistical software (v 3.6.0) was used for statistical analyses with a significance level of a ≤0.05 [74].

Decennial Yield and Partial Harvest Success
Net stand yield ten years after harvesting ranged from −10.6 to 6.7 m 2 ha −1 decade −1 , with an average of −0.8 m 2 h −1 decade −1 and a median of −0.4 m 2 h −1 decade −1 . One-third (33%) of the sampling plots had a reduction in basal area ten years after partial harvesting, whereas a gain in basal area was recorded for 60% of the plots (Figure 2).

Importance of Stand Components on Decennial Stand Yield Following Partial Harvesting
The best model selected with multiple regression analyses for the decennial net stand yield (Appendix A, AIC = 304.2) included the three yield components, i.e., residual tree growth, recruitment into the merchantable tree class, and tree mortality (Table 3). Net stand yield significantly increased with the growth of residual trees and the post-harvest recruitment into the merchantable tree class, whereas it significantly decreased with increasing residual tree mortality ( Figure 3). However, based on the slopes of these relationships, mortality had the greatest influence among the three components on stand yield (Figure 3).
Based on the three yield components, the first regression tree partitioned the data and determined how these stand components influence the positive decennial stand yield (Figure 4). The first node split the plots based on tree mortality with a threshold of 2.7 m 2 ha −1 (which is the equivalent of 44% of basal area occupied by dead trees) into two groups. Of the total number of plots, 62% were below this threshold with an average gain in stand yield of 0.97 m 2 ha −1 , whereas 38% of plots were above this threshold with an average loss of stand yield of 1.6 m 2 ha −1 . Then the group with a gain in basal area was divided in terms of residual tree recruitment into two sub-groups. In the high (≥263 trees ha −1 ) recruitment group, a gain in decennial stand yield was observed in 12% of the plots, whereas a positive stand yield was observed in 50% of plots if tree recruitment was <263 trees ha −1 . In the low (<263 trees ha −1 ) recruitment sub-group, when tree recruitment was <88 trees ha −1 but tree mortality basal area was <1.8 m 2 ha −1 (which is the equivalent of 30% of basal area occupied by dead trees), a positive decennial stand yield was observed in 13% of partially harvested plots. In contrast, a residual tree mortality ≥1.8 m 2 ha −1 with tree recruitment <88 trees ha −1 led to a negative decennial stand yield in 9% of the plots. Moreover, a gain in decennial stand yield was observed in 19% of partially harvested plots when tree recruitment was ≥88 trees ha −1 and residual tree growth was >2.9 cm 2 yrs −1 . In contrast, a tree growth ≤2.9 cm 2 yrs −1 led to a gain in decennial stand yield in only 10% of the plots. This first regression tree indicated that tree mortality and recruitment were the most influential components explaining variations in stand yield ten years after harvesting, whereas the influence of residual trees growth appeared marginal (Figure 4).  . Regression tree predicting decennial stand yield (m 2 h −1 ) after partial harvesting according to yield components, i.e., recruitment into merchantable tree class (RTR), residual tree mortality (RTM), and growth (RTG) over the ten years since partial harvesting.

Factors Influencing Merchantable Tree Recruitment
The nonlinear mixed-effect models showed that merchantable tree recruitment ten years following harvest was significantly influenced by several factors such as organic layer thickness, residual tree basal area, residual sapling density, mean DBH of residual trees, and harvesting intensity (Table 4).
Ten years following partial harvesting, residual tree recruitment significantly increased with increasing residual sapling density (Figure 5a), whereas it significantly decreased with increasing organic layer thickness, residual tree basal area, mean DBH of residual tree, and harvesting intensity (Figure 5b-e). Specifically, based on the slopes of these relationships between residual tree recruitment and the explanatory variables ( Figure 5), we observed that merchantable tree recruitment is mostly influenced by sapling density and harvesting intensity (Figure 5a,c), followed by organic layer thickness and residual tree basal area (Figure 5b,d), and to a lesser extent by mean residual diameter at breast height (Figure 5e). Figure 5. Variation in merchantable tree class recruitment ten years after harvesting as a function of (a) sapling density, (b) organic layer thickness, (c) harvesting intensity, (d) residual basal area, and (e) residual mean DBH as predicted by nonlinear mixed-effect models ( Table 4).
The second regression tree indicated that residual sapling density, organic layer thickness (OLT), residual mean DBH (DBHr), mean residual basal area, and harvesting intensity influenced merchantable tree recruitment ten years after harvesting ( Figure 6). A first node split the plots into two sub-groups based on the DBHr. Merchantable tree recruitment was 181 trees ha −1 in 33% of the stands where DBHr <118 mm, whereas forest stands where DBHr was ≥118 mm showed a recruitment of only 71 trees ha −1 . In both cases, recruitment into the merchantable tree class was important when OLT was <16 cm (22% of the plots, n = 20) or <17 cm (41% of the plots, n = 38) with respectively 231 and 102 trees ha −1 ten years after harvesting, whereas merchantable tree recruitment was less important (<100 trees ha −1 ) when OLT ≥16 cm (12% of the plots, n = 11). When DBHr <118 mm and OLT <16 cm, the regression tree showed an increase in merchantable tree recruitment from 162 to 300 trees ha −1 with respect to the density of saplings (RSD). In contrast, when OLT <17 cm, a net increase in merchantable tree recruitment was favored where the DBHr was <150 mm (33%, n = 31), and residual basal area (RBA) was ≥4.2 m 2 ha −1 (26%, n = 24), and finally when RSD ≥815 saplings ha −1 (16%, n = 15). Figure 6. Regression tree predicting recruitment into merchantable tree class (trees ha −1 decade −1 ) ten years after harvesting according to significant explanatory factors obtained using a nonlinear mixed-effect model. DBHr: mean residual diameter at breast height; OLT: organic layer thickness; RBA: residual basal area; RSD: residual sapling density.

Factors Influencing Tree Mortality
Residual tree basal area, mean residual tree DBH, residual sapling density, and harvesting intensity were the four major factors explaining residual tree mortality ten years following harvest (Table 5).
Ten years following partial harvesting, residual tree mortality significantly increased with increasing residual tree basal area, mean DBH of residual tree, and harvesting intensity, whereas it significantly decreased with increasing residual sapling density (Figure 7). Specifically, based on the slopes of these relationships, we observed that residual tree mortality results mostly from residual tree basal area and residual mean DBH (Figure 7a,b), followed by harvesting density (Figure 7d), and to a lesser extent by residual sapling density (Figure 5e). Variation in the residual tree mortality ten years following partial harvesting as a function of (a) residual basal area, (b) residual mean DBH, (c) sapling density, and (d) harvesting intensity as predicted by nonlinear mixed-effect models ( Table 5).

Factors Influencing the Residual Tree Growth
The nonlinear mixed-effect models showed that only sapling density had a significant positive relationship with residual tree growth ten years following partial harvesting (0.0 > CI < 0.02; Table 6; Figure 9).
The forth regression tree showed how residual sapling density influenced the variation in the decennial tree growth after harvesting ( Figure 10). This regression tree showed that the residual sapling density (RSD) was determinant in explaining the probability of the decennial tree growth after harvesting. Tree growth ten years after harvesting was 5.9 cm 2 yrs −1 (17% of the plots; n = 16) where RSD was ≥1875 saplings ha −1 . When RSD was <1875 saplings ha −1 , the tree growth varied between 2.4 and 4.1 cm 2 yrs −1 (Figure 10). Figure 9. Variation in tree growth (cm 2 yrs −1 decade −1 ) ten years following partial harvesting as a function of sapling density as predicted by nonlinear mixed-effect models ( Table 6). Figure 10. Regression tree predicting residual tree growth ten years after harvesting according to significant explanatory factors obtained using a nonlinear mixed-effect model. RSD: residual sapling density.

Discussion
Quantifying the response of stands to partial harvesting treatments is a critical step towards assessing their potential from either a conservation or timber supply perspective in boreal forests [15,23,32,[53][54][55][56]. In previously unmanaged black spruce forest stands, this study shows that depending on pre-harvest structure and site characteristics, especially initial basal area, sapling density, tree diameter, and organic layer thickness, partial harvesting ten years following experimental treatments can result in either an increase in post-harvest tree recruitment and growth or a loss of stand volume because of standing tree mortality. However, our results also show that residual tree growth plays a relatively minor role in stand yield ten years following partial harvesting treatments compared to both tree mortality and recruitment into the merchantable tree class.
First, merchantable tree recruitment ten years following harvest was a complex phenomenon significantly influenced by several factors. Merchantable tree recruitment significantly decreased with increasing thickness of the organic layer, residual tree basal area, and mean DBHr, whereas it increased with increasing residual sapling density. The thickness of the organic layer has often been identified as an important factor influencing the structure and the productivity of black spruce forests [75][76][77]. In this study, in plots where the thickness of the organic layer was ≥17 cm, black spruce establishment was limited (i.e., ≤22 stems ha −1 ), leading to low productive stands, likely due to relatively low soil nutrient availability locked up in the poorly decomposed organic matter that makes up the forest floor [78,79]. The accumulation of a thick organic layer-or paludification [51,52]-causes both anoxic conditions, i.e., soil moisture and water table level to rise, and low access to nutrients, i.e., anchoring and intertwining of the black spruce root system to decrease, which leads in turn to decreased aboveground biomass and productivity [79][80][81]. Moreover, it has recently been suggested that trembling aspen establishment and growth are also limited on sites where soil organic layer thickness was ≥26 cm [81]. This suggests that in boreal forests, several tree species respond similarly to soil organic layer thickness likely due to its negative effects on soil physical and chemical properties associated with the rise of the water table, especially when it is located above the mineral soil surface [82]. In contrast, in plots where the thickness of the organic layer was <17 cm, tree establishment rate and growth were generally greater, probably due to the high availability of fine roots in the mineral soil than in the organic layer [76,83]. Our results suggest that site productivity is one of the main factors limiting the ability of roots to reach the mineral soil. In fact, the activity of the microbial biomass generally is influenced by the anaerobic conditions induced by the shallow water table rise [44,76,77,84].
Pre-harvest structure characteristics including the diameter of the residual trees, residual tree basal area, and sapling density are also factors influencing merchantable tree recruitment mainly in sites where the organic layer is <17 cm thick. Indeed, when the mean DBHr <15 cm, forest stands with a residual basal area ≥4.2 m 2 ha −1 are more favorable to merchantable tree recruitment, especially if ≥815 saplings per hectare are preserved after harvesting. Generally, on thin organic layers, small merchantable stems (9 < DBH ≤ 14 cm) present a better survival rate [85] regardless of harvesting intensity and initial tree density, because in coniferous forests, following canopy opening by partial harvesting, residual saplings tend to invest more in mechanical support, and to allocate more resources to root growth to counterbalance increased wind and evaporative demand of the exposed canopy [86,87]. Knowing that dense stands generally have a lower water table than less dense sites because of evapotranspiration, partial harvesting promotes sapling growth and their recruitment into the tree class. This suggests that in sites where organic layer is thinner than 17 cm (i.e., rich soil conditions), the way to ensure a tree recruitment into merchantable class will also depend on pre-existing stand structure, especially on stand initial basal area and presence and preservation of saplings.
Second, as observed in others coniferous boreal forests, our results indicate that residual basal area and mean diameter of the residual trees, and the density of residual saplings, are the main factors influencing stem mortality among residual trees ten years after partial harvesting [25,32,88]. Generally, after natural or anthropic disturbances (e.g., harvest operations), new edges and gaps are created and stand density is reduced, leading to an increased risk of tree mortality or windthrow [24,31].
Our results indicate that ten years post-harvest, tree losses were largely dependent on pre-harvest structure characteristics and harvest treatment intensity, which is in agreement with previous study results in unmanaged black spruce stands [32]. Low mortality was observed when residual tree basal area was <9 m 2 ha −1 (open stands) and DBHr was <12.2 cm. When DBHr was >12.2 cm, low mortality was observed only on plots where initial basal area was harvested at more than 48% and the density of residual saplings was >1250 saplings ha −1 . This could result from the important density of saplings which can act as buffer zone to the risk of windthrow [75,89]. The fact that small black spruce stems generally tend to allocate more resources to root growth makes them more resistant to post-harvest mortality [33,89,90]. However, a high basal area occupied by dead trees (6 m 2 ha −1 ) was observed on sites where more than 48% of initial basal area was harvested. Partial harvesting increases wind penetration into the managed forest stand, which increases the risk of windthrow [75,91,92] and can then induce mechanical damage to individual stems and result in uprooting or windthrow and, ultimately, tree mortality [15]. Accordingly, these results suggest that during silvicultural treatment, because large trees are more vulnerable to windthrow than small trees [75], especially during the first five years after treatment [24,93], large standing stems (DBH >12.2 cm in our case) should be left behind only if two conditions were met: harvesting intensity is more than 48% of initial basal area; and the conservation of more than 1250 of saplings per hectare, which could act as buffer to the windthrow exposition. In the last two decades, several partial harvesting treatments have been developed in the boreal forest to maintain and reproduce the irregular diameter structure of natural stands in the eastern Canadian boreal forests [23]. Studies conducted in Canadian boreal forests reported that the mortality after silvicultural treatments was proportional to harvesting treatment intensity [17,30,31,94,95]. In the present study, our results suggest that high intensity of partial harvesting should be avoided in poorly drained sites, likely because of their vulnerability to paludification [51,52,80,84], which can be identified by an organic layer thickness on the ground approaching 17 cm or more. That is, partial harvest could be more efficient than conventional management methods to mitigate negative effects of a high water table while limiting paludification [96]. However, knowing that mortality is one of the most complex ecological processes to monitor in post-harvesting forest stand dynamics, long-term monitoring considering other variables (e.g., topography, exposition, wind-speed, competition, and the growth of seedlings) will be necessary to better assess the efficiency of partial harvest systems from either a conservation or wood supply perspective in boreal stands. Moreover, our results indicate that harvesting intensity was one of the main factors involved in both the recruitment into merchantable tree class and the mortality of residual trees.
Finally, in terms of tree growth, our results indicated that although residual tree growth showed a significant effect on decennial stand yield after partial harvesting, this effect seemed marginal. This is likely because in case of tree loss or recruitment of new merchantable trees into the stand, all of the volume of the tree matters rather than only the addition to the pre-existing volume, which is the case of tree growth. Moreover, residual sapling density appeared to be the main factor involved in tree growth dynamics, confirming the fact that an increase in growth depends on an increased residual sapling density [23,54]. This can be explained by the fact that the canopy reopening due to harvesting treatment generally stimulates tree growth. Canopy opening by partial harvesting would increase solar radiation within forest stands, thus favoring the growth of residual trees and the recruitment into the merchantable tree class [85][86][87]. Previous boreal forest studies conducted with different species and partial harvesting treatments in both North American [23] and Fennoscandinavian boreal forests [97] have shown an increase in radial growth some years following partial harvesting. The results of our study showed that regardless of the harvesting intensity, both the CCCC and HARP treatments increased radial growth of residual trees ten years after partial harvesting in all studied sites. However, previous studies showed that the intensity of partial harvesting was relevant in explaining black spruce growth response, especially when harvesting was more than 50% [98,99]. Therefore, our results suggest that both the CCCC and HARP treatments induce a similar growth response in black spruce boreal forest stands.
Partial harvesting can be a viable silvicultural option for black spruce boreal forest stands. However, for this sustainable potential to be realized, it is important to devote more attention to understand the factors of success and failure of different existing partial harvest treatments in these conditions, especially in black spruce forests prone to paludification. For example, in this study, no cases of high mortality (basal area occupied by dead trees) were observed when treatment intensity, i.e., the percentage of harvested basal area, was ≤48%. Fenton et al. (2013) stated that basal area retention needed to maintain at least pre-harvest biodiversity for a large variety of taxa is between 40 and 60%. To ensure an increase in decennial stand yield after harvesting in black spruce forests, partial harvesting should be avoided in sites where the organic layer thickness on the ground approaches 17 cm or more. Second, to achieve successful partial harvesting, a stand selection before cutting should be done based on the structure and habitat conditions criteria. Uneven and low-density forest stands with moderate basal area (<9 m 2 ha −1 ) should be favored. Finally, during silvicultural treatment, large standing stems (DBH >12 cm) should be left on site only when they can be useful as snag legacy for biodiversity, as they are too vulnerable to windthrow.

Conclusions and Implications for Forest Management
Our results evidenced that partial harvesting treatments maintain stand structural attributes for a longer period in the short term in the context of sustainable boreal forest management, and are effective to promote residual tree growth and recruitment into merchantable tree class in black spruce boreal forests. This research shows that depending on pre-harvest stand structures and site conditions, partial harvesting could result in either an increase in post-harvest tree recruitment and growth or a loss of stand basal area because of a high rate of standing tree mortality. An increase in decennial stand yield after harvesting in black spruce forests results mostly from the combination of recruitment into the merchantable tree class, and tree mortality compared to residual tree growth, which appears to play a minor role in stand yield development ten years following partial harvesting treatments.
Because few studies have addressed the factors controlling post-harvest stand dynamics after partial harvesting practices in black spruce boreal forests, before undertaking large-scale partial harvest, it is important to devote more attention to understanding the factors that determine success of partial harvest, in particular treatment intensity (percentage of harvested basal area), the site conditions in terms of drainage class and degree of paludification, and the initial stand characteristics in terms of diametrical structure. This could not only limit forest harvest impacts on biodiversity, but it could also maintain the resilience of the forest by limiting both post-harvest failure of regeneration, and growth and residual tree mortality rate.

Conflicts of Interest:
The authors declare that they have no competing interests.
Appendix A Table A1. The results of step-wise multiple regression models exploring the effect of yield components on decennial stand yield after harvesting in black spruce forests. Elements in bold indicate that the effect of the model on the response variable was significant and had a low AIC. AIC: Akaike's information criterion.

Fixed Effect
Model