Evaluating the Performance of a Forest Succession Model to Predict the Long-Term Dynamics of Tree Species in Mixed Boreal Forests Using Historical Data in Northern Ontario, Canada

: Environmental concerns and economic pressures on forest ecosystems have led to the development of sustainable forest management practices. As a consequence, forest managers must evaluate the long-term effects of their management decisions on potential forest successional pathways. As changes in forest ecosystems occur very slowly, simulation models are logical and efﬁcient tools to predict the patterns of forest growth and succession. However, as models are an imperfect representation of reality, it is desirable to evaluate them with historical long-term forest data. Using remeasured tree and stand data from three data sets from two ecoregions in northern Ontario, the succession gap model ZELIG-CFS was evaluated for mixed boreal forests composed of black spruce ( Picea mariana [Mill.] B.S.P.), balsam ﬁr ( Abies balsamea [L.] Mill.), jack pine ( Pinus banksiana L.), white spruce ( Picea glauca [Moench] Voss), trembling aspen ( Populus tremuloides Michx.), white birch ( Betula papyrifera Marsh.), northern white cedar ( Thuja occidentalis L.), American larch ( Larix laricina [Du Roi] K. Koch), and balsam poplar ( Populus balsamefera L.). The comparison of observed and predicted basal areas and stand densities indicated that ZELIG-CFS predicted the dynamics of most species consistently for periods varying between 5 and 57 simulation years. The patterns of forest succession observed in this study support gap phase dynamics at the plot scale and shade-tolerance complementarity hypotheses at the regional scale.


Introduction
Emulating natural disturbances in managed forests has been recommended as a means to sustainably manage boreal forests for timber production, climate regulation and climate change mitigation [1], biodiversity [2], and many other ecosystem services that are essential to human well-being [3][4][5]. These multiple social and environmental demands have led to increasing pressure on forest ecosystems. Consequently, forest managers must be concerned with the application of the basic principles of ecosystem management and integrity, which requires the evaluation and prediction of potential successional pathways [6,7]. Here, succession is defined as the directional change through time of species composition and vegetation physiognomy of a single site where climate remains relatively constant [8]. Bazzaz [9] introduced the notion of scale into the discussion of succession. In boreal forests, evidence points to gap dynamics at the local scale and species succession at the landscape scale.
For boreal forests, a growing body of evidence [10][11][12][13][14][15] suggests that emulating natural disturbances in managed forests should be performed in the context of forest gap dynamics theory [16]. Forest gap dynamics is a process that occurs when one or a few canopy trees die

Study Area
This study was conducted using historical data obtained from the boreal forest region of northern Ontario in ecoregions 3E and 3W in the middle section of the boreal forest of the Ontario Shield Ecozone (Figure 1). Details on ecological characteristics may be found in [38]. A brief summary follows. The area of this ecozone is approximately 65 millions hectares (ha) with variable topography and its climate is characterized by long, cold and moist winters and short warm summers. Variations in temperature and precipitations are considerable. In average, daily temperatures are approximately −15 and 17 • C in January and July, respectively. As the ecoregions 3E and 3W differ in basic ecological characteristics, results are presented separately. Ecoregion 3E is in the Humid Mid-Boreal ecoclimatic region. Precipitation varies between 652 and 1029 mm year −1 . Its deposits consist of morainal clays in the northeastern part, while ground moraine is generally found in most other areas. Soils consist mostly of podzols. The ecoclimatic region of ecoregion 3W is moist mid-boreal with lower annual precipitations than ecoregion 3E, varying from 654 to 879 mm. The geology of this ecoregion consists of granitic bedrock, but also includes basalt and volcanic rocks, greenstone, siltstone and shale. Soils consist of brunisols and podzols. The following tree species were found in both ecoregions: black spruce (Picea  Table 1). The importance of species composition differed between ecoregions 3E and 3W. Mean basal area was higher in Ecoregion 3E than 3W for black spruce, trembling aspen, American larch, and balsam poplar, but the opposite for jack pine and white birch. Differences in mean basal area between both ecoregions were relatively small for balsam fir, white spruce and northern white cedar. For the range of variation between minimum and maximum values, differences between the ecoregions were pronounced for white spruce, trembling aspen, white birch, northern white cedar, American larch and balsam poplar.

Historical Forest Data Sets
The individual-tree sample plot data used in this study to compare model predictions against observations were obtained from three data sets maintained by the Ontario Ministry of Northern Development, Mines, Natural Resources and Forestry (NDMNRF): American Can, Kimberly Clark and Spruce Falls Power and Paper Co. (SFPP), which included 148, 115 and 113 sample plots, respectively. Measurements were conducted between 1951 and 2000 for American Can, between 1952 and 2000 for Kimberly Clark and between 1930 and 1989 for SFPP. All the trees in each sample plot were measured for diameter at breast height (dbh). For all three data sets, the combination of sample plots and remeasurement periods resulted in tree measurements that exceeded 1 million records (Table 2). For all three data sets, stand age varied between 5 and 214 years. Most sample plots were remeasured every 5 years (Table A1). Remeasurement intervals between 1 and 4 years and between 6 and 49 years were less frequent. The number of sample plots in different age classes varied substantially (Table 3). Relatively high numbers of sample plots were found in the age classes from 30 to 39 to 90 to 99 years. The duration of remeasurements generally varied from a minimum of 5 years to a maximum of 56 years. Most of the maximums were close to 50 years. Sample plot size varied between 400 and 4047 m 2 (Table A2). The majority of them measured 809 m 2 in size (30.5%). While the sample plots for the American Can data set were in ecoregions 3E and 3W, the sample plots for Kimberly Clark and SFPP were in 3W and 3E, respectively ( Table 4). The sample plots were established on sites of natural origin and soils that had sandy, sandy-loam, loam, silty-loam or clayey-loam textures. However, for the SFPP data set, 21 sample plots originated from stands that were cut or cut and burned.   Maximum  5 to 19  10  10  56  20 to 29  19  5  51  30 to 39  50  20  55  40 to 49  27  5  49  50 to 59  44  5  57  60 to 69  37  10  56  70 to 79  39  5  48  80 to 89  70  5  56  90 to 99  34  5  47  100 to 109  15  5  48  110 to 119  8  15  48  120 to 129  5  15  48  130 to 139  7  18  48  160 to 169  10  10  48  ≥170  1  15 15 Variation in stand density and dbh was considerable for most species (Table 2). For stand density, substantial variations were evident for black spruce, balsam fir, trembling aspen, jack pine and white birch in the three data sets, but were less pronounced for the other species. Low species density values observed may be related to the mixed nature of the stands. Only 2.3% of the stands contained 1 species, which consisted of black spruce, but 68.2% of the stands included between 4 and 6 species. Approximately one-quarter of the stands were composed of 2 or 3 species (22%) and 7.4% had 7 or 8 species. Thus, variations in species proportions were important among the stands, varying from small to large presence.

