Implications of future disturbance regimes on the carbon balance of Canada’s managed forest (2010–2100)

Recent increases in fire and insect disturbances have contributed to a transition of Canada’s managed forest carbon balance from sink to source. Further increases in area burned could contribute positive feedback to climate change. We made probabilistic forecasts of the recovery of C sinks in Canada’s managed forest between 2010 and 2100 under two assumptions about future area burned by wildfire: (1) no increase relative to levels observed in the last half of the 20th century and (2) linear increases by a factor of two or four (depending on region) from 2010 to 2100. Recovery of strong C sinks in Canada’s managed forest will be delayed until at least the 2030s because of insect outbreaks, even if predicted increases in area annually burned do not occur. After 2050, our simulations project an annual probability of a sink near 70% with no increase in area burned and 35% with increasing area burned. All simulations project a cumulative C source from 2010–2100, even if annual area burned does not increase. If the sink strength of terrestrial ecosystems is reduced because of increasing natural disturbances, then it will become more difficult to achieve global atmospheric CO2 stabilization targets.


Introduction
The strength and dynamics of carbon (C) sinks in terrestrial and ocean systems affect the residual airborne fraction of fossil fuel emissions remaining in the atmosphere. On average, 57% of annual CO 2 emissions were absorbed by these sinks between 1959(Le Quéré et al., 2009), but the residual airborne fraction appears to have increased over time because rates of emissions have increased faster than rates of uptake by global sinks (Canadell et al., 2007;Le Quéré et al., 2009). In many of the world's high-latitude forest ecosystems, it is not clear yet if the combined effects of climatic and global change impacts will tend to drive these ecosystems towards a sink or source Luo, 2007). In Canada, the managed forest (230 Mha) has recently transitioned to near C neutral or a small C source due to an unprecedented outbreak of mountain pine beetle in British Columbia (Kurz et al., 2008a) and is expected to remain so at least until 2020 under the continued risk of wildfire and the spruce budworm outbreak that is anticipated to occur in eastern Canada over the coming decade (Kurz et al., 2008b). Model projections and measurements both suggest post-disturbance recovery of C sinks in regrowing forests at the stand level (Amiro et al., 2006;Grant et al., 2007;Grant et al., 2009) and landscape level (Kurz et al., 2008a). Possible climatic and global change feedbacks Luo, 2007;Kurz et al., 2008c), could affect the time frame over which this recovery will occur. An important climate change feedback in Canada is a forecast increase in area burned by wildfire (Flannigan et al., 2005;Balshi et al., 2009a;Flannigan et al., 2009;Krawchuk et al., 2009;Podur and Wotton, 2010), leading to increases in direct greenhouse gas emissions from biomass combustion (Amiro et al., 2009;Balshi et al., 2009b). In this paper, we examine the potential impact of these forecast increases in area burned on the C balance of Canada's managed forest over the course of the 21st century. Determining if such an increase will delay or prevent the recovery of the C sink in Canada's managed forest is of both scientific and policy interest because it will affect the level of mitigation effort needed to reach future atmospheric CO 2 stabilization targets.

Carbon budget model of the Canadian forest sector (CBM-CFS3)
We conducted our simulations with the CBM-CFS3 (Kurz et al., 2009), an empirical model of stand-and landscape-level C dynamics that implements Tier 3 standards of the Intergovernmental Panel on Climate Change (IPCC) good practice guidance for reporting on C stocks and stock changes resulting from land use, land-use change and forestry (Penman et al., 2003).
It is the core model of Canada's National Forest Carbon Monitoring, Accounting and Reporting System (NFCMARS - Kurz and Apps, 2006), which estimates greenhouse gas (GHG) emissions and removals from managed forests in Canada for reporting under the United Nations framework convention on climate change and the Kyoto protocol. The model simulates above-and belowground biomass, detrital and soil C dynamics at the stand and landscape levels and accounts explicitly for human and natural disturbances. It tracks five categories of living biomass (foliage, merchantable wood, other wood, and coarse and fine roots) by species group (softwoods and hardwoods) and eleven detrital and soil C pools categorized by the type of material they contain and by their anticipated rate of decay. The pools represent snags (standing dead trees) and separate detrital and soil C into aboveground (woody litter and organic soil horizons) and belowground (mineral soil horizons) components. The pools decompose by pool-specific temperaturedependent decay rates per year (Smyth et al., 2010). Decay releases some C to the atmosphere and transfers the remainder to a stable, slowly decomposing pool from which decay releases all C to the atmosphere. At the stand level, aboveground biomass C is estimated from merchantable stem wood volume (Boudewyn et al., 2007); belowground biomass C from aboveground biomass C and species group (Li et al., 2003). The model simulates events that transfer C between pools or change the annual growth increment, including natural disturbances (e.g. fire, insects and windthrow) or forest management (e.g. harvesting, slash burning and planting). Disturbance matrices defining proportional C transfers between pools and fluxes to the atmosphere or forest products are used to simulate disturbance impacts in an event-driven framework. The model and its development are described in several publications (Kurz et al., 1992;Kurz and Apps, 1999;Kurz et al., 2009; and references therein) and is freely available with a user guide (Kull et al., 2006) and training material (http://carbon.cfs.nrcan.gc.ca).

