Modelling Post-Disturbance Successional Dynamics of the Canadian Boreal Mixedwoods

Abstract

Natural disturbances, such as fire and insect outbreaks, play important roles in natural forest dynamics, which are characterized over long time scales by changes in stand composition and structure. Individual-based forest simulators could help explain and predict the response of forest ecosystems to different disturbances, silvicultural treatments, or environmental stressors. This study evaluated the ability of the SORTIE-ND simulator to reproduce post-disturbance dynamics of the boreal mixedwoods of eastern Canada. In 1991 and 2009, we sampled all trees (including seedlings and saplings) in 431 (256 m²) plots located in the Lake Duparquet Research and Teaching Forest (western Quebec). These plots were distributed in stands originating from seven wildfires that occurred between 1760 and 1944, and which represented a chronosequence of post-disturbance stand development. We used the 1991 inventory data to parameterize the model, and simulated short- to long-term natural dynamics of post-fire stands in both the absence and presence of a spruce budworm outbreak. We compared short-term simulated stand composition and structure with those observed in 2009 using a chronosequence approach. The model successfully generated the composition and structure of empirical observations. In long-term simulations, species dominance of old-growth forests was not accurately estimated, due to possible differences in stand compositions following wildfires and to differences in stand disturbance histories. Mid- to long-term simulations showed that the secondary disturbance incurred by spruce budworm did not cause substantial changes in early successional stages while setting back the successional dynamics of middle-aged stands and accelerating the dominance of white cedar in late-successional post-fire stands. We conclude that constructing a model with appropriate information regarding stand composition and disturbance history considerably increases the strength and accuracy of the model to reproduce the natural dynamics of post-disturbance boreal mixedwoods.

1. Introduction

Simulation of post-fire succession in forest ecosystems is a tool for processing and optimizing knowledge contained within empirical data that are collected through field measurements, and for providing a better vision of future dynamics of forest systems. In North American boreal forests, boreal mixedwoods are an important source of timber [1]. Among boreal forest ecosystems, mixedwoods are recognized as being structurally heterogeneous and the most productive [2]. In addition to traditionally recognized effects of site characteristics and climate on forest stand distributions at different scales through the regulation of the recruitment and growth of tree species [3], disturbances, which are partly driven by climate [4], also modify the composition and structure of forest ecosystems [5]. Specifically, the structure and composition of boreal forest ecosystems are widely influenced by insect outbreaks [6] and wildfires [7,8], which spatially and temporally vary in size, frequency, and severity.

Mixedwoods of the southern Canadian boreal forest consist of stands that are dominated by shade-tolerant conifers mixed with intolerant hardwoods, particularly trembling aspen (Populus tremuloides Michx.) [9]. Generally, after a stand-replacing disturbance, such as fire, the structural pattern and species composition of boreal mixedwoods are the consequence of stand dynamics [2,8], during which a gradual transition occurs, from early successional hardwoods to admixtures of hardwood–conifer species to conifer-dominated stands [10]. These three successional developmental stages of the boreal mixedwoods are at the core of the three-cohort model that was developed by Bergeron et al. [10] to depict forest successional dynamics. Through processes, such as mortality of early successional species and recruitment of mid- and late-successional species, each of these development stages co-exists, which results in more complex multi-cohort forest structures at the landscape level [11,12]. If intervals between fires (fire cycles) are short, then changes in species dominance are limited, which leads to pre-fire and post-fire stands with similar compositions. More specifically, this phenomenon is observed in stands that are dominated by trees bearing serotinous cones, such as jack pine (Pinus banksiana Lamb.) and black spruce (Picea mariana [Mill.] B.S.P.), and stands that are dominated by species possessing root systems that can survive fire, such as aspen and white or paper birch (Betula papyrifera Marsh.) [8,13,14]. During the development of mixedwoods, gaps that are generated by partial disturbances, such as insect outbreaks and windthrow, together with tree- and stand-level processes, such as competition and senescence, facilitate the establishment and growth of shade-tolerant conifers [8]. These species might stay suppressed in the understory for long periods before the formation of canopy gaps (and the concomitant increase in light availability) induces an acceleration of their growth rate. Thus, gap dynamics act as modulators of stand composition and structure in late-successional conifer stands [15,16]. Moreover, changes in canopy composition are the consequences of differential growth rates among species [17] and not just a replacement of species that is driven by numerous inhibition or facilitation processes. That explains the occurrence of some late-successional species immediately following fire, although species, such as balsam fir (Abies balsamea (L.) Miller) and eastern white cedar (Thuja occidentalis L.), are favored by the increased availability of suitable seedbeds in older mixed or conifer stands [8].

Stand replacement dynamics following disturbances depend upon the nature of the disturbance, the composition of the disturbed stands, and site conditions [8,18,19]. For instance, seedbed limitations [20,21], lack of seed sources of shade-tolerant conifers [8,22], and their poor establishment and survival beneath hardwood canopies [23,24] prolong the dominance of shade-intolerant or pioneer species. In addition, insect outbreaks can increase the regeneration of hardwood species if conifer regeneration is not sufficient [12]. Defoliation affects tree regeneration and the soil organic layer less severely than wildfire does, and exerts a species-specific influence on forest stands [2,25]. Insect defoliation often extends over several years and has varying effects on species and stand development dynamics. For example, repeated and severe defoliation of aspen stands by forest tent caterpillar (Malacosoma disstria Hübner; FTC) accelerates the mortality of mature aspen trees and facilitates the release and rapid growth of shade-tolerant conifers, such as spruce (Picea spp.) and balsam fir, into the canopy [26,27]. In fir-dominated stands, outbreaks of spruce budworm (Choristoneura fumiferana [Clemens]), a defoliator of both balsam fir and spruce, may help to maintain the continuous dominance of balsam fir, through the growth release of suppressed saplings that results from the formation of canopy openings [28,29]. Budworm outbreak may also induce the conversion of conifer-dominated stands to mixedwood stands by favoring the recruitment of shade-intolerant tree species, such as trembling aspen and white birch [30,31]. Moreover, the effects of other insects that selectively kill trees, such as bark beetles (e.g., Dendroctonus ponderosae Hopkins), may considerably affect soil organic content due to additional litter inputs from killed trees [32]. The successional dynamics of boreal mixedwoods from simple structured stands to more complex multi-cohort forest structures are thus the result of complex interactions that are controlled by numerous allogenic and autogenic factors [12]. Several factors, such as stand attributes (e.g., species composition and age), prior to and after disturbances, as well as disturbance severity, intensity and spatial distribution, and time since disturbance define the level of complexity in successional stand dynamics [33].

