Emulating Succession of Boreal Mixedwood Forests in Alberta Using Understory Protection Harvesting

Understory protection harvesting is a form of partial cutting that can be used in aspen (Populus tremuloides Michx.)-dominated stands that have understories of white spruce (Picea glauca (Moench) Voss). This practice involves removing 75% to 85% of the merchantable aspen while minimizing damage to the advance spruce regeneration, in addition to leaving 15% to 25% of the aspen standing to reduce potential windthrow of the spruce understory. In this paper, we summarize results from 18 stands measured 10 to 12 years after understory protection harvest. Diameter growth of spruce increased during the first five years after harvest while height growth increased during the second five-year period (5 to 10 or 7 to 12 years after release). Consistent with other studies, mortality rates of aspen trees ≥7.1 cm DBH (diameter breast height, 1.3 m) averaged 45.0% over the 10–12 year period following harvesting. Spruce mortality averaged 27.5% over the same 10–12 year period. Substantial aspen regeneration was evident across most harvested blocks, with aspen sapling densities 10–12 years from harvest being higher in removal (14,637 stems·ha−1) than in buffer areas (6686 stems·ha−1) and in extraction trails (7654 stems·ha−1). Spruce sapling (>1.3 m height and <4 cm DBH) densities averaged 1140 stems·ha−1 in removal areas at ages 10–12, with these trees likely being present as seedlings at the time of harvest. Mixedwood Growth Model projections indicate merchantable volumes averaging 168 m3·ha−1 (conifer) and 106 m3·ha−1 (deciduous) 70 years from harvest, resulting in MAI (mean annual increment) for this period averaging 2.0 m3·ha−1·y−1 with MAI for a full 150-year rotation of approximately 2.5 m3·ha−1·y−1.


Introduction
Advance regeneration of white spruce (Picea glauca (Moench) Voss) is common in the understory of aspen (Populus tremuloides Michx.)-dominated stands in the boreal forests of western Canada. In addition to regenerating at the same time as aspen following wildfire or other stand-destroying disturbances, spruce may also regenerate over an extended period of time [1,2]. Due to faster growth of young aspen, which regenerates primarily from root suckers, the spruce, which regenerate from seed or through planting, remain in the understory and begin to enter the main canopy after age 60.
In boreal mixedwood stands dominated by aspen and white spruce, growth of aspen peaks between ages 40 and 80, while spruce trees are generally too small to be utilized at this age [3,4]. Harvesting aspen before age 80 is ideal since the incidence of stem decay increases and becomes prevalent in many stands after age 80 [5,6]. Applying partial cutting to remove some or all of the aspen at this age emulates natural succession while allowing harvest of the aspen when it is still sound and mature [7,8]. Protection and release of advanced regeneration can provide for more rapid stocking of stands and shortening of

Experimental Design and Study Sites
The SCUP Permanent Sample Plot (PSP) network was established between 2005 and 2007 with 92 PSPs established in 18 recently harvested SCUP blocks within the Central Mixedwood Natural Subregion of Alberta ( Figure 1). Blocks selected for use in the network all exceeded 120 m × 300 m in size. Within each block, between 2 and 6 PSPs were randomly distributed to capture within-block variability and were remeasured at ages 5-7 and 10-12. For this analysis, the second re-measurement (10-12 years following harvest) was available for the entire SCUP Permanent Sample Plot network. This paper presents results based on measurements collected over 10 to 12 years following strip cut understory protection harvesting in 18 blocks in Alberta. In this paper, we examine spruce response to release, mortality of residual spruce and aspen, and ingress of spruce and aspen, and we explore growth and yield implications of understory protection using the Mixedwood Growth Model (MGM) [16].

