How have past fire disturbances contributed to the current carbon balance of boreal ecosystems?

. Boreal ﬁres have immediate effects on regional carbon budgets by emitting CO 2 into the atmosphere at the time of burning, but they also have legacy effects by initi-ating a long-term carbon sink during post-ﬁre vegetation recovery. Quantifying these different effects on the current-day pan-boreal (44–84 ◦ N) carbon balance and quantifying relative contributions of legacy sinks by past ﬁres is important for understanding and predicting the carbon dynamics in this region. Here we used the global dynamic vegetation model ORCHIDEE–SPITFIRE (Organising Carbon and Hydrology In Dynamic Ecosystems – SPread and InTensity of FIRE) to attribute the contributions by ﬁres in different decades between 1850 and 2009 to the carbon balance of 2000–2009, taking into account the atmospheric CO 2 off ﬁres 0.17 ﬁres

are suspected to have collectively enhanced the vegetation production and carbon fixing. However, as climate change continues, carbon stocks in boreal forest may become more vulnerable, as indicated by (1) deceleration of "greening" over this biome as seen by satellites (Xu et al., 2013), (2) locally observed decreased vegetation productivity (Beck and Goetz, 2011), and (3) evidence for large climate-related disturbances such as insect outbreaks  and catastrophic fires (Kasischke and Hoy, 2012) that cause CO 2 losses to the atmosphere.
Fire has always been a natural disturbance in boreal ecosystems (Anderson et al., 2006), and it has multiple impacts on vegetation dynamics, carbon cycling, soil processes, atmospheric chemistry and permafrost dynamics. Fire plays an important role in the evolution of ecosystem species composition in this region through complex fire-climatevegetation feedbacks on different timescales (Kelly et al., 2013;Schulze et al., 2012). The carbon balance of boreal forests is modified immediately by fire through fire carbon emissions, but fires also lead to successional post-fire carbon accumulation as the ecosystem recovers -a long-term process of CO 2 removal from the atmosphere (Amiro et al., 2010;Goulden et al., 2011). Additionally, fires impact soil carbon dynamics, primarily by direct combustion of the organic layer at the soil surface but also through the creation and deposition of recalcitrant charcoal (Santín et al., 2015). Furthermore, organic soil carbon is also restored as post-fire vegetation carbon recovers , though the extent of restoration may depend on factors like post-fire vegetation type and regenerating forest stand density (Kashian et al., 2006). Lastly, soil carbon dynamics are also changed by altered soil temperature and moisture conditions after fire (Harden et al., 2006).
Many factors contribute to the currently observed boreal carbon sink, including the fertilization effect of increasing CO 2 concentration (Balshi et al., 2007), nitrogen deposition (DeLuca et al., 2008), forest management (Kauppi et al., 2010), climate change (Wang et al., 2011), and the balance between ecosystem (mainly forest) recovery from past disturbances (Pan et al., 2011b) and emissions from current fires. However, the relative contributions of these factors and their interactions are still poorly known, although a large part of the carbon sink in boreal forests has been attributed to forest recovering from past disturbance or degradation (Kauppi et al., 2010;Pan et al., 2011a). Given the role of fire in driving the demography and carbon balance of boreal forests, several studies used biogeochemical models to examine the carbon balance of boreal ecosystems and the related impacts from fires (Balshi et al., 2007;Hayes et al., 2011;Yuan et al., 2012). These studies conducted simulations with fire and without fire (or with a stationary fire regime) and examined the total-sum impacts of all preceding fires on the boreal carbon balance for a particular target time period. However, the immediate-source impacts of current fires through emissions and the sink legacies by previous fires were not formally sep-arated. Consequently, the contributions of fires that occurred before the current time (and associated post-fire vegetation recovery) to the current carbon balance, i.e., the legacy sink effects of past fire, remained largely unknown.
In the current study, we focus on the contributions of fires during different past periods to the carbon balance in boreal ecosystems. Theoretically, assuming stable environmental conditions, fires would have a close-to-zero net effect on the vegetation carbon storage over the fire cycle as the ecosystems are at a dynamic equilibrium state: fire emissions would be compensated for by post-fire vegetation regrowth (Kashian et al., 2006;Odum, 1969), as illustrated by the black curve in Fig. 1a. In this case, the forest net ecosystem production (NEP, which is photosynthesis minus respiration) may follow the classical temporal pattern, being negative in young forest, peaking in intermediately aged forest and declining in old forest. The temporal integration of NEP should be equal to the pulse of fire emissions, as the carbon balance over the entire fire cycle is expected to be zero.
However, when anthropogenic perturbations, especially those since preindustrial times as a result of intensive use of fossil fuels, come into play, this equilibrium state in which emissions are balanced by cumulative NEP may be disturbed. Of the anthropogenic perturbations affecting the environment, three prominent changes could exert a strong influence on the carbon dynamics related to disturbances. Climate change, predominantly temperature rise, could increase the growing-season length of Northern Hemisphere vegetation, strengthening plant physiological activities such as photosynthesis (Saxe et al., 2001). Atmospheric CO 2 increase could further enhance vegetation productivity, directly as a resource for photosynthesis but also indirectly by alleviating plant water stress (Franks et al., 2013). Nitrogen availability is considered as one limiting factor for boreal forest growth, and nitrogen deposition has been found to enhance vegetation productivity (Magnani et al., 2007). These three factors are abbreviated as CCN (climate, CO 2 , nitrogen) perturbations hereafter in this paper and are intended to represent the perturbations that collectively enhance the growth of vegetation regenerating after stand-replacing fires. As a result, the CCN perturbations could cause the curve of forest NEP against time since disturbance to shift toward higher carbon uptake, and the integration of NEP over time would probably exceed the fire emission pulse, making the vegetation a CO 2 sink (Fig. 1b, blue curve). Note here that, as fires are an agent leading to forest regeneration, the contributions of fires to the carbon balance are closely related to post-fire forest carbon dynamics and include the CCN perturbation effects that modify forest carbon uptake.
Based on this understanding, past fires must have contributed to the current boreal carbon balance through the enhanced post-fire forest regrowth as a result of CCN perturbations, termed the fire legacy carbon sink in this paper. The central aim of our study is to develop a conceptual framework to quantify the decadal contributions of past fires dur- , and carbon fluxes from undisturbed mature forests (with area being S − S). The nature (sink or source; blue or red arrow) and size (the width of arrows) of carbon balance of fire cohorts of different ages are shown quantitatively in the figure. The mathematical symbols for the carbon fluxes of the fire cohorts of the decades 2000-2009 and 1970-1979 and those from undisturbed mature forests are indicated; the symbols are the same as in Eq. (2) in the text. Note that, for clarity, the flux under preindustrial conditions (f c (g, b)) and the additional flux caused by CCN perturbations ( f c (g, b)) are not separated for all (red and blue) arrows that represent carbon fluxes.
ing 1850-2009 to the current carbon balance (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) in the pan-boreal region (44-84 • N). The tool used is the global dynamic vegetation model ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) with the prognostic fire module SPITFIRE (SPread and InTensity of FIRE). Fire occurrences are simulated in a prognostic way, with the dynamic vegetation module being activated. Our objectives are (1) to compare the simulated versus observed distribution of tree cover and tree groups, given fire disturbance; (2) to separate the contribution of legacy sink of past fires from emissions of current fires to the pan-boreal carbon balance and to further quantify the relative sink contributions by fires in different decades of the past. Being a preliminary effort, the different driving factors influencing fire contributions (such as CCN) are not individually separated; rather, their effects are included in the decadal fire contributions.

