Skip to main content

Historical and future carbon stocks in forests of northern Ontario, Canada



Forests in the Far North of Ontario (FNO), Canada, are likely the least studied in North America, and quantifying their current and future carbon (C) stocks is the first step in assessing their potential role in climate change mitigation. Although the FNO forests are unmanaged, the latter task is made more important by growing interest in developing the region’s natural resources, primarily for timber harvesting. In this study, we used a combination of field and remotely sensed observations with a land surface model to estimate forest C stocks in the FNO forests and to project their future dynamics. The specific objective was to simulate historical C stocks for 1901–2014 and future C stocks for 2015–2100 for five shared socioeconomic pathway (SSP) scenarios selected as high priority scenarios for the 6th Assessment Report on Climate Change.


Carbon stocks in live vegetation in the FNO forests remained relatively stable between 1901 and 2014 while soil organic carbon (SOC) stocks steadily declined, losing about 16% of their initial value. At the end of the historical simulation (in 2014), the stocks were estimated at 19.8, 46.4, and 66.2 tCha−1 in live vegetation, SOC, and total ecosystem pools, respectively. Projections for 2015–2100 indicated effectively no substantial change in SOC stocks, while live vegetation C stocks increased, accelerating their growth in the second half of the twenty-first century. These results were consistent among all simulated SSP scenarios. Consequently, increase in total forest ecosystem C stocks by 2100 ranged from 16.7 to 20.7% of their value in 2015. Simulations with and without wildfires showed the strong effect of fire on forest C stock dynamics during 2015–2100: inclusion of wildfires reduced the live vegetation increase by half while increasing the SOC pool due to higher turnover of vegetation C to SOC.


Forest ecosystem C stock estimates at the end of historical simulation period were at the lower end but within the range of values reported in the literature for northern boreal forests. These estimates may be treated as conservatively low since the area included in the estimates is poorly studied and some of the forests may be on peat deposits rather than mineral soils. Future C stocks were projected to increase in all simulated SSP scenarios, especially in the second half of the twenty-first century. Thus, during the projected period forest ecosystems of the FNO are likely to act as a C sink. In light of growing interest in developing natural resources in the FNO, collecting more data on the status and dynamics of its forests is needed to verify the above-presented estimates and design management activities that would maintain their projected C sink status.


Forests occupy over 40 million km2 globally, accounting for about 30% of total land area, and store about 860 billion tonnes of C in live biomass and dead organic matter pools [51]. The amount of C stored in forests continues to increase, partly offsetting increasing atmospheric CO2; forests are the main component of the terrestrial C sink, estimated as of 2018 at 3.5 billion tonnes of C annually or about 30% of total annual emissions [21]. Not surprisingly, forests feature prominently in discussions about climate change and its mitigation [62]. Of particular importance is estimating forest’s future C stocks and whether they will continue to act as a C sink given that climate change is projected to substantially affect forest condition and growth. On the one hand, rising air temperature and atmospheric CO2 concentrations can stimulate forest productivity, especially in areas where soil moisture, nitrogen (N), and phosphorous are not a limiting factor [17]. On the other hand, positive effects on productivity may be counteracted by heat-related stress in plants, worsening drought conditions, and changes in disturbance regime [17]. Therefore, quantifying the effects of climate change on forests is needed to develop climate change mitigation strategies.

Canadian forests occupy about 3.56 million km2, of which 3.09 million km2 are in the boreal zone and the remainder in the temperate zone [53]. Of the total area, 2.3 million km2 are covered by managed forests, as defined for the purposes of annual reporting on greenhouse gas (GHG) emissions and removals, while another 1.26 million km2 of forested area in northern latitudes are unmanaged [37]. Managed forests are much better studied in terms of their C stocks, due to international reporting obligations as well as historical interest in collecting data and developing tools for timber supply needs, with various components serving as inputs to C balance estimation [35]. The current state of C stocks is much less known for unmanaged forests, and despite recent efforts in data collection (e.g., National Forest Inventory plot grid), major gaps and uncertainties remain unaddressed [35].

In Ontario, the unmanaged portion of the forest is in the northern part of the province, usually referred to as the Far North of Ontario (FNO), that lies roughly above the 51st parallel (51°N; see Fig. 1). This area is likely the least studied part of northern forests in North America compared with, for example, Alaska and Quebec, where considerable effort has been focused on the assessment of forest resources and C stocks [41, 45]. Meanwhile, the effects of climate change in the FNO are projected to be more pronounced relative to other parts of Canada; for example, increases in the average winter temperature by the end of twenty-first century are projected to double those in southern Ontario [14].

Fig. 1
figure 1

Boundaries of the Far North of Ontario (FNO), Canada, with area in green classified as non-wetland forest. Delineated area in bottom left corner of the FNO is the Whitefeather Forest Management Unit

The first published attempt at estimating current and predicting future C stocks in forests of the FNO was undertaken by [26]. In that study, historical and future C stocks were estimated from 1900 to 2100 for two climate change scenarios (representative concentration pathways RCP4.5 and RCP8.5). The current study builds on [26] by using newly available data sources, correcting winter temperature bias, and expanding future projections to cover the broader spectrum of climate change scenarios. More specifically, the objective was to simulate changes in vegetation and SOC stocks for historical and five shared socioeconomic pathway (SSP) scenarios from the Sixth Assessment Report on Climate Change by the IPCC [42]. The latter projections would provide forest managers and policy makers with benchmark estimates of forest C stocks in the unmanaged forests of Ontario and indicate whether, in the twenty-first century, these forests will act as a C sink or source under possible climate change scenarios.


Carbon stocks were simulated for 120,634 km2 of area classified in the land cover map of the FNO as non-wetland coniferous, deciduous, or mixedwood forest (Fig. 1). Hereafter this area is referred to as the FNO forested area or the FNO forests; see "Methods" for study area description. Simulations were performed for historical data for 1901–2014 and five SSP scenarios for 2015–2100 (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) Mean annual temperature, total annual precipitation, and atmospheric CO2 concentration for the simulation periods are shown in Fig. 2, while summary statistics are presented in Table 1.

Fig. 2
figure 2

Mean annual temperature (a), total annual precipitation (b) in the Far North of Ontario, and surface CO2 concentration in Northern Hemisphere (c) from historical (Climate Research Union observations, [42]) for 1901–2014 and shared socioeconomic pathway scenarios for 2015–2100 (CanESM5 projections, [42])

Table 1 Summary statistics on annual air temperature and annual total precipitation during 2015–2100 for the Far North of Ontario for five shared socioeconomic pathway (SSP) scenarios

Simulated historical (1901–2014) C stocks in live vegetation, soil, and total ecosystem are presented in Fig. 3, along with projected (2015–2100) C stocks for SSP2-4.5 scenarios with and without wildfires; maps of spatial distribution of C stocks density across forested area of the FNO for 2015 and 2100 simulated for SSP2-4.5 scenario wildfires are shown in Fig. 4. The latter scenario was selected for presentation of the results because it is a “middle-of-the-road” one among the five scenarios used in this study, with the mean 2015–2100 annual air temperature and atmospheric CO2 concentration of 3.62 °C and 522 ppm, respectively (see Table 1 and Fig. 2). As seen in Fig. 3, between 1901 and 2014, C stocks in live vegetation in the FNO forests remained relatively stable, ranging from a minimum of 18.3 tCha−1 in 1918 to a maximum of 20.2 tCha−1 in 1901. Soil organic carbon stocks, on the other hand, steadily declined from 54.9 tCha−1 in 1901 to 46.3 tCha−1 in 2014; here and throughout the text SOC includes organic C in the upper 1 m of soil and dead organic matter (DOM) pools such as litter, downed dead wood, and standing dead trees. Consequently, total forest ecosystem C stocks also declined from 75.2 tCha−1 in 1901 to the minimum of 65.9 tCha−1 in 2002 and then rebounded slightly to 66.1 tCha−1 in 2014 (Fig. 3).

Fig. 3
figure 3