Overview of ZELIG-CFS
ZELIG-CFS, a gap model that belongs to the JABOWA [21] type model, simulates individual-tree growth and mortality and regeneration on an annual time step. Its core structure is based on the integration of species-specific ecological properties ( Table 5) and representation of local environmental growth conditions, including monthly mean temperatures and precipitations, site fertility, and light interception at the individual-crown and canopy scales at different height strata on tree growth and seedling establishment. The basic approach in modelling tree growth consists in computing species-specific potential diameter increment (calculated from Dmax in Table 5), which is then modified by light available growing factor and species-specific shade tolerance class (Toler in Table 5), site growing degree-days parabolic response functions (estimated from temperature values and DDmin and DDmax in Table 5), soil moisture (estimated from precipitation values), and a soil fertility factor (estimated from site-specific potential biomass production). Regeneration establishment is modelled based on potential number of seedlings that can germinate on a site modified by understory light conditions adjusted by species-specific shade tolerance class, growing degree-days and soil moisture and fertility response classes. Mage is the maximum expected age for a species (year), Dmax is maximum diameter at breast height (cm) that a species can reach in its area of distribution, Hmax is the maximum height (cm) that a species can reach in its area of distribution, DDmin and DDmax are the minimum and maximum growing degree-days needed for a species to grow, Toler is species-specific shade tolerance class (1 = shade tolerant; 5 = shade intolerant), Drought is the species-specific drought tolerance class (1 = very drought-intolerant; 5 = very drought tolerant), and Nutri represents a soil fertility response class (1 = nutrient stress intolerant; 3 = nutrient stress tolerant).
ZELIG-CFS originates from the ZELIG model [39][40][41], but was modified to improve the representation of some processes. The major modification was made to the modelling of the available light growing factor (ALGF) and its effect on species-specific dbh growth rate. In ZELIG, individual-tree leaf area is assumed to be distributed uniformly with crown depth. In ZELIG-CFS, a sigmoidal cumulative leaf area distribution function is computed to more realistically represent the pattern of changes in light conditions with crown depth that integrates species-specific shade tolerance classes, which influences how crown recession rate occurs and, ultimately, dbh growth rate. Additionally, as crown recession takes place, change in crown width is adjusted based on species-specific crown shape functions: (1) a parabolic function for white spruce, northern white cedar, trembling aspen, jack pine and balsam poplar, (2) an ellipsoid function for American larch and white birch, (3) a conical function for balsam fir and (4) a cylindrical function for black spruce. This approach is consistent with findings reported by [42][43][44][45]. Other modifications were made to the mortality and regeneration model components. Mortality is handled by computing species-specific survival probability functions from dbh and dbh growth rate. For regeneration, a stocking factor representing space occupancy was added to the model. ZELIG-CFS was used to successfully predict the development of North American mixed forest types, including hardwood and softwood species [37,46] and subtropical dry forests [47].