Model introduction
This study uses the process-based dynamic global vegetation model (DGVM) ORCHIDEE (Krinner et al., 2005). The ORCHIDEE model has three sub-modules. The SECHIBA sub-module simulates the fast exchange of water and energy between the land and the atmosphere. The STOMATE submodule simulates the vegetation carbon cycle processes including photosynthesis, photosynthate allocation, litter fall, litter and soil organic matter decomposition. The third submodule simulates vegetation dynamics. The equations of vegetation dynamics are mainly taken from the LPJ (Lund-Potsdam-Jena) model (Sitch et al., 2003), with modifications described by Krinner et al. (2005).
For this study, the prognostic fire module SPITFIRE as originally developed by Thonicke et al. (2010) was incorporated into ORCHIDEE, from here on referred to as ORCHIDEE-SPITFIRE. Global validation of simulated burned area and fire carbon emissions were described by Yue et al. (2014) and Yue et al. (2015). Notably, ORCHIDEE-SPITFIRE is able to capture the decadal variations of burned area in boreal Russia when compared to the historical reconstruction data by Mouillot and Field (2005) and the interannual variations of burned area in boreal North America when compared with the fire agency data. All fire processes are the same as described in Yue et al. (2014), except that the human suppression of lightning-ignited fires is introduced, as a function of human population density, following Li et al. (2012): where, D p is the population density (individuals per square kilometers), and F s a multiplicative coefficient applied to lightning ignitions to account for human suppression at a given D p . This corresponds to a suppression fraction of 0.01 in sparsely inhabited regions and of 0.99 in highly populated regions (i.e., D p → +∞). Within SPITFIRE, fire occurrence depends on vegetation and climate conditions and has feedbacks on forest mortality through crown scorching and cambial damage, which reduces forest stem density (Thonicke et al., 2010). Thus, in ORCHIDEE-SPITFIRE, vegetation dynamics are affected by both climatic factors, as simulated by the dynamic vegetation module, and fire disturbances, as simulated by SPIT-FIRE. In addition to the climatic limits that give the adaptation or extinction for different tree vegetation types under specific climate and climate variability conditions (Krinner et al., 2005;Sitch et al., 2003), fires further impact the treegrassland competition and the competition within woody vegetation types.
The ORCHIDEE-SPITFIRE used here includes the DGVM improvements made by Zhu et al. (2015), which improved the simulation of northern vegetation distribution. The improved DGVM processes include (1) tree mortality dependence on growth efficiency, defined as the ratio of net annual biomass increment to the preceding-year maximum leaf area index (LAI); (2) tree mortality induced by winter extreme coldness for all tree plant functional types (PFTs), except boreal deciduous needleleaf, and by spring frost in broadleaf forests only; (3) the definition of the tree line limit as an isotherm of a growing-season mean soil temperature of 6.7 • C. A threshold of a mean monthly temperature of 22 • C is used to limit the distribution of C4 grass, following Still et al. (2003). Maximum carboxylation rates (V c max , µmol m −2 s −1 ) were adjusted based on the results of parameter optimization for ORCHIDEE against flux tower measurements (Kuppel, 2012).