Historical (1901–2014) and projected (2015–2100) for SSP2-4.5 scenario carbon density in live vegetation (green lines), soil organic carbon (SOC; brown lines), and total ecosystem (blue lines) in the forests of the Far North of Ontario (FNO). Solid and dashed lines correspond to projections with and without fire, respectively; vertical dashed line indicates the end of historical simulation

Fig. 4
figure 4

Map of total ecosystem carbon density for the forested area of the Far North of Ontario (FNO) for a 2015 and b 2100 for SSP2-4.5 scenario with fires

As expected, model runs showed strong effects of wildfires on the projected C stocks. In the absence of wildfires, live vegetation C stocks doubled from 19.9 tCha−1 in 2014 to 40.7 tCha−1 in 2100; simulation of wildfire activity reduced this increase almost by half, to 31.7 tCha−1 in 2100 (Fig. 3). Simulating wildfires had the opposite effect on SOC, which effectively showed no change from 46.3 tCha−1 in 2014 to 46.9 tCha−1 in 2100 with wildfires, and decreasing to 43.7 tCha−1 in 2100 in the absence of wildfires. Summation of live vegetation C and SOC resulted in forest ecosystem stocks increasing from 66.2 tCha−1 in 2014 to 84.3 tCha−1 and 78.6 tCha−1 in 2100 without and with wildfires, respectively (Fig. 3).

The patterns presented in Fig. 3 were consistent among all five climate change scenarios; changes in total forest ecosystem C stocks for the five scenarios, with and without simulation of wildfires, are presented in Fig. 5. In all scenarios, total forest ecosystem C stocks dropped slightly in 2020–2025 losing the maximum of 0.14 in SSP1-2.6 scenario. Subsequently, the stocks steadily increased until 2100, reaching values between 77.3 tCha−1 and 79.9 tCha−1 (Fig. 5). Thus, by the end of the century, total forest ecosystem C stocks in the FNO are projected to increase by 16.7–20.7% relative to 2015. The results of simulation with and without wildfires for all climate change scenarios are shown in the Additional file 1.

Fig. 5
figure 5

Projected (2015–2100) total forest ecosystem carbon stock density in forests of the Far North of Ontario (FNO) for five climate change scenarios: SSP1-1.9 (blue); SSP1-2.6 (brown); SSP2-4.5 (green); SSP3-7.0 (yellow); SSP5-8.5 (red). Solid and dashed lines correspond to projections with and without fire, respectively

Figure 6 shows forest ecosystem, vegetation, and soil C fluxes (i.e., annual changes in C stocks) simulated for the historical period 1991–2005 in this study and by three Earth System Models (ESMs): CanESM2, HadGEM2, and IPSL. Results are shown for 1991–2005 because it is the overlapping period between InTEC and the three ESM simulations forced by observed changes in atmospheric composition and land cover. Figure 7 shows C fluxes for 2015–2100 projected in this study for SSP2-4.5 scenario and the averages of five-year moving minima and maxima over six scenarios projected by the three ESMs and calculated as follows. For a given year and for each model-scenario combination, a minimum and maximum are estimated over a five-year period centering on the given year, and the average minimum and maximum are then calculated for the six scenarios. The results for scenario SSP2-4.5 from this study are shown in Fig. 7 since this scenario is most similar to RCP4.5 scenarios used in the ESMs. The fluxes projected by the ESMs for all six scenarios are shown in Fig. S1.

Fig. 6
figure 6

Historical (1901–2005) carbon stock changes in a total ecosystem, b soil and c live vegetation in forests of the Far North of Ontario (FNO) simulated by InTEC (this study) and three Earth System Models (CanESM2, HadGEM2, and IPSL)

Fig. 7
figure 7

Projected (2015–2100) carbon stock changes in a total ecosystem, b soil and c live vegetation in forests of the Far North of Ontario (FNO) simulated by InTEC (this study). Highlighted in grey is the area between five-year moving average of maxima and minima from six scenarios simulated by three Earth System Models (CanESM2, HadGEM2, and IPSL); see text for more detail


Live vegetation C stocks estimated for 2014 were within the range of values reported for northern boreal forests, albeit at the lower end. Globally, live vegetation C in the boreal forest biome is estimated at 46.7 tCha−1 as of 2007 [51]. However, it is not uncommon to see much lower estimates, especially for the northern boreal forest. For example, Gower et al. (2001) reported values of 12.8, 16.1, 99.9, and 15.3 tCha−1 in black spruce (Picea mariana (Mill.) BSP)-dominated stands in Alaska, USA, aged 51, 55, 130, and 200 years old, respectively. Other estimates of mean live tree biomass for the same region include 26.0 tCha−1 for three black spruce stands ranging in age from 160 to 200 years [56] and 25.5 tCha−1 for four black spruce-dominated stands with an average age of 95 years [13]. Schvidenko and Nilsson [60] estimated live vegetation C stocks at 23.9 tCha−1 for forest tundra in the northern and sparse taiga bio-climatic zone in Russia while [45] used two different methods to arrive at estimates of 18.8 and 25.4 tCha−1 for total aboveground vegetation in the northern and southern boreal zones of Quebec, respectively.

Similarly, our estimate of SOC is at the low end of values reported in the literature. For example, Schulze et al. [58] estimated the average SOC content of a chronosequence (28–383 years) of six Scots pine (Pinus sylvestris L.) stands in Central Siberia, Russia, at 37.8 tCha−1, of which 18.2 tCha−1 were in mineral soil and 19.6 tCha−1 in DOM pools (forest floor, downed dead wood, and standing dead trees). In a similar study in Central Siberia, the average mineral soil and DOM C stocks in four chronosequences of Scots pine stands (total of 22 stands) was estimated at 15.6 tCha−1 and 31.8 tCha−1, respectively [73]. In a 68-year-old Norway spruce (Picea abies L.H.Karst) stand growing in the northern boreal zone in the European part of Russia, C stocks in mineral soil and dead organic surface layers were estimated at 43.5 tCha−1 and 7.8 tCha−1, respectively [32]. In the latter three studies [32, 58, 73], C in mineral soil was measured to 0.5 m depth. However, based on the analysis of nearly 17,000 forest soil measurements, De Vos et al. [16] concluded that the top 0.5 m accounts for 77% of the total C content in 1 m of mineral soil. Thus, if we expand the estimates from [32, 58, 73] by applying a correction factor of 1.3 estimated in [16], the resulting C content in the DOM pool is about the same as that estimated in this study.

It is worth noting that, in 2014, C stocks in live vegetation in the FNO forest accounted for 34.4% of total forest ecosystem C stocks. This number is similar to the estimated fraction of live vegetation C stocks in several other studies, e.g., 30.1% in [34] or 33.7% in [9]. In other words, our simulated estimate of the ratio of live vegetation C to SOC is within the range of values reported in the literature.

Live vegetation C stocks simulated in this study remained relatively stable during 1900–2015. These results are consistent with the conclusions reached by [8]; the authors of the latter study analysed more than 20,000 sample plots in five states across the United States and did not find substantial increases in forest growth rates during the twentieth century. In a large-scale study that included permanent sample plots in five provinces, no significant trends were observed in the rate of biomass change in eastern Canada (Ontario and Quebec) from mid-seventies of the twentieth century to the early years of this century [39]. In an unrelated large-scale study based on the analysis of tree-ring width data covering about 60 years starting in 1950, significant changes were not found in productivity of black spruce and jack pine (Pinus banksiana Lamb.) forests growing in the Boreal Shield West region of Canada, which includes the FNO [25].

During the same historical period, SOC stocks steadily declined, diminishing nearly 16% by 2014 relative to 2001. This decline is likely caused by two factors. The first one is an increase in soil respiration and rate of organic matter decomposition caused by rising mean air temperature that, during 1901–2014, on average increased by 1.5 °C in the FNO [26]. For example, a heating experiment in a red spruce (Picea rubens Sarg.)-dominated forest in Maine, USA, showed an increase in litter decay as a result of warming, and consequently a significant loss of litter mass and C content [57]. In another study in mixed hardwood forest in Massachusetts, USA, soil C losses were recorded as a result of increased temperature, primarily due to microbial respiration associated with soil organic matter decay that was responsible for over three quarters of total annual CO2 efflux from soil [43]. The second factor is the above-discussed stability of live vegetation stocks that is a source of C transfers to SOC pool. The combination of increasing efflux and constant influx resulted in the steady loss of SOC simulated in our study.

