We use eddy covariance measurements of net ecosystem productivity (NEP) from 21 FLUXNET sites (153 site-years of data) to investigate relationships between phenology and productivity (in terms of both NEP and gross ecosystem photosynthesis, GEP) in temperate and boreal forests. Results are used to evaluate the plausibility of four different conceptual models. Phenological indicators were derived from the eddy covariance time series, and from remote sensing and models. We examine spatial patterns (across sites) and temporal patterns (across years); an important conclusion is that it is likely that neither of these accurately represents how productivity will respond to future phenological shifts resulting from ongoing climate change. In spring and autumn, increased GEP resulting from an ‘extra’ day tends to be offset by concurrent, but smaller, increases in ecosystem respiration, and thus the effect on NEP is still positive. Spring productivity anomalies appear to have carry-over effects that translate to productivity anomalies in the following autumn, but it is not clear that these result directly from phenological anomalies. Finally, the productivity of evergreen needleleaf forests is less sensitive to phenology than is productivity of deciduous broadleaf forests. This has implications for how climate change may drive shifts in competition within mixed-species stands.
Numerous studies have documented the effects of recent climate change on the phenology of plant and animal species across a wide range of taxa (Peñuelas et al. 2002; Badeck et al. 2004; Schwartz et al. 2006; Cleland et al. 2007; Parmesan 2007). There is also increasing interest in how these phenological shifts may affect ecosystem processes and services (Gu et al. 2003; Morisette et al. 2009; Noormets 2009). Related to this are important questions about how phenology mediates feedbacks of terrestrial ecosystems to the climate system (Schwartz 1992; Moore et al. 1996; Hogg et al. 2000), as for example through seasonal variation in albedo, partitioning of available energy to latent and sensible heat fluxes, and vegetation-atmosphere exchange of carbon dioxide and other greenhouse gases (Hayden 1998; Pielke et al. 1998; Chapin et al. 2000). However, there are still large uncertainties about the magnitude of these feedbacks, and how they may be influenced by ongoing climate change (Peñuelas et al. 2009).
It is generally assumed that warming will increase the length of the growing season (or ‘active season’) in boreal and temperate ecosystems, where wintertime temperatures are metabolically limiting. Increased primary productivity is expected to result simply from the fact that more days are available for carbon assimilation and biomass growth. Indeed, two previous modelling studies have presented similar estimates of the slope of the relationship between growing season length and net productivity in deciduous forests (5.2 g C m−2 d−1, Lieth 1975; 5.9 g C m−2 d−1, Baldocchi & Wilson 2001). Empirical analysis of spatial patterns across sites (Baldocchi et al. 2001; Churkina et al. 2005; Baldocchi 2008) and interannual variation at individual sites (e.g. Chen et al. 1999; Black et al. 2000; Carrara et al. 2003; Barr et al. 2004, 2009; Hollinger et al. 2004; Niemand et al. 2005; Richardson et al. 2009a) also give support for this standard model. However, other studies have indicated more surprising and somewhat contradictory results. For example, Dunn et al. (2007) could not discern a significant relationship between growing season length and net ecosystem productivity (NEP) in a boreal conifer stand, likely because increased gross ecosystem photosynthesis (GEP) in years with a longer growing season was offset by concurrent increases in ecosystem respiration (ER). This can be attributed to coupling between photosynthesis and autotrophic respiration, and the sensitivity of both photosynthesis and respiration to temperature. Sacks et al. (2007) found that annual NEP of a subalpine conifer forest decreased in years with an early spring because earlier onset resulted from a shallow snowpack, which translated to reduced soil water content, and ultimately drought, later in the growing season (see also Hu et al. 2010). More recently, Piao et al. (2008) reported that while recent autumn warming tends to prolong the growing season in higher latitude ecosystems, the resulting late season respiratory losses far exceed the gains in GEP at that time of year, reducing the impact of concurrent early spring gains on net carbon sequestration. Note, however, that carbon assimilation of many boreal coniferous forest ecosystems is light-limited in the autumn (e.g. Suni et al. 2003) and therefore warmer autumn temperatures may only marginally increase photosynthesis. In fact, in autumn and winter at high latitudes, temperature and solar radiation tend to be negatively correlated, and warming would likely bring more cloudy weather and less radiation (Vesala et al. 2009). By comparison, in spring, air temperature triggers the onset of photosynthetic uptake, and by the time this occurs, solar radiation is typically not limiting (Tanja et al. 2003).
Both indirect and lagged effects of ‘early spring’ or ‘late autumn’ on productivity can be envisioned. In addition to triggering earlier spring onset, warmer temperatures might increase N mineralization rates, potentially resulting in increased plant uptake and higher levels of foliar N, which could be expected to stimulate photosynthetic rates over the course of the entire growing season (hypothesized by Richardson et al. 2009a). Alternatively, earlier spring onset could potentially lead to greater display of leaf area, thereby enhancing light interception and canopy-level photosynthetic potential (Jolly et al. 2004; Luyssaert et al. 2007). Or, earlier spring could increase transpiration rates, leaving less moisture in the soil in summer and reducing summer (and annual) productivity (Leuzinger et al. 2005; Kljun et al. 2006). At the same time, during a cool, dry autumn with early senescence, decomposition of fresh litter could be inhibited, which would lead to a transient sink of labile carbon that might be quickly respired after snowmelt the following spring (hypothesized by Vesala et al. 2009). This last scenario would directly result in an increase in net C sequestration in autumn, but because of lagged effects this would be offset by a decrease in net sequestration the following spring.
Some of these ideas are summarized in the conceptual model shown in figure 1. Relative to a year where the timing of spring onset (which can be defined in a number of ways; see below) is ‘normal’, in figure 1a there is a direct effect of earlier spring: productivity increases are driven by the ‘extra’ days for assimilation and growth. This is the standard model described above. In figure 1b, lagged effects of earlier spring result in an increase in productivity during the peak of the growing season (e.g. owing to higher foliar N or leaf area index), but are not associated with a change in the timing of autumn senescence. In figure 1c, a shift towards earlier spring onset is accompanied by a shift towards earlier autumn senescence, and thus increased productivity in spring is offset by decreased productivity in autumn (e.g. owing to leaf ageing effects). In figure 1d, earlier spring onset is associated with decreased, rather than enhanced, productivity later in the growing season (e.g. owing to water limitation), and thus early-season gains are offset by sustained late-season reductions in physiological activity. Variations on these scenarios are of course possible, but given the uncertainty about which (if any) of these models is correct, it is clear that a more detailed exploration of relationships between phenology and productivity is justified (Norby et al. 2003; Piao et al. 2007).
Continuous (integrated over a 30-min time step, 24 h a day, 365 days a year) eddy covariance measurements of surface-atmosphere exchanges of CO2, water and energy are being made at hundreds of research sites around the world, many of which contribute data to the FLUXNET database (Baldocchi et al. 2001, 2008). With these data, the NEP (in grams of carbon per square-metre per unit time) at time scales from hours to years can be quantified, and relationships between productivity and environmental drivers characterized. The measured net exchange can be decomposed to the constituent fluxes GEP and ER, where, by definition, NEP = GEP − ER. This partitioning yields information about the ecosystem-scale processes controlling spatio-temporal variation in the uptake and release of CO2.
Here, we conduct an analysis of relationships between forest productivity and the timing of spring and autumn phenological events, with respect to spatial patterns across sites, and interannual variation at individual sites. We use long-term records from 21 temperate and boreal sites within FLUXNET to quantify gross and net productivity as GEP and NEP, respectively. Traditional ground observations of ‘budburst’ and ‘senescence’ are not available for many of these sites, and thus as phenological indicators, we instead use (i) dates derived from the fluxes themselves, representing the phenology of ecosystem processes (e.g. Gu et al. 2009), and (ii) dates derived from remotely sensed canopy greenness, representing land surface phenology (Reed et al. 2009; White et al. 2009; Xiao et al. 2009), and dates derived from phenological models (Jolly et al. 2005; Stöckli et al. 2008).
By leveraging the high-frequency eddy covariance time series aggregated to a daily time step, we are able to obtain process-level insights at a fine temporal scale, which would simply not be possible using biometric inventories, conventional tree growth measurements or other integrated measures of productivity, e.g. as inferred from remote sensing or atmospheric CO2 concentrations (Myneni et al. 1997; Randerson et al. 1999). The empirical analysis also improves on model-based studies (e.g. Lieth 1975; White et al. 1999; Baldocchi & Wilson 2001; Jolly et al. 2004; Loustau et al. 2005; Euskirchen et al. 2006), which generally do not adequately describe indirect or lagged effects of phenological variation (Schwartz 1992; Richardson et al. 2009a). This is further exacerbated by the fact that plant phenology is itself poorly represented in most biogeochemical, ecosystem and land surface models (e.g. Levis & Bonan 2004; Kucharik et al. 2006; Ryu et al. 2008).
In the following analysis, we address four key questions:
— How does the timing of spring and autumn phenological transitions control ecosystem productivity, and do the effects of phenology on gross productivity, inferred from estimates of GEP, differ from those on net productivity? Based on previous work, we expect that respiration, while not as directly controlled by phenology as is photosynthesis, will be increased both by earlier spring onset and later autumn senescence, and this will to some degree reduce the impact on net productivity.
— How do relationships between phenology and ecosystem productivity across sites compare with those within sites (i.e. with respect to interannual variation)? Because variation in weather events throughout the year (and not just in spring and autumn) contributes to the interannual variability in productivity, we predict stronger relationships across sites and climatic gradients compared to within sites.
— Do deciduous broadleaf (DBF) and evergreen needleleaf forests (ENF) vary in sensitivity to indicators of phenological events, and how do these patterns compare in spring and autumn? Physiological differences between deciduous and coniferous species, reflecting divergent ecological strategies, provide a strong a priori basis to expect dissimilar sensitivities.
— Which (if any) of the four conceptual models presented in figure 1 is best supported by the data? Important to answering this question will be whether results give support for sustained indirect or lagged effects of spring phenology on summer and autumn productivity.
2. Data and methods
(a) Site selection and data processing
We used data from V2 of the FLUXNET ‘La Thuile’ dataset (http://www.fluxdata.org), in which the 30 min eddy covariance measurements (net ecosystem exchange, NEE, of CO2) have been standardized, gap-filled and partitioned to the component fluxes of ER and GEP using a set of common algorithms (Reichstein et al. 2005; Papale et al. 2006; Moffat et al. 2007). From the 253 sites with available data, we identified those temperate or boreal ecosystems with DBF or ENF vegetation that had not been recently disturbed or heavily managed. In these ‘summer-active’ ecosystems, temperature and day length, rather than water availability, are considered the main drivers for seasonal variation in phenology.
Working with the daily data files for each site, we further retained only those years for which 75 per cent or more of the 30 min periods consisted of actual measurements or were gap-filled with ‘high confidence’ (Category A in Reichstein et al. 2005), and the fraction of missing daily NEE values did not exceed 5 per cent. Sites for which a minimum of 5 years of data had passed these criteria were included in our analyses, yielding the 21 sites and 153 site-years given in table 1.
Negative NEE is defined as net CO2 uptake by the ecosystem. So that the meaning of increases in carbon storage is clear, we use the term NEP (NEP = −NEE), and consider the component fluxes ER and GEP to both be positive fluxes (thus NEP = GEP − ER).
(b) Phenological indicators
On the ground, observer-based, phenological observations are unavailable for most FLUXNET sites, and thus we relied on phenological indicators (i) derived from the flux time series themselves, and (ii) derived from remote sensing observations or modelling analyses.
(i) Indicators derived from flux measurements
A number of previous studies (including, for example, Gu et al. 2003, 2009; Churkina et al. 2005; Barr et al. 2007; Richardson et al. 2009a) have used seasonal inflection points in the time series of NEP and GEP to characterize transition dates associated with the phenology of ecosystem processes. We used smoothing splines to determine the first and last dates at which the NEP CO2 source–sink transitions occurred, and the corresponding dates at which NEP reached a threshold value of 1 g C m−2 d−1 net uptake (figure 2). Two types of indicators were derived from the GEP time series: (i) first and last dates at which specific absolute thresholds of daily GEP (=2, 4 and 6 g C m−2 d−1) were reached, and (ii) first and last dates at which relative thresholds (defined in terms of GEPmax, the peak value of daily GEP) of daily GEP were reached—i.e. 25%, 50% and 75% of GEPmax.
Our objective with these different thresholds is to identify dates with which the phase of the seasonal cycles of GEP and NEP may be characterized. We apply both relative and absolute GEP thresholds specifically because GEPmax varies among sites, and thus the functional relevance of any particular absolute threshold may likewise vary among sites. Similarly, we include a range of different NEP and GEP thresholds because we do not have strong a priori beliefs about which indicator will prove to be ‘best’ or most useful.
(ii) Indicators derived from remote sensing and models
Additional phenological indicators were obtained from several remote sensing platforms. For deciduous forests, the dominant signal seen by satellite-based sensors is the seasonal pattern of canopy development and abscission; for evergreen forests, changes in ‘greenness’ can be attributed to seasonal variation in canopy biochemistry, the production of new foliage by canopy species, and, particularly where the overstory is sparse, the phenology of understory vegetation.
The MOD12Q2 (Collection 4) global land cover dynamics product (Zhang et al. 2006), based on data from NASA's MODIS instrument aboard the Terra and Aqua spacecraft, is produced at 1 km spatial resolution using 16 day nadir bi-directional reflectance distribution function (BRDF)-adjusted reflectance. Two transition dates, corresponding to the onset of spring ‘green-up’ and ‘maturity’, and autumn ‘senescence’ and ‘dormancy’, were estimated in each of spring and autumn from inflection points in sigmoid functions fit to the MODIS EVI (enhanced vegetation index) time series for the pixel containing each tower.
The European Commission Joint Research Center (JRC) fraction of absorbed photosynthetically active radiation (fAPAR) product (Gobron et al. 2001, 2006; data from http://fapar.jrc.ec.europa.eu/) derived from SeaWiFS sensor data, is a 10-day time composite product spatially averaged over a 3 × 3 pixel matrix (approx. 6 × 6 km), centred on each eddy flux tower. The first and last dates of half-maximum JRC fAPAR were determined using a logistic function for interpolation (e.g. as in Richardson et al. 2009b).
Following procedures designed for boreal forests (Delbart et al. 2005, 2006), 10-day satellite pour l'observation de la terre (SPOT) and advanced very high resolution radiometer (AVHRR) data were used in conjunction with the normalized difference water index (NDWI) and a pixel-specific threshold (PST) approach to estimate the onset of spring green-up at both 0.1° (NDWI-PST composite) and 1 km (SPOT NDWI only) spatial resolutions around each tower.
As a model-based alternative to remote-sensing approaches, we calculated the growing season index (GSI) of Jolly et al. (2005). GSI uses photoperiod, evaporative demand and minimum temperature to predict seasonal variation in vegetation greenness. GSI values fall on a scale between 0 and 1; we identified the first and last dates at which a half-maximum threshold value of GSI was predicted to occur.
A final indicator applied data-model fusion techniques to merge remote-sensing and model-based approaches. As described more fully by Stockli et al. (2008), MODIS data products were used to constrain parametrization of the PROGNOSTIC phenology model, which is based on Jolly et al.'s (2005) GSI. We identified the first and last dates of half-maximum fAPAR predicted by the optimized model.
Missing values (e.g. MODIS data are not available before 2001; JRC fAPAR data are not available before 1997) were filled using a linear model conditioned on the non-missing observations for each site.
(c) Statistical analyses
We calculated pairwise Pearson correlations, separately for ENF and DBF, between phenological indicators derived from CO2 flux measurements and indicators derived from remote sensing and models. This analysis was conducted using both mean dates for each site, and year-to-year anomalies at a given site, to examine agreement among metrics with respect to spatial and temporal patterns, respectively.
We used analysis of covariance to quantify the degree to which the different ‘start of season’ and ‘end of season’ phenological indicators explained spatial (across site means) variation in early- (January–June) and late-season (July–December) GEP and NEP integrals, and to test whether the slopes of these relationships differed between ENF (n = 12 sites) and DBF (n = 9 sites) forests. Our model included forest ‘vegetation type’ as a class variable, ‘phenological indicator’ (day of year) as a covariate and a ‘vegetation type’ × ‘phenological indicator’ interaction effect. Thus the y-axis intercept, and the slope of the relationship between CO2 flux integral and phenological indicator, could vary between DBF and ENF.
Similarly, analysis of covariance was used to evaluate temporal patterns, i.e. the response of CO2 fluxes to inter-annual variation in phenology. This analysis was based on year-to-year anomalies in both fluxes and phenological indicators for each site. Again, we included a vegetation type × phenological indicator interaction effect in the predictors to test for differences between ENF (n = 81 site-years) and DBF (n = 62 site-years) forests. An additional interaction effect, with sites nested within vegetation type, allowed the slope of relationships between phenological anomalies and flux anomalies to vary among sites. In the analysis that follows, we focus on the overall pattern across all sites within each vegetation type.
Finally, to quantify lagged effects, we used correlation analysis to investigate relationships between anomalies in flux integrals over successive six month periods, and whether these anomalies were also correlated with antecedent phenological anomalies.
(a) Spatial relationships among different phenological indicators
Spatial patterns indicated by phenological indicators derived from remote sensing and modelling were, in some but not all instances, in good agreement with the spatial patterns of variation in phenological indicators derived from flux measurements. We provide some examples here; the full correlation matrix is shown in the electronic supplementary material, table S1.
Across sites, the mean date at which the JRC fAPAR product reached its half-maximum in spring was well-correlated with the mean NEP springtime source–sink transition date across both DBF (r = 0.81, p < 0.01, n = 9) and ENF (r = 0.92, p < 0.01, n = 12) forests (figure 3a). However, in spite of this strong correlation, ENF half-maximum fAPAR was, on average, not reached until 50 days after the NEP source–sink transition. Furthermore, these correlations were weak for the corresponding autumn mean half-maximum and sink–source transition dates (r = 0.68, p = 0.04 and r = 0.29, p = 0.35, respectively; figure 3b).
The mean MODIS greenup date in spring was well-correlated with mean NEP source–sink transition date across DBF (r = 0.88, p < 0.01) and ENF (r = 0.84, p < 0.01), whereas mean MODIS maturity date was better correlated with indicators related to GEP, e.g. the mean date at which half-maximum daily GEP was achieved (r = 0.87, p < 0.01, and r = 0.95, p < 0.01, respectively).
For DBF, the mean spring and autumn dates of GSI half-maximum were correlated with mean first (r = 0.86, p < 0.01) and last (r = 0.74, p = 0.02) dates, respectively, at which daily NEP = 1 g C m−2. By comparison, for ENF, mean GSI half-maximum correlated better with mean first (r = 0.93, p < 0.01) and last (r = 0.86, p < 0.01) dates at which daily GEP = 6 g C m−2.
(b) Interannual variation across different phenological indicators
With regard to interannual variation at a given site, there was generally a strong overall coherence among different phenological indicators; examples for two sites are shown in figure 4. Consequently, temporal anomalies in flux-derived phenological indicators tended to be reasonably well correlated with temporal anomalies in remote sensing and model-derived indicators. Again, we provide some examples here, and present the full correlation matrix in electronic supplementary material, table S2.
For DBF, anomalies in all remote sensing-derived spring indicators were significantly correlated with anomalies in the first date at which daily GEP = 2 g C m−2 (mean r = 0.52 ± 0.08, all p < 0.01) and with anomalies in the NEP source–sink transition date (mean r = 0.47±0.07, all p < 0.01). The corresponding autumn indicators were, in most cases, not significantly correlated (figure 3c,d and electronic supplementary material, table S2).
For ENF, although there were some statistically significant correlations between anomalies in remotely sensed phenological indicators and anomalies in flux-derived phenological indicators, MODIS greenup was the only spring indicator for which any such r > 0.5; the best correlation was with the date at which daily GEP = 4 g C m−2 d−1 (r = 0.58, p < 0.01). In autumn, anomalies in MODIS senescence and dormancy were not significantly correlated with anomalies in any of the flux-derived phenological indicators (all p > 0.10). With the exception of JRC fAPAR and GSI half-maximum dates, for which several correlations were significant at p < 0.01, results for other remotely sensed indicators were similar (figure 3d and electronic supplementary material, table S2).
Results here and in the previous section demonstrate that the phenological indicators derived from remote sensing were better able to capture spatial (across-site), compared with temporal (interannual), patterns of variability in phenological indicators derived from CO2 flux time series. With respect to interannual variability, there was better agreement among independent phenological indicators in spring than in autumn, and better agreement for DBF when compared with ENF.
(c) Relationships between phenological indicators and CO2 fluxes across sites
We now focus on quantifying the degree to which different ‘start of season’ and ‘end of season’ phenological indicators explained across-site variation in January–June and July–December integrals of GEP and NEP, and whether this varied between ENF and DBF. Our analysis was conducted using these six month periods so that ‘seasonal’ CO2 flux integrals are calculated over the same time period for each site. Given the wide range (across sites) of spring and autumn transition dates, these half-year periods are essentially the lowest common denominator (cf. quarterly integrals in Richardson et al. 2009a).
The mean JRC fAPAR half-maximum date and mean MODIS greenup date consistently explained more of the spatial variation in mean January–June flux integrals than other remote-sensing or model-derived indicators (selected results shown in table 2; full results in electronic supplementary material, table S3). With respect to flux-derived indicators, absolute GEP threshold dates (e.g. daily GEP = 4 g C m−2) were the best predictors of mean January–June GEP, whereas mean NEP source–sink transition dates were the best predictors for January–June NEP. Overall, some relatively consistent patterns emerge when we focus on these indicators (table 2 and figure 5). First, the slopes of the relationships between mean phenological indicators and mean flux integrals generally did not differ significantly between ENF and DBF sites. Second, the slopes of the relationships were more or less consistent across the different phenological indicators, with coefficients of variation of 30 per cent or less. From this analysis, we estimate that, across sites, a 1 day advancement of spring increased total GEP over the January–June period by roughly 12 g C m−2. This was offset by increased ER losses of 8 g C m−2, resulting in a net increase in carbon uptake of 4 g C m−2 for every 1 day advancement of spring onset.
Compared with the spatial patterns in spring, mean autumn phenological indicators poorly explained the spatial patterns of July–December fluxes, and there were inconsistencies among indicators in the retrieved slope coefficients (table 2 and electronic supplementary material, table S3). Statistically significant differences were rarely observed between DBF and ENF. As in the spring, autumn mean GEP threshold dates best explained the spatial variability in July–December GEP integrals, whereas the autumn mean NEP sink–source transition date best explained the spatial variability in NEP integrals. The remote sensing-derived phenological indicators were typically unable to explain spatial patterns of autumn GEP or NEP. Across the selection of indicators in table 2, and both forest types, later autumn was spatially associated with about a 6 g C m−2 d−1 increase in July–December integrals for both GEP and NEP, although slope coefficients were, as noted, highly variable among phenological indicators.
(d) Relationships between interannual phenological anomalies and CO2 flux anomalies
The overall patterns between interannual phenological anomalies and corresponding GEP and NEP flux anomalies were relatively consistent, although the retrieved slope coefficients differed in magnitude among phenological indicators (electronic supplementary material, table S4). Phenological indicators derived from remote sensing or models were unable to explain more than 40 per cent of the interannual variability in January–June GEP or NEP, whereas some indicators derived from the flux time series explained 50–80% of this variability (table 3 and figure 6). As was the case across sites (table 2), it was more difficult to phenologically explain interannual variability in CO2 fluxes from July–December compared with January–June (table 3). Estimated slope coefficients in autumn were, in many instances, not significantly different from zero for most indicators derived from remote sensing (table 3 and electronic supplementary material, table S4).
When the difference between DBF and ENF sites was statistically significant, GEP and NEP of DBF sites were consistently more sensitive to phenological anomalies. Across the indicators in table 3, a 1 day phenological anomaly (towards earlier spring onset) for DBF was associated with a 7 g C m−2 increase in January–June GEP, and a 4 g C m−2 increase in January–June NEP. Corresponding values for ENF sites were increases of 5 and 3 g C m−2 for January–June GEP and NEP, respectively. In autumn, a 1 day phenological anomaly (towards later senescence), was associated with an 8 g C m−2 increase in DBF GEP, compared with just 3 g C m−2 in ENF GEP. Concurrent increases in ER losses resulted in an overall increase in total July–December NEP by 5 g C m−2 d−1 for DBF, but no increase in NEP of ENF, in response to a 1 day autumn phenological anomaly.
We conducted a similar analysis on annual GEP and NEP integrals, including both spring onset and autumn senescence indicators as covariates (electronic supplementary material, table S5). A much higher proportion of interannual variability in GEP, compared with NEP, could be explained by the phenological indicators; similarly, indicators derived from fluxes explained more of the variability than those derived from remote sensing or models. Estimated GEP sensitivities to spring phenological indicators were markedly larger for the annual sums compared with the January–June integrals, offering indirect evidence of lagged effects of spring phenological anomalies.
(e) Productivity and phenological anomalies in relation to environmental drivers
Anomalies in spring (defined as the ±30 days surrounding mean onset date, for each site, across all phenological indicators) air temperature, but not precipitation or solar radiation, explained a considerable proportion of the interannual variability in both spring and January–June productivity. The observed sensitivity of GEP to air temperature (i.e. increase in spring GEP resulting from a +1°C air temperature anomaly) was greater for DBF (35 ± 5 g C m−2 °C−1) than ENF (20 ± 3 g C m−2 °C−1) forests (figure 7a,b). However, these gains were offset by concurrent increases in ER (15 ± 3 and 11±2 g C m−2 °C−1, respectively). Thus, a +1°C spring air temperature anomaly was associated with a 20 ± 3 g C m−2 increase in DBF NEP, and a 9 ± 2 g C m−2 increase in ENF NEP. In autumn, no statistically significant relationships (all p > 0.10) between autumn air temperature anomalies and autumn or July–December GEP, ER or NEP anomalies were observed (figure 7c,d).
Anomalies in spring air temperature were also associated with anomalies in the timing of spring phenological events. For most phenological indicators, these correlations were significant (p < 0.05). Generally, a +1°C spring air temperature anomaly led to an approximately 3 day advancement in spring onset date for both DBF and ENF. By comparison, in autumn, correlations between temperature anomalies and phenological anomalies were much weaker, and for many phenological indicators, not statistically significant; overall, a +1°C autumn air temperature anomaly led to only an approximately 1 day delay in autumn senescence date.
(f) Lagged effects of anomalies in GEP and ER
Correlation analyses indicated potential lagged effects of GEP and ER anomalies on ecosystem function. A 1 g C m−2 increase in July–December GEP was associated with a GEP increase of roughly 0.8 g C m−2 during the following January–June period (r = 0.34, p < 0.01; figure 8a). Autumn ER anomalies were also positively correlated with ER anomalies during the following January–June (r = 0.22, p < 0.01). However, there was only weak evidence for relationships between July and December GEP or ER anomalies and subsequent NEP anomalies.
Similarly, a 1 g C m−2 increase in January–June GEP was associated with a GEP increase of 1.2 g C m−2 during the following July–December period (r = 0.41, p < 0.01; figure 8b); a comparable increase in ER was associated with an ER increase of 1.5 g C m−2 during the following July–December (r = 0.36, p < 0.01). There was no consistent evidence for lagged effects of spring GEP or ER anomalies on subsequent NEP.
(g) Lagged effects of phenological anomalies
We did not see any consistent evidence suggesting that interannual variation in the timing of spring phenological transitions had any effect on subsequent dates of autumn phenological transitions (i.e. lagged effects as suggested in figure 1c). There were no statistically significant correlations observed between anomalies in the first and last dates at which (i) the NEP source–sink transition occurred, (ii) daily GEP reached a threshold value of 2 g C m−2, or (iii) either JRC fAPAR or GSI reached half maximum. Only for half-maximum PROGNOSTIC fAPAR was there a significant correlation between anomalies in the first and last dates (r = −0.27, p = 0.03 for DBF; r = −0.36, p < 0.01 for ENF).
At the same time, we did not generally observe statistically significant correlations between spring phenological anomalies and rates of physiological activity during summer months, or summer or autumn GEP anomalies. One notable exception was that for both forest types, anomalies in the spring date at which daily GEP reached 50 per cent of GEPmax were negatively correlated (all r = −0.3 to −0.4, p < 0.01) with anomalies in July photosynthetic light use efficiency and maximum photosynthetic rate (at the ecosystem scale), as well as anomalies in integrated July and July–December GEP. Thus, for this phenological indicator, earlier spring onset was associated with increases in gross productivity later in the summer (and hence the observed negative correlation coefficients). In parallel to this and building on results shown in figure 8b, anomalies in January–June GEP were positively correlated (r ≈ 0.5, p < 0.01), for both forest types, with anomalies in July GEP, as well as July photosynthetic light use efficiency and maximum photosynthetic rate. Thus, increased GEP in spring was also associated with increases in productivity later in the summer and autumn (and hence the observed positive correlation coefficient).
To frame the discussion, we return to the four questions posed in §1.
(a) Gross versus net productivity
At the ecosystem scale, differences between gross and net productivity are mediated by respiration, which, like photosynthesis, is temperature-sensitive. Recent photosynthetic products also provide substrate for the autotrophic component of ER, and in the autumn, heterotrophic respiration is stimulated by leaf litter inputs, representing a pulse of highly labile C. Thus, there is a strong but indirect connection between phenology and respiration.
In both spring and autumn, and with respect to both spatial patterns across sites and interannual variation at individual sites, increased GEP resulting from an ‘extra’ day during the active season tended to be offset by a concurrent, but smaller, increase in ER; thus, in most of the cases examined, NEP also increased in response to either earlier spring onset or delayed autumn senescence (tables 2 and 3). Temperature anomalies appeared to be the main driver of phenological anomalies, particularly in spring. About 50 per cent of the increase in GEP in years with a warmer than average spring tended to be offset by respiratory losses, resulting in a reduced, but still substantial, increase in NEP (figure 7a,b). However, we were unable to document statistically significant impacts of autumn temperature anomalies on net or gross productivity (figure 7c,d), although warmer temperatures appeared to be correlated with delayed senescence for many autumn phenological indicators.
Because the ecosystems studied are all temperature-limited during the dormant season, these patterns are consistent with expectations: increases in metabolic activity, whether related to photosynthetic or respiratory processes, should result from temperature-driven advances in spring onset or delays in autumn senescence. Patterns similar to these have been reported in previous analyses using subsamples of the present dataset (Barr et al. 2009; Richardson et al. 2009a). However, unlike Dunn et al. (2007) or Piao et al. (2008), we did not generally observe instances where the increased respiratory losses were sufficient to offset concurrent GEP gains.
Vesala et al. (2009) suggested that observed increases in autumn respiration owing to warming could be the result of (i) more rapid turnover of ‘new’ leaf litter inputs (i.e. a transient C source which would be offset by a compensating decrease in respiration the following spring), or (ii) an overall increase in respiratory losses from the ecosystem, perhaps resulting from enhanced decomposition of ‘older’ soil C. We observed positive, rather than negative, correlations between respiratory anomalies in spring and the following autumn, and between anomalies in autumn and the following spring. These patterns do not support scenario (i) of Vesala et al. (2009). This, plus the fact that the geometric mean slope of the relationship between autumn GEP anomalies (x) and autumn ER anomalies (y) is close to unity (compared with >1.5 in spring), suggest an alternative interpretation, whereby increased autotrophic respiration is driven by recent photosynthesis, and is entirely accounted for (i.e. does not represent a change in the overall C balance of the system) by concurrent increases in GEP during warmer autumns. However, discriminating between this and scenario (ii) of Vesala et al. (2009) would require additional information about the quality and age of the substrate being used to support the increased ER. In general, the measurements needed for this kind of analysis (e.g. isotopic signatures of respiratory components; Trumbore 2006; Carbone et al. 2008) are not being made at most FLUXNET sites.
(b) Spatial versus temporal patterns
For almost all phenological indicators considered, the slope of relationships between early-season (January–June) fluxes and phenology was about twice as steep across sites when compared with within sites; this pattern held for GEP, ER, and for ENF sites, NEP (tables 2a and 3a). In autumn, ENF similarly exhibited greater sensitivity to phenological variation across sites compared with over time. Churkina et al. (2005) similarly predicted (but did not test) that spatial patterns of productivity in relation to growing season length across sites should be stronger than interannual responses to variation in growing season length at a single site. Our data give some support for this hypothesis (cf. Barr et al. 2009). However, complicating the interpretation of spatial patterns in NEP is the strong influence of forest age on net carbon balance (Piao et al. 2009). Even if phenological transition dates in spring and autumn were similar for two adjacent stands, younger forests are expected to have high NEP and older forests, low NEP.
Our results highlight the challenge of trying to equate spatial and temporal variability; as suggested by Burke et al. (1997), the ‘structure’ (a combination of vegetation structure and biogeochemical structure, as for example species composition and nutrient pools, respectively) of systems constrains temporal variability in ways that are different from spatial patterns of ecosystem function. Thus, spatial patterns mostly represent equilibrium responses and the interaction of covarying factors (climate, nutrient availability, species composition, etc.). By comparison, temporal patterns mostly reflect transient responses to external forcing (weather, disturbances, etc.), which are dependent on ecosystem structure. Comparing two previous studies, Burke et al. (1997) concluded that the spatial patterns of aboveground net primary productivity in relation to precipitation (Sala et al. 1988) were much stronger than the interannual patterns of productivity in relation to rainfall at individual sites (Lauenroth & Sala 1992). In the context of the present study, interannual patterns of carbon uptake and release are likely to be driven by internal system dynamics, particularly related to carbon pools with intermediate (months–years) turnover times and high sensitivity to climatic drivers (e.g. decomposition of fresh, and hence labile, leaf litter; Richardson et al. 2007). Thus, transient sinks or sources of CO2 are a major determinant of patterns of interannual variation (Vesala et al. 2009), whereas these just contribute noise to patterns of spatial variation. In all likelihood, neither of these is directly analogous to ecosystem responses to a gradually changing climate, especially as structural shifts in response to climate are likely to be slow. A natural conclusion to draw from this is that caution should be exercised when attempts are made to estimate responses to future climate change from either short-term or space-for-time studies.
(c) Differences between ENF and DBF
Across sites, slopes of spatial relationships between mean flux integrals and mean dates of phenological transitions tended not to be significantly different between DBF and ENF forests. In the very few autumn cases where the difference was significant, ENF forests were more sensitive than DBF forests (table 2). With respect to interannual variation in fluxes and phenology within sites, patterns were reversed. Statistically significant differences in slopes of temporal relationships (DBF > ENF) were commonly observed for flux-derived phenological indicators for both GEP and NEP (table 3). The productivity of DBF was also more sensitive to interannual spring air temperature anomalies when compared with ENF (figure 7). Previous empirical (e.g. Welp et al. 2007; Barr et al. 2009; Richardson et al. 2009a), and other studies cited in §1 (note that in many cases some of the same data were used in the present study, and thus these results are not fully independent of those presented here) and model-based (Loustau et al. 2005; Piao et al. 2007) analyses have reported similar differences. For example, based on patterns across sites, Churkina et al. (2005) found that the sensitivity of measured NEP to growing season length (‘carbon uptake period’) was 5.8 ± 0.7 g C m−2 d−1 in deciduous forests, compared with 3.4 ± 0.3 g C m−2 d−1 in ENF. Similarly, Piao et al. (2007) reported that in DBF, the sensitivity of modelled GEP to growing season length (in days) was 9.8 ± 2.6 g C m−2 d−1, compared with just 4.9 ± 2.5 g C m−2 d−1 in ENF. These patterns reflect not only the fact that DBF tend to have higher productivity per day during the growing season compared with ENF (e.g. Schulze et al. 1977; Falge et al. 2002; Piao et al. 2007; Richardson et al. 2009a), but also fundamental differences in phenological strategy between the two forest types (Barr et al. 2009). The production of new foliage is a prerequisite for photosynthesis in deciduous stands, and while early leaf-out or delayed senescence could be advantageous in terms of carbon uptake, there are tradeoffs between increasing growing season length versus increases in the probability of early-spring or late-autumn frost damage (Saxe et al. 2001; Norby et al. 2003; Gu et al. 2008). By comparison, the recovery of photosynthesis in ENF stands is reversible, proceeds slowly and involves multiple steps (Monson et al. 2005). While the evergreen strategy is more conservative, it appears that a consequence of this is reduced potential for productivity gains in response to earlier or warmer spring onset.
The observed differences between ENF and DBF in phenology–productivity relationships have potential ecological consequences, especially in light of recent (and anticipated future) effects of climate change on phenology (Kramer 1994). The climatic determinants of phenology are known to differ among species (e.g. Richardson & O'Keefe 2009). Models predict that because of this, species will respond in different ways to future climate change (Morin et al. 2008). Our results suggest that gross productivity of DBF is somewhat more sensitive to phenological variability (at least over time, if not across space) than that of ENF. This has implications for the effects of future phenological shifts on competition for light and other resources, including space, in mixed species stands (Kramer et al. 2000; Rötzer et al. 2004; Welp et al. 2007). This is in addition to ways in which climate change will directly affect productivity, and presumably fitness, through temperature and moisture effects on growth conditions and habitat suitability (e.g. Mohan et al. 2009).
(d) Selecting among competing conceptual models
Based on the fact that the influence of earlier springs tended to be about twice as large when annual, rather than springtime, flux integrals were considered, Richardson et al. (2009a) concluded that there must be substantial lagged or indirect effects associated with the timing of spring onset. Vaganov et al. (2009) also found evidence in the autocorrelation structure of measured tree-ring widths for ‘carry-over’ effects of growing conditions in one year that tended to influence tree growth rates in subsequent years as well. Results presented here clearly indicate large direct effects of the timing of spring on early-season flux integrals (tables 2a and 3a), and also seem to indicate potential lagged effects as well. Hypotheses of lagged effects (figure 1b) are supported by observed correlations between January–June GEP anomalies and subsequent July–December GEP and ER anomalies. However, given that only for one phenological indicator (the date at which relative GEP reached 50% of GEPmax) were spring onset anomalies significantly correlated with July or July–December flux anomalies, or measures of July photosynthetic capacity, the evidence for phenological control of these lagged anomalies is not very strong.
We developed the four scenarios shown in figure 1 as a conceptual model with which to explore the potential effects of variability in spring onset. This analysis does not support scenario (c), as there was no correlation between anomalies in spring and autumn transition dates. Similarly, the analysis does not support scenario (d), as there was no general evidence of reduced productivity in late summer being correlated with earlier spring onset, or increased spring productivity. As discussed in the preceding paragraph, we see some evidence for selecting scenario (b) over scenario (a), although it is not clear that increased productivity in summer and autumn can be attributed to ‘early spring’ phenological anomalies. Rather, phenological anomalies are just one source of variability in spring productivity (explaining, for example, no more than 62% of the variability in spring NEP; see table 3), and it is spring productivity anomalies, not spring phenological anomalies per se, that are best correlated with lagged productivity anomalies. As discussed by Gu et al. (2003), other factors, such as the rate of green-up and development (or recovery) of photosynthetic capacity, are also important regulators of spring productivity anomalies, but are not captured by the indicators we use here.
5. Summary and conclusions
Using continuous measurements of surface-atmosphere CO2 exchange from 21 FLUXNET sites, in conjunction with a range of phenological indicators for spring and autumn transition dates, we have evaluated the degree to which phenological transitions control productivity of DBF and ENF in the temperate–boreal zone. We have quantified relationships between phenology and productivity, both across space in relation to phenological gradients, and with respect to interannual variation in phenology. However, we caution against attempts to predict future climate change responses from either spatial or temporal patterns (see also Piao et al. 2009). That being said, the results presented here represent a novel opportunity to evaluate predictions of state-of-the art land surface models, and to test the way in which lagged and indirect effects of phenology are thus represented (e.g. Kucharik et al. 2006). There is also the potential to test, using experimental studies, some of the ecological hypotheses raised here—in particular, investigation of the mechanisms underlying the observed lagged effects is important so that these can be more accurately modelled.
Our analysis essentially presumed that specific dates representing, for example, ‘spring onset’ or ‘autumn senescence’ could be used to represent the timing of the start or end of the entire seasonal trajectory of development from dormancy through activity, senescence and back to dormancy. Even ignoring the fact that the velocity of developmental changes is perhaps as important as the timing of those changes (Gu et al. 2003), it is clear that unequivocal determination of dates representing the start and end of the growing season—or even the length of the growing season—is difficult (Gu et al. 2009), as different definitions (with different physiological connotations) abound and results may be sensitive to the specific definition used (White & Nemani 2003; Piao et al. 2007; also tables 2 and 3). In the present study, although the slopes of relationships between phenological indicators and flux integrals varied depending on the phenological indicator in question, we observed patterns that were generally consistent among indicators, and also with the results of previous empirical studies (e.g. Baldocchi 2008; Barr et al. 2009; Richardson et al. 2009a).
Future studies could improve on the present analysis in a number of ways. While we have tried to characterize relationships between remotely sensed phenology and physiological changes on the ground, additional data—such as could be provided by ground observations, tower-based radiometric data or webcam imagery (Richardson et al. 2009a,b)—would be extremely valuable in this regard. These would also contribute to improved understanding of spatial and temporal patterns of seasonal variation in canopy structure, and how this relates to the phenology of ecosystem processes. Further improvements to remote sensing algorithms for phenology are also needed; important targets include distinguishing phenologies of overstory and understory vegetation, and resolving issues related to spatial heterogeneity and disaggregation of phenological signals in pixels with mixed landcover.
We have focused here on patterns across a range of DBF and ENF ecosystems; we have not attempted to investigate how (or why) patterns differ among sites with similar vegetation—for example, why some sites are more sensitive, and others less sensitive, to phenological variability. Stand age, disturbance history and species composition are factors likely to be important.
Given the rich data contained in the FLUXNET database, our study could easily be expanded to a wider range of sites, including ecosystems where both phenology and productivity are more tightly regulated by moisture supply. Such an analysis would logically link phenology to latent heat fluxes (evapotranspiration), and emphasize the coupled nature of the carbon and water cycles. A similar analysis could also focus on phenological control of other vegetation feedbacks to the climate system, including, for example, surface energy balance and biogenic emission of volatile organic compounds.
A.D.R. acknowledges support from both the Northeastern States Research Cooperative and the Office of Science (BER), US Department of Energy, through the Terrestrial Carbon Programme under Interagency Agreement No. DE-AI02-07ER6455, and through the Northeastern Regional Center of the National Institute for Climatic Change Research. A.D.R. also thanks the participants in the 2009 Kittery Point workshop for their feedback and suggestions.
We thank the FLUXNET Site PIs for contributing data, the agencies and institutions that funded long-term measurements at these sites, and the networks to which the sites belong, particularly AmeriFlux, Fluxnet-Canada and CarboEurope-IP. This project resulted from the 2007 La Thuile FLUXNET workshop, which would not have been possible without the financial support provided by CarboEuropeIP, FAO-GTOS-TCO, iLEAPS, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, US Department of Energy. Moreover, we acknowledge technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, University of California-Berkeley and University of Virginia.
One contribution of 11 to a Theme Issue ‘The role of phenology in ecology and evolution’.
- advanced very high resolution radiometer
- AVHRR PST
- spring transition date estimated from pixel-specific threshold from AVHRR data
- bi-directional reflectance distribution function
- deciduous broadleaf forest
- evergreen needleleaf forest
- ecosystem respiration
- enhanced vegetation index
- fraction of absorbed photosynthetically active radiation
- gross ecosystem photosynthesis (integrated to daily, half-yearly and annual sums)
- peak value of GEP, estimated via smoothing spline
- GEP = 2, 4, 6 g C m−2 d−1
- first/last dates of spring/autumn absolute threshold (2, 4, 6 g C m−2 d−1) values for GEP
- growing season index
- GSI half-maximum
- date of spring/autumn half-maximum growing season index
- European Commission Joint Research Council
- JRC fAPAR half-maximum
- date of spring/autumn half-maximum fAPAR from SeaWiFS sensor
- moderate resolution imaging spectroradiometer
- MODIS dormancy
- second autumn transition date estimated from MODIS EVI
- MODIS greenup
- first spring transition date estimated from MODIS EVI
- MODIS maturity
- second spring transition date estimated from MODIS EVI
- MODIS senescence
- first autumn transition date estimated from MODIS EVI
- normalized difference water index
- composite metric, combining AVHRR PST (1990–1997, 1999) and SPOT NDWI (1998, 2000–present), 0.1° resolution
- net ecosystem exchange
- net ecosystem productivity (= GEP−ER=−1×NEE)
- NEP = 1 g C m−2 d−1
- date of spring/autumn NEP = 1 g C m−2 d−1 threshold
- NEP source/sink
- date of spring/autumn source/sink transition in net ecosystem productivity
- PROGNOSTIC half-maximum
- date of spring/autumn half-maximum fAPAR from PROGNOSTIC model tuned to MODIS data
- pixel-specific threshold
- daily GEP relative to GEPmax (e.g. rGEP=25%)
- rGEP = 25, 50, 75% GEPmax
- first/last dates of spring/autumn relative threshold (25, 50, 75% of GEPmax) values for daily GEP
- satellite pour l'observation de la terre
- SPOT NDWI
- spring transition data estimated from normalized difference water index from SPOT data
- © 2010 The Royal Society