The conceptual framework
In this section we develop a conceptual framework which forms the basis of our simulation protocol and allows us to separate legacy carbon sinks from past fires for the carbon balance for the 2000-2009 decade from emissions by current fires. This conceptual framework was inspired by the theoretical attribution framework for the role of land-use change in carbon balance by Gasser and Ciais (2013). The influence of CCN perturbations on the carbon balance of regenerating forests as compared to a case without CCN is introduced in Sect. 1. Further, one should note that CCN perturbations also tend to increase carbon sinks in otherwise carbon-neutral old forests, i.e., land that is not disturbed by fires during the time of the CCN perturbation. Likewise, as the CCN perturbation increases forest carbon stock, when forests are burned, carbon emissions will also increase compared with the case without CCN perturbation. Consequently, for the decade of 2000-2009, the carbon balance of a grid cell is the sum of (1) fire emissions during 2000-2009, (2) legacy sinks caused by fires that occurred since 1850 and are impacted by CCN to various degrees (shown as the blue curve in Fig. 1a), and (3) source or sink of the tracts of forests that have not burned since 1850 but are influenced by CCN (i.e., which are considered undisturbed mature ecosystems). The composition of the carbon balance of 2000-2009 is illustrated in Fig. 1b.
The carbon balance of a geographical area covered by a given biome (g, b) for the 2000-2009 decade, under the CCN perturbation and taking into account decadal fire disturbances since 1850, can be expressed as where F ON (g, b) is the total carbon balance of the area S(g, b), typically expressed in grams of carbon per year, with presence of fire, and all lowercase f functions indicate the area-based carbon balance expressed as grams of carbon per square meter per year for various cases: f * u (g, b) is the undisturbed land impacted by the CCN perturbation (thus not equal to zero); f c (g, b) is the fire-generated cohort carbon flux density without the CCN perturbation; and f c (g, b) is the deviation of carbon flux from a cohort under steady environment conditions because of the CCN perturbation ( Fig. 1a blue curve). δS i represents the fire-disturbed land cohorts within the ith decade, with i ranging from the 1850s (1850-1859) to the 2000s (2000-2009); S(g, b) is the sum of disturbed land areas from fires of all decades since 1850. Note that, in Eq. (2), we separated the total carbon flux into lands undisturbed and those disturbed by fire. Further, we assume that fires also occurred before 1850, but their influence on the 2000-2009 carbon flux is included in the undisturbed land flux, given the observed very small net ecosystem produc-Biogeosciences, 13, 675-690, 2016 www.biogeosciences.net/13/675/2016/ tivity in boreal forests older than 150 years (Goulden et al., 2011).
In studies using numerical biogeochemical models, Eq.
(2) represents a case in which fire-generated forest cohorts are explicitly simulated: the 2nd part on the right-hand side of the equation gives the contributions of different decadal fires to the carbon balance for the 2000-2009 decade. However, for models that do not explicitly simulate forest cohorts (which is the case for the version of ORCHIDEE used here), a workaround is possible by manually suppressing fires in the model within a particular decade to allow quantifying the contribution of fires from this decade by the difference between the two simulations. Similar to Eq. (2), the carbon flux for the 2000-2009 decade if fires are suppressed in a particular decade D can be written as is the carbon balance for the 2000-2009 decade but with fires suppressed in the D decade and with the contribution by fires of the D decade being simultaneously removed from the right-hand side of the equation. Thus, the contribution by fires of the D decade is the difference between F ON (g, b) and F OFF, D (g, b): where Cont D is the contribution of fires within the D decade to the carbon balance of the 2000-2009 decade. In contrast with explicit cohort simulation, this factorial approach quantifies the past-fire-generated cohort contribution, taking as a baseline the carbon flux of otherwise undisturbed land but as influenced by the CCN perturbation. Finally, one could vary D from the 1850s to the 2000s to derive the contribution by fires within each decade between 1850 and 2009. This conceptual framework remains valid when integrating all the variables in Eqs.
(2)-(4) over the geographical extent and different vegetation types to attribute carbon fluxes on a regional scale. Note that, in this framework, the effects of different factors of the CCN perturbation are not individually separated, but rather their impact is embedded as a whole in the fire contribution.