Trends in historical C stocks and their fluxes simulated in this study are similar to those simulated for the FNO using three ESMs (Fig. 6 and Additional file 1). Estimated as a slope of a simple linear regression, the trend in historical vegetation C stocks (0.0025) and C fluxes (0.0009) are within the ranges for trends in stocks (0.0011–0.0159) and fluxes (0.0009–0.0029) of those simulated using the ESMs. The same is true for the trends in soil C fluxes (0.0007 in this study vs. 0.0–0.0019 for ESMs) and in ecosystem C fluxes (0.0015 in this study vs. 0.0010–0.0035 for ESMs). The only historical trends projected in this study that fall out of those simulated by the ESMs are those for soil C stocks (− 0.0685 in this study vs. 0.0–0.0014 for ESMs) that in turn dictate the negative trend for ecosystem C stocks (-0.0662 vs. 0.0012–0.0173 for ESMs). However, vegetation C stocks are better studied and quantified than those for soil, and vegetation C stocks and fluxes simulated in this study closely agree with both the empirical observations [25, 39] and simulations using the ESMs. This finding, combined with the above-discussed interplay between vegetation and soil C stocks serves as a verification of historical C stocks and fluxes simulated in this study for the FNO.

Simulations of C stocks during 2015–2100 were performed for five climate change scenarios based on the so-called shared socioeconomic pathways (SSPs) combined with various mitigation efforts [24]. The SSPs describe how the world might evolve based on socioeconomic factors such as population, economic growth, education, urbanisation, and the rate of technological development [55]. The five scenarios used in this study were selected as “high-priority scenarios” for the Sixth Assessment Report on Climate Change by the IPCC and included: three scenarios—SSP1-2.6, SSP2-4.5, and SSP5-8.5—approximately corresponding to the previous generation scenarios RCP 2.6, RCP 4.5, and RCP 8.5, respectively; “medium–high” scenario SSP3-7.0; and SSP1-1.9, which most closely reflects a 1.5 °C target under the Paris Agreement [42]. For more detail on the development of SSPs and for the so-called narratives (i.e., textual description of how the future might unfold in terms of broad societal trends) of SSPs selected for this study the reader is referred to [55].

Vegetation C stocks in the twenty-first century increased in all climate change scenarios. Similar trends have been simulated in other studies. For example, simulations of total land C (vegetation plus SOC) in British Columbia, Canada, projected an increase in C stocks from 2005 to 2050 by 12% and 16% (compared to 2005) for the lowest (RCP2.6) and highest (RCP8.5) scenarios, respectively [3]. Changes modelled in our study for the same period are lower, likely because SOC stocks are projected to remain relatively stable throughout the twenty-first century. In another North American study, changes in terrestrial C stocks in Alaska, USA, were simulated by two general circulation models for three SRES scenarios (A1B, A2, and B1) [23]. The authors of the latter study projected an increase in upland vegetation C stocks in all six scenarios ranging from 12 to 24%. Modelling studies for boreal forests in Europe also indicated positive trends in C stocks and productivity. Increases in net primary productivity of the European boreal zone are projected by the end of the twenty-first century for A1B and B1 scenarios generated by three general circulation models [54]. Similar results were produced by [33] who used a variety of climate change scenarios to evaluate sources of uncertainty in projections of the primary productivity in boreal forests of Finland and concluded that productivity was likely to increase by the end of the century.

The changes in C stocks projected in this study are consistent with those simulated using the three ESMs. For example, forest ecosystem C stocks are projected to increase by 18.7% during 2015–2100 in the SSP2-4.5 scenario; respective changes simulated by the ESMs range from − 7.3% (for the CanESM2-Climate scenario) to 28.8% (for the HadGEM2-CO2 scenario). A comparison of future linear trends in C stocks (as was done above for the historical period) would not be informative due to the non-linear nature of C stock dynamics (Fig. 5). The C stock fluxes, however, also indicate good agreement between the results of this study and the CMIP5 project: as seen from Fig. 7, for both SSP2-4.5 scenarios (with and without wildfires) the C stock fluxes for all three pools (ecosystem, soil, and vegetation) are within the range of minima and maxima simulated using the three ESMs.

Results from several ESMs indicate that rising atmospheric CO2 will promote C uptake by both plant biomass and soil organic matter [65, 67]. These results seem to be based on the concept of increased vegetation growth resulting in higher input to and therefore an increase in soil C. However, the meta-analysis of 53 experiments undertaken by [67] indicated that atmospheric CO2 enrichment stimulates both the input and the turnover of C in soil, and therefore does not significantly affect soil C content. In other words, increased CO2 stimulates plant growth and C input to soil; however, the response is counteracted by the microbial response, reducing the net effect on soil C stocks [67]. This finding is echoed by [65] who state “there is both empirical and theoretical evidence that increases in soil inputs, especially under elevated CO2, may have little effect on soil organic carbon stocks.” Results of another meta-analysis of 49 field experiments, located across North America, Europe, and Asia, suggest that the effect of warming depends on the amount of initial soil C stock: effects are negligible for sites with small initial C stock but occur beyond the stock values of 200–400 tCha−1 and are considerable in soil with stocks higher than 600 tCha−1 [15]. The area simulated in our study belongs to the former category of low initial SOC stocks; hence, stocks are projected to remain relatively stable throughout the simulation period. As stated by [15], in ecosystems with low initial SOC losses resulting from increased decomposition rate due to warming, changes are relatively minor and can be offset by concurrent increases in vegetation growth and soil C stabilization; however, in areas with large initial SOC, increased decomposition occurs faster than C accumulation from enhanced vegetation growth and therefore leads to substantial C losses to the atmosphere [15].

One of the predicted outcomes of the future C stock dynamics is the change in allocation of C between live vegetation and SOC pools (Fig. 3). By the end of the historical simulation (i.e., in 2014), live vegetation accounted for 30% of total ecosystem C. That fraction increased by 2100 to about 40%, ranging from 39.3% to 41.8% among scenarios. This shift in C allocation with the increasing temperature and atmospheric concentration of CO2 has been observed in empirical studies. For example, in a CO2-enrichment experiment in Massachusetts, USA, a deciduous forest stand increased the share of C stored in live vegetation while losing some of the SOC stock [43]. The shift in relative distribution of C between forest floor and aboveground vegetation with increasing soil temperature has been reported in a comparative study of C allocation in nine black spruce stands in boreal forests of Canada and the United States [69]; similar change in relative distribution of C with increasing air temperature was observed in Scots pine forests in Europe [69].

Previous simulations demonstrated the effects of strong winter warmth biases inherent in ESMs [26]. The latter biases may lengthen the growing season, which in turn would result in an increase in vegetation C stocks and consequently SOC. For the current simulations, the winter biases were corrected for all climate variables at a monthly time scale with the linear scaling bias correction method using observed and simulated data from overlapping years (2015–2018). Comparison between the closest scenarios in this study and [26] (SSP5-8.5 and RCP8.5 + historical fire scenarios, respectively) showed that without winter bias correction, by the end of the twenty-first century, live vegetation, SOC stocks, and total ecosystem C stocks were overestimated by 15.5%, 3.0%, and 7.9%, respectively. This finding is consistent with the effect of bias correction on changes in C pools reported by [2] and emphasizes the importance of bias correction for more realistic prediction of C stocks.

Total C stocks in forests of the FNO in 2014 were estimated at 798.6 million tonnes, of which live vegetation and DOM contained 239.2 and 559.4 million tonnes, respectively. For comparison, the total amount of C in 29.4 million hectares of Ontario’s forests managed for timber production has been estimated at 4,719.2 million tonnes, for the average density of 160.5 tCha−1 [11]. The lower C density in FNO forests may be caused by differences in species composition, climatic conditions, and the lack of active wildfire suppression measures applied in the managed forest. Also, since we excluded areas classified as treed wetland, our results may underestimate C storage in FNO forests. It is also possible that some of the forests included in our study are on peat deposits rather than mineral soils [26]. Overall, this uncertainty underscores the need for collecting data on the condition and dynamics of the FNO forest; primarily forest inventory information, which is lacking for almost the entire FNO, but also sample plots that would allow assessment of changes in forest composition and characteristics over time.