Initialization Files
The development of each sample plot of the three data sets was simulated. For each sample plot, basic site information, including plot area and shape, latitude, longitude, altitude, soil field capacity and wilting points, soil fertility factor and monthly mean temperature and precipitation along with their standard deviations were stored in initialization files. Standard deviations for climatic data are used to randomly select monthly mean temperature and precipitation variations over the course of simulations. Site-specific monthly mean temperature and precipitation values were obtained from BIOSIM, an application that is used interpolate temperature and precipitation among regional weather stations [48].
Individual-tree data provided for each sample plot in the initialization files included tree number, species code and dbh. The first tree dbh measurements for each sample plot of the data sets were used as initialization data to initiate the simulations. They were considered as the state of the stands at simulated year 0. For each plot, the periodic observed data obtained from remeasurements over time were relative to simulated year 0, which facilitated the matching of observed and predicted basal areas and stand densities at specific simulated years.
Initial regeneration data provided as input in the initialization files included number of seedlings per square meter and stocking. The number of seedlings is, in fact, a potential number of seedlings per species that may germinate every year during the simulations. This number is then adjusted in every simulation year by a random value that represents the natural variability, understory light conditions, soil moisture and fertility and air temperature [37]. Stocking represents the proportion of the area that seedlings of a species may occupy within a forest ecosystem.
As regeneration data were not available in the data set, a procedure was developed to derive realistic estimates for both potential number of seedlings and stocking by assuming that regeneration was related to species-specific density in the tree layer at the individualplot scale. A literature review on regeneration studies conducted in the boreal forest of northern Ontario and Quebec allowed us to identify several studies providing results on both seedling abundance and stocking [49][50][51][52][53][54][55][56][57]. As no regeneration study was found for northern white cedar in northern Ontario or Quebec, values documented by [58] in the Ottawa valley, near the southern limit of the boreal forest were used.
For each sample plot and species within each data set, the potential number of seedlings (PNS) provided in the initialization files to run the simulations was estimated as: where MNS is an estimate of the maximum number of seedlings per m 2 that can be found in northern Ontario, as determined from the literature, STD the species-specific tree density in a sample plot and MAD the maximum species-specific tree density observed in the data set.
For each species present in the sample plots, stocking values (STOCK) were estimated as: where MINS is the minimum stocking value that can be found in northern Ontario, MIND the minimum species-specific tree density observed in the data set, and MAXS the maximum stocking value determined from the literature. For the initialization files, balsam fir had the highest estimated average potential number of seedlings per m 2 , followed by trembling aspen, black spruce and northern white cedar ( Table 6). The potential numbers of seedlings were much lower for jack pine, white spruce, white birch, American larch and balsam poplar. The large numbers of seedlings per m 2 for balsam fir and black spruce were also associated with relatively high average stockings ( Table 6). The highest average stocking was observed for white birch, but was not associated with a relatively high potential number of seedlings. Trembling aspen was among the species with high potential number of seedlings, but had a relatively low stocking. Relatively lower stocking values were observed for jack pine, white spruce, northern white cedar, American larch and balsam poplar. Relatively large standard deviations indicated much variation in the potential number of seedlings and stocking among the sample plots for each species. As the presence of northern white cedar was very marginal, simulation results for this species were not included here. Table 6. Average potential number of seedlings per m 2 and stocking values for each species provided as input for the sample plots in the initialization files.

Species
Number

Data Analysis
Three statistics were computed to compare model predictions against observed plot attributes, including basal area and stand density over time. Both basal area and stand density were computed for each plot and by species and year and then averaged for each year. The relative model bias was computed to compare predicted and observed basal areas for each sample plot as: This statistic was discussed in several studies that evaluated the performance of models [23,27,59,60]. It was averaged among plots to calculate an agreement indicator between simulated and observation data. High agreement between simulations and observations is associated with a low bias value. Negative values indicate that a model underestimates, while positive values indicate that it overestimates.
Model efficiency (ME) was computed to evaluate the capacity of a model to reproduce observations on a relative scale with values below 0 indicating a poor performance and 1 perfect performance [59]. It is computed as: where y i is observation i,ŷ i prediction i, and y the average of observations.
Observed and predicted diameter distributions were compared by computing a similarity coefficient (PS) [25] for every observation year to evaluate the extent to which observed and predicted dbh distributions differed: where x i and y i are the predicted and observed numbers of trees in dbh class i, respectively. This coefficient varies between 0 and 1 and indicates the proportion of values common to predicted and observed data sets [61].

