Simulating the Effects of Intensifying Silviculture on Desired Species Yields across a Broad Environmental Gradient

In the past two decades, forest management has undergone major paradigm shifts that are challenging the current forest modelling architecture. New silvicultural systems, guidelines for natural disturbance emulation, a desire to enhance structural complexity, major advances in successional theory, and climate change have all highlighted the limitations of current empirical models in covering this range of conditions. Mechanistic models, which focus on modelling underlying ecological processes rather than specific forest conditions, have the potential to meet these new paradigm shifts in a consistent framework, thereby streamlining the planning process. Here we use the NEBIE (a silvicultural intervention scale that classifies management intensities as natural, extensive, basic, intensive, and elite) plot network, from across Ontario, Canada, to examine the applicability of a mechanistic model, ZELIG-CFS (a version of the ZELIG tree growth model developed by the Canadian Forest Service), to simulate yields and species compositions. As silvicultural intensity increased, overall yield generally increased. Species compositions met the desired outcomes when specific silvicultural treatments were implemented and otherwise generally moved from more shade-intolerant to more shade-tolerant species through time. Our results indicated that a mechanistic model can simulate complex stands across a range of forest types and silvicultural systems while accounting for climate change. Finally, we highlight the need to improve the modelling of regeneration processes in ZELIG-CFS to better represent regeneration dynamics in plantations. While fine-tuning is needed, mechanistic models present an option to incorporate adaptive complexity into modelling forest management outcomes.


Introduction
Forest management has undergone several major paradigm shifts in the past two decades. While forest managers continue to aim to maintain or improve productivity by intensifying silviculture [1][2][3], guidelines such as those advocating the emulation of natural disturbance patterns at landscape scale [4][5][6] and a desire to enhance forest complexity at stand scale [7,8] have created a need to revise forest productivity and stand dynamics models. These needs, coupled with improvements in successional theory [9,10], computing power, increasing availability of long-term forest inventories [11], and the ongoing acceleration of climate change effects on forest dynamics [12] present both challenges and opportunities to forest management modelling, especially as global wood supplies become plete [9]. Few long-term experiments address the dynamics of mixed-species planted forests after treatments such as releasing conifer regeneration from overtopping hardwood cover or thinning of hardwood species [47]. Forest managers in boreal and northern temperate forests of Canada commonly predict stand development of planted forests using classical succession theory [17,18], based on which plant communities develop along a linear pathway from initial occupancy by shade-intolerant pioneers to late successional shade-tolerant species [48]. However, recent developments in succession theory describe forest stands as moving through multiple potential pathways, depending on pre-and post-disturbance species compositions and traits [9,49,50]. These multiple potential developmental pathways are difficult to project using empirically-based models currently available to forest managers in Canada [21,22,25], which are traditionally calibrated for even-aged monospecific or hyper-dominant stands. The number of possible pathways increases further when type, frequency, duration, and intensity of disturbance or the suppression of these disturbances are considered. For instance, in eastern boreal forests, where fire cycles are long and climates are relatively humid, a combination of fire suppression and harvesting without applying any silvicultural interventions has greatly increased the proportion of balsam fir (Abies balsamea [L.] Mill), a typically late-successional shadetolerant conifer, in young stands relative to pre-management conditions [51]. However, in central boreal forests, where fire cycles are shorter and climates are drier, harvesting without further silvicultural intervention has led to a proliferation of mixed wood conditions due to the rapid ingress of shade-intolerant hardwood species that typically existed in small proportions before harvesting [52,53]. Effective forest management simulation models must effectively incorporate these divergent pathways but even with the increasing availability of data from permanent sample plot networks, many forest conditions remain undersampled [11].
Climate change poses a real challenge for forest ecosystems and for forest management planning [12,54]. Rapid and non-linear changes in environmental conditions also challenge the use of traditional growth and yield models and tables, which were developed based on historical climates [20]. For example, reductions in tree longevity due to climate change would necessitate a decrease in rotation age that is not captured by yield-based models [55,56]. Climate change can affect stand dynamics if certain species are favoured over others and thus further disrupt already complex successional pathways, potentially reducing the yield of desired species and possibly triggering biome shifts [57][58][59]. Given the dearth of longer-term operational studies designed to determine the effects of climate change on forest productivity and stand dynamics [60], forest managers need to rely on models to predict climate change effects. Such models need to explicitly capture the effects of a range of climate conditions (i.e., increases in atmospheric CO 2 and annual temperatures and changes in precipitation regimes).
Finally, robust stand dynamic models also need to be applicable to a range of silvicultural systems (i.e., the process of tending, removing, and replacing crop trees; [61]), including even-and uneven-age silvicultural systems [62]. For instance, in forests dominated by jack pine (Pinus banksiana Lamb.), clearcuts, high mineral soil exposure, and high seeding or replanting densities are typically needed to regenerate to pre-harvest conditions [63]. Even then, vegetation management to control the ingress of shade-intolerant hardwoods may be needed if initial site preparation is insufficient to favour jack pine regeneration [64,65]. In temperate forests, uneven-aged selection systems are applied to regenerate tolerant hardwoods [66]. These differences in silvicultural systems further increase complexity in forest modelling.
While the development of empirical models for each paradigm shift is possible, a modelling platform that supports all paradigm shifts would be substantially more relevant [17,18]. Gap models are useful to fill this need: their emphasis on modelling underlying ecological processes means they can be used to understand dynamics for a variety of forest types, silvicultural systems, and climate scenarios [21,62,67,68]. Recent studies have demonstrated that various gap models can account for partial cuts in boreal mixedwoods as well as natural disturbances in temperate forests [26,27].
Here we use an experimental plot network encompassing much of the range of managed forests in Ontario, Canada, to examine how increasing silviculture intensity can maintain high commercial yields while managing for complexity. The experimental platform is implemented across four fire regimes and encompasses four silvicultural intensities as well as monitoring plots in undisturbed forest conditions [37]. To predict future yield and species composition, we used the ZELIG-CFS forest gap model (a version of the ZELIG tree growth model adapted for use in Canada by the Canadian Forest Service) [26,68] to simulate stand dynamics for 70 years after harvesting, which represents rotation age for most of the study sites. We further examined whether climate scenarios could be included in the ZELIG-CFS modelling framework by examining a subset of study sites expected to undergo the most significant changes in climate.