Simulations with and without wildfires illustrated the effects of wildfires on C stocks in the FNO forests (Fig. 3). In the absence of wildfires, live vegetation C stocks in SSP2-4.5 scenario increased twofold by the end of twenty-first century; the increase was much smaller when wildfires were simulated, reflecting losses of C to combustion emissions and transfers to SOC pool [35]. The latter transfers result in wildfires having the opposite effect on SOC stocks that are higher in the scenario with wildfires (Fig. 3). Simulations without wildfires also corroborated the previously discussed effect of increasing air temperature on soil respiration. In the absence of wildfires and other disturbances, soil C stocks are primarily controlled by C inputs from the vegetation pool and C losses to the atmosphere due to soil respiration. As vegetation C stocks increase, C inputs into the soil pool increase [12] so to maintain stable soil C stocks, these inputs must be offset by increases in soil respiration (Fig. 3). Thus, simulations without fire indicate that soil respiration increases with air temperature which is consistent with empirical observations [43, 57].

Simulation of wildfires also sheds light on their relative contribution versus that of climate change to changes in C stocks of the FNO forests (Figs. 3 and 5). The difference in projected live vegetation C stocks between scenarios with and without wildfires may be attributable to changes in climatic variables and CO2 fertilisation. This attribution is supported by the pattern of C stock dynamics during the first half of the twenty-first century. Without wildfires, live vegetation (Fig. 3) and total ecosystem (Fig. 5) C stocks start increasing immediately after 2015 because wildfire simulation stops at the end of historical period, leading to increases in forest age. However, in the scenario with wildfires forest age is unlikely to increase, and changes in C stocks become noticeable much later since it takes longer time for the effects of gradual changes in climate and atmospheric CO2 concentration to become pronounced. It is worth noting that, by the end of simulation period, the rate of increase in ecosystem C stocks in simulations with wildfires fire matches or exceeds that in scenarios without wildfire (Fig. 5). This effect is similar to previously reported simulations in which increases in the average annual fire-burned area do not reduce forest ecosystem C stocks when the effect of CO2 fertilization is accounted for [4, 26]. While the latter result may seem counterintuitive, the proposed explanation lies in wildfire-caused changes to forest age structure: fires replace older slow growing forest stands with young ones that can take advantage of the increasing atmospheric CO2 concentration. The resulting accelerated biomass growth combined with slow decomposition of fire-transferred matter to SOC pools is sufficient to compensate for C losses due to direct combustion emissions.

One of the limitations of this study is the assumption of equal probability of burning for all pixels classified as forest (with the exception of forests younger than 11 years). Meanwhile, indications are that individual wildfire events may exhibit selectivity for species composition and forest age. For example, Bernier et al. [6] in a Canada-wide study of wildfires during 2002–2011 concluded that probability of burning was significantly higher in forests older than 30 years as well as in forests dominated by conifers (vs. broadleaved stands). Another limitation is the lack of dynamic representation of species composition; albeit not unique in the modelling literature (e.g., see [3]), the assumption of fixed spatial distribution of forest types does not reflect possible effects of climate change on this distribution. For example, studies in Alaska (another large area of unmanaged northern boreal forests in North America) suggested that changes in the wildlife frequency and intensity have the potential to increase the extent of deciduous forest cover at the expense of the conifers [23, 30]. More research is needed to assess how total forest C stocks would be affected by these two potentially counterbalancing phenomena, i.e., higher probability of wildfires targeting older conifer stands that contain higher C stocks vs. reduced probability of wildfires resulting from the increasing fraction of deciduous forests.

Climate change is projected to increase the burned area in the FNO. For example, [20] estimated the annual burned area for a scenario corresponding to a threefold increase in the atmospheric concentration of CO2, with 1.5–2.0 increases projected in the western part of the FNO (containing most of the FNO forests) and 2.1 increase in the eastern part relative to levels registered during 1959–1997. Boulanger et al. [7] projected the annual burned area to increase less than twofold and fourfold in the western and eastern parts of the FNO, respectively; the latter increases were projected to occur during 2071–2100 in the A2 scenario [7]. These predictions are likely to be at the upper limit of possible changes in the fire-burned area since triple-CO2 and A2 scenarios are close to the most aggressive RCP8.5 and SSP3-7.0-SSP5-8.5 scenarios [42, 68]. Unfortunately, neither paper [7, 20] contained the models used in their predictions. In the precursor to this study [26], future fires were simulated using equations from [5]. However, this meant using equations with relatively low explanatory power (R2 < 0.5) and extrapolating them to the parts of the FNO for which no equations have been estimated. As noted in [5], this is one of the areas in eastern Canada where developing fire models driven by fuel moisture and temperature is more challenging due to the influence of large water bodies (Hudson and James bays). In the absence of reliable models predicting future fire regime in the FNO for various climate change scenarios, we decided on the conservative approach of simulating the historical rate of annual burned area.

Finally, future emissions from wildfires may also be affected by increases in proportion of combusted soil material. For example, a study in black spruce-dominated stands in Alaska indicated that the depth of burning in the ground layer of biomass (including moss, litter, organic soils, and peat) increased with the length of the wildfire season [66]. Similarly, studies conducted in boreal forests of the Northwest Territories (Canada) showed that C soil combustion depended on the forest stand’s landscape position and age [70, 71]. As with the above-discussed ambiguity of wildfires’ forest type selectivity and possible shifts to more deciduous-dominated forests, more studies are needed to assess whether higher growth rate in young stands due to CO2 fertilization can overcome C losses from the increased spatial extent and depth of burning of future wildfires. The latter is important as the province of Ontario, among many other jurisdictions, attempts to reduce its GHG emissions [44]. Currently, the FNO forests fall in the category of unmanaged and therefore, according to international rules, their C balance would not count towards GHG emission reduction targets. However, despite the current surplus in harvestable timber in the managed forests of Ontario [50], the interest in developing natural resources of the FNO, primarily timber harvesting, is growing [48]. Harvesting is the main forest management activity affecting forest C dynamics in North American forests by reducing their C stocks and emitting C to the atmosphere [17]. Although various management activities could increase forest uptake of atmospheric C and decrease emissions in the forest sector [37], it would be challenging to develop resource management in the FNO while maintaining the projected C sink status of its forests. However, regardless of possible management activities, the principal need is the above-mentioned data collection in the FNO forests that would allow verification of their current state and reduce uncertainty about their future dynamics.


We simulated historical (1901–2014) and future (2015–2100) C stocks in forests of the FNO. These simulations improved on previous results by using a better-defined land base and a larger number of more up-to-date climate change scenarios. Simulated live vegetation C stocks were relatively stable during the historical period, while SOC stocks steadily declined, losing about 16% of their initial amount by the end of the historical period. Consequently, total forest ecosystems C stocks were also reduced by about 12%. Carbon stock estimates at the end of historical simulation period were at the lower end, but within the range, of values reported in the literature for northern boreal forests. These estimates may be treated as conservatively low; the area included in the estimates is poorly studied and it is possible that some of the forests are on peat deposits rather than mineral soils. Future C stocks were projected for five SSP scenarios selected as “high-priority scenarios” for the Sixth Assessment Report on Climate Change by the IPCC. Soil organic carbon stocks were projected to remain relatively constant from 2015 to 2100; live vegetation C stocks, however, were projected to increase, especially in the second half of the twenty-first century. These results were consistent among all five simulated climate change scenarios. The results were improved by applying a correction to winter warmth bias; in the absence of this correction, vegetation, SOC, and total ecosystem C stocks would be overestimated by up to 15.5%, 3.0%, and 7.9%, respectively. By and large, during the projected period forest ecosystems of the FNO are likely to act as a C sink. In light of growing interest in developing natural resources in the FNO, collecting more data on the status and dynamics of its forests is needed to be able to verify the above-presented estimates and design management activities that would maintain their projected C sink status.