Basal Area
Agreement between mean predicted and observed basal areas was high for black spruce in the two ecoregions for nearly every year and the standard deviations between observations and predictions overlapped considerably (Figure 2a,b). For Ecoregion 3E, percentages of differences between observations and predictions smaller than 2 and 5 m 2 ha −1 in absolute values were 54 and 78%, respectively ( Figure 3a). The maximum amplitudes of the majority of the differences were ±10 m 2 ha −1 over years and the number of differences smaller or higher than −10 or 10 m 2 ha −1 were marginal, representing only 6% of the differences. The same statistics for differences smaller than 2 and 5 m 2 ha −1 in absolute values for Ecoregion 3W were 73 and 90%, with a pronounced proportion of negative differences (Figure 3b). Only 2.9% of the differences between observations predictions were smaller or higher than 10 m 2 ha −1 over years and were observed in the past 30 years. Fluctuations were important for both ecoregions, but most of them occurred in the same years for predictions and observations (Figure 3a,b). Mean basal area increased slightly over time in Ecoregion 3E. For Ecoregion 3W, there was a relatively stable pattern of change over time, despite fluctuations.
The pattern of change over years in mean basal area for balsam fir was relatively stable for both ecoregions, but fluctuations were more pronounced than those for black spruce and did not always occur in the same years, particularly in Ecoregion 3W (Figure 2c,d). For Ecoregion 3W, agreement between predictions and observations was reasonably high from years 0 to 15, and for years 19, 20, 30, 41, 41 and 50, with large overlaps between the standard deviations of observations and predictions. At the plot scale and by year, 82 and 94% of the differences between predictions and observations in Ecoregion 3W were less than 2 and 5 m 2 ha −1 , respectively, and were within the same amplitude of differences over years ( Figure 3d). Only 6% of the differences between observations and predictions were higher than 5 m 2 ha −1 . A close examination of the sample plot data in observation years 25 and 26 indicated that the large differences between observations and predictions were caused by extreme values in nine sample plots. The pronounced increases in basal area that was observed from years 25 to 26 and from years 44 to 46, which increased the differences between observations and predictions, were probably caused by measurement or record errors in the historical data. Differences between mean predictions and observations for Ecoregion 3E were generally less pronounced than in Ecoregion 3W, but comparable percentages of differences among sample plots and years were observed: 83 and 95% were less than 2 and 5 m 2 ha −1 in absolute values, respectively (Figure 3c). The range of variation was the same over years, except for some large differences at years 10, 25, 26 and 47, which represented 2.5% of the total number of differences. Substantial overlaps in the standard deviations of observations and predictions and the occurrences of fluctuations were evident at the same time for most years. Agreement between mean observations and predictions was high for jack pine in Ecoregion 3E, which was characterized by a slight pattern of decrease over time (Figure 2e). Differences between predictions and observations smaller than 2 and 5 m 2 ha −1 represented 77 and 87% of the differences at the plot scale and by year. The amplitude of variation was relatively homogeneous over years, except for some large negative differences from year 10 to 60 (Figure 3e), representing 13% of the total number of differences. For Ecoregion 3W, mean observed and predicted basal areas for jack pine more than doubled in 48 years, despite fluctuations (Figure 2f). Differences between mean predictions and observations were generally higher than those in Ecoregion 3E (Figure 3f), but their standard deviations overlapped considerably for most years: 50 and 63% of the differences between observations and predictions were less than 2 and 5 m 2 ha −1 in absolute values, respectively, over years.
Mean predicted and observed basal areas for white birch in Ecoregion 3E fluctuated considerably and generally did not occur the same years (Figure 2g). Agreement was high between mean observations and predictions for the first 25 years. Then, differences generally increased until year 50, but were all not pronounced: 80 and 91% of the differences between observations and predictions were smaller than 2 and 5 m 2 ha −1 in absolute values, respectively (Figure 3g). The greatest differences between observations and predictions were observed from years 27 to 51 (8%). Mean observed and predicted basal areas in Ecoregion 3W were characterized by a pattern of decrease in the first 40 years, followed by a sharp increase (Figure 2h). Agreement was high between predictions and observations for the first 15 years. Then, relatively large negative differences were observed between years 15 and 43, with no strong overlap between the standard deviations of predictions and observations. When considered at the plot scale and by year, these relatively large differences can be attributed to sharp decreases in observations for some years only (Figure 3h), considerably decreasing the mean values. Despite these discrepancies, 86 and 96% of the differences between predicted and observed basal areas were less than 2 and 5 m 2 ha −1 over years, respectively. Compared to the other species, white spruce occupied small proportions of the stands, as indicated by the relatively small basal area values (Figure 4a,b). Mean predictions and observations were similar until year 10 and their standard deviations overlapped considerably for Ecoregion 3E (Figure 4a). However, compared to other species, there was a higher agreement between observations and predictions: 72% of the differences between predictions and observationsover years were smaller than 1 m 2 ha −1 and only 2 differences were greater than 2 m 2 ha −1 (Figure 5a). Even though differences between mean predictions and observations increased over years (Figure 4a), they remained relatively small (Figure 5a). Agreement was high between mean predicted and observed basal areas in Ecoregion 3W for all years (Figure 4b).
The standard deviations for both mean predictions and observations overlapped considerably and fluctuations occurred in the same years. Deviations between observations and predictions were similar to those in Ecoregion 3E: 95% of the differences between predictions and observations were smaller than 1 m 2 ha −1 in absolute values (Figure 5b).  Mean predictions and observations for balsam poplar in Ecoregion 3E fluctuated, with no obvious pattern of increase or decrease over time (Figure 4c). Despite apparent large discrepancies, the percentages of differences smaller than 2 and 5 m 2 ha −1 between predictions and observations were 62 and 81%, respectively, over years (Figure 5c). Very large differences between mean predictions and observations at years 11, 21, 30, 35 and 45 also resulted from sudden decreases in observed basal area relative to previous years in some sample plots. However, the number of large differences between observations and predictions represented only 12% of these values. Mean predictions and observations for balsam poplar in Ecoregion 3W were very close and their standard deviations overlapped considerably ( Figure 4d) and 76% of the differences between predictions and observations over years were less than 1 m 2 ha −1 (Figure 5d).
Results for trembling aspen were not as consistent as those for the other species. Large differences between mean predictions and observations were observed for nearly every year in Ecoregion 3E and standard deviations did not much overlap (Figure 4e). Large positive differences between predictions and observations were observed nearly every year (Figure 5e). For Ecoregion 3W, mean predicted basal area decreased over years and mean observed basal area fluctuated considerably (Figure 4f). However, agreement was relatively high for years 29, 30 and 42. In years 5, 10, 15, 25 and 32, the standard deviations for both predictions and observations overlapped. Despite large differences between mean predictions and observations for several years, 52 and 66% of the differences between observations and predictions were less than 2 and 5 m 2 ha −1 , respectively (Figure 5f). The large differences for years 26, 35 to 41 and 43 were caused by high mortality rates predicted by ZELIG-CFS for some sample plots.
American larch occupied a small proportion of the forests compared to the other species (Figure 4g,h). For the two ecoregions, agreement was high between mean predictions and observations and their standard deviations largely overlapped. Over years, 57 and 81% of the differences between observations and predictions were less than 1 m 2 ha −1 for ecoregions 3E and #w, respectively (Figure 5g,h).
Bias values indicated that, on average, ZELIG-CFS slightly underestimated basal area for black spruce in ecoregions 3E and 3W, but slightly overestimated basal area for jack pine in Ecoregion 3E (Table 7). In Ecoregion 3E, high positive bias values were observed for balsam fir, white birch, balsam poplar and American larch, but relatively high negative bias for white spruce and trembling aspen. Balsam fir and balsam poplar had high average bias values in Ecoregion 3E, but lower average bias values in Ecoregion 3W. The reverse was observed for jack pine between Ecoregions. In absolute values, lower bias values were observed for white birch, trembling aspen and American larch in Ecoregion 3W than in 3E. While a negative bias was observed for white spruce in Ecoregion 3E, a positive bias was observed in Ecoregion 3W. In many cases, high bias values were caused by extremes that were not necessarily associated with many substantial differences between predictions and observations. For instance, 25% of the bias values at the combined plot scale and by year for balsam fir in Ecoregion 3E were less than 10% and most of the high bias values occurred in simulation results past 20 years. Additionally, high bias values were observed in many sample plots that did not contain many trees of particular species. For instance, predicted basal area for balsam fir in a sample plot was 0.011 m 2 ha −1 , but observed basal area was 0.056 m 2 ha −1 , which resulted in a bias of 80% in absolute value. The predicted basal area represented 15 small trees with average dbh of 2.5 cm over 184 trees in a 809 m 2 sample plot. For the other five species in the sample plot, average dbh was 9 cm, with 42 trees larger than 15 cm in dbh. This pattern was observed in several sample plots for species with high average bias values (data now shown). Model efficiency was fairly high for black spruce and white spruce in both ecoregions and for white birch and balsam poplar in Ecoregion 3W (Table 7). It was greater than 0.50 for American larch in both ecoregions, and trembling aspen in Ecoregion 3W, between 0.40 and 0.50 for balsam fir in both ecoregions and jack pine in Ecoregion 3E, and very low for trembling aspen in Ecoregion 3E. Only two negative values were observed: balsam poplar in Ecoregion 3E and jack pine in Ecoregion 3W. Low or negative values of model efficiency observed for some species can be associated with the occurrence of substantial differences between observations and predictions, but they did not affect the general patterns of similarity of changes in basal area over time for most species (Figures 2 and 4).