National forest database
The NFCMARS and CBM-CFS3 use timber supply modelling data, including a forest inventory and merchantable volume yield curves, to describe and forecast C dynamics. The default source for this information in the NFCMARS is the Canadian National Forest Inventory (CanFI2001 - Power and Gillis, 2006), which has been replaced for many provinces with more recent detailed data used in provincial wood supply analyses. Monte-Carlo projections for 100-yr time periods with these comprehensive inventory data cannot presently be conducted due to computational restrictions. To overcome this constraint, we developed a compressed database that retains only the four forest inventory classifiers that are absolutely required by the model: (1) administrative boundary, (2) terrestrial ecozone, (3) a code distinguishing lands eligible for timber harvest and (4) leading tree species. Administrative boundary, terrestrial ecozone and management code are required to assign forest management activities, natural disturbances and ecological parameters to spatial strata. The leading species code is required to attach merchantable volume yield curves and volume to biomass conversion algorithms to inventory records and to specify eligibility for disturbance types. Area-weighted average yield curves were generated from the full database for records within these classifiers, based on the inventory area in 1990. This compression reduced the number of inventory records by one order of magnitude and assimilates data from over 20 project databases in the NFCMARS into a single database that can be executed as one CBM-CFS3 project. The results obtained from simulations with the compressed database mimic the Canada-wide results obtained for the full databases for several key indicators when evaluated over the period 1990-2007.