Simulation protocol and input data sets
Following the conceptual framework, we conducted factorial simulations to quantify the decadal contributions of past "fire cohorts" to the simulated carbon balance of 2000-2009. The carbon balance is defined as the net biome production (NBP): where NPP is net primary production (i.e., the net biomass accumulation by plants after accounting for their own use), RH is the ecosystem heterotrophic respiration, and EMI is carbon released by fire. A positive NBP indicates a net carbon flux from the atmosphere to land, i.e., a land carbon sink.
In the following, we use the terms "carbon sink" and "NBP" interchangeably, unless otherwise specified (e.g. if stated as a negative NBP, it is a carbon source releasing carbon to the atmosphere). We conducted a reference simulation (SIM fireON ) from 1850 until 2011, accounting for climate change, atmospheric CO 2 concentration change and prognostically simulated fire disturbance. We then conducted a series of other simulations (named SIM OFF ), which branch off from the SIM fireON simulation from the beginning year of each decade between 1850 and 2009. In the SIM OFF simulations, the fire module was switched off sequentially from the decade of the 1850s (1850-1859) to the 2000s (2000-2009) and switched on afterwards, with all remaining parameter settings and input data sets the same as in the reference simulation. Following Eq. (4), the contribution by fires within a specific decade to the carbon balance of each year for the time after this decade would be quantified as the difference between the reference simulation and the decadal SIM OFF simulation. In all simulations, the vegetation dynamics module of ORCHIDEE was switched on to allow the vegetation distribution to respond to climate variations and fire disturbances.
The spatial domain of our simulation covers the land pixels of 44-84 • N at a 2 • resolution. The land north of 84 • was excluded as it is covered mainly by ice and snow. The model was forced by the CRUNCEP climate data at a 2 • resolution, regridded from its original resolution of 0.5 • . The CRUNCEP consists of 6-hourly gridded climate data generated by combining CRU TS 3.1 0.5 • monthly climate data and NCEP 6-hourly 2.5 • reanalysis data (thus the name CRUNCEP). Rainfall, cloudiness, relative humidity and temperature are from the CRU data set and interpolated at a 6-hourly time step following the temporal variability of NCEP. Pressure, longwave radiation, and wind speed are from NCEP, reinterpolated on a 0.5 • scale. The values for these variables for 1948 were also used for the period before 1948. For more details, see http://dods.extra.cea.fr/ store/p529viov/cruncep/V4_1901_2012/readme.htm. A single global annual atmospheric CO 2 concentration time series since 1850 was applied everywhere in the spatial domain of the model, which is a combination of ice core and NOAA station measurements. The fire module needs additional input data for lightning flashes and human population density. Lightning flashes were retrieved from the High Resolution Monthly Climatology of lightning flashes by the Lightning Imaging Sensor-Optical Transient Detector (LIS/OTD) (http://gcmd.nasa.gov/records/GCMD_ lohrmc.html). The LIS/OTD data set provides annual mean flash rates over the period of 1995-2000 on a 0.5 • scale with monthly time step, which was cycled each year throughout www.biogeosciences.net/13/675/2016/ Biogeosciences, 13, 675-690, 2016 the simulation. An annual historical population density map was retrieved from the Netherlands Environmental Assessment Agency (http://themasites.pbl.nl/tridion/en/themasites/ hyde/download/index-2.html). Both lightning and population density data sets were regridded at a 2 • resolution before being fed into the model. The reference simulation SIM fireON consists of a spin-up run from bare soil and a transient run, with the fire module being activated. For the spin-up, climate data for the period 1901-1930 were cycled and atmospheric CO 2 concentration (285 ppm) and population density were prescribed at the 1850 level. The spin-up run lasted for 400 years but contained three runs of soil-only processes each lasting 1000 years to speed up reaching equilibrium for slow and passive soil carbon pools. We verified that the average annual NBP during the last 30 years of the spin-up run was −0.003 Pg C yr −1 (a negative value as the model recovers from fast accumulation of soil carbon in the soil-only runs) and that no significant trend exists for annual NBP, indicating that the model had approximately reached an equilibrium state. The spin-up was followed by a transient simulation for 1850-2011, in which transient climate data, atmospheric CO 2 concentration and population density data were used. For 1850-1900, cycling climate data of 1901-1930 continue to be used.
As our focus is the carbon dynamics of natural vegetation in response to fires within the boreal region, croplands were not simulated in the model. This is acceptable given that land-use change during the 20th century in this region was small (Hurtt et al., 2006). Cropland fractions within grid cells were prescribed according to a current-day vegetation map (the IGBP-DIS 1 km global land-cover map; Loveland et al., 2000), and fractions of natural vegetation (i.e., trees and grasses) were simulated. Tundra in the high-arctic regions is simulated as C3 grassland.

Comparison of simulated forest distribution and fires to observations
We compared the spatial distribution of three morphological and phenological tree groups between the model simulation and MODIS land-cover data for the year 2010: broadleaf (including evergreen and deciduous), evergreen needleleaf and deciduous needleleaf trees, corresponding to the three boreal tree PFTs in ORCHIDEE. The MCD12Q1 version 5 landcover data (Friedl et al., 2010) were used (http:glcf.umd. edu/data/lc, with a northern limit of 84 • N). Fractions of the 17 different land-cover types in the IGBP land classification scheme were calculated at a 2 • resolution based on the 500 m original resolution data. Further, the 2 • land-cover fractions were cross-walked to PFT fractions using the approach developed by Poulter et al. (2011), in which the mixed treegrass land-cover types such as shrublands are assumed to be composed of different fractions of trees and grasses (see Table 6 in Poulter et al., 2011, for more details). The simulated maximum foliage projective cover for each of the three tree groups was compared with the corresponding MODIS observation, with the sum of the three groups being compared as tree cover. Simulated burned area and fire carbon emissions were compared with GFED3.1 burned area data , and carbon emission estimates were simulated by the CASA biosphere model . Burned areas and fire carbon emissions from agricultural fires were excluded from GFED3.1 data before comparison because these fires are not included in the model. Northern peatland fires were not simulated due to a lack of peatland PFT in the model nor are they included in the GFED3.1 emission data.

Simulated forest distribution
The simulated spatial extent of forest distribution is broadly similar to that of MODIS land-cover data over the region north of 44 • N for the year 2010, with the forest biome extending from eastern Canada northwestward to Alaska in boreal North America, and to that in northern and northeastern Europe, as well as most of Siberia (Fig. 2). The magnitude of foliage projective tree cover between ORCHIDEE and MODIS land-cover data is generally comparable, except at the southern and northern fringes of the study region (mainly Asia and America), where tree cover is overestimated by approximately 30-50 % in ORCHIDEE (hatched areas in Fig. 2). When considering the uncertainties in different observation data sets (by comparing different land-cover data sets of ESA-CCI, GLC2000 and VCF; see the Supplement for more details on data sources and their treatment), the error in simulated tree cover is less prominent (Supplement Fig. S1). The over-or underestimation of tree cover by ORCHIDEE in central and northern Siberia disappears; however, the overestimation of tree cover in southern Asian and North American boreal forests remains. In central Alaska and western Canada, tree cover is also underestimated by 10-30 % of ground area. Figure 3 presents simulated and observed spatial distribution of three tree groups: broadleaf (including evergreen and deciduous), evergreen needleleaf and deciduous needleleaf. There is a widespread presence of broadleaf forest, but with generally low fractional cover, across the study region, which is reproduced fairly by ORCHIDEE (Fig. 3, panel 1a and b). Both MODIS land-cover data and ORCHIDEE simulation indicate the dominance of evergreen needleleaf forest in North America, in western Siberia, and in northern and eastern Europe (Fig. 3, panel 2a and b). In contrast, MODIS data show that central and eastern Siberia is dominated by deciduous needleleaf forests (Fig. 3, panel 3b). ORCHIDEE successfully captures this, but the spatial extent and magnitude of tree cover are overestimated (Fig. 3, panel 3a). In addi-  tion, ORCHIDEE also erroneously allocates more deciduous needleleaf forests in Alaska and northwestern Canada than the MODIS data. We also extend the comparison of different tree group extents by including more land-cover data sets (see Figs. S2, S3 and S4). Again, when considering other land-cover maps (ESA-CCI, GLC2000 and VCF), the model error is less than when using the MODIS data set. Notably, both ESA-CCI and GLC2000 data sets indicate a larger extent of deciduous needleleaf forest in eastern Siberia compared to MODIS, resulting in much lower errors in the OR-CHIDEE simulation (nevertheless, a model overestimation of 20-50 % of ground area persists in western Siberia).