Study Site and Plot Measurements
The NEBIE plot network was used to estimate the effects of silvicultural intensities and fire regimes on managed forest function and productivity. The plot network is fully replicated in six locations across the province of Ontario, Canada: Sioux Lookout, Dryden, Timmins, Kapuskasing, Petawawa, and North Bay (in ascending order of fire return interval). The locations of the sites are presented in Figure S1. Sioux Lookout and Kapuskasing are typical boreal conifer sites, Dryden and Timmins are boreal mixedwoods, and Petawawa and North Bay are temperate forests. At each site, a randomized complete block design was established consisting of four blocks with complete replication of all silvicultural interventions in each block, except for North Bay and Timmins that had only three blocks remeasured consistently post-treatment. Silvicultural intensities in increasing order of extent were: natural (i.e., left unharvested), extensive, basic, intensive, and elite. Exact interventions varied by site but generally stocking of desired species regeneration was 40% for extensive, 60% for basic, and 80% for intensive and elite treatments. Desired species were: jack pine in Sioux Lookout, white spruce (Picea glauca (Moench) Voss) and black spruce (Picea mariana [Mill.] B.S.P.) in Kapuskasing, white spruce in both Dryden and Timmins, white pine in Petawawa, and sugar maple (Acer saccharum Marsh.) with an increasing component of yellow birch (Betula alleghaniensis Britt.) in North Bay. Complete information about each treatment at each site as well as pre-harvest stand conditions and site information is presented in Bell et al. (2017).
Measurements from before harvesting and second, fifth-, and tenth-year post-harvesting are available from the plot network. In each treatment and block combination, four growth and yield permanent sample plots were established following Ontario's Forest Growth and Yield Program specifications [69,70]. Sample plots were circular with an area of 400 m 2 . All trees more than 2.5 cm in diameter at breast height (DBH) were located, tagged, and measured; however, in the 10th year remeasurement, new recruits were only measured in the northeast quadrant in Sioux Lookout, Dryden, Timmins, and Kapuskasing due to high stem densities. To properly match the recruits to the veteran trees (measured on the full plot area of 400 m 2 ), the tree list for recruits was replicated three times (to scale one quadrant measurements to the full plot). All veteran trees (i.e., those that were present pre-harvest) and all trees in the natural plots were measured in the whole plot area and were not included in the replication list. Finally, 10 stocking plots measuring 2 × 2 m were used to assess the stocking of seedlings in each sample plot. In total 110 experimental units were assessed (five treatments per block, four blocks per site, with three blocks in Timmins and North Bay, and six sites total), including 440 growth and yield permanent sample plots and 4440 stocking plots ( Figure S2).