Stand Density
Agreement for dbh class distributions for black spruce in Ecoregion 3E was relatively high between mean observed and predicted stand densities for several years and dbh classes (Figure 6a). Differences between mean observations and predictions were generally more pronounced in older years and large dbh classes. However, for most dbh classes and years, the standard deviations between observations and predictions largely overlapped. For the 45 cm dbh class, the absence of standard deviations was the result of lack of variation in the observed and predicted numbers of trees. These results are consistent with the similarity coefficients computed at the plot scale, which indicated a gradual decrease over year, but with relatively high coefficients in the first 40 years (Figure 7). Differences between mean observed and predicted stand densities in Ecoregion 3W were larger than those in Ecoregion 3E for the 5 cm dbh class, but agreement was relatively high for the 15, 25 and 35 cm dbh classes for most years (Figure 6a). For the 45 and 55 cm dbh classes, observed and predicted mean stand densities at years 45 and 55 were the same, which explains the absence of standard deviations. The similarity coefficients were generally smaller than those in Ecoregion 3E, particularly in the first few years, and were also characterized by a pattern of decrease over time (Figure 7). When analyzed by combining plots, years and dbh classes, results indicated that 47.3 and 61.2% of the differences between observed and predicted stand densities were less than 100 trees ha −1 for ecoregions 3E and 3W, respectively (Tables 8 and 9). Differences greater than 100 trees ha −1 may appear considerable, but, in reality, they are less important when adjusted to sample plot area That is, a difference of 100 trees ha −1 means a difference of 4 trees in a 400 m 2 sample plot. Differences greater than 100 trees ha −1 were observed for sample plots combined with years and dbh classes with high densities. Table 8. Percentages of plots in different classes (20 trees ha −1 ) of differences between observed and predicted stand densities (trees ha −1 ) in Ecoregion 3E, northern Ontario, Canada.  Table 9. Percentages of plots in different classes (20 trees ha −1 ) of differences between observed and predicted stand densities (trees ha −1 ) in Ecoregion 3W in northern Ontario, Canada.  For balsam fir in Ecoregion 3E, relatively good agreements were observed between mean observed and predicted densities for some years in the 15, 25 and 35 cm dbh classes (Figure 6b). Low agreement was evident for all years for the 5 cm dbh class. For year 50 for the 35 cm dbh class, ZELIG-CFS predicted the presence of 7 trees ha −1 , but none were present in the observations. However, 7 years later (year 57), trees were present in this dbh class, suggesting an error in the observation data set. No differences were observed between 22% of the observed and predicted stand densities and 37% were greater than 200 trees ha −1 ( Table 8). Percent differences decreased gradually from the lowest (1-20 trees ha −1 ) to the highest (181-200 trees per plot) density classes. Except for 11 years, the similarity coefficients were higher than 0.50 (Figure 7). Differences between observed and predicted mean stand densities were more pronounced in Ecoregion 3W than 3E, particularly for the 5 and 15 cm dbh classes (Figure 6b), which resulted in lower similarity coefficients (Figure 7). ZELIG-CFS predicted high densities for both dbh classes, but observations indicated that few trees were present. It is possible that errors in the observation data set may be the cause of the disparity. Despite the large differences between observations and predictions, agreement was high in many instances. The absence of differences between observations and predictions were evident for 33.8% of the combinations of sample plots, years and dbh classes, 33.7% of the differences varied between 1 and 200 trees ha −1 and 32.5% of the differences were greater than 200 trees ha −1 (Table 8). Larges differences between mean observed and predicted stand densities were observed for jack pine in Ecoregion 3E for most years in each dbh class (Figure 6c). However, the percentages of differences between observed and predicted stand densities suggest that the large differences were caused by only a relatively small number of combinations of sample plots, years and dbh classes in the results ( Table 8). As much as 32.7% of the differences between observed and predicted stand densities were less than 1 tree ha −1 , 56.8% of them were between 1 and 200 trees ha −1 and only 10.6% of them were greater than 200 trees ha −1 ( Table 8). The relatively high percentage of the absence of differences between observations and predictions may explain the occurrence of relatively high similarity coefficients (Figure 7). Except for the 5 cm dbh class, differences between observed and predicted stand densities were smaller for Ecoregion 3W than for 3E (Figure 6c). Substantial overlaps between the standard deviations of observations and predictions suggested that agreement was high in several instances. As much as 27.5% of the differences between observations and predictions were negligible, 47.2% were between 1 and 200 stems ha −1 , and 25.3% were greater than 200 stems ha −1 (Table 9). Similarity coefficients were greater than 0.50 for 13 years over a total of 40 years (Figure 7).

