Simulating boreal forest carbon dyna ics after stand-replacing fire disturbance : insights from a global process-based vegetation model

Introduction Conclusions References

due to increased fire and insect disturbances Stinson et al., 2011). Stand-replacing fires are the major natural disturbance in the North American boreal forest and thus fire frequency and severity, and their changes, play a key role in controlling large-scale boreal forest carbon dynamics (Balshi et al., 2007;Bond-Lamberty et al., 2007b;Harden et al., 2000;Hayes et al., 2011).
A long-term ecological consequence of stand-replacing fire is that trees are generally completely killed (although with a possible delay in time) and a large amount of snags are left to decompose (Amiro et al., 2006;Liu et al., 2005;Manies et al., 2005), while a secondary ecosystem succession is initiated, with a new forest cohort being established (Kasischke et al., 1995). Time-since-disturbance is an ecosystem state variable which is often important in determining the carbon flux and carbon stock (Litvak et al., 2003;Pregitzer and Euskirchen, 2004;Amiro et al., 2006Amiro et al., , 2010. Forest age is recognized as a key factor in explaining the sign and magnitude of the carbon balance of North American boreal forests (Kurz and Apps, 1999;Kurz et al., 2009;Stinson et al., 2011). Besides fluxes, fire and forest age may also impact the pattern of land-atmosphere water and energy exchange (Liu and Randerson, 2008;Liu et al., 2005). Thus, to simulate the effects of fire disturbance on the regional carbon balance and biophysical changes, these stand-age effects must be included in general purpose, large-scale vegetation models.
Since the 1970s, boreal forests have experienced a significant temperature rise and probably an intensifying fire regime (Euskirchen et al., 2007;Kasischke and Turetsky, 2006). A plethora of modeling studies have investigated how the historical carbon stock of boreal forest has changed in response to changes in climate, atmospheric CO 2 and the fire regime (Balshi et al., 2007;Bond-Lamberty et al., 2007b;Hayes et al., 2011;Yuan et al., 2012). There is general agreement that increasing atmospheric CO 2 tends to favor carbon accumulation in the system, while fires exert a negative effect on the net terrestrial carbon sink, but the role of temperature rise is uncertain. Satellite-derived evidence indicates that rising temperature generally increases the length of the thermal growing season; however, the carbon uptake resulting from the earlier growing season onset is partly offset by warming in the fall due to enhanced respiration (Barichivich et al., 2012;Piao et al., 2008). Despite all these efforts, the effects of climate variability and increased atmospheric CO 2 on postfire forest carbon flux trajectories have not been examined at the site level, using a model-observation combined approach within the context of postfire forest succession.
Here, we use a general process-based global vegetation model to simulate postfire forest carbon dynamics. After fire, the model establishes a new forest cohort with explicit stand structure, and simulates forest dynamics such as selfthinning and snag decomposition. Carbon flux and carbon stock data from 13 flux sites across North American boreal forests are used for comprehensive model calibration and er-ror characterization. The role of historical climate and atmospheric CO 2 trend in driving the postfire carbon flux trajectory is explored by model manipulation in combination with in situ measurements. Specifically, the objectives of this study are: -to perform a calibration of the model by adjusting its parameters simultaneously against CO 2 fluxes (gross primary production or GPP, net ecosystem production or NEP, total ecosystem respiration or TER), and on carbon stocks (total and aboveground biomass carbon stock, forest floor carbon stock, woody debris carbon, mineral soil carbon), and forest stand structure (basal area, stand density, mean diameter at breast height, or DBH); -to evaluate the model performance at different sites across different soil drainage conditions; -to attribute the effects of climate and atmospheric CO 2 trends in driving the postfire trajectory carbon fluxes at the sites examined; -to quantify the uncertainty in CO 2 fluxes and the residual model biases which cannot be reduced by calibration, in order to assess the modeled carbon balance uncertainty for regional-scale applications.

Site descriptions
Details of the three clusters of chronosequence sites distributed across the boreal forest biome in North America are given in Table 1. The sites include three flux towers (US-Bn1 to US-Bn3) in Alaska, USA, seven flux towers (CA-NS1 to CA-NS7) in western Manitoba, Canada, and another three flux towers (CA-SF1 to CA-SF3) in Saskatchewan, Canada. The climate conditions differ among the three clusters, with mean annual temperatures (MAT) of −2.1 • C for Alaska,−3.2 • C for Manitoba and 0.4 • C for Saskatchewan. Mean annual precipitation (MAP) is 290 mm for Alaska, 536 mm for Manitoba and 470 mm for Saskatchewan.
The year of the most recent fire event is available for all these sites. The eddy-covariance (EC) method was used to measure CO 2 fluxes for different times after fire disturbance (Table 1). For clarity, the fire events listed in Table 1 are hereinafter referred to as "the most recent fire event", and the periods during which EC observations have been made are referred to as "the EC observation period". The individual EC measurement locations with multiple towers (US-Bn1-3/CA-NS1-7/CA-SF1-3) will be referred to as "evaluation site(s)".
All the study sites are documented to have experienced stand-replacing fires, with all aboveground biomass being killed in the fire, resulting in complete forest regeneration  (Amiro et al., 2006;Goulden et al., 2006;Liu and Randerson, 2008). Vegetation before burning was dominated by black spruce (Picea mariana (Mill.) BSP) or jack pine (Pinus banksiana Lamb.), and was in various stages of forest recovery when measurements were collected. Dead tree boles were found to remain intact and upright until approximately 5 yr after fire, and most of these snags then fell within 10-15 yr after burning Liu and Randerson, 2008). Soils at the three evaluation site clusters are moderately well-drained to well-drained (Gower et al., 1997;Liu et al., 2005;Goulden et al., 2006). Permafrost occurs at the Alaska and Manitoba sites but is absent for the top 2 m of soil for all the Saskatchewan sites (Gower et al., 1997). The model was evaluated against various measurements collected from a range of published studies, or retrieved by the authors. Variables evaluated comprise GPP, NEP, TER, leaf area index (LAI), total (aboveground) biomass carbon, woody debris carbon, forest floor carbon, mineral soil carbon, forest DBH, individual density, and basal area. The detailed data source information is provided in Table S1.