Methods and several input data sets repeat those used in the precursor to this study [26]. Therefore, here we briefly describe models and input data used in the study, with more attention to the novel elements. The latter include new landcover data, updated version of historical climate data, and new projections of future climate data corresponding to five SSPs (as opposed to only two scenarios of an older vintage used in [27]). Future climate data were corrected for winter biases the effects of which were evident in the previous simulations [27].

The FNO is an area of 439,756 km2 located approximately to the north of 51 °N and bound from the west and east by Manitoba and Hudson and James bays, respectively [47] (Fig. 1). The FNO overlaps two ecozones: Hudson Plains and Boreal Shield. Forested areas are located mostly in the Boreal Shield since flat terrain and poor drainage of the Hudson Plains has resulted in the largest contiguous wetlands in the world. Tree species composition in the FNO is typical of northern boreal forests of Canada, with black spruce being the dominant species, particularly on lowland sites, along with white spruce (Picea glauca (Moench) Voss), jack pine (Pinus banksiana Lamb.), trembling aspen (Populus tremuloides Michx.), tamarack (Larix laricina (Du Roi) K. Koch), and white birch (Betula papyrifera Marsh.) [26]. This study includes land cover cells classified as non-wetland coniferous (104,426 km2), deciduous (5,770 km2), and mixedwood (10,438 km2) forests totalling 120,634 km2 [27]; areas classified as treed wetland were not included. The long-term (1901–2015) average mean annual air temperature and total annual precipitation are − 1.3 °C and 640 mmyear−1, respectively. Forests in the study area are classified as unmanaged. In 2009, a fraction of the study area along the southwest boundary of the FNO was designated an addition to the forests managed for timber production (Fig. 1). A forest management plan for 6,270 km2 of forested area in the newly formed Whitefeather Forest Management Unit (FMU) was developed in 2012. However, to date no harvesting or other management activity has occurred in the Whitefeather FMU [49]; therefore, the entire forested area included in this study is considered unmanaged.

The forest C cycle was modelled using the Integrated Terrestrial Ecosystem Carbon Cycle (InTEC) model developed to simulate C balance in Canada’s forests by integrating effects of disturbance and non-disturbance factors such as climate, CO2 concentration, and N deposition [12]. Complete theories and formulations of the InTEC model are documented by [10, 12, 31]. In the model, carbon balance is simulated using 13 C pools, four of them describing vegetation (foliage, stem, fine root, and coarse root), and nine for SOC (surface structural litter, soil structural litter, woody litter, surface metabolic litter, soil metabolic litter, surface microbial, soil microbial, slow soil organic matter C pool, and passive soil organic matter C pool); simulated SOC pools account for C content to one m depth. The InTEC model consists of three components: a canopy-level photosynthesis module for simulating net primary productivity (NPP), a module for simulating soil C and N dynamics, and a hydrological module for simulating soil moisture and temperature.

The soil C and N dynamics module is based on the CENTURY model [52] modified to account for multiple soil C pools; temporally and spatially varying N depositions; the effects of drainage, soil temperature, and moisture on decomposition rate; and climatic and C pool size effects on N fixation. The hydrological module, parameterized based on fractions of sand, clay, silt, and organic matter and vegetation properties, simulates soil water content and temperatures of three soil layers at monthly time steps needed to quantify decomposition rates of soil C pools and soil water stress effects on photosynthesis. The photosynthesis module was developed from a canopy-level Farquhar’s leaf biochemical model using a temporal and spatial scaling scheme [12, 19]. The module quantifies the integrated effects of changes in stand age, climate, and CO2 and N deposition since the preindustrial period on the interannual variability of NPP to progressively calculate annual NPP from an initial NPP value. The NPP value in a reference year (we used 2004 because that was the year all spatial data sets were available), simulated at daily time steps using the Boreal Ecosystem Productivity Simulator (BEPS) [38], was the benchmark used to tune the initial NPP value. For each pixel, the initial NPP value was repeatedly adjusted until the difference between NPP simulated by InTEC and the benchmark output from BEPS in the reference year (i.e., 2004) was less than 1%. Collection of the field data used to validate the BEPS reference NPP validations and development of NPP-age curves are described in detail in [27].

To initialize various C pools, InTEC assumes them to be in a steady state before 1901 for stands disturbed after 1901 or in the year before the most recent disturbance for stands undisturbed after 1901. The initialization was run until the C pools reached a steady state in which the absolute value of the ecosystem net C balance became smaller than 2% of the initial NPP, using as the latter the NPP and stand age in 1901 along with mean climatic conditions in 1901–1910 to initialize C pools in vegetation and SOC. The sizes of the various C pools were estimated by solving a set of differential equations that consider the interaction among pools under steady state C dynamics [10].

All wildfires simulated by InTEC are assumed to be stand-replacing, i.e., causing complete stand mortality. For both historical and projected simulations, the constant number of burned cells is randomly allocated to all forest stands older than 11 years; the probably of getting burned by a wildfire does not depend on forest type and age. Once burned, a fraction of biomass and soil C is emitted to the atmosphere while the rest of biomass C is transferred to the soil pools. The amount of directly emitted C includes 100% of foliage, 100% of surface structural and surface metabolic litter, and 25% of stem wood; the latter fraction accounts for the losses of primary and secondary branches and part of the boles. Dead biomass C remaining after fire is transferred to woody litter and surface metabolic and structural litter pools, and its decomposition is assumed to start the year after fire disturbance. The model does not consider fire intensity and hence does not account for the effect of the latter on C losses. The forest disturbed by wildfires regenerates without a forest cover type change in the second year after a disturbance.

Inputs to the InTEC model included spatial data sets of climatic variables, N deposition, soil texture, drainage, digital elevation model, land cover, leaf area index (LAI), forest stand age, reference NPP, and fire disturbance; the spatial data sets covered the entire FNO forest area. Prior to model runs, all spatial data sets were interpolated to 500 m resolution, the highest spatial resolution of the remote sensing input data sets. Four new data sets were compiled for the current study (compared to [26]), namely: new land cover data, new version of historic climate data (1901–2014), future climate data (2015–2100), and future fire data (2015–2100). A brief description of input data sets is given below.

Land cover data were obtained from the Far North Land Cover, version 1.4 [46]. The initial 24 classes included in the data set were aggregated into six classes (i.e., water, wetland, treed wetland, tundra heath, coniferous forest, deciduous forest, and mixedwood forest). Only non-wetland coniferous, deciduous, and mixedwood forests were considered. The area classified as wetland and treed wetland was excluded because processes controlling C stock dynamics in these ecosystems differ from those in forests [22] and are not simulated by the InTEC model. Gridded monthly climate data, including mean air temperature, water vapour pressure, and precipitation at 500 m resolution for the period 1901–2014 were interpolated bilinearly from the 0.5° global data set produced by the Climate Research Unit, version CRU TS v. 4.03 [29]. The monthly incoming shortwave radiation data for 1901–2014 were obtained from 20th Century Reanalysis V3 data set provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA [61].

The future climate data were from the Canadian Earth System Model (CanESM5) prepared for the sixth phase of the Coupled Model Intercomparison Project (CMIP6) that is under the auspices of the World Climate Research Programme. Future climate data for five scenarios (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5) were downloaded from [72]. The CanESM5 data were interpolated bilinearly from the 2.8o × 2.8° original spatial resolution to 500 m. Previous simulations demonstrated the effects of strong winter warmth biases inherent in all Earth System Models [26]. For the current simulations, the winter biases are corrected using the so-called delta change and linear scaling bias correction methods [40]; albeit relatively simple, they have been shown to perform well at monthly scale [59]. In our study, the delta change and linear scaling bias correction methods were applied using the observed and simulated data for 2015–2018. For each variable, the overlapping period (2015–2018) provided four years or 48 monthly data points that were used to derive a correction factor. The corrections, based on the difference in the mean monthly value of observed and simulated temperature and on the ratio of the mean monthly value of observed divided by simulated climate data, were applied to monthly simulated data throughout the 2015–2100 period. The ratio was derived for each climate variable separately. The CO2 and total CO2 equivalent GHG concentrations for the selected SSP scenarios needed to estimate future N deposition were obtained from [42].