Forest Simulation Model
ZELIG-CFS is a gap model that simulates individual tree growth, recruitment, and mortality, and has proven useful for estimating successional trajectories of northern forests and simulating the effect of silvicultural interventions [26,68]. It originates from the model ZELIG [71] but, based on a previous study [72], several modifications were introduced including the modelling of crown interaction effects, survival rate, and regeneration [68]. The ZELIG-CFS model is initialized from a list of tree variables, including the tag number of the tree, species, and DBH and requires a series of site-specific inputs influencing the simulated forest dynamics: plot location, elevation, soil texture type (used to calculate wilting point and field capacity), a fertility factor (an estimate of site productivity, included as the maximum productivity in megagrams per hectare per year), average monthly climatic conditions (temperature and precipitation and their standard deviations), seeds per square metre (to indicate the potential number of seedlings of each species that may germinate), and seedling stocking (to indicate the proportion of plot area covered by seedlings of a given species). ZELIG-CFS has been evaluated in North American mixed forests and each species is characterized by a set of species-specific allometric and ecological parameters (e.g., maximum age, shade, and drought tolerance class; [68]). For this study, the set of parameters from Laroque et al. (2011) was maintained, except for the parameter controlling growing season temperature (minimum degree days) for trembling aspen (Populus tremuloides Michx.), which was lowered (from 889 to 581) to be consistent with the North American version of the FORCLIM model [73], and for white pine, which was lowered (from 1500 to 1204) to equal that of sugar maple thus allowing its growth in northwestern Ontario at the current northern edge of its range in Ontario.
Site-specific input data such as plot location, elevation, and soil texture type were collected during plot establishment before treatments were applied. The fertility factor was calculated using information from the natural treatment plots and the ZELIG-CFS fertility factor calculator. These values were referenced against productivity values from Ontario's growth and yield permanent plot network and lowered slightly for Sioux Lookout to better reflect observed growth rates of larger jack pine.
To investigate differences among silvicultural treatments in all locations, simulations were run under a scenario of current climate based on the 1980-2010 period. Climate series were derived using BioSIM version 11 (Laurentian Forestry Centre, Quebec City, QC, Canada), which uses historical weather data from the nearest weather stations to predict past, current, or future weather based on adjustments for differences in elevation, latitude, and longitude [74]. To assess the effect of changing climates, two climate change scenarios were simulated in Sioux Lookout, Dryden, and Petawawa (i.e., one of each forest type); treatments simulated were natural, extensive, and elite. Climate change scenarios were based on projected climates using projected normals from 2020 to 2050 (i.e., roughly the midpoint of stand rotation age) based on two standard representative concentration pathways (RCP), 2.6 and 4.5, from the Canadian Earth System Model version 2 (CanESM2; [75]). While temperatures were generally higher in the climate change scenarios, only a marginal difference was found between RCP scenarios 2.5 and 4.6; precipitation was generally consistent among all three climate scenarios due to a large amount of variability in precipitation predictions (Tables S1-S3).
Seeds per square meter and seedling stocking were estimated for each plot using the stocking quadrat information. Seeds per square meter were initially estimated as the number of seedlings (individuals per tree species with sizes below 1.3 m in height, i.e., without a DBH) present in the quadrats divided by 160 m 2 (i.e., the total area of the quadrats in each experimental unit) and stocking was estimated as the proportion of quadrats in each experimental unit including at least one individual of each species below the 2.5 cm DBH threshold. However, ZELIG-CFS assumes the number of potential seedlings is constant throughout stand development. This number of potential regenerating seedlings is then modified by understory light conditions, soil moisture and fertility, degree days, and random variation. Since our stands were very young at initialization, the regeneration plots had substantial numbers of seedlings, leading to the model estimating a high number of potential regenerating seedlings and subsequently estimating extremely high recruitment levels. To correct for this overestimate, both seeds per square meter and stocking were adjusted based on species autecology as follows: seeds per square meter were adjusted to the number of quadrats with a seedling present divided by the total area measured. This approach removed quadrats where some species, particularly balsam fir, black spruce, and sugar maple, had extremely concentrated seedling densities that did not reflect average conditions or expected regeneration across the site. Furthermore, for all species below a shade tolerance value of 5 (i.e., the maximum shade tolerance in ZELIG-CFS, or extremely shade-tolerant) all stems less than 50 cm in height were removed from the seeds per square meter analysis. This adjustment better reflects the reduced ability of these species to regenerate under closed canopy conditions while maintaining some regeneration for early gaps and maintenance of secondary cohorts in suppressed canopy positions. Personal communication with area foresters noted quick die-off of birch (Betula spp.) seedlings in particular. Stocking was estimated as the proportion of quadrats with a recorded percent cover of at least one individual of each species smaller than the 2.5 cm DBH threshold. Both seeds per square meter and stocking values were multiplied by the relative shade tolerance of the species (0.2 reflecting least shade-tolerant and 1 reflecting most shade-tolerant, or species shade tolerance divided by the maximum of 5) to better reflect successional dynamics.
The tree list used to initialize ZELIG-CFS was gathered in each experimental unit during the tenth-year post-treatment survey. When recruits were only measured in one quadrant, the tree list of recruits was replicated three times to produce an estimated tree list for the full 400 m 2 plot. The tree lists generated from the four growth and yield plots per experimental unit were combined. Each experimental unit on each site was simulated for 70 years. Simulations ran on annual time steps and outputs were collected for each year.