In the current study, we took advantage of empirical data from 431 forest inventory plots that were sampled in 1991 and resampled in 2009 in the Abitibi region of Quebec, Canada, to simulate the natural successional dynamics of stands originating from seven wildfires in the eastern Canadian boreal mixedwoods. To do so, we used SORTIE-ND, a spatially explicit individual-based forest stand dynamics model [34]. This mechanistic model has been previously used to explore natural forest dynamics in a number of forest systems in Canada [35,36,37,38] and elsewhere in the world [39,40,41,42]. Spatially explicit individual-based models have proved to be satisfactory tools in predicting stand succession based upon species autoecology [43]; we hypothesized that SORTIE-ND would be able to reproduce the compositional and structural dynamics of boreal post-fire stands from early to late-successional stages. Our first objective was to parametrize SORTIE-ND for simulating the dynamics of post-fire mixedwoods, which were sampled in 1991 over 18 years. Benefiting from the empirical data that were sampled in 2009, we evaluated model performance in terms of the species composition, stem density, and size distribution of trees within stands. Our second objective was to evaluate the contribution of the recent spruce budworm outbreak to current post-fire stand composition and structure in 2019. The recent outbreak extended from 1972 to 1987 and was severe, defoliating ca. 56% (of the total number of stems per hectare) of host tree species, i.e., fir and spruce [44]. In taking advantage of empirical data in 1991 that account for the status of stands both with and without including the effect of the recent spruce budworm outbreak, we carried out two 28-year simulations that culminated in the current age of stands in 2019. Comparing the two 28-year simulation outputs for data with and without the outbreak effect enabled us to evaluate the short-term influence of spruce budworm outbreak on the post-fire composition and structure of boreal mixedwoods stands that were located at the Lake Duparquet Research and Teaching Forest (LDRTF). Finally, conducting long-term simulations, our third objective was to determine how the reoccurrence of spruce budworm outbreaks would contribute to changes in the natural successional dynamics of post-fire mixedwoods over a century.

2. Materials and Methods

2.1. Study Area

Our study landscape was located in the boreal mixedwoods of eastern Canada, more specifically at the LDRTF (79°19’ W–79°30’ W, 48°86’ N–48°32’ N), where a dozen large wildfires occurred between 1760 and 1944. The study landscape is located within the balsam fir-white birch bioclimatic domain [45]. For the 1981–2010 period, the mean annual temperature was 1.0 °C, mean annual total precipitation was 985 mm, of which 30% falls during the growing season, and the average number of degree days (5 °C base temperature) was 1350 (Mont Brun weather station) [46]. Forests within the study area are characterized by a mixed composition of boreal conifers and shade-intolerant hardwood species. More specifically, early successional stands are dominated by trembling aspen, white birch, and jack pine, whereas balsam fir and eastern white cedar, in association with white birch, dominate late-successional stands [8].

The combined influences of catastrophic wildfires and defoliating insect outbreaks are the main natural disturbances of the region. The fire regime of the past 300 years has been reconstructed and is characterized by high-intensity fires covering vast areas, particularly in terrain with flat topography [47]. During that period, fire return intervals increased from 83 years prior to 1850 to 325 years since 1920 [48]. During the 20th century, three outbreaks of eastern spruce budworm were reported in the territory, for the periods 1919–1929, 1930–1950, and 1972–1987 [49,50]. Bergeron et al. [44] qualified the last spruce outbreak as being very severe, i.e., killing more than 56% (percentage of the total number of stems per hectare) of host species (i.e., balsam fir and spruces), especially balsam fir. The forest tent caterpillar, with lesser effects on host species mortality [51] and shorter outbreak cycles than spruce budworm [52], defoliates hardwood species, particularly trembling aspen. In Quebec, areas that were affected by FTC outbreaks have been reported since 1985 and in 2018, 38,903 ha of forest were affected across the province [53].

2.2. Sampling Design

We used a network of 431 sampling plots, which were distributed systematically within the study area. Originally, the plots were set up to observe changes in stand composition and structure in stands originating from seven fires representing a chronosequence of 249 years of forest succession [54,55]. The study plots (16 m × 16 m, or 256 m²) were located every 50 m along transects that were established in each burned area and which were sampled in 1991 and 2009 (i.e., 18 years apart).

Table 1. Characteristics of the seven study fires in the measurement year, 1991.

Fire Date	Fire Area * (ha)	Stand Age (year)	Number of Transects	Plots Count	Mean Basal Area (m2ha−1, ±SE) of Standing Live Trees (DBH > 1 cm)
					Fir	Birch	Spruce	Aspen	Cedar
1760	>7760	231	6	54	2.60 (0.55)	6.69 (0.86)	2.04 (0.40)	2.44 (0.84)	10.89 (1.46)
1797	178	194	4	50	2.77 (0.42)	7.92 (0.83)	2.08 (0.39)	0.50 (0.20)	3.84 (0.62)
1823	288	168	5	66	4.32 (0.47)	3.88 (0.37)	2.46 (0.38)	8.03 (1.03)	0.53 (0.26)
1847	122	144	4	74	6.84 (0.50)	4.62 (0.45)	5.40 (0.58)	11.88 (1.00)	0.77 (0.29)
1870	555	121	6	64	3.46 (0.36)	5.86 0.60)	4.46 (0.56)	12.78 (1.51)	0.02 (0.02)
1916	35	75	3	52	4.54 (0.40)	7.63 (0.83)	1.99 (0.57)	16.74 (2.01)	0.97 (0.36)
1944	>298	47	4	71	3.81 (0.48)	9.34 (0.76)	0.81 (0.19)	11.57 (1.72)	0.22 (0.14)

* The fire area was extracted from Bergeron (2000), where the extent of fires that occurred in 1760 and 1944 is a minimum estimate since the fire burned a larger area than the investigated area. The species are balsam fir (Fir), white birch (Birch), white spruce (Spruce), trembling aspen (Aspen), and eastern white cedar (Cedar).

summarizes the study plots according to the fire-related characteristics (year of occurrence, burned area, and stand age at first sampling), number of transects, and sampling plots per fire and basal area (m² ha⁻¹) of the main tree species at the first sampling (i.e., 1991). In each plot, all trees greater than 5 cm in diameter at breast height (DBH), alive or dead (standing or fallen with bark and branches intact), were identified and classified in 5-cm DBH classes. In addition, within each plot, a 64 m² (8 m × 8 m) sub-plot and 12 one square meter (1 m × 1 m) micro-plots were installed to record saplings (stems ≤5 cm and height >1 m) and seedlings (stems ≤1 cm and height <1 m), respectively. All plots were located primarily on mesic or sub-hydric clay soils and occasionally on mesic loams.

2.3. SORTIE-ND Simulator

SORTIE-ND is a spatially explicit, individual-based forest dynamics model that simulates changes in tree populations over time [34]. Nothing in the model is pre-defined, default, or automatic, and the model acts, based upon processes, which are also termed behaviors. Behaviors control whatever happens during the simulation and are a combination of empirical and mechanistic processes. Empirical behaviors correspond to biological and environmental processes, such as seed dispersal, individual tree growth, and mortality, whereas mechanistic behaviors perform as helper functions to measure, calculate, and record forest metrics [34,38]. The model is designed to study neighborhood processes, in which trees are modelled individually. In SORTIE-ND, a large community of individual trees, which are classified according to their development stages (i.e., seedling, sapling, mature trees, or snags), represents the forest. Every single tree has a location and is described as a discrete object with several attributes (e.g., dimensions, age, species identify, growth rate; see Appendix A for attributes used by the model). During the simulation of stand development, SORTIE-ND sums up the spatio-temporal interactions of each tree with its surroundings to simulate the population-level dynamics of forest stands. These interactions include the effect of neighboring species and distances to other individuals on the growth of each individual, together with environmental conditions, such as light availability and substrate type.