Simulated burned area and fire carbon emissions
The spatial distribution of simulated mean annual burned fraction for 1997-2009 is compared with GFED3.1 data in Fig. 4, with non-modeled agricultural fires being excluded from GFED data. The comparisons of cumulative latitudinal distribution of burned area and fire carbon emissions are shown in Fig  area exist, ORCHIDEE-SPITFIRE simulates an annual total burned area of 11.9 Mha yr −1 and fire carbon emissions of 0.20 Pg C yr −1 , which are close to GFED3.1 estimates giving an annual burned area of 16.9 Mha yr −1 and fire carbon emissions of 0.20 Pg C yr −1 . Spatially, burned area is underestimated within the latitude band 44-54 • N in Eurasia, concurrent with an overestimation of tree cover in the same region ( Figs. 2 and 3). On the other hand, there is an overestimation of burned area in the regions north of 54 • N covered by forest, shrubland and tundra according to the MCD12Q1 landcover map. Over North America, the spatial distribution of simulated burned area is in fair agreement with the GFED3.1 data, with burned area being dominated by the northwest-tosoutheast boreal forest fires.

Decadal contributions of fire to the simulated carbon sink
The simulated annual NBP for 1850-2011 for the study region in non-agricultural land and contributions of decadal fire cohorts to the carbon balance after the fire occurrence are shown in Fig. 6. The annual carbon sink of the reference simulation for 1990-2011 is 0.91 Pg C yr −1 (Fig. 6a), which falls within the range of forest-inventory-based estimates (∼ 0.7 Pg C yr −1 by Pan et al., 2011b) and the mean value of the terrestrial carbon cycle models (∼ 1.1 Pg C yr −1 ) as assessed by IPCC AR5 . Figure 6b shows how each decadal fire cohort contributes to the NBP of the study domain. For example, the curve labeled "1910s" shows the annual contribution of the cohort of the decade 1910-1919, which produced a net carbon source, followed by a long-term carbon sink whose magnitude decreases with time.  sizes, whereas fires within the 2000-2009 decade contribute as a source term. Figure 7 shows the contributions of fires within each decade to the annual NBP of the study region for 2000-2009. All decades before 2000 cause a fire legacy sink, collectively having a total sink of 0.23 Pg C yr −1 . These legacy sinks are compensated for by a carbon source of 0.17 Pg C yr −1 from fires in 2000-2009, leaving a net fire effect of 0.06 Pg C yr −1 . This net sink fire effect represents only a very small fraction (6.3 %) of the annual carbon sink of the reference simulation (0.95 Pg C yr −1 ), indicating that most of this sink occurs in unburned natural ecosystems for which the model produces enhanced carbon storage due to climate warming (e.g., longer growing seasons) and the CO 2 fertilization effect. The sink contributions of different decadal fire cohorts (1850-1999) exhibit a general decaying trend as the cohort ages, with the variations being affected by changes in climate, atmospheric CO 2 concentration and fire disturbance. Fires in the 4 decades prior to 2000-2009 (1960-1999, i.e., corresponding to a "cohort age" of 10-40 years) collectively contribute 0.14 Pg C yr −1 , accounting for 61 % of total legacy sink effect. Fires in the past century (1900-1999) contribute 0.19 Pg C yr −1 , or 83 %, of the total legacy sink.
The whole study region can be classified into six fire groups according to their different fire return intervals (FRIs, here quantified as the inverse of burned fraction) as simulated by the model, with the shortest FRI of 2-10 years and the longest of more than 500 years. This classification was done for each decade of 1850-1999 (i.e., decades having a carbon sink effect for 2000-2009), using a simulated mean decadal burned fraction, followed by partitioning the decadal sink contribution into these fire groups. Figure 8 shows the relative contributions of each fire group by summing together the partitioning results of all the decades. The fire group with an FRI of 10-50 years emerges as the biggest contributor, contributing a carbon sink of 0.1 Pg C yr −1 or 42.7 % of the total sink effect. Fires with intermediate FRIs (50-200 years) contribute 0.06 Pg C yr −1 (26.1 % of the total sink effect), while very rare fires (with an FRI > 500 years) or very frequent fires (with an FRI of 2-10 years) contribute least to the total sink effect (collectively contributing 0.04 Pg C yr −1 or 15.6 % of the total sink effect).

Discussion
We first describe in general the fire-climate-vegetation feedbacks in boreal regions and the role of fires in the regional carbon balance to put our findings in a more appropriate context (Sect. 4.1). Section 4.2 discusses some general model performance issues, with Sect. 4.3 presenting more detailed comparisons of our results with similar studies. Section 4.4 discusses uncertainties and future perspectives.