Silvicultural Interventions
Silvicultural interventions were implemented in our simulations as described in Bell et al. (2017) and were consistent with recommendations in the provincial silviculture guide [76]. A summary of all interventions is provided in Table 1. At Sioux Lookout, intensive and elite treatments were thinned from below to about 2500 stems ha −1 in the 20th year post-harvest and thinned from above to remove 25% of jack pine basal area and all pin cherry (Prunus pensylvanica L.f.) at 35 years post-harvest. At Kapuskasing, 30% of the basal area of both white spruce and black spruce were thinned from below in the 45th year in both intensive and elite treatments. At Dryden, 30% of white spruce were thinned from below and 75% of white birch (Betula papyrifera Marsh.) were thinned from above in the 40th year in both intensive and elite treatments. At Timmins, 75% of trembling aspen and 50% of white birch were thinned from above in the 40th year post-harvest in both intensive and elite treatments. At Petawawa, between 40% and 60% of white pine (depending on the size of individuals in each experimental unit; all individuals greater than 12 m in height were removed) and 50% of red maple (Acer rubrum L.) were thinned from above to simulate canopy removal in the basic, intensive, and elite treatments at the 40th year post-harvest. Finally, at North Bay, 50% of sugar maple and other hardwoods except yellow birch were thinned from above in the 40th year in intensive and elite treatments to simulate overstory removal and 50% of sugar maple was thinned from below at the 60th year post-harvest in both treatments. Table 1. Simulated silvicultural interventions for sites in a silvicultural intensity study in Ontario, Canada. Parameters included in the intervention files for the PartialCut module in ZELIG-CFS. All treatments differed in initial planting density, planting stock quality, site preparation, and initial silvicultural treatments. Extensive treatments had no simulated silvicultural interventions.

Site Basic Treatment Intensive and Elite Treatments
Sioux Lookout None Jack pine thinned from below to 2500 stems per hectare at 20th year; 25% thinned from above at 35th year post-harvest Kapuskasing None 30% basal area of white and black spruce thinned from below at 45th year post-harvest Dryden None 30% of white spruce thinned from below and 75% of white birch thinned from above in the 40th year post-harvest Timmins None 75% of trembling aspen and 50% of white birch were thinned from above in the 40th year post-harvest Petawawa All white pine taller than 12 m removed and 50% of red maple thinned from above to simulate canopy removal in 40th year post-harvest All white pine taller than 12 m removed and 50% of red maple thinned from above to simulate canopy removal in 40th year post-harvest North Bay None 50% of the sugar maple and other hardwoods, except yellow birch, thinned from above at the 40th year post-harvest to simulate overstory removal and 50% of sugar maple thinned from below at the 60th year post-harvest

Analysis
The model output includes a list of individual trees at each time step. Volume yield was estimated from the model output using Honer's volume equations for gross total volume by species [77]. Standing volume was the total volume of live individuals in each experimental unit, harvested volume was the total removed volume over the course of the simulation, and the total volume was the summation of final standing volume and harvested volume. The relative abundance was calculated from the model output as the relative proportion of each species basal area to the total basal area in the experimental unit for each time step. For visualization purposes, we chose standing volume and species composition because these two variables are often the most pertinent to forest managers in Ontario. However, other important characteristics for informing management decisions, such as diameter distributions or height of stems, can also be derived from the outputs of ZELIG-CFS. Changes in standing volume and species composition were examined graphically. Descriptive analyses were chosen because estimating the statistical significance of annual changes in volume or species composition was not the goal of the study.