Model description of ORCHIDEE-FM
For our research objectives, the global general purpose, process-based vegetation model ORCHIDEE (Krinner et al., 2005) is used. Brought to steady equilibrium state for vegetation and soil carbon pools after a long spinup, ORCHIDEE originally represents an "average" mature forest in a bigleaf approximation. To explicitly account for forest stand structure and processes related with forest age, a forest management module (FM) was developed (Bellassen et al., 2010). The version of ORCHDEE with FM module is hereinafter referred to as "ORCHIDEE-FM", and the original ORCHDEE (without FM) will be referred to as "the standard ORCHIDEE".
ORCHIDEE-FM includes several processes that are important in simulating forest stand development: (1) agerelated stand dynamics, including the decline in NPP in old forests, age limitation of LAI growth in young stands, and age-dependent allocation of woody NPP among stems and coarse roots; (2) a woody litter pool to account for the slow decomposition of that material; and (3) an explicit description of forest stand structure (stand density, tree diameter distribution, etc.) with the model able to simulate the selfthinning process, given a maximum density-biomass relationship for different types of forests.
The decomposition for above-and belowground litter and mineral soil organic carbon is modeled to follow a first-order kinetic equation in ORCHIDEE-FM (Parton et al., 1988). Acknowledging the different turnover rates of its components, the litter pool is sub-divided into metabolic, structural and woody litter. Plant live biomass is divided into eight compartments: leaf, above-and belowground sapwood, above-and belowground heartwood, root, fruit, and carbon reserve pool. In ORCHIDEE-FM, when carbon is transferred from the live biomass to the litter pools, the carbon from live biomass does not go directly to the three litter pools, it goes first to a temporary litter buffer, and is then allocated to the three litter pools according to prescribed ratios (Krinner et al., 2005).

8236
C. Yue et al.: Modeling boreal forest carbon dynamics after fire disturbance simulate boreal crown fire processes. A crown fire kills all the aboveground biomass and initiates a new forest cohort, and is similar to a clear-cut event from a successional modeling perspective, although the aboveground carbon is exported in the case of clear-cut. The "clear-cut" routine of the FM module was thus adapted to mimic fire burning effects and the establishment of new forest after fire. When crown fire occurs, part of the ground litter and tree biomass are removed from existing pools as carbon emissions into the atmosphere, and the unburned biomass is simulated as standing dead wood (snags), which gradually moves into the litter pool over time.

Age-related changes of LAI and productivity
Previous studies have reported that LAI increases with stand age in young boreal forests until ∼ 25 yr; after that LAI tends to saturate (Wang et al., 2003;Goulden et al., 2011). This age dependence of LAI is empirically modeled by scaling the maximum LAI (a parameter specified for each biome in ORCHIDEE which sets an achievable climax LAI, not necessarily reached at a given site) to increase with the square root of stand age until 25 yr. The model's default maximum LAI of 4.5 m 2 m −2 for the boreal forest biome was used.
To account for the productivity decrease in aging boreal forests (Bond-Lamberty et al., 2004;Mack et al., 2008;Goulden et al., 2011), a decrease in the optimal carboxylation rate in old-growth forest was added to the model. For stands aged between 100 and 200 years old, the optimal carboxylation rate is modeled to decrease by up to 10 % over the 100 yr in a linear way, to mimic the NPP reduction between the 74 and 154 year old forests in central Canada observed by Goulden et al. (2011). To promote agreement between the simulated and field-based estimates of productivity, the optimal carboxylation rate was adjusted from the default value of 35 µmol m −2 s −1 to 24 µmol m −2 s −1 .