SORTIE-ND functions are based upon a parameter file (see details in Appendix A), which compiles field data, local conditions, and a selected list of individual trees’ behaviors to make up the simulation in a defined period set as the number of timesteps. Timestep is the unit of time in SORTIE-ND and can have a length of one or more years. Once per timestep, each behavior, with a clearly defined action in a run, performs in a pre-defined order to structure the forest dynamics correctly. In the late 2000s, a parameter file was developed, tested, and modified for the LDRTF [38]. Several field experiments and studies attempted to either parameterize specific functions as various behaviors (e.g., light, tree allometry growth, and mortality) or modify the existing model parameters, for reconstructing the natural stand dynamics in the boreal mixedwoods of the LDRTF in a more realistic fashion [37,38,55,56]. To date, the LDRTF parametrized model is available for the five species that are presented in Table 1, jack pine and mountain maple (Acer spicatum Lamb.).

2.4. Simulation Runs

On the one hand, the observed post-fire stand composition of the study area was highly variable. On the other hand, the model assumes that only species that are present in each plot can grow, develop, and provide seeds. In order to account for the observed variability in the sampled plots, we created 431 starting conditions representing the empirical data that were collected from the 431 plots in 1991 (see

Table 2. Overview of simulation scenarios.

Source of Empirical Data Used for Setting Initial Conditions	Number of Initial Conditions *	Simulation Period	Purpose of Simulation
Short-Term Simulations
Empirical data in 1991 following spruce budworm outbreak	431	18 years	Model parameterization to reproduce empirical data in 2009.
Mid-Term Simulations
Empirical data in 1991 following spruce budworm outbreak	431	28 years	(1) Evaluating the mid-term effect of the recent budworm outbreak on stand dynamics. (2) Estimating stand composition and species proportion to define mortality rate of spruce budworm event for the long-term simulations.
Empirical data in 1991 prior to spruce budworm outbreak	431	28 years
Long-Term Simulations
Empirical data in 1991 following spruce budworm outbreak	431	60 years	(1) Evaluating the long-term performance of the model. (2) Evaluating the long-term effect of the recent budworm outbreak on stand dynamics.
Empirical data in 1991 prior to spruce budworm outbreak	431	60 years

* The tree population in the LDRTF parameter file was adjusted separately for each plot, giving a total of 431 starting conditions for simulation runs.

). The plot size (stem map) for simulation runs was set to 4 ha (200 m × 200 m), as the realistic minimum size that is recommended for SORTIE-ND simulations [34]. Since some behaviors that are included in the parameter file (e.g., competition mortality) could not handle timesteps longer than one year, the timestep was set to one year.

2.4.1. Model Parameterizations and Evaluations through Short-Term Simulations

Using empirical data that were collected in 1991, we adjusted the initial densities (number per hectare) of sapling and adult trees in 5-cm DBH classes and number of seedlings per hectare in the LDRTF parameter file (see Appendix A). We then carried out short-term simulations over 18 years, to reach the age of the empirical data in 2009. Consequently, we were able to evaluate the performance of the model by comparing simulation outputs with empirical data that were collected in 2009 in terms of stem size distribution and live stem basal area of tree species within each plot. We employed polynomial regression analyses to describe the general trends of species abundance, which were explained by their absolute basal area, in relation to the time since the fire. We also used non-parametric kernel density estimation to visualize the underlying distribution of the species basal area within 431 inventory plots and their corresponding simulated outputs, and controlled their equality of probability distributions with nonparametric two-sample Kolmogorov–Smirnov (K–S) tests. If the simulated outputs differed significantly (p-value < 0.05) from the empirical data for any species, we tried to identify various parameters that would have caused the discrepancies and modified the behavior functions accordingly. To adjust the LDRTF parameter file, we went through previous reports and published work for the LDRTF or similar sites and extracted new values for reparametrizing the behaviors that we found problematic. Simulation runs, model parameterizations, and validation analyses were repeated with the LDRTF-modified parameter file until the model could reconstruct the 18-year dynamics of our empirical data with good accuracy in terms of species compositions and diameter distributions (see Appendix B,

Table A1. The final list of modified parameters during short-term simulations.

Evaluations *	Parameter	Old Value	New Value	Source of New Parameters
Overestimation for birch >5 cm DBH	NCI maximum potential growth	0.692	0.593	Figure 7 in [8]
	DBH at onset of senescence	50 cm	30 cm	[100]
	Competition mortality shape parameter	4.33E-8	4.33E-6	[56]
	Competition mortality maximum parameter	0.1	0.25
	Adult stochastic mortality	0.00004	0.009
Underestimation for aspen 5–15 cm DBH	Seedling stochastic mortality	0.1	0.01
Underestimation for aspen >50 cm DBH	DBH at onset of senescence	30	40	[100]

* Evaluations show the significant differences (p-value >0.05) between 2009 empirical data and short-term simulation output. Old values are the model parameters at the beginning of simulations and model parameterizations. New values are the final adjusted values for model parameters.

, for parameter modifications).

2.4.2. Mid-Term Simulations

As previously mentioned, Bergeron et al. [44] reported a very severe spruce budworm outbreak for the period 1972–1987, which killed more than 56% (percentage of the total number of stems per hectare) of balsam fir and spruce trees. The first inventory in 1991 was conducted immediately after the maximum mortality event that was caused by the spruce budworm outbreak had taken place [49]. Thus, most of the trees that died during this outbreak could be easily identified since they were still standing (75%, on average) or present on the forest floor with bark and branches intact [44]. The inventory that was conducted in 1991 tallied dead trees (standing trees and fallen trees with bark and branches intact) greater than 5 cm DBH. We are confident that this approximation of stand density prior to spruce budworm outbreak is reliable. On the one hand, the natural mortality rate is generally low in the absence of spruce budworm [57]. On the other hand, the odds of including trees that were killed by intense suppression was decreased by considering trees greater than 5 cm DBH [44].

Table 3. Mean number (stems ha⁻¹, ±SE), basal area (m² ha⁻¹, ±SE), and mortality rate (Mort, %) of balsam fir and spruce prior to and following the recent spruce budworm outbreak (1972–1987).

Fire Date	Balsam Fir						White Spruce
	Number			Basal Area			Number			Basal Area
	Prior	Post	Mort	Prior	Post	Mort	Prior	Post	Mort	Prior	Post	Mort
1760	1560 (179)	489 (109)	68.65	15.16 (0.88)	2.60 (0.55)	82.85	67 (11)	43 (8)	35.82	2.77 (0.50)	2.04 (0.40)	26.35
1797	2303 (221)	1117 (192)	51.50	14.55 (1.04)	2.77 (0.42)	80.96	159 (24)	118 (20)	25.79	3.03 (0.42)	2.08 (0.39)	31.35
1823	2434 (211)	1485 (206)	38.99	13.87 (0.70)	4.32 (0.47)	68.85	209 (38)	188 (36)	10.05	2.90 (0.41)	2.46 (0.38)	15.17
1847	3289 (231)	1523 (174)	53.69	17.45 (1.06)	6.84 (0.50)	60.80	159 (16)	132 (14)	16.98	6.06 (0.65)	5.40 (0.58)	10.89
1870	1311 (132)	801 (95)	38.90	7.94 (0.61)	3.46 (0.36)	56.42	242 (27)	218 (25)	9.92	5.02 (0.61)	4.46 (0.56)	11.16
1916	1757 (144)	1319 (138)	24.93	7.72 (0.56)	4.54 (0.40)	41.19	254 (40)	246 (40)	3.15	2.04 (0.58)	1.99 (0.57)	2.45
1944	1272 (198)	1096 (174)	13.84	4.69 (0.66)	3.81 (0.48)	18.76	179 (50)	167 (47)	6.70	0.95 (0.23)	0.81 (0.19)	14.74