Results
At the end of the simulation period, the basic treatments resulted in the highest standing gross volume among silvicultural treatments at both boreal conifer sites, Sioux Lookout and Kapuskasing ( Figure 1). However, when considering the amount of volume removed, intensive and elite treatments supported comparable final volumes at Sioux Lookout and higher final gross volumes at Kapuskasing (Table 2). Extensive treatments resulted in slightly lower average volume, but, at Kapuskasing, simulation outcomes varied more than those of other treatments (Figure 1). On natural (unharvested and untreated) plots, volume increased initially at both sites, but quickly reached a maximum and generally remained constant thereafter ( Figure 1). The four post-harvesting treatments reached volumes comparable to those in the natural treatment within the 70-year simulation period. Boreal mixedwoods demonstrated opposite outcomes. The extensive treatments resulted in, on average, higher volumes across the entire simulation period at both Dryden and Timmins (Figure 2). Including removed volume brought the intensive and elite treatments in line with the extensive treatment at Dryden, but these treatments produced considerably more volume than the extensive treatments at Timmins (Table 2). Natural treatments showed disparate responses: Over the 70-year simulation period, volumes increased moderately at Dryden, while they decreased by almost half at Timmins ( Figure  2). By the end of the simulation period, all harvested treatments reached the initial starting volumes of the natural stands at Dryden, but only intensive and elite treatments had volumes similar to initial volumes in the natural treatments at Timmins (Figure 2).  Boreal mixedwoods demonstrated opposite outcomes. The extensive treatments resulted in, on average, higher volumes across the entire simulation period at both Dryden and Timmins (Figure 2). Including removed volume brought the intensive and elite treatments in line with the extensive treatment at Dryden, but these treatments produced considerably more volume than the extensive treatments at Timmins (Table 2). Natural treatments showed disparate responses: Over the 70-year simulation period, volumes increased moderately at Dryden, while they decreased by almost half at Timmins (Figure 2). By the end of the simulation period, all harvested treatments reached the initial starting volumes of the natural stands at Dryden, but only intensive and elite treatments had volumes similar to initial volumes in the natural treatments at Timmins (Figure 2). At the temperate sites, North Bay and Petawawa, the effect of silvicultural intensity was noticeable. At both sites, by the end of the simulation period extensive treatments resulted in the highest average standing volume among post-harvesting treatments but showed very little gain; in fact, at North Bay, the extensive treatment resulted in a small loss in standing volume (Figure 3). Basic treatments at both sites resulted in a small increase in volume, more so at North Bay where stands were not thinned (Figure 3). Although final standing volumes were relatively low, the highest total volumes occurred in the intensive and elite treatments when the volume removed during overstory harvesting 40 years after the first harvest was factored in (Table 2, Figure 3). Over the 70-year simulation period, volume in natural treatments increased at Petawawa whereas it decreased somewhat at North Bay; at the end of the simulation period, standing volume in the treated stands differed from that in natural stands (Figure 3). At the temperate sites, North Bay and Petawawa, the effect of silvicultural intensity was noticeable. At both sites, by the end of the simulation period extensive treatments resulted in the highest average standing volume among post-harvesting treatments but showed very little gain; in fact, at North Bay, the extensive treatment resulted in a small loss in standing volume (Figure 3). Basic treatments at both sites resulted in a small increase in volume, more so at North Bay where stands were not thinned (Figure 3). Although final standing volumes were relatively low, the highest total volumes occurred in the intensive and elite treatments when the volume removed during overstory harvesting 40 years after the first harvest was factored in (Table 2, Figure 3). Over the 70-year simulation period, volume in natural treatments increased at Petawawa whereas it decreased somewhat at North Bay; at the end of the simulation period, standing volume in the treated stands differed from that in natural stands (Figure 3). At the boreal conifer sites, changes in species composition through time and along silvicultural intensities occurred as expected. Interestingly, at Sioux Lookout the extensive treatment was sufficient to regenerate predominately jack pine; however, at Kapuskasing, extensive treatments regenerated mostly to trembling aspen (Figure 4). In the boreal mixedwoods, only the intensive and elite treatments regenerated mostly to spruce; extensive and basic treatments regenerated to poplar and birch, which were the dominant tree species at the 10th year plot surveys ( Figure 5). The species composition of the natural stands in the boreal mixedwoods remained mostly stable, with a slight shift towards a spruce-fir stand type at Timmins near the end of the simulation period ( Figure 5). Finally, composition at both temperate sites remained relatively stable through time and across silviculture intensities, except at North Bay, where yellow birch dominated after the maple was removed via thinning ( Figure 6). The relative species composition at these sites was comparable to that in the natural stands ( Figure 6). At the boreal conifer sites, changes in species composition through time and along silvicultural intensities occurred as expected. Interestingly, at Sioux Lookout the extensive treatment was sufficient to regenerate predominately jack pine; however, at Kapuskasing, extensive treatments regenerated mostly to trembling aspen (Figure 4). In the boreal mixedwoods, only the intensive and elite treatments regenerated mostly to spruce; extensive and basic treatments regenerated to poplar and birch, which were the dominant tree species at the 10th year plot surveys ( Figure 5). The species composition of the natural stands in the boreal mixedwoods remained mostly stable, with a slight shift towards a spruce-fir stand type at Timmins near the end of the simulation period ( Figure 5). Finally, composition at both temperate sites remained relatively stable through time and across silviculture intensities, except at North Bay, where yellow birch dominated after the maple was removed via thinning ( Figure 6). The relative species composition at these sites was comparable to that in the natural stands ( Figure 6).    Stand dynamics and predicted volumes were sensitive to climate change ( Figure 7). Notably, at Sioux Lookout, a substantial shift was evident away from mature, closed, and high yield jack pine stands towards a lower-volume open canopy state. Harvested stands were predicted to fail to reestablish and natural stands were projected to decline over time, which contrasts sharply with current climate scenarios that project they will remain relatively stable (Figure 1, Figure 7). At Dryden and Petawawa, on the other hand, a general increase in volume was projected for the two climate change scenarios (Figure 7). Stand dynamics and predicted volumes were sensitive to climate change (Figure 7). Notably, at Sioux Lookout, a substantial shift was evident away from mature, closed, and high yield jack pine stands towards a lower-volume open canopy state. Harvested stands were predicted to fail to reestablish and natural stands were projected to decline over time, which contrasts sharply with current climate scenarios that project they will remain relatively stable (Figures 1 and 7). At Dryden and Petawawa, on the other hand, a general increase in volume was projected for the two climate change scenarios (Figure 7).