Projections
2.3.1. Fire. The Canadian large fires database (Amiro et al., 2001;Stocks et al., 2002) contains spatially explicit observations of annual area burned from 1959-1999. Annual area burned is highly variable, with about 97% of the area burned in Canada occurring in fires >200 ha. These large fires represent about 3.1% of the number of fires during that period. To simulate future fire (Kurz et al., 2008b) we developed from the observations in the large fires database, regionally calibrated probability distributions of annual areas burned for modelling regions derived by the intersection of provincial boundaries with terrestrial ecozones, for a total of 29 fire-modelling regions (Fig. 1). We considered any probability distribution for positive random values as a candidate (Evans et al., 2004) and chose the distribution with the lowest root mean squared error when distributions were fit using cumulative-distribution function regression (Cao, 2004), also fitting the criterion of a probability of an annual area burned greater than the largest in the recent historical record of between 1 and 2%. We also constrained the projected maximum burnable area by truncating values at twice the maximum observed in a modelling region during the recent historical period. Assumptions about the frequency and magnitude of extreme fire years are highly uncertain, but these rules typically generate fire cycles approximating those in the period of observation in long-term stochastic simulations. We assumed no temporal autocorrelation DISTURBANCES AND CANADA'S FUTURE FOREST CARBON BALANCE 721 Fig. 1. Map of Canada showing the managed forest and areas where fire was assumed to increase 2× and 4× in our simulations. (Armstrong, 1999) and also no spatial autocorrelation between regions. Process-based analysis of atmospheric circulation patterns suggests a negative correlation between area burned east and west of the Rocky Mountains (Macias Fauria and Johnson, 2006;Macias Fauria and Johnson, 2008), but otherwise it is uncertain if correlations observed in the last half of the 20th century (Magnussen, 2008) will endure in a future climate. Not accounting for this correlation likely reduced the estimated variance of our future projections (Magnussen, 2009).
Future increases in wildfire were implemented with a multiplier that we applied to the randomly drawn area burned. The multiplier increased linearly from one in 2010 to either two or four in 2100 to represent either a doubling or quadrupling of average annual area burned over that time. We assumed a doubling of fire for eastern Canada and British Columbia and a quadrupling of fire in western and northern Canada ( Fig. 1), based on projections in Flannigan et al. (2005); Balshi et al. (2009a) and Flannigan et al. (2009). Studies examining climate change impacts on future fire exist only for a small region in southern British Columbia (Nitschke and Innes, 2008a) and so we used a more conservative estimate (2×) for this jurisdiction, even though it is in western Canada. We assumed that the amount of fuel consumed by fire remained unchanged over time in both scenarios, as most of the increase in future emissions will likely be caused by increased area burned (Amiro et al., 2009). More sophisticated approaches for forecasting future area burned that account for interactions between climate, vegetation and fire are currently under development.
2.3.2. Harvest. To estimate future harvest levels, we started with forecasts from Kurz et al. (2008b), developed through consultations with provincial and territorial timber supply planners. We projected these forward to 2100 using the simple assumption that harvest rates would increase at a rate of 1% every 20 yr after 2040 in simple step functions (Fig. 2). Simple extrapolations usually outperform more complex methods for forecasting socio-economic time series when compared on the basis of predictive accuracy (Makridakis and Hibon, 2000). The drop in harvest rates in the first decade is the result in an anticipated reduction in salvage logging in British Columbia towards the end of the outbreak of mountain pine beetle. Future harvest assumptions were provided to CBM-CFS3 as area targets (Fig. 2). All harvest was simulated as clearcutting, using the C transfer parameters found in Kurz et al. (2009). In the simulations, oldest stands were selected for harvest first and 90% of a given database record was allowed to be harvested at one time. These rules most closely approximate those used in operational forest management in Canada.

Insects.
We included a single outbreak of mountain pine beetle (Dendroctonus ponderosae Hopkins) in British Columbia ending in 2018 and a single outbreak of spruce budworm (Choristoneura fumiferana Clemens) in eastern Canada (provinces from Ontario to the East) ending in 2022 in all simulations (Fig. 2). These represent average area forecasts from Kurz et al. (2008a). Recent monitoring results (Westfall and Ebata, 2010) confirm that observed area attacked by mountain pine beetle in BC to 2009 is well within the range of outbreak area forecasts in Kurz et al. (2008a). Spruce budworm was simulated using a series of annual area targets causing 15% stand mortality in stands between 40-and 250-yr of age (Fig. 2). We also examined a second scenario with a less severe spruce budworm outbreak, affecting the same area but causing only 9% stand mortality to evaluate the sensitivity of our results to assumptions about projected future spruce budworm impacts (Supplementary Material). A more detailed assessment of the potential regional impacts of spruce budworm is reported in Dymond et al. (In press). We did not simulate additional insect outbreaks after 2022, which is a conservative assumption.

Simulations and assumptions
We conducted two sets of 100 Monte Carlo simulations. In the first, we assumed that future area burned would follow the same statistical distribution as that observed in the last half of the 20th century. In the second, we assumed that annual area burned would increase between 2010 and 2100 by either a factor two or four depending on the modelling region (Fig. 1). The input data for each simulation in the second set were identical to the corresponding simulation in the first, with the exception that the annual area burned for each modelling region was increased by the annual multiplier derived through linear interpolation between one (in 2010) and two or four (in 2100). Here we report only changes in ecosystem C stocks, assuming implicitly and in accordance with current international GHG accounting guidelines that all C removed from the forest in harvested wood represents a reduction in ecosystem C stocks that was immediately emitted to the atmosphere.

Distribution of forest area by age class
The forest age-class structure in 2010 reflects the disturbance history of past decades. By 2100, both scenarios of future fire project more forest area younger than 60 yr and less forest area aged 60-160 yr, relative to 2010 (Fig. 3). The differences are greater for the simulations with increased area burned, but the shifts are in the same direction. Both scenarios of future fire show more forest area in age classes from 160 to 200 yr in 2100 than in 2010, due to the aging of forests in the 40-80 yr age classes and no difference in the area of forest older than 200 yr.