summarizes the number, basal area, and mortality rate of balsam fir and spruce trees prior to and following the 1972–1987 spruce budworm outbreak.

In order to account for episodic mortality that was caused by spruce budworm, we created another set of 431 parameter files by including tallied dead balsam fir and spruce trees with the living tree population. Eventually, having created two sets of parameter files for stands prior to and stands following the spruce budworm outbreak, we carried out two 28-year simulations (see Table 2) to reach the current age of post-fire stands in 2019.

2.4.3. Long-Term Simulations

For our long-term evaluations, we simulated the inventory data in 1991 for a 60-year period, i.e., for the data prior to the recent spruce budworm outbreak, by counting trees that were dead due to the recent spruce budworm outbreak as live trees. The recent outbreak occurred during 1972–1987; indeed, it was the only spruce budworm event between 1951 and 2019 [49,50]. Therefore, by simulating the stands with tree populations prior to the recent outbreak for a period of 60 years, we eventually considered a century of post-fire stand development during the 1951–2051 period, during which no spruce budworm outbreak occurred (see Figure 1). Note that this period of post-fire development without insect disturbance was based upon the initial conditions of post-fire stands, in which some stands had been affected by insect outbreaks prior to the 1972–1987 outbreak. We carried out another long-term (60 years) simulation considering the mortality caused by spruce budworm. To do so, we used the tree population of stands in 1991, following the recent (1972–1987) outbreak, to account for the spruce budworm outbreak during 1951–2019. Since the frequency of budworm outbreaks for the study region is about every 30 to 35 years [49], we expected the occurrence of another outbreak during 2019–2051. Thus, we added a spruce budworm mortality episode for fir and spruce (episodic mortality in SORTIE-ND) in 2020, i.e., 33 years after the end of the previous outbreak in 1987. The budworm-induced mortality rate was set based on the research of Blais [58] for spruce in the boreal region of eastern Canada and Bergeron et al. [44] for fir in the study region. In both aforementioned studies, the magnitude of the fir and spruce mortality rates depends upon the relative proportion of coniferous species to deciduous species; we used the model output after 28 years of simulation for stands following the outbreak to estimate the basal area proportion of species in the year prior to the budworm incident in 2020. Accordingly, depending upon the coniferous and deciduous relative (to total) basal area proportions, we set the mortality rate for fir in 5–10-, 10–15-, and ≥15-cm DBH classes and for spruce for trees greater than 10 cm DBH.

2.4.4. Assessing the Mid- to Long-Term Performance of the Model

For evaluating the mid- to long-term model performance, we used the simulation output over a 60-year period, accounting for spruce budworm incidents, at different lengths (timestep) of simulations, in which the age of younger stands after some years of simulation (up to 60 years) attained the age of older empirical stands in 1991 (see Table 1). Since the simulation accounted for two spruce budworm incidents, for the model evaluation (see Figure 1), we used empirical data prior to the outbreak counting for two outbreaks of three during the last century of stand development [49,50]. We then compared the species compositions and stem diameter distributions of empirical data and simulated stands representing stands of similar ages. Using the species relative (to total) basal area, we conducted principal component analysis (PCA) to visualize the potential links between fire year (stand age) and the composition of stands originating from each fire, in terms of the presence and dominance of species measured by their relative (to total) basal area. For each set of empirical and simulated data of the same age, we performed separate ordinations. The ordination axes were standardized to the first two components explaining the greatest percentages of total variance. Gaussian bivariate confidence ellipses enclosed 95% of the range of variability in stand composition and structure on the two axes. Furthermore, two-way mixed ANOVA that included a data (empirical and simulation) × species interaction was used to test whether the abundance (basal area m² ha⁻¹) of species in post-fire and simulated stands of the same age was different. We treated the species nested within the quadrats in each post-fire stand as random effects. The empirical stand originating from the 1944 wildfire was excluded from the analyses, since there was no available simulation of younger stands (than stand 1944) representing a similar age. However, the simulation output for the 1944 post-fire stand, for various simulation years, was used as simulated data for older stands.

2.4.5. Mid- to Long-Term Influence of Spruce Budworm Outbreak

We used two-way mixed ANOVA, including a data (prior to and following the outbreak) × species interaction, to test how the abundance (basal area m² ha⁻¹) of species in post-fire stands has changed over 28 years (mid-term) because of the recent outbreak. For the long-term effect of episodic coniferous mortality that was induced by spruce budworm, we compared the species composition and structure of the two simulated stands over a 60-year period. To do so, we performed polynomial regressions to describe the general trends of species abundance in relation to the time since the fire. The succession starts from 47, i.e., the age of the most recent post-fire stand (1944) at the time of first inventory in 1991, and ends at 291, the age of the oldest post-fire stand (1760) after 60 years of simulation. Given the species identities, basal area, and age of stands, we were able to illustrate the chronosequence of 291 years of successional dynamics of post-fire stands at the LDRTF.

Given that we were dealing with a large dataset, we utilized the Compute Canada (www.computecanada.ca) computing platform to implement the simulation runs and accelerate the process. We then transferred outputs to the R environment (version 3.5.0) for further data preparation and analyses [59]. We used the FactoMineR package to perform the ordinations, and ggpmsic and ggplot2 packages to perform the polynomial regressions and visualize the successional dynamics of species following the occurrence of wildfire.

3. Results

3.1. Model Evaluations

3.1.1. Short-Term Model Evaluation

The output of short-term simulation runs showed good agreement (p-value >0.05) with the empirical data that were collected in 2009 in terms of the live stem density, live stem basal area of the study species, and the distribution of the live stem basal area in the 5-cm diameter classes.

Table 4. Comparison of species basal area (m² ha⁻¹, ±SE) for empirical data in 2009 and short-term simulation output.

Fire Date	Species Mean Basal Area for Standing Live Sapling and Adult Trees
	Fir		Birch		Spruce		Aspen		Cedar
	Emp.	Sim.	Emp.	Sim.	Emp.	Sim.	Emp.	Sim.	Emp.	Sim.