Figure 7.
Simulated gross standing volume from the 10th to the 80th year post-harvest for the three selected sites in the silvicultural intensity study in Ontario, Canada, including climate forcing. Lines represent a LOESS smooth with a span of four years of gross standing volume of each treatment at each site. Ribbons represent the minimum and maximum gross standing volume at each year, for each treatment, at each site. Climate forcings followed Relative Concentration Pathways 2.6 (low) and 4.5 (moderate) climate change scenarios. Line type represents the two climate scenarios. Basic and Intensive treatments are omitted for readability.

Discussion
We simulated forest dynamics for 70 years using a gap model and a plot network specifically designed to address the paradigm shifts facing forest management today. We simulated a range of silvicultural systems and interventions: from unharvested stands to clearcuts and shelterwood cuts with no interventions to clearcuts and shelterwood cuts with substantial silvicultural intervention. Stand complexity ranged widely as well: from even-aged stands with one hyper-dominant species (i.e., jack pine clearcuts at Sioux Lookout) to multiple age cohort complex stands with more than 10 co-occurring species (i.e., temperate shelterwood systems at Petawawa). Stand-level volumes generally increased as silvicultural intensities increased and simulated species compositional changes aligned with expectations. Finally, we were able to assess the effects of climate change on our sites. During our simulations, the model itself and the initialization file and intervention file structures were consistent across all sites. Our results suggest mechanistic models can be used to project a range of forest conditions and silvicultural interventions. Thus, mechanistic models are an underused asset for forest management and have the potential to both streamline and enhance forest management planning [17]. Simulated gross standing volume from the 10th to the 80th year post-harvest for the three selected sites in the silvicultural intensity study in Ontario, Canada, including climate forcing. Lines represent a LOESS smooth with a span of four years of gross standing volume of each treatment at each site. Ribbons represent the minimum and maximum gross standing volume at each year, for each treatment, at each site. Climate forcings followed Relative Concentration Pathways 2.6 (low) and 4.5 (moderate) climate change scenarios. Line type represents the two climate scenarios. Basic and Intensive treatments are omitted for readability.

Discussion
We simulated forest dynamics for 70 years using a gap model and a plot network specifically designed to address the paradigm shifts facing forest management today. We simulated a range of silvicultural systems and interventions: from unharvested stands to clearcuts and shelterwood cuts with no interventions to clearcuts and shelterwood cuts with substantial silvicultural intervention. Stand complexity ranged widely as well: from even-aged stands with one hyper-dominant species (i.e., jack pine clearcuts at Sioux Lookout) to multiple age cohort complex stands with more than 10 co-occurring species (i.e., temperate shelterwood systems at Petawawa). Stand-level volumes generally increased as silvicultural intensities increased and simulated species compositional changes aligned with expectations. Finally, we were able to assess the effects of climate change on our sites. During our simulations, the model itself and the initialization file and intervention file structures were consistent across all sites. Our results suggest mechanistic models can be used to project a range of forest conditions and silvicultural interventions. Thus, mechanistic models are an underused asset for forest management and have the potential to both streamline and enhance forest management planning [17].

Simulation of Silvicultural Intensities
Simulation of silvicultural intensities generally produced outcomes that aligned with the intent of each treatment. While intensive and elite treatments often ended with similar, or slightly lower, standing volumes than those of the extensive and basic treatments, the total volume (i.e., final standing volume plus removed volume) was often much higher for the former. This finding is comparable to results of previous studies on the effects of commercial thinning in other Canadian forests, including for jack pine [78], spruce [79], and temperate pine stands [80]. Our results combined with those from previous investigations suggest that commercial thinning may reflect an opportunity to improve overall yields while producing high-quality lumber at rotation age and preserving desired species composition. Regarding harvest modelling, ZELIG-CFS only supports thinning based on the percentage of standing basal area at a given time step. To better incorporate forest management planning and new forest management tools such as density management diagrams [81], future versions of the model should consider supporting thinning down to target stem counts per hectare or target specific stand basal areas rather than percentages.

