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

C. Yue, Philippe Ciais, S. Luyssaert, P. Cadule, J. Harden, J. Randerson, Valentin Bellassen, T. Wang, S. L. Piao, B. Poulter, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

C. Yue, Philippe Ciais, S. Luyssaert, P. Cadule, J. Harden, et al. Simulating boreal forest carbon dynamics after stand-replacing fire disturbance: insights from a global process-based vegetation model: Simulating boreal forest carbon dynamics after stand-replacing fire disturbance. Biogeosciences, 2013, 10 (12), pp.8233. 10.5194/bg-10-8233-2013. [https://hal.science/hal-01190133](https://hal.science/hal-01190133)

HAL Id: hal-01190133

https://hal.science/hal-01190133

Submitted on 28 May 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

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

C. Yue<fnref n="1" />, P. Ciais<fnref n="1" />, S. Luyssaert<fnref n="1" />, P. Cadule<fnref n="1" />, J. Harden<fnref n="2" />, J. Randerson<fnref n="3" />, V. Bellassen<fnref n="4" />, T. Wang<fnref n="1" />, S. L. Piao<fnref n="5" />, B. Poulter<fnref n="1" />, and N. Viovy<fnref n="1" />

*Correspondence to:* C. Yue (chao.yue@lsce.ipsl.fr)

Received: 12 March 2013 – Published in Biogeosciences Discuss.: 24 April 2013

Revised: 30 October 2013 – Accepted: 9 November 2013 – Published: 13 December 2013

**Abstract.** Stand-replacing fires are the dominant fire type in North American boreal forests. They leave a historical legacy of a mosaic landscape of different aged forest cohorts. This forest age dynamics must be included in vegetation models to accurately quantify the role of fire in the historical and current regional forest carbon balance. The present study adapted the global process-based vegetation model ORCHIDEE to simulate the CO₂ emissions from boreal forest fire and the subsequent recovery after a stand-replacing fire; the model represents postfire new cohort establishment, forest stand structure and the self-thinning process. Simulation results are evaluated against observations of three clusters of postfire forest chronosequences in Canada and Alaska. The variables evaluated include: fire carbon emissions, CO₂ fluxes (gross primary production, total ecosystem respiration and net ecosystem exchange), leaf area index, and biometric measurements (aboveground biomass carbon, forest floor carbon, woody debris carbon, stand individual density, stand basal area, and mean diameter at breast height). When forced by local climate and the atmospheric CO₂ history at each chronosequence site, the model simulations generally match the observed CO₂ fluxes and carbon stock data well, with model-measurement mean square root of deviation comparable with the measurement accuracy (for CO₂ flux ~100 g C m⁻<fnref n="2" /> yr⁻<fnref n="1" />, for biomass carbon ~1000 g C m⁻<fnref n="2" /> and for soil carbon ~2000 g C m⁻<fnref n="2" />). We find that the current postfire forest carbon sink at the evaluation sites, as observed

by chronosequence methods, is mainly due to a combination of historical CO₂ increase and forest succession. Climate change and variability during this period offsets some of these expected carbon gains. The negative impacts of climate were a likely consequence of increasing water stress caused by significant temperature increases that were not matched by concurrent increases in precipitation. Our simulation results demonstrate that a global vegetation model such as ORCHIDEE is able to capture the essential ecosystem processes in fire-disturbed boreal forests and produces satisfactory results in terms of both carbon fluxes and carbon-stock evolution after fire. This makes the model suitable for regional simulations in boreal regions where fire regimes play a key role in the ecosystem carbon balance.

## 1 Introduction

Boreal forests store a large amount of the global terrestrial-ecosystem carbon, with 78 Pg C in biomass and ~230 Pg C in soil (Kasischke, 2000). This forest biome has been estimated to be a carbon sink over the late 20th century (Kurz and Apps, 1999; McGuire et al., 2009; Pan et al., 2011), but the carbon stock and carbon sink are highly sensitive to fire disturbance (Balshi et al., 2009; Bond-Lamberty et al., 2007b). Particularly, forests in boreal North America appeared likely to become a carbon source in the past decade

<footnote n="1">¹Laboratoire des Sciences du Climat et de l'Environnement, UMR8212, LSCE CEA CNRS UVSQ, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">²US Geological Survey, Menlo Park, CA, USA</footnote>

<footnote n="3">³Department of Earth System Science, University of California, Irvine, California, USA</footnote>

<footnote n="4">⁴CDC Climat, 47 rue de la Victoire, 75009 Paris, France</footnote>

<footnote n="5">⁵Department of Ecology, College of Urban and Environmental Science, Peking University, Beijing 100871, China</footnote>
due to increased fire and insect disturbances (Hayes et al., 2011; 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., 2006, 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₂ 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₂ 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₂ 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 self-thinning 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 error characterization.

The role of historical climate and atmospheric CO₂ 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₂ 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₂ trends in driving the postfire trajectory carbon fluxes at the sites examined;

- to quantify the uncertainty in CO₂ 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.

# 2 Materials and methods

## 2.1 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₂ 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
<table><caption>**Table 1.** Measurement sites used in this study for model evaluation, their geographical coordinates, soil texture, pre-fire dominant vegetation species, year of the most recent fire event, and the period of eddy-covariance observation.</caption><thead><tr><th rowspan="2">Site Name</th><th rowspan="2"></th><th rowspan="2">Lat</th><th rowspan="2">Lon</th><th colspan="3">Soil texture<sup>a</sup></th><th rowspan="2">Pre-fire dominant species<sup>b</sup></th><th rowspan="2">Year of burn</th><th rowspan="2">OSY<sup>c</sup></th><th rowspan="2">OEY<sup>c</sup></th></tr><tr><th>Sand</th><th>Silt</th><th>Clay</th></tr></thead><tbody><tr><td rowspan="3">Saskatchewan</td><td>CA-SF1</td><td>54.5</td><td>-105.8</td><td>0.58</td><td>0.32</td><td>0.1</td><td>Jack pine</td><td>1977</td><td>2003</td><td>2005</td></tr><tr><td>CA-SF2</td><td>54.3</td><td>-105.9</td><td>0.58</td><td>0.32</td><td>0.1</td><td>Jack pine</td><td>1989</td><td>2003</td><td>2005</td></tr><tr><td>CA-SF3</td><td>54.1</td><td>-106.0</td><td>0.58</td><td>0.32</td><td>0.1</td><td>Jack pine</td><td>1998</td><td>2003</td><td>2005</td></tr><tr><td rowspan="8">Manitoba</td><td>CA-NS1</td><td>55.9</td><td>-98.5</td><td>0.02</td><td>0.13</td><td>0.86</td><td>Black spruce</td><td>1850</td><td>2002</td><td>2005</td></tr><tr><td>CA-NS2</td><td>55.9</td><td>-98.5</td><td>0.27</td><td>0.31</td><td>0.42</td><td>Black spruce</td><td>1930</td><td>2001</td><td>2005</td></tr><tr><td>CA-NS3</td><td>55.9</td><td>-98.4</td><td>0.27</td><td>0.31</td><td>0.42</td><td>Black spruce</td><td>1964</td><td>2001</td><td>2005</td></tr><tr><td>CA-NS4</td><td>55.9</td><td>-98.4</td><td>0.27</td><td>0.31</td><td>0.42</td><td>Black spruce</td><td>1964</td><td>2002</td><td>2004</td></tr><tr><td>CA-NS5</td><td>55.9</td><td>-98.5</td><td>0.27</td><td>0.31</td><td>0.42</td><td>Black spruce</td><td>1981</td><td>2001</td><td>2005</td></tr><tr><td>CA-NS6</td><td>55.9</td><td>-99.0</td><td>0.27</td><td>0.31</td><td>0.42</td><td>Black spruce</td><td>1989</td><td>2001</td><td>2005</td></tr><tr><td>CA-NS7</td><td>56.6</td><td>-100.0</td><td>0.34</td><td>0.29</td><td>0.37</td><td>Black spruce</td><td>1998</td><td>2002</td><td>2005</td></tr><tr><td>US-Bn1</td><td>63.9</td><td>-145.4</td><td>0.82</td><td>0.12</td><td>0.06</td><td>Black spruce</td><td>1920</td><td>2003</td><td>2003</td></tr><tr><td rowspan="2">Alaska</td><td>US-Bn2</td><td>63.9</td><td>-145.4</td><td>0.82</td><td>0.12</td><td>0.06</td><td>Black spruce</td><td>1987</td><td>2003</td><td>2003</td></tr><tr><td>US-Bn3</td><td>63.9</td><td>-145.7</td><td>0.82</td><td>0.12</td><td>0.06</td><td>Black spruce</td><td>1999</td><td>2003</td><td>2003</td></tr></tbody></table>

(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 (Goulden et al., 2006; 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.

## 2.2 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 big-leaf 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 ORCHIDEE with FM module is hereinafter referred to as “ORCHIDEE-FM”, and the original

ORCHIDEE (without FM) will be referred to as “the standard ORCHIDEE”.

ORCHIDEE-FM includes several processes that are important in simulating forest stand development: (1) age-related 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 self-thinning 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).

## 2.3 Modifications of ORCHIDEE-FM into ORCHIDEE-FM-BF for boreal fires

ORCHIDEE-FM was adapted into ORCHIDEE-FM-BF (ORCHIDEE-FM for Boreal Fires) to make it able to

<footnote><sup>a</sup>Soil texture information, where available, is provided by site PIs, otherwise it is completed from the soil map of Zobler (1986) and translated into sand/silt/clay fractions by the model default values that correspond to soil types in the Zobler map. <sup>b</sup>Black spruce: *Picea mariana* (Mill.) BSP; Jack pine: *Pinus banksiana* Lamb. <sup>c</sup>OSY, Eddy covariance observation period start year; OEY, Eddy covariance observation period end year.</footnote>
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.

### 2.3.1 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² m⁻² 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⁻² s⁻¹ to 24 µmol m⁻² s⁻¹.

### 2.3.2 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, 1989; Kasischke 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).

### 2.3.3 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):

$$ \frac{dD}{dt} = -kD, \qquad (1) $$

where $\frac{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⁻¹ 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 \text{ 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 \text{ 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).

### 2.3.4 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 (Trumbore and Harden, 1997) 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.
## 2.4 Simulation set-up

### 2.4.1 Climate forcing data

Two sets of climate data are used for climate forcing. Composite Monthly Climate Data (hereafter referred to as CMCD) include monthly temperature and precipitation from meteorological stations close to the evaluation sites with the missing fields being filled by CRU3.1 data. Half-hourly climate data (hereafter referred to as HHCD) include gap-filled in situ meteorological measurements from the eddy-covariance sites (see Supplement Sect. 2 for more detailed description).

### 2.4.2 Vegetation, soil and other input data

ORCHIDEE-FM-BF is built on the concept of plant functional types (PFTs), with each PFT being associated with a specific vector of model parameters. Boreal evergreen needleleaf forests (black spruce or jack pine) are the climax vegetation type on the evaluation sites, although broadleaf trees such as trembling aspen often dominate at the early succession stage. The dynamic vegetation mechanisms in ORCHIDEE-FM-BF do not allow the realistic representation of this species shift at the intermediate forest stage on a single simulation pixel, thus each site is prescribed to be fully covered by the boreal needleleaf evergreen forest PFT (see 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.

### 2.4.3 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 \pm 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.
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Illustration of the simulation protocol, climate forcing and atmospheric CO₂ forcing data for various simulations. First spinup: the model is first run for 400 yr starting from bare soil without fire; second spinup: after first spinup, the model is run for 3200 yr, which consists of 20 successive “fire rotations” with assumed crown fires occurring every 160 yr (i.e., a FRI of 160 yr). Postfire simulation: the most recent fire event is simulated during the occurrence year, and the model is driven with actual observed climate forcing data for postfire regrowth. For clarity, CNT-HHCD, GPPCAL-HHCD and CO₂FIX-CLIMFIX simulations were not shown. Refer to Sect. 2.4.4 and Table 2 for more detailed descriptions of different simulations.</figcaption></figure>

Two scenarios of atmospheric CO₂ concentration were used in the simulations, namely fixed CO₂ (CO₂FIX) and variable CO₂ (CO₂VAR). For the CO₂FIX scenario, the atmospheric CO₂ concentration was fixed at the 1850 level (285 ppm) throughout the simulation. For the CO₂VAR scenario, for the “first spinup” until the beginning of the 20th fire rotation of the second spinup, CO₂ concentration was fixed at the 1850 level, then during the 20th fire rotation, beginning from some point the atmospheric CO₂ was prescribed to increase (transient CO₂ concentration was used), and corresponded exactly to the CO₂ concentration at the year of burning for the most recent fire at each site. This was done to reflect the atmospheric CO₂ history experienced at each individual site. The difference between CO₂VAR and CO₂FIX allowed us to separate the effect of CO₂ fertilization on postfire CO₂ fluxes.

### 2.4.4 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₂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 simulation,

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₂FIX-CLIMVAR:** no CO₂ fertilization simulation. The CO₂FIX scenario was used in the simulation, with varying climate data.

**CO₂FIX-CLIMFIX:** no CO₂ fertilization simulation. The CO₂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₂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₂ 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₂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₂FIX
<table><caption>**Table 2.** List of simulations, and their climate and atmospheric CO₂ forcing.</caption><thead><tr><th>Name</th><th>Model used</th><th>Climate forcing</th><th>CO₂ forcing<sup>a</sup></th></tr></thead><tbody><tr><td>CNT-CMCD</td><td>ORCHIDEE-FM-BF</td><td>CMCD</td><td>CO₂VAR</td></tr><tr><td>CNT-HHCD</td><td>ORCHIDEE-FM-BF</td><td>HHCD for only EC obs. period</td><td>CO₂VAR</td></tr><tr><td>GPPCAL-CMCD</td><td>ORCHIDEE-FM-BF</td><td>CMCD</td><td>–</td></tr><tr><td>GPPCAL-HHCD</td><td>ORCHIDEE-FM-BF</td><td>HHCD for only EC obs. period</td><td>–</td></tr><tr><td>CO₂FIX-CLIMVAR (only for sites CA-NS1, CA-SF1, US-Bn1)</td><td>ORCHIDEE-FM-BF</td><td>CMCD</td><td>CO₂FIX</td></tr><tr><td>CO₂FIX-CLIMFIX (only for sites CA-NS1, CA-SF1, US-Bn1)</td><td>ORCHIDEE-FM-BF</td><td>The same average CMCD as used for spinup runs<sup>b</sup></td><td>CO₂FIX</td></tr><tr><td>ORC-STD</td><td>Standard ORCHIDEE with fire and without snag</td><td>CMCD</td><td>CO₂VAR</td></tr><tr><td>ORC-FM-NOSNAG</td><td>ORCHIDEE-FM with fire but without snag</td><td>CMCD</td><td>CO₂VAR</td></tr></tbody></table>

scenario observation data available to derive the site-specific ratios.

ORC-STD: simulation with standard version of ORCHIDEE. 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.

## 2.5 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

$$Y_i = \text{slope} \cdot X_i + \varepsilon_i, \qquad (2)$$

where $X_i$ are measurement data, $Y_i$ are simulation data, and $\varepsilon_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:

$$RMSD = \sqrt{\frac{\sum (X_i - Y_i)^2}{n}}, \qquad (3)$$

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:

$$RMSD\_sys = \sqrt{\frac{\sum (X_i - \hat{Y}_i)^2}{n}}, \qquad (4)$$

where $\hat{Y}_i$ is the predicted value by RTO regression, $X_i$ is measurement value.

$$RMSD\_unbias = \sqrt{\frac{\sum (Y_i - \hat{Y}_i)^2}{n}}, \qquad (5)$$

where $\hat{Y}_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

<footnote><sup>a</sup>Atmospheric CO₂ forcing for GPPCAL-CMCD and GPPCAL-HHCD simulations is not applicable because GPP is externally forced to the model in the simulation, by nudging mean multi-year observed annual GPP into the model. <sup>b</sup>For the CO₂FIX-CLIMFIX simulation, the average monthly climate data used are: 1968–1980(avg) for CA-NS1, 1943–1976(avg) for site CA-SF1, and 1930–1959(avg) for site US-Bn1 (see also Table S2).</footnote>
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 (Bond-Lamberty et al., 2006). 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 poor-drainage 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.

# 3 Results

## 3.1 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

<table><caption>**Table 3.** Simulated fire carbon emissions (kg C m⁻²) from live biomass, aboveground litter, and total carbon emissions for the most recent fire event for the three site clusters by GPPCAL-CMCD simulation. The mean values for each site cluster are shown in the table with standard deviation being shown within brackets.</caption><thead><tr><th>Site cluster</th><th>Live biomass</th><th>Aboveground litter</th><th>Total carbon emissions</th></tr></thead><tbody><tr><td>Saskatchewan</td><td>0.60 (0.15)</td><td>2.52 (1.22)</td><td>3.12 (1.35)</td></tr><tr><td>Manitoba</td><td>0.31 (0.06)</td><td>0.72 (0.12)</td><td>1.03 (0.17)</td></tr><tr><td>Alaska</td><td>0.12 (0.02)</td><td>0.41 (0.02)</td><td>0.53 (0.03)</td></tr><tr><td>Average</td><td>0.33 (0.18)</td><td>1.06 (0.98)</td><td>1.40 (1.15)</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Simulated vs. observed annual gross primary production (sky blue), total ecosystem respiration (orange) and net ecosystem production (green), before (a, c) and after (b, d) nudging observed mean multi-year GPP. Observed annual carbon fluxes are from Amiro et al. (2010). A 1:1 line is shown as the dashed gray line. Colored dashed lines indicate RTO regression lines. Distinctions are made among Manitoba (circles), Saskatchewan (cross symbol, “+”) and Alaska (stars) data.</figcaption></figure>

simulation are: for Saskatchewan 3.12 ± 1.3 kg C m⁻², for Manitoba 1.03 ± 0.17 kg C m⁻², and for Alaska 0.53 ± 0.03 kg C m⁻² (Table 3). The mean fire-carbon emissions among all the 13 study sites are 1.40 ± 1.15 kg C m⁻², with 0.33 kg C m⁻² being emitted from crown burning and 1.06 kg C m⁻² from surface burning.
<table><caption>**Table 4.** Model-measurement comparison metrics for annual GPP, TER and NEP. The left-hand (right-hand) column shows values before (after) nudging the observed multi-year GPP into the model (see Sect. 2.4.4 for more details).</caption><thead><tr><th rowspan="2">Metrics considered</th><th colspan="3">CNT-CMCD</th><th colspan="3">GPPCAL-CMCD</th></tr><tr><th>GPP</th><th>TER</th><th>NEP</th><th>GPP</th><th>TER</th><th>NEP</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">CNT-CMCD</td></tr><tr><td>Slope</td><td>1.28</td><td>1.29</td><td>0.79</td><td>1.00</td><td>1.05</td><td>0.41</td></tr><tr><td>Adjusted R<sup>2</sup></td><td>0.95</td><td>0.96</td><td>0.62</td><td>0.98</td><td>0.98</td><td>0.48</td></tr><tr><td><i>p</i> value</td><td><b>0.00</b></td><td><b>0.00</b></td><td>0.06</td><td>1.00</td><td>0.05</td><td><b>0.00</b></td></tr><tr><td>RMSD</td><td>247</td><td>215</td><td>61</td><td>76</td><td>87</td><td>70</td></tr><tr><td>RMSD_sys</td><td>175</td><td>160</td><td>21</td><td>0</td><td>31</td><td>57</td></tr><tr><td>RMSD_unbias</td><td>174</td><td>143</td><td>57</td><td>76</td><td>81</td><td>39</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">CNT-HHCD</td></tr><tr><td>Slope</td><td>0.91</td><td>1.02</td><td>0.13</td><td>1.00</td><td>0.98</td><td>0.85</td></tr><tr><td>Adjusted R<sup>2</sup></td><td>0.84</td><td>0.92</td><td>-0.02</td><td>0.98</td><td>0.98</td><td>0.58</td></tr><tr><td><i>p</i> value</td><td>0.19</td><td>0.76</td><td><b>0.00</b></td><td>1.00</td><td>0.31</td><td>0.24</td></tr><tr><td>RMSD</td><td>244</td><td>166</td><td>138</td><td>76</td><td>76</td><td>69</td></tr><tr><td>RMSD_sys</td><td>57</td><td>9</td><td>84</td><td>0</td><td>14</td><td>15</td></tr><tr><td>RMSD_unbias</td><td>237</td><td>166</td><td>110</td><td>76</td><td>74</td><td>67</td></tr></tbody></table>

Linear regression (of form $y = slope \cdot 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.

## 3.2 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.

<figure><img src="image_7.png" /></figure>

**Fig. 3.** Simulated and measured leaf area index (LAI) as a function of time after fire. Model results (Manitoba: sky blue, Saskatchewan: green, Alaska: orange) are presented by pooling together outputs for all evaluation sites of the same cluster, with the solid line indicating the mean value, and shaded area showing between-site minimum-maximum range. Measurements from different sources are shown separately for Manitoba (circles), Saskatchewan (diamonds) and Alaska (triangles), with wet (dry) site measurements as filled (open) signs. Error bars on the measurement points indicate 90% confidence interval measurement uncertainty. The inset panel shows the overall model-measurement agreement along a 1:1 line for dry (small open circles) and wet (small cross symbol, "+") site measurements separately.

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. Model-measurement 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.

## 3.3 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 (Table 5). Note that the two GPP calibration simulations only differ for the EC observation period and carbon stocks are not expected to change greatly during these
<figure><img src="image_8.png" /></figure>

**Fig. 4.** Simulated vs. observed total biomass carbon (for Manitoba) and aboveground biomass carbon (for Saskatchewan and Alaska) as a function of time after fire. Model results (Manitoba: sky blue, Saskatchewan: green, Alaska: orange) are presented by pooling together outputs for all evaluation sites of the same cluster, with the solid line indicating the mean value, and shaded area showing between-site minimum–maximum range. Measurements from different sources are shown separately for Manitoba (circles), Saskatchewan (diamonds) and Alaska (triangles), with wet (dry) site measurements as filled (open) signs. Error bars on the measurement points indicate 90% confidence interval measurement uncertainty. The inset panel shows the overall model-measurement agreement along a 1:1 line for dry (small open circles) and wet (small cross symbol, “+”) site measurements separately.

few years, thus the simulation results from the GPPCAL-CMCD simulation are used for comparison.

### 3.3.1 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 old-aged 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 broadleaf trees, which often dominate the early succession stage and have higher LAI and productivity.

<figure><img src="image_9.png" /></figure>

**Fig. 5.** Simulated vs. observed forest floor carbon stock as a function of time after fire. Model results (Manitoba: sky blue, Saskatchewan: green, Alaska: orange) are presented by pooling together outputs for all evaluation sites of the same cluster, with the solid line indicating the mean value, and shaded area showing between-site minimum–maximum range. Measurements from different sources are shown separately for Manitoba (circles), Saskatchewan (diamonds) and Alaska (triangles), with wet (dry) site measurements as filled (open) signs. Error bars on the measurement points indicate 90% confidence interval measurement uncertainty. The inset panel shows the overall model-measurement agreement along a 1:1 line for dry (small open circles) and wet (small cross symbol, “+”) site measurements separately.

### 3.3.2 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.

### 3.3.3 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><caption>**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).</caption><thead><tr><th>Items</th><th>Drainage</th><th>LAI</th><th>Biomass carbon</th><th>Forest floor carbon</th><th>Woody debris</th><th>DBH</th><th>Stand individual density</th><th>Basal Area</th></tr></thead><tbody><tr><td rowspan="3">N</td><td>All</td><td>33</td><td>37</td><td>35</td><td>60</td><td>12</td><td>11</td><td>12</td></tr><tr><td>Dry</td><td>26</td><td>27</td><td>26</td><td>53</td><td>7</td><td>6</td><td>7</td></tr><tr><td>Wet</td><td>7</td><td>10</td><td>9</td><td>7</td><td>5</td><td>5</td><td>5</td></tr><tr><td rowspan="3">Slope</td><td>All</td><td>0.76</td><td>1.03</td><td>0.51</td><td>0.31</td><td>0.79</td><td>1.14</td><td>0.89</td></tr><tr><td>Dry</td><td>0.71</td><td>0.94</td><td>0.54</td><td>0.30</td><td>0.74</td><td>1.06</td><td>0.76</td></tr><tr><td>Wet</td><td>1.15</td><td>1.55</td><td>0.47</td><td>0.50</td><td>0.89</td><td>1.32</td><td>1.50</td></tr><tr><td rowspan="3">p value</td><td>All</td><td><strong>0.00</strong></td><td>0.69</td><td><strong>0.00</strong></td><td><strong>0.00</strong></td><td>0.06</td><td>0.54</td><td>0.43</td></tr><tr><td>Dry</td><td><strong>0.00</strong></td><td>*0.37</td><td><strong>0.00</strong></td><td><strong>0.00</strong></td><td><strong>0.02</strong></td><td>*0.80</td><td><strong>0.03</strong></td></tr><tr><td>Wet</td><td>0.48</td><td><strong>0.02</strong></td><td><strong>0.00</strong></td><td>0.24</td><td>0.70</td><td>0.56</td><td>0.25</td></tr><tr><td rowspan="3">Adjusted R<sup>2</sup></td><td>All</td><td>0.78</td><td>0.86</td><td>0.59</td><td>0.34</td><td>0.85</td><td>0.72</td><td>0.80</td></tr><tr><td>Dry</td><td>0.79</td><td>*0.89</td><td>0.57</td><td>*0.37</td><td>*0.94</td><td>*0.81</td><td>*0.94</td></tr><tr><td>Wet</td><td>0.85</td><td>0.89</td><td>0.66</td><td>0.23</td><td>0.77</td><td>0.64</td><td>0.82</td></tr><tr><td rowspan="3">Overlap ratio</td><td>All</td><td>0.43</td><td>0.38</td><td>0.43</td><td>0.43</td><td>0.50</td><td>0.27</td><td>0.33</td></tr><tr><td>Dry</td><td>*0.47</td><td>*0.38</td><td>0.35</td><td>0.43</td><td>*0.57</td><td>*0.33</td><td>*0.43</td></tr><tr><td>Wet</td><td>0.25</td><td>0.38</td><td>0.67</td><td>—</td><td>0.40</td><td>0.20</td><td>0.20</td></tr><tr><td rowspan="3">RMSD</td><td>All</td><td>1.13</td><td>1018</td><td>1884</td><td>2165</td><td>1.74</td><td>6595</td><td>8.57</td></tr><tr><td>Dry</td><td>1.20</td><td>*865</td><td>*1810</td><td>2199</td><td>*1.52</td><td>*5481</td><td>*6.87</td></tr><tr><td>Wet</td><td>0.86</td><td>1347</td><td>2083</td><td>1890</td><td>2.00</td><td>7723</td><td>10.50</td></tr><tr><td rowspan="3">RMSD_sys</td><td>All</td><td>0.58</td><td>69</td><td><strong>1416</strong></td><td><strong>1839</strong></td><td>0.93</td><td>1292</td><td>2.08</td></tr><tr><td>Dry</td><td>0.74</td><td>152</td><td>1279</td><td><strong>1933</strong></td><td><strong>1.25</strong></td><td>656</td><td><strong>5.36</strong></td></tr><tr><td>Wet</td><td>0.26</td><td><strong>965</strong></td><td><strong>1763</strong></td><td>900</td><td>0.42</td><td>2401</td><td>6.05</td></tr><tr><td rowspan="3">RMSD_unbias</td><td>All</td><td>0.98</td><td>1015</td><td>1242</td><td>1142</td><td>1.46</td><td>6467</td><td>8.31</td></tr><tr><td>Dry</td><td>0.94</td><td>851</td><td>1281</td><td>1047</td><td>0.87</td><td>5441</td><td>4.30</td></tr><tr><td>Wet</td><td>0.83</td><td>939</td><td>1111</td><td>1662</td><td>1.96</td><td>7340</td><td>8.58</td></tr></tbody></table>

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.

### 3.3.4 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.

### 3.3.5 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).
<figure><img src="image_10.png" /></figure>

**Fig. 6.** Simulated vs. observed woody debris as a function of time after fire for (a) standing dead wood (SDW); (b) downed woody debris (DWD); and (c) total woody debris (TWD). Data for Manitoba and Saskatchewan sites are shown for SDW, and only Manitoba sites are shown for DWD and TWD. Model results are presented by pooling together outputs for all the evaluation sites in the same site cluster, with the thick black (and dashed) line indicating the mean value, and shaded area showing between-site minimum-maximum range. All black lines in the three panels show the result for Manitoba sites, the dashed line in panel (a) shows the result for Saskatchewan sites. Measurements from different sources are shown for both dry (open circles and triangles) and wet (filled dots) measurements. Error bars on the measurement points indicate 90% confidence interval measurement uncertainty. The inset panel shows the overall model-measurement agreement along a 1:1 line for dry (small open circles) and wet (small cross symbol, +) measurements separately.

## 3.4 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²), 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.

## 3.5 Attributing the role of past climate and CO₂ trends in the postfire evolution of carbon fluxes

To attribute the roles of the atmospheric CO₂ and varying climate in postfire forest CO₂ fluxes trajectory, CO₂FIX-CLIMFIX and CO₂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₂ and climate for the site CA-NS1 (Manitoba), we assume that the simulation result by CO₂FIX-CLIMFIX for the time before 1968 (the starting year of meteorological station observations) followed the same curve as CO₂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 post-fire CO₂ 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₂. The EC-measured magnitudes of postfire GPP for different ages of forest after fire are much higher than the simulation result with fixed CO₂ (CO₂FIX-CLIMVAR). The GPP measurements can only be reproduced by the model if the effect of increasing CO₂ 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₂FIX-CLIMVAR simulations, accounting for the effect of rising
<table><caption>**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.</caption><thead><tr><th>Variable</th><th>RMSD</th><th>RTO regression slope</th><th>Overlapping ratio</th><th>Measurement Accuracy*</th></tr></thead><tbody><tr><td>Gross primary production (g C m<sup>-2</sup> yr<sup>-1</sup>)</td><td>76</td><td>–</td><td>–</td><td>100</td></tr><tr><td>Total ecosystem respiration (g C m<sup>-2</sup> yr<sup>-1</sup>)</td><td>76</td><td>–</td><td>–</td><td>200</td></tr><tr><td>Net ecosystem production (g C m<sup>-2</sup> yr<sup>-1</sup>)</td><td>69</td><td>–</td><td>–</td><td>50</td></tr><tr><td>Leaf area index (m<sup>2</sup> m<sup>-2</sup>)</td><td>1.13</td><td>0.76</td><td>0.43</td><td>2</td></tr><tr><td>Biomass carbon (g C m<sup>-2</sup>)</td><td>1018</td><td>1.03</td><td>0.38</td><td>1000</td></tr><tr><td>Woody debris (g C m<sup>-2</sup>)</td><td>1885</td><td>0.35</td><td>0.43</td><td>2000</td></tr><tr><td>Forest floor carbon (g C m<sup>-2</sup>)</td><td>1884</td><td>0.51</td><td>0.43</td><td>1000</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**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₂ with varying climate (GPPCAL-CMCD) are shown for all seven sites in colored lines, with the scenario of fixed CO₂ with varying climate (CO₂FIX-CLIMVA) as a gray line, and the scenario of fixed CO₂ with fixed climate (CO₂FIX-CLIMFIX) as black line for the site CA-NS1. Corresponding eddy-covariance CO₂ 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₂ for each site in the GPPCAL-CMCD simulation.</figcaption></figure>

CO₂ 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₂FIX-CLIMFIX simulation than CO₂FIX-CLIMVAR, indicating that recent climate trends alone decrease the post-fire carbon sequestration in our model simulation. GPP is decreased to a greater extent than TER when climate varies but CO₂ 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₂ caused an increase of 30.5 g C m⁻² yr⁻¹ in mean annual NEP but climate trends alone a decrease by 3.5 g C m⁻² yr⁻¹; their combined effect is an increase of mean annual NEP by 26.9 g C m⁻² yr⁻¹ (Table 7). For the period after the most recent fire event, the GPPCAL-CMCD (varying CO₂ 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₂ 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₂ caused an increase in NEP by 54.1 g C m⁻² yr⁻¹, while varying climate resulted in a decrease by 22.8 g C m⁻² yr⁻¹; their combined effect is an increase in mean annual NEP of 31.2 g C m⁻² yr⁻¹. For US-Bn1 (83 yr), increasing CO₂ caused an increase of 7.0 g C m⁻² yr⁻¹ in mean annual NEP while varying climate a decrease by 7.5 g C m⁻² yr⁻¹, leading to a decrease in mean annual NEP of 0.6 g C m⁻² yr⁻¹.

## 4 Discussion

In the preceding sections, the ORCHIDEE-FM-BF model is evaluated against chronosequence measurements of carbon

<footnote>*Measurement accuracy information is from Goulden et al. (2011). For GPP, TER and NEP, measurement and aggregation accuracy is used, and for other variables, across-landscape sampling accuracy is used.</footnote>
<table><caption>**Table 7.** Effects of varying atmospheric CO₂ and climate and their combined effect on mean annual carbon fluxes (g C m⁻² yr⁻¹) 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.</caption><thead><tr><th rowspan="2">Carbon flux</th><th rowspan="2"></th><th colspan="3">Effect on mean annual carbon flux over the chronosequence period*</th><th colspan="3">Simulated mean annual carbon flux over the chronosequence period</th></tr><tr><th>Climate trends</th><th>Rising CO₂</th><th>Combined</th><th>CO₂FIX-CLIMFIX</th><th>CO₂FIX-CLIMVAR</th><th>GPPCAL-CMCD</th></tr></thead><tbody><tr><td rowspan="4">CA-NS1</td><td>GPP</td><td>-9.3</td><td>91.4</td><td>82.1</td><td>461.4</td><td>452.1</td><td>543.5</td></tr><tr><td>TER</td><td>-5.8</td><td>61.0</td><td>55.2</td><td>454.0</td><td>448.2</td><td>509.2</td></tr><tr><td>NEP</td><td>-3.5</td><td>30.5</td><td>26.9</td><td>7.4</td><td>3.9</td><td>34.3</td></tr><tr><td>GPP</td><td>-35.8</td><td>201.7</td><td>165.9</td><td>681.0</td><td>645.2</td><td>846.9</td></tr><tr><td rowspan="3">CA-SF1</td><td>TER</td><td>-13.0</td><td>147.7</td><td>134.7</td><td>672.7</td><td>659.8</td><td>807.4</td></tr><tr><td>NEP</td><td>-22.8</td><td>54.1</td><td>31.2</td><td>8.3</td><td>-14.6</td><td>39.5</td></tr><tr><td>GPP</td><td>-31.8</td><td>35.1</td><td>3.3</td><td>230.7</td><td>198.9</td><td>234.0</td></tr><tr><td rowspan="2">US-Bn1</td><td>TER</td><td>-24.3</td><td>28.1</td><td>3.9</td><td>225.1</td><td>200.9</td><td>229.0</td></tr><tr><td>NEP</td><td>-7.5</td><td>7.0</td><td>-0.6</td><td>5.6</td><td>-2.0</td><td>5.0</td></tr></tbody></table>

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.

## 4.1 Effect of nudging observed multi-annual average GPP on improving other CO₂ 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 calibration,

which has also been noticed by Zhao et al. (2012) and Lin et al. (2011).

## 4.2 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⁻²) than in Canadian fires (1.2–1.9 kg C m⁻²), 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⁻² and are therefore within the range given by French et al. (2011). However, the simulated emissions at Alaskan sites (0.53 kg C m⁻²) are far below their reported value, and are also lower than those reported by Randerson et al. (2006, 1.56 kg C m⁻² for the site US-Bn3), and Kasischke and Hoy (2012, 3.01 kg C m⁻² during large fire years and 1.69 kg C m⁻² 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⁻² vs. observed ~3.3 kg C m⁻²). Correspondingly, fire carbon emissions are underestimated by 65–80% (simulated 0.5 kg C m⁻² vs. reported 1.5–3 kg C m⁻² from local studies). Besides, the simple combustion fraction scheme used in the model does not

<footnote>*Climate effect was calculated as difference between simulations of CO₂FIX-CLIMVAR and CO₂FIX-CLIMFIX, and CO₂ effect was calculated as difference between simulations of GPPCAL-CMCD and CO₂FIX-CLIMVAR, and combined effect as difference between simulations of GPPCAL-CMCD and CO₂FIX-CLIMFIX.</footnote>
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).

## 4.3 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., 1997, 2001; Wang 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 (Bond-Lamberty et al., 2006).

ORCHIDEE-FM-BF is found to perform better in well-drained 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 (Harden et al., 2006). 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 ORCHIDEE, the poor drainage related hydrological and eco-physiological 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 well-drained. 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.

## 4.4 The role of past climate and CO₂ trends in the postfire evolution of carbon fluxes

### 4.4.1 The role of increasing atmospheric CO₂

We find a significantly positive role of increasing atmospheric CO₂ in explaining the postfire carbon uptake as observed by the chronosequence method. Other studies also report similar positive effects of historical CO₂ in increasing the boreal regional carbon stock (Balshi et al., 2007; Hayes et al., 2011). This is reasonable as the ambient CO₂ 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₂ 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⁻² yr⁻¹). This is because the CO₂ 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–2005 is lower than the site CA-NS5 (615 vs. 741 g C m⁻² yr⁻¹) even though both forests share the same ambient CO₂ concentration, as the forest age at the latter site is much older (155 vs. 25 yr).

This finding might lead to the hypothesis that, for the same CO₂ 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₂. If correct, this hypothesis implies that the current forest carbon sink should not be ascribed only to either forest succession or CO₂ increase, but must be due to a combination of both, i.e., the forest succession accompanied by concurrent atmospheric CO₂ 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₂ fertilization effect might be overestimated.

### 4.4.2 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₂ 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₂ is prescribed to increase. Our results also support a negative role of varying climate in the postfire CO₂ 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 post-fire CO₂ 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₂FIX-CLIMVAR than CO₂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 (Harden et al., 2006; O'Donnell et al., 2011; Yi et al., 2009b, 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.

### 4.4.3 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₂ increase, nutrient availability (which is not included in ORCHIDEE-FM-BF), permafrost and soil dynamics, and forest succession.

## 5 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 struc-
ture, with self-thinning and postfire snag decomposition also
being simulated. By using a simulation protocol with recur-
rent 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 CO2 fertilization in explaining the current post-fire 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 CO2 fertilization but must be due to a combination of both, i.e., the successional sink has been amplified by the CO2 fertilization.

Despite lack of some important features in boreal forest
ecosystems (which may include poor drainage processes, nu-
trient 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 con-
ditions. The model calibration in this study will allow the re-
gional 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.

**Supplementary material related to this article is**
**available online at http://www.biogeosciences.net/10/
8233/2013 bg-10-8233-2013-supplement.pdf.**

Acknowledgements. This study was supported by the fireCCI project (http://www.esa-fire-cci.org/), funded by the European Space Agency. The authors would like to thank S. Yi for providing the climate data for the meteorological stations of Delta Junction and Fairbanks, Alaska. We thank B. D. Amiro for providing the CO2 flux data and the constructive comments and M. L. Goulden for providing carbon flux data from the Manitoba sites. We also thank the two anonymous reviewers for providing constructive reviews of previous versions of this paper. This work used eddy covariance data acquired by the FLUXNET community and in particular by the following networks: AmeriFlux (US Department of Energy, Biological and Environmental Research, Terrestrial Carbon Program (DE-FG02-04ER63917 and DE-FG02-04ER63911)), AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada (supported by CFCAS, NSERC, BIOCAP, Environment Canada, and NRCan), GreenGrass, KoFlux, LBA, NECC, OzFlux, TCOS-Siberia, USCCC. We acknowledge the financial support to the eddy covariance data harmonization provided by CarboEuropeIP, FAO-GTOS-TCO, iLEAPS, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, Université Laval and Environment Canada and US Department of Energy and the database

development and technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, University of California – Berkeley, University of Virginia.

Edited by: P. Stoy

<figure><img src="image_12.png" /></figure>

The publication of this article is
financed by CNRS-INSU.

**References**

Amiro, B. D., Barr, A. G., Black, T. A., Iwashita, H., Kljun, N., McCaughey, J. H., Morgenstern, K., Murayama, S., Nesic, Z., Orchansky, A. L., and Saigusa, N.: Carbon, energy and water fluxes at mature and disturbed forest sites, Saskatchewan, Canada, Agr. Forest Meteorol., 136, 237–251, doi:10.1016/j.agrformet.2004.11.012, 2006.

Amiro, B. D., Barr, A. G., Barr, J. G., Black, T. A., Bracho, R., Brown, M., Chen, J., Clark, K. L., Davis, K. J., Desai, A. R., Dore, S., Engel V., Fuentes J. D., Goldstein A. H., Goulden M. L., Kolb T. E., Lavigne M. B., Law B. E., Margolis H. A., Martin T., McCaughey J. H., Misson L., Montes-Helu M., Noormets A., Randerson J. T., Starr G., and Xiao J.: Ecosystem carbon dioxide fluxes after disturbance in forests of North America, J. Geophys. Res., 115, G00K02, doi:10.1029/2010JG001390, 2010.

Anderson, R. S., Hallett, D. J., Berg, E., Jass, R. B., Toney, J. L., De Fontaine, C. S., and DeVolder, A.: Holocene development of Boreal forests and fire regimes on the Kenai Lowlands of Alaska, The Holocene, 16, 791–803, doi:10.1191/0959683606hol966rp, 2006.

Balshi, M. S., McGuire, A. D., Zhuang, Q., Melillo, J., Kick-
lighter, D. W., Kasischke, E., Wirth, C., Flannigan, M., Harden,
J., Clein, J. S., Burnside, T. J., McAllister, J., Kurz, W. A.,
Apps, M., and Shvidenko, A.: The role of historical fire dis-
turbance in the carbon dynamics of the pan-boreal region:
A process-based analysis, J. Geophy. Res.-Biogeosci., 112,
G02029, doi:10.1029/2006JG000380, 2007.

Balshi, M. S., McGuire, A. D., Duffy, P., Flannigan, M., Kick-
lighter, D. W., and Melillo, J.: Vulnerability of carbon storage in
North American boreal forests to wildfires during the 21st cen-
tury, Glob. Change Biol., 15, 1491–1510, 2009.

Barichivich, J., Briffa, K. R., Osborn, T. J., Melvin, T. M., and Cae-
sar, J.: Thermal growing season and timing of biospheric carbon
uptake across the Northern Hemisphere, Global Biogeochem.
Cy., 26, GB4015, doi:10.1029/2012GB004312, 2012.

Bellassen, V., Le Maire, G., Dhôte, J. F., Ciais, P., and Viovy, N.: Modelling forest management within a global vegetation model – Part 1: Model structure and general behaviour, Ecol. Model., 221, 2458–2474, doi:10.1016/j.ecolmodel.2010.07.008, 2010.

Bond-Lamberty, B. and Gower, S. T.: Decomposition and Frag-
mentation of Coarse Woody Debris: Re-visiting a Boreal
Black Spruce Chronosequence, Ecosystems, 11, 831–840,
doi:10.1007/s10021-008-9163-y, 2008.
Bond-Lamberty, B., Wang, C., and Gower, S. T.: Net primary production and net ecosystem production of a boreal black spruce wildfire chronosequence, Glob. Change Biol., 10, 473–487, 2004.

Bond-Lamberty, B., Gower, S. T., Goulden, M. L., and McMillan, A.: Simulation of boreal black spruce chronosequences: Comparison to field measurements and model evaluation, J. Geophys. Res., 111, G02014, doi:10.1029/2005JG000123, 2006.

Bond-Lamberty, B., Gower, S. T., and Ahl, D. E.: Improved simulation of poorly drained forests using Biome-BGC, Tree Physiol., 27, 703–715, 2007a.

Bond-Lamberty, B., Peckham, S. D., Ahl, D. E., and Gower, S. T.: Fire as the dominant driver of central Canadian boreal forest carbon balance, Nature, 450, 89–92, doi:10.1038/nature06272, 2007b.

Chapin III, F. S., Pamela, A. M., and Harold, A. M.: Principles of Terrestrial Ecosystem Ecology, Springer-Verlag New York, Inc., 2002.

De Groot, W. J., Pritchard, J. M., and Lynham, T. J.: Forest floor fuel consumption and carbon emissions in Canadian boreal forest fires, Can. J. Forest Res., 39, 367–382, doi:10.1139/X08-192, 2009.

Ducoudré, N. I., Laval, K., and Perrier, A.: SECHIBA, a New Set of Parameterizations of the Hydrologic Exchanges at the Land-Atmosphere Interface within the LMD Atmospheric General Circulation Model, J. Climate, 6, 248–273, doi:10.1175/1520-0442(1993)006<0248:SANSOP>2.0.CO;2, 1993.

Euskirchen, E. S., MCGUIRE, A. D., and Chapin, F. S.: Energy feedbacks of northern high-latitude ecosystems to the climate system due to reduced snow cover during 20th century warming, Glob. Change Biol., 13, 2425–2438, doi:10.1111/j.1365-2486.2007.01450.x, 2007.

French, N. H. F., Goovaerts, P., and Kasischke, E. S.: Uncertainty in estimating carbon emissions from boreal forest fires, J. Geophys. Res., 109, D14S08, doi:10.1029/2003JD003635, 2004.

French, N. H. F., de Groot, W. J., Jenkins, L. K., Rogers, B. M., Alvarado, E., Amiro, B., de Jong, B., Goetz, S., Hoy, E., Hyer, E., Keane, R., Law, B. E., McKenzie, D., McNulty, S. G., Ottmar, R., Pérez-Salicrup, D. R., Randerson, J., Robertson, K. M., and Turetsky, M.: Model comparisons for estimating carbon emissions from North American wildland fire, J. Geophys. Res., 116, G00K05, doi:10.1029/2010JG001469, 2011.

Goulden, M. L., Winston, G. C., McMillan, A. M. S., Litvak, M. E., Read, E. L., Rocha, A. V., and Rob Elliot, J.: An eddy covariance mesonet to measure the effect of forest age on land–atmosphere exchange, Glob. Change Biol., 12, 2146–2162, 2006.

Goulden, M., McMillan, A., Winston, G., Rocha, A., Marnies, K., Harden, J., and Bond-Lamberty, B.: Patterns of NPP, GPP, respiration, and NEP during boreal forest succession, Glob. Change Biol., 17, 855–871, doi:10.1111/j.1365-2486.2010.02274.x, 2011.

Gower, S. T., Vogel, J. G., Norman, J. M., Kucharik, C. J., Steele, S. J., and Stow, T. K.: Carbon distribution and aboveground net primary production in aspen, jack pine, and black spruce stands in Saskatchewan and Manitoba, Canada, J. Geophys. Res., 102, 29029–29041, doi:10.1029/97JD02317, 1997.

Harden, J. W., O'Neill, K. P., Trumbore, S. E., Veldhuis, H., and Stocks, B. J.: Moss and soil contributions to the annual net car-

bon flux of a maturing boreal forest, J. Geophys. Res.-Atmos., 102, 28805–28816, doi:10.1029/97JD02237, 1997.

Harden, J. W., Trumbore, S. E., Stocks, B. J., Hirsch, A., Gower, S. T., O'Neill, K. P., and Kasischke, E. S.: The role of fire in the boreal carbon budget, Glob. Change Biol., 6, 174–184, 2000.

Harden, J. W., Meier, R., Silapaswan, C., Swanson, D. K., and McGuire, A. D.: Soil drainage and its potential for influencing wildfires in Alaska, US Geological Survey, Professional Paper, 139–144, 2001.

Harden, J. W., Manies, K. L., Turetsky, M. R., and Neff, J. C.: Effects of wildfire and permafrost on soil organic matter and soil climate in interior Alaska, Glob. Change Biol., 12, 2391–2403, doi:10.1111/j.1365-2486.2006.01255.x, 2006.

Hayes, D. J., McGuire, A. D., Kicklighter, D. W., Gurney, K. R., Burnside, T. J., and Melillo, J. M.: Is the northern high-latitude land-based CO₂ sink weakening?, Global Biogeochem. Cy., 25, GB3018, doi:10.1029/2010GB003813, 2011.

Hayes, D. J., Turner, D. P., Stinson, G., McGuire, A. D., Wei, Y., West, T. O., Heath, L. S., de Jong, B., McConkey, B. G., Birdsey, R. A., Kurz, W. A., Jacobson, A. R., Huntzinger, D. N., Pan, Y., Post, W. M., and Cook, R. B.: Reconciling estimates of the contemporary North American carbon balance among terrestrial biosphere models, atmospheric inversions, and a new approach for estimating net ecosystem exchange from inventory-based data, Glob. Change Biol., 18, 1282–1299, doi:10.1111/j.1365-2486.2011.02627.x, 2012.

Kasischke, E. S., Christensen, N. L., and Stocks, B. J.: Fire, Global Warming, and the Carbon Balance of Boreal Forests, Ecol. Appl., 5, 437–451, doi:10.2307/1942034, 1995.

Kasischke, E. S.: Boreal ecosystems in the global carbon cycle, in: Fire, Climate Change, and Carbon Cycling in the Boreal Forest, edited by: Kasischke, E. S. and Stocks, B. J., 19–30, Springer New York, 2000.

Kasischke, E. S. and Hoy, E. E.: Controls on carbon consumption during Alaskan wildland fires, Glob. Change Biol., 18, 685–699, doi:10.1111/j.1365-2486.2011.02573.x, 2012.

Kasischke, E. S. and Turetsky, M. R.: Recent changes in the fire regime across the North American boreal region – Spatial and temporal patterns of burning across Canada and Alaska, Geophys. Res. Lett., 33, L09703, doi:10.1029/2006GL025677, 2006.

Kasischke, E. S., O'Neill K. P., French, N. H. F., and Bourgeau-Chavez, L. L.: Controls on patterns of biomass burning in Alaskan boreal forests, in: Fire, Climate Change, and Carbon Cycling in the Boreal Forest, edited by: Kasischke, E. S. and Stocks, B. J., 173–196, Springer New York, 2000.

Koven, C. D., Ringeval, B., Friedlingstein, P., Ciais, P., Cadule, P., Khvorostyanov, D., Krinner, G., and Tarnocai, C.: Permafrost carbon-climate feedbacks accelerate global warming, PNAS, 108, 14769–14774, doi:10.1073/pnas.1103910108, 2011.

Kreuzwieser, J., Papadopoulou, E., and Rennenberg, H.: Interaction of flooding with carbon metabolism of forest trees, Plant. Biol., 6, 299–306, doi:10.1055/s-2004-817882, 2004.

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S., and Prentice, I. C.: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Global Biogeochem. Cy., 19, GB1015, doi:10.1029/2003GB002199, 2005.
Kurz, W. A. and Apps, M. J.: A 70-YEAR RETROSPECTIVE ANALYSIS OF CARBON FLUXES IN THE CANADIAN FOREST SECTOR, Ecol. Appl., 9, 526–547, 1999.

Kurz, W. A., Dymond, C. C., White, T. M., Stinson, G., Shaw, C. H., Rampley, G. J., Smyth, C., Simpson, B. N., Neilson, E. T., Trofymow, J. A., Metsaranta, J., and Apps, M. J.: CBM-CFS3: A model of carbon-dynamics in forestry and land-use change implementing IPCC standards, Ecol. Model., 220, 480–504, doi:10.1016/j.ecolmodel.2008.10.018, 2009.

Lin, J. C., Pejam, M. R., Chan, E., Wofsy, S. C., Gottlieb, E. W., Margolis, H. A., and McCaughey, J. H.: Attributing uncertainties in simulated biospheric carbon fluxes to different error sources, Global Biogeochem. Cy., 25, GB2018, doi:10.1029/2010GB003884, 2011.

Litvak, M., Miller, S., Wofsy, S. C., and Goulden, M.: Effect of stand age on whole ecosystem CO₂ exchange in the Canadian boreal forest, J. Geophys. Res., 108, 8225, doi:10.1029/2001JD000854, 2003.

Liu, H. and Randerson, J. T.: Interannual variability of surface energy exchange depends on stand age in a boreal forest fire chronosequence, J. Geophys. Res., 113, G01006, doi:10.1029/2007JG000483, 2008.

Liu, H., Randerson, J. T., Lindfors, J., and Chapin, F. S.: Changes in the surface energy budget after fire in boreal ecosystems of interior Alaska: An annual perspective, J. Geophys. Res., 110, D13101, doi:10.1029/2004JD005158, 2005.

Ma, Z., Peng, C., Zhu, Q., Chen, H., Yu, G., Li, W., Zhou, X., Wang, W., and Zhang, W.: Regional drought-induced reduction in the biomass carbon sink of Canada's boreal forests, PNAS, doi:10.1073/pnas.1111576109, 2012.

Mack, Treseder, K. K., Manies, K. L., Harden, J. W., Schuur, E. A. G., Vogel, J., Randerson, J. T., and Chapin III, F. S.: Recovery of aboveground plant biomass and productivity after fire in mesic and dry black spruce forests of Interior Alaska, Ecosystems, 11, 209–225, doi:10.1007/s10021-007-9117-9, 2008.

Magnani, F., Mencuccini, M., Borghetti, M., Berbigier, P., Berninger, F., Delzon, S., Grelle, A., Hari, P., Jarvis, P. G., Kolari, P., Kowalski, A. S., Lankreijer, H., Law, B. E., Lindroth, A., Loustau, D., Manca, G., Moncrieff, J. B., Rayment, M., Tedeschi, V., Valentini, R., and Grace, J.: The human footprint in the carbon cycle of temperate and boreal forests, Nature, 447, 849–851, doi:10.1038/nature05847, 2007.

Manies, K. L., Harden, J. W., Bond-Lamberty, B. P., and O'Neill, K. P.: Woody debris along an upland chronosequence in boreal Manitoba and its impact on long-term carbon storage, Can. J. Forest Res., 35, 472–482, doi:10.1139/x04-179, 2005.

McGuire, A. D., Anderson, L. G., Christensen, T. R., Dallimore, S., Guo, L., Hayes, D. J., Heimann, M., Lorenson, T. D., Macdonald, R. W., and Roulet, N.: Sensitivity of the carbon cycle in the Arctic to climate change, Ecol. Monogr., 79, 523–555, 2009.

Michaelian, M., Hogg, E. H., Hall, R. J., and Arsenault, E.: Massive mortality of aspen following severe drought along the southern edge of the Canadian boreal forest, Glob. Change Biol., 17, 2084–2094, doi:10.1111/j.1365-2486.2010.02357.x, 2011.

O'Donnell, J. A., Harden, J. W., McGUIRE, A. D., Kanevskiy, M. Z., Jorgenson, M. T., and Xu, X.: The effect of fire and permafrost interactions on soil carbon accumulation in an upland black spruce ecosystem of interior Alaska: implications

for post-thaw carbon loss, Glob. Change Biol., 17, 1461–1474, doi:10.1111/j.1365-2486.2010.02358.x, 2011.

Pan, Y., Birdsey, R. A., Fang, J., Houghton, R., Kauppi, P. E., Kurz, W. A., Phillips, O. L., Shvidenko, A., Lewis, S. L., Canadell, J. G., Ciais, P., Jackson, R. B., Pacala, S. W., McGuire, A. D., Piao, S., Rautiainen, A., Sitch, S., and Hayes, D.: A Large and Persistent Carbon Sink in the World's Forests, Science, 333, 988–993, doi:10.1126/science.1201609, 2011.

Parton, W., Stewart, J., and Cole, C.: Dynamics of C, N, P and S in grassland soils: a model, Biogeochemistry, 5, 109–131, doi:10.1007/BF02180320, 1988.

Piao, S., Ciais, P., Friedlingstein, P., Peylin, P., Reichstein, M., Luyssaert, S., Margolis, H., Fang, J., Barr, A., Chen, A., Grelle, A., Hollinger, D. Y., Laurila, T., Lindroth, A., Richardson, A. D., and Vesala, T.: Net carbon dioxide losses of northern ecosystems in response to autumn warming, Nature, 451, 49–52, doi:10.1038/nature06444, 2008.

Pietsch, S. A., Hasenauer, H., Kucera, J., and Cermák, J.: Modeling effects of hydrological changes on the carbon and nitrogen balance of oak in floodplains, Tree Physiol., 23, 735–746, 2003.

Pregitzer, K. and Euskirchen, E.: Carbon cycling and storage in world forests: biome patterns related to forest age, Glob. Change Biol., 10, 2052–2077, 2004.

Randerson, J. T., Liu, H., Flanner, M. G., Chambers, S. D., Jin, Y., Hess, P. G., Pfister, G., Mack, M. C., Treseder, K. K., Welp, L. R., Chapin, F. S., Harden, J. W., Goulden, M. L., Lyons, E., Neff, J. C., Schuur, E. A. G., and Zender, C. S.: The Impact of Boreal Forest Fire on Climate Warming, Science, 314, 1130–1132, doi:10.1126/science.1132075, 2006.

Ryan, M. G., Binkley, D., and Fownes, J. H.: Age-Related Decline in Forest Productivity: Pattern and Process, in: Advances in Ecological Research, Volume 27, edited by: Begon, M. and Fitter, A. H., 213–262, Academic Press, available at: http://www.sciencedirect.com/science/article/pii/S0065250408600094 (last access: 27 October 2013), 1997.

Soil Landscapes of Canada Working Group, Soil Landscapes of Canada, version 3.2, Agriculture and Agri-Food Canada, digital map and database at 1:1 million scale, http://sis.agr.gc.ca/cansis/nsdb/slc/v3.2/index.html (last access: 07 February 2013), 2010.

Stinson, G., Kurz, W. A., Smyth, C. E., Neilson, E. T., Dymond, C. C., Metsaranta, J. M., Boisvenue, C., Rampley, G. J., Li, Q., White, T. M., and Blain, D.: An inventory-based analysis of Canada's managed forest carbon dynamics, 1990 to 2008, Glob. Change Biol., 17, 2227–2244, doi:10.1111/j.1365-2486.2010.02369.x, 2011.

Stocks, B. J.: Fire behavior in immature jack pine, Can. J. Forest Res., 17, 80–86, doi:10.1139/x87-014, 1987.

Stocks, B. J.: Fire behavior in mature jack pine, Can. J. Forest Res., 19, 783–790, doi:10.1139/x89-119, 1989.

Stocks, B. J., Mason, J. A., Todd, J. B., Bosch, E. M., Wotton, B. M., Amiro, B. D., Flannigan, M. D., Hirsch, K. G., Logan, K. A., Martell, D. L., and Skinner, W. R.: Large forest fires in Canada, 1959–1997, J. Geophys. Res.-Atmos., 107, 8149, doi:10.1029/2001JD000484, 2002.

Trumbore, S. E. and Harden, J. W.: Accumulation and turnover of carbon in organic and mineral soils of the BOREAS northern study area, J. Geophys. Res., 102, 28817–28830, doi:10.1029/97JD02231, 1997.
Turetsky, M. R., Kane, E. S., Harden, J. W., Ottmar, R. D., Manies, K. L., Hoy, E., and Kasischke, E. S.: Recent acceleration of biomass burning and carbon losses in Alaskan forests and peatlands, Nat. Geosci., 4, 27–31, doi:10.1038/ngeo1027, 2011.

Wang, C., Bond-Lamberty, B., and Gower, S.: Carbon distribution of a well-and poorly-drained black spruce fire chronosequence, Glob. Change Biol., 9, 1066–1079–1066–1079, 2003.

Wickland, K. P. and Neff, J. C.: Decomposition of soil organic matter from boreal black spruce forest: environmental and chemical controls, Biogeochemistry, 87, 29–47, doi:10.1007/s10533-007-9166-3, 2007.

Yi, S., Manies, K., Harden, J., and McGuire, D.: Characteristics of organic soil in black spruce forests: implications for the application of land surface and ecosystem models in cold regions, Geophys. Res. Lett., 36, L05501, doi:10.1029/2008GL037014, 2009a.

Yi, S., McGuire, A. D., Harden, J., Kasischke, E., Manies, K., Hinzman, L., Liljedahl, A., Randerson, J., Liu, H., Romanovsky, V., Marchenko, S., and Kim, Y.: Interactions between soil thermal and hydrological dynamics in the response of Alaska ecosystems to fire disturbance, J. Geophys. Res.-Biogeo., 114, G02015, doi:10.1029/2008JG000841, 2009b.

Yi, S., McGuire, A. D., Kasischke, E., Harden, J., Manies, K., Mack, M., and Turetsky, M.: A dynamic organic soil biogeochemical model for simulating the effects of wildfire on soil environmental conditions and carbon dynamics of black spruce forests, J. Geophys. Res., 115, G04015, doi:10.1029/2010JG001302, 2010.

Yuan, F., Yi, S., McGuire, A. D., Johnson, K. D., Liang, J., Harden, J., Kasischke, E. S., and Kurz, W.: Assessment of Historical Boreal Forest C Dynamics in Yukon River Basin: Relative Roles of Warming and Fire Regime Change, Ecol. Appl., 22, 2091–2109, doi:10.1890/11-1957.1, 2012.

Zhao, Y., Ciais, P., Peylin, P., Viovy, N., Longdoz, B., Bonnefond, J. M., Rambal, S., Klumpp, K., Olioso, A., Cellier, P., Maignan, F., Eglin, T., and Calvet, J. C.: How errors on meteorological variables impact simulated ecosystem fluxes: a case study for six French sites, Biogeosciences, 9, 2537–2564, doi:10.5194/bg-9-2537-2012, 2012.