Age‐dependent variation in the biophysical properties of boreal forests

The changes in boreal forest hydrology, biogeochemistry, and biophysics during succession have critical implications for the sign and magnitude of the vegetation‐climate feedbacks that might occur with a change in fire frequency, and also for the identification and attribution of changes in boreal forest to climate. We combined in situ measurements from eddy covariance sites located along an age transect in a Canadian boreal forest with spectral vegetation indices (SVIs) derived from Landsat and MODIS imagery. We found tight spatial relationships between Landsat SVIs and in situ measurements of three important biophysical properties: albedo, maximum daily uptake of CO2 (FCO2‐max), and leaf area index (LAI). The tasseled cap indices were particularly well suited for tracking biophysical variation along an age transect. Trends in brightness, greenness, and wetness from 1984 to 2005 indicated how succession drives temporal trends in biophysical properties. Albedo and FCO2‐max increased rapidly in the decade following fire and then decreased for the remainder of succession, while LAI continued to increase until ∼135 years and may decrease thereafter. The ratio of greenness to wetness indicated that photosynthesis was limited by leaf area before 10–12 years and by reduced leaf‐level photosynthetic rates thereafter, coinciding with the successional replacement of broadleaf deciduous species by evergreen conifer species. The timing of phenological events was also strongly age‐dependent, but the normalized difference vegetation index (NDVI) confounded the disappearance of the snowpack in spring for the onset of photosynthesis. Secondary succession was the dominant source of temporal variability in the biophysical properties we examined.


Introduction
[2] The biophysical and biogeochemical properties of boreal forest change markedly during secondary succession. Albedo is low immediately following fire, increases rapidly from 2 to 30 years, and then decreases again in mature forests [Amiro et al., 2006b;Randerson et al., 2006]. Rates of photosynthesis recover to maximum values 24 to 36 years following fire, and decrease thereafter [Litvak et al., 2003;Amiro et al., 2006a;Goulden et al., 2006;Welp et al., 2006]. Leaf area index (LAI) increases steadily through the first $70 years of succession and decreases thereafter [Bond-Lamberty et al., 2002b]. The length of the growing season responds to interannual variations in weather according to the fraction of evergreen vegetation in the canopy, a quantity that also varies as a function of stand age [McMillan et al., 2008].
[3] The changes in boreal forest properties during succession have critical implications for the sign and magnitude of the vegetation-climate feedbacks that might occur with a change in fire frequency , and also for the identification and attribution of changes in boreal forest to climate [Goetz et al., 2005]. Warmer temperatures are predicted to decrease the fire return interval, resulting in larger areas of young stands. Young stands have higher albedos than mature stands [Amiro et al., 2006b]. The cooling effect of the elevated albedo in young stands may exceed the combined warming effects of greenhouse gases and black carbon released by fire, and the reduced CO 2 uptake capacity of very young stands ]. An increase in fire frequency and a shift to younger stands could therefore attenuate, rather than amplify, a temperature-driven increase in forest fire frequency. However, the balance between the factors that contribute to warming (especially direct fire emissions) and those that counter warming (especially increased vegetation albedo, increased snow cover exposure, and enhanced photosynthetic rates with regeneration) is strongly dependent on an accurate determination of the trajectories of these properties during secondary succession.
[4] In situ studies along age transects provide direct, unambiguous observations of biophysical properties at specific locations, but are labor intensive and limited in their scope for spatial or temporal extrapolation. Remote sensing observations offer great scope for spatial and temporal extrapolation, but are indirect, and biophysical properties must be inferred by proxies that are often ambiguous. Several recent developments have expanded options for understanding the progression of biophysical properties during succession. First, increased computer performance (both CPU speed and disk space) and the wide availability of user-friendly image analysis software have facilitated the processing of multiple Landsat images and thus the construction of Landsat time series that extend multiple decades, and contain one or two images per year. The high spatial and spectral resolution of Landsat TM/ ETM+ imagery allows better characterization of the spectral properties of individual burn scars. The Landsat TM/ETM+ record, from 1984 until present, is sufficiently long that both decadal and successional trends can be detected. Second, the availability of the 8-day MODIS land products, such as the MOD09A1 Surface Reflectance product, opens the possibility of direct comparison with in situ observations of the seasonal variation in biophysical properties. Third, the recent collection of multiyear continuous eddy covariance measurements at multiple sites across ecological gradients allow the extraction of key biophysical properties, such as the maximum net daily uptake of CO 2 . Goulden et al. [2006] described an eddy covariance data set from a boreal chronosequence comprising 2 to 4 years of measurements from seven forests aged between 0 and 153 years (since last burn).
[5] In this study, we take advantage of these recent developments to better quantify and understand the spatial and temporal dynamics of ecosystem functioning during boreal forest succession. We combine eddy covariance data from a boreal forest age transect with remote sensing data to address two central questions. First, how do albedo, maximum daily CO 2 uptake (F CO2-max ) and LAI change during boreal forest succession? Second, how does the length of the growing season change during succession? We discuss the implications of age-dependent changes in biophysical properties for carbon and energy feedbacks in boreal forests, and for the inference of climate-driven decadal-scale trends in ecosystem function.