Simulation of Mixed-Species and Multiple Cohort Stands
While our simulations generally returned expected outcomes, model outputs aligned with expectations more when the initial overstory composition was similar to that of the desired end goal (e.g., Dryden, Sioux Lookout, Kapuskasing). Sites with significant changes in species composition (e.g., Timmins) or an overstory removal (e.g., Petawawa, North Bay) had somewhat lower than expected final volumes. While the secondary cohort (i.e., individuals in a suppressed canopy position) supported the expected species where the canopy was removed, the growth response of this cohort to release was likely underestimated. For Timmins, for instance, final total stand volumes (i.e., including removals) averaged between 250 m 3 ha −1 and 350 m 3 ha −1 , somewhat lower than the 400 m 3 ha −1 predicted from similar aspen-spruce plantations in Alberta [82]. This result suggests that the regeneration submodel in ZELIG-CFS should be refined to improve the representation of regeneration processes, particularly for planted seedlings. While regeneration is modelled, the number of potential seedlings and stocking values are held constant through time. The number of seedlings that regenerate are then calculated annually from the number of potential seedlings and stocking values modified by understory light conditions, soil moisture and fertility, degree-days, and random variation [68]. As a result, after a thinning release, recruitment pulses occur multiple times rather than just after the thinning, which can reduce available growing space and cause more competition to be modelled than would actually occur. Regardless, trends in yields aligned with expectations; generally, increasing silvicultural intervention led to similar or higher total volumes relative to no intervention.

Simulation of Stand Dynamics Related to Changes in Species Composition
Changes in species composition generally followed the intended outcomes of the silvicultural treatments. At all boreal sites, harvesting without subsequent vegetation management increased the hardwood component as expected [52,53]. The only site where a significant number of conifers were established without vegetation management was at Sioux Lookout. Here, the previous forest type was majority jack pine with no detectable hardwood component other than a minimal amount of pin cherry. The silviculture clearcut system and sandy soils clearly favoured jack pine establishment [63]. Increases in the hardwood proportion were far higher at sites where poplar and other hardwoods were present prior to harvesting. At Kapuskasing and both boreal mixedwoods, the hardwood proportion increased in all extensive treatments relative to that in the natural stands. However, as silvicultural intensities increased the proportion of hardwoods decreased, indicating that the simulation model was able to capture the long-term effects of harvesting on species composition both with and without silvicultural interventions.
One aspect where the model performed unexpectedly was the maintenance of a fir component. This issue is perhaps most notable at Kapuskasing, where fir was 25% to 50% of the initial species composition in basic, intensive, and elite treatments but decreased to less than 15% on average. Similar, although less drastic, decreases in fir occurred at all sites. This result was unexpected, particularly for older stands, as fir is a gapspecialist that is able to quickly colonize patches and can shade out other regenerating species. In managed boreal forests with longer fire return intervals, it has been found to increase in proportion [51]. It is possible that the model overestimated the mortality of smaller individuals of shade-tolerant species that are able to survive with low growth rates in suppressed positions over long periods. Despite this, the model performed well in simulating species composition of the aspen-spruce intimate mixture at Timmins. While total volumes were likely underestimated, removing the aspen canopy immediately shifted the stand type to mostly spruce, providing further evidence that the overtopping of spruce by poplar (hybrid or otherwise) for the first few decades is a viable strategy that can lead to a significant volume of spruce [41,82].
At the temperate sites, the same overstory species maintained dominance throughout the rotation, which is expected from shelterwood systems. At North Bay, increasing silviculture intensity was necessary to regenerate a substantial quantity of yellow birch, likely because the yellow birch was shaded and unable to attain merchantable sizes under the dense sugar maple canopy, further reinforcing that silvicultural interventions are required to maintain this important, but often minor, species [83]. Overstory removals did not affect the dominance of white pine at Petawawa since most regeneration was also white pine. However, intensive and elite treatments supported higher proportions of white pine at the end of the simulation period, possibly because of the more intense vegetation control, which aids white pine regeneration [84] leading to a lower proportion of red maple in the stands.