1760	6.49 (0.81)	5.61 (0.88)	4.88 (0.67)	4.03 (0.68)	3.28 (0.61)	2.37 (0.43)	3.45 (1.21)	2.21 (0.68)	15.73 (1.89)	13.80 (1.77)
1797	10.38 (0.79)	8.43 (0.86)	7.40 (0.79)	5.88 (0.54)	4.53 (0.78)	3.04 (0.52)	0.73 (0.34)	0.67 (0.24)	5.11 (0.93)	5.09 (0.81)
1823	11.93 (0.74)	10.54 (0.84)	3.07 (0.34)	3.63 (0.35)	4.35 (0.60)	3.84 (0.50)	9.99 (1.30)	10.74 (1.25)	0.86 (0.40)	0.71 (0.35)
1847	11.36 (0.60)	11.97 (0.58)	2.58 (0.32)	3.88 (0.35)	6.47 (0.73)	6.85 (0.67)	13.67 (1.15)	11.56 (0.80)	0.63 (0.22)	0.85 (0.32)
1870	6.69 (0.64)	6.42 (0.62)	3.72 (0.53)	4.33 (0.46)	5.52 (0.74)	6.10 (0.71)	11.56 (1.33)	15.33 (1.72)	0.04 (0.04)	0.03 (0.03)
1916	9.13 (0.64)	7.98 (0.59)	5.51 (0.64)	6.52 (0.72)	3.33 (0.56)	2.56 (0.59)	17.80 (2.28)	16.20 (1.86)	1.37 (0.77)	1.22 (0.46)
1944	6.60 (0.62)	6.26 (0.71)	8.04 (0.78)	8.01 (0.68)	1.71 (0.36)	1.28 (0.31)	11.50 (1.66)	12.20 (1.78)	0.33 (0.25)	0.27 (0.17)

Species are balsam fir (Fir), white birch (Birch), white spruce (Spruce), trembling aspen (Aspen), and white cedar (Cedar). Bold values indicate that basal area in the empirical dataset is significantly (p-value <0.05; two-sample K–S test) different from the basal area in the simulated data set.

summarizes the live basal area for the empirical data that were collected in 2009 and the 18-year simulated output. As presented, the empirical and simulated basal areas of the study tree species are statistically similar (p-value >0.05), except for birch in stands originating from wildfire in 1847, with significant basal area overestimation (p-value <0.05). Since the DBH limit for saplings in the LDRTF parameter file (stems ≤3–10 cm depending on tree species) differed from the field data (stems ≤5 cm), saplings and adult trees were considered as one group for further analyses that are described throughout the results section.

Figure 2 illustrates the dynamics of the species basal area over 184 years, from the youngest post-fire stand with an age of 65 years to the oldest post-fire stand with an age of 249 years, for both empirical data in 2009 and model output at the end of the 18-year simulation. As demonstrated here (Figure 2), the model could successfully capture the changes in the stem basal area of the study species within the 95% confidence interval.

In addition (Table 4), the distribution of the species basal area among the empirical plots is not significantly different (p-values >0.05; K–S test) from the simulated plots, except for birch trees in stands originating from the fire in 1847 (p-value 0.025), which was slightly overestimated (see also the kernel density curves in Appendix C, Figure A1).

To account for tree size distributions in terms of the stem number and basal area per hectare, we also examined the distribution of the basal area in 5-cm diameter classes for standing live trees for the empirical data and model output after an 18-year period of simulation (Figure 3). The 18-year simulation proved efficient in reconstructing the tree species that were distributed in different size classes. Yet, in some fire areas, the simulation projected a higher basal area/stem number of birch (DBH < 10 cm), aspen (DBH < 20 cm), and cedar trees (DBH < 5 cm) (p-values <0.05).

3.1.2. Mid- to Long-Term Model Evaluation

When comparing the simulated stands with the empirical data of the same ages, the mid- to long-term simulation of post-fire stands could faithfully reconstruct the species composition and diameter distributions of the empirical data (Figure 4); however, in some cases, the simulated basal area (m² ha⁻¹) was significantly different from the empirical basal area (p-values <0.05 in

Table 5. Basal area comparisons of empirical stands that were surveyed in 1991 (prior to the spruce budworm outbreak) with the simulated stands of similar ages (with spruce budworm outbreak events).

Empirical Data Fire (Stand Age)	Simulation Data
	Fire Year + Length of Simulation (years)
1760 (231 y)	1797 + 37 (Total−)
	Fir− Birch Spruce Aspen Cedar
1797 (194 y)	1823 + 26	1847 + 50 (Total)
	Fir Birch− Spruce+ Aspen+ Cedar−	Fir Birch− Spruce+ Aspen+ Cedar−
1823 (168 y)	1847 + 24 (Total+)	1870 + 47 (Total)
	Fir Birch Spruce+ Aspen+ Cedar	Fir−Birch−Spruce+Aspen+ Cedar
1847 (144 y)	1870 + 23 (Total)
	Fir− Birch Spruce Aspen Cedar−
1870 (121 y)	1916 + 46 (Total)
	Fir+ Birch Spruce− Aspen Cedar+
1916 (75 y)	1944 + 28 (Total)
	Fir Birch Spruce Aspen Cedar

The significant differences detected by two-way ANOVA (p-value <0.05) are in bold, where + and − superscripts refer to significant (p-value <0.05) overestimation and underestimation of the model, respectively. Total refers to the total basal area of stands, regardless of tree species.

).

As presented in Table 5, the total basal area of all tree species trees within the plots was underestimated (p-value <0.05) for stands that are aged 231 years, and overestimated for stands at an age of 168 years (post-fire stand 1847 after 24 years of simulation). In some stands that were older than 121 years, the model underestimated (p-value <0.05) the basal area proportion of fir, birch, and cedar, and overestimated (p-value <0.05) the basal area proportion of spruce and aspen. In younger stands, there is some evidence for overestimation (p-value <0.05) for the fir and cedar basal area and underestimation for the spruce basal area (Table 5).

Figure 5 illustrates the variability and composition of both empirical and simulated stands of equal ages. The first and second axes of the PCA explained 51% to 70.7% and 16.7% to 25.9% of the total variance, respectively. The model appeared to reproduce the variability that was found in post-fire stands’ composition and structure with a good convergence. The ordinations seem to be driven by the higher abundance of hardwoods in early successional post-fire stands, rather equal abundances of hardwoods and softwoods in mid-successional post-fire stands, and the higher abundance of coniferous species in late-successional post-fire stands. Simulated data show, more or less, greater variability (see the relative size of the ellipses in Figure 5) than the empirical data of the same ages, more specifically in the middle-aged to old post-fire stands.

3.2. Mid- to Long-Term Influence of Spruce Budworm