Site Description
[6] The study sites were within 100 km of the northern BOREAS Study Area [Sellers et al., 1997] (Figure 1) and were selected to represent a black spruce chronosequence at various stages of recovery from stand-replacing wildfire. The main features of these sites are described by Goulden et al. [2006]. The day of maximum snowmelt (t snowmelt ) in the spring was defined as the day with the most rapid decrease in snowpack depth. Ideally, the timing of snowmelt would be measured in situ, but because of the remote location of the field sites and the inherently high spatial heterogeneity of snowpack depth beneath forest canopies, we obtained snowpack depth data from a nearby weather station at Thompson Airport, Thompson, MB, which was 43 km east of UCI-1930 (http://www.climate.weatheroffice.ec.gc.ca).

Eddy Covariance Measurements
[7] Net exchanges of CO 2 , water vapor and sensible energy were measured at seven sites using tower-based eddy covariance. Each site was equipped with an eddy covariance  UCI-1850, UCI-1930, UCI-1964, UCI-1981, UCI-1989, and UCI-2003. UCI-1998 is approximately 100 km NE of UCI-1989. The geographical coordinates of the upper left corner of the image are 56°4 0 41 00 N, 99°11 0 52 00 W (northern Manitoba, Canada). tower that measured the fluxes of CO 2 , H 2 O and sensible heat, net radiation, relative humidity, air temperature and the incoming and reflected solar and photosynthetically active radiation at 30-min intervals. We calculated the maximum daily net CO 2 exchange, F CO2-max , from the 30-min data and followed a sign convention whereby positive values indicate movement of CO 2 from the atmosphere to the land surface (uptake). Incoming and reflected solar radiation were measured with thermopile pyranometers (CM3, Kipp and Zonen, Delft, the Netherlands). Daily values of albedo were calculated as the ratio of reflected radiation to incoming radiation averaged between 10AM and 2PM local time. Daily values of albedo and F CO2-max were averaged from 15 June to 15 August a period between the major phenological events when these properties were generally consistent from day to day. Other experimental details of the eddy covariance measurements are given by Goulden et al. [2006].

Total LAI and Fraction of Evergreen Needleleaf
Area ( f ENLA ) [8] The herbaceous, understory, and overstory LAIs at each site were determined in September 2004 by harvest or allometry using site-specific allometric equations and specific leaf areas [Goulden et al., 2006]. For woody species, foliage biomass was calculated using site-specific allometric equations for each species [Bond-Lamberty et al., 2002a] and leaf area index (LAI) was calculated using specific leaf area values given by Bond-Lamberty et al. [2002b]. LAI of herbaceous species was calculated from peak biomass harvests and multiplied by a specific leaf area of 11.17 ± 2.28 m 2 kg À1 obtained for Epibolium augustifolium by Bond-Lamberty et al. [2002b]. Deciduous and evergreen components were separated visually and the fraction of evergreen needleleaf area (f ENLA ) from total harvested biomass was estimated. In all cases, LAI is defined as hemisurface leaf area, or one-half total leaf area per unit ground area [Chen et al., 1997]. . Twenty-four images were selected on the basis of (1) location: 6 images from Path 34/Row 21 and 18 images from Path 33/Row 21; (2) fraction of cloud-free coverage (<50% clouds); and (3) temporal proximity to midsummer. The images covered 21 years from 22 June 1984 to 7 June 2005. The LEDAPS image was a surface reflectance image from 17 September 2001 (Path 33/Row 21) that had been corrected using the 6S atmospheric correction model [Vermote et al., 1997b] and the Dark Dense Vegetation method [Kaufman et al., 1997].
[10] The Landsat digital numbers (DN) representing radiance at sensor were linearly transformed to physical units (W m À2 steradians À1 mm À1 ) using gains and offsets provided in the image documentation. Next, we converted images to top of atmosphere reflectance (R-TOA) using formulae given by Chander and Markham [2003]. The images were stacked by band if needed, and then georeferenced to a single base image that was projected in the UTM system (Zone 14, north, Datum: WGS84). Eight or more ground control points were used to georeference the images to a common projection using first degree polynomial warping and nearest-neighbor resampling. All preliminary processing was done using the Environment for Visualizing Images (ENVI v4.0) software (RSI, Boulder, CO).  Rouse et al. [1974] EVI 2:5ÃðNIR850ÀREDÞ NIR850þ6ÃREDÀ7:5BLUEþ1 2:5ðr 4 Àr 3 Þ r 4 þ6r 3 À7:5r 1 þ1 2:5ðr 2 Àr 1 Þ r 2 þ6r 1 À7:5r 3 þ1 Huete et al. [2002] VIg GREENÀRED GREENþRED r 2 Àr 3 r 2 þr 3 r 4 Àr 1 r 4 þr 1 Gitelson et al. [2002] VARI GREENÀRED GREENþREDÀBLUE r 2 Àr 3 r 2 þr 3 Àr 1 r 4 Àr 1 r 4 þr 1 Àr 3 Brown et al. [2000] min(r 5 ) = 0.064 min(r 6 ) = 0.02 max(r 5 ) = 0.270 max(r 6 ) = 0.40 [11] The georeferenced images were then radiometrically normalized to the atmospherically corrected LEDAPS image. We used a set of 39 light and dark pseudoinvariant regions (PIR) which included lakes and unsealed roads. Each image was normalized to the base image by determining the linear relationship between the reflectances of the PIRs in the two images on a band-by-band basis [Hall et al., 1991]. For each band, PIR reflectances were always tightly correlated among images, and r 2 values exceeded 0.9 in 90% of the comparisons, giving us high confidence in the normalization procedure. A 3 by 3 square of 30-m pixels centered on the eddy covariance tower was selected for each site. Pixels contaminated with cloud were identified by high radiance in both the red and NIR bands. Sixteen bit DN exceeding 6000 in both red and NIR were rejected from subsequent analyses. Eight bit DN exceeding 120 in red and 120 in NIR were rejected from subsequent analyses. We ensured that our target regions were topographically flat from analyses of a 90-m spatial resolution Shuttle Radar Topographic Mission Digital Terrain Elevation Data. Tower sites had a slope of 0 to 1.7% and PIR targets had a slope of 0 to 2.7%.

Spectral Vegetation Indices
[12] We calculated ten ratio-based spectral vegetation indices (SVIs) from the surface reflectance of individual Landsat bands (Table 1). We also employed the tasseled cap approach described by Kauth and Thomas [1976]. Three sets of six coefficients were used to linearly transform the six spectral channels from the Landsat reflectance data to three new indices, brightness, greenness and wetness. The coefficients for brightness, greenness and wetness indices were given by Crist [1985]. While these coefficients were originally intended for the TM sensor, we assumed they were also appropriate for the ETM+, since all band reflectances were normalized to an atmospherically corrected image.

MODIS NDVI
[13] Normalized difference vegetation index (NDVI) was extracted from the 16-day, 250-m MOD13Q1 vegetation indices product and the 8-day, 500 m MOD09A1 Surface Reflectance product [Vermote et al., 1997a]. We obtained 200 Â 200 km MOD09A1 subsets distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at the U.S. Geological Survey (USGS) Center for Earth Resources Observation and Science (EROS) (http:// LPDAAC.usgs.gov). Images for all available dates coinciding with tower measurements (July 2001 to October 2005) were downloaded. We then extracted 3 by 3 pixel matrices of individual band reflectance centered on each tower. The data were screened by the MODLAND Mandatory bits (the two rightmost bits of the Quality Assurance (QA) science data set). Only values of the highest quality (zeros for the QA bits) were included in the following analysis. We also rejected MOD13Q1 NDVI on dates when the QA bits indicated snow contamination.

Detection of the Growing Season Start by Eddy
Covariance, t GSS-tower , and From MOD09A1-Derived NDVI, t GSS-NDVI [14] We defined the start of the growing season (t GSS-tower ) as the first day in which F CO2-max exceeded 35% of the annual maximum value of F CO2-max . We used MOD09A1-derived NDVI to remotely detect the start (t GSS-NDVI ) of the growing season. We found t GSS-NDVI by identifying the first day of each year that exceeded a threshold value. We required that this day be found in the first 200 days of each year and that it occur after the maximum NDWI (Table 1) was observed.
[15] We systematically searched for a single NDVI threshold that best indicated the start of the growing season across all sites and years. We used the root mean square error (RMSE), bias and dispersion as statistics to evaluate the agreement between the satellite and in situ estimates for the growing season start [Delbart et al., 2005]. The best agreement (when the RMSE was smallest) was found with a NDVI threshold of 0.58. We also attempted to find the end of the growing season by a similar approach. However, the RMSE was too high (>26 days) for predictive purposes.

In Situ Patterns of Land Surface Biophysical Properties
[16] The biophysical properties of boreal forests changed markedly during successional recovery from fire, and large variations in albedo, CO 2 exchange and LAI were observed ( Figure 2). Albedo peaked at 15 years and decreased thereafter (Figure 2a). Although albedo was not measured at the youngest stand (UCI-2003), other studies have shown that albedo is small (between 0.05 and 0.10) in the first 2 years following fire because of the presence of soot and The r 2 value for each regression is shown at the lower left of each panel. Asterisks indicate the significance level of the regression: two asterisks = P < 0.01; one asterisk = P < 0.05; no asterisks = P > 0.05. charred surfaces [Amiro et al., 2006b;Randerson et al., 2006]. Summer-F CO2-max was smallest at the youngest site, peaked at 23 years (8 years after albedo), and decreased slightly thereafter (Figure 2b). LAI was smallest at the youngest stand and continued to increase until 74 years ( Figure 2c). The fraction of evergreen leaf area increased sharply from 23% at 15 years to 85% at 23 years ( Figure 2c).

Relationships Between Landsat Spectral Vegetation Indices and Stand-Level Biophysical Properties
[17] We explored relationships between site-to-site variation in albedo, summer-F CO2-max and LAI with site-to-site variation in SVIs calculated from the Landsat TM image acquired on 6 July 2004. Short-wave albedo was best correlated with tasseled cap brightness (r 2 = 0.86) (Figure 3). Two of the ratio-based SVIs, RSR and SR, and the tasseled cap wetness, were also well correlated (r 2 > 0.60). Brightness was unique in that it depicted the progressively smaller reductions in albedo that occur as stands age beyond 15 years. The tasseled cap coefficients for brightness are all positive with the highest weighting on band 3 (Red) and 4 (NIR 850 ). The tight relationship between albedo and brightness on level terrain is unsurprising since brightness is the sum of the individual band reflectances after each has been multiplied by a positive coefficient [Zhang et al., 2002].
[18] Strong correlations were found between the SVIs and summer-F CO2-max , even though site to site differences in CO 2 uptake were not a simple monotonic function of stand age (Figure 4). The SVIs that were most strongly related to summer-F CO2-max were enhanced vegetation index (EVI) (r 2 = 0.82), NDII6 (r 2 = 0.90), NDII7 (r 2 = 0.90) and tasseled cap greenness (r 2 = 0.86). We separately tested how well these four SVIs depicted summer-F CO2-max for young (<23 years since burn) and old stands (> = 23 years since burn). All four SVIs were well correlated with summer-F CO2-max at young stands (r 2 > 0.9). EVI and greenness were also well correlated with summer-F CO2-max at old stands (r 2 > 0.9) whereas NDII6 and NDII7 were poorly correlated with summer-F CO2-max at old stands (r 2 < 0.2) (plot not shown).
[19] NDVI, VIg, VARI, NDII6, NDII7, SR and RSR all performed well (r 2 > 0.85) as proxies for LAI ( Figure 5). The tasseled cap wetness also performed well (r 2 = 0.84) (Figure 5l). Brown et al. [2000] recommended RSR (the modification of the SR with NIR 1240 ) as a proxy for LAI because of its ability to unify deciduous and evergreen species. Our analysis suggests that SR performs equally as well as RSR when considering variations across an entire chronosequence, though we included fewer sites than Brown et al. [2000].
[20] Boreal forests vary markedly with stand age in the three important biophysical parameters addressed here. Remarkably, each of the three tasseled cap transformations of Landsat surface reflectance was strongly related to one of these biophysical properties. Brightness was well correlated with albedo, greenness was well correlated with summer-F CO2-max , and wetness was well correlated with LAI (Table 2). Further work is needed to test whether these relationships are universal.

Changes in the Landsat Tasseled Cap Indices as a Function of Stand Age
[21] The tasseled cap transform was applied to the entire Landsat time series. Brightness, greenness and wetness followed orderly progressions as stands recovered from fire ( Figure 6). The time-dependent slope within each individual stand followed the overall trend with age between stands.
[22] All three indices recovered rapidly following fire. Brightness increased from a minimum value of 0.1 immediately following fire to a maximum value of 0.25 7 years after the fire (Figure 6a). Greenness increased from a minimum value of 0.01 immediately following fire to a maximum value of 0.16 after 11 years (Figure 6b). After reaching maximum values, brightness and greenness decreased for the remainder of succession. The initial rapid  [23] After reaching maximum values, the subsequent decrease in the values of brightness and greenness could be approximated with exponential decay functions fitted to all stands older than 12 years: brightness ¼ 0:1123 þ 0:1938e À0:0271t ; t > 12 years; r 2 ¼ 0:93 ð3Þ greenness ¼ 0:0733 þ 0:084e À0:024t ; t > 12 years; r 2 ¼ 0:73 [24] Wetness increased rapidly following fire, but unlike brightness and greenness, wetness continued to increase (at diminishing rates) until $135 years (Figure 6c). Up until this age, wetness was well approximated by a negative exponential function: wetness ¼ À0:1525 þ 0:1272 ð1 À e À0:033t Þ; r 2 ¼ 0:93 ð5Þ [25] These equations can be combined with the linear coefficients in Table 2 to estimate the time series of albedo, summer-F CO2-max and LAI ( Figure 6).
[26] The ratio of greenness-derived summer-F CO2-max to wetness-derived LAI represented a canopy average of net carbon assimilation per unit leaf area with units equivalent to the average leaf-level rate of CO 2 uptake (mmol m 2 leaf s À1 ) (Figure 6d). The summer-F CO2-max : LAI ratio followed an orderly decrease with age during middle to late succession. Values ranged between 4 and 10 at sites dominated by early successional species, between 1.5 and 3.6 at sites dominated by midsuccessional species (Jack Pine and Aspen at UCI-1981 andUCI-1964) and were below 1.6 for sites dominated by late successional species (Black Spruce at UCI-1930 and UCI-1850; Figure 6d).

Seasonal Patterns in MODIS NDVI: Correspondence to Measured CO 2 Exchange
[27] Remotely sensing the seasonal pattern of F CO2-max posed a different set of challenges to the sensing of the spatial pattern of biophysical properties. The Landsat time series offered high spatial resolution, but the availability of images was limited to one to three per year. We required higher temporal resolution to resolve the seasonal changes in the physical environment (e.g., snow) and the vegetation phenology. We therefore compared a time series of NDVI from the 16-day, 250-m MODIS MOD13Q1 vegetation indices product with F CO2-max measured by eddy covariance.  Table 2.
[28] The correspondence between the MOD13Q1-derived NDVI and daily F CO2-max was examined for comparatively warm (2003) and cool (2004) years at three study sites (Figure 7). NDVI captured the broad variations in the seasonal amplitude of F CO2-max between forests of different age, and in both years. Both NDVI and F CO2-max decreased in winter and increased in summer. Both NDVI and F CO2-max were comparatively small at the youngest site ( Figure 7c) and large at the oldest site (Figure 7a).
[29] However, the variations in F CO2-max that occurred at the beginning and end of the growing seasons were poorly depicted by NDVI. NDVI increased sharply following the day of maximum snowmelt (t snowmelt ), and remained elevated after the period of active photosynthesis had ended. Immediately after t snowmelt , NDVI reached between 48% and 82% of the maximum annual NDVI by t GSS-tower , whereas F CO2-max remained below 35% (Table 3). The absolute values of NDVI at t GSS-tower at the three sites during the 2 years varied between 0.504 and 0.663.
[30] At the oldest stand (Figure 7a), the interval between t snowmelt and t GSS-tower was short during a year with a warm spring (6 days in 2003) but longer during a year with a cool spring (47 days in 2004). NDVI increased sharply, from $0.22 to $0.58, following t snowmelt and then increased more gradually, from $0.58 to $0.80, following t GSS-tower .
[31] At the youngest stand (Figure 7c), the interval between t snowmelt and t GSS-tower was long during a year with a warm spring (62 days in 2003) and even longer during a year with a cool spring (92 days in 2004). NDVI increased sharply, from $À0.02 to $0.58, following t snowmelt and then increased more gradually (from $0.58 to $0.69) following t GSS-tower .
[32] Seasonal increases in NDVI were numerator-driven when the absolute rate of change in the numerator of NDVI (jd(r NIR À r red )/dtj) exceeded the absolute rate of change in the denominator (jd(r NIR + r red )/dtj). Conversely, increases in NDVI were denominator-driven when the rate of change in the denominator of NDVI exceeded the rate of change in the numerator. The first increase in NDVI, associated with snowmelt, was denominator-driven at all sites because of simultaneous decreases in NIR reflectance and red reflectance, which increased the value of the NDVI denominator relative to the numerator. The second NDVI increase, associated with the onset of photosynthesis, was also denominator-driven at the oldest stand, because of simultaneous decreases in NIR reflectance and red reflectance, but was numerator-driven at younger stands because of a decrease in red reflectance and an increase in NIR reflectance.

Age-Dependent Variations in the Start of the Growing Season (t GSS-NDVI )
[33] The consistency of NDVI values following snowmelt suggested that an absolute value of NDVI could be used as a threshold to identify the start of the growing season. To provide greater temporal resolution, we calculated a NDVI time series from the 8-day MOD09A1-derived NDVI (instead of using the 16-day MOD13Q1 NDVI product) and calculated t GSS-NDVI for each year across all sites, using an absolute threshold value of 0.58. We compared t GSS-NDVI with t GSS-tower and obtained a RMSE of 14 days. At the youngest sites, the t GSS were systematically overestimated. However, the method was only moderately biased and this simple approach captured the age-dependent trends in the start of the growing season. We were unable to reliably predict the end of the growing season using this approach. Previous studies have also encountered difficulty in the remote sensing of the end of the growing season [Delbart et al., 2005], possibly because the cues for the onset of dormancy are different to those for the onset of the growing season [Tanja et al., 2003].
[34] The day when the growing season began (t GSS ) at each of our study sites was age-dependent and increased sharply during the first 24 years of succession, but was relatively constant thereafter ( Figure 8). We retrieved NDVI from other nearby (<200 km) forests with ages determined  Figure 8. Growing season start, t GSS , detected by finding the first day when NDVI exceeded a threshold of 0.58. Large circles, triangles, and squares represent t GSS at the six main study sites: UCI-1850, UCI-1930, UCI-1964, UCI-1981, UCI-1989, and UCI-1998. Smaller, open circles (''Other Burns'') represent GS start and end days from 14 other nearby forests with known fire history between 1959 and 2005 [Stocks et al., 2002]. Plus signs (''Other Burns BD'') represent erroneous t GSS that resulted from anomalous seasonal patterns of NDVI.
from the Large Canadian Fire Database [Stocks et al., 2002]. The age-dependent trends in the starts of the growing season at these forests broadly followed the trends observed at the main tower sites. From 4 to 24 years, the start of the growing season advanced by 2 to 3 days for each year of forest age. From 24 to 155 years, the start of the growing season did not advance substantially with age. The eddy covariance data indicated that the start of the growing season was always earlier in the years with warm springs (2003 and 2005) but the NDVI-based t GSS detection captured this effect only at the younger stands.

How Might the Presence of Snow Confound the Detection of Decadal-Scale Trends in Phenology?
[35] Several analyses of the AVHRR NDVI have suggested that recent warming has led to changes in either the photosynthetic rate of boreal forest and/or an extension of the growing season [Myneni et al., 1997;Hicke et al., 2002]. These changes have major implications for the global carbon cycle; expansion of the growing season length has been correlated with changes in atmospheric CO 2 in several studies [e.g., Lucht et al., 2002]. By comparing our in situ measurements of CO 2 exchange with the seasonal course of NDVI, we found that the largest changes in NDVI during spring were a response to the melting of the snowpack rather than to phenological changes in vegetation (Figure 7). The seasonal course of NDVI must be cautiously interpreted before attributing spring increases to the onset of photosynthesis. Previous work has recognized the influence of snow cover on NDVI [Shabanov et al., 2002;Dye and Tucker, 2003], and others have hypothesized that the lack of correlation between spring NDVI and the seasonal cycle of CO 2 is related to snow cover [Angert et al., 2005]. Here, we demonstrate that spring increases in NDVI can occur up to several weeks before the spring increase in photosynthesis.
[36] In younger forests, the canopy is dominated by deciduous species, and the onset of photosynthesis must be preceded by budburst, an event that may occur several weeks later than the melting of the snowpack. In older forests, dominated by evergreen species, the onset of the growing season is triggered by increasing air temperatures during the spring [Tanja et al., 2003] and is therefore more likely to coincide with snowmelt. Therefore, confounding the NDVI-based detection of the onset of photosynthesis with the snowmelt signal will lead to a greater underestimation of t GSS in younger forests.
[37] A ''greening'' trend inferred from decadal changes in NDVI will be qualitatively correct in the sense that years with extended periods of elevated NDVI occur in years with extended growing seasons. However, in a quantitative sense, the growing season and concomitant CO 2 fluxes will be overestimated. Deciduous plants may still be dormant following the disappearance of the snowpack, but respiratory processes may have been stimulated by temperature increases that would increase, not decrease, atmospheric CO 2 concentrations. The ability to differentiate snowmelt from phenological changes in deciduous stands is particu-larly important in light of recent reports of increasing forest fire frequency resulting in changes to the age distribution of boreal forests [Flannigan et al., 2005;Kasischke and Turetsky, 2006]. Larger areas of young forests mean that the error in regional CO 2 exchange resulting from the confusion of snowmelt and phenology will increase. Previous studies have asserted that changes to the seasonal cycle of atmospheric CO 2 are caused by increased spring photosynthesis inferred from NDVI at high latitudes [Randerson et al., 1999]. As long-term (>10 years) eddy covariance data sets become available, the fidelity of NDVI, and other SVIs, as proxies for the onset of photosynthesis can be properly established. As an interim measure, applying the threshold approach described above will increase confidence that spring increases in NDVI are driven by changes in the dynamics of vegetation, rather than by changes in the dynamics of snowmelt. The lack of correspondence between CO 2 uptake and NDVI caused by melting snow has the greatest implications for hypotheses regarding trends in spring photosynthesis, rather than those regarding trends in summer or fall photosynthesis.

Remote Sensing of Biophysical Properties Along Age Transects in Boreal Forests
[38] Spatial variability in boreal landscapes was determined largely by stand age and was greatest during early succession because of a rapid increase in canopy development, followed by a transition from deciduous to evergreen vegetation ( Figure 2). The Landsat TM/ETM+ surface reflectance imagery successfully captured this variability, and certain Landsat SVIs were strongly correlated with biophysical properties measured at the stand level (Figures 3, 4, and 5).
[39] Albedo increased rapidly over the first 15 years of succession and then decreased gradually for the next $65 years. SR, RSR and brightness were all significantly correlated with albedo ( Figure 3). Brightness was the most tightly correlated with albedo, and it serves as a particularly useful proxy for the sensing of surface albedo on level terrain.
[40] Age-dependent variation in summer-F CO2-max was best predicted by EVI, NDII6, NDII7, and greenness ( Figure 4). These SVIs increased linearly together with summer-F CO2-max over the first 23 years of succession. Summer-F CO2-max decreased after 23 years. This decrease was captured by EVI and greenness but not NDII6 or NDII7. In the first 23 years of succession, NDII6 and NDII7 were driven by the increasing reflectance of NIR and increasing absorptance of SWIR. As the canopy became increasingly dominated by evergreen vegetation, NIR reflectance decreased sharply whereas SWIR absorptance continued to increase. These opposing effects caused NDII6 and NDII7 to saturate at older stands. In contrast, EVI and greenness responded linearly to rates of F CO2-max at older stands, and they serve as particularly useful SVIs for detecting spatial variations in F CO2-max across all stand ages.
[41] LAI varied by over an order of magnitude among the seven stands (Figure 2c). Eight of the 12 SVIs were significantly correlated with LAI ( Figure 5). The underlying basis for the high correlations was either due to increased leaf area absorbing more strongly in the red and SWIR and/ or reflecting more strongly in the NIR. Red absorptance is associated with chlorophyll abundance whereas SWIR absorptance is associated with plant cell water content [Nemani et al., 1993]. NIR reflectance increases with increasing leaf area because of photon scattering at air-cell interfaces in the spongy mesophyll [Woolley, 1971]. SR was ideally suited to the detection of LAI across the wide range of biomass encountered in the age transect (Figure 5h). It responded to both the increased reflectance of NIR (in the numerator) and the decreased reflectance of red (in the denominator). We found SR to be the most highly correlated to LAI, and it serves as a particularly useful SVI for detecting spatial variations across all stand ages. Of the tasseled cap SVIs, wetness was most correlated with LAI. In a single tasseled cap transform of Landsat imagery, using wetness as a proxy for LAI would strongly complement the use of brightness as a proxy for albedo and greenness as a proxy for photosynthesis.

How Do Biophysical Properties Change During Succession?
[42] Strong correlations existed between Landsat SVIs and biophysical properties among sites of different age. The Landsat time series of tasseled cap SVIs allowed us to reconstruct biophysical history over 21 years along an age transect ( Figure 6). Sources of variability in the SVI time series included age-related effects, atmospheric effects and interannual variability in climate and plant growth. The Landsat record is of sufficient length that younger stands reach the same successional stage that older stands passed through at earlier stages of the record. The similarity in the tasseled cap SVI values during the period when the ages of two stands overlapped indicated that the dominant source of the variability was stand age.
[43] Albedo is an important component of the surface energy balance [Bonan et al., 1992;Amiro et al., 2006b]. The land surface darkens immediately after fire, but then rapidly brightens, reaching maximum values 10 to 15 years after fire, and subsequently decreases. Randerson et al. [2006] found that the negative radiative forcing (cooling) due to the increased albedo of young stands was greater than positive radiative forcing (warming) due to the combined effects of emissions from the fire (black carbon and greenhouse gases) and the reduced net ecosystem production (NEP) of the early successional vegetation. Using the Landsat record of brightness allowed us to better determine the implications of age-dependent albedo on the energy budget. Amiro et al. [2006b] suggested that postmaximum albedo reduces linearly by 5% per decade. We found that the decreasing albedo follows an exponential decay function rather than a linear function (equation (3)). The rate of darkening is more rapid in the 65 years following peak albedo than in the subsequent 75 years. We propose that the period of time that stands have elevated (>0.11) albedo lasts $75 years.
[44] The initial darkening is caused by the deposition of soot and charred surfaces. The subsequent brightening over the next 12 years was due to increases in surface reflectance, resulting from changes in both the soil and the overlying vegetation. Over the course of succession, the vegetation species composition changed from predominantly deciduous to predominantly evergreen. The majority of this transition occurred between 15 and 23 years (Figure 2c). Evergreen needleleaf (ENL) species have much lower albedo than broadleaf deciduous species [Fuentes et al., 2001]. We propose that the nonlinear nature of their inclusion into the canopy is responsible for the exponential decay in albedo between 24 and $150 years. Differences in the reflectance spectra among the various ENL species are comparatively minor [Fuentes et al., 2001].
[45] Randerson et al. [2006] estimated that a reduced fire return interval would increase areas of younger stands and cause a negative rather than a positive radiative forcing. The Landsat record of albedo allowed us to interpolate accurately between stand-level measurements made along a chronosequence and further extend the range of forest ages that were sampled. We found that albedo had fallen to within 14% of full range by 65 years of succession, and was still decreasing (at an instantaneous rate) by $0.6% per year. Both the period of postfire albedo elevation and the magnitude of the elevation determine the radiative impact of these changes. Randerson et al. [2006] found that the period of albedo elevation lasted $55 years. We found that albedo was still elevated at the UCI-1930 site ($75 years). This suggests that only slight changes to the fire return interval could cause albedo-related cooling, though a quantitative analysis is beyond the scope of this paper.
[46] Researchers have hypothesized that NEP follows an orderly progression through succession [e.g., Odum, 1969]. Eddy covariance measurements of CO 2 exchange in boreal forests along age transects (chronosequence studies) have quantified this progression [Litvak et al., 2003;Goulden et al., 2006;Welp et al., 2006]. The time series of tasseled cap greenness from the Landsat imagery allowed us to infer summer-F CO2-max at stand ages not sampled by the eddy covariance towers. We observed rapid recovery of primary productivity following fire. Stand age was the dominant source of variation in summer-F CO2-max . At the sites where maximum values of summer-F CO2-max occurred, the canopy composition was in transition from predominantly deciduous broadleaf species to predominantly evergreen needleleaf species, in particular, jack pine (Pinus banksiana) (Figure 2c). At younger sites, summer-F CO2-max was limited by LAI (Figures 2 and 6c) whereas at older sites it was limited by the lower leaf level photosynthetic rates (Figure 6d). Leaf-level rates of photosynthesis in the dominant species at the most productive stands is higher than that of the dominant species at the older stands [Dang et al., 1997]. The age-dependent trajectory of biophysical properties reported here are from a single chronosequence. Under disturbance regimes where burn severity or fire return intervals are different, we would also expect differences in age-dependent trajectories. Similar analyses of Landsat tasseled cap time series indices will allow characterization of these other trajectories.
[47] Production efficiency models (PEMs), such as CA-SA, predict photosynthesis from the product of intercepted photosynthetically active radiation (IPAR) and the production efficiency, e, which specifies carbon gain per unit of IPAR [Potter et al., 1993]. LAI is the principal determinant of the fraction of incoming PAR intercepted by vegetation (fPAR) and increased sharply with age ( Figure 2). Summer-F CO2-max increased more gradually with age indicating that the production efficiency is lower at older stands. In boreal landscapes that are spatially heterogeneous because of wildfire, estimating gross photosynthesis using a PEM is complicated by the need to specify appropriate values of e for variously aged stands. The linear relationship that exists between summer-F CO2-max and SVIs, such as EVI and greenness, would preclude the need to specify e, and a value for F CO2-max could be obtained directly from remote sensing data. At nine eddy covariance sites in North America, in situ gross primary production (GPP) was more tightly correlated to EVI than to the production efficiencybased MODIS MOD17 GPP product [Rahman et al., 2005] although the relationship between EVI and GPP was stronger at deciduous stands than evergreen stands [Sims et al., 2006]. An alternative approach for the remote sensing of GPP has used EVI as a measure of fPAR on the basis of its sensitivity only to the portion of fPAR that is due photosynthetically active vegetation [Xiao et al., 2005].
[48] In many situations, summer-F CO2-max is likely to be well correlated with annual NEP. Under these conditions, brightness and greenness will act as proxies for two separate radiative forcing processes: brightness is inversely related to the absorption of incoming solar radiation and hence the net radiation available for heating the troposphere; greenness may indicate the capacity of ecosystems to sequester atmospheric CO 2 . The age-dependent trajectory of albedo and summer-F CO2-max reported here will allow a detailed analysis of the interplay between direct radiative forcing (via albedo) and indirect radiative forcing (via the photosynthetic uptake of CO 2 ). The utility of the tasseled cap SVIs as proxies for the radiative impact of disturbance deserves further exploration.

Conclusions
[49] Four implications arise from this exploration of spectral vegetation indices and in situ biophysical measurements within a boreal chronosequence: [50] 1. Each of the three tasseled cap indices derived from Landsat surface reflectance were tightly correlated with one of three important biogeophysical or biogeochemical properties.
[51] 2. The melting of the snowpack in spring confounded NDVI-based detection of the seasonal onset of photosynthesis in boreal forests.
[52] 3. Age-dependent changes in the tasseled cap indices indicated that albedo and F CO2 increase rapidly in the initial 10 -12 years of succession, and decrease thereafter, whereas LAI continues to increase until $135 years.
[53] 4. The magnitude of change in biophysical properties due to age is much greater than year-to-year variations, or decadal-scale trends.