Simulation of Climate Change Effects
Climate change is expected to alter temperature, precipitation, and seasonality with potentially acute effects on Canada's forests [54,85]. When climate variables are not included in growth models, the implicit assumption is that the rate of tree growth observed when the data used to calibrate the model was collected will remain constant into the future [29]. To assess potential climate change effects, we simulated stand dynamics under changing climate scenarios for three forest types, observing responses ranging from near-collapse of boreal jack pine to enhanced growth of temperate hardwoods. Simulated response of forest stands to climate change was quite harsh in Sioux Lookout; while the combination of increased temperatures and reduced precipitation could shift the closed canopy boreal conifer biome to a more woodland biome type [86], using this model with even mild climate change scenarios led to complete collapse. This outcome is likely because of the quadratic relationship between temperature and growth in ZELIG-CFS. Jack pine was growing near its growing degree limit specified in the model (particularly in the RCP 4.5 scenario) and had lowered projected growth due to temperature that differed from the current climate. This lowered growth increased the mortality rates in the model and led to jack pine failing to establish. While a threshold of decreased growth in boreal tree species due to higher temperatures has been shown for a 3 • C temperature increase in experimental warming treatments [87], space-for-time modelling suggests a moderate increase in growth for temperature increases below that level [88]. Thus, the reduction in growth and subsequent mortality was likely higher than what would be likely to occur based on the increase in temperature alone. Other forest growth models assume asymptotic growth with temperature [89] which might lead to a more robust simulation of the effects of climate change when boreal species are not competing with other temperate species.
At Dryden and Petawawa, increased growth was evident in climate change simulations compared to that under the current climate, likely because while projections for total precipitation values were generally constant, temperature increased steadily. These increases in temperature benefited the more temperate species mix (particularly the minor white pine component at Dryden) and the poplar species. Both of these groups had temperature limited growth in the current climate scenario at Dryden and had very high degree day maximums in the model (i.e., no reduction with increasing temperatures as at Sioux Lookout). The outcome of the climate change simulations for Dryden suggests that attempts to increase temperate species mix in advance of climate change may be a way to prevent ecosystem collapse such as that projected for the Sioux Lookout site. It is important to note that the current version of ZELIG-CFS does not include effects of potential CO 2 fertilization. Although other process-based models have incorporated CO 2 fertilization as a potential variable of interest [90], the lack of CO 2 fertilization in our models is likely not a concern because the effects of increasing CO 2 on tree growth in Canada are equivocal with reports of both increases [91] and no changes in growth with increasing atmospheric CO 2 [92].

Management Implications
The use of mechanistic models in the framework of the forest management process has been suggested for some time [11,[17][18][19]; however, their use in planning has been limited, particularly in Canada [19]. Here we demonstrate that a currently available mechanistic model can be used to provide estimates of rotation age yield and species composition across a range of forest conditions. This approach has the potential to streamline the forest management process and allow for the exploration of more management and climate scenarios. Furthermore, the availability of new remote sensing data sets such as broadscale high-resolution aerial light detection and ranging (LIDAR), coupled with advances in computer processing power and speed, will allow for more of the forest to be explicitly and accurately modelled than in the past. Building on previous modelling studies using both photo interpretation to examine stand development [9,93] and aerial LIDAR to examine yields [94] will improve the applicability and adaptability of these mechanistic models. Advanced forest models and computer processing capabilities will enable forest managers to capitalize on the full potential of enhanced forest resource inventories and airborne LIDAR data.
However, as with any other mathematical model that represents a simplification of complex ecological processes, models like ZELIG-CFS have limitations. In our study, the limitations mostly related to a disconnect in stocking estimates and regeneration modelling. For mechanistic models such as ZELIG-CFS to be fully incorporated into the planning process, they need to be modified to allow for changes in regeneration through time and to account for initial planting densities. While we focused on the implementation of silviculture treatments and climate change scenarios, the incorporation of natural disasters such as spruce budworm outbreaks and wildfire would greatly increase the accuracy of ZELIG-CFS outcomes [26]. Disturbances are important considerations, particularly at the landscape scale and under expected climate change, which may result in both boreal and temperate forests being more susceptible to natural disturbances [95,96]. Model applicability to forest managers is highly dependent on ease of use. ZELIG-CFS is managed by a user friendly graphical tool, AMSIMOD (Application for the Management of Simulation MODels). Inputs to ZELIG-CFS include diameter at breast height and species codes from sample plots and climate information. AMSIMOD can generate outputs of ZELIG-CFS in the form of size distributions, line and bar charts, and for export to Geographic Information System mapping software. Future quality of life improvements will include integration of climate generation software into the AMSIMOD environment. Still, as currently constructed, this mechanistic model can provide a quicker estimate of expected stand outcomes than those empirical models that require stands to develop to a point where a site index can be estimated for their use or those that use a prior stand index which further compounds issues related to using past conditions to predict future outcomes in an uncertain and complex future [19].

Conclusions
This study demonstrated the flexibility of the gap model ZELIG-CFS when applied to forest management scenarios. Mechanistic models can enable forest managers to examine the effects of intensified silviculture, understand trade-offs between increasing complexity and expected yields, estimate end of rotation species compositions, and incorporate the effects of climate change. This flexibility could improve the planning process and enable forest managers to examine a much broader range of prescriptions than is possible with traditional yield tables and allow them to incorporate the paradigm shifts in forest management into their modelling frameworks. Future improvements are suggested for ZELIG-CFS in particular, especially if it is applied to forest management: namely, modifications in the regeneration submodel to better model changes in regeneration through time and incorporate planted seedlings, and an option for highly shade-tolerant individuals to survive long periods with low growth rates. However, these modifications will require better understanding of the dynamics of seedling establishment and survival in the early years by establishing appropriate field experiments. Mechanistic models should also be constructed to support the use of new data sources from advanced forest inventories.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.