Except for the youngest post-fire stand, the basal area proportion of balsam fir during 28 years of development in post-fire stands was significantly (p-value <0.05) less for simulated outputs accounting for the most recent spruce budworm outbreak during 1972–1987, compared to simulated outputs without the spruce budworm effect (

Table 6. Comparison of species basal area (m² ha⁻¹, ±SE) for 1991 empirical data including or not mortality induced by the spruce budworm outbreak after mid- and long-term simulations.

Fire Date	Simulation Length (years)	Species Mean Basal Area for Standing Live Sapling and Adult Trees
		Fir		Birch		Spruce		Aspen		Cedar
		No-SBW	SBW	No-SBW	SBW	No-SBW	SBW	No-SBW	SBW	No-SBW	SBW
1760	28	23.69 (0.96)	7.98 (1.09)	2.92 (0.61)	3.64 (0.65)	3.36 (0.56)	2.68 (0.48)	1.57 (0.52)	1.94 (0.58)	15.06 (1.91)	16.87 (2.07)
	60	32.97 (1.26)	9.25 (1.06)	1.26 (0.4)	2.92 (0.56)	3.54 (0.66)	3.93 (0.78)	0.53 (0.17)	1.26 (0.35)	19.46 (2.22)	25.88 (2.63)
1797	28	25.73 (1.01)	12.65 (1.15)	4.52 (0.52)	5.26 (0.53)	4.89 (0.65)	3.72 (0.61)	0.56 (0.20)	0.70 (0.26)	5.71 (0.92)	6.52 (1.02)
	60	35.74 (1.05)	15.02 (1.24)	1.93 (0.39)	3.53 (0.55)	6.7 (0.89)	5.57 (0.93)	0.49 (0.2)	0.82 (0.37)	8.24 (1.27)	12.04 (1.74)
1823	28	23.85 (0.96)	14.77 (1.04)	3.16 (0.32)	3.35 (0.33)	5.18 (0.60)	4.73 (0.57)	10.71 (1.31)	11.35 (1.29)	0.80 (0.39)	0.88 (0.42)
	60	32.68 (1.18)	18.72 (1.09)	1.59 (0.26)	1.91 (0.28)	7.17 (0.80)	6.07 (0.68)	7.40 (1.06)	8.76 (1.03)	1.13 (0.52)	1.52 (0.66)
1847	28	23.60 (0.73)	15.56 (0.68)	3.13 (0.29)	3.47 (0.32)	8.47 (0.78)	7.71 (0.72)	10.76 (0.86)	11.14 (0.84)	0.90 (0.33)	0.97 (0.35)
	60	31.79 (1.03)	20.46 (0.86)	1.32 (0.17)	1.82 (0.27)	9.31 (0.85)	7.24 (0.66)	6.25 (0.71)	7.33 (0.7)	1.28 (0.46)	1.72 (0.59)
1870	28	14.11 (0.95)	8.68 (0.83)	3.58 (0.39)	3.72 (0.40)	7.71 (0.86)	7.20 (0.83)	15.72 (1.76)	15.80 (1.74)	0.03 (0.03)	0.03 (0.03)
	60	23.29 (1.32)	12.68 (1.13)	1.60 (0.2)	1.73 (0.23)	11.3 (1.15)	9.05 (1.01)	11.41 (1.32)	12.41 (1.36)	0.04 (0.04)	0.05 (0.05)
1916	28	14.29 (0.85)	10.70 (0.72)	6.42 (0.71)	6.47 (0.71)	3.12 (0.63)	3.09 (0.63)	14.67 (1.77)	14.73 (1.76)	1.43 (0.54)	1.60 (0.61)
	60	24.13 (1.09)	15.04 (0.91)	4.90 (0.57)	5.12 (0.58)	5.67 (0.84)	4.62 (0.74)	6.99 (1.15)	7.40 (1.15)	2.33 (0.82)	3.01 (1.07)
1944	28	9.08 (1.02)	8.22 (0.88)	7.86 (0.64)	7.93 (0.65)	1.86 (0.45)	1.65 (0.41)	12.11 (1.73)	12.23 (1.73)	0.34 (0.22)	0.34 (0.22)
	60	16.69 (1.50)	11.75 (1.17)	6.19 (0.54)	6.40 (0.55)	3.47 (0.79)	2.49 (0.61)	8.26 (1.27)	8.64 (1.28)	0.57 (0.36)	0.67 (0.41)

Species are balsam fir (Fir), white birch (Birch), white spruce (Spruce), trembling aspen (Aspen), and white cedar (Cedar). The significant differences that were detected by two-way ANOVA (p-value <0.05) are in bold. SBW refers to spruce budworm outbreak.

). The reduction in the balsam fir basal area varied from 9% in the youngest stand, which originated from the wildfire in 1944, to 66% in the oldest stand, which originated from the wildfire in 1760. Despite the decreased basal area of white spruce, the influence of the recent spruce budworm outbreaks did not appear statistically significant (Table 6).

To summarize the successional dynamics of species in post-fire stands of the study area, with and without spruce budworm outbreaks, we regressed the changes in the species basal area against the time since the fire for the seven post-fire stands (Figure 6). The proportion of the species basal area significantly changed through a gradual transition of post-fire stands from early stages to late-successional stages (p-value <0.05). As illustrated in Figure 6, the relative proportion of hardwoods (aspen and birch) is higher in early years following a wildfire and it decreases with the time since the fire. The occurrence of a spruce budworm outbreak, as the secondary disturbance for post-fire stands, increases hardwood proportions, specifically in old-growth post-fire stands (see significant increases in Table 6). Throughout succession, the basal area proportion of coniferous species (fir, spruce, and cedar) increases. Around the mid-successional stages, they dominate post-fire stands (see Table 6). The mortality that was induced by two spruce budworm incidents significantly (p-value <0.05) decreased the basal area proportion of fir in all post-fire stands. In contrast, spruce was not meaningfully affected, given that it subsequently recovered its growth and dominance status prior to the outbreak (60-year simulation output in Table 6). The spruce budworm outbreak set back the successional dynamics of the middle-aged post-fire stands, by decreasing the proportion of conifers and consequently increasing the relative abundance of hardwoods. Likewise, due to the high mortality of balsam fir, the dynamics in stands in which cedar were already established accelerated the dominance of cedar in post-fire stands around the age of 250. Within undisturbed stands, balsam fir is the most important species in terms of the basal area.

4. Discussion

We verified whether model parameterization allows reconstruction of natural post-fire successional dynamics of our study area over various time steps by comparing model outputs to field observations over three simulation horizons. As we discuss below, the short-term simulation proved competent in reproducing the species composition and size distribution that was observed in the field data. Moreover, using the mid- and long-term simulations while accounting for episodic spruce budworm outbreaks made it possible to examine the model’s potential to generate the species compositions and size distributions of empirical data that characterized old-growth forest. It also enabled us to assess the potential changes in the successional dynamics of post-fire stands through coniferous mortality that were induced by spruce budworm as a secondary disturbance, which reoccurs episodically in the study region.

4.1. Short-Term Evaluation

Short-term simulation outputs showed good agreement with the empirical data in terms of the species composition, species basal area, and species density. The empirical and simulated data were generally similar, despite a slight basal area overestimation of white birch for the stands originating from the 1847 fire, higher basal area/stem number of trembling aspen with diameters smaller than 20 cm, and high ingrowth of small cedar trees (DBH < 5 cm). Various factors that affected the recruitment and survival of species, but which were not incorporated into the model, may explain the significant differences we detected. For instance, the model neither considered the reported defoliation of aspen and birch trees in 2001, nor the dry summers in 2001 and 2002 [38,60]. The 2001 drought was severe, which increased tree mortality with varying effects on species, depending upon their vulnerability to drought [61]. Consequently, it can be reasonably expected that simulations deviated from empirical data for some species, especially if the model is not simulating such drought events. The model also did not account for any other sources of disturbance (i.e., major windthrow events) that might be important, although less frequent than insect outbreaks [62,63]. These aforementioned disturbances could possibly create small gaps, which are more suitable for birch regeneration compared to aspen, which requires full sunlight for establishment. Nevertheless, the emergence of a few aspen suckers may help maintain aspen (co)dominance [32]. Although both aspen and birch exist along with conifers in the mid-successional stands (i.e., 1823, 1847, and 1870 post-fire stands) [31], birch is also expected to have lower recruitment than aspen [12,64,65]. Specifically, birch has the lowest proportion of basal area in the 1847 post-fire stands (16% of the total basal area in the 1847 post-fire stand). The collective effects of low proportions of birch trees, defoliation by forest tent caterpillar, and probable windthrow likely reduced the recruitment of small birch trees into the larger diameter classes in the 1847 post-fire stands, which were not addressed by the model. Likewise, the higher proportion of small birch trees (DBH < 10 cm) that were observed in the field, when compared to simulation outputs, might be explained by the aspen mortality that was induced by FTC defoliation, dry summers, or windthrow. In the absence of aspen, birch possibly has enhanced sprouting or seeding [8], and stronger recruitment in gaps that were created by aspen mortality [16,27].