Experimental Design and Study Sites
The SCUP Permanent Sample Plot (PSP) network was established between 2005 and 2007 with 92 PSPs established in 18 recently harvested SCUP blocks within the Central Mixedwood Natural Subregion of Alberta ( Figure 1). Blocks selected for use in the network all exceeded 120 m × 300 m in size. Within each block, between 2 and 6 PSPs were randomly distributed to capture within-block variability and were remeasured at ages 5-7 and 10-12. For this analysis, the second re-measurement (10-12 years following harvest) was available for the entire SCUP Permanent Sample Plot network. Sites were located across a climatic gradient with mean annual temperature ranging from −0.3 to 2.5 °C, mean annual precipitation from 432 to 506 mm, and climate moisture index from −0.4 to 3.3 (Table 1). Climatic conditions for the 1981-2010 period were obtained from ClimateNAv6.11 [17] with the Climate Moisture Index (CMI) calculated following methods described by Hogg et al. (2013) [18]. Elevation of these 18 blocks ranges from 362 to 766 m a.s.l. and slope 0% to 20% (85 plots within 10% slope), with all topographic positions represented in the dataset. According to field-determined ecosites [19], 83 plot clusters are situated in ecosite d (low-bush cranberry) and 9 in ecosite e Sites were located across a climatic gradient with mean annual temperature ranging from −0.3 to 2.5 • C, mean annual precipitation from 432 to 506 mm, and climate moisture index from −0.4 to 3.3 (Table 1). Climatic conditions for the 1981-2010 period were obtained from ClimateNAv6.11 [17] with the Climate Moisture Index (CMI) calculated following methods described by Hogg et al. (2013) [18]. Elevation of these 18 blocks ranges from 362 to 766 m a.s.l. and slope 0% to 20% (85 plots within 10% slope), with all topographic positions represented in the dataset. According to field-determined ecosites [19], 83 plot clusters are situated in ecosite d (low-bush cranberry) and 9 in ecosite e (dogwood) of the Central Mixedwood Natural Subregion of Alberta [20]. Soils assessed for a subset of the 60 plots show typical conditions for boreal mixedwoods with mesic to subhygric soil moisture regime and medium to rich soil nutrient regime. Table 1. Information for Strip Cut Understory Protection (SCUP) blocks used in this study. Aspen age (AGEAw) and site index (SI Aw ) were calculated using data from aspen in unharvested buffers, and spruce site index (SI sw ) was obtained using an aspen-to-spruce SI conversion model [21]; both represent tree height at a reference age of 50 years total age. Climate data were from ClimateNA version 6.11 for the 1981-2010 normal period. MAT = mean annual temperature; MAP = mean annual precipitation. Climate Moisture Index was calculated from model outputs of monthly temperature and precipitation following Hogg et al. (2013) [18].

Plot Cluster Design and Measurements
The SCUP PSP network uses a plot cluster design consisting of a set of adjacent plots sampling the three treatment areas, as illustrated in Figure 2. The Extraction Trail is the treatment area within which all trees have been harvested to provide access for harvesting equipment. The Removal Strip is the treatment area subjected to overstory removal (generally aspen) for the purpose of releasing understory conifers (generally white spruce). The Buffer Strip is a "leave" area within which the overstory is retained to reduce windthrow effects on released trees in the Removal Strip.
Plot cluster dimensions vary depending on the width of extraction, removal, and buffer strips. In order to sample sufficient area in each plot, plot length is 8 when buffer strips are 11-20 m wide and 16 m when buffer strips are 5-10 m wide, and this is the same for each subplot (wind buffer (B), retention area (R1/R2), and extraction trail (E) treatments), while sapling plots (S1-S7) are 2 m in length. Subplot widths, including nested sapling plots, vary with the width of each SCUP treatment, as illustrated in Figure 2. Total plot cluster sizes range between 150 m 2 and 463 m 2 . Adjacent to the buffer subplot, an age plot was also established to provide age information for aspen site index determination.
Initial measurements were planned for the year of establishment after harvesting (Measurement 1, establishment year, 2005-2007) followed by Measurement 2 at year 5 to 7 (2010-2014) and Measurement 3 at year 10 to 12 (2015-2019). For measurements conducted after 15 July, one growing season was added for analysis. The first response period, between Measurements 1 and 2, was 5 to 7 years (growing seasons), and the second (between Measurements 2 and 3) was 5 years for all plots, with Measurement 3 occurring 10 to 12 years following harvest.
In each subplot (treatment area), all trees ≥7.1 cm DBH were measured. Saplings (trees ≥1.3 m height and <7.1 cm DBH) were generally measured in corresponding sapling subplots (S1-S7), but in some cases, they were measured in the whole subplot (treatment area). For each measured tree, the following were recorded: species, total height (nearest 0.1 m), diameter at breast height (nearest 0.1 cm), height to live crown (nearest 0.1 m), lean (%), crown class, and up to three condition codes, in order of priority.
period, between Measurements 1 and 2, was 5 to 7 years (growing seasons), and the second (between Measurements 2 and 3) was 5 years for all plots, with Measurement 3 occurring 10 to 12 years following harvest.
In each subplot (treatment area), all trees ≥7.1 cm DBH were measured. Saplings (trees ≥1.3 m height and <7.1 cm DBH) were generally measured in corresponding sapling subplots (S1-S7), but in some cases, they were measured in the whole subplot (treatment area). For each measured tree, the following were recorded: species, total height (nearest 0.1 m), diameter at breast height (nearest 0.1 cm), height to live crown (nearest 0.1 m), lean (%), crown class, and up to three condition codes, in order of priority.

