Agricultural and Forest Meteorology Inﬂuence of stand age on the magnitude and seasonality of carbon ﬂuxes in Canadian forests

Proper management and accounting of forest carbon requires good knowledge of how disturbances and climate affect the carbon dynamics of different stand types. We have investigated such relationships by measuring, over a 5-year period (2003–2007), the net ecosystem productivity (NEP), gross ecosystem productivity (GEP) and ecosystem respiration (ER) of 26 forest sites in Canada using the eddy covariance technique. The study included black spruce, jack pine, Douglas-ﬁr, aspen, boreal mixedwood and white pine forest ecosystems ranging in age from 1- to 153-years. The dataset included six chronosequences (one afforested plantation, three harvested and two burned). Following planting, the afforested white pine stands quickly became carbon sinks and offset initial carbon losses after 4 years. Depending on forest type, the other forest stands were carbon sources for 10–18 years following a disturbance, offset initial carbon losses after 19–47 years, and showed net total gains ranging from 38 to 86 Mg C ha − 1 at 80 years. Peak NEP ranged from 0.9 to 2.9 Mg C ha − 1 year − 1 at ages of 35–55 years except for the afforested white pine where it was 6.9 Mg C ha − 1 year − 1 at 15–20 years. Stepwise regression and Pearson correlation analyses indicated that the GEP and ER of mature stands (>70 years old) were driven mainly by climate, while ﬂuxes of young stands (<19 years old) were driven by both leaf area index and climate. Although stand age of the afforested white pine plantations did not affect the GEP growing season lengths, the growing season length of the other forests increased with age until about 20 years and this coincided with the switch from carbon source to sink. With the exception of the afforested white pine, peak GEP/ER ratios of the youngest sites occurred later in the growing season compared to older sites. The strong inﬂuence of stand age on the seasonal dynamics of GEP ﬂuxes needs to be considered to avoid confounding the impacts of climate change with those of disturbance. These age-related seasonality effects are continental in scope and should be important in interpreting the time series of atmospheric CO 2 concentration measurements at regional and global scales.


Introduction
Anthropogenic greenhouse gas emissions have been increasing since the beginning of the industrial age and the resulting increase in atmospheric CO 2 concentrations are now widely believed to be changing the planet's climate. The Earth's terrestrial and ocean ecosystems have been sequestering a large portion of these CO 2 emissions (Canadell et al., 2007a;Sarmiento and Gruber, 2002). By constraining the increases in atmospheric CO 2 concentrations with the observations of a weak north-south concentration gradient, Tans et al. (1990) concluded that there was a northern extratropical net carbon sink of ∼2 Pg C year −1 . Later studies suggest that the efficiency of the terrestrial carbon sink has been decreasing due to changes in the Earth's climate (Canadell et al., 2007b).
The spatial distribution of the terrestrial sink has been the subject of significant debate over recent years (Canadell and Raupach, 2008;Ciais et al., 2010;Le Quéré et al., 2007Rödenbeck et al., 2003;Stephens et al., 2007) and the current and future role of northern forests in the terrestrial sink is still an active research question (Ciais et al., 2010). Using forest inventory data and long-term ecosystem carbon studies, Pan et al. (2011) concluded that boreal forests accounted for 21% (0.50 ± 0.08 Pg C year −1 ) of the overall global carbon sink from established forests (2.41 ± 0.42 Pg C year −1 ) between 1990 and 2007. Overall, the boreal carbon sink over this time period was found to be the sum of a reduction in Canadian forest carbon stocks due to disturbance that was offset by an increasing biomass sink in other boreal regions. Clearly ecological disturbance and forest management are having a major influence on the contribution of boreal forests to global carbon sequestration (Kurz et al., 2008a,b,c). It is also well recognized that northern forests contain large amounts of carbon (C) in both biomass and soils Apps, 1995, 1999;Tarnocai et al., 2009) and these reservoirs may be even more vulnerable to future changes in climate. Since Canada contains 10% of the world's forests, proper accounting and management of these large C stocks requires a solid scientific understanding of how disturbance and climate variability impact the emission and sequestration of carbon by these forests and how we might separate the effects of these two factors on regional C budgets.
Eddy covariance (EC) flux towers are one of the main tools used to monitor, understand and quantify CO 2 exchange between forests and the atmosphere (Baldocchi, 2008) because they provide nearcontinuous half-hourly time series of carbon, water, and energy exchange at the ecosystem scale. Regional networks of flux towers have been developed over the last two decades in different parts of the world. In Canada, the Fluxnet-Canada Research Network (2002Network ( -2007 and the follow-on Canadian Carbon Program (2007)(2008)(2009)(2010)(2011) Research Network have had a leadership role in the study of the C cycle of Canada's forests and peatlands, as affected by disturbances and climate variability . Amiro et al. (2010) presented an overview of how fire, harvest, windthrow, insects, and silvicultural treatments influence C exchange in North American forest ecosystems.
In this study, we focus on a continental-scale transect of forest flux towers across Canada to address how disturbance types, stand age and climate variability may impact the C exchange of northern and coastal forests. We use six different chronosequences to determine: (1) the age at which forests reach the C compensation point when they switch from a C source to a C sink, and the offset point when carbon sequestration equals carbon loss following disturbance or afforestation, (2) the influence of stand age on the seasonal dynamics of the fluxes, and (3) the extent to which climate versus stand structure are the main drivers of interannual variability following disturbance.

Fluxnet-Canada and the Canadian Carbon Program
The Fluxnet-Canada Research Network (FCRN) was established in 2002 and consisted of a series of eddy-covariance flux towers installed in mature and disturbed forests as well as peatlands across the country. Continuous year-round measurements of carbon, water and energy exchanges between the land surface and the atmosphere were made using standard measurement protocols. The FCRN was the next logical step in the study of the interactions between Canadian forests and the atmosphere, following the Boreal Ecosystem-Atmosphere Study (BOREAS) which was carried out between 1994 and 1996 (Sellers et al., 1997). The FCRN was established using the main flux tower sites from BOREAS in Saskatchewan (which had become the Boreal Ecosystems Research and Monitoring Sites (BERMS) program) and Manitoba, as well as a few other existing sites, to which new sites were added so as to produce an east-west transect across Canada. Several chronosequences were included in the transect, facilitating the study of disturbance effects and stand age on the carbon cycle. The FCRN was succeeded in 2007 by the Canadian Carbon Program (CCP).

Study sites
The description and characteristics of the 26 forests sites, ranging in age from 1-to 153-years old, used in this study are given in Tables 1 and 2. Sixteen of these sites were supported by the FCRN/CCP Margolis et al., 2006), four (WP39-ON, WP74-ON, WP89-ON, WP02-ON) were FCRN/CCP associated sites (Peichl and Arain, 2006), and six (BS1850-MB, BS30-MB, BS64-MB, BS81-MB, BS89-MB, BS98-MB) were part of the Ameriflux network (Goulden et al., 2006). The study, which included two burn chronosequences (BBS-MB, BJP-SK), three harvested chronosequences (HBS-QC, HDF-BC, HJP-SK) and an afforested plantation chronosequence (PWP-ON), was comprised of 10 black spruce (BS), seven jack pine (JP), three Douglas-fir (DF) and four white pine stands (WP), and one each of aspen (ASP) and a mixedwood stand (MW). Sites within a chronosequence were similar in climate, soil characteristics, disturbance history and topography (Goulden et al., 2006). Measurements were available for periods ranging from 2.5 to 5 years, with 50% of the sites having a full 5 year measurement record.

Measurements and data analysis
Meteorological and flux data for 105 site-years were obtained from either the FCRN Data Information System (http://fluxnet.ccrp.ec.gc.ca) or directly from Principal Investigators (PIs). Half-hourly C fluxes were measured using the EC technique (Aubinet et al., 2000;Morgenstern et al., 2004) and, with the exception of the BBS-MB chronosequence, were processed following FCRN protocols Zha et al., 2009). Net ecosystem exchange (NEE) was computed as NEE = F c + S c , where F c is the measured C flux and S c is the rate of change in CO 2 storage between the measurement height and the ground. Depending on the site, F c was measured using either a closedor open-path infrared gas analyzer (IRGA). S c was computed using a multi-level CO 2 concentration measurement system when available, or as a one-level storage term calculated using the CO 2 concentration measured by the EC IRGA. Initial quality control of all half-hourly data was carried out by each site's PI. With the exception of the BBS-MB chronosequenc, cross-validation of meteorological and eddy-covariance equipment was carried out as part of FCRN quality assurance activities using a standard set of roving flux and meteorological calibration equipment.
Night-time measurements were eliminated when atmospheric conditions were calm. These conditions were identified when the measured friction velocity (u * , m s −1 ) was lower than the sitespecific threshold (u *th ). The threshold was estimated for each site by plotting nighttime NEE against u * and selecting the minimum u * value (aggregated in 0.05 m s −1 bins) when NEE was no longer dependent on u * (Humphreys et al., 2006). Moreover, NEE outliers (30-min fluxes farther than four standard deviations away from the monthly mean) were eliminated except for the SK-HJP02 site where most of the outliers were retained because they are believed to be caused by atmospheric gravity waves. Net ecosystem productivity (NEP), which was calculated as NEP = −NEE, is positive when CO 2 is sequestered by the ecosystem (sink) and negative when CO 2 is emitted to the atmosphere (source). Data gaps in CO 2 fluxes were filled and ecosystem respiration (ER) and gross ecosystem productivity (GEP) were estimated following the standard FCRN gap-filling procedure (Barr et al., 2004;Amiro et al., 2006) for all sites with the exception of WP74-ON, WP89-ON, WP02-ON, F77-SK, F89-SK and F98-SK. The FCRN procedure uses a three-parameter logistic equation relating ER to soil temperature, the Michaelis-Menten equation to relate GEP to PPFD and a 100-point moving window to allow certain parameters (e.g. A max ) to vary over time. The gap-filling procedure used for WP74, WP89, WP02 sites is described in Peichl et al. (2010a,b). The use of open-path analysers at the F77, F89 and F98 sites posed challenges during cold weather. Although there is a known heating issue with these analysers (Burba et al., 2008), Amiro (2010) demonstrated that the heating corrections were not a satisfactory solution for this data set, and that gap-filling winter data with respiration functions was the best solution. Open-path analysers were also used at the HBS00-QC and the HJP02-SK (for July 2003-November 2004 period only) but sensor heating issues were not observed at either of these cutover sites.
The photosynthetic growing season was identified using nongapfilled GEP, except for the WP74-ON, WP89-ON and WP02-ON sites where gapfilled GEP was used instead because of the presence of large data gaps (several weeks in length). Five-day running means of GEP were calculated for each data set and a yearly maximum was identified for each site. The beginning of each growing season was identified as the first occurrence of five consecutive days with a GEP running mean greater than 15% of the yearly daily maximum GEP. The end of the growing season was identified as the first occurrence when the GEP running mean fell below 15% of the maximum value for five consecutive days during the calendar year. This method is similar to that used by McMillan et al. (2008) and Zha et al. (2009). Using non-gapfilled GEP data was the preferred method in our case so as to avoid any influence of the gapfilling algorithm which set GEP to 0 when site temperatures fell below specific thresholds and we used 5-day running means to ensure that each half-hour of the day was represented in the calculation of mean daily GEP. The 15% threshold was selected after analyzing the results obtained using various thresholds (5-20%) and best identified the period of continuous photosynthesis at all sites as well as spring ramp-up, and corresponded well with the end of the snow cover period.
To identify periods with snow cover, we first calculated the mean value of the ratio of upward to downward photosynthetic photon flux density (PPFD) for each day using only daytime data. A Table 2 Summary description of flux sites used in this study.

TSD/age (years) a
Ta  threshold specific to each site was identified above which the day was considered to have snow on the ground. Snow cover periods were identified when there were at least five consecutive days with the ratio greater than the threshold. When multiple periods existed within a calendar year (e.g. at the beginning and end of the snow season), the periods were merged together if they were separated by less than 15 days. When no upward PPFD data were available, we used near-surface soil temperature measurements. When possible, the snow cover periods identified using these criteria were validated using other data such as snow depth. A robust nonlinear regression analysis conducted in MATLAB (The MathWorks Inc., Natick, MA) was used to predict annual net ecosystem production (NEP p = a * exp (b * age) + c * exp (d * age) ), C losses and net gains at age 80 years, and the carbon compensation (C pt ) and offset (C off ) points for specific chronosequences and ecosystems. The BS-QC chronosequence was not used for this analysis because it was composed of only two sites. Following Dragoni et al. (2007), the uncertainties for the total C loss (cumulative NEP p from age 1 to C pt ) and net C gain (cumulative NEP p at age 80 years) estimates were calculated using a Monte Carlo simulation. An estimate of the random error (ε s ) was generated using a normal distribution with a mean of 0 and a standard deviation equal to the standard deviation of each NEP p as generated by the regression model. Simulated values of NEP were then calculated as NEP s = NEP p + ε s . The simulation was repeated 10,000 times and cumulative C losses and gains were then calculated for each simulation and the uncertainty of the net C loss and net C gains was estimated by calculating the standard deviation of the simulated summations.
Forward stepwise linear regressions (using p = 0.05 for entry into the model) and Pearson correlations performed by SAS (SAS Institute, Cary, NC) were used to identify variables having the greatest effect on annual GEP, ER and NEP. All other regressions were performed using SigmaPlot (Systat Software, Inc., Chicago, IL).

Results and discussion
3.1. Source/sink transition age and carbon offset point Fig. 1a shows the composite picture of annual NEP as a function of age since disturbance or afforestation, for all stands in this study. Estimates of NEP obtained from exponential curves fit to the data from subsets of this full dataset indicate that the number of years needed to reach the C compensation point, the amount of carbon lost during this period and the time it takes to offset initial C losses, all depend on the stand type (Table 3). This curve-fitting exercise indicates that, generally, boreal forests in Canada are C sources for the first 9 years of growth and offset initial losses at Table 3 Results from nonlinear regression analyses of annual NEP versus age since disturbance data in Fig. 1 Total C lost as predicted by the model (cumulative NEP from age 0 to carbon compensation point).
f Age when C accumulation is predicted to offset initial carbon loss (age when cumulative NEP = 0). g Age at which maximum NEP occurs, and maximum NEP as predicted by the model.  Table 1 for site definitions.
26 years of age. An analysis of specific chronosequences indicates that the burned black spruce (BBS-MB), harvested jack pine (HJP-SK) and harvested Douglas-fir chronosequences (HDF-BC) become C sinks 10, 14 and 18 years after disturbance, respectively, and offset initial C losses caused by the disturbance after 19, 34 and 47 years of regrowth, respectively (Table 3). Giasson et al. (2006) suggested that the speed of post-disturbance recovery is controlled by the speed of GEP recovery against a background level of heterotrophic respiration, which seems to be confirmed by our results (Fig. 1). Using various methods, previous studies of these same forest chronosequences have found that they would become C sinks at ages ranging from 11 to 20 years (Goulden et al., 2011;Grant et al., 2007a,b;Schwalm et al., 2007;Zha et al., 2009). For a set of chronosequences located throughout North America and originating from various disturbance types, Amiro et al. (2010) found that disturbed ecosystems became C sinks at between 10 and 20 years of age. Needle-leaf forests in the U.S have been predicted to become sinks at ages ranging from 4 to 16 years (Thornton et al., 2002), while Kolari et al. (2004) reported that a Scots pine chronosequence transitioned from a source to a sink at 12 years of age. In comparison to the boreal and coastal forest sites, the afforested white pine plantations (PWP-ON) could be C sources for the first 2 years after planting and could offset initial losses the following year (Table 3). The much shorter time to reach the C compensation point and the higher NEP rates for the  Table 1 for site definitions. PWP-ON chronosequence compared to other chronosequences in this study are likely due to the low residual soil carbon of the abandoned agricultural land on which the PWP-ON chronosequence was established (Peichl et al., 2010a), which would have resulted in significantly lower respiration rates. In particular, this low carbon stock legacy enabled these plantations to become C sinks after only 2 years, in contrast to the very productive Douglas-fir stands which took 18 years to become sinks. Clark et al. (2004) reported similarly high rates of NEP for a chronosequence of slash pine plantations established on previous plantations (with tree stumps present) in Florida, although the youngest stand in the chronosequence was a strong C source.
According to our regression model, the NEP of the HJP-SK chronosequence would peak at 0.98 ± 0.30 Mg C ha −1 year −1 in the 38-43 year age range, whereas the NEP of BBS-MB and HDF-BC would peak at 1.49 ± 0.40 and 2.92 ± 1.0 Mg C ha −1 year −1 in the age ranges of 53-68 and 46-49 years, respectively (Table 3). The NEP of the afforested WP-ON chronosequence, on the other hand, would peak at rates of 6.94 ± 0.66 Mg C ha −1 year −1 at 16-18 years of age. Furthermore, peak NEP would coincide with peak leaf area index (Fig. 2). These results are comparable with those of Pregitzer and Euskirchen (2004) who reported that the peak NEP of the boreal forest biome occurred in the age range of 31-120 years and Zha et al. (2009) who found that the HJP-SK chronosequence would reach peak NEP values at 50 years of age. Thornton et al. (2002), on the other hand reported that evergreen needleleaf forests (harvest and burn origin) in the U.S. could reach peak NEP values much earlier (8-19 years of age) and only 2-7 years after becoming sinks.
Our analysis indicates that boreal forest stands lose between 3.7 ± 1.8 and 12.2 ± 0.8 Mg C ha −1 during their initial stages of growth after harvesting or fire and before switching to C sinks, while the temperate HDF-BC stands lose 60.3 ± 2.5 Mg C ha −1 (Table 3). Initial losses by the afforested WP-ON stands are only on the order of 2.1 ± 1.5 Mg C ha −1 for reasons explained above. In comparison, Thornton et al. (2002) predicted that the evergreen needle leaf forests located in the warmer climates included in their study typically would lose between 14.9 and 30.5 Mg C ha −1 , and as much as, 85 Mg C ha −1 . It is important to note that our figures do not include the initial C losses due to combustion and wood harvesting which can be quite variable depending on fire intensity and stand density. For instance, it has been estimated that the average amount of carbon lost due to fires in the region containing the BBS-MB chronosequence is 12.4 Mg ha −1 (Bond-Lamberty et al., 2004), which, according to our analysis, would take an additional 13 years to offset and reduce the net C gain of this chronosequence from 81.3 ± 3.4 Mg C ha −1 at 80 years of age to 68.9 Mg C ha −1 . Canada-wide estimates of C losses from fire combustion are about 8-20 Mg C ha −1 , depending on ecoregion, so larger combustion losses from some fires would likely take a longer period for NEP recovery.
We estimate that at 80 years of age, forest stands in this study would have net C gains ranging from 38.1 ± 2.1 to 86.3 ± 12.1 Mg C ha −1 whereas the afforested white pine plantations could have gained as much as 270 ± 6.5 Mg C ha −1 at the same age (Table 3). Furthermore, annual NEP of mature forest stands at 80 years of age would range from 0.5 ± 0.2 to 2.0 ± 3.3 Mg C ha −1 and decrease to values ranging from 0.3 ± 0.2 to 1.4 ± 3.3 Mg C ha −1 (results not shown) at 100 years of age. We recognize that the projections for net C gain at 80 years of age for the HDF-BC and PWP-ON chronosequences may be somewhat uncertain, because of a lack of data points in that age range and beyond. Trofymow et al. (2008) report a total ecosystem C content estimate of 750 Mg C ha −1 for an old-growth Douglas-fir landscape, an amount that can only be reached through sustained net C uptake. Schwalm et al. (2007) predicted a net C gain for the HDF-BC chronosequence of 21.0 Mg ha −1 by 56 years of age, whereas Grant et al. (2007b) used a process model, Ecosys, to estimate an annual NEP of 2-4 Mg C ha −1 year −1 at 50 years of age for the same chronosequence. For comparison, our model predicts a net C gain of 27.8 Mg ha −1 at 56 years and an annual NEP of 2.9 Mg C ha −1 year −1 at 50 years of age for this ecosystem. Furthermore, fertilization of the HDF-BC chronosequence in early 2007 increased the NEP of the DF49-BC and HDF88-BC by 1.68 and 1.82 Mg C ha −1 year −1 , respectively, while the NEP of the youngest site (HDF00-BC) was likely reduced by 0.03 Mg C ha −1 year −1 for the last year of data in our analysis (Jassal et al., 2010).

Stand age and the seasonality of carbon fluxes
GEP season lengths (SL GEP ), the contiguous period when GEP was at least 15% of the maximum annual rate, varied with respect to age, plant functional type and climate/geographic region (Table 4). Generally, SL GEP of the forested stands increased with stand age. The deciduous OA-SK stand had considerably shorter growing seasons compared to coniferous stands from the same region, while the Douglas-fir stand and white pine plantations located in temperate regions had longer growing seasons than stands of similar age in other regions. However, stand age did not have an effect on the SL GEP of the afforested white pine plantations (Table 4).
For a 1-year long period of measurements in a chronosequence study of burned stands in Alaska, USA, Welp et al. (2006) reported that the onset of the growing season in a 15-year old stand was delayed by 3 weeks compared to an 80-year old stand and the length of the carbon uptake period (CUP) was reduced at 3and 15-year old stands. Previous studies using some of the same chronosequences included in this analysis have reported that, over a 2-to 3-year period, recently burned and harvested black spruce stands (BS98-MB, HBS00-QC) had shorter growing seasons with later onsets and earlier endings (Bergeron et al., 2008;McMillan et al., 2008) compared to a mature stand, and a seven to 11-year Table 4 Day of year start and end dates for GEP growing seasons (GSGEP) and GEP growing season lengths (SLGEP: mean daily GEP > 15% of maximum mean daily GEP for the year) for 2003-2007. Data within parentheses are standard deviations for the means of all years. old harvested jack pine stand (HJP94-SK) tended to have a shorter growing season and later onset ) than a mature stand.
The methods used to identify the growing season in these previous studies varied. Consequently, to better determine the influence of disturbance on SL GEP and the timing of the onset and end of the growing season, we analyzed data from the five forest chronosequences in our data set. We did not include the white pine plantation chronosequence in the analysis because its growing season length did not vary with stand age. The oldest stand in each chronosequence was compared to the younger stands in the same chronosequence. The differences in SL GEP and the timing of the onset and end of the growing season between the oldest stand and younger stands in the chronosequence decreased with increasing age of the younger stands and became minimal when stands reached approximately 15-20 years of age (Fig. 3a,  b and Table 4). The youngest stand in each chronosequence had GEP growing seasons that were on average between 9 and 58 days shorter, and as much as 109 days shorter for the HDF00-BC site, than the oldest stand in the same chronosequence ( Fig. 3a and Table 4). These shorter growing seasons are attributable to either delays in the start of photosynthesis in the spring (Fig. 3b) or an earlier end to the growing season, or both (Table 4), but the effect of age appears to be more pronounced in the spring than in the autumn.
The shorter GEP growing seasons exhibited by the younger disturbed stands likely result from the transition from the evergreen conifer dominated canopies of older stands to the deciduous and semi-persistent shrubs and herbaceous plants that dominate the canopies of recently disturbed stands (Goulden et al., 2006;Welp et al., 2006). The age at which there are no longer any differences in SL GEP likely indicates the transition back to an evergreendominated canopy, which coincides with the switch from C source to C sink for the chronosequences in this study, and very likely, with canopy closure. Using the Ecosys model, Grant et al. (2010) reported that pioneer deciduous plant functional types were mostly responsible for GEP during the first 10 years after clearcutting at the HBS-QC, HJP-SK and HDF-BC chronosequences while the dominant coniferous species were mostly responsible for GEP thereafter. The absence of an age effect on SL GEP for the white pine plantation chronosequence (PWP-ON) is probably because of the very different ground cover properties and dynamics of these intensively managed sites.
Differences in annual NEP among stands in the same chronosequence tended to increase as differences in GEP season length increased (Fig. 3c), with annual NEP decreasing by approximately 40 g C m −2 year −1 with each 10 day decrease in season length for the Douglas-fir chronosequence and 25 g C m −2 year −1 for the jack pine and black spruce chronosequences. These rates are similar to carbon uptake rates ranging from 2 to 4 g C m −2 d −1 reported for coniferous forests in Europe and North America (Churkina et al., 2005) and for mature coniferous and deciduous forests in eastern North America (Richardson et al., 2009).  3. (a and b) Differences in GEP season length and number of days to the start of GEP growing season, respectively, between the oldest stand in the chronosequence and each of the younger stands in the same chronosequence. The x-axis indicates the age of the younger stand of the pair. (c) Differences in annual NEP between the oldest stand in the chronosequence and each of the younger stands in the same chronosequence, with respect to the difference in season length for the same pair. The afforested white pine plantation chronosequence (PWP-ON) was not included in the analysis because we did not detect an effect of stand age on season length. See Table 1 for site definitions.
The carbon balance of northern ecosystems is sensitive to changes in the timing of spring and autumn (Goulden et al., 1998;Piao et al., 2008;Randerson et al., 1999). Atmospheric CO 2 measurements have shown that warmer and wetter springs are associated with greater ecosystem CO 2 uptake in boreal North America (Goetz et al., 2007;Randerson et al., 1999). Furthermore, satellite-based analyses have examined the effects of increasing temperatures in northern latitudes on vegetation structure and phenology with, at times, conflicting results (Goetz et al., 2005;Myneni et al., 1997;Wang et al., 2011;Zhang et al., 2007). Our results suggest that the influence of stand age on the timing of ecosystem fluxes is continental in scope and needs to be considered when predicting the C balance of forests. This is supported by results from Welp et al. (2006) who reported an increase in the seasonal amplitude of CO 2 over a 15-year old burned stand compared to an 80-year old  Table 1 for  site definitions. stand, along with a delayed onset of the growing season. Furthermore, Zimov et al. (1999) reported that disturbed Siberian forest stands had greater seasonal amplitudes of CO 2 compared to undisturbed forest stands, these greater amplitudes being accompanied by shorter CUPs, and these differences among stands were linked to shifts from coniferous species to grasses, herbs and deciduous species.
With the exception of the burned jack pine (BJP-SK) and plantation white pine chronosequences (PWP-ON), the delay in the springtime start of the GEP growing season for the youngest stands leads them to reach maximum monthly GEP/ER ratios 1-3 months later compared to older stands (Fig. 4a-d). The mature deciduous aspen stand (ASP-SK) showed a similar pattern, where GEP/ER ratios peaked in June or July, 1-2 months after the mature coniferous stands (Fig. 4e). These results indicate that recently harvested stands can only attain maximum sink or minimum source strength once their deciduous dominated canopy is fully developed.

Controls of annual GEP, NEP and ER
Given the strong relationship between C fluxes and stand age, stands were divided into three age classes according to time since disturbance/age (mature: 73-153 years, intermediateaged: 23-64 years and young: 1-15 years) to identify the meteorological and/or structural variables having the most significant effects on annual GEP, ER and NEP. A stepwise regression analysis indicated that the combination of mean growing season volumetric soil water content (SWC-GS) and near-surface soil temperature (T s ) explained 85% and 86% of the variability in annual GEP and ER, respectively, of mature stands, and growing season mean air temperature explained 44% of the variation in annual NEP (Table 5). Total above-ground biomass of mature stands was also strongly and positively correlated to all three C fluxes (Tables 5 and 6). Interestingly, SL GEP of mature stands was strongly and negatively correlated to annual GEP and ER (Table 6). Although negative correlations between growing season length and NEP have been reported for sub-alpine forests and linked to decreased water availability due to an earlier snow-melt (Hu et al., 2010;Sacks et al., 2007), we believe the negative correlations we found between SL GEP and GEP and ER were caused by the presence of different plant functional types in the analysis of mature stands. The mature deciduous aspen (ASP-SK) and mixedwood (MW-ON) stands tended to have higher annual GEP and ER rates (Fig. 1b) and shorter growing seasons (Table 4) compared to the boreal mature stands.
Previous multi-site regional and global analyses of climatic controls on C fluxes have yielded somewhat similar results to those we found for mature boreal stands in Canada. Using global datasets, Luyssaert et al. (2007) reported that 71% of the variability in annual GEP was explained by air temperature and precipitation, whereas Law et al. (2002) reported that 64% of the variation in GEP was explained by air temperature and site water balance. Janssens et al. (2001) and Law et al. (2002) both found poor correlations between annual ER and mean annual air temperatures. Reichstein et al. (2007), on the other hand, reported that annual ER and GEP of European sites located above 52 • N were correlated with mean annual air temperature and sites located below 52 • N were correlated with available soil moisture. Finally, annual NEP has been reported by Yi et al. (2010) to be strongly correlated with mean annual air temperature at mid-and high-latitude sites, whereas Law et al. (2002) and Piao et al. (2008) reported only weak correlations between annual NEP and water balance  and mean annual air temperature.
The regression and correlation analyses indicate that annual C fluxes of the intermediate-aged sites are related to both climate and stand structure. SL GEP (partial R 2 = 0.91 and 0.89) combined with SWC-GS explained 95% and 94% of the variation in the annual GEP and ER of these stands, while 74% of annual NEP was explained by a combination of the number of days the soil was frozen (T s0 ) and T s (Table 5). LAI and total rainfall during the GEP growing season (RAIN-GS) were also strongly correlated to annual GEP of the intermediate-aged stands while LAI, temperature (soil and/or air), snow cover and rainfall (ER only) were also strongly correlated to annual ER and NEP (Table 6). We believe that the strong correlation between LAI and the C fluxes actually reflect regional climate differences among the stands within the intermediate-aged group. The DF49-BC and WP39-ON located in temperate regions had larger fluxes (Fig. 1) and greater LAI (Fig. 2a) compared to the other stands in this age class and LAI was also strongly correlated to both annual soil and air temperatures (r = 0.81 and 0.78, respectively, results not shown). This would indicate that, as with mature stands, the fluxes of the intermediate-aged stands in this study are mainly affected by climate through SL GEP and air and soil temperatures.
Our results indicate that the combination of LAI (partial R 2 = 0.87) and SL GEP explained 90% of the variability in the annual GEP of young, recently disturbed stands (Table 5). The very strong positive correlation between LAI and GEP confirms that stand Annual NEP -young f ns ns ns ns ns ns ns ns ns ns n/a ns a Ta: mean annual air temperature; Ts: mean annual near-surface soil temperature; Ta-GS: mean GEP growing season air temperature; Ts-GS: mean GEP growing season near surface soil temperature; Ts0: number of days when Ts < 0 • C on annual basis; SNOW: number of days with snow cover on annual basis; RAIN-GS: total growing season rainfall; SWC-GS: mean growing season soil water content; PPFD-GS: growing season mean PPFD; SLGEP: GEP growing season length; AGB: total above-ground biomass; LAI: leaf area index. b n/a: r not calculated for the pair of variables, no r was calculated between AGB and carbon fluxes of intermediate-aged and young recently disturbed stands because there was only a one-time AGB measurement available for each stand and AGB varies as stands age to maturity. c ns: r not significant at the 0.05 level. d BS1850-MB not included in analysis because at least on meteorological variable was not available. e WP74-ON not included in analysis because at least one meteorological variable was not available. f WP89-ON not included in analysis because at least one meteorological variable was not available.
structure is more important than climate in determining GEP fluxes of young stands in the years immediately following a disturbance. Annual ER, on the other hand, was correlated with both climate and stand structure, with the length of the snow-cover period, LAI and soil and air temperatures combining to explain 90% of the variation in ER from these stands (Table 5). No significant regressions or correlations were found for annual NEP of the young stands which is a similar result to Law et al. (2002) and Piao et al. (2008), who found only weak correlations between annual NEP and mean annual air temperature.

Conclusion
Coordinated measurements of ecosystem fluxes are important for understanding the influence of disturbance and climate variability on the C cycle of Canadian forests at the continental scale. Our results have shown that, depending on species/ecosystem type, harvested and burned forest stands are sources of carbon for the first 9-17 years following the disturbance and may lose between 3.7 ± 1.8 Mg C ha −1 and 60.3 ± 2.5 Mg C ha −1 during this period. These initial losses would be offset when the stands had reached 19-47 years of age (excluding the initial C lost through either combustion or harvest removals) and stands would attain peak NEP values of 0.89 ± 0.19-2.92 ± 1.0 Mg C ha −1 at ages ranging from 35 to 55 years. We also estimate that forest stands in Canada achieve net total gains of 38-86 Mg C ha −1 at 80 years of age and are carbon neutral or weak sinks at 100 years of age. Forest plantations established on abandoned agricultural land in the temperate zone, on the other hand, only lose 2.1 ± 1.5 Mg C ha −1 over the initial few years following establishment and then become a C sink. Overall, these data should help to more accurately determine the influences of forest management on the Canadian C budget by serving as data benchmarks for C cycle models.
There is little information currently available in the literature concerning the effects of disturbance on the onset, end, and duration of the growing season. Our results indicate that, with the exception of the afforested plantation stands, the youngest stands in each chronosequence had GEP growing seasons that were between 9 and 109 days shorter than the oldest stand in the same chronosequence. The differences in GEP season length decreased with increasing stand age, becoming minimal at 15-20 years of age and coinciding with the transition from C source to sink. The shorter growing seasons of the younger stands were associated with changes in the seasonal pattern of the C fluxes. Thus, disturbance influences not only stand structure and species composition, but also the seasonality of the GEP and NEP fluxes. Atmospheric CO 2 concentrations in the northern hemisphere also have a very strong seasonality that has a large influence on global patterns of atmospheric CO 2 . The changes in the seasonality of ecosystem fluxes following disturbance across the range of ecosystems described in this study improve our understanding of how ecological disturbance might influence the time series of atmospheric CO 2 concentrations. This could be important for separating the influences of climate variability versus disturbance at larger spatial scales.
Our results indicate that climate (soil moisture, air and soil temperatures and growing season length) is the primary factor affecting annual GEP, ER and NEP of mature stands and stand structure, through LAI, is the main factor affecting the annual GEP of young, recently disturbed stands. These differences in C flux drivers need to be considered when we try to detect a climate change signal that is mixed into a disturbance signal. While annual ER of young, recently disturbed stands is correlated with both climate and LAI, none of the variables tested were found to be related to annual NEP of these stands.