The model simulated high ingrowth of small cedar trees (DBH < 5 cm), mainly in old stands where white cedar establishment dominates. Yet, on the one hand, the model did not account for the competitive effects of herbs and woody shrubs, except for mountain maple, on the growth and survival of seedlings. On the other hand, the death of companion species (aspen, birch, spruce, and fir) in old-growth stands and their fast decay rates following their death [66] would provide suitable conditions for cedar recruitment. Thus, a small degree of cedar overestimation can be the consequence of a high number of cedar trees in old stands that were dispersing seeds on decaying wood, as their most favorable germination substrate [23,67,68], and growing in small canopy openings that would facilitate their establishment and growth [69] in a less competitive environment. However, constructing a model with proper behaviors defining the growth and mortality of adult white cedar trees would properly control the recruitment of cedar trees into the larger diameter classes. Despite the aforementioned contradictions between the field observations and simulation outputs in terms of the recruitment and survival of birch, aspen, and cedar, the overall simulated dynamics of stands were consistent with the observed dynamics of the empirical data, which we discussed in the following section.

Individual Species Dynamics within the Post-Fire Stands

The post-fire stands in the study represent a chronosequence of 249 years of forest succession originating from seven wildfires. Within each fire area, the successional dynamics of species can be appropriately explained by their life-history traits, such as longevity, mode of regeneration, growth rate, and shade tolerance [8], together with inhibition or facilitation processes [8,70]. In utilizing suitable parameters for species life-history traits and other individual- and stand-level processes, we could successfully reproduce post-fire species and stand dynamics in the boreal mixedwoods of LDRTF that were presented earlier by Bergeron [8]. In his study, Bergeron [8] pointed out three important developmental stages that characterize the post-fire stand dynamics of the study region: The post-fire cohort is dominated by hardwoods with a conifer understory; first, the aspen cohort’s breakup is followed by subsequent cohorts of aspen mixed with conifers that are recruited into the canopy; and finally, stands are dominated by coniferous species undergoing spruce budworm outbreaks. The relative importance of each species in post-fire stands is controlled more by the initial numbers than by different niches [71] and the model properly benefited from including the observed initial numbers of species in the field data to simulate the stand dynamics following the wildfires.

The young stands originating from wildfires in 1944 and 1916 formed the early successional post-fire hardwood cohort (also called the aspen cohort), which are dominated by fast-growing trembling aspen and white birch [2,72]. Following the wildfires, massive aspen root suckers colonize burned areas, together with the resprouting of white birch from the stem collars or from seeds dispersed over long distances [2,8]. Similar to Bergeron [8], simulations displayed a decreasing importance of trembling aspen and white birch in all post-fire stands throughout succession: The highest basal area values for aspen were found in the first two cohorts, and for birch, in the first and last cohorts. Indeed, due to its great longevity (exceeding 230 years of age), birch can remain suppressed under the aspen canopy and grow slowly during the period of aspen dominance [8]. Although less abundant, where their seeds are available, slow-growing balsam fir and spruce appear in the post-fire cohort with a 5- to 10-year delay compared to hardwood species [8,73]. Moreover, the very low constancy of white cedar that was observed in post-fire cohorts is possibly the result of its low seed dispersal capability [74] or the absence of suitable seedbeds following fire [8]. As highlighted by Bergeron [8] and seen in both the field and simulated data, variation in the species composition of the post-fire cohort has great influence on the successional dynamics of stands from hardwood trees to conifer regeneration and establishment.

The breakup of the aspen cohort comes in the middle-aged stands, originating from the 1823, 1847, and 1870 wildfires, where declining proportions of hardwoods, particularly aspen, and an increasing occurrence of balsam fir and spruce form mid-successional mixedwoods. Gaps that are created in the canopy due to the breakup of the aspen cohort are filled either through self-replacement, in which a new aspen cohort re-establishes from successful root suckering [8,14,73], or by birch, fir, and spruce, which have already established beneath the aspen canopy [8,30,75]. The continuous occurrence of fir throughout the succession is the consequence of fir’s shade tolerance, which allows for abundant regeneration under closed canopies or within small gaps [72], which was reflected appropriately in simulations of the process. As Bergeron [8] suggested, the limited recruitment of spruce in undisturbed soil following canopy disturbances is responsible for the low occurrence of spruce during succession after the initial post-fire cohort, which was also adequately addressed by the model. Indeed, spruce establishment on organic soil or on substrates other than mineral soil is likely limited because of its small seeds [76]. Since the model used the number of surviving spruce and fir after spruce budworm outbreak, the short-term simulation precisely reflected the status of these two species in the 2009 field observations.

Eventually, the old growth stands originating from wildfires in 1797 and 1760 form the old successional coniferous cohort with a noticeable occurrence of white cedar trees. The simulated data reveal the dominance of eastern white cedar in old-growth forests some decades after aspen cohort breakup, which accords with the findings of Bergeron [8]. Some studies have related a higher cedar abundance to the decreased inhibitory effects of hardwood litter, especially aspen litter [67,76], to the increased amount of moss layer during succession [67,77,78,79], and to intense layering under the conifer canopy in more humid conditions [8]. The model successfully reproduced cedar longevity [80] and shade tolerance [74], as well as favorable conditions for cedar establishment. These conditions insured the maintenance and increasing abundance of white cedar with time since the fire. In the coniferous cohort, aspen abundance and its relative proportion decreased noticeably. Bergeron [8] related this decline to reduced root sucker production with time. Aspen root sucker emergence or growth might be weakened in small gaps receiving insufficient heat and sunlight [81], an increased humus layer inhibiting aspen suckering [82], and in stands with high competition from other strongly established species [8]. Through the actions of various behaviors, such as light, neighborhood competition, and substrate (see Appendix A), the simulator managed to obtain recruitment of aspen that was appropriate in different developmental stages of post-fire stands.

4.2. Mid- to Long-Term Evaluation—Can We Reconstruct the Post-Fire Forest Succession?