Annual C balance (2010-2100)
Our forecasts project that, on average, a C source will continue in Canada's managed forest until at least the 2030s for both scenarios of future area burned ( Fig. 4a and b). Recovery towards a sink begins to occur after that point and by about 2050 the annual C balance projections appear to stabilize. After 2050, the projections predict a median annual C sink of about 9 Tg C yr −1 (Fig. 4a) with an annual probability of a sink near 70% (Fig. 4c) if area burned does not increase. Typically, 80% of the simulations were between an annual source of 18 Tg C yr −1 and an annual sink of 21 Tg C yr −1 after 2050 if area burned does not increase.  If area burned increases, then the projections predict a median annual C source of about 8 Tg C yr −1 after 2050 (Fig. 4b), with an annual probability of a sink near 35% (Fig. 4c). Typically, 80% of the simulations were between a source of 66 Tg C yr −1 and a sink of 13 Tg C yr −1 after 2050 for an increase in area burned.

Cumulative C balance (2010-2100)
If area burned is not increased, then our projections forecast a median cumulative C source of 650 Tg C between 2010 and 2100, with 80% of the simulations showing a C source between 450 and 800 Tg C (Fig. 5). If area burned increases, then our projections forecast a median cumulative C source of 2200 Tg between 2010 and 2100, with 80% of the simulations showing a C source of between 1900 and 2500 Tg C (Fig. 5).

Total ecosystem C (2010-2100)
Both scenarios forecast a decline in total ecosystem C from 2010 to 2100. With no change in area burned, the forecast median total ecosystem C stock was 46.2 Pg C in 2100 and 80% of simulations had total ecosystem C stocks between 46.0 and 46.4 Pg C (Fig. 6a). With an increase in annual area burned, the forecast median total ecosystem C stock in 2100 was 44.6 Pg C and 80% of the simulations had total ecosystem C stocks between 44.3 and 44.9 Pg C. The difference between the median forecasts of total ecosystem C in 2100 between the two scenarios is 1.6 Pg C. The median increase in net ecosystem production (NEP) required to make up for the decline in total ecosystem C in 2100 due to increased annual area burned would be 7.8 g C m −2 yr −1 , with 80% of simulations requiring increases of between 6.7 and 8.8 g C m −2 yr −1 (Fig. 6b). This increase would need to occur each year between 2010 and 2100, for all 230 million hectares of managed forest in Canada. This represents an enhancement of 25%, relative to the average NEP (31 g C m −2 yr −1 ) observed in Canada's managed forest between 1990 and 2008 by the NFCMARS (Canadian Forest Service, unpublished data).