Dataset Screening and Processing
The dataset was screened for errors and growth anomalies prior to analysis following standard procedures for quality checking of permanent sample plot measurements, including outlier analysis and visual inspection of growth trajectories (diameter change, and height change) and height-diameter ratio, consistency in species and origin determination, and tree factor calculation. Missing or anomalous diameter and height values were corrected using height-diameter equations in use in Alberta (Huang 1994 [22]), or through interpolation between values from preceding and succeeding measurements.
Stand growth parameters were aggregated at the subplot level (buffer, removal, extraction) and per hectare by applying appropriate weighting according to the different sizes of trees and sapling subplots. Slenderness was calculated by dividing tree height by diameter at breast height for each tree. Periodic annual diameter and height increment were calculated for the period between measurements and then divided by the number of growing seasons included in the measurement to provide average annual increment values. Survival, mortality, and ingress of each tree were identified (based on tracking tree ID) and used to compute periodic density-related variables (i.e., Mortality rate1-2 =

Dataset Screening and Processing
The dataset was screened for errors and growth anomalies prior to analysis following standard procedures for quality checking of permanent sample plot measurements, including outlier analysis and visual inspection of growth trajectories (diameter change, and height change) and height-diameter ratio, consistency in species and origin determination, and tree factor calculation. Missing or anomalous diameter and height values were corrected using height-diameter equations in use in Alberta (Huang 1994 [22]), or through interpolation between values from preceding and succeeding measurements.
Stand growth parameters were aggregated at the subplot level (buffer, removal, extraction) and per hectare by applying appropriate weighting according to the different sizes of trees and sapling subplots. Slenderness was calculated by dividing tree height by diameter at breast height for each tree. Periodic annual diameter and height increment were calculated for the period between measurements and then divided by the number of growing seasons included in the measurement to provide average annual increment values. Survival, mortality, and ingress of each tree were identified (based on tracking tree ID) and used to compute periodic density-related variables (i.e., Mortality rate 1-2 = (Mortality 1-2 /Density 1 ) × 100; Ingress rate 1-3 = (Ingress 1-3 /Density 3 ) × 100; Survival rate 1-2 = (Density 2 /Density 1 ) × 100).
Aspen top height information was collected in aspen age plots (100 m 2 with one aspen tree selected or 300 m 2 where three aspen trees were selected for measurement) located in the buffer of each plot cluster. The most recent measurement was used to calculate aspen site index for each plot cluster, at a reference age of 50 years total age, using Huang et al.  [23] top height models. Age corrections [23] were applied to calculate total age from ring counts at breast height for each tree. Where more than one tree was available in an age plot, total age and height were averaged before calculating site index (SI). Results show aspen site index ranging from 17.7 m to 26.6 m ( Figure 3), with the majority of plots falling between 19 and 22 m. Three blocks were very productive (SI > 26 m).   Due to a lack of measured spruce in age plots, as well as concerns about using spruce that had spent decades growing in the understory, an SI conversion model (aspen-to-white spruce) was used to calculate spruce SI [21]. Mean spruce SI varied between 14.6 m and 26.5 m across SCUP blocks, with 5 blocks having spruce SI above 23.5 m ( Figure 3).

Assessing Growth Response to Treatments
Growth, mortality, and ingress of spruce and aspen were analyzed using data for all measured trees above the common tagging limit (≥7.1 cm in DBH) in removal and buffer areas, while on extraction trails, trees were, in most cases, smaller than the common tagging limit and, since they were not remeasured, were not included in the analysis of tree growth responses. Analysis was performed in two ways: (a) for the first response period between Measurement 1 and Measurement 2, and for the second response period between Measurements 2 and 3; and (b) for the entire 10-year response period between Measurements 1 and 3. The same subset of trees was used for both (a) and (b) analyses. Buffer subplots from 5 blocks (7012, 27131, 19191, 29691, 16751) were not measured at Measurement 2 and therefore could not be used in this analysis. In addition, plot 5 from block 36551 was excluded from analysis because it lacked complete Measurement 2 data. As a result, 61 plots were available for analysis of buffer areas and 91 plots were available for both removal and extraction treatments (by species).
Saplings, defined as deciduous or coniferous trees ≥1.3 m height and <7.1 cm DBH, were analyzed separately. Plots which had incomplete or no sapling measurements were excluded from analysis, resulting in a total of 58 plots available for the buffer areas and 91 plots available for both removal and extraction treatments. Aspen saplings were not measured at Measurement 1. Sapling ingress analysis was limited to tagged and measured saplings, where the lower tagging limit was 1.3 m in height and the upper tagging limit was 7.1 cm in DBH.
To analyze spruce response to release, we tested the effects of treatments (removal and buffer) on spruce diameter and height increment at both the plot and tree levels. Spruce growth response and aspen and spruce mortality rates were examined using a linear mixed-effect model fitted using the lme function from the nlme package [24] in R statistical software [25]. The hierarchical data structure was represented by making treatment a fixed effect and block a random effect, while for tree level analysis, plot was included as random effect. For both plot and tree level analysis, spruce SI and initial size were included in the model as covariates. For all statistical analyses, a significance level of α = 0.05 was used. Homogeneity of slopes between covariates and treatment was confirmed in order to meet assumptions for ANCOVA. Assumptions of homogeneity of variances and normality of residuals were examined using diagnostic plots of the residuals. The response variables were transformed with natural logarithms to meet assumptions of normality and homoscedasticity, where necessary. To calculate marginal and conditional coefficients of determination (R 2 ) for mixed-effects models, MuMIn package [26] was used.

Yield Projections
Yield predictions were computed for the entire network of 18 SCUP blocks (91 plot clusters) using the Mixedwood Growth Model (MGM21) [16,27]. Full tree lists were created for each strip (extraction, removal, and buffer) in each plot cluster based on the latest measurements (Measurement 3) collected 10 to 12 years after the SCUP harvest (2016-2019) and which included all trees taller than 1.3 m height. Aspen and spruce SI was estimated at the block level for MGM, while a default SI of 15 m was used for black spruce.
Starting age for growth simulations was set to 0 for all subplots (buffer, removal, extraction) and simulations were run to age 200 (i.e., 200 years since Measurement 3). All MGM simulations used the GYPSY site index curves for Alberta [24], and the Alberta Central Mixedwood taper equations for volume calculations [22]. Merchantable volumes were calculated using minimum DBH of 13.51 cm and 13.67 cm for spruce and aspen, MGM projections were made with MGM's Light/Adjacency Submodel (MSLight) enabled. The Light/Adjacency Submodel reduces the growth of trees based on the shading effect from adjacent strata (most notably from tall trees in the buffer). To run MGM's Light/Adjacency Submodel, subplot boundaries were created with georeferenced plot centers, corridor azimuths, and plot dimensions from the SCUP database. Since the SCUP measurement plot design only includes one buffer (Figure 2), the buffer polygon was replicated next to the R2 subplot in order to properly simulate shading effects from both sides of the harvest pattern in the SCUP treatment. MGM defaults were used to define the beginning of the growing season (Julian Day 106), the end of the growing season (Julian Day 274), and the diffuse radiation fraction (0.5).

Spruce Growth Response to Release
Size and growth of measured spruce (≥7.1 cm DBH) within buffer and removal subplots are summarized in Table 2. Quadratic mean diameter (QMD, the diameter of the tree of average basal area) and mean height were larger in buffer than removal plots, but differences were small (ranging from 0.5 to 2.0 cm for QMD and 0.6 to 1.1 m for height). Across the three measurements, mean diameter increased while mean height decreased in both buffer and removal areas. This is the result of simultaneous tree mortality, damage, and ingress (trees passing the tagging limit), which reflects the irregular (structured) stand conditions created after the SCUP harvest. Slenderness (ratio of height to DBH) also changed as a reflection of declining average height and increasing diameter. Spruce density increased, which, in combination with increased DBH, led to an increase in basal area from 8 to 11.7 m 2 ·ha −1 and from 5.7 to 7.8 m 2 ·ha −1 for buffer and removal areas, respectively. Table 2. Treatment means and standard deviations for spruce size by measurement. Tagging limit is ≥7.1 cm DBH for all measurements. Note: Plots with no spruce trees were excluded from this analysis; the period between Measurements 1 and 2 ranged from 5 to 7 years; between Measurements 2 and 3, it was 5 years. Overall, between Measurements 1 and 3, the period ranges from 10 to 12 years. Diameter increment and height increment were analyzed for all trees with at least two measurements (Table 3). Sample size varied between measurements due to mortality and ingress. Mean annual diameter increment during the first response period was 0.44-0.47 cm·y −1 , and 0.58-0.60 cm·y −1 during the second period, averaging 0.50-0.54 cm·y −1 over the full 10-12-year period. Mean annual height increment at the plot level was 0.07-0.09 m·y −1 during the first 5-7 years and increased to 0.24-0.27 m·y −1 during the next 5 years, while for the whole period, height increment ranged between 0.14 m·y −1 and 0.16 m·y −1 . Table 3. Summary of spruce growth (trees ≥7.1 cm DBH) response at plot level and at tree level. Id = diameter increment. Ih = height increment. Note: Plots with no spruce trees were excluded from this analysis; the first response period between Measurements 1 and 2 ranged from 5 to 7 years; the second between Measurements 2 and 3 was 5 years. The overall response period between Measurements 1 and 3 ranged from 10 to 12 years. We compared growth in removal to growth in buffer areas using spruce SI and initial size as covariates in the mixed model ANOVA to account for effects of site and tree size. None of the models tested for each period and response (diameter and height increment) revealed significant differences between removal and buffer for spruce growth (Table 4). SI was only significant for height increment during the first period at the plot level, while initial size was significant in all models except for the plot-level diameter increment in the second response period (which was only marginally non-significant). Table 4. Results from linear mixed-effect models for diameter and height increment of white spruce by measurement period. R 2 m is marginal coefficient of determination; R 2 c is the conditional coefficient of determination. Id is diameter increment, Ih is height increment. Bold numbers indicate a statistically significant difference (p-value less than 0.05).

Growth of Spruce and Aspen Saplings
Annual diameter increment and height increment were analyzed using data for all tagged and remeasured saplings (Figure 4). Due to a lack of suitable data for aspen at Measurement 1, analysis of aspen growth was limited to the second growth period. Diameter increment of aspen saplings during the second response period ranged between 0.14 and 0.16 cm·y −1 , and height increment ranged between 0.23 and 0.25 m·y −1 with similar growth responses for aspen in buffer, removal, and extraction strata. Due to the small numbers of spruce saplings present, we were unable to analyze the growth of spruce saplings in the extraction trails. Spruce sapling diameter increment averaged 0.26 and 0.30 cm·y −1 for buffer and removal, respectively, with second period values averaging 0.22 and 0.28 cm·y −1 for buffer and removal, respectively. Mean annual height increment of the spruce saplings was 0.12-0.15 m·y −1 during the first period, and 0.14-0.21 m·y −1 during the second period. Height increment of spruce saplings does appear to be increasing with time since removal.
Forests 2022, 13, x FOR PEER REVIEW 11 of 20 0.28 cm⋅y −1 for buffer and removal, respectively. Mean annual height increment of the spruce saplings was 0.12-0.15 m⋅y −1 during the first period, and 0.14-0.21 m⋅y −1 during the second period. Height increment of spruce saplings does appear to be increasing with time since removal.

Mortality of Spruce and Aspen Trees
The common tagging limit of 7.1 cm DBH was applied in the analysis of mortality, as that was the tagging limit applied across all measurements. Results (Table 5) show that mean stand density in the buffer at initial measurements was similar for both species (363 stems⋅ha −1 of spruce and 392 stems⋅ha −1 of aspen). In removal areas, density of both species was lower (268 stems⋅ha −1 for spruce and 39 stems⋅ha −1 of residual aspen trees, on average), reflecting the expected effect of the partial harvest. However, high variability in stand density across sample plots was found for both species.

Mortality of Spruce and Aspen Trees
The common tagging limit of 7.1 cm DBH was applied in the analysis of mortality, as that was the tagging limit applied across all measurements. Results (Table 5) show that mean stand density in the buffer at initial measurements was similar for both species (363 stems·ha −1 of spruce and 392 stems·ha −1 of aspen). In removal areas, density of both species was lower (268 stems·ha −1 for spruce and 39 stems·ha −1 of residual aspen trees, on average), reflecting the expected effect of the partial harvest. However, high variability in stand density across sample plots was found for both species.
Stand density in buffer areas increased with time for spruce from 363 stems·ha −1 to 479 stems·ha −1 as ingress was higher than mortality (Table 6). In contrast, aspen density decreased from 392 stems·ha −1 to 277 stems·ha −1 because mortality was almost three times higher than ingress over the entire observed period. As expected, ingress caused an increase in stand density for both species in removal areas. However, while ingress of spruce (92-95 stems·ha −1 ) was fairly constant over the first two periods, a delay in aspen ingress was evident. Mortality was highest for aspen in buffer areas as old aspen trees were exposed to extreme weather conditions after SCUP harvest. Mortality rate of aspen was similar during the first two periods for both treatments, ranging between 24% and 31%, with cumulative mortality over the whole period of 43% to 47%. For spruce, mortality was higher during the first period, averaging 50 stems·ha −1 and 77 stems·ha −1 for buffer and removal, respectively, but declined to 21-23 stems·ha −1 during the next five years. Overall, mortality rate was higher for aspen than for spruce. Spruce mortality rate decreased from 21-24% in the first period to 5-9% in the second period for both treatments. Mixed-model ANOVA, conducted separately for each species and period, did not reveal significant differences in mortality rate between buffer and removal areas for either species.  Table 6. Summary of mortality and ingress for aspen and spruce trees (≥7.1 cm DBH) over the three measurement periods. Note: There were different lengths of the first response period (between Measurements 1 and 2) of 5-7 years for different blocks, and the same length of the second response period (between Measurements 2 and 3) of 5 years across all blocks and plots. Subplots with no recorded trees were included in averaging (as zero density), but were not used in computation of the mortality rate.  Figure 5 shows sapling density per hectare by species. Due to a lack of suitable data for aspen at Measurement 1, analysis of aspen density was limited to Measurements 2 and 3. Highest numbers of saplings were found in removal areas, then in extraction, and the lowest in buffer areas. However, trees crossing tagging limits during the period affect density-related variables. Between Measurements 2 and 3, there was an increase in aspen sapling density in removal areas from 12,334 to 14,637 stems·ha −1 , while aspen sapling density decreased in the buffer from 7771 to 6686 stems·ha −1 , and decreased in the extraction areas from 9697 to 7654 stems·ha −1 . Based on these densities, aspen ingress was higher than mortality in removal areas, while mortality was higher than ingress in buffers and extraction trails. The aspen ingress rate was the same in removal and extraction areas, and these were about two times higher than in buffer subplots.

Density of Aspen and Spruce Saplings
Forests 2022, 13, x FOR PEER REVIEW 13 of 20 Figure 5. Sapling density and ingress by species and by treatment. Lower tagging limit is 1.3 m height and upper tagging limit is 7.1 cm DBH for both species. Means and standard errors (error bars) are presented. n is sample size (number of plots). Note: First response period (between Measurement 1 and 2) ranged from 5 to 7 years; for the second response period (between Measurements 2 and 3), the time interval was 5 years; and between Measurements 1 and 3, the overall time interval was between 10 and 12 years. Subplots with no recorded desired trees were accounted in averaging (as zero density), but were not used in the computation of ingress rate. Due to a lack of data for Measurement 1, aspen density at Measurement 1 and ingress for the first response period could not be calculated.
In the buffer areas, spruce saplings had similar density over the first two measurements with 427-440 stems⋅ha −1 and decreased at Measurement 3 to 346 stems⋅ha −1 ( Figure 5). In the removal areas, spruce sapling density increased from 449 to 624 stems⋅ha −1 from Measurement 1 to Measurement 2 and increased further to 1140 stems⋅ha −1 at Measurement 3. For the extraction trails, spruce sapling density increased from 0 to 19 stems⋅ha −1 from Measurement 1 to Measurement 2 and increased further to 55 stems⋅ha −1 at Measurement 3. Ingress of spruce was higher in the removal area than in other areas and increased over time, while in the buffer, it showed a decrease. Spruce ingress rate in both treatments decreased from 46% to 37% and from 34% to 23% for removal and buffer, respectively. Figure 5. Sapling density and ingress by species and by treatment. Lower tagging limit is 1.3 m height and upper tagging limit is 7.1 cm DBH for both species. Means and standard errors (error bars) are presented. n is sample size (number of plots). Note: First response period (between Measurement 1 and 2) ranged from 5 to 7 years; for the second response period (between Measurements 2 and 3), the time interval was 5 years; and between Measurements 1 and 3, the overall time interval was between 10 and 12 years. Subplots with no recorded desired trees were accounted in averaging (as zero density), but were not used in the computation of ingress rate. Due to a lack of data for Measurement 1, aspen density at Measurement 1 and ingress for the first response period could not be calculated.
In the buffer areas, spruce saplings had similar density over the first two measurements with 427-440 stems·ha −1 and decreased at Measurement 3 to 346 stems·ha −1 ( Figure 5). In the removal areas, spruce sapling density increased from 449 to 624 stems·ha −1 from Measurement 1 to Measurement 2 and increased further to 1140 stems·ha −1 at Measurement 3. For the extraction trails, spruce sapling density increased from 0 to 19 stems·ha −1 from Measurement 1 to Measurement 2 and increased further to 55 stems·ha −1 at Measurement 3. Ingress of spruce was higher in the removal area than in other areas and increased over time, while in the buffer, it showed a decrease. Spruce ingress rate in both treatments decreased from 46% to 37% and from 34% to 23% for removal and buffer, respectively.

Growth and Yield Implications
Stratum and block (whole plot) level trends in conifer, deciduous, and stand level merchantable volume estimated using MGM21 are shown in Figure 6. Predictions were initiated for each plot at age 0 using Measurement 3 tree lists (10 years following SCUP harvest). Deciduous dominated in the extraction trails, and conifers dominated in the removal and buffer areas. Deciduous volume began to decline after 20 years in buffer areas for all plots, while conifer volume showed a steep rise in the first 50 years and continued to increase up to 85 years. In removal areas, deciduous volume increased slowly up to 72 years, and conifer volume increased gradually and reached culmination values at 96 years. Deciduous merchantable volume for the whole plot (weighted by treatment areas) decreased after 60 years. Conifer merchantable volume for the whole plot increased over time up to 93 years; however, the rate of increase was low after 60 years. Similarly, total stand volume at the whole plot level (all strata included) reached a maximum of 293 m 3 ·ha −1 at 81 years. In year 60 (Figure 7), deciduous volume was higher for extraction (173 m 3 ·ha −1 ) than for buffer (124 m 3 ·ha −1 ) or removal (58 m 3 ·ha −1 ) areas. Coniferous volume was higher in buffer (240 m 3 ·ha −1 ) than in removal areas (194 m 3 ·ha −1 ) and was negligible in extraction areas (5 m 3 ·ha −1 ). Merchantable total (stand) volume (weighted by treatment areas) in year 60 ( Figure 7  On the X-axis, age represents years after initiation of the simulations (10 to 12 years after SCUP harvest) for all treatments.

Spruce Growth Response to Release
Trees >7.1 cm in DBH showed only moderate increases in diameter and height of spruce following understory protection harvest. A lack of differences between spruce in buffer and removal areas is not surprising, given that spruce in the buffers also experienced improved growing conditions due to the removal of trees in the adjacent removal areas. Height increment during the second 5-year response period increased approximately threefold compared to the first 5-7 year response period after harvesting, while that increase was about 25% for diameter increment. Krebs (2016) [13] reported a 3-5-year delay in height growth response of spruce advanced regeneration after removal of aspen canopies but an immediate diameter response, while Kneeshaw et al. (2002) [28] observed a 2-3-year delay in height growth response. Based on a meta-analysis using 17 data sources, Man and Greenway (2004) [12] found that for trees ranging in initial height between 1 m and 19 m, spruce height growth increases were largest for 8 m tall trees, and diameter responses were immediate. Many other partial cutting studies show increased growth of spruce following reductions in the amount of aspen, with growth generally decreasing with increasing levels of canopy retention [14,[29][30][31]. Smith et al. (2016) [32] report stronger responses of suppressed trees compared to codominant and dominant trees, and, similar to Krebs (2016) [13], report that conifer neighbors were the primary competitors affecting spruce growth. Studies also indicate that aspen growth increases following partial canopy removal, with growth increasing mostly at the lowest levels of retention [8,33].

Spruce Growth Response to Release
Trees >7.1 cm in DBH showed only moderate increases in diameter and height of spruce following understory protection harvest. A lack of differences between spruce in buffer and removal areas is not surprising, given that spruce in the buffers also experienced improved growing conditions due to the removal of trees in the adjacent removal areas. Height increment during the second 5-year response period increased approximately threefold compared to the first 5-7 year response period after harvesting, while that increase was about 25% for diameter increment. Krebs (2016) [13] reported a 3-5-year delay in height growth response of spruce advanced regeneration after removal of aspen canopies but an immediate diameter response, while Kneeshaw et al. (2002) [28] observed a 2-3-year delay in height growth response. Based on a meta-analysis using 17 data sources, Man and Greenway (2004) [12] found that for trees ranging in initial height between 1 m and 19 m, spruce height growth increases were largest for 8 m tall trees, and diameter responses were immediate. Many other partial cutting studies show increased growth of spruce following reductions in the amount of aspen, with growth generally decreasing with increasing levels of canopy retention [14,[29][30][31]. Smith et al. (2016) [32] report stronger responses of suppressed trees compared to codominant and dominant trees, and, similar to Krebs (2016) [13], report that conifer neighbors were the primary competitors affecting spruce growth. Studies also indicate that aspen growth increases following partial canopy removal, with growth increasing mostly at the lowest levels of retention [8,33].

Mortality and Ingress of Aspen and Spruce Trees
Windthrow is a common cause of loss of trees over 7.5 m tall following understory protection but was not documented in this study. MacIsaac and Krygier (2017) [9] recommend that buffers should be no more than 75 m apart (i.e., not exceeding 2.5 times the height of dominant/codominant aspen), with preference for narrow buffers being located less than 35 m apart. Given that distances between buffers for the plots used in this study were generally less than 25 m, only limited windthrow might be expected. Solarik et al. (2012) [34] observed 30% mortality of aspen during the first 5 years and 50% mortality of aspen over the first 10 years following removal of 80% of the overstory canopy at the EMEND (Ecosystem Management Emulating Natural Disturbance) site in northwestern Alberta, while mortality in the unharvested was 10% and 20% over these time periods. Xing et al. (2016) [31] found that mortality of aspen increased to 6.9% during the 5 to 10 year period after removals with 20% overstory retention compared to unharvested stands, where it was 2.4%. In contrast,   [35] observed that mortality of aspen trees >10 cm DBH was lower under 33% retention than in unharvested areas. In general, mortality rates and windthrow appear to decline with time after harvesting [9,31,32]. Our results indicate aspen mortality rates of 27.5%, 26.6%, and 45.0% over the first, second, and overall response periods, which are consistent with ranges reported by these other studies.
Man and Greenway (2004) [12] report that release had little influence on survival/mortality of white spruce. Xing et al. (2016) [30] report that spruce mortality following 20% retention was approximately 5% in the 1-5-year period compared to approximately 2.5% in non-harvested areas and was lower in the 20% retention (1.8%) than in the non-harvested (2.1%) in the 5-10 year period. Solarik et al. (2012) [34] report cumulative 10-year mortality of spruce residual trees averaging 30%, with taller spruce with live crown percent below 50% having the highest mortality rate. Our results indicate spruce mortality rates of 22.4%, 7.0%, and 27.5% over the first, second, and overall response periods, respectively, following understory protection harvest, which are substantially larger than values reported by Xing et al. (2016) [30] but are consistent with Solarik et al. (2012) [34].

Density of Aspen and Spruce Saplings
While loss of apical dominance promotes suckering of aspen [36], when more than 35% of the original canopy is retained following partial cutting there is usually little improvement in suckering or sucker survival due to the effects of shade from the overstory canopy [33,34,37,38]. Our observations of aspen sucker densities being higher in removal (14,637 stems·ha −1 ) than in buffer areas (7654 stems·ha −1 ) at age 10 are consistent with the suggestion that hormonal control is likely the dominant factor controlling aspen suckering [33].
Conifer sapling recruitment generally increases with time following harvesting and tends to be higher in partial cuts than in clearcuts [34], presumably due to the presence of seed sources and reduced amounts of competition. We observed recruitment of spruce saplings in the removal areas with the numbers increasing by 691 stems·ha −1 over 10 years to reach a total of 1140 stems·ha −1 taller than 1.3 m at age 10. Very few spruce saplings were observed in extraction trails at this age, due to limited amounts of regeneration and slow growth of natural spruce regeneration on the extraction trails.

Growth and Yield Implications of Understory Protection
Strip cut understory protection using a two-stage harvest design can optimize the yield of boreal mixedwoods and overcome challenges caused by different timing of maximum mean annual increment (MAI) for aspen and white spruce [4]. Seventy years after the first harvest, block-level conifer volumes projected by MGM averaged 168 m 3 ·ha −1 and deciduous volumes averaged 106 m 3 ·ha −1 , with MAI for this period averaging 4.5 m 3 ·ha −1 ·y −1 . Given that aspen volume removed in the initial harvest was likely to have been approximately 200 m 3 ·ha −1 , this would result in a total aspen MAI of 2.0 m 3 ·ha −1 ·y −1 over 150 years, resulting in total stand MAI over 150 years of 2.5 m 3 ·ha −1 ·y −1 .
When spruce regenerates naturally in these stands, as was the case in our study, costs of planting spruce are avoided. Comparing understory protection to growing pure aspen stands or pure spruce plantations indicates that, while all three approaches can provide similar levels of mean annual increment (MAI) over a full rotation, the pure stand approaches yield primarily individual species while understory protection results in yield of both species [3].

Study Limitations
A lack of paired non-harvested control plots precluded direct comparison of outcomes from understory protection to non-harvested stands. Establishing additional control plots in similar stands near each SCUP block could improve this study and future analysis and outcomes. Changes in measurement protocols, a lack of full measurement of regeneration (missing plots), and failure to tag aspen saplings at the time of the initial measurement limited analysis and interpretation. Follow-up with a fully replicated designed experiment that includes a non-harvested control, an understory protection treatment (perhaps including both strip and dispersed aspen retention), and a single-pass clearcut, with measurements initiated prior to harvest, would provide valuable additional information on the merits of understory protection relative to other options. To provide data on the merits of understory protection under a changing climate, we recommend supplemental collection of data on soil and understory microclimate linked to annual growth responses of trees in all layers to treatment and climate.

Conclusions
We found no differences in spruce growth between removal and adjacent buffer areas 10-12 years after harvest. However, height increment during the second 5-year response period increased approximately threefold compared to the first period after the SCUP harvest, while that increase was about 25% for the diameter increment of spruce trees. Aspen saplings were more abundant in removal areas (14,637 stems·ha −1 ) than in extraction trails (7654 stems·ha −1 ) at Measurement 3, although there was substantial variation in aspen densities between blocks. Aspen sapling density in buffer areas averaged 6686 stems·ha −1 . At age 10-12 (Measurement 3), spruce sapling density averaged 1140 stems·ha −1 in removal areas, 346 stems·ha −1 the buffers, and 55 stems·ha −1 in extraction trails. Overall, the mortality rate was higher for aspen than for spruce. Ingress rates for aspen saplings were the same in removal areas and extraction trails, which were about two times higher than in buffer areas. Ingress rates for spruce saplings in removal and extraction sub-plots decreased slightly over time.
Seventy years after understory protection harvest, MGM estimated yields of 105 m 3 ·ha −1 for deciduous and 168 m 3 ·ha −1 for conifer, in addition to an estimated 200 m 3 ·ha −1 removed during the first entry of the understory protection harvest. These results indicate the potential for understory protection to increase total stand yields through the removal of aspen when it is mature, followed by later removal of spruce and any merchantable aspen in a second pass approximately 70 years later.
The results shown here demonstrate the potential to apply understory protection harvesting in aspen-dominated stands that have a white spruce understory and suggest an opportunity to increase aspen yields through the mid-rotation harvesting of aspen, while also increasing spruce yields in these stands. As well as providing an opportunity to provide mid-term yields of spruce through application of understory protection harvesting in existing stands, the planned application of understory protection when growing aspenspruce mixtures provides an opportunity to emulate succession in these stands, enables the use of aspen as a nurse crop during the establishment of spruce, and could be used effectively in the maintenance of a diversity of stand types in a forest management unit.
Author Contributions: I.B. and P.C., analysis and writing of the manuscript; S.M., secured funding for and coordinated measurements, secured funding to support data analysis; S.M. and B.R., contributed to writing of the final manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: Funding to support data analysis, preparation and publication of this paper was provided by the Forest Growth Organization of Western Canada (FGrOW) Mixedwood Practices Team. Data Availability Statement: Data from this study may be made available via application to FGrOW. Any request for data should include a detailed explanation of the purpose of the request and plans for use of these data.