Species Classes of Ranges of Differences between Observations and Predictions in Stand
The patterns of mean dbh class distribution for white birch in Ecoregion 3E indicated high agreement between observations and predictions for some years and dbh classes (Figure 8a). Substantial differences were observed for the older years in the 5, 15 and 35 cm dbh classes and only for year 57 in the 25 cm dbh class. Large overlaps in the standard deviations between observations and predictions indicated agreement was high in several cases: differences between observations and predictions did not differ for 29.3% of the combinations of sample plots, years and dbh classes and only 16.2% of them were greater than 200 trees ha −1 (Table 8). Similarity coefficients were relatively low for most years (Figure 7). Results for Ecoregion 3W differed considerably from those for Ecoregion 3E (Figure 8a). Most of the large differences between mean observed and predicted stand densities were observed in the older years for most dbh classes, except for the 25 cm dbh class. As suggested by the large overlaps between the standard deviations of observations and predictions, agreement was high for many combinations of sample plots, years and dbh classes: no difference between observed and predicted stand densities were observed in 37.8% of the cases and only 12.9% had differences greater than 200 trees ha −1 ( Table 9). Most of the similarity coefficients were smaller than 0.50 (Figure 7). For white spruce in Ecoregion 3E, agreement between observed and predicted mean stand densities were high for most dbh classes and years (Figure 8b). There were no differences between observed and predicted stand densities for 36.7% of the combinations of sample plots, years and dbh class distributions and only 1.7% of them were greater than 200 trees ha −1 ( Table 8). The similarity coefficients were all greater than 0.60 ( Figure 7). Except for the 5 cm dbh class, agreement between observed and predicted stand densities was high for Ecoregion 3W (Figure 8b). In general, differences between observed and predicted stand densities were smaller than those in Ecoregion 3E (Figure 8b): 58.6% of the differences between observed and predicted stand densities were negligible and only 0.2% of them were higher than 200 trees ha −1 . The similarity coefficients were all greater than 0.60.
Differences between observed and predicted mean stand densities were large for balsam poplar in Ecoregion 3E for most of the dbh classes and years (Figure 8c). Relatively high agreement was observed only for year 10 in the 5, 15 and 25 cm dbh classes. However, despite the pronounced differences in the mean densities, the results at the sample plot scale and by dbh class and year indicated that 31.9% of observed stand densities coincided with the predicted densities and 14.5% of them were higher than 200 trees ha −1 ( Table 8). The similarity coefficients were low beyond year 20 (Figure 7). The comparison of mean stand densities for Ecoregion 3W indicated better agreement than that for Ecoregion 3E between observations and predictions ( Figure 8c). This can be seen in the percent differences between observed and predicted stand densities (Table 9). Only 0.3 % of the differences between observed and predicted stand densities were higher than 200 trees ha −1 and as much as 56.6% of them were negligible. The similarity coefficients were all greater than 0.60, except for years 31 and 33 (Figure 7).
In Ecoregion 3E, observed and predicted mean stand densities for trembling aspen differed considerably in the 5, 15, 25 and 55 cm dbh classes for nearly all years (Figure 9a). Relatively high agreements between observations and predictions were observed for years 10 and 20 in the 5 cm class and for most years in the 35 and the 45 cm dbh classes. Similar to the other species, there were large overlaps between the standard deviations for most combinations of sample plots, years and dbh classes that can be associated with several high agreements between observations and predictions. As much of 28.2% of the differences between observed and predicted stand densities did not differ and 54.9% of the percent differences were between 1 and 200 trees ha −1 ( Table 8). Most of the similarity coefficients were smaller than 0.50 ( Figure 7). Higher agreements between observed and predicted mean stand densities were observed for Ecoregion 3W than that for 3E for most years and dbh classes, except for year 30 in the 15, 25 and 45 cm dbh classes, year 40 in the 5, 15, 35 and 45 cm dbh classes (Figure 9a). Despite the large departures, agreement was high between observed and predicted stand densities in many cases (Table 9). There were no differences between observed and predicted stand densities for 36.9% of the combinations of sample plots, years and dbh classes and only 12.4 of them were greater than 200 trees ha −1 . Similarity coefficients were higher than 0.40 until year 35, except for years 29 and 33, then gradually declined, except for year 46 (Figure 7). For American larch, agreement was high between observed and predicted mean stand densities in Ecoregion 3E, except for years 20 to 57 in the 5 cm dbh class and years 10, 20, 40 and 50 in the 25 cm dbh class (Figure 9b). As much as 29.8% of the differences between observed and predicted stand densities were negligible and only 7.9% of them were greater than 200 trees ha −1 (Table 8). Most similarity coefficients were greater than 0.40 for most years (Figure 5d). American larch was less present in Ecoregion 3W (Figure 9b). Agreement was high between mean observations and predictions only in the 25 cm dbh class. As much of 66.7% of the differences between observations and predictions were negligible and 11.1% of them were higher than 21 trees ha −1 (Table 9). These results can be associated with relatively high similarity coefficients (Figure 7).