Boreal fire-climate-vegetation feedbacks and fire contribution to the regional carbon balance
In boreal regions the climate, vegetation dynamics and fire disturbances are intrinsically linked with each other (Campbell and Flannigan, 2000). Given the long time of exposure under insolation during summer days, fuels (e.g., litter on the ground) could get dry enough for fires to start if there are enough consecutive days with little precipitation.  In turn, plant traits adapt for fires, and fire adaption is used as a strategy to maintain competitiveness by different tree species (Wirth, 2005). For example, the gradual rising of black spruce (Picea mariana) in place of Betula in Alaskan forests during the Holocene has been aided by increased fire activities as a result of climate warming since the last glacial maximum (Kelly et al., 2013), since spruce trees keep their dead branches to promote fires and have serotinous cones that geminate after fire, making them more competitive against Betula under increasing fire disturbances.
Given a stable fire regime (fire return interval, fire severity, etc.), spruce forests form stable self-replacement succession cycles: carbon stored in fuels (litter and crown fuel) is released into atmosphere during fire; young forest stand is regenerated, and surface organic litter and biomass carbon stock are restored during forest growth until the next fire event . At the early successional stage, deciduous broadleaf trees (aspen, birch) often occur as pioneer species and are outcompeted at the late-successional stage due to their shade intolerance (Johnstone et al., 2010b). As such, fire cycles are internally coupled with vegetation carbon dynamics (and hydrological and energetic dynamics). As most carbon in boreal ecosystems is stored in organic soil, which is the dominant source of fire carbon emissions, fires have a comparatively big impact on the vegetation carbon cycling . However, evidence shows that more intense fires could sustain the dominance of broadleaf trees for a longer time and had the potential to alter the regional vegetation composition (Johnstone et al., 2010a).
With growing atmospheric concentrations of greenhouse gases and anthropogenic warming of the climate during past decades, there is increasing interest in examining boreal ecosystems as a potential carbon sink and, especially, in how likely it is that increasing fire activities would impact the carbon dynamics of this region. Research foci include quantifying contemporary regional fire carbon emissions (French et al., 2011), site-level post-fire carbon dynamics (Goulden et al., 2011), and regional carbon balance analysis using largescale biogeochemical models (Balshi et al., 2007;Hayes et al., 2011). The large-scale biogeochemical models have the particular advantage of evaluating the carbon balance on the regional scale and separating the impacts of different environmental factors such as climate, atmospheric CO 2 and disturbances. Most modeling studies examined the impacts of a changed fire regime or the collective impact of past fires Biogeosciences, 13, 675-690, 2016 www.biogeosciences.net/13/675/2016/ on the carbon balance for a target period. Bond-Lamberty et al. (2007) found that the central Canadian boreal forest is a small carbon sink (9.9 ± 11.8 g C m −2 yr −1 ) for 1958-2005, and, compared to a stable fire regime of the mid-20th century, fire disturbances have reduced the sink by 8.5 g C m −2 yr −1 . Balshi et al. (2007) and Hayes et al. (2011) used additive biogeochemical model simulations (i.e., simulations with and without fire) and quantified the collective impact of past fires on the pan-boreal carbon balance for different decades of the second half of 20th century, with fire contribution varying from small source to sink effects (around 0.1 Pg C yr −1 ) depending on different time periods. Nevertheless, given increasing fire frequency during the second half of the 20th century in this region (Stocks et al., 2003) and the important post-fire vegetation carbon dynamics linked with anthropogenic perturbations (such as the CCN perturbations as introduced in Sect. 1), few studies have tried to examine the potentially different impacts from fires occurring at different times in the past and elucidate how the current pan-boreal carbon balance is determined by past fire legacy sinks and current-day fire carbon emissions. Using a factorial simulation protocol, we found that fires during 2000-2009 have a net source contribution of −0.17 Pg C yr −1 to the carbon balance of the decade 2000-2009. However, this source effect is compensated for by legacy sinks (in total 0.23 Pg C yr −1 ) in land recovering from fires prior to the 2000s (1850-1999). These legacy sinks are ameliorated by climate warming and CO 2 fertilization. We further found that more than 60 % of the sink effects are contributed by fires during 1960-1999. Our finding is unique in that it separates the effects of previous fire legacy sinks and current-day fire emissions.