Measurements of N deposition for Canadian forests during 1983–1994 made by the Canadian Air and Precipitation Monitoring Network [18] were used to produce a N deposition map for Canada, and then extrapolate it temporally based on the relationship between deposition rate and historical national greenhouse gas emissions [10]. Drainage, soil water characteristics, and soil texture (fractions of sand, silt, clay, and organic matter) were compiled from the Soil Landscape of Canada database [1]. Missing soil data values were interpolated based on a nearest neighbour scheme. The 30-arcsec digital elevation model data used to calculate the wetness index and local slopes for simulating the horizontal redistribution of soil water was derived from the Canadian Digital Elevation Data at 1:250,000 scale by the Canadian Forest Service [36].

The LAI, NPP, and wildfire fire disturbance data used in the study are described in detail by [26]. The LAI was derived from the moderate resolution imaging spectroradiometer (MODIS) 8-day reflectance observations and validated using field measurements collected in 2010 in 32 sample plots established in unmanaged forest in the southwestern part of the FNO. Reference NPP for 2004 was simulated using the BEPS model. Two sets of field measurements were used (1) to validate the BEPS reference NPP simulations and (2) to develop NPP-age curves used in simulations by InTEC model. Finally, data on historical fires for 1959–1999 were obtained from the Large Fire Database [63]. These data were supplemented with remote sensing observations to verify and expand the historical records to 2004. The historical fire data were used to assign forest age to the stands regenerating in the areas burned during 1959–2004 and to estimate the average annual wildfire rate that was used in simulations for 2005–2100. The rationale for choosing the historical wildfire rate in simulations of future C stocks is presented in "Discussion" section.

Carbon stocks for the FNO forests were simulated using historical data for 1901–2014 and using five SSP scenarios for 2015–2100. For each SSP scenario, two simulations were performed: with and without wildfires. The results for simulated C pools are combined in live vegetation, SOC, and total ecosystem pools and presented in the form of C density (amount of C per unit area). For cross-model comparison we used the following three models within the Coupled Model Intercomparison Project (CMIP5) project for which soil and vegetation C stocks were available in fully coupled model simulations [64]: the Canadian Earth System Model version 2 (CanESM2), the Hadley Centre Global Environmental Model version 2 (HadGEM2-ES) and the Institut Pierre Simon Laplace Climate Model 5A (IPSL-CM5A-LR). We used three experiments simulated by these models [64]: the ‘historical’ experiment for 1901 − 2005 is driven by prescribed CO2, aerosol, solar forcing, climate and land use change forcing; the ‘esmFdbk2’ experiment for 2006 − 2100 (referred to as “Climate”) with constant pre-industrial CO2 for physiological effects on vegetation and RCP4.5 greenhouse gas concentrations for radiative calculations; and the ‘esmFixClim2’ experiment for 2006 − 2100 (referred to as “CO2”) with RCP4.5 greenhouse gas concentrations for physiological effects on vegetation and pre-industrial greenhouse gas concentrations for radiative calculations. Data simulated for the above-listed scenarios were acquired for a geographic area spanning 50 − 55°N and 82 − 95°W, encompassing the entire forested area of the FNO.

Availability of data and materials

Results of the analysis and references to the data sources are included in the Additional file 1 accompanying this paper.



Boreal Ecosystem Productivity Simulator




Coupled Model Intercomparison Project


Dead organic matter


Earth System Model


Forest management unit


Far North of Ontario


Greenhouse gas


Leaf area index




Net primary productivity


Soil organic carbon