Discussion
Canada's forest was a C sink for much of the 20th century , but the managed forest has recently transitioned to near C neutral or a small C source and is expected to remain so at least until 2020 (Kurz et al., 2008b). In this paper, we used data from Canada's NFCMARS   and the CBM-CFS3 (Kurz et al., 2009) to forecast the potential recovery of a C sink in Canada's managed forest over the 21st century (2010-2100) under two potential scenarios of future annual area burned: (1) no change relative to observations in the last half of the 20th century and (2) climate-change induced linear increases by a factor of either two or four in different regions of Canada (Fig. 1). Our results suggest that recovery of strong C sinks in Canada's managed forest may be delayed until the 2030s, even if predicted increases in area annually burned do not occur. Our simulations project that the probability of an annual C sink is around 70% in the last half of the 21st century if there is no increase in area burned and only 35% if area burned increases. However, all simulations project that Canada's managed forest will be a cumulative C source from 2010-2100, even if annual areas burned does not increase. This occurs because a large source in one extreme fire year can reverse small sinks in several low fire years. The magnitude of the cumulative source is about three times larger if annual area burned increases than if it does not.
Our analysis did not consider all possible feedbacks resulting from climatic and global change effects. It is not yet certain if the balance of climate change effects on net primary production (NPP) and heterotrophic respiration (Rh) will act as a negative or positive feedback (Friedlingstein et al., 2006;Luo, 2007). We did not account for possible increases in NPP due to elevated atmospheric CO 2 concentration, increased atmospheric N deposition and longer growing seasons. At present, agreement among process-based simulation models about the magnitude (and direction) of forest responses to these factors remains limited (Friedlingstein et al., 2006), but observations over the last half of the 20th century suggest that these effects may have been enhancing C sinks where water is not strongly limiting (Boisvenue and Running, 2006). In short-term experiments, warming increases Rh. (Kirschbaum, 2000;Rustad et al., 2001). However, the long-term persistence of this enhancement is uncertain (Luo et al., 2001;Mellilo et al., 2002) and respiration may decline in drier soils (Allison and Treseder, 2008). Our analysis did not account for these effects, which remain a subject of inquiry (Davidson and Janssens, 2006;Bronson et al., 2008;Vanhala et al., 2008;Karhu et al., 2010). Our analysis may also be optimistic because we did not account for several potential source-side risks, including drought-induced forest mortality (Allen et al., 2009;van Mantgem et al., 2009), melting of permafrost (Schuur et al., 2009), increased emissions from peat fires  or insect outbreaks beyond 2022. In all likelihood additional insect outbreaks will occur over the 21st century because insects such as the spruce budworm have cyclical outbreak patterns (Blais, 1983;Royama, 1984;Candau et al., 1998). We also did not account for changes in postfire succession (Landhaeusser and Wein, 1993;Johnstone and Kasischke, 2005;Johnstone and Chapin, 2006) or regeneration failure (Nitschke and Innes, 2008b).
A previous study of a hypothetical boreal landscape suggested that long-term increases in NEP, over large geographic areas would be needed to offset the predicted impacts of increases in annual area burned and that this could only occur if NPP was enhanced relative to Rh (Kurz et al., 2008c). For increased NEP to compensate for C lost to increased area burned by wildfire (1.6 Pg C difference in 2100 between the median of total ecosystem C in our two simulations), the net effect of changes in NPP and Rh must result in a 7.8 g C m −2 yr −1 increase in NEP each year between 2010 and 2100, across all 230 million hectares of managed forest in Canada. This represents a 25% increase in NEP relative to the period 1990 to 2008. As both scenarios forecast a decline in C stocks, the increase would need to be even greater for total ecosystem C stocks to recover to 2010 levels. If we assume, unrealistically, that Rh does not change over time, then this also approximates the minimum increase in NPP required. The maximum increase in NPP required to offset the projected changes in C stocks cannot be evaluated using our DISTURBANCES AND CANADA'S FUTURE FOREST CARBON BALANCE 725 simulations because increased NPP also increases Rh and disturbance losses and because increased future temperature may also enhance Rh.
No simulation predicted an annual C sink of the magnitude that we estimate to have occurred in Canada during the 1990s, based on the NFCMARS. We believe that this is due to the legacy of age-class structure effects. In 1990, forests in Canada had a skewed age-class structure, with less forest in young age classes (0-40 yr) than in middle age classes (60-100 yr) . This age-class structure was the result of a period of high disturbances followed by a period of lower disturbances. The productivity of forests with such an age-class structure is high, but it is a transient and unsustainable condition (Böttcher et al., 2008). Productivity will decline as stands mature into older (less productive) age classes and as older forests are disturbed and returned to very young (again less productive) age classes (Böttcher et al., 2008). By 2100, even the simulations with no increase in area burned predict that total ecosystem C stocks will not have recovered to 2010 levels. This suggests that the residual effects of age-class structure on the forest C balance in Canada will continue over the 21st century. As the age-class structure of 2100 (skewed toward a predominance of younger forests) ages, forests may again temporarily recover to a strong sink status, particularly if there is a period of low disturbance. However, if further increases in the severity of disturbance regimes occur after 2100 and the age-class structure shifts to even younger forests, a return to a net sink may be delayed past that time.
We only included projection of insect disturbance impacts in the early portion of our simulation period -including the remainder of the ongoing mountain pine beetle outbreak in BC and one outbreak of eastern spruce budworm, an insect that has historically exhibited cyclical outbreak dynamics (Blais, 1983). We conducted a scenario analysis to explore the sensitivity of our results to the severity of spruce budworm outbreak simulated and found very limited sensitivity (Supplementary Material). Having no projection of insect outbreak dynamics beyond these two outbreaks, however, does make our results optimistic relative to outcomes that might be realized if there are severe insect outbreaks in the future that result in further reductions in productive capacity of the forest. Even greater increases in NEP (and hence NPP) would be required to offset these additional losses that were not accounted for in our projections.
Our results describe only C stored in the forest ecosystem and do not consider C stored in long-lived wood products or the effects of substituting wood-based products for more C intensive alternatives (Nabuurs et al., 2007;NCASI, 2007;Hennigar et al., 2008;Neilson et al., 2008). In the 20th century, about half of harvested wood C in Canada entered long-term storage . Forests remove C from the atmosphere, but they also provide ecosystem services and supply products such as timber, wood fibre and energy to meet society's demands. Atmospheric benefits of forest management in the 20th century may often have been a secondary benefit of our desire to maximize the yield of these products (Kauppi et al., 2010). The IPCC concluded that 'a sustainable forest management strategy, aimed at maintaining or increasing forest C stocks, while providing an annual sustained yield of timber, fibre and energy from the forest will generate the largest sustained mitigation benefit' (Nabuurs et al., 2007). Our projected C stock decreases suggest that the legacy effects of age-class structure may possibly overwhelm effects of forest management activities in the short term, even if annual area burned does not increase above levels observed in the 20th century. Forest sector mitigation options designed to provide atmospheric benefits need to be identified and quantified, and the atmospheric benefits accounted for relative to business-as-usual baselines (Böttcher et al., 2008). Recent economic conditions have significantly reduced the demand for wood products. Consequently, the volume of wood harvested in Canada has declined from 2004 to 2008, the last year for which statistics were available, to its lowest amount since 1990 1 and below what was forecast for 2008 in Kurz et al. (2008b) and in our simulations here. Future harvest levels are uncertain because future economic conditions are unpredictable. It is not known to what extent future demands on wood fibre for bioenergy may increase or substitute for present demands for products such as lumber and pulp and paper or if additional salvage logging of areas killed by insects will actually take place as we have assumed.