General model performance, simulated vegetation dynamics and burned area
ORCHIDEE-SPITFIRE successfully captured the largescale spatial pattern of tree cover distribution and the distribution of broadleaf versus needleleaf and evergreen versus deciduous forests in different continents, with the presence of fire disturbances being prognostically simulated. The larger spatial extent of deciduous needleleaf forests in Siberia and northern regions of America in ORCHIDEE may be due to our DGVM parameterization according to which winter extreme coldness leads to elevated mortality of all forests except deciduous needleleaf ones; this expands their presence within the tree line limit as represented by an isotherm of growing-season soil temperature (Zhu et al., 2015). Schulze et al. (2012) found that in a transitional zone (61-64 • N, 90-107 • E) in central Siberia, where the species Picea obovata and Abies sibirica (evergreen conifers) are natural late-successional species, frequent surface fires are the major factor explaining the dominance of Larix over the evergreen climax tree species. Infrequent crown fires initiate new Larix cohorts, while surface fires thin them and prevent ev-ergreen needleleaf saplings from reaching the canopy. Even though our model does not account explicitly for these two different fire impacts, on a broad scale, the dominance of evergreen coniferous forests in northern Europe and western Siberia coincides with slightly lower fire frequencies (Figs. 3 and 4). This is consistent with the observed pattern that more frequent fires in eastern Siberia are associated with the dominance of Larix deciduous needleleaf trees.
For the majority of the pan-boreal region, ORCHIDEE-SPITFIRE simulates a fire return interval of 10-200 years (Fig. 4, corresponding to burned fraction of 0.5-10 %), which is consistent with the evidence from various observational data sets Stocks et al., 2003). The simulated fire frequency (0.2-2 % yr −1 ) in Canada agrees with that reported by Stocks et al. (2003) using the Canadian Large Fire Database. The general spatial extent and magnitude of fires in northern Eurasia (> 54 • N) roughly agrees with GFED3.1 data, although burned fractions in northern tundra and shrubland are overestimated. This may be because tundra is treated as generic C3 grass in the model and thus assigned a low fuel bulk density (Thonicke et al., 2010) that promotes fast fire propagation. In reality tundra has a more dense growth form than temperate grasslands and therefore has a much higher bulk density (Pfeiffer et al., 2013). Fires are greatly underestimated by the model at the southern edge of the study area in Eurasia, with a simulated burned fraction of 0.2-2 % compared to values of 1-30 % in GFED3.1 data. This underestimation, especially in central Asian grasslands over Kazakhstan and Mongolia, is accompanied by an overestimation of tree cover (Fig. 2). This indicates that the role of fires in promoting grasslands over forests as shown by other modeling studies (e.g., Bond et al., 2005;Poulter et al., 2015) in these semiarid regions is underestimated in ORCHIDEE-SPITFIRE, probably due to excessive tree sapling recruitment. Despite this, our simulated boreal carbon sink for the 1990-1999 and 2000-2009 decades is comparable with other independent approaches, with simulated fire carbon emissions being close to GFED3.1 data. Therefore, though spatial model errors exist, we believe that the quantified total carbon fluxes on the regional scale remain valid. Balshi et al. (2007) and Hayes et al. (2011) used an additive simulation protocol to examine fire impact on the carbon balance, i.e., the contribution of fire to the carbon balance of a target decade (e.g., the 2000s) is given by the difference between two simulations, with and without fires. Note that this approach examines the total sum effect of all fires occurring before but also within the target decade, i.e., equivalent to the effect of all fires of 1850-2009 and termed net fire effect in our analysis. Balshi et al. (2007) further conducted parallel simulations with and without CO 2 fertilization for all additive runs. They found that during 1996-2002, the sum effect of fires in the pan-boreal region (north of 45 • N) increased the ecosystem carbon storage (ranging from 0.08 to 0.5 Pg C yr −1 ) for all years except 2002, according to a simulation that includes the CO 2 fertilization effect. When the CO 2 fertilization effect is excluded, the role of fires is more varied, leading to a close to zero sum fire effect for the same period. We also found that the net fire effect during the 2000-2009 decade to be a carbon sink of 0.06 Pg C yr −1 (i.e., the equivalent of the sum fire effect in Balshi et al., 2007), a value smaller than that reported by Balshi et al. (2007). However, we noticed that in their study, the contribution of fires varied greatly in magnitude from year to year, and it was sometimes even 3 times higher than the sink term due to the CO 2 fertilization effect, which may indicate the great uncertainty in their results (Fig. 6 in Balshi et al., 2007). Hayes et al. (2011) also used the additive approach to find a net carbon sink fire effect on the pan-boreal carbon balance for the decades of 1960-1969 to 1990-1999 with a similar magnitude to that in our study (0.03-0.08 Pg C yr −1 ). They argue that fires have changed from a carbon sink to source term for the 2000-2009 decade (ca. −0.13 Pg C yr −1 ) due to increased fire activities (Fig. 3 in Hayes et al., 2011), which is different from our conclusion. However, it should be noted that their estimated pan-boreal carbon sink for 1997-2006 (0.04 Pg C yr −1 ) was much lower than those based on atmospheric inversion or inventory approaches . On the other hand, their estimated fire carbon emissions (0.3 Pg C yr −1 for north of 45 • N) are 50 % higher than GFED3.1 data. Thus, it is likely that the biases in their estimated carbon fluxes (overestimation of emissions and underestimation of carbon sink) could lead to an overestimation of the carbon source effect by fires in the 2000-2009 decade. Finally, Yuan et al. (2012) examined the effect of changes in fire regime on the carbon balance of the Yukon River basin forests in Alaska from 1960 to 2006 by comparing simulations with time-varying and fixed fire regimes. They found that increased fires, compared with a stationary fire regime, have reduced the total ecosystem carbon storage by 185 Tg C, or 4 Tg C yr −1 . Despite not using exactly the same simulation approach, we also found a net carbon source fire effect of 1.5 Tg C yr −1 for the 2000-2009 decade carbon balance for Alaska, similar to Yuan et al. (2012) but with a smaller magnitude.