Shared socioeconomic pathway


  1. [AAFC] Agriculture and Agri-Food Canada. Soil Landscapes of Canada Working Group, Soil Landscapes of Canada, version 3.2 (digital map and database at 1:1 million scale). 2010. Accessed 15 Dec 2020.

  2. Ahlström A, Schurgers G, Smith B. The large influence of climate model bias on terrestrial carbon cycle simulations. Environmental Research Letters. 2017;12(1):014004.

  3. Arora VK, Peng Y, Kurz WA, Fyfe JC, Hawkins B, Werner AT. 2016. Potential near‐future carbon uptake overcomes losses from a large insect outbreak in British Columbia, Canada. Geophysical Research Letters. 2016;43(6):2590–2598.

  4. Balshi MS, McGuire AD, Duffy P, Flannigan M, Kicklighter DW, Melillo J. Vulnerability of carbon storage in North American boreal forests to wildfires during the 21st century. Glob Change Biol. 2009;15(6):1491–510.

    Article  Google Scholar 

  5. Balshi MS, McGuire AD, Duffy P, Flannigan M, Walsh J, Melillo J. Assessing the response of area burned to changing climate in western boreal North America using a Multivariate Adaptive Regression Splines (MARS) approach. Glob Change Biol. 2009;15(3):578–600.

    Article  Google Scholar 

  6. Bernier PY, Gauthier S, Jean PO, Manka F, Boulanger Y, Beaudoin A, Guindon L. Mapping local effects of forest properties on fire risk across Canada. Forests. 2016;7(8):157.

    Article  Google Scholar 

  7. Boulanger Y, Gauthier S, Gray DR, Le Goff H, Lefort P, Morissette J. Fire regime zonation under current and future climate over eastern Canada. Ecol Appl. 2013;23(4):904–23.

    Article  Google Scholar 

  8. Caspersen JP, Pacala SW, Jenkins JC, Hurtt GC, Moorcroft PR, Birdsey RA. Contributions of land-use history to carbon accumulation in US forests. Science. 2000;290(5494):1148–51.

    Article  CAS  Google Scholar 

  9. Chen J, Colombo SJ, Ter-Mikaelian MT, Heath LS. Carbon budget of Ontario’s managed forests and harvested wood products, 2001–2100. For Ecol Manage. 2010;259(8):1385–98.

    Article  Google Scholar 

  10. Chen JM, Ju W, Cihlar J, Price D, Liu J, Chen W, Pan J, Black A, Barr A. Spatial distribution of carbon sources and sinks in Canada’s forests. Tellus B Chem Phys Meteorol. 2003;55(2):622–41.

    Google Scholar 

  11. Chen J, Ter-Mikaelian MT, Ng PQ, Colombo SJ. Ontario’s managed forests and harvested wood products contribute to greenhouse gas mitigation from 2020 to 2100. For Chron. 2018;43(3):269–82.

    Google Scholar 

  12. Chen W, Chen J, Cihlar J. An integrated terrestrial ecosystem carbon-budget model based on changes in disturbance, climate, and atmospheric chemistry. Ecol Model. 2000;135(1):55–79.

    Article  CAS  Google Scholar 

  13. Cleve KV, Oliver L, Schlentner R, Viereck LA, Dyrness CT. Productivity and nutrient cycling in taiga forest ecosystems. Can J For Res. 1983;13(5):747–66.

    Article  Google Scholar 

  14. Colombo SJ, Gray PA, Partington PJ, Pearson D. Beyond 450 parts per million: Climate change hazards in a 4°C warmer world and how Ontario can help avoid them. Ontario Centre for Climate Impacts and Adaptation Resources, Sudbury, ON/Ontario Ministry of Natural Resources and Forestry, Peterborough, ON. 2015. Accessed 15 Dec 2020.

  15. Crowther TW, Todd-Brown KE, Rowe CW, Wieder WR, Carey JC, Machmuller MB, Snoek BL, Fang S, Zhou G, Allison SD, Blair JM. Quantifying global soil carbon losses in response to warming. Nature. 2016;540(7631):104–8.

    Article  CAS  Google Scholar 

  16. De Vos B, Cools N, Ilvesniemi H, Vesterdal L, Vanguelova E, Carnicelli S. Benchmark values for forest soil carbon stocks in Europe: Results from a large scale forest soil survey. Geoderma. 2015;251:33–46.

    Article  CAS  Google Scholar 

  17. Domke G, Williams CA, Birdsey R, Coulston J, Finzi A, Gough C, Haight B, Hicke J, Janowiak M, de Jong B, Kurz WA. Forests. In Second State of the Carbon Cycle Report (SOCCR2): A Sustained Assessment Report. Cavallaro, N.; Shrestha, G.; Birdsey, R.; Mayes, MA; Najjar, RG; Reed, SC; Romero-Lankao, P.; Zhu, Z., eds. Washington, DC, USA: US Global Change Research Program. 2018:365–98.

  18. [ECCC] Environment and Climate Change Canada, Canadian Air and Precipitation Monitoring Network. 2017. Accessed 15 Dec 2020.

  19. Farquhar GD, von Caemmerer SV, Berry JA. A biochemical model of photosynthetic CO2 assimilation in leaves of C 3 species. Planta. 1980;149(1):78–90.

    Article  CAS  Google Scholar 

  20. Flannigan MD, Logan KA, Amiro BD, Skinner WR, Stocks BJ. Future area burned in Canada. Clim Change. 2005;72(1–2):1–6.

    Article  CAS  Google Scholar 

  21. Friedlingstein P, Jones M, O'sullivan M, Andrew R, Hauck J, Peters G, Peters W, Pongratz J, Sitch S, Le Quéré C, Bakker DCE. Global carbon budget 2019. Earth System Science Data. 2019;11(4):1783–838.

  22. Frolking S, Talbot J, Jones MC, Treat CC, Kauffman JB, Tuittila ES, Roulet N. Peatlands in the Earth’s 21st century climate system. Environ Rev. 2011;19:371–96.

    Article  CAS  Google Scholar 

  23. Genet H, He Y, Lyu Z, McGuire AD, Zhuang Q, Clein J, D’Amore D, Bennett A, Breen A, Biles F, Euskirchen ES. The role of driving factors in historical and projected carbon dynamics of upland ecosystems in Alaska. Ecol Appl. 2018;28(1):5–27.

    Article  Google Scholar 

  24. Gidden M, Riahi K, Smith S, Fujimori S, Luderer G, Kriegler E, van Vuuren DP, van den Berg M, Feng L, Klein D, Calvin K. Global emissions pathways under different socioeconomic scenarios for use in CMIP6: a dataset of harmonized emissions trajectories through the end of the century. Geoscientific Model Development Discussions. 2019;12(4):1443–75.

    Article  CAS  Google Scholar 

  25. Girardin MP, Bouriaud O, Hogg EH, Kurz W, Zimmermann NE, Metsaranta JM, de Jong R, Frank DC, Esper J, Büntgen U, Guo XJ. No growth stimulation of Canada’s boreal forest under half-century of combined warming and CO2 fertilization. Proc Natl Acad Sci. 2016;113(52):E8406–14.

    Article  CAS  Google Scholar 

  26. Gonsamo A, Chen JM, Colombo SJ, Ter-Mikaelian MT, Chen J. Global change induced biomass growth offsets carbon released via increased forest fire and respiration of the central Canadian boreal forest. J Geophys Res Biogeosci. 2017;122(5):1275–93.

    Article  Google Scholar 

  27. Gonsamo A, Ter-Mikaelian MT, Chen JM, Chen J. Does earlier and increased spring plant growth lead to reduced summer soil moisture and plant growth on landscapes typical of Tundra-Taiga interface? Remote Sensing. 2019;11(17):1989.

    Article  Google Scholar 

  28. Gower ST, Krankina O, Olson RJ, Apps M, Linder S, Wang C. Net primary production and carbon allocation patterns of boreal forest ecosystems. Ecol Appl. 2001;11(5):1395–411.

    Article  Google Scholar 

  29. Harris I, Osborn TJ, Jones P, Lister D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. Scientific Data. 2020;7(1):1–8.

    Article  Google Scholar 

  30. Johnstone JF, Hollingsworth TN, Chapin FS III, Mack MC. Changes in fire regime break the legacy lock on successional trajectories in Alaskan boreal forest. Glob Change Biol. 2010;16(4):1281–95.

    Article  Google Scholar 

  31. Ju W, Chen JM, Black TA, Barr AG, Liu J, Chen B. Modelling multi-year coupled carbon and water fluxes in a boreal aspen forest. Agric For Meteorol. 2006;140(1–4):136–51.

    Article  Google Scholar 

  32. Kalinina O, Chertov O, Dolgikh AV, Goryachkin SV, Lyuri DI, Vormstein S, Giani L. Self-restoration of post-agrogenic Albeluvisols: Soil development, carbon stocks and dynamics of carbon pools. Geoderma. 2013;207:221–33.

    Article  CAS  Google Scholar 

  33. Kalliokoski T, Mäkelä A, Fronzek S, Minunno F, Peltoniemi M. Decomposing sources of uncertainty in climate change projections of boreal forest primary production. Agric For Meteorol. 2018;262:192–205.

    Article  Google Scholar 

  34. Kauppi PE, Posch M, Hänninen P, Henttonen HM, Ihalainen A, Lappalainen E, Starr M, Tamminen P. Carbon reservoirs in peatlands and forests in the boreal regions of Finland. Accessed 15 Dec 2020.

  35. Kurz WA, Shaw CH, Boisvenue C, Stinson G, Metsaranta J, Leckie D, Dyk A, Smyth C, Neilson ET. Carbon in Canada’s boreal forest—a synthesis. Environ Rev. 2013;21(4):260–92.

    Article  CAS  Google Scholar 

  36. Lawrence K, Hutchinson M, McKenney D. Multi-scale digital elevation models for Canada. Frontline, Forestry Research Applications. 2008. Accessed 15 Dec 2020.

  37. Lemprière TC, Krcmar E, Rampley GJ, Beatch A, Smyth CE, Hafer M, Kurz WA. Cost of climate change mitigation in Canada’s forest sector. Can J For Res. 2017;47(5):604–14.

    Article  Google Scholar 

  38. Liu J, Chen JM, Cihlar J, Chen W. Net primary productivity mapped for Canada at 1-km resolution. Glob Ecol Biogeogr. 2002;11(2):115–29.

    Article  Google Scholar 

  39. Ma Z, Peng C, Zhu Q, Chen H, Yu G, Li W, Zhou X, Wang W, Zhang W. Regional drought-induced reduction in the biomass carbon sink of Canada’s boreal forests. Proc Natl Acad Sci. 2012;109(7):2423–7.

    Article  CAS  Google Scholar 

  40. Maraun D. Bias correcting climate change simulations-a critical review. Current Climate Change Reports. 2016;2(4):211–20.

    Article  Google Scholar 

  41. McGuire AD, Genet H, Lyu Z, Pastick N, Stackpoole S, Birdsey R, D’Amore D, He Y, Rupp TS, Striegl R, Wylie BK. Assessing historical and projected carbon balance of Alaska: A synthesis of results and policy/management implications. Ecol Appl. 2018;28(6):1396–412.

    Article  Google Scholar 

  42. Meinshausen M, Nicholls ZR, Lewis J, Gidden MJ, Vogel E, Freund M, Beyerle U, Gessner C, Nauels A, Bauer N, Canadell JG. The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500. Geoscientific Model Develop. 2020;13(8):3571–605.

    Article  CAS  Google Scholar 

  43. Melillo JM, Butler S, Johnson J, Mohan J, Steudler P, Lux H, Burrows E, Bowles F, Smith R, Scott L, Vario C. Soil warming, carbon–nitrogen interactions, and forest carbon budgets. Proc Natl Acad Sci. 2011;108(23):9508–12.

    Article  CAS  Google Scholar 

  44. [MECP] Ministry of the Environment, Conservation and Parks. Preserving and protecting our environment for future generations: A made-in-Ontario environment plan. 2018. Accessed 15 Dec 2020.

  45. Nelson R, Boudreau J, Gregoire TG, Margolis H, Næsset E, Gobakken T, Ståhl G. Estimating Quebec provincial forest resources using ICESat/GLAS. Can J For Res. 2009;39(4):862–81.

    Article  Google Scholar 

  46. [OMNRF] Ontario Ministry of Natural Resources and Forestry. Far North Land Cover Data Specifications Version 1.4. 2014. Accessed 15 Dec 2020.

  47. [OMNRF] Ontario Ministry of Natural Resources and Forestry. Forest Resources of Ontario. 2016. Accessed 15 Dec 2020.

  48. [OMNRF] Ontario Ministry of Natural Resources and Forestry. Far North Land Use Strategy: Discussion paper. 2019. Accessed 15 Dec 2020.

  49. [OMNRF] Ontario Ministry of Natural Resources and Forestry. Natural Resources Information Portal. 2020. Accessed 15 Dec 2020.

  50. [OMNRF] Ontario Ministry of Natural Resources and Forestry. Sustainable Growth: Ontario’s Forest Sector Strategy. 2020. Accessed 14 Jun 2021.

  51. Pan Y, Birdsey RA, Phillips OL, Jackson RB. The structure, distribution, and biomass of the world’s forests. Annu Rev Ecol Evol Syst. 2013;44:593–622.

    Article  Google Scholar 

  52. Parton WJ, Scurlock JM, Ojima DS, Gilmanov TG, Scholes RJ, Schimel DS, Kirchner T, Menaut JC, Seastedt T, Garcia Moya E, Kamnalrut A. Observations and modeling of biomass and soil organic matter dynamics for the grassland biome worldwide. Global Biogeochem Cycles. 1993;7(4):785–809.

    Article  CAS  Google Scholar 

  53. Price DT, Alfaro RI, Brown KJ, Flannigan MD, Fleming RA, Hogg EH, Girardin MP, Lakusta T, Johnston M, McKenney DW, Pedlar JH. Anticipating the consequences of climate change for Canada’s boreal forest ecosystems. Environ Rev. 2013;21(4):322–65.

    Article  Google Scholar 

  54. Reyer C, Lasch-Born P, Suckow F, Gutsch M, Murawski A, Pilz T. Projections of regional changes in forest net primary productivity for different tree species in Europe driven by climate change and carbon dioxide. Ann For Sci. 2014;71(2):211–25.

    Article  Google Scholar 

  55. Riahi K, Van Vuuren DP, Kriegler E, Edmonds J, Oneill BC, Fujimori S, Bauer N, Calvin K, Dellink R, Fricko O, Lutz W. The shared socioeconomic pathways and their energy, land use, and greenhouse gas emissions implications: an overview. Global Environ Change. 2017;42:153–68.

    Article  Google Scholar 

  56. Ruess RW, Hendrick RL, Burton AJ, Pregitzer KS, Sveinbjornssön B, Allen MF, Maurer GE. Coupling fine root dynamics with ecosystem carbon cycling in black spruce forests of interior Alaska. Ecol Monogr. 2003;73(4):643–62.

    Article  Google Scholar 

  57. Rustad LE, Fernandez IJ. Soil warming: Consequences for foliar litter decay in a spruce‐fir forest in Maine, USA. Soil Science Society of America Journal. 1998:1072–80.

  58. Schulze ED, Lloyd J, Kelliher FM, Wirth C, Rebmann C, Lühker B, Mund M, Knohl A, Milyukova IM, Schulze W, Ziegler W. Productivity of forests in the Eurosiberian boreal region and their potential to act as a carbon sink–-a synthesis. Glob Change Biol. 1999;5(6):703–22.

    Article  Google Scholar 

  59. Shrestha M, Acharya SC, Shrestha PK. Bias correction of climate models for hydrological modelling–are simple methods still useful? Meteorol Appl. 2017;24(3):531–9.

    Article  Google Scholar 

  60. Shvidenko A, Nilsson S. A synthesis of the impact of Russian forests on the global carbon budget for 1961–1998. Tellus B Chem Phys Meteorol. 2003;55(2):391–415.

    Google Scholar 

  61. Slivinski LC, Compo GP, Whitaker JS, Sardeshmukh PD, Giese BS, McColl C, Allan R, Yin X, Vose R, Titchner H, Kennedy J. Towards a more reliable historical reanalysis: Improvements for version 3 of the Twentieth Century Reanalysis system. Q J R Meteorol Soc. 2019;145(724):2876–908.

    Article  Google Scholar 

  62. Smith P, Bustamante M, Ahammad H, Clark H, Dong H, Elsiddig EA, Haberl H, Harper R, House J, Jafari M, Masera O. Agriculture, forestry and other land use (AFOLU). Climate change 2014: mitigation of climate change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Chapter. 2014. 11:811–922.

  63. Stocks BJ, Mason JA, Todd JB, Bosch EM, Wotton BM, Amiro BD, Flannigan MD, Hirsch KG, Logan KA, Martell DL, Skinner WR. Large forest fires in Canada, 1959–1997. Journal of Geophysical Research: Atmospheres. 2002;107(D1):FFR-5.

  64. Taylor KE, Stouffer RJ, Meehl GA. An overview of CMIP5 and the experiment design. Bull Am Meteor Soc. 2012;93(4):485–98.

    Article  Google Scholar 

  65. Todd-Brown KE, Randerson JT, Hopkins F, Arora V, Hajima T, Jones C, Shevliakova E, Tjiputra J, Volodin E, Wu T, Zhang Q. Changes in soil organic carbon storage predicted by Earth system models during the 21st century. Biogeosciences. 2014;11(8):2341–56.

    Article  CAS  Google Scholar 

  66. Turetsky MR, Kane ES, Harden JW, Ottmar RD, Manies KL, Hoy E, Kasischke ES. Recent acceleration of biomass burning and carbon losses in Alaskan forests and peatlands. Nat Geosci. 2011;4(1):27–31.

    Article  CAS  Google Scholar 

  67. Van Groenigen KJ, Qi X, Osenberg CW, Luo Y, Hungate BA. Faster decomposition under increased atmospheric CO2 limits soil carbon storage. Science. 2014;344(6183):508–9.

    Article  CAS  Google Scholar 

  68. Van Vuuren DP, Carter TR. Climate and socio-economic scenarios for climate change research and assessment: reconciling the new with the old. Clim Change. 2014;122(3):415–29.

    Article  Google Scholar 

  69. Vogel JG, Bond-Lamberty BP, Schuur EA, Gower ST, Mack MC, O’Connell KE, Valentine DW, Ruess RW. Carbon allocation in boreal black spruce forests across regions varying in soil temperature and precipitation. Glob Change Biol. 2008;14(7):1503–16.

    Article  Google Scholar 

  70. Walker XJ, Baltzer JL, Cumming SG, Day NJ, Ebert C, Goetz S, Johnstone JF, Potter S, Rogers BM, Schuur EA, Turetsky MR. Increasing wildfires threaten historic carbon sink of boreal forest soils. Nature. 2019;572(7770):520–3.

    Article  CAS  Google Scholar 

  71. Walker XJ, Rogers BM, Baltzer JL, Cumming SG, Day NJ, Goetz SJ, Johnstone JF, Schuur EA, Turetsky MR, Mack MC. Cross-scale controls on carbon emissions from boreal forest megafires. Glob Change Biol. 2018;24(9):4251–65.

    Article  Google Scholar 

  72. [WCRP] World Climate Research Programme, Coupled Model Intercomparison Project (Phase 6). 2019. Accessed 15 Dec 2020.

  73. Wirth C, Schulze ED, Lühker B, Grigoriev S, Siry M, Hardes G, Ziegler W, Backor M, Bauer G, Vygodskaya NN. Fire and site type effects on the long-term carbon and nitrogen balance in pristine Siberian Scots pine forests. Plant Soil. 2002;242(1):41–63.

    Article  CAS  Google Scholar 

Download references


We thank Lisa Buse (with the OMNRF) for editing an earlier version of this manuscript. We also thank three anonymous reviewers, whose thoughtful comments helped us improve the manuscript.


Funding for this study was provided by the Ontario Ministry of Natural Resources and Forestry (OMNRF) through Modelling Forest Biomass in the Far North project.

Author information

Authors and Affiliations



All authors were involved in the design of the project and discussed and refined the results. G. Mo and A. Gonsamo performed the data analyses and model simulations. All authors contributed to writing the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Michael T. Ter-Mikaelian.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Histotical and projected C stocks and fluxes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ter-Mikaelian, M.T., Gonsamo, A., Chen, J.M. et al. Historical and future carbon stocks in forests of northern Ontario, Canada. Carbon Balance Manage 16, 21 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: