Evaluating the impacts of climate variability and cutting and insect defoliation on the historical carbon dynamics of a boreal black spruce forest landscape in eastern Canada

Ecological Modelling 321 (2016) 98–109

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Evaluating the impacts of climate variability and cutting and insect defoliation on the historical carbon dynamics of a boreal black spruce forest landscape in eastern Canada Bin Chen a,c,∗ , M. Altaf Arain b , Jing M. Chen c , Holly Croft c , Robert F. Grant d , Werner A. Kurz e , Pierre Bernier f , Luc Guindon f , David Price g , Ziyu Wang d a

International Institute for Earth System Science, Nanjing University, Nanjing, China McMaster Centre for Climate Change and School of Geography and Earth Sciences, McMaster University, Hamilton, ON, Canada c Department of Geography and Program in Planning, University of Toronto, Toronto, ON, Canada d Department of Renewable Resources, University of Alberta, Edmonton, AB, Canada e Canadian Forest Service, Pacific Forest Centre, Natural Resources Canada, Victoria, BC, Canada f Canadian Forest Service, Laurentian Forestry Centre, Natural Resources Canada, Quebec, QC, Canada g Canadian Forest Service, Northern Forestry Centre, Natural Resources Canada, Edmonton, AB, Canada b

a r t i c l e

i n f o

Article history: Received 24 August 2015 Received in revised form 20 November 2015 Accepted 23 November 2015 Available online 11 December 2015 Keywords: Net biome productivity Climate sensitivity analysis Scenario analysis Canadian land surface scheme

a b s t r a c t In this study, the Carbon and Nitrogen coupled Canadian Land Surface Scheme (CN-CLASS) was used to investigate the impact of climate variability, seasonal weather effects, disturbance, and CO2 fertilization effects on the historical carbon (C) dynamics of an eastern Canadian boreal forest landscape (6275 ha) from 1928 to 2008. The model was parameterized with ecological, soil texture, forest inventory and historical disturbance data and driven by hourly meteorological data constructed from the historical climate records. Before performing the landscape-level simulation, model results were evaluated against sitelevel eddy covariance (EC) measurements. Landscape-level simulated C fluxes showed that the forest ecosystem was a small C sink in all of the years prior to cutting and insect defoliation in 1963, which resulted in the removal of 23849 Mg C from the forest landscape. As a consequence, the study area was a large C source in 1963 (net biome productivity, NBP = −537 g C m−2 yr−1 ). After that, the forest landscape was mainly a net annual C sink, with total ecosystem C stocks increasing from 14.8 to 16.0 kg C m−2 by 2008, during which total biomass increased from 3.1 to 4.2 kg C m−2 . Analysis of landscape-level, agedetrended, simulated C fluxes for the undisturbed forest landscape from 1928 to 2002 indicated that summer temperature was the dominant control on C fluxes with higher temperature causing a much faster increase in landscape-level annual Re than that of GPP (i.e. 12.3 vs. 1.3 g C m−2 yr−1 ◦ C−1 , respectively). Scenario analysis suggested that forest disturbances had a less profound impact on landscape-level C fluxes and stocks compared to inter-annual climate variability in this landscape. Climate sensitivity analysis revealed that landscape-level simulated C fluxes and stocks were sensitive to the change of air temperature, while only dead organic matter (DOM) and soil organic matter (SOM) were sensitive to the change of precipitation. This study will help to explore the impact of future climate change scenarios and forest management on boreal forest landscapes. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Boreal forests cover circumpolar areas in the Northern Hemisphere in both North America and Eurasia, occupying

∗ Corresponding author at: International Institute for Earth System Science, Nanjing University, 22 Hankou Road, Gulou, Nanjing, Jiangsu, China. E-mail address: [email protected] (B. Chen). http://dx.doi.org/10.1016/j.ecolmodel.2015.11.011 0304-3800/© 2015 Elsevier B.V. All rights reserved.

16.6 million km2 area worldwide (Canadian Forest Service, 2005) and accounting for half of the terrestrial biosphere’s carbon (C) stocks (Dixon et al., 1994). In Canada, boreal forests occupy 3.3 million km2 and around 37% of the country’s landmass, which is about 77% of Canada’s forested land (Brandt, 2009). In 1989, Canada’s boreal forest contained approximately 335.5 Mg C ha−1 (Kurz and Apps, 1999). The ability of boreal forests to sequester C is governed by both climate and disturbance regimes (Chertov et al., 2009; Kang et al., 2006). The productivity of boreal forests is influenced by short-term (Bergeron et al., 2007) and long-term

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

variability in climate (Kang et al., 2006; Wang et al., 2013), as warmer springs tend to increase the yearly CO2 uptake due to the earlier beginning of the growing season (Barr et al., 2004; Bergeron et al., 2007; Tanja et al., 2003). In contrast, autumn warming can enhance respiration more than photosynthesis and thus decrease the yearly CO2 uptake of the forest ecosystems (Piao et al., 2008). In addition, higher summer air temperature (Ta ) may reduce net ecosystem productivity (NEP) (Brooks et al., 1998; Dang and Lieffers, 1989; Tang et al., 2010) through concurrent declines in gross primary productivity (GPP) and increases in ecosystem respiration (Re ) (Grant et al., 2009; Griffis et al., 2003). The increases in NEP from higher spring or summer Ta are likely the more important responses in cooler or wetter climates, while the decreases in NEP resulting from higher summer Ta are likely to be more important in warmer or drier climates (Hofgaard et al., 1999; Nishimura and Laroque, 2011). In Canada, black spruce (Picea mariana (Mill.) B.S.P.) stands constitute a large portion of Canadian boreal forest landscapes (Pavlic et al., 2007), and have a greater total ecosystem C content than any other major forest ecosystems in this biome (Gower et al., 1997). Over the last few decades, boreal black spruce forests have been impacted by both natural and anthropogenic disturbances, which have an effect on the forests’ C cycle. Bergeron et al. (2008) compared carbon dioxide (CO2 ) fluxes of an old black spruce site in eastern Canada with a black spruce site that was harvested in 2000 (HBS00) near Chibougamau, Quebec, Canada. They found that the C budget of boreal black spruce forests was much more affected by the disturbance regime, which consequently affected stand age, than the inter-annual climate variability (Bergeron et al., 2008). Bond-Lamberty et al. (2004) measured all the components of net primary productivity (NPP) and NEP in seven black spruce sites which comprised a boreal forest wildfire chronosequence near Thompson, Manitoba, Canada. This study indicated that site-level total NPP was lowest (50–100 g C m−2 yr−1 ) immediately after fire, highest (332 and 521 g C m−2 yr−1 in the dry and wet stands, respectively) 12–20 years after fire, but around 50% lower than this level in the oldest stands (Bond-Lamberty et al., 2004). This study demonstrated the profound impact of wildfire on the rate of C exchange between forest and atmosphere and the need to account for soil drainage, bryophyte production and species succession when modeling boreal C fluxes (Bond-Lamberty et al., 2004). However, most studies in literature are stand specific and, with very few studies investigating the impact of disturbance on the historical C dynamics of boreal forests at landscape scale, particularly for the black spruce landscape that dominates the Canadian boreal forests (Bernier et al., 2010; Wang et al., 2013). In this study, we incorporated the disturbance matrix of the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3) carbon accounting model into the Carbon and Nitrogen coupled Canadian Land Surface Scheme (CN-CLASS) and applied it to a 6275 ha forest landscape in Chibougamau, Quebec from 1928 to 2008. The specific objectives of this study are to (1) compare model predictions of C fluxes with contemporary flux tower measurements at the Eastern Old Black Spruce (EOBS) site within this study area; (2) to investigate the impacts of climate variability, disturbance events and long-term trends in atmospheric CO2 concentrations on the landscape-level historical C dynamics, and (3) to explore how the interaction of climate variability with disturbance events and CO2 fertilization affect landscape-level C fluxes and stocks.

99

Quebec, Canada and is dominated by black spruce forests. Most of the forests in Chibougamau area were originated from major forest fires in late 1800s, although about 22 ha of 120-year old balsam fir (Abies balsamea (L.) Mill) stands are still present in the study area. The study landscape covers 6275 ha area with UTM co-ordinates (Zone 18, NAD 83). The elevation of the study area varies from 368 m in the southwest to 444 m in the north. The major soil types in the study area include sand, sandy loam and organic. Five surface deposit types (fluvio-glacial deposits, deep glacial till, shallow glacial till, deep organic and shallow organic) were indentified in the Chibougamau study area from a vector map of surface deposit type and soil texture properties (Wang et al., 2013). The drainage of these soil ranges from well drained sand to very poor drained organic soils. The 30-year (1971–2000) mean annual air temperature (MAT) and precipitation are 0.0 ◦ C and 961 mm, respectively, at Chapais which is less than 30 km away from EOBS Flux Station main tower site located in this landscape (Bergeron et al., 2007, 2008; Coursolle et al., 2006). The historical disturbance data of Chibougamau forest landscape in Quebec was created and validated at Laurentian Forestry Centre of Canadian Forest Service (Bernier et al., 2010). Six intermediate mosaics of black and white aerial photos (taken during survey flights in 1953–1954, 1959, 1965, 1967–1968, 1969–1970 and 1982) were interpreted to retrieve the disturbance layers (clear cut, partial cut, insect defoliation and infrastructure) for the period of 1928–2003. Moreover, three decennial provincial forest inventory maps (1970–2005) and Quickbird Pansharped Multispectral images (2003) with 0.6 m resolution were also used to obtain more precise disturbance information after the mid 1970s. Currently, the majority of this disturbed landscape (4601 ha) is dominated by different-aged black spruce stands, among which there are 3857 ha pure black spruce stands and approximately 744 ha black spruce stands are mixed with trembling aspen (Populus tremuloides). In the rest of the landscape, there are 123 ha pure trembling aspen forest stands, 984 ha of non-productive area and 455 ha is covered by water. In 2005, the area-weighted mean stand age was 121 and 43 years for undisturbed, and disturbed stands, respectively. This was determined from the photo-interpretation of the 1998 provincial forest cover map at 1:20,000 scale (Bernier et al., 2010). Fig. 1 shows the levels and types of forest disturbance that was recreated from inventory and remote sensing data. There was no wood exploitation in this area until 1950s, after which the study area was partially disturbed. Major disturbance events include partial cuts (25–75% forest crown density removal) in 1950, 1953, 1957 and 1963 in about 130, 14, 11 and 75 ha areas, respectively (Fig. 1). In the 1950s, about 155 ha of the forest were partially cut (<75% removal) and 11 ha forest was clear cut (> 75% removal) (Fig. 1). In 1963, around 683 ha forest was harvested, of which 608 ha was clear cut and 75 ha was partial cut (Fig. 1). In the same year, approximately 6 ha forest was affected by mild insect defoliation (Fig. 1). From 1966 to 1969, about 95 ha forest was clear cut (Fig. 1). From 1982 to 1991, 12 ha forest was clear cut, of which 10 ha of the forest had regeneration protection (Fig. 1). Approximately, 222 ha stands were converted into infrastructure development (Bernier et al., 2010), causing the permanent reduction of forest area (deforestation). These grid cells were not included in the simulations. 2.2. Model overview

2. Materials and methods 2.1. Study landscape and its disturbance history The Chibougamau study area (centered at 49◦ 41 32.9” N, 74◦ 20 31.3” W) is located about 30 km south of Chibougamau,

CN-CLASS is the carbon and nitrogen coupled version of Canadian Land Surface Scheme (CLASS), which was originally developed for coupling with the Canadian Global Climate Model, CGCM (Verseghy et al., 1993; Verseghy, 1991, 2000). CLASS may have up to five surface/vegetation types in each grid cell including

100

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

Fig. 1. Disturbed area in Chibougamau forest landscape (6275 ha) by the disturbance event type during the simulation period. There were no disturbance events before 1950.

needleleaf trees, broadleaf trees, crops, grasses and urban areas. The soil column is divided into three layers (0.10, 0.25 and 3.75 m thick) and snow is treated as an analogous fourth layer with variable depth. In order to address the interactions between the climate and biogeochemical systems, the ecosystem carbon and nitrogen algorithms were implemented in CLASS version 3.0, yielding the CN-CLASS model (Arain et al., 2002, 2006). CN-CLASS includes modules for photosynthesis (sunlit and shaded leaves), canopy conductance, autotrophic respiration, live biomass allocation (i.e. photosynthate, leaf, wood and root pools), short-lived and stable soil C pools and heterotrophic respiration (Rh ) which is estimated as sum of ground surface litter decomposition, dead wood decomposition, dead roots decomposition, short-lived and stable soil organic matter (SOM) decomposition (Arain et al., 2006). A simple tree allometry module with dynamic leaf phenology was also incorporated in CN-CLASS when the model was applied to seven boreal and temperate conifer forests across Canadian continental transect (Arora and Boer, 2005; Yuan et al., 2008). In this study, a disturbance matrix was used to describe the C flux between the donating (pre-disturbance) and receiving (postdisturbance) C pools as the ratios of the C pools after the logging (Bernier et al., 2010).

2.3. Model initialization and parameterization The total area (3825 ha) for model simulation was determined by overlaying the 1928 and 1998 forest cover maps and extracting the area that was forested in both 1928 and 1998 (Bernier et al., 2010). Totally, four forest cover types were identified in this area (Fig. 2a): 2949 ha were black spruce; 325 ha were mixed forest stands of balsam fir, black spruce and jack pine; 117 ha were jack pine and 434 ha were trembling aspen. For surface deposit maps, fluvio-glacial deposits occupied 20%, deep organic occupies 11%; shallow organic occupies 14%; deep glacial till occupies 51% and shallow glacial till occupies 3% of the total studied area (Fig. 2b). Above ground biomass (AGB) from the 1928 inventory was used to initialize the CN-CLASS model (Fig. 2c). The initial wood biomass, heartwood biomass, root biomass and fine root biomass in the model was assumed as 80%, 60%, 20% and 10% of AGB in 1928, respectively. The initial total ecosystem carbon (TEC) was prescribed from CBM-CFS3’s output (Bernier et al., 2010) (Fig. 2d). The initial forest floor litter, dead wood, dead roots, fast mineralized SOM and stable SOM were assumed as 50%, 10%, 10%, 5% and 25% of necromass. Initial canopy and soil temperatures and soil moisture were derived from a 5-year model spin-up run, using

Fig. 2. Maps of initial data as the input to the CN-CLASS model. (a) Forest cover types. NF: non forest; BS: black spruce; BF: balsam fir; JP: jack pine; TA: trembling aspen. (b) Surface deposit types. Fluvio: fluvio-glacial deposits; Organicd: deep organic; Organics: shallow organic; Tilld: deep glacial till; Tills: shallow glacial till. (c) Aboveground tree biomass in 1928. (d) Total ecosystem carbon in 1928.

B. Chen et al. / Ecological Modelling 321 (2016) 98–109 Table 1 Initial carbon/nitrogen characteristics of foliage used to run the CN-CLASS model. Species Black spruce Jack pine Trembling aspen

SLW g m−2

C/N ratio

230.8 296.1 57.09

67.1 65.3 16.0

SLW: Specific leaf weight. Most of these values were assigned from Middleton et al. (1997).

observed meteorology data (2004 and 2005). The spin-up run to 1928 was also conducted to stabilize ecosystem C stocks, at which time AGB was re-initialized with values from CBM-CFS3. Initial Carbon/Nitrogen characteristics of foliage used in the model are given in Table 1. Soil data and disturbance data (in GIS format) were provided by the Historical Carbon Modeling Project of the Fluxnet Canada Research Network (FCRN) (Bernier et al., 2010; Wang et al., 2013). The weather data used to drive the time series of CN-CLASS model from 1928 to 2003 were derived from climate archives in Great Lakes Forestry Centre, Canadian Forest Service, Natural Resources Canada which had been interpolated to the EOBS eddy covariance flux tower site (Wang et al., 2013). From 2004 to 2008, the observed half-hourly meteorology data from EOBS flux tower site were used (Wang et al., 2013). Historical daily weather records from 1928 to 2003 were interpolated to calculate hourly values for use in CNCLASS model. Specifically, daily shortwave radiation was generated from maximum and minimum air temperature and then daily values were interpolated over calculated day length to obtain hourly values using a sine function. Hourly air temperature was also interpolated using sine function with the daily minimum temperature at dawn and maximum temperature occurring three hours after solar noon. Key parameters used to run the CN-CLASS model over the Chibougamau forest landscape are shown in Table 2. The forest landscape was resolved into 6275 grid cells (1 ha each). The soils, inventory, forest cover type and disturbance data in these grid cells were then aggregated into 319 unique combinations of soils, inventory, forest cover types and disturbances, each of which was represented by one model run. The results from these model runs were then allocated to corresponding individual 1 ha grid cells. Spatial variations in radiation, temperature and precipitation, as well interactions between neighboring grid cells, were not accounted for in our study. The model was spun-up from the year 1923 using the same meteorology forcing data (2004 and 2005 meteorology data) repetitively and resetting the wood biomass and root biomass each year to their initial values to keep model pre-running in a steady forest before starting formal simulations on January 1st, 1928 and completing on December 31st, 2008. 2.4. Climate sensitivity analysis and scenario runs To evaluate the relative influences of climatic variables (i.e. air temperature and precipitation) on historical C dynamics in the landscape, model sensitivity tests were conducted for the whole simulation period (1928–2008). For climate variables, the sensitivity was assessed by increasing/decreasing hourly values of (i) Ta by ±1.0 ◦ C and ±2.0 ◦ C; (ii) precipitation by ±10% and ±25%, respectively. To isolate the compounding effects of climatic forcing on C budget, we changed one variable at a time, while all others were kept unchanged during the sensitivity runs. It should be pointed out that this approach may not fully account for the complexities of real world climate changes and their impacts. Eight scenarios were simulated from 1928 to 2008 over the whole landscape to evaluate the impact of disturbance,

101

Table 2 Photosynthesis and respiration parameters used in the model simulations. Parameters

Value

Maximum carboxylation rate, Vcmax (␮mol m−2 s−1 ) Empirical constant related to intercellular CO2 concentration, a1 Empirical parameter for stomatal sensitivity to VPD, D0 (Pa) Extinction coefficient for canopy nitrogen Maximum transfer rate of assimilate to leaf structure, kp (s−1 ) Quantum efficiency of RuBP regeneration (mol e mol−1 quanta) Leaf respiration rate at reference condition, Rref (␮mol m−2 s−1 ) Living wood respiration rate at reference condition (␮mol m−3 s−1 ) Fine root respiration rate at reference condition (␮mol kg−1 C s−1 ) Coarse root respiration rate at reference condition (␮mol kg−1 C s−1 ) Base litter respiration rate (␮mol kg−1 C s−1 ) Base dead wood respiration rate (␮mol kg−1 C s−1 ) Base dead root respiration rate (␮mol kg−1 C s−1 ) Base short-lived SOM respiration rate (␮mol kg−1 C s−1 ) Base stable SOM respiration rate (␮mol kg−1 C s−1 ) Relative change in stem/coarse root respiration for 10 ◦ C (Q10 ) Relative change in leaf/fine root respiration for 10 ◦ C (Q10 ) Q10 in soil/litter respiration when temperature in 10–30 ◦ C Q10 in soil/litter respiration when temperature is under 10 ◦ C Q10 in soil/litter respiration when temperature is over 30 ◦ C One-half field capacity water content One-half saturated water content (porosity) * Wood turnover rate (␮mol kg−1 C s−1 ) Root turnover rate (␮mol kg−1 C s−1 )

20.0 6.0

*

2000.0 0.14 5.0 × 10−6 0.2 0.5 (15 ◦ C) 3.0 (15 ◦ C) 0.5 (15 ◦ C) 0.3 (15 ◦ C) 0.3 (10 ◦ C) 0.5 (10 ◦ C) 1.5 (10 ◦ C) 1.0 (10 ◦ C) 0.20 (10 ◦ C) 2.0 2.0 2.0 6.0 2.0 0.20 0.23 0.05 (1.89% yr−1 ) 0.30 (11.35% yr−1 )

Wood includes both stem and branch wood.

environmental change (CO2 concentration) and climate variability on the historical dynamics of C fluxes and stocks. These scenarios were: (1) Full: base simulation (i.e. all three effects included); (2) No–CO2 : anthropogenic-driven increase in atmospheric CO2 concentration with time was excluded (i.e. used CO2 concentration in 1828 to 1908 for the simulation period of 1928 to 2008); (3) No–dis: disturbance effects excluded; (4) No–dis/CO2 : both CO2 fertilization and disturbance excluded (climate only effects); (5) No–Clim: inter-annual variability in climate variables excluded; (6) No–Clim/CO2 : both climate variability and CO2 fertilization effects excluded (i.e. disturbance only effects); (7) No–Clim/Dis: both climate variability and disturbance effects excluded (i.e. CO2 fertilization only); (8) No–Clim/CO2 /Dis: all three effects were excluded. 3. Results 3.1. Overall CN-CLASS model performance Before conducting landscape-level analysis of the impacts of climate variability on C fluxes, we first evaluated model performance by comparing modeled and measured C fluxes (GPP, Re and NEP) for the EOBS flux tower site. The regression statistics are summarized in Table 3. The root mean square errors (RMSE) of simulated C

102

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

Table 3 Intercepts (a), slopes (b), coefficients of determination (R2 ), root mean square of errors (RMSE) from regressions of hourly CO2 fluxes by the CN-CLASS model vs. hourly-averaged CO2 fluxes measured by EC from 2004 to 2008 in the fetch area of the EOBS flux tower (grid cell ID 1868). Year

2004 2005 2006 2007 2008

a (g C m−2 h−1 )

RMSE (g C m−2 h−1 )

R2

b

n

GPP

Re

NEP

GPP

Re

NEP

GPP

Re

NEP

GPP

Re

NEP

GPP

Re

NEP

0.01 0.01 0.01 0.01 0.01

0.00 0.00 0.00 0.00 0.00

0.01 0.01 0.01 0.01 0.01

0.86 0.86 0.93 0.96 0.93

0.93 0.92 0.92 0.92 0.93

0.74 0.59 0.66 0.68 0.67

0.82 0.81 0.83 0.77 0.81

0.79 0.80 0.80 0.80 0.80

0.74 0.68 0.71 0.66 0.70

0.06 0.07 0.07 0.07 0.07

0.07 0.07 0.07 0.07 0.07

0.07 0.07 0.07 0.07 0.07

7311 7845 7840 7606 7993

2447 2971 2946 2769 2799

4688 5743 5727 5501 5644

fluxes range from 0.06 to 0.07 g C m−2 hr−1 (Table 3). RMSE between modeled and measured fluxes were similar to a RMSE value of 0.08 g C m−2 hr−1 estimated for CO2 flux measurements at EOBS site (Richardson et al., 2006). Thus, variation in observed fluxes unexplained by the CN-CLASS model could be largely attributed to measurement uncertainty. CN-CLASS is able to explain 77% to 83% of the variance of GPP, 79% to 80% of the variance of Re and 66% to 74% of the variance of NEP (Table 3). The RMSE and R2 values indicated that the CN-CLASS model can generally capture the hourly variability in measured CO2 fluxes. 3.2. Contemporary climate variability effects on forest productivity 3.2.1. Modeled carbon fluxes for two years with contrasting weather patterns CO2 fluxes measured at EOBS site for the coolest year (2004, mean annual temperature (MAT) = −0.3 ◦ C) vs. the warmest year

(2006, MAT = 2.3 ◦ C) during the experimental period (2004–2008) provided an opportunity to test how CN-CLASS model simulated CO2 fluxes under contrasting weather. Different seasonal temperatures caused daily NEP, calculated as the sum of hourly CO2 fluxes, to follow different time courses in 2004 and 2006 (Fig. 3). Earlier spring warming in 2006 vs. 2004 hastened the onset of net C uptake from DOY 129 in 2004 (Fig. 3b) to DOY 104 in 2006 (Fig. 3c). Net C uptake attained maximum values between late June to late July in both years (Fig. 3b and c). The EC measured CO2 fluxes indicated that the onset of net C uptake was quite rapid in both of the years 2004 and 2006, suggesting the activation of CO2 fixation by a thermal signal in spring. Net C uptake generally declined after late July in both years, particularly when Ta exceeded 25 ◦ C, indicating adverse effects of higher Ta on NEP during summer (Fig. 3). Figs. 4 and 5 examine more closely the impact of early season warming events on modeled and measured C fluxes. In 2004, measured CO2 fixation did not respond to the first (Apr. 19–24th) of four successive warming events (Fig. 4d), very weakly to the

Fig. 3. Hourly air temperature (a) and daily net ecosystem productivity (NEP) calculated from CO2 fluxes measured by eddy covariance and modeled by CN-CLASS in flux tower fetch area for 2004 and 2006. Lines for modeled values and dots for measured values.

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

103

Fig. 4. Hourly air temperature, gross primary productivity (GPP), ecosystem respiration (Re ) and net ecosystem productivity (NEP) measured by eddy covariance and modeled by CN-CLASS over flux tower fetch area in spring 2004. Lines for modeled values and dots for measured values.

Fig. 5. Hourly air temperature, gross primary productivity (GPP), ecosystem respiration (Re ) and net ecosystem productivity (NEP) measured by eddy covariance and modeled by CN-CLASS over flux tower fetch area in spring 2006. Lines for modeled values and dots for measured values.

104

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

Fig. 6. Hourly air temperature (a) and hourly gross primary productivity (GPP) (b), ecosystem respiration (Re ) (c) and net ecosystem productivity (NEP) (d) measured by eddy covariance and modeled by CN-CLASS in flux tower fetch area in summer 2006. Lines for modeled values and dots for measured values.

second (Apr. 25–30th) and the third (May 1–6th), but strongly to the fourth (May 7–13th). In 2006 (Fig. 5), EC measured CO2 fixation did not respond to the first (Mar. 27th to Apr. 3rd) and the third (Apr. 8–12th) of four successive warming events (Fig. 5), strongly to the second (Apr. 4–5th) and the fourth (Apr. 13–24th). The timing of this thermal signal was simulated with some accuracy by CN-CLASS allowing the earlier onset of the net C uptake to be simulated in 2006 than in 2004 (Figs. 4d and 5d). Given the rapid increase in CO2 fixation following activation, accurate timing of this signal is important to accurately simulate NEP in boreal conifer forest ecosystem. A continued warming event in July 2006 provided an opportunity to test the simulation of these adverse effects by CN-CLASS. During this event (Fig. 6a), higher Ta and VPD in CN-CLASS forced lower stomatal conductivity and hence sharp declines in midafternoon CO2 influxes when Ta > 25 ◦ C on Jul. 7th, Jul. 12th and Jul. 13th, as also apparent in the EC fluxes (Fig. 6b and d). These declines in mid-afternoon CO2 influxes caused the sharp declines in net daily C uptake modeled with warming in 2004 and 2006 (Fig. 3b and c). 3.2.2. Relationships between simulated landscape-level fluxes and meteorological variables The correlation between simulated seven-year mean C flux deviations and the meteorological variables for the undisturbed landscape from 1928 to 2002 was explored. A seven-year average was taken to remove the age effect on the C fluxes. Over this period, mean daily minimum air temperature (Tmin ) ranged from −8 to 1 ◦ C in spring (April and May), while the mean daily maximum air temperature (Tmax ) was around 2 to 14 ◦ C in spring. Summertime (June, July and August) mean daily Tmin and mean daily Tmax was between 6–12 ◦ C, and 17–23 ◦ C, respectively.

The regression of GPP deviations with Tmax in spring and summer indicated that landscape-level annual photosynthesis rate increased with air temperature in spring when the air temperature was relatively low while in summer, the weather is warm and air temperature is not a major factor that limits the photosynthesis rate (Table 4). The Re deviations appear to be positively related to air temperature (Tmin and Tmax ) in both spring and summer (Table 4). In contrast to GPP, NEP deviations showed stronger negative relationships with air temperature in summer due to the increase of Re with temperature in summer. However, NEP deviations had no correlation with temperature in spring because both GPP and Re increase with air temperature in spring at approximately the same rate (Table 4). This highlights the importance of summertime conditions to the C budget of this forest landscape. Similar results were found by Wang et al. (2013). 3.3. Historical disturbance impacts on forest productivity and C stocks 3.3.1. Disturbance impacts on C fluxes and stocks In 1928, CN-CLASS simulated landscape-level GPP, NPP and NEP values were 643, 331 and 105 g C m−2 yr−1 , respectively (Fig. 7a). The modeled results show that disturbance events in 1950, 1953 and 1957 (155 ha partial cut and 11 ha clear-cut, in total), had insignificant impacts on the simulated landscape-level C fluxes (Fig. 7a), and the forest ecosystems in the study area were small C sinks. In 1963, around 18% of the forest landscape was disturbed, which includes 683 ha forests logged (608 ha clear-cut and 75 ha partial cut) and around 6 ha forests influenced by light insect defoliation. The disturbance events totally caused 23849 Mg C removed from the whole forest landscape in 1963, resulting in the study area becoming a large C source (NBP = −537 g C m−2 yr−1 ).

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

105

Table 4 Intercept (a), slope (b), coefficient of determination (R2 ) from regressions of CN-CLASS simulated annual C flux deviations to daily minimum and maximum air temperature in spring and summer for Chibougamau undisturbed forest landscape from 1928 to 2002. Temperature (◦ C)

a (g C m−2 yr−1 )

◦

Daily Tmin in spring (−8 to 1 C) Daily Tmin in summer (6 to 12 ◦ C) Daily Tmax in spring (2 to 14 ◦ C) Daily Tmax in summer (17 to 23 ◦ C)

b (g C m−2 yr−1 ◦ C−1 )

R2

GPP

Re

NEP

GPP

Re

NEP

GPP

Re

NEP

36.50 −8.83 −101.46 −59.14

35.90 −107.72 −90.36 −243.03

0.61 98.89 −11.10 183.88

11.22 1.32 11.96 3.04

11.54 12.28 10.51 11.99

−0.33 −10.97 1.45 −8.95

0.33 0.00 0.53 0.01

0.34 0.17 0.41 0.19

0.00 0.42 0.02 0.33

In 1966, around 59 ha forests were clear cut which caused 2509 Mg C removed from the forest landscape (3825 ha) and the landscape was a small C sink in 1966 (NBP = 10 g C m−2 yr−1 ). In 1969, 28 ha forests were clear cut which caused 1148 Mg C removed from the forest landscape and the landscape was a small C sink in 1969 (NBP = 17 g C m−2 yr−1 ). After 1969, the forest landscape had very few disturbance events acting as a net C sink for most of the years (Fig. 7a). In 1928, CN-CLASS simulated initial landscape-level aboveground tree biomass (AGB) (leaf, wood and bark biomass), below ground tree biomass (BGB) (live fine and coarse root biomass) were 2.3 and 0.7 kg C m−2 , respectively, derived from the inventory data (Fig. 7b). AGB and BGB gradually increased to 3.0 and 0.9 kg C m−2 , respectively, in 1962 (Fig. 7b). In 1963, AGB and BGB were reduced to 2.4 and 0.7 kg C m−2 due to the harvest of 683 ha forest stands in the landscape (Fig. 1). After 1963, both

AGB and BGB gradually increased to 3.3 and 0.9 kg C m−2 , respectively, in 2008 (Fig. 7b). Overall, CN-CLASS simulated an increase of 39040 Mg C (1.0 kg C m−2 or 0.013 kg C m−2 yr−1 ) in AGB over the 3825 ha forested area in the landscape during the 81 years’ simulation period (Fig. 7b). Initial landscape-level dead organic matter (DOM) (forest floor litter, dead wood and belowground dead roots) and soil organic matter (SOM) (short-lived soil organic matter and stable soil organic matter) were 8.4 and 3.8 kg C m−2 , respectively, derived from the CBM-CFS3 model (Fig. 7b). DOM gradually decreased to 7.9 kg C m−2 in 2008 (Fig. 7b). SOM was quite stable during the simulation period and was 3.9 kg C m−2 in 2008. CNCLASS simulated 5% increase in total ecosystem carbon (TEC) during the simulation period of 1928–2008. The mass balance calculation indicated that this 3825 ha of forest landscape gained 0.8 kg C m−2 from 1928 through 2008. The change in TEC stock matched with the cumulative NBP over the study period (1928–2008).

Fig. 7. (a) Simulated landscape-level mean annual C fluxes: gross primary productivity (GPP), net primary productivity (NPP), net ecosystem productivity (NEP) and net biome productivity (NBP). Positive values are net additions of carbon to the ecosystem and negative values are net removal of carbon from the ecosystem; (b) simulated landscape-level mean annual C stocks.

106

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

Fig. 8. Standard deviations of simulated landscape-level aboveground biomass (AGB) and net biome productivity (NBP) from 1928 to 2008.

3.3.2. Spatial heterogeneity of forest productivity and above ground biomass Standard deviation of landscape-level NBP and AGB represents the spatial heterogeneity of the respective variables in the study area. Our analysis showed that from 1928 to 1949, standard deviations of NBP ranged from 32 to 74 g C m−2 yr−1 (Fig. 8). In this period, no disturbance events took place and the spatial heterogeneity of the forest landscape productivity was predominantly determined by other factors, such as stand age-class structure and site fertility. Standard deviation of NBP increased to 122 g C m−2 yr−1 in 1950 (Fig. 8), when the first forest stands (130 ha) were partial cut (Fig. 1). Standard deviation of NBP reached the highest value (1123 g C m−2 yr−1 ) in 1963 (Fig. 8), when around 683 ha forest was harvested in this landscape (Fig. 1). Standard deviation of

NBP decreased to low value (33 g C m−2 yr−1 ) in 1970 and onwards (Fig. 8), when no disturbance events occurred in the study area (Fig. 1). The AGB standard deviation gradually increased from 1.4 to 2.1 kg C m−2 , indicating that the impact of disturbance events to spatial heterogeneity of AGB was insignificant for the study period (1928–2008). The small increase of AGB standard deviation may be attributed to the difference of forest age structure in the landscape between 2008 and 1928. 3.4. Climate sensitivity and scenario analysis A sensitivity analysis was performed to investigate the response of C fluxes and stocks to air temperature and precipitation

Fig. 9. Sensitivity of simulated landscape-level carbon fluxes and stocks to air temperature. (a) gross primary productivity, GPP; (b) autotrophic respiration, Ra; (c) heterotrophic respiration, Rh ; (d) net ecosystem productivity, NEP; (e) aboveground biomass, AGB; (f) belowground biomass, BGB; (g) dead organic matter, DOM; and (h) soil organic matter, SOM.

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

107

Fig. 10. Sensitivity of simulated landscape-level carbon fluxes and stocks to precipitation. (a) gross primary productivity, GPP; (b) autotrophic respiration, Ra; (c) heterotrophic respiration, Rh ; (d) net ecosystem productivity, NEP; (e) aboveground biomass, AGB; (f) belowground biomass, BGB; (g) dead organic matter, DOM; and (h) soil organic matter, SOM.

(Figs. 9 and 10, respectively). Fig. 9a–d demonstrates that the CN-CLASS simulated C fluxes (GPP, Ra , Rh and NEP) and C stocks (Fig. 9e–h) are sensitive to the air temperature. For C stocks, AGB, BGB, DOM and SOM were sensitive to the change of air temperature due to the temperature effects on the foliage, wood, root turnover rates and C losses through the decomposition (Fig. 9). In contrast to the air temperature simulations, landscape-level C fluxes and biomass were not sensitive to the change in precipitation (Fig. 10). DOM and SOM were sensitive to the change of precipitation (Fig. 10) due to the soil moisture effects on the C losses through the decomposition. Our scenario simulations revealed long-term impacts of all three effects: disturbance, rising atmospheric CO2 concentrations (i.e. the CO2 fertilization effect), and interannual climate variability, on C stock accumulation and C fluxes, but the degree of their impacts was variable (Table 5). When comparing the impact of the three different effects individually, with respect to the null model simulation (none of the three effects were included when comparing the

absolute percentage difference in the variable size), we observed that climate variability effects had the most significant impacts on AGB, BGB and SOM (Table 5). In our simulation runs, inter-annual climate variability alone produced a decrease of 11.6% in AGB and a decrease of 15.7% in BGB and a decrease of 2.2% in SOM. In contrast, all the three effects had insignificant effects on DOM. For the C fluxes, climate variability had the most significant effects on simulated C fluxes (GPP, Re , NEP and NPP). Interannual climate variability alone decreased GPP, Re , NEP and NPP by 20.8%, 14.2%, 59.5% and 26.0%, respectively. 4. Discussion 4.1. Climate variability The response of CN-CLASS simulated NEP to warming between 2004 and 2006 was the net result of two opposing responses: a rise with warming at Ta < 15 ◦ C and a decline with warming at Ta > 20 ◦ C

Table 5 Simulated landscape-level total carbon stocks and fluxes calculated for eight different model scenarios, where model sensitivity to the three effects (i.e. climate variability, CO2 fertilization and disturbances) was tested, separately and in combination. For the stocks, model scenarios were run from 1928 to 2008 and the final values for 2008 are shown, while for the fluxes means of annual values from 2004 to 2008 are shown. In all model runs, mean values (1928 to 2005) of the climate variables were used. Numbers in parenthesis show the percentage change in the stock or flux value, as a result of including the given effect(s) in the model, with respect to the null model simulation (i.e. none of the three effects were included). The sign in the brackets indicate an increase (+) or a decrease (−) with respect to the null model simulation. Stocks (kg C m−2 )

Fluxes (g C m−2 yr−1 )

Scenarios

AGB (%)

BGB (%)

DOM (%)

SOM (%)

GPP (%)

Ra (%)

Rh (%)

NEP (%)

None Clim only CO2 only Dis only Clim & CO2 only Clim & dis only Dis & CO2 only All three effects

3.4(0.0) 3.0(−11.6) 3.4(0.0) 3.4(0.0) 3.0(−11.6) 3.0(−11.6) 3.4(0.0) 2.8(−18.1)

1.0(0.0) 0.9(−15.7) 1.0(0.0) 1.0(0.0) 0.9(−15.7) 0.9(−15.7) 1.0(0.0) 0.9(−15.3)

7.7(0.0) 7.8(0.3) 7.7(0.0) 7.7(0.0) 7.8(0.3) 7.8(0.3) 7.7(0.0) 7.8(0.9)

3.9(0.0) 3.8(−2.2) 3.9(0.0) 3.9(0.0) 3.8(−2.2) 3.8(−2.2) 3.9(0.0) 3.8(−1.9)

713.7(0.0) 565.5(−20.8) 713.7(0.0) 713.7(0.0) 565.5(−20.8) 565.5(−20.8) 713.7(0.0) 558.5(−21.7)

356.2(0.0) 300.9(−15.5) 356.2(0.0) 356.2(0.0) 300.9(−15.5) 300.9(−15.5) 356.2(0.0) 300.0(−15.8)

254.5(0.0) 222.8(−12.5) 254.5(0.0) 254.5(0.0) 222.8(−12.5) 222.8(−12.5) 254.5(0.0) 216.7(−14.9)

103.0(0.0) 41.7(−59.5) 103.0(0.0) 103.0(0.0) 41.7(−59.5) 41.7(−59.5) 103.0(0.0) 41.8(−59.4)

Legend: clim—interannual variability in climate; dis—disturbance; CO2 —CO2 fertilization; AGB—aboveground tree biomass; BGB—belowground tree biomass; DOM—dead organic matter, including dead wood, dead root and forest floor litter; SOM—soil organic matter; GPP—gross primary productivity; Re —ecosystem respiration; NEP—net ecosystem productivity; NPP—net primary productivity.

108

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

(Fig. 3). The rise in NEP with warming at lower Ta in spring can be attributed to three sources in the model: (1) A lengthening of the net C uptake period by 25 days in spring (modeled vs. measured lengthening from 2004 to 2006 in days: 25 vs. 25 at EOBS site in Figs. 4 and 5) with warming of ca. 2 ◦ C in 2005 and 2006 vs. 2004 (Table 3). Lengthening arose from seasonal temperature effects on photosynthetic capacity (Arain et al., 2002; Bergh and Linder, 1999; Jarvis and Linder, 2000; Krishnan et al., 2008; Medhurst et al., 2006). Thus, the earlier spring warming would increase the length of growing season. (2) A rise in carboxylation activity based on parameters (i.e. Tmin and Topt ) for the temperature response function in the CNCLASS model. The minimum and optimum temperatures were set as 0 and 15 ◦ C for Rubisco activity. Thus, Vcmax will increase with the air temperature when the air temperature was between 0 and 15 ◦ C. (3) When the air temperature is above zero for several days, the soil started to thaw and the snow above the soil would begin to melt resulting to the increase of the soil moisture. The increase of soil moisture in early spring would trigger the root nutrient uptake from the soil. More rapid nitrogen uptake by roots maintained or lowered foliar C:N ratio and increased the leaf Rubisco-N and chlorophyll content during spring warming so that N constraints on CO2 fixation remained constant or declined with rising Ta . The decline in NEP with warming at higher Ta in summer was attributed to the following sources in the model: (1) A sharp decline in mid-afternoon stomatal conductance under higher VPD modeled from the parameter D0 which represents the sensitivity of stomatal conductivity to VPD. (2) A rise in Ra and Rh driven by the Q10 function of Ta and Ts . The net effect of these contrasting responses to warming was to increase NEP modeled with warming of 1–2 ◦ C from 2004 to 2008. This effect may be continued with further warming up to mean annual temperature (MAT) of 6–8 ◦ C above which water may limit NEP at mid-continental sites (SOBS and NOBS) with lower precipitation (Grant et al., 2009). Rises in NEP with Ta modeled here were consistent with the experimental findings of Medhurst et al. (2006) that an increase in Ta of 3 ◦ C raised NEP by around one-third in Norway spruce after 3 years of growth in whole tree chambers in northern Sweden (MAT 2.3 ◦ C). However, raising Ta by more than 3 ◦ C above a mean annual value of 7.5 ◦ C caused modeled NEP of a temperate conifer to decline from water stress (Thornley and Cannell, 1996) as found by Grant et al. (2009) for mid-continental sites at SOBS and NOBS. The forest C stocks are being affected by both climate change and disturbances. Atmospheric composition (i.e. CO2 fertilization) and climate variability directly affect the forest productivity by increasing the length of growing season, increasing the gradient between ambient and intercellular CO2 concentration and the closure of leaf stomata due to higher air temperature and larger VPD in summer and so on. In general, there is a negative relationship between temperature and annual C sequestration rates, with higher summer temperatures causing a faster increase in Re than that of GPP (Arain et al., 2002). The results of our study are consistent with this previous finding, indicating that in the absence of disturbances, air temperature in summer was the dominant control on the boreal forest productivity. Furthermore, our climate sensitivity analysis indicated that temperature impacts on C fluxes (i.e. GPP, Ra , Rh and NEP) and C stocks (i.e. AGB, BGB, DOM and SOM) were highly nonlinear. Within the prescribed Ta range (±1.0 and +2.0 ◦ C), all C fluxes and stocks appeared sensitive to Ta , because temperature is the most significant constraint factor on photosynthetic and respiratory biochemical processes. 4.2. Disturbance impacts A large portion of the forest in this study area originated from major forest fires in the late 1800s, although some small pockets of “old growth” forests dominated by balsam fir are still present

since the late 1800s. Approximately, 25% of the study area’s 3825 ha were further disturbed at some time between 1928 and 2008. The photointerpretation of forest cover map revealed that in 1928 this landscape was covered with different aged forests, which included 478, 1564, 1627, 476, 577 and 22 ha stands that were 10, 30, 50, 70 and 90 years old, respectively. The area-weighted mean stand age was 47 years in 1928 indicating relatively younger age of forests in the Chibougamau landscape. Younger forests were underrepresented in the sample plot due to the operational sampling bias toward older stands and the low prevalence of fires in Quebec’s boreal forests during most of the 20th century (Bergeron et al., 2006; Bernier et al., 2010). This underrepresentation has resulted in a possible bias in the strata-level estimates of aboveground biomass in 1928 (Bernier et al., 2010). Results of our study suggest that disturbances had a less profound impact on landscape-level C fluxes and stocks compared to interannual climate variability in this landscape (Table 4). Our results are supported by the findings of (Bernier et al., 2010), who used CBM-CFS3 to reconstruct and model 71 years of forest growth in the same forest landscape as in this study. Bergeron et al. (2008) investigated the impact of forest harvest on CO2 fluxes of black spruce ecosystems in eastern North America by comparing CO2 fluxes over a mature black spruce stand (EOBS) with a site that was harvested in 2000 (HBS00) which has similar soil parent material, site fertility, climate and pre-harvest species composition (Bergeron et al., 2008). Earlier, it was pointed out that natural disturbance events, land use history, atmospheric CO2 concentration and inter-annual climate variability should be considered for the model simulation of biosphere-atmosphere CO2 fluxes (Desai et al., 2007). Our study is a step forward in this direction, indicating that once adequately parameterized, a process-based ecosystem model such as CN-CLASS can become a valuable tool to study the climate variability, forest management, land use and disturbance impacts on the C cycle in boreal forest landscapes.

5. Conclusions In this study, the CN-CLASS model was used to investigate seasonal weather effects and the impacts of climate variability and the disturbance events on the historical C dynamics of a boreal black spruce forest landscape (6275 ha) in eastern Canada from 1928 to 2008. The study results indicated that: (1) The forest landscape in this study was a small C sink of 51 g C m−2 yr−1 in 2008 while it was a larger C sink (105 g C m−2 yr−1 ) at the beginning of the simulation period in 1928. (2) Analysis of landscape-level, age-detrended, simulated C fluxes for the undisturbed forest landscape from 1928 to 2002 indicated that summer temperature was the dominant control on C fluxes with higher temperature causing a much faster increase in landscape-level annual Re than that of GPP (i.e. 12.3 vs. 1.3 g C m−2 yr−1 ◦ C−1 , respectively). Thus, the decrease in landscape-level annual NEP with summer temperature (mean daily Tmin or Tmax in summer) was attributed to the faster increase in landscape-level annual Re than that of GPP with the mean daily Tmin or Tmax in summer. (3) Climate sensitivity analysis indicated that landscape-level C fluxes and stocks were sensitive to the change of air temperature while they were not sensitive to the change of precipitation. The scenario analysis revealed that the impacts of inter-annual climate variability were significant to the simulated landscape-level C fluxes and stocks while other factors (i.e. disturbance and CO2 fertilization effects) were insignificant. Our study highlights the importance of historical climate, forest inventory, disturbance data and process-based land surface model

B. Chen et al. / Ecological Modelling 321 (2016) 98–109

in investigating the impacts of climate variability and disturbance on the boreal forest landscape. Acknowledgements This research is supported by National Natural Science Foundation of China (Grant number: 41503070). This study was conducted as a part of the historical carbon modeling project of the Canadian Carbon Program (CCP). CCP was supported through funding from Canadian Foundation for Climate and Atmospheric Sciences (CFCAS) and funding from National Science and Engineering Research Council (NSERC). We also acknowledge support from Natural Resources Canada (NRC) and the Québec Ministère des Ressources naturelles et de la Faune for the provision of the 1928 inventory data, the 1998 forest cover map and historical climate records. CN-CLASS model runs were made at the SHARCNET at the McMaster University. We want to thank the two anonymous reviewers whose constructive comments have helped us improve our text. References Arain, M.A., Black, T.A., Barr, A.G., Jarvis, P.G., Massheder, J.M., Verseghy, D.L., Nesic, Z., 2002. Effects of seasonal and interannual climate variability on net ecosystem productivity of boreal deciduous and conifer forests. Can. J. For. Res. 32 (5), 878–891. Arain, M., Altaf Yuan, F., Andrew Black, T., 2006. Soil-plant nitrogen cycling modulated carbon exchanges in a western temperate conifer forest in Canada. Agric. For. Meteorol. 140 (1-4), 171–192. Arora, V.K., Boer, G.J., 2005. A parameterization of leaf phenology for the terrestrial ecosystem component of climate models. Global Change Biol. 11 (1), 39–59. Barr, A.G., Black, T.A., Hogg, E.H., Kljun, N., Morgenstern, K., Nesic, Z., 2004. Interannual variability in the leaf area index of a boreal aspen-hazelnut forest in relation to net ecosystem production. Agric. For. Meteorol. 126 (3-4), 237–255. Bergeron, O., Hank, A.M., T. Andrew, B., Carole, C., Allison, L.D., Alan, G.B., Steven, C.W., 2007. Comparison of carbon dioxide fluxes over three boreal black spruce forests in Canada. Global Change Biol. 13 (1), 89–107. Bergeron, O., Margolis, H.A., Coursolle, C., Giasson, M.-A., 2008. How does forest harvest influence carbon dioxide fluxes of black spruce ecosystems in eastern North America? Agric. For. Meteorol. 148 (4), 537–548. Bergeron, Y., Cyr, D., Drever, C.R., Flannigan, M., Gauthier, S., Kneeshaw, D., et al., 2006. Past, current, and future fire frequencies in Quebec’s commercial forests: implications for the cumulative effects of harvesting and fire on age-class structure and natural disturbance-based management. Can. J. For. Res. 36 (11), 2737–2744. Bergh, J., Linder, S., 1999. Effects of soil warming during spring on photosynthetic recovery in boreal Norway spruce stands. Global Change Biol. 5 (3), 245–253, http://dx.doi.org/10.1046/j. 1365-2486.1999.00205.x. Bernier, P.Y., Guindon, L., Kurz, W.A., Stinson, G., 2010. Reconstructing and modelling 71 years of forest growth in a Canadian boreal landscape: a test of the CBM-CFS3 carbon accounting model. Can. J. For. Res. 40 (1), 109–118. Bond-Lamberty, B., Wang, C., Gower, T.S., 2004. Net primary production and net ecosystem production of a boreal black spruce wildfire chronosequence. Glob. Change Biol. 10 (4), 473–487. Brandt, J.P., 2009. The extent of the North American boreal zone. Environ. Rev. 17, 101–161. Brooks, J.R., Flanagan, L.B., Ehleringer, J.R., 1998. Responses of boreal conifers to climate fluctuations: indications from tree-ring widths and carbon isotope analyses. Can. J. For. Res. 28 (4), 524–533, http://dx.doi.org/10.1139/x98-018. Canadian Forest Service, 2005. The State of Canada’s Forests 2004–2005: The Boreal Forest. C. F. S. Natural Resources Canada Headquarters, Planning, Operations and Information Branch, Ottawa, pp. 96 p. Chertov, O., Bhatti, J.S., Komarov, A., Mikhailov, A., Bykhovets, S., 2009. Influence of climate change, fire and harvest on the carbon dynamics of black spruce in Central Canada. For. Ecol. Manage. 257 (3), 941–950. Coursolle, C., Margolis, H.A., Barr, A.G., Black, T.A., Amiro, B.D., McCaughey, J.H., et al., 2006. Late-summer carbon fluxes from Canadian forests and peatlands along an east-west continental transect. Can. J. For. Res. 36 (3), 783–800. Dang, Q.L., Lieffers, V.J., 1989. Climate and annual ring growth of black spruce in some Alberta peatlands. Can. J. Bot. 67 (6), 1885–1889, http://dx.doi.org/10.1139/b89239.

109

Desai, A.R., Moorcroft, P.R., Bolstad, P.V., Davis, K.J., 2007. Regional carbon fluxes from an observationally constrained dynamic ecosystem model: Impacts of disturbance, CO2 fertilization, and heterogeneous land cover. J. Geophys. Res. 112 (G1), G01017. Dixon, R.K., Solomon, A.M., Brown, S., Houghton, R.A., Trexier, M.C., Wisniewski, J., 1994. Carbon pools and flux of global forest ecosystems. Science 263 (5144), 185–190. Gower, S.T., Vogel, J.G., Norman, J.M., Kucharik, C.J., Steele, S.J., Stow, T.K., 1997. Carbon distribution and aboveground net primary production in aspen, jack pine, and black spruce stands in Saskatchewan and Manitoba, Canada. J. Geophys. Res. 102 (D24), 29029–29041. Grant, R.F., Margolis, H.A., Barr, A.G., Black, T.A., Dunn, A.L., Bernier, P.Y., Bergeron, O., 2009. Changes in net ecosystem productivity of boreal black spruce stands in response to changes in temperature at diurnal and seasonal time scales. Tree Physiol. 29 (1), 1–17, http://dx.doi.org/10.1093/treephys/tpn004. Griffis, T.J., Black, T.A., Morgenstern, K., Barr, A.G., Nesic, Z., Drewitt, G.B., et al., 2003. Ecophysiological controls on the carbon balances of three southern boreal forests. Agric. For. Meteorol. 117 (1–2), 53–71, http://dx.doi.org/10.1016/ S0168-1923(03)00023-6. Hofgaard, A., Tardif, J., Bergeron, Y., 1999. Dendroclimatic response of Picea mariana and Pinus banksiana along a latitudinal gradient in the eastern Canadian boreal forest. Can. J. For. Res. 29 (9), 1333–1346. Jarvis, P.G., Linder, S., 2000. Botany: constraints to growth of boreal forests. Nature 405 (6789), 904–905. Kang, S., Kimball, J.S., Running, S.W., 2006. Simulating effects of fire disturbance and climate change on boreal forest productivity and evapotranspiration. Sci. Total Environ. 362, 85–102. Krishnan, P., Black, T.A., Barr, A.G., Grant, N.J., Gaumont-Guay, D., Nesic, Z., 2008. Factors controlling the interannual variability in the carbon balance of a southern boreal black spruce forest. J. Geophys. Res. Atmos. 113 (D9), D09109, http://dx. doi.org/10.1029/2007jd008965. Kurz, W.A., Apps, M.J., 1999. A 70-year retrospective analysis of carbon fluxes in the canadian forest sector. Ecol. Appl. 9 (2), 526–547. Medhurst, J., Parsby, J.A.N., Linder, S., Wallin Gör Ceschia, E., Slaney, M., 2006. A whole-tree chamber system for examining tree-level physiological responses of field-grown trees to environmental variation and climate change. Plant, Cell Environ. 29 (9), 1853–1869, http://dx.doi.org/10.1111/j. 1365-3040.2006. 01553.x. Middleton, E.M., Sullivan, J.H., Bovard, B.D., Deluca, A.J., Chan, S.S., Cannon, T.A., 1997. Seasonal variability in foliar characteristics and physiology for boreal forest species at the five Saskatchewan tower sites during the 1994 Boreal Ecosystem-Atmosphere Study. J. Geophys. Res. Atmos. 102 (D24), 28831–28844 (1984–2012). Nishimura, P.H., Laroque, C.P., 2011. Observed continentality in radial growth—climate relationships in a twelve site network in western Labrador, Canada. Dendrochronologia 29 (1), 17–23. Pavlic, G., Chen, W., Fernandes, R., Cihlar, J., Price, D.T., Latifovic, R., et al., 2007. Canada-wide maps of dominant tree species from remotely sensed and ground data. Geocarto Int. 22 (3), 185–204, http://dx.doi.org/10.1080/ 10106040701201798. Piao, S., Ciais, P., Friedlingstein, P., Peylin, P., Reichstein, M., Luyssaert, S., et al., 2008. Net carbon dioxide losses of northern ecosystems in response to autumn warming. Nature 451 (7174), 49–52. Richardson, A.D., Hollinger, D.Y., Burba, G.G., Davis, K.J., Flanagan, L.B., Katul, G.G., et al., 2006. A multi-site analysis of random error in tower-based measurements of carbon and energy fluxes. Agric. For. Meteorol. 136 (1–2), 1–18. Tang, G., Beckage, B., Smith, B., Miller, P.A., 2010. Estimating potential forest NPP, biomass and their climatic sensitivity in New England using a dynamic ecosystem model. Ecosphere 1 (6), art18, http://dx.doi.org/10.1890/ES10-00087.1. Tanja, S., Berninger, F., Vesala, T., Markkanen, T., Hari, P., Mäkelä, A., et al., 2003. Air temperature triggers the recovery of evergreen boreal forest photosynthesis in spring. Global Change Biol. 9 (10), 1410–1426. Thornley, J.H.M., Cannell, M.G.R., 1996. Temperate forest responses to carbon dioxide, temperature and nitrogen: a model analysis. Plant, Cell Environ. 19 (12), 1331–1348, http://dx.doi.org/10.1111/j. 1365-3040.1996.tb00012.x. Verseghy, D.L., 1991. CLASS—a Canadian land surface scheme for GCMs. I. Soil model. Int. J. Climatol. 11, 111–113. Verseghy, D.L., 2000. The Canadian land surface scheme (CLASS): its history and future. Atmos. Ocean 38 (1), 1–13. Verseghy, D.L., McFarlane, N.A., Lazare, M., 1993. CLASS—a Canadian land surface scheme for GCMs. II. Vegetation model and coupled runs. Int. J. Climatol. 13, 347–370. Wang, Z., Grant, R.F., Arain, M.A., Bernier, P.Y., Chen, B., Chen, J.M., et al., 2013. Incorporating weather sensitivity in inventory-based estimates of boreal forest productivity: a meta-analysis of process model results. Ecol. Modell. 260 (0), 25–35, http://dx.doi.org/10.1016/j.ecolmodel.2013.03.016. Yuan, F., Arain, M.A., Alan, G.B., Black, T.A., Charles, P.A.B., Carole, C., et al., 2008. Modeling analysis of primary controls on net ecosystem productivity of seven boreal and temperate coniferous forests across a continental transect. Global Change Biol. 14 (8), 1765–1784.