Fire combustion fraction for different carbon pools
In fires of boreal North America, carbon emissions come mainly from ground fuel (or organic soil horizons) and live biomass combustion is limited to tree leaves and small branches. The combustion fraction depends on a range of factors and is not necessarily constant; however, the present study was not seeking an accurate simulation of the combustion process and thus simple fixed combustion fractions were adopted for different types of fuel. Considering the data in available field reports (Stocks, 1987(Stocks, , 1989Kasischke and Stocks, 2000;Kasischke and Hoy, 2012), the combustion fraction was set to 7 % for aboveground live biomass and 30 % for ground litter. The unburned biomass was transferred to the snag pool (see Supplement Sect. 1 for more detailed description).

Snag decay after a fire event
At the three study sites, fire-formed snags and downed deadwood represent a significant component of the forest carbon stocks. In reality, the decrease of snags (either standing or downed) occurs through three processes: decomposition by microbes, fragmentation and leaching. Field measurements show that the change of snag and woody debris as a function of time after fire follows a first-order kinetic equation (Bond-Lamberty and Gower, 2008): where dD dt is the change of woody debris, D is the existing woody debris and k is the annual decomposition constant. Bond-Lamberty and Gower (2008) reported a k value of 0.05-0.07 yr −1 for standing snag and downed woody debris using three independent methods: direct respiration measurement, repeated field surveys and chronosequence observation.
To represent this process, a snag pool was added in ORCHIDEE-FM-BF, which includes the aboveground unburned heartwood biomass and belowground thick pivotal woody roots. To simplify the model, it was assumed that no respiration occurs for standing snags and that the reduced snag mass enters the litter pool completely. This is consistent with the results of Manies et al. (2005), who found standing snags (stems) respire very slowly until they fall and come into contact with moss and soil.
During the first 20 yr after fire, we assumed an annual decay rate of 8 % (k = 0.08 yr −1 ) of snag being transferred to the litter pool, and during the next 20-40 yr an annual rate of 5 % (k = 0.05 yr −1 ). Forty years after fire, all the remaining snag becomes litter. These parameters are derived from the field measurements by Bond-Lamberty and Gower (2008).

Modifying litter and soil carbon turnover time
The default values of litter and soil carbon turnover time in ORCHIDEE were originally defined from the CENTURY soil carbon decomposition model (Parton et al., 1988), which was based on grassland ecosystems. We modified the residence time parameters in ORCHIDEE-FM-BF according to observations of radiocarbon content  and mass balance constraint studies (Harden et al., 2000;O'Donnell et al., 2011) for boreal soils. When driven by local climate, the modified mass weighted mean turnover time is ∼ 63 yr for aboveground litter (against 15.7 yr in the standard version of ORCHIDEE) and ∼ 2100 yr (against 1555 yr) for the mineral soil carbon pool.  Krinner et al., 2005). The PFT-specific parameters in ORCHIDEE-FM-BF are set to the model default values given by Krinner et al. (2005) except the modifications described in Sect. 2.3. Soil texture is prescribed at each site as shown in Table 1.

Simulation protocol
Our simulations are conducted in a way that does not require specific biometric measurement inputs (such as LAI, initial biomass, mineral soil carbon etc.). Rather, these variables are determined from climate forcing data and the model equations in a prognostic way, by starting with a spinup simulation followed with transient simulations. Boreal forests are known to experience recurrent fire disturbances, and similarly, our simulation protocol is designed to allow a gradual ecosystem carbon pool accumulation under recurrent fires until equilibrium state, before the most recent fire event is simulated at each site ( Fig. 1). As we did not know the exact history of fire disturbance for each evaluation site, the fire return interval (FRI) for all the sites was set uniformly to be 160 yr. This value was chosen for three reasons. First, the oldest stand in the evaluation sites is 155 yr after fire (CA-NS1 in Manitoba). Second, the FRI for Canadian central boreal forest was reported to range between 66 and 200 yr by Stocks et al. (2002). Third, Anderson et al. (2006) reported that during the past 4600 yr, the FRI for a lowland boreal forest in south central Alaska was 142 ± 70 yr. We acknowledge that assuming a uniform FRI of 160 yr is a simplification, as fire occurrence depends on several factors including climate (which varied during the Holocene) and the availability of fuel -thus not necessarily occurred at the equal intervals in the past.
According to the defined simulation protocol (Fig. 1), the model was initially run for a "first spinup" period starting from bare soil and without fires for 400 yr, followed by a "second spinup" of 3200 yr, i.e., 20 successive "fire rotations" with assumed stand-replacing fires occurring every 160 yr (each 160 yr period is called a "fire rotation"). This second spinup allows all forest ecosystem carbon pools (especially the mineral soil carbon pool) to reach a long-term equilibrium state in the presence of recurrent fire disturbance. Finally, the most recent fire event is simulated for the year of burning. The model is then driven with actual climate forcing data during the postfire period of forest regrowth, so that the model output can be compared to and evaluated against field measurements.
We define a "climate forcing history" which was used by most of the simulations. The monthly CMCD climate data were used in this climate forcing history. The rationale of the climate forcing history was to ensure that, for the historical spinup, the average state of historical climate was used, and for the period before and after the most recent fire event, actual observed climate data were used. Specifically, the climate forcing data used in different stages of the simulation protocol are as below (see also Table S2): a. For the first spinup and the first 19 fire rotations of the second spinup, we use the multi-year mean CMCD climate data of each site from the first year of meteorological station measurements until the year of the most recent fire (see Table S2 for details). This stable climate forcing was used with the goal of driving the model into an equilibrium with the average state of historical climate conditions. b. For the last (20th) fire rotation of the second spinup, at the beginning the average CMCD were used, but when entering the period of meteorological station measurement, the climate forcing shifted to actual year data. The time of shift depended on the year of the most recent fire and the duration of meteorological station observation at each site. By doing so, we made sure that before the most recent fire, actual monthly climate data were used to reflect the historical climate trends and variability that the ecosystem experienced at each site (see Table S2 for details).
c. For the postfire simulation after the most recent fire event, actual CMCD climate data continued to be used. But if the most recent fire occurred earlier than the start year of the meteorological station measurement, the average CMCD data were repeated until the year when meteorological station observations started and then real CMCD data began to be used. Two scenarios of atmospheric CO 2 concentration were used in the simulations, namely fixed CO 2 (CO 2 FIX) and variable CO 2 (CO 2 VAR). For the CO 2 FIX scenario, the atmospheric CO 2 concentration was fixed at the 1850 level (285 ppm) throughout the simulation. For the CO 2 VAR scenario, for the "first spinup" until the beginning of the 20th fire rotation of the second spinup, CO 2 concentration was fixed at the 1850 level, then during the 20th fire rotation, beginning from some point the atmospheric CO 2 was prescribed to increase (transient CO 2 concentration was used), and corresponded exactly to the CO 2 concentration at the year of burning for the most recent fire at each site. This was done to reflect the atmospheric CO 2 history experienced at each individual site. The difference between CO 2 VAR and CO 2 FIX allowed us to separate the effect of CO 2 fertilization on postfire CO 2 fluxes.

List of simulations in this study
Based on the defined simulation protocol, several simulations were performed; these are listed in Table 2 and described below: CNT-CMCD: control simulation. The simulation was performed using ORCHIDEE-FM-BF with all the modified features, processes, and parameters described in Sect. 2.3. The CO 2 VAR scenario was used.
GPPCAL-CMCD: GPP calibration simulation. This simulation used the same forcing as with the CNT-CMCD simulation, but the EC-observed multi-year mean GPP was assimilated (nudged) into the model. GPP assimilation was done by first calculating the ratio of average simulated to observed GPP at each study site for the EC observation period, and then applying this site-specific ratio throughout all the simulation stages (first spinup, second spinup and postfire simu-lation, see Fig. 1) to correct the simulated GPP for each run step. In this manner, the mean modeled GPP is tuned to be equal to the multi-year mean value observed by eddy covariance.
CNT-HHCD and GPPCAL-HHCD: these simulations are the same as CNT-CMCD and GPPCAL-CMCD, except that HHCD (half-hourly climate data) rather than CMCD were used during the EC observation period. Note that GPPCAL-CMCD and GPPCAL-HHCD simulations will have exactly the same GPP in every time step of the simulation (30 min) and the change from CMCD data to HHCD for the EC observation period will only affect respiration, and thus net ecosystem production.
CO 2 FIX-CLIMVAR: no CO 2 fertilization simulation. The CO 2 FIX scenario was used in the simulation, with varying climate data.
CO 2 FIX-CLIMFIX: no CO 2 fertilization simulation. The CO 2 FIX scenario was used in the simulation. But for the postfire simulation after the most recent fire, input climate forcing was fixed as the average monthly climate data that was used in the spinup runs. This simulation, together with CO 2 FIX-CLIMVAR simulation, was done only for sites CA-NS1, CA-SF1 and US-Bn1, in order to be compared with the GPPCAL-CMCD simulation to separate the roles of varying climate and CO 2 in the postfire carbon flux trajectory. To make the results comparable with the GPPCAL-CMCD simulation, the same site-specific GPP correction ratio used in GPPCAL-CMCD simulation was applied to correct modeled GPP in each step of the CO 2 FIX runs. This was done for two reasons. First, this site-specific ratio is considered to reflect the model internal structural error that could not be resolved by parameterization and is therefore independent of forcing factors. Second, there are no corresponding CO 2 FIX  1968-1980(avg) for CA-NS1, 1943-1976 for site CA-SF1, and1930-1959(avg) for site US-Bn1 (see also Table S2).
scenario observation data available to derive the site-specific ratios. ORC-STD: simulation with standard version of OR-CHIDEE. This simulation was done using the standard version of ORCHIDEE, with the same crown fire process as in Sect. 2.3 being implemented. No snag pool was represented and all unburned live biomass was sent to the litter pool (via the litter buffer) immediately after fire.
ORC-FM-NOSNAG: simulation with ORCHIDEE-FM with the same crown fire process implemented. Again, no fire caused snag process was represented in the model. This simulation and ORC-STD simulation were mainly for comparing with the CNT-CMCD simulation to demonstrate the "improvement chain" in simulating postfire forest regrowth by moving from standard ORCHIDEE to ORCHIDEE-FM and further to ORCHIDEE-FM-BF.

Method for simulation-measurement comparison
Due to the differences in the scope of the forest floor, woody debris and mineral soil carbon as measured in the field and as modeled, a scheme was developed to match the model output with field measurements for these variables (refer to Supplement Sect. 3 for more details).
Three metrics were used to quantitatively evaluate the simulation-measurement agreement in a comprehensive way: 1. Metric 1. Linear regression was used to examine the overall model-measurement agreement. Simulation data were regressed against measurements with the regression line forced through the origin (regression through origin; RTO) as it was assumed when the measurement value is zero, the simulation value should be zero as well. The RTO model used is where X i are measurement data, Y i are simulation data, and ε i is random error. If the regression slope was not significantly different from unity, then we considered it as fairly good agreement.
2. Metric 2. The root mean square deviation (RMSD) was used to quantify the model-measurement agreement in absolute terms. RMSD is the average quadratic distance between simulation and measurement values: where Y i and X i are simulated and measured data, respectively. We further calculated the systematic RMSD (RMSD_sys), which describes the error caused by systematic difference between simulation and measurement data, and unbiased RMSD (RMSD_unbias), which describes the error caused by internal variation among simulation values: whereŶ i is the predicted value by RTO regression, X i is measurement value.
whereŶ i is the predicted value by RTO regression, and Y i is simulation value. If RMSD is close to the field measurement error (instrument error, aggregation error, site-to-site and year-to-year precision in measurement) and RMSD is dominated by RMSD_unbias, then we considered that the modeling error was acceptable, and that the model realistically reproduced measurement data.
3. Metric 3. Measurement-simulation uncertainty overlap ratio was used to characterize measurement-simulation agreement with uncertainties in both being considered. First, we collected for each evaluation variable the measurement standard error or 90 % confidence interval, when they are reported in the source literature, and treated them as the measurement uncertainty. Then we calculated the number of data points where the simulation uncertainty and the measurement uncertainty overlapped. Finally we calculated the ratio of this overlapped number to the total number of measurement data points . We considered one overlap between model and measurement uncertainty to indicate that the model was able to simulate the measurement correctly. The model simulation uncertainty was constructed by pooling together simulation data for all the evaluation sites at the same site cluster, and the minimum-maximum range among model outputs was treated as simulation uncertainty.
Soil drainage condition is recognized as an important site-specific feature affecting forest ecosystem processes in North American boreal forests (Wang et al., 2003;Yi et al., 2009a) and may contribute to model-data misfit. Therefore, all the measurement data have been classified as either "dry" (with good soil drainage) or "wet" (with poor soil drainage) according to the information provided in the data source literature. As noted in Sect. 2.1, all the evaluation sites simulated could be considered as dry sites. Nevertheless, to have an idea of the model performance for the poordrainage sites as well, the model-measurement comparison statistics by metrics (1), (2), and (3) were calculated separately for "dry" and "wet" sites, and for all dry and wet sites combined.

Model improvement in ORCHIDEE-FM-BF and simulated fire carbon emissions
Moving from standard ORCHIDEE to ORCHIDEE-FM (or ORCHIDEE-FM-BF) significantly improved simulation results. GPP, total biomass carbon and heterotrophic respiration simulated by ORCHIDEE-FM-BF are in good agreement with the observation data (Refer to Supplement Sect. 4 for more detailed information). The simulated fire carbon emissions in the most recent fire events for the three site clusters for the GPPCAL-CMCD simulation are: for Saskatchewan 3.12 ± 1.3 kg C m −2 , for Manitoba 1.03 ± 0.17 kg C m −2 , and for Alaska 0.53 ± 0.03 kg C m −2 ( Table 3). The mean fire-carbon emissions among all the 13 study sites are 1.40 ± 1.15 kg C m −2 , with 0.33 kg C m −2 being emitted from crown burning and 1.06 kg C m −2 from surface burning. Linear regression (of form y = slope · x) is fitted between simulated and observed annual data across all evaluation sites. Sample size is 33 for GPP, and 31 for TER and NEP. Data reported here include regression slope, and the probability of the slope being significantly not different from 1 (p value), adjusted goodness of fit (adjusted R 2 , a modification of R 2 that adjusts for the number of explanatory variables), root mean square deviation (RMSD), systematic and unbiased RMSD (see Sect. 2.5.2 for detailed description). Tests with p value < 0.05 (i.e., slope = 1) are bold for emphasis.

Evaluation of simulated annual GPP, TER and NEP
Simulated annual GPP, TER and NEP is compared to EC measurements in Fig. 2. Model outputs are shown for simulations both before (CNT-CMCD and CNT-HHCD) and after (GPPCAL-CMCD and GPPCAL-HHCD) multi-year average GPP assimilation. Note, in the GPPCAL simulations, only the average multi-year GPP from the model is optimized to the measured mean value, so that the simulated inter-annual variability of GPP can still be compared to EC data. The CNT-CMCD simulation without GPP nudging, overestimates both GPP and TER across all study sites (Table 4, Fig. 2a) by approximately 30 %. No overestimation can be seen for NEP in CNT-CMCD simulation (Table 4, Fig. 2a), suggesting compensation of GPP and TER biases in NEP. In the CNT-HHCD simulation, the RTO regression slopes for GPP and TER are not significantly different from unity, indicating no systematic bias in simulation of these two carbon fluxes. Whereas NEP in CNT-HHCD is significantly lower than measurement values (∼ 85 % by RTO slope). This shows that the choice of climate forcing data (low vs. high frequency) strongly affects the model-data misfit.
For CNT-CMCD simulation, the simulated to observed multi-year average GPP ratios are different for each of the three clusters of evaluation sites. This ratio is 1.15 ± 0.27 for Saskatchewan sites, 1.42 ± 0.17 for Manitoba sites, and 2.13 ± 0.13 for Alaska sites. The overestimate of simulated GPP tends to be larger at higher latitudes. Because the GPPCAL-CMCD and GPPCAL-HHCD runs use the same nudged GPP, the model-measurement metrics for annual GPP are the same for these two simulations (Table 4). As expected, nudging multi-year average GPP greatly improves the model-measurement agreement for annual GPP, with RTO regression slope equal to unity, and RMSD being reduced by more than ∼ 70 % compared to the non-assimilated runs (Table 4). Nudging GPP simultaneously improves the TER simulation in both GPPCAL-CMCD and GPPCAL-HHCD simulations. Modelmeasurement agreement for NEP in GPPCAL-HHCD is also improved after assimilation (RTO regression slope changes from 0.13 to 0.85, RMSD is reduced by more than half; Table 4), but NEP remains underestimated (too small a modeled carbon sink) in GPPCAL-CMCD.

Evaluation of postfire evolution of LAI, biomass, forest floor, woody debris, and mineral soil organic carbon
In this section, the evaluations of ORCHIDEE-FM-BF output against biometric measurements are presented for LAI, biomass carbon, forest floor carbon, woody debris and mineral soil carbon ( few years, thus the simulation results from the GPPCAL-CMCD simulation are used for comparison.

Leaf area index
RTO regression slope between simulated and measured LAI indicates that the model underestimates LAI by an overall fraction of 24 % (RTO regression equal to 0.76, see Table 5) when all sites are considered together. The LAI in dry sites is underestimated by ∼ 30 % by the model and wet sites overestimated by 15 % (Table 5). The overlap ratio between simulated and measured data is 0.43, indicating 43 % of measurements are well simulated by the model. A close look at the model-and measurement data shows that LAI in the Manitoba dry sites have been underestimated by 50 % in the middle-aged forest (∼ 80 yr) and by 30 % in the oldaged forest (∼ 150 yr) (Fig. 3). The underestimation in the intermediate-aged forest is partly because the whole simulation pixel is prescribed to be fully covered by the boreal needleleaf trees, and this precludes the occurrence of broadleaf trees, which often dominate the early succession stage and have higher LAI and productivity.

Total and aboveground biomass
With the re-growing of the forest stand after fire, forest biomass carbon is modeled to increase continuously with forest age (Fig. 4). The RTO regression slope for all dry and wet sites combined is not significantly different from unity, indicating a very good overall model-measurement agreement. Yet for wet sites, modeled biomass carbon is overestimated by about 50 %, with a larger systematic RMSD than unbiased RMSD, likely because ORCHIDEE-FM-BF does not limit growth sufficiently at wet sites.

Forest floor carbon
RTO regression analysis indicates that forest floor carbon is underestimated by ∼ 50 %, when pooling all study sites together (Table 5), explaining the bigger systematic RMSD than unbiased RMSD. But due to the large uncertainty in the measured data (Fig. 5), the overlap ratio between simulated and measured data is 0.43, which is moderately good (43 % of the measurements are reproduced by the model).
If the three clusters of sites are examined separately, some additional biases are revealed. For the Manitoba sites, forest floor carbon is found to be underestimated by the model, mainly in very young forest (< 10 yr) and forest older than ∼ 70 yr. In contrast, for Saskatchewan sites the model overestimates forest floor carbon very strongly (by a factor of 3), Table 5. Model-measurements comparison metrics for LAI, biomass carbon, forest floor carbon, total woody debris, diameter at breast height (DBH), stand individual density, and basal area (BA) (see Table S1 for the definition of the variables; see Supplement Sect. 3 for how modeled and field observation data were matched against each other for woody debris and forest floor). See Sect. 2.5 for explanation for all items except N, which means the sample size. Cases with p value < 0.05 are bold for emphasis (which means regression slope is significantly different from 1). Bold RMSD_sys indicates that the value of RMSD_sys is bigger than RMSD_unbias, which means RMSD is dominated by systematic error and poor model-measurement agreement. Stars ( * ) indicate a better model-measurement agreement in dry sites than in wet sites. Woody debris includes all standing dead wood (STD), downed woody debris (DWD) and total woody debris (TWD).
although this result could well be biased by the scarcity of available measurements.

Woody debris
Model-measurement comparisons for woody debris, including standing dead wood (SDW), downed woody debris (DWD) and total woody debris (TWD), are shown in Fig. 6. The RTO regression slopes between simulated and measured woody debris for dry, wet and for all sites combined are 0.31, 0.30 and 0.50 respectively, suggesting underestimation by the model. The overlap ratio is 0.43, indicating moderate model-measurement agreement. For dry sites and all sites combined, the systematic RMSD is bigger than unbiased RMSD, indicating a persistent model bias. But for wet sites, this bias is reduced and becomes smaller than unbiased RMSD. At Manitoba sites, simulated SDW carbon is higher than field measurements for young forests (< 20 yr) (Fig. 6a). Note, for both downed woody debris and total woody debris, the measurements are extremely high in forests around 20 yr old. This is because tree mortality is delayed after fire (Bond-Lamberty and Gower, 2008), whereas in the model, unburned tree biomass is considered to become standing dead wood immediately after fire.

Mineral soil carbon and stand structure
The simulated mineral soil carbon is compared with the observation in a qualitative way. Summarizing, the simulated mineral soil carbon agrees best with the well-drained measurements, and is generally smaller than the poorly drained measurements. ORCHIDEE-FM-BF only has moderate to low capability at reproducing soil carbon. However the differences to observation are mainly found in the carbon pool with a small turnover rate, and thus the overall contribution to the flux error is relatively small (see Supplement Sect. 5 for more detailed discussion). The simulated forest stand individual density, mean DBH and basal area generally agree well with observation (see Supplement Sect. 6 for more detailed description).

Overall model performance across different soil drainage conditions
When examining model-measurement agreement for different soil drainage conditions, in terms of RTO regression goodness of fit (adjusted R 2 ), overlap ratio, and RMSD, the model is found to perform better at dry sites than at wet sites for LAI, biomass carbon stock, DBH, stand density and basal area (5 out of 6 variables for which the comparison is applicable). This indicates that the model has a generally better performance at dry sites than wet ones. Model-measurement agreement metrics across all evaluation sites for carbon fluxes and carbon stocks are summarized in Table 6, together with measurement accuracy. The average model-measurement overlap ratio for LAI, biomass, forest floor and total woody debris carbon was 42 %. Except for forest floor carbon, the model-measurement RMSDs for all variables in Table 6 are lower than or equal to the uncertainty of field measurement, indicating that the model performance is on average satisfactory, given the uncertainty in the observed data.

Attributing the role of past climate and CO 2 trends in the postfire evolution of carbon fluxes
To attribute the roles of the atmospheric CO 2 and varying climate in postfire forest CO 2 fluxes trajectory, CO 2 FIX-CLIMFIX and CO 2 FIX-CLIMVAR simulations were done for sites CA-NS1, CA-SF1, and US-Bn1 to be compared with GPPCAL-CMCD simulations. The postfire GPP, TER and NEP trajectory for Manitoba sites are shown in Fig. 7 (see Fig. S6 for sites in Saskatchewan and Alaska). The attribution analysis is done only for the period after the most recent fire events, the same period as the chronosequence study on each cluster of sites (for CA-NS1 155 yr, for CA-SF1 28 yr and for US-Bn1 83 yr). When attributing the effects of varying CO 2 and climate for the site CA-NS1 (Manitoba), we assume that the simulation result by CO 2 FIX-CLIMFIX for the time before 1968 (the starting year of meteorological station observations) followed the same curve as CO 2 FIX-CLIMVAR. This is mainly due to the restriction in the availability of historical climate data from the meteorological station. This assumption may lead to underestimated varying climate effects as the climate trend before 1968 is not taken into account. For clarity, we focus on the site CA-NS1 and then briefly discuss the sites CA-SF1 and US-Bn1.
We find that the temporal pattern and magnitudes of postfire CO 2 fluxes at the Manitoba sites (CA-NS1 to CA-NS7) over the past 150 yr are greatly driven by the fertilization effect of increasing atmospheric CO 2 . The EC-measured magnitudes of postfire GPP for different ages of forest after fire are much higher than the simulation result with fixed CO 2 (CO 2 FIX-CLIMVAR). The GPP measurements can only be reproduced by the model if the effect of increasing CO 2 is accounted for (GPPCAL-CMCD, Fig. 7a). The same also applies for TER and NEP, although the slight underestimation of NEP (underestimated carbon sink) by the GPPCAL-CMCD simulation is again shown (as shown in Fig. 2c). When comparing GPPCAL-CMCD and CO 2 FIX-CLIMVAR simulations, accounting for the effect of rising  Table 6. Model-measurement RMSD, RTO regression slope, overlap ratio for all sites combined, and measurement accuracy for GPP, TER, NEP, LAI, biomass carbon, total woody debris and forest floor carbon.
Climate effect CO 2 effect Combined effect Fig. 7. Simulated GPP (a), TER (b) and NEP (c) trajectory for the time of the 20th fire rotation (since the first vertical red dashed line) and for the chronosequence period (since the second vertical red dashed line) for Manitoba sites. Simulation results for the scenario of varying CO 2 with varying climate (GPPCAL-CMCD) are shown for all seven sites in colored lines, with the scenario of fixed CO 2 with varying climate (CO 2 FIX-CLIMVA) as a gray line, and the scenario of fixed CO 2 with fixed climate (CO 2 FIX-CLIMFIX) as black line for the site CA-NS1. Corresponding eddycovariance CO 2 flux measurements (Amiro et al., 2010;Goulden et al., 2011) at each site are shown as colored dots, with the colors corresponding to the colors of GPPCAL-simulation results for each site. For GPPCAL-CMCD simulation results, the numbers after the site names in the legend indicate the year of most recent fire event. The colored small arrows at the bottom of the panel (a) indicate the time to increase atmospheric CO 2 for each site in the GPPCAL-CMCD simulation. CO 2 leads to both an increase in GPP and TER, and a net increase in NEP (postfire carbon sink) ( Table 7).
All three fluxes of GPP, TER and NEP are higher for the CO 2 FIX-CLIMFIX simulation than CO 2 FIX-CLIMVAR, indicating that recent climate trends alone decrease the postfire carbon sequestration in our model simulation. GPP is decreased to a greater extent than TER when climate varies but CO 2 remains fixed, causing a net decrease in NEP (Table 7). In summary, in terms of mean annual NEP over the entire postfire period for the site CA-NS1 (155 yr), increasing CO 2 caused an increase of 30.5 g C m −2 yr −1 in mean annual NEP but climate trends alone a decrease by 3.5 g C m −2 yr −1 ; their combined effect is an increase of mean annual NEP by 26.9 g C m −2 yr −1 (Table 7). For the period after the most recent fire event, the GPPCAL-CMCD (varying CO 2 with varying climate) simulates a strong carbon sink in all the carbon stock compartments (total biomass, aboveground litter, belowground litter and mineral soil organic carbon) in the forest and accumulated NEP is 4.6 times higher than the carbon emissions in the fire (see Tables S3, S4 and S5 in the Supplement for the carbon budget for the three simulations at the three sites).
The role of varying CO 2 and climate in postfire carbon fluxes trajectory for the sites CA-SF1 and US-Bn1 is similar to that for CA-NS1. In terms of mean annual NEP over the period after the most recent fire event for CA-SF1 (28 yr), increasing CO 2 caused an increase in NEP by 54.1 g C m −2 yr −1 , while varying climate resulted in a decrease by 22.8 g C m −2 yr −1 ; their combined effect is an increase in mean annual NEP of 31.2 g C m −2 yr −1 . For US-Bn1 (83 yr), increasing CO 2 caused an increase of 7.0 g C m −2 yr −1 in mean annual NEP while varying climate a decrease by 7.5 g C m −2 yr −1 , leading to a decrease in mean annual NEP of 0.6 g C m −2 yr −1 .

Discussion
In the preceding sections, the ORCHIDEE-FM-BF model is evaluated against chronosequence measurements of carbon Table 7. Effects of varying atmospheric CO 2 and climate and their combined effect on mean annual carbon fluxes (g C m −2 yr −1 ) over the chronosequence study period after the most recent fire event. The respective postfire period length for the evaluation sites CA-NS1, CA-SF1 and US-Bn1 are 155, 28 and 83 yr.

Carbon flux
Effect on mean annual carbon flux Simulated mean annual carbon flux over the chronosequence period * over the chronosequence period 229.0 NEP −7.5 7.0 −0.6 5.6 −2.0 5.0 * Climate effect was calculated as difference between simulations of CO 2 FIX-CLIMVAR and CO 2 FIX-CLIMFIX, and CO 2 effect was calculated as difference between simulations of GPPCAL-CMCD and CO 2 FIX-CLIMVAR, and combined effect as difference between simulations of GPPCAL-CMCD and CO 2 FIX-CLIMFIX. fluxes, ecosystem carbon pools and stand structure, with focus on the capability of the model to reproduce the temporal evolution of each variable as a function of time-after-fire. To our knowledge, this is the first study that tries to evaluate a global process-based vegetation model for fire disturbances against multiple sites and multiple observation data sets. Although we believe the chronosequence stands are usually carefully selected to make them comparable and thus represent the process of forest development, there are still site-specific conditions which lead to spatial heterogeneity among sites and complicate model-measurement comparison. So the discussion below will focus on the more general issues on model-measurement comparison rather than explaining why some specific site or variable is overestimated or underestimated by the model.

Effect of nudging observed multi-annual average GPP on improving other CO 2 fluxes
As shown in Sect. 3.2, nudging observed GPP has improved the model-measurement agreement for carbon fluxes. Besides, the RMSDs for the carbon pool and stand structure variables have been reduced by 1-40 % in the GPPCAL-CMCD simulation compared with CNT-CMCD simulations (RMSDs for CNT-CMCD simulations not shown). The GPPCAL-HHCD simulation has better agreement with observation in terms of TER and NEP than the GPPCAL-CMCD simulation. This is because the 30 min climate fields used in the GPPCAL-HHCD simulation are assumed to be more realistic than the monthly ones in the GPPCAL-CMCD simulation. This forcing-related model output bias could not be completely solved by model calibra-tion, which has also been noticed by Zhao et al. (2012) and Lin et al. (2011).

Fire carbon emissions
The simulated fire carbon emissions are highest at Saskatchewan sites and lowest at Alaskan sites. As the fire combustion fraction is fixed in the model, the amount of simulated carbon emissions is directly related to the simulated forest floor carbon stock, which is the major source of the carbon emissions in boreal fires.
In a synthesis study of fire carbon emissions in North America, French et al. (2011) reported higher average carbon emissions in Alaska (2-5 kg C m −2 ) than in Canadian fires (1.2-1.9 kg C m −2 ), probably due to a higher level of organic soil carbon and deeper burning in Alaskan forests. Our simulated average fire carbon emissions for the Manitoba and Saskatchewan sites are 1.65 kg C m −2 and are therefore within the range given by French et al. (2011). However, the simulated emissions at Alaskan sites (0.53 kg C m −2 ) are far below their reported value, and are also lower than those reported by Randerson et al. (2006, 1.56 kg C m −2 for the site US-Bn3), and Kasischke and Hoy (2012, 3.01 kg C m −2 during large fire years and 1.69 kg C m −2 during small fire years).
The underestimation of fire carbon emissions at Alaskan sites could be partly explained by the underestimation in forest floor carbon stock. The forest floor carbon stock is underestimated by ∼ 60 % (simulated ∼ 1.3 kg C m −2 vs. observed ∼ 3.3 kg C m −2 ). Correspondingly, fire carbon emissions are underestimated by 65-80 % (simulated 0.5 kg C m −2 vs. reported 1.5-3 kg C m −2 from local studies). Besides, the simple combustion fraction scheme used in the model does not allow combustion fraction to vary with fire conditions, and this might also contribute to the errors in carbon emission estimation, as several studies point out that surface fuel combustion fraction contributes to the biggest uncertainty in estimating fire carbon emissions from boreal forest fires (French et al., 2004;de Groot et al., 2009).

Model performance across different soil drainage conditions
One characteristic of boreal forest ecosystems is that ecosystem processes are greatly modulated by soil drainage conditions (Wang et al., 2003;Bond-Lamberty et al., 2006). Well-drained stands occur on flat upland or south-facing slopes and are often not underlain by permafrost. Poorly drained stands occur on flat lowland, or north-facing slopes and are generally underlain by continuous or discontinuous permafrost (Harden et al., , 2001Wang et al., 2003;Turetsky et al., 2011). Stands with poor soil drainage are often found to be associated with open canopy forests with relatively poor tree growth and a low biomass, abundant bryophyte layer such as sphagnum (Sphagnum spp.) which typically grows in wet environments (Wang et al., 2003), frequently flooded soil, and a massive amount of organic soil carbon due to the slow decomposition in the anaerobic environment . ORCHIDEE-FM-BF is found to perform better in welldrained sites than in poorly drained ones. This is an expected result which could be explained by several key processes in the model. First, the soil hydrological processes in ORCHIDEE-FM-BF are simulated in a way that allows soil water to drain away as runoff or infiltration when excessive precipitation occurs (Ducoudré et al., 1993), this procedure does not allow soil flooding. In reality, soils on poorly drained sites (either underlain by permafrost or due to topographic effects) tend to be saturated with water. In addition, the thick surface organic layer acts to maintain moisture . These processes are however not included in ORCHIDEE-FM-BF. Second, in current versions of ORCHIDEE, the soil moisture always has a positive effect on photosynthesis (Krinner et al., 2005); the model fails to represent the detrimental effect of excessive soil water on plant roots and the negative effect on photosynthesis.
To improve model performance on poorly drained sites of a general process biogeochemical model such as OR-CHIDEE, the poor drainage related hydrological and ecophysiological processes need to be incorporated into the model. These processes may include, for example, frequent soil flooding, detrimental effects of excessive soil water on root function and photosynthesis (Kreuzwieser et al., 2004), and reduced soil organic matter decomposition and nutrient mineralization (Wickland and Neff, 2007) in case of excessive soil moisture. Despite some valuable attempts (Bond-Lamberty et al., 2007a;Pietsch et al., 2003;Yi et al., 2010), the explicit process-based representation of a forest with poor soil drainage or forested wetlands or peatland remains a big challenge in general process biogeochemical models.
Nevertheless, to examine the potential errors for regional application of ORCHIDEE-FM-BF on carbon fluxes and biomass carbon stocks, we tried to upscale the site-level simulation errors from both good and poor drainage conditions (Table 5) to regional scale, by using the soil drainage distribution information in both Alaska and Canada. The soil drainage map for Alaska by Harden et al. (2001) shows that ∼ 60 % of the soils were well-drained to moderately welldrained. For Canada, 65-75 % of the soils are well-drained or moderately well-drained (Soil Landscapes of Canada Working Group, 2010). To upscale the site-level error to regional scale in a rather simple way, the RTO regression slopes for dry and wet sites (Table 5) were used together with dry/wet soil distribution to derive an area-weighted error. Using this method, the model will probably generate an overestimation of total/aboveground biomass carbon stock by 12 % (Canada) to 18 % (Alaska), which is still within or comparable to the uncertainty of inventory-based net land-atmosphere carbon fluxes at national (Stinson et al., 2011) or regional scales (Hayes et al., 2012).
In summary, the model performance is found to be generally acceptable if all dry and wet sites are considered together. A process-based generic model like ORCHIDEE-FM-BF should not be fine-tuned at each study site for further reducing each error. Contrarily, only a multi-site agreement can be expected to assess for instance the model's capability to make regional simulations. Based on the results in Table 6, the key information is that the model-measurement error across multiple sites is comparable with the measurement accuracy, which justifies using the model for regional applications. We find a significantly positive role of increasing atmospheric CO 2 in explaining the postfire carbon uptake as observed by the chronosequence method. Other studies also report similar positive effects of historical CO 2 in increasing the boreal regional carbon stock (Balshi et al., 2007;Hayes et al., 2011). This is reasonable as the ambient CO 2 concentration is one of the factors that limit plant photosynthetic activity (Chapin III et al., 2002). Our results also reveal some more interesting aspects of the CO 2 fertilization effect on plant carbon uptake in the context of forest succession. As shown in Fig. 7, the magnitude of postfire GPP on each site, as measured by the chronosequence method, could only be reproduced at exactly the same simulated site, but not at others (for easy identification, the simulation data for the GPPCAL-CMCD simulation and measurement data for a single site are plotted in the same color). For example, when site CA-NS5 reaches 25 yr old, its annual GPP is much higher than site CA-NS1 at the time when it reached 25 yr old (∼ 700 vs. ∼ 450 g C m −2 yr −1 ). This is because the CO 2 concentration for 25-year-old CA-NS1 (289 ppm, in the year of 1875) is much lower than that for 25-year-old CA-NS5 (377 ppm, in the year of 2005). Similarly, the mean annual GPP for site CA-NS1 for 2003 is lower than the site CA-NS5 (615 vs. 741 g C m −2 yr −1 ) even though both forests share the same ambient CO 2 concentration, as the forest age at the latter site is much older (155 vs. 25 yr).

The role of past climate and CO
This finding might lead to the hypothesis that, for the same CO 2 concentration (or its increment), the positive fertilization effect is greater in young forests than in old ones. Because young forests have higher productivity than old forests (Chapin III et al., 2002;Goulden et al., 2011;Ryan et al., 1997), and therefore are likely to have a better capability of making use of the increased ambient CO 2 . If correct, this hypothesis implies that the current forest carbon sink should not be ascribed only to either forest succession or CO 2 increase, but must be due to a combination of both, i.e., the forest succession accompanied by concurrent atmospheric CO 2 increase. However, caution should be taken here because the increased nutrient availability also contributes to the carbon uptake in boreal forest (Magnani et al., 2007). As ORCHIDEE-FM-BF does not include nutrient dynamics, the CO 2 fertilization effect might be overestimated.

The role of varying climate
The results of previous studies of the overall effect of past climate variability on forest net carbon uptake are variable. Balshi et al. (2007) argue for a consistent positive climate effect (sink term) irrespective of the atmospheric CO 2 variation, whereas Hayes et al. (2011) find the climate effect shifts from positive to negative around 1970-1980s. Using factorial model simulation, Yuan et al. (2012) show that warming temperature alone would have increased the vegetation carbon stock but reduced the total soil carbon, leaving a moderate net sink term for the boreal forest in the Alaskan Yukon river basin, given that the atmospheric CO 2 is prescribed to increase. Our results also support a negative role of varying climate in the postfire CO 2 uptake. However, both photosynthetic carbon uptake and carbon release by respiration have been reduced, but with the former term being reduced more.
The negative effect of climate trends in reducing postfire CO 2 uptake might be due to increasing water stress, probably caused by the increase in temperature with the absence of sufficient increase in precipitation. Over the meteorological station observation period, mean annual temperature increased on each of the three site clusters, with 0.6 • C per decade in Alaska, 0.4 • C per decade in Manitoba and 0.3 • C per decade in Saskatchewan. The decadal trends in annual precipitation over the same period were 2.8 mm in Alaska, −14 mm in Manitoba and 7.7 mm in Saskatchewan.
For all the three regions, simulated soil moisture is lower in CO 2 FIX-CLIMVAR than CO 2 FIX-CLIMFIX (data not shown), indicating that with the climate trend, plants tend to suffer more water stress. Thus the negative effect on forest NEP of varying climate indicated by model simulation might be due to increasing drought, which has also been reported by, for example, Michaelian et al. (2011) and Ma et al. (2012).
Another important feature that is lacking in ORCHIDEE-FM-BF is the interactions among fire, organic soil dynamics, soil temperature and permafrost dynamics. In boreal forest, the organic soil layer thickness decreases after fire disturbance, leading to further increase of soil active layer depth and soil temperature and deepening of thawing front, this stimulates the decomposition of mineral soil carbon O'Donnell et al., 2011;Yi et al., 2009bYi et al., , 2010. The ORCHIDEE-FM-BF does not represent this varying soil organic layer thickness (or ground litter depth) nor its thermal and hydrological role during the cycles of fire disturbance. Thus the model fails to reproduce the seasonal amplitude of soil temperature change for different periods after fire (Fig. S7). The lack of these processes in the model might bias the effect of past temperature change on the carbon flux as presented here. However, the permafrost processes are included in other versions of ORCHIDEE (Koven et al., 2011) and future work will combine these different developments to take full account of the complex interactions discussed here.

The role of fire disturbance
Fires are generally recognized as a negative agent in the regional carbon balance (source term) due to the large carbon emissions released during the burning process (Balshi et al., 2007;Bond-Lamberty et al., 2007b;Hayes et al., 2011;Yuan et al., 2012). Besides, the soil active layer depth is also deepened after burning, leaving more soil carbon exposed to decomposition, as discussed above. However, fires also initiate new forest stands that will absorb carbon during their whole successional cycle, and the forest age effect is an important factor in explaining the contemporary carbon sink in the boreal forests (Kurz and Apps, 1999;Pan et al., 2011;Stinson et al., 2011). Thus the ultimate role of fire in the regional carbon balance will depend on the complex interactions among climate, atmospheric CO 2 increase, nutrient availability (which is not included in ORCHIDEE-FM-BF), permafrost and soil dynamics, and forest succession.

Summary and conclusion
In this study, we adapted a general purpose process-based global vegetation model ORCHIDEE to ORCHIDEE-FM-BF, to simulate stand-replacing fires in boreal forests and investigate postfire carbon dynamics during forest regrowth. The carbon dynamics in relation to stand age are explicitly simulated by initiating a new forest cohort with stand structure, with self-thinning and postfire snag decomposition also being simulated. By using a simulation protocol with recurrent crown fire processes in a completely prognostic way, both realistic fire carbon emissions and postfire carbon fluxes could be reproduced by the model.
Our results also highlight the role of forest succession and atmospheric CO 2 fertilization in explaining the current postfire forest carbon sink as generally observed by the chronosequence method. We find that this stand-age-dependent sink could not simply be ascribed to either the forest succession or the CO 2 fertilization but must be due to a combination of both, i.e., the successional sink has been amplified by the CO 2 fertilization.
Despite lack of some important features in boreal forest ecosystems (which may include poor drainage processes, nutrient dynamics, permafrost dynamics, organic soil horizon and succession dynamics etc.), the model is generally able to reproduce postfire forest carbon dynamics when errors are upscaled across multiple sites and different soil drainage conditions. The model calibration in this study will allow the regional carbon balance analysis for boreal forest by using the approach that, the demographic effects of standing-replacing fires are accounted for with a mosaic of different aged forest cohorts being established and simulated. And this progress will help to more accurately quantify the contribution of fires to the historical and current carbon balance of boreal forest.