Discussion
As far as could be evaluated from the literature covering recent decades, support is high for the existence of forest gap dynamics in boreal forests [10][11][12][13][14][15]. The present study is, however, one of the rare efforts, such as the study by [62], that have evaluated the performance of a gap model using extensive long-term historical data for forest types varying substantially in species composition. This issue was raised by [26] who highlighted the need for using extensive data sets for evaluating the performance of gap models to avoid bias toward some species or ages. The factors that must be considered in the comparison of the predictions of gap models with observed stand data can be classified in three groups [29]: (1) degree to which the environmental site variables used in the simulations are accurate and representative of the field conditions; (2) extent to which observed data, which may originate from inventory surveys, represent the forests under examination; and (3) extent to which potential successional changes and management activities are taken into account. As previously mentioned, ZELIG-CFS is a semi-mechanistic model that includes basic environmental site variables for the simulations, including monthly mean temperatures and precipitations and basic soil texture, in the initialization files. For climatic and soil texture data, the estimates were fairly accurate and representative of field conditions at the sample plot scale. The application BioSIM [48], which was used to estimate monthly mean temperatures and precipitations, is based on interpolations from regional weather stations by taking each sample plot's latitudes, longitudes and altitudes into account to the coordinates of each sample plot. Soil texture data were obtained from soil maps and basic description of soil characteristics for each sample plot. The historical long-term sample plot data set used in this study originated from a databank developed and maintained by the NDMNRF Growth and Yield Program and the measurements were conducted by ministry staff or forest companies mentioned above, which provides a strong confidence in the validity of the measurements. Of course, measurement errors were possibly made at the time of measurements or introduced in inventory reports, but it is unlikely that they were sufficient to invalidate a large proportion of the data set. The third factor was less relevant for the present study, as very few sample plots were subject to management activities, such as partial cuts. However, since remeasurement data spanned relatively long periods of time, the successional pathways were likely well captured.
The patterns of change in predicted mean basal area over relatively long periods were consistent with observations for most species in both ecoregions, except for trembling aspen in Ecoregion 3E. The simultaneous occurrences of fluctuations in observations and predictions for several years suggest that the patterns of predicted mortality and regeneration were generally well captured by ZELIG-CFS. Relatively small fluctuations from years to years are normal in the course of stand dynamics when observations and predictions from several sample plots are assembled and analyzed. However, large fluctuations were observed in mean basal area in several cases within very short periods, such as those in Ecoregion 3W from years 25 to 27 for trembling aspen and from years 43 to 50 for balsam fir and, in Ecoregion 3W, from years 42 to 44 for black spruce. These rapid changes might have been caused by some errors in the records.
ZELIG-CFS generally overpredicted stand density, particularly in the small dbh classes, but, similarly to the results by [63,64], predicted that the differences between observed and predicted stand densities tended to decrease with increasing dbh class for most species. This pattern can be explained by high regeneration rates predicted by ZELIG-CFS for several sample plots, which increased predicted stand density in the smallest dbh classes. Then, differences between observed and predicted densities were reduced over time, as young trees died in the years following their establishment. Thus, our results suggest that ZELIG-CFS consistently predicted small tree mortality in the years after their establishment. However, overprediction of regeneration did not occur consistently in all sample plots. As indicated above, large overlaps in the standard deviations were evident between observed and predicted mean densities, which could be associated with many cases of agreement between observed and predicted stand densities and high similarity coefficients. Nevertheless, these patterns do not indicate that ZELIG-CFS consistently failed to realistically predict regeneration establishment. Large differences between observations and predictions do not necessarily indicate that a model does not adequately represent mechanisms or failed to make realistic predictions [29]. A good example of this statement is the modelling of regeneration [37]. ZELIG-CFS requires as input the potential numbers of seedlings that might germinate every year and then randomly selects the number of seedlings that will appear by considering the site conditions, including understory light conditions and soil moisture. This modelling procedure is biologically consistent and is an optimal solution given the lack of experimental data. Many surveys have been conducted to estimate regeneration establishment in forest ecosystems, but they are just a picture at one point in time and, as such, do not capture potential annual variations and early seedling survival. This may partially explain some of the differences between observed and predicted stand densities in small dbh classes. Despite consistent predictions of regeneration, these results suggest that there is a need to gain a better understanding on the different stages of regeneration, including seed germination, seedling growth and mortality and transition from seedlings to trees to improve the modelling of regeneration. However, the fact that differences between observed and predicted stand densities diminished over time for some species suggests that ZELIG-CFS consistently and realistically predicted tree mortality, which is a complex process involving interactive effects of several factors [65][66][67]. The prediction of mortality may be more difficult for some species. For instance, mortality in paper birch stands may be important throughout all its life, unless suppressed trees are released early [68].
Black spruce is a dominant species in both ecoregions. It was present in 94 and 97% of the sample plots in ecoregions 3E and 3W, respectively, and was characterized by high basal areas and stand densities for nearly all years and dbh classes. These results can be related to the ecological properties of black spruce. Black spruce grows well in a wide range of ecological conditions, including organic and mineral soils, cold and high temperatures and low and high precipitation [69,70]. Additionally, as a shade-tolerant species, it grows well in pure and mixed forests and its seedlings may grow in as little as one-tenth of understory sunlight [71]. Black spruce occupies a special niche as a late succession species where most succession pathways converge in the boreal forest of Ontario [6].
Balsam fir was less present in terms of basal area and stand density than black spruce in both ecoregions, but was found in 65 and 48% of sample plots in ecoregions 3E and 3W, respectively. It has growing requirements similar to those of black spruce with respect to temperature and precipitation [71,72] and is a late successional species that may maintain itself in mixed forests with black spruce, white spruce and white birch for long periods in the absence of disturbances, particularly fires [71]. Additionally, its seedlings can grow well under poor understory light conditions for long periods, but respond vigorously to increases in understory light conditions following gap openings [71]. Its development in mixed stands and longevity was consistently represented by ZELIG-CFS.
Basal area for jack pine was similar in both ecoregions for approximately the first 20 years of the simulations. While there was a pattern of slight increase in basal area in Ecoregion 3W, it decreased in Ecoregion 3E. Jack pine occurred more often in Ecoregion 3W than in 3E. The decline in Ecoregion 3E may be explained by the fact that jack pine is a transient species that is eventually replaced by late successional species, such as black spruce or balsam fir, that better tolerate understory shade conditions [69]. Additionally, as jack pine is a shade-intolerant species, its seedlings do not grow well in low-light conditions [70,71]. Two possible reasons for its more frequent occurrence in Ecoregion 3W than 3E are: first, more stands were in a state of successional transition in Ecoregion 3E. Second, soil conditions were more favorable for jack pine in Ecoregion 3W. The best soil conditions for jack pine are well drained or dry loamy acidic sands, such as those Ecoregion 3W [38,69]. Large areas in Ecoregion 3E have calcareous fine texture soils, including silts or clay [38], which are less favorable for jack pine.
White birch was present in 58 and 50% of the sample plots in ecoregions 3E and 3W, respectively. It is a shade-intolerant species that can grow as a minor component in openings of mixed late successional forests [18,69]. White birch presented patterns similar to those of jack pine with respect to basal area in both ecoregions: a decline in Ecoregion 3E, but an increase in 3W, despite fluctuations. Again, more stands were in a state of transition in Ecoregion 3E than 3W, reducing white birch occurrence. Another explanation is that it may grow as codominant with trembling aspen, black spruce, white spruce, jack pine and balsam fir, but rarely is a late successional species [69]. The presence of white birch in mixed old forests may be associated with its robust capacity to produce and propagate seedlings on different soil conditions [68,72]. The fact that its seedlings can develop under 50% of full sunlight probably contributes to its higher densities in Ecoregion 3W [65], contrary to jack pine. The effect of differences in soil properties and temperature and precipitation conditions between both ecoregions may not explain differences in basal area development in both ecoregions. White birch usually grows on well-drained sandy loams and is adapted to wide variations in temperature and precipitation [70].
White spruce was present in 4 and 33% of the sample plots in ecoregions 3E and 3W, respectively, but stand density was generally higher in Ecoregion 3E than 3W. Its ecological properties do not explain the differences of its development between ecoregions. White spruce can grow in large variations of climatic conditions and soil properties, including acid or calcareous loams and clays and wet or dry sites [69]. As suggested by the basal area and stand density values relative to the other species, it is a minor component in most stands surveyed in the present study. White spruce can grow in successional transitions with trembling aspen and paper birch that may eventually become pure stands, but may also reach a late successional status with other boreal species, such as black spruce, which explains why it was found in so many mixed stands for both observations and predictions [69]. In particular, white spruce may benefit from the presence of trembling aspen [73].
Despite fluctuations, basal area and stand density for balsam poplar were generally higher in Ecoregion 3E than 3W and this species was present in 37 and 22% of the sample plots in the respective ecoregions. Balsam poplar is a low shade-tolerant species that grows in a wide range of climatic conditions and, although it grows best on moist sites, it can tolerate dry sites [69]. It is an early successional species that colonizes disturbed wet sites that will generally lead to mixed forests [69]. The fine textured soils with silts or clay likely explains its greater presence and development of basal area in Ecoregion 3E. These soils retain more water than the generally well drained soils in Ecoregion 3W.
Trembling aspen was the species with the lowest agreement between observed and predicted basal areas and stand densities in both ecoregions. It is a pioneer species that grows well in various cold boreal climatic conditions and can colonize disturbed sites with poor soils, including bare soils, but requires moist conditions to grow well [69,70,74]. It is widely distributed in North America [70], but was found in only 26 and 57% of the sample plots in this study, which may indicate that sample plots were not evenly disturbed, considering its capacity to invade disturbed sites. It is also possible that, as a shadeintolerant species, its seedlings normally cannot grow below its own canopy [74], unless there are gaps caused by disturbances, such as windthrow, as observed in aspen-white spruce forests [74,75]. An explanation for the large differences between observations and predictions may be that ZELIG-CFS did not well capture the complexity of the mortality process for this species in the region. It seems that mortality does not occur on a regular basis for trembling aspen, as there may be episodes of high mortality rates in young stands caused by competition for light from shade-tolerant dominant species that will replace them in the course of succession, drought conditions or intensive browsing [70,76].
American larch was a minor species in the data set. It was present in 6 and 5% of the sample plots in ecoregions 3E and 3W, respectively, and its basal area and stand density were smaller than those for most species. The simulation results agree with what has been observed for this species. In the boreal forest of Ontario, American larch may grow alone or with other species, including black spruce, balsam fir, white spruce, white birch or trembling aspen, but remains in small proportions in mixed stands [69]. The reasons that may explain why it was found in only a few sample plots include its specific site requirements and its shade intolerance to the extent that its seedlings cannot grow under its own canopy. Its best growth conditions occur on moist to poorly drained soils, including acid swamps, bogs, or muskegs, but it can also be found in moist areas along streams or lakes [65].

Conclusions
The evaluation of the performance of the succession model ZELIG-CFS for various boreal mixed forest types in northern Ontario indicated that it generally predicted the development of most species well when compared with historical data sets and also captured the effects of differences in species-specific ecological properties. Some important discrepancies between observations and predictions could be partly explained by potential errors in the data sets. However, as with any model, the representations of some mechanisms could be improved. The major one involves the dynamics of the transitions from seedling to sapling to tree stages. For most species, understanding about these transitions remains sparse as forest ecosystem studies focus on young or mature forests. Even with current knowledge or information, the modelling of some processes, such as regeneration, competitive effects or temperature and precipitation variation, will always include some degree of uncertainty. Most available models provide deterministic results. Including uncertainty algorithms to provide uncertainty estimates in model predictions might be helpful.