Simulating the development of post-fire stands of younger ages to reach the age of older post-fire stands illustrated the ability of the model to reconstruct the species composition and stem size distributions of not only young- to mid-successional stands but also old-growth post-fire stands. The differences that were seen in the stand variability and species abundances in the observed and simulated data might be the cumulative effect of the secondary disturbance history, such as insect outbreaks and climatic disturbances. Different disturbance regimes with various patterns fluctuating in space and time have been shown to be the result of different forest compositions, thereby promoting species coexistence [83]. The empirical post-fire stands have passed through various disturbances for which the model was not intended (i.e., not constructed to include the relevant information), or those that it could not address accurately. For instance, older stands have suffered several conifer-defoliator spruce budworm [49] and hardwood-defoliator tent caterpillar [53] outbreaks, which reduced the abundance of host species during the natural successional dynamics of stands. We simulated the spruce budworm mortality of balsam fir and white spruce, but the adjusted mortality rates were based upon the most recent documented outbreaks [44,58]. The simulations were probably different from the real situations, when post-fire stands experienced episodic outbreaks during their developments. During actual stand development, a spruce budworm outbreak and its consequences (i.e., mortality and weakened growth of the host species) usually occur over the course of time [49,50] while the model only simulated mortality that was induced by spruce budworm, immediately and within one year of post-fire stand development. On the one hand, comparing the stands originating from different wildfires of the same age at various spatial locations increases the risk of incurring different seed dispersions from the neighboring stands beyond the plots, because the model does not allow for interactions among plots and their neighborhoods [34]. On the other hand, the uncertainty is increased regarding having similar initial stand compositions, which is a decisive factor in successional post-fire stand development [8,71] and for simulation runs [34]. Dissimilar initial stand compositions are the result of pre-fire stand compositions that frequently control post-fire stand composition [84,85].

4.3. Spruce Budworm Contribution to the Successional Dynamics of Post-Fire Stands

Budworm-induced fir mortality largely governs post-fire stand dynamics when the abundance of host species defines the severity of a budworm outbreak [44,58]. Using mid- to long-term simulations, we tested the effects of spruce budworm outbreak on species and stand dynamics over a century of post-fire stand development. The model illustrated a loss in balsam fir basal area, which is the lowest in young post-fire stands and the greatest in the oldest post-fire stands, both in the mid-term and over the long term. Since budworm-induced mortality is supposed to be more severe in mature than in immature stands [57,58], and the abundance of mature balsam fir gradually rises in mid- and late-successional post-fire forests, increasing fir mortality caused by budworm with time since the fire is expected. Following outbreaks, the basal area proportion of white birch and trembling aspen gradually increases, which is most evident in old-growth post-fire stands. Several factors define the recruitment of shade-intolerant birch and aspen into the older post-fire stands following the budworm outbreak. Among these factors are defoliation severity and the degree of subsequent canopy mortality, canopy composition, and density of coniferous advanced regeneration [8,27,86,87]. The simulation produced a greater post-disturbance increase in the aspen basal area compared to white birch. Indeed, fast-growing aspen root suckers can rapidly fill larger gaps that are created by budworm outbreaks, whereas birch and fir recruit into smaller gaps in the absence of competition from aspen [30,72].

Compared to the high mortality rates of balsam fir that are incurred through spruce budworm outbreaks, mortality was low for white spruce [44]. Spruce mortality mostly occurs in the two last cohorts, more specifically in middle-aged post-fire stands. This difference between the two species could be the consequence of a smaller proportion of spruce relative to fir, the faster emergence of buds in fir [88], the greater ability of spruce to tolerate defoliation, and budworm preferences for balsam fir [89]. Accordingly, despite the limited canopy recruitment of white spruce, mainly in the post-fire cohort, and the occurrence of spruce budworm outbreaks, spruce survives to an old age [8] and exists throughout the simulated succession of post-fire stands for over 291 years. Similar to observations that were made by Bergeron [8], white spruce’s relative proportions in empirical post-fire stands decreased throughout the succession following the outbreak. Yet, the model produced a higher abundance of white spruce in late-successional post-fire stands compared to the observed data. Likewise, an increase in white spruce abundance, following spruce budworm outbreaks, has been reported in other studies in Ontario and Quebec [90]. In his study, Bergeron [8] stated that old stands had yet to recover from large gaps that had formed in the fir canopy, and which were created by the 1972–1987 budworm mortality, Further, the length of the chronosequence was not sufficient to predict future stand development. The increased basal area proportion of white spruce in old-growth stands could be the consequence of several factors: Low seedling mortality, longevity greater than the maximum age of all post-fire stands, no signs of higher mortality in the oldest stands, and an ability to attain the highest canopy position [8]. Moreover, increased abundance of well-decomposed deadwood in old stands that is due to budworm mortality can facilitate white spruce establishment [67,68]. These factors, together with the increased growth of residual trees following the reduction in stand density and competition [35,91], not only secure the continuous recruitment of white spruce in old growth post-fire stands but also cause an increase in the white spruce basal area following the outbreaks.

The long-term simulation also suggested an increased basal area for eastern white cedar in post-fire stands disturbed by budworm outbreaks. Compared to balsam fir and white spruce, white cedar is generally a sub-canopy tree [8]. Like white spruce, cedar has a great longevity and an increase in the deadwood component of the forest floor favors its establishment [8,23]. The greater persistence of white cedar post-disturbance explains the exceptional cedar longevity, increased deadwood components due to budworm mortality, and release from the surrounding fir and spruce, occurring mostly in small canopy openings [69]. Openings in the canopy that result from defoliation and tree mortality also increase the growth of advanced balsam fir regeneration, resulting in cyclical fir replacement [28,92]. Altogether, as shown by the long-term simulations and highlighted by Bergeron [8], budworm outbreak did not reset the succession back to the hardwood stage, although it permitted some recruitment of trembling aspen and white birch.

Simulations for both disturbed and undisturbed old-growth post-fire stands had projected dominance of balsam fir and white cedar, as was observed by Kneeshaw and Bergeron [93]. As discussed earlier, however, the high mortality of balsam fir and increased abundance of white cedar resulted in a more rapid dominance of balsam fir by cedar. It must be noted that the study area did not experience any windthrow immediately following the spruce budworm outbreaks [8].

5. Conclusions

The current study took advantage of short- to long-term simulations to reconstruct the successional dynamics of post-fire stands, and examined the influence of defoliation by spruce budworm outbreaks, over a century-long period of post-fire stand development. Using SORTIE-ND, we were able to reasonably reproduce the compositional variability of LDRTF post-fire stands, which had been previously identified by chronosequence analysis [94] and that were confirmed by a stand-reconstruction approach [8]. This approach afforded us the opportunity to use SORTIE to predict responses to specific silvicultural treatments, especially over the short term, where the model appeared to make good predictions. In other words, the model apparently has great potential to evaluate alternative management or disturbance scenarios. Despite the importance of understory vegetation, which may compete with tree species regeneration for resources, the model only took into account the effects of mountain maple on the simulation processes [93]. Therefore, we were unable to verify the influence of herbs and woody shrubs, particularly pin cherry (Prunus pensylvanica L.f.) and willows (Salix spp.), which were considered in the study conducted by Bergeron [8]. The simulator did not correctly predict the density of seedlings for the study species; however, it showed its ability for reconstructing successional dynamics in the eastern Canadian boreal mixedwood forest, following some model adjustments (i.e., survival and recruitment of seedlings to saplings and adult trees). As understory competition (e.g., grasses, forbs, or shrubs) may limit seedling recruitment [95,96], implementing SORTIE-ND with some information about their competitive abilities could improve simulations of seedling abundance. Moreover, the model’s efficiency could be improved through appropriate quantification of seed production, dispersal, and germination, as well as seedling survival rates on different seedbeds in LDRTF and mortality at different stages.