Conclusions
Recent global analyses suggest that the residual airborne fraction of global CO 2 emissions has been increasing (Canadell et al., 2007;Le Quéré et al., 2009). In the early part of the 21st century, Canada's managed forest shifted from a C sink to near C neutral or a small C source and is expected to experience an increase in the area burned by wildfire over the coming century. We forecast the recovery of C sinks in the 21st century under two assumptions about future area burned by wildfire: (1) no increase relative to levels observed in the last half of the 20th century and (2) gradual increases by a factor of two or four between 2010 and 2100. Strong C sinks in Canada's managed forest may be delayed until at least the 2030s and possibly later. We project that the probability of an annual C sink is around 70% in the last half of the 21st century, but only 35% if area burned increases. All simulations project that Canada's managed forest will be a cumulative C source from 2010-2100, even if annual area burned does not increase. Our results suggest that, relative to its size, Canada's managed forest is projected to make a disproportionately small contribution to reducing the airborne fraction of global fossil fuel emissions over the 21st century. In part, this result is due to the legacy effects of the age-class structure of Canada's forest in the late 20th century, when a period of relatively high disturbances followed by a period of relatively low disturbance resulted in a right-shifted age-class structure with high, but ultimately transiently enhanced, productivity and correspondingly high C sink strength. In regions where natural disturbance regimes are projected to increase, the resulting reduction in sink strength for these terrestrial ecosystems will make it increasingly more difficult to achieve global atmospheric CO 2 stabilization targets.

Acknowledgments
Funding for this project was provided by the Canadian Forest Service and JMM was supported by the Visiting Fellowships in Canadian Government Laboratories Program. We thank Drs. Mike Flannigan, Mike Wotton and Bill de Groot for discussions about representation of future area burned and in particular Mike Flannigan for helpful comments on an earlier version of this manuscript. We also thank Greg Rampley for maintenance of the simulation software and Michael Magnan and Scott Morken for assistance with the computing infrastructure used to conduct the simulations. Three anonymous reviewers provided valuable feedback on our manuscript, which we greatly appreciated.