Comparison of simulated fire impacts with other studies and fire contributions linked with burned area and fire frequency
The sink contributions by different decadal "fire cohorts" show a general decreasing trend in the past, with more than half of the total sink effect contributed by the 4 decades before 2000 . This pattern may be partly explained by the strong carbon uptake in the young to medium-aged forests, as shown by site-level measurement (Goulden et al., 2011) and partly reflected in the model (Fig. 6b). One may consider whether the sink magnitude could be related to the amount of burned area, as suppressing of strong fire may lead to strong recovery (and thus a strong legacy sink). As shown in Fig. S5, the variation in decadal sink contribution magnitude does not echo that of burned area exactly, despite the fact that the correlation does exist (r = 0.54, p<0.05). Thus, we suspect that the variation in decadal fire legacy sinks may be related with both the known temporal pattern of post-fire forest carbon uptake and the fire extent. The CCN perturbations (represented in the model by applying transient climate forcing and increasing atmospheric CO 2 ) must also exert some control, but the full separation of their impacts is beyond the scope here.
We also found the highest legacy sink is contributed by the fire group with a fire return interval of 10-50 years (0.10 Pg C yr −1 , or 43 % of the total sink effect), followed by the 100-200-year fire group (0.04 Pg C yr −1 ) and 50-100 years (0.03 Pg C yr −1 ). In fact, the highest contribution by 10-50-year fire group is related to its dominance in total burned area (58 % of the total burned area of all fire groups) (Table S1 in Supplement). When examining the ratio of legacy sink effect to burned area (somewhat like fire sink efficiency), the 100-200-year and 200-500-year fire groups emerge to have the highest ratio (0.037 Pg C Mha −1 ). This ratio is reasonable as fires with this long return interval often occur in forest (or tundra, but more rarely) that has a strong and long-term recovery carbon uptake. The ratio of sink to burned area decreases as the fire return interval increases, indicating more frequent fires leading to weaker sink recovery, probably because increasing fire frequency is associated with increasing grassland fraction (Yue et al., 2014), which has a weaker sink recovery than forest. It is hard to conclude that more frequent fires will necessarily lead to a stronger sink effect. However, in general, if the same vegetation type could be maintained (e.g., forest regenerating after fire) rather than more intense fire leading to the replacement of forest by grassland, then, combined with the CCN perturbations and the strong carbon uptake of young to mediumaged forest, vegetation carbon uptake may increase with increasing fire frequency.
We highlight important contributions of past fire disturbances to the current ecosystem carbon sink, thanks to postfire vegetation recovery being enhanced by CO 2 fertilization and climate warming. These two factors, in spite of their roles not having been entirely separated out in the current study, may also influence the occurrence of fires and their emissions in the 2000-2009 decade, which partially counteract the sink effects by previous fires. In the long term, change in ecosystem structure and species will also affect fuel load and combustion completeness and modify fire emissions as well. Therefore, the future role of fires in the carbon balance of boreal regions remains rather uncertain and depends on how the post-fire recovery sink and fire carbon emissions respond to the changes in climate and atmospheric CO 2 concentration. , 13, 675-690, 2016 www.biogeosciences.net/13/675/2016/

Uncertainties and future perspective
As the version of ORCHIDEE used here does not include explicit forest stand structure and successional dynamics (age classes) within grid cells, we are unable to distinguish between the ecosystem effects of surface and crown fires. Instead, simulated fire effects (e.g., fuel combustion completeness, tree mortality) are applied to the whole grid cell in proportion to the burned fraction, as is done in most other fire models (Kloster et al., 2010;Li et al., 2012;Pfeiffer et al., 2013). Due to this inability to characterize the sub-grid level fire regime, fires seldom lead to the complete destruction of the whole forest stand and the re-establishment of a new cohort at the grid cell level (because the burned fraction seldom approaches unity). Instead, live biomass is removed in proportion to the simulated mortality multiplied by the simulated burned fraction. As forest is never completely killed, this approach may lead to a faster post-fire recovery in the model compared with that after a crown fire in reality. Our finding that the legacy sink peaked in the decade of 1990-1999 may be biased by this model behavior. Due to lack of explicit forest structure and vertical profile, the model is not able to simulate the thinning effects of surface fires. However, the evolution of fire impacts the simulated NBP with time, since disturbance on the regional scale ( Fig. 6) generally resembles the temporal pattern of post-fire forest NEP observed at site level (e.g., Fig. 1 in Amiro et al., 2010), that is, a carbon source effect at the time of and for a few years after fire occurrence, followed by long-term decaying sink effect. Besides the uncertainties introduced by the model's inability to distinguish crown fire versus surface fire, the underestimation of burned area in central Asian grasslands and eastern Siberian boreal forests is another source of uncertainty in our results. We expect the underestimation of grassland burned area to make little impact on the estimated fire legacy sink effects, as grasslands quickly recover from fires; thus, over a centennial timescale, their fire legacy impact on NBP would be close to zero. The underestimation of forest-fire-burned area in eastern Siberia, on the other hand, may lead to an underestimation of the fire legacy sink effect, as it is clear that crown fires create a long-term sink and surface fires also result in enhanced forest growth due to a short-term increase in available resources (Schulze et al., 2012).
However, it is difficult to quantify the uncertainties in our results by comparing them with observational data. For one thing, as forest age is not explicitly simulated within each grid cell, no forest age map could be derived from our model simulation; this precludes evaluating our results against inventory-based forest age maps. Despite the fact that a current-day forest age map has been compiled for boreal North America (Pan et al., 2011a;Stinson et al., 2011), such maps are still scarce for boreal Eurasia. Furthermore, the reconstruction of historical forest age dynamics will need a hindcast of the current forest age map by combining it with known disturbance histories. Geospatially explicit burned area data sets are available for Alaska, the USA and Canada, starting from the 1950s Stocks et al., 2003); those for Russia are only available starting satellitebased mapping of burned area (Giglio et al., 2013), and existing reconstructed data were based on simple assumptions and are subject to great uncertainties (Balshi et al., 2007;Mouillot and Field, 2005). To derive a better estimate of the role of fire in the boreal carbon cycle, a two-pronged approach is required: collecting historical fire data for the Eurasian boreal region and developing models further to include forest age groups in ORCHIDEE (Naudts et al., 2014).

Data availability
All data used in this manuscript will be made available upon request to the corresponding author (email address: chaoyue-joy@gmail.com).
The Supplement related to this article is available online at doi:10.5194/bg-13-675-2016-supplement.