Edinburgh Explorer Biosphere-atmosphere exchange of CO2 in relation to climate: a cross-biome analysis across multiple time scales

. The net ecosystem exchange of CO 2 (NEE) varies at time scales from seconds to years and longer via the response of its components, gross ecosystem productivity (GEP) and ecosystem respiration (RE), to physical and biological drivers. Quantifying the relationship between ﬂux and climate at multiple time scales is necessary for a comprehensive understanding of the role of climate in the terrestrial carbon cycle. Orthonormal wavelet transformation (OWT) can quantify the strength of the interactions between gappy eddy covariance ﬂux and micrometeorological measurements at multiple frequencies while expressing time series variance in few energetic wavelet coefﬁcients, offering a low-dimensional view of the response of terrestrial carbon ﬂux to climatic variability. The variability of NEE, GEP and RE, and their co-variability with dominant climatic drivers, are explored with nearly one thousand site-years of data from the FLUXNET global dataset consisting of 253 eddy co-Correspondence variance research sites. The NEE and GEP wavelet spectra were similar among plant functional types (PFT) at weekly and shorter time scales, but signiﬁcant divergence appeared among PFT at the biweekly and longer time scales, at which NEE and GEP were relatively less variable than climate. The RE spectra rarely differed among PFT across time scales as expected. On average, RE spectra had greater low frequency (monthly to interannual) variability than NEE, GEP and climate. CANOAK ecosystem model simulations demonstrate that “multi-annual” spectral peaks in ﬂux may emerge at low (4+ years) time scales. Biological responses


Introduction
Variability in the global carbon cycle is dominated by terrestrial ecosystem metabolism (Houghton, 2000;Canadell et al., 2007) and it is critical to understand how and why the terrestrial carbon cycle varies to advance in our knowledge of the Earth system. The net ecosystem exchange of CO 2 (NEE) between the biosphere and atmosphere is determined by the response of its components, gross ecosystem productivity (GEP) and ecosystem respiration (RE) to physical and biological drivers which vary at multiple time scales. Modeling and time series analysis of eddy covariance (EC) measurements of NEE have demonstrated that its variability is largely determined by physical (i.e. climatic) controls at short controls at short (e.g. daily) time scales and by biological responses to climatic varaibility at longer (e.g. seasonal and interannual) time scales Richardson et al., 2008). Process-based ecosystem models often experience difficulties in replicating interannual variability in NEE (Siqueira et al., 2006;Urbanski et al., 2007), suggesting that our knowledge of the carbon cycle at multiple time scales must improve.
Models synthesize our knowledge of terrestrial carbon exchange and thus represent an explicit hypothesis about how ecosystem function transfers variability in climatic forcing to an ecological response, namely the flux of mass or energy. These ecosystem transfer properties can now be explored using long-term meteorological and carbon flux measurements from the international FLUXNET project, which consists of eddy covariance tower flux measurements from regional networks [CarboeuropeIP, AmeriFlux, Fluxnet-Canada, LBA, Asiaflux, Chinaflux, USCCC, Ozflux, Carboafrica, Koflux, NECC, TCOS-Siberia, (e.g., Aubinet et al., 2000;Baldocchi et al., 2001a;Baldocchi, 2008)]. The FLUXNET project offers the unprecedented opportunity to relate directly the measured variability in mass and energy flux to the measured meteorological variability at time scales from hours to years, and in some cases over a decade Urbanski et al., 2007;Granier et al., 2008).
Information on ecosystem properties (e.g. canopy N, leaf area index, vegetation growth) in the FLUXNET database is increasingly available, but currently lacking are measurements of how ecosystem properties change over time. Such measurements are needed to synthesize how climate and ecology interact to dictate NEE across time scales in global ecosystems. In spite of these limitations, our knowledge of ecosystem carbon cycling and its variability can be improved by investigating the global relationships between observed climate and plot-level biosphere-atmosphere fluxes across time scales from hours to years.
Complicating efforts to understand the relationship between flux and climate at multiple frequencies is the size of the FLUXNET database (Table 1), inherent gaps in the eddy covariance measurement records (Falge et al., 2001), and nonstationarity in flux time series Scanlon and Albertson, 2001). Extracting patterns from the FLUXNET database to investigate multi-scale relationships between flux and climate requires a methodology that is robust to data gaps and time-varying statistics that can simultaneously decrease the magnitude of flux time series in the time scale domain. Orthornormal wavelet transformation (OWT) with the Haar wavelet basis function  is ideally suited to mine (extract pattern from) flux data given: (i) the finite support and time locality of the Haar basis function for removing the effects of data gaps on spectral and cospectral calculations on nonstationary time series and (ii) the data reduction properties of OWT.
Our objective is to quantify the spectral properties of flux time series and the cospectral relationships among flux and climate in the FLUXNET database at time scales from hours to multiple years to investigate the degree to which biosphere-atmosphere C flux can be explained by instantaneous climate variability. This investigation contributes to major challenges in carbon cycle science including Ameri-Flux science question 3: "What is the causal link between climate and the exchanges of energy, CO 2 and water vapor for major vegetation types, and how do seasonal and interannual climate variability and anomalies influence fluxes?" Given previous findings in temperate forest ecosystems Richardson et al., 2008), we anticipate that the instantaneous effects of climate variability on flux variability will diminish at seasonal to interannual time scales, when biological changes that result from ecosystem development and ecosystem response to abrupt and slow variations in local climate dominate flux. Consequently, we hypothesize that the magnitude of the flux-climate cospectra will be lower at seasonal and annual time scales than daily time scales.
To accomplish our objectives and address the hypothesis, we quantify statistical differences among the spectra of net ecosystem exchange (NEE), gross ecosystem productivity (GEP) and ecosystem respiration (RE) using 999 siteyears of eddy covariance data from 253 ecosystems across the globe in the FLUXNET database using OWT . After exploring the inherent variability in NEE, GEP and RE across time scales, we quantify: (1) ecosystem spectral transfer (EST), defined here as the simple ratio between the orthonormal wavelet power spectra of ecosystem response variables (in this case NEE, GEP and RE), and the OWT spectra of different meteorological forcings, here photosynthetic photon flux density (PPFD), air temperature (T a ), vapor pressure deficit (VPD), and precipitation (P); (see (2) the wavelet co-spectra, which can be used to investigate the scale-wise correlations between climate and flux (see Fig. 1: Analysis II). In an appendix, we further investigate (3) 18 years of CANOAK model outputs to quantify the role of climate in driving interannual variability in NEE, GEP and RE using Fourier decomposition .
Analysis I assesses the frequencies at which an ecosystem response modulates a given meteorological forcing via Table 1. The number of FLUXNET research sites in version 2 of the LaThuille FLUXNET dataset with ecosystem type information, the total number of site years of data, the number of potential 1/2 h eddy covariance flux measurements, and resulting number of wavelet coefficients after orthonormal wavelet transformation (OWT), with respect to the climate and vegetation classes investigated here (after Williams et al., 2009 7  8  0  0  53  213  3 734 300  761  Tropical  1  1  0  9  0  1  0  3  0  15  49  858 960  214  Dry  0  1  0  1  1  3  0  4  0  10  27  473 280  153  Boreal  0  4  2  0  22  4  2  0  4  38  161  2 822 782  549  Arctic  0  0  2  0  0  1  0  0  3  6  22  385 678  86  Subtrop.-Medit.  5  6  11  5  17  11  2  6  0  63  272  4768551  941   Total  30  16  32  18  69  45  16  13  11  250  999  17 514  changes in state variables or functional parameters, and how the frequencies of this modulation shift among climate and ecosystem type. Analysis II is intended to unveil if scalewise correlations between ecosystem responses and climatic variables shift with respect to ecosystem type and climate regime to test the hypothesis. Analysis III explores low frequency flux variability. Results are discussed in the context of the hypotheses and the implications for multi-scale ecological modeling of the terrestrial carbon cycle (Williams et al., 2009).

FLUXNET database
Flux and meteorological data from version 2 of the LaThuile FLUXNET database (www.fluxdata.org) was used. Data were collected at individual sites according to network specific protocols (e.g., Aubinet et al., 2000), although deviations in methodology cannot be fully excluded for all sites. Half-hourly flux data were then processed according to the FLUXNET protocols. These include a filter for periods of insufficient friction velocity (u * ), the despiking of halfhourly flux data (Papale et al., 2006), and the partitioning of measured NEE into GEP and RE . FLUXNET products that fill the data gaps that result from missing measurements are not used in the present study, which investigates only measured flux and meteorological data that have passed quality control checks. The version 2 database includes 253 research sites encompassing 7 climate types and 11 vegetation types (Table 1) Agarwal et al., 2007). To obtain sufficient sample sizes of different ecosystem types for the statistical analysis, the vegetation classes "savanna" and "woody savanna" were combined to create one class called savanna, and the classes "open shrubland" and "closed shrubland" were combined to create one class called shrub. Data records extended from several months to over a decade. The database currently contains over 17.5 million half-hourly data points for each variable, and over 18 billion total cells of information.

Orthonormal Wavelet Transformation (OWT)
A brief, qualitative description of wavelet methodology is presented here; we refer the reader to Torrence and Compo (1998) for a basic treatment for geo-scientific applications, and Katul et al. (1995Katul et al. ( , 2001 for a detailed discussion of wavelet analysis for flux applications. Scanlon and Albertson (2001) and Stoy et al. (2005) present conceptual descriptions of wavelet techniques and further examples of wavelet analysis for biosphere-atmosphere flux research.
Wavelet decomposition differs from standard Fourier techniques in that it employs a finite basis function, called a "mother wavelet", that is translated (shifted) and dilated (expanded and contracted) across a signal to quantify, for the case of a time series, signal variance across both time and time scale. An infinite number of wavelet basis functions exist given the admissibility criteria that its integral is zero (Daubechies, 1992). The choice between multiple wavelet basis functions for a given application may appear subjective, but basis functions that are optimal for given time series properties can and should be selected (Torrence and Compo, 1998  Response to climatic variability also depends on ecosystem state through, for example, leaf area index (LAI) or biomass (B), both of which vary across time and space. This study introduces the ecosystem spectral transfer function (EST), defined as the ratio between the wavelet spectra of NEE and climatic drivers (Analysis I). A co-spectral analysis to identify scale-wise correlations between climate and flux is performed in Analysis II.
Haar wavelet basis function (a square wave) is the most logical choice given its localization in the temporal domain and consequent ability to control for the effects of the sharp discontinuities created by the inherent gaps (Falge et al., 2001;Katul et al., 2001). In this way, scale-wise contributions to the variance of the measured flux data can be quantified without considering the contributions of models to gapfill missing flux data, which are characterized by the frequency dynamics of their underlying driving variables. The wavelet transform can be continuous, but discretizations that avoid redundant information to return a tractable number of coefficients are advantageous. Orthonormal wavelet transformation (OWT) features a log-spaced discretization that quantifies the total variance of a time series of length 2 N in only N coefficients. For the case of FLUXNET data, the response between climatic drivers and the surface exchange of mass and energy for each measurement site across temporal scales can be expressed in terms of tens of wavelet coefficients that represent variability at different scales in time (the "time scale domain") rather than tens to hundreds of thousands of data points in the temporal domain Braswell et al., 2005;Stoy et al., 2005;. The wavelet cospectra between driver and flux can also be quantified to explore scale-wise climate-flux correlations Stoy et al., 2005) (see also Saito and Asunama, 2008).
We note that differences in the spectral estimate between Fourier and wavelet based methods are expected given that the two techniques decompose data using fundamentally different basis functions and algorithms. However, it should also be emphasized there is no "true" spectrum for complex, finite time series; Fourier decomposition assumes that the time series is composed of a combination of sinusoidal curves, which need not be the case. There is no reason why two spectra obtained with different kernels should be identical (Torrence and Compo, 1998). The only necessary requirement for spectral decomposition is conservation of spectral energy when summed across all frequencies (i.e. Parseval's Identity, Dunn and Morrison, 2005), which OWT satisfies. This is one reason why spectral scaling exponents inferred from OWT and Fourier methods agree reasonably well and are insensitive to the choice of the basis function when the time series exhibits an extensive scaling law across a wide range of scales (Katul and Parlange, 1994).
To compute OWT coefficients, all half-hourly measurements of NEE, GEP, RE, T a , VPD, PPFD, P and latent heat flux (LE) in the FLUXNET database were multiplied by their quality control flag (1 for raw data, 0 for missing data) to exclude missing or gapfilled measurements . This treatment ensured that data gaps do not contribute to the total series variance after Haar wavelet decomposition. All time series were then normalized to have zero mean (with zeroes in place of gaps) and unit variance for comparison. The time series for each site then underwent a standard "zeropadding" (Torrence and Compo, 1998) by adding zero values to both ends of the time series in order to make the length of each time series equal to a power of 2 for fast wavelet decomposition. The resulting zero-padded series was then again normalized to ensure that the normalized time series have unit variance and that all periods without valid measurements were set to zero. Differing lengths of the flux data records resulted in 15-18 dyadic scale representations per site and data record. The lowest-frequency wavelet coefficient for each normalized time series is poorly constrained and was dropped from the analysis, resulting in a maximum time scale of 2 17 half hours, or 7.48 years for the 16 year Harvard Forest time series.
For simplicity, orthonormal wavelet coefficients for all flux time series (NEE, GEP and RE) are abbreviated OWT FLUX , all meteorological time series are abbreviated OWT MET and some combination of flux and meteorological drivers are abbreviated OWT X . For the purposes of this investigation, LE is often considered alongside the meteorological drivers to investigate the coupling between carbon and water fluxes.
The ecosystem spectral transfer (EST) can be defined as The flux signal is said to be "amplified" ("dampened") compared to the respective climatic input if EST FLUX,MET is positive (negative) ( Fig. 1: Analysis I). We note that concepts and terminology from the signal processing literature, "amplifying", "dampening", "modulating", and "resonating" find a natural application when discussing eddy covariance measurements because the flux and meteorological variables are time-varying signals and ecosystems can be thought to process this signal in a corresponding, time-varying response ( Fig. 1). Amplification or dampening need not imply causality in systems, like ecosystems, that respond to a range of factors; rather, the EST is investigated to ascertain the simplest possible description of the instantaneous response of ecosystems to climate across scales of time.
To investigate how climatic inputs and flux outputs coresonate and to test the hypothesis, we quantified the wavelet covariance (Torrence and Compo, 1998) for all flux and meteorology combinations, and abbreviate the resulting coefficients OWT FLUX,MET Fig. 1: Analysis II). It is important to note that OWT NEE,MET and OWT GEP,MET will be negative if the relationship between flux and meteorological variability is positive due to the micrometeorological sign convention where flux from atmosphere to biosphere is denoted as negative.

Statistical analyses
A mixed model (PROC MIXED in SAS 9.1) was implemented to test for differences among OWT FLUX , EST FLUX,MET , and OWT FLUX,MET for different climate and vegetation classes for the 250 of 253 FLUXNET sites with climate and vegetation information. Wavelet variances were log-transformed prior to analysis when necessary to ensure normality of the response variables, noting that EST P. C. Stoy et al.: CO 2 flux and climate spectra is computed from log-transformed coefficients (Eq. 1). Coefficients at the 3.74 year and 7.48 year time scales were not included in the statistical analysis because the variable length of the flux data records resulted in a smaller sample size at low frequencies (see Fig. 2). Time scale, climate type, and PFT were specified as fixed main effects, and time scale×climate and time scale×vegetation were included as interaction effects. The vegetation×climate interaction effect was not included because of the large number of unreplicated combinations (Table 1). Observations for each site at different time scales were considered repeated measures, and a first-order autoregressive covariance structure was employed. Overall least-square means were calculated for the main effects climate and vegetation. The Tukey-Kramer multiple comparison test was used to evaluate overall differences among effect levels. We also tested for simple interaction effects among climate and vegetation ("sliced") per time scale using the PROC MIXED procedure in SAS 9.1. Test statistics that are significant at the 95% confidence level are reported.

Variability in NEE, GEP and RE
Medians and ranges of OWT FLUX for all FLUXNET sites are presented in Fig. 2. Abbreviations ("hourly", "daily", "annual", etc.) are used to approximate the discrete time scales for simplicity; for example, the annual coefficients represent variability at 0.935 years, daily coefficients represent 16 h variability as 24 is not a power of 2, and "multi-day" represents 1.33 day and 2.67 day variability because these time scales are longer than 1 day yet shorter than 1 week (Fig. 2, see also Table 2). These descriptions are meant to present a simple and intuitive means of communicating temporal scale for interpreting results and are used subsequently. It should be noted that sampling every third data point (1.5 h) would have resulted in OWT coefficients at exactly 24 h, but the annual variability would be poorly resolved given that coefficients at 256 and 512 days would also be returned.
Local reductions in OWT FLUX power spectra (the socalled "spectral gaps") exist at hourly time scales, and at multi-day to bi-monthly time scales. A low frequency spectral gap begins to appear at interannual time scales in OWT NEE and OWT GEP ; less so in OWT RE (Fig. 2). NEE and GEP variability were, on average, relatively lower at these time scales than the corresponding spectral peaks at daily and seasonal/annual time scales. The relative difference in OWT RE among sites increased at the longest time scale for which there were replicates (3.74 y), while inter-site differences in OWT NEE and OWT GEP decreased at low frequencies (Fig. 2). In other words, RE variability among sites was more different at interannual time scales, while seasonal and annual variability in NEE and GEP were more different among sites than RE. The single site with sufficient time series length to identify 7.48 y variability, Harvard Forest, had relatively high flux variability in GEP and RE at these "multiannual" time scales compared to the average 3.74 y interannual variability of other ecosystems.

Flux variability by climate zone and plant functional type
OWT FLUX (Fig. 2) was separated by climate and vegetation type (Fig. 3). Time scales for which there were significant differences among climate or vegetation types at the 95% confidence level, as determined by the multiple comparison tests, are listed in Table 2 and also demarcated by the horizontal bars in Fig. 3. (We note again that significance at the two lowest frequency time scales, 3.74 and 7.48 years, could not be determined due to lack of replicates.) There were no significant differences in OWT NEE among climate zones at hourly to weekly time scales, and among PFT at hourly (2 h) to bi-weekly time scales. This implies universality in the total relative variability, but not magnitude, of NEE at these time scales in the FLUXNET database. OWT GEP differed among climate types at bi-monthly to annual time scales, and among PFT at high-frequency (hourly) time scales and also monthly to seasonal time scales. High frequency (hourly to daily) OWT RE differed among climate types, but lower-frequency differences were only found at the seasonal time scale. PFT-related differences in OWT RE emerged at seasonal to interannual frequencies (Table 2, Fig. 3).

Analysis I: spectra of meteorological drivers and ecosystem responses
The median normalized power spectra of the flux and meteorological drivers varied by orders of magnitude across time scales (Fig. 4), highlighting the temporal dynamics of microclimate at characteristic frequencies and subsequent vegetation/ecosystem response via carbon flux Baldocchi et al., 2001b;Katul et al., 2001). Median OWT T a and OWT LE were notably energetic at seasonal and annual time scales. OWT P tended to be roughly equal across time scales larger than the isolated storm time scale, approximating white-noise, consistent with other studies Mahecha et al., 2008); P is known to possess a composite spectrum with multiple scaling laws across various frequencies (Gilman et al., 1963;Fraedrich and Larnder, 1993;Peters et al., 2002). The magnitude of OWT GEP , OWT RE and OWT LE at Harvard Forest at the 7.48 y time scale was large in comparison to FLUXNET median flux interannual variability at the 1.87 and 3.74 y time scales (Figs. 3 and 4).   Fig. 3. Mean wavelet spectra of the normalized net ecosystem exchange of CO 2 (OWT NEE ), gross ecosystem productivity (OWT GEP ), and ecosystem respiration (OWT RE ) per climate (left-hand panels) and vegetation type (right-hand panels) for the 250 eddy covariance measurement sites with ecosystem-level information in version 2 of the LaThuille FLUXNET database. Horizontal lines demarcate time scales at which there are significant differences among climate (top of subplot) or vegetation type (bottom of subplot) at the 95% significance level (see Table 2). Significance at the longest time scales (3.74 and 7.48 y) could not be determined among climate types due to lack of replicates. Circles at the 7.48 year time scales represent flux variability at Harvard Forest.
The tendency of ecosystems to dampen or amplify the variability of climatic drivers was explored via EST (Figs. 1 and 5). Figure 5 displays the median EST for each flux/meteorological driver comparison. At hourly time scales where flux footprint variability is known to dominate the measurement signal in spatially heterogeneous ecosystems , the variability of NEE tended to be greater than that of T a and VPD, but not PPFD. NEE and GEP were on average less variable than PPFD at time scales greater than an hour and shorter than a day, consistent with the observed saturating response of GEP to PPFD. NEE and GEP were damped with respect to T a and VPD variability at weekly and longer time scales, on average ( Fig. 5a and b), Table 2. Results from the mixed model analysis for significant interaction effects between meteorology (MET) by time scale (TS), and plant functional type (PFT) by TS on net ecosystem exchange (NEE) gross ecosystem productivity (GEP) and ecosystem respiration (RE) variability. Time scales for which spectra are significantly different at the 95% confidence level are denoted by the flux term abbreviation. These significant interaction effects are presented as bars in Fig. 3  but variability in RE was on average similar to or greater than T a variability at time scales of weeks and longer (Fig. 5c). Median OWT NEE was most similar to OWT PPFD across time scales as evidenced by the near-zero scale-wise EST NEE,PPFD (see also Fig. 4), but at lower (monthly to annual) frequencies, OWT GEP was most similar to OWT VPD (Figs. 4, 5a and 5b). OWT RE was most similar to OWT T a across time scales, but was also surprisingly similar to the median OWT LE spectra.
As a whole, EST GEP,MET followed similar scale-wise patterns to EST NEE,MET ( Fig. 5a and b), but median EST NEE,MET was dampened to a greater degree, on average, at lower frequencies. Conversely, EST RE,MET was amplified compared climatic variability, on average, at the lower frequencies.
EST FLUX,MET was significantly different among climate type and PFT across a complex array of time scales with distinct patterns, as demonstrated by the horizontal bars in Fig. 5, noting again that the bars in the upper sections of the subplot demarcate significant differences among climate type and the bars in the lower sections of the subplots signify significant differences among PFT. Slicing by climate type reveals significant differences at both short and long time scales in the relationship between all flux and meteorological variability. Significant differences among PFT emerged at longer (seasonal to annual) time scales in EST NEE,MET and EST GEP,MET , less so EST RE,MET .

Analysis II: wavelet cospectra between NEE and climatic drivers
The covariance between flux and climate was dominated by variability at daily and seasonal/annual time scales, but the magnitude of this covariance tended to be greater at daily rather than seasonal/annual time scales (Fig. 6). Mean OWT NEE,T a (i.e. the wavelet covariance between NEE and T a ) at the seasonal and annual time scales for different climate types were often of the same magnitude as those at the daily time scale, in contrast to OWT NEE,LE , OWT NEE,VPD , and OWT NEE,PPFD , which had a greater magnitude, on average, at daily time scales than seasonal or annual time scales (Fig. 6a). OWT GEP,LE was on average larger at seasonal/annual time scales than at daily time scales, and the median OWT RE,LE was unexpectedly large across monthly to annual frequencies.
Despite the relatively low OWT FLUX,MET magnitude at interannual time scales, these cospectra are in almost all cases significantly greater than that of a random signal at the 95% confidence level, as determined by a Monte Carlo analysis with 1000 synthetic autoregressive pink noise time series of the sort commonly applied in geoscientific studies (e.g., Allen and Smith (1996), see also Richardson et al. (2008), data not shown).
Most of the statistically-significant differences in OWT FLUX,MET among climate type and PFT were clustered around the daily and seasonal/annual spectral peaks, where  most of the energy in the cospectra resides (Fig. 6). No significant differences in the interaction effects among cospectra sliced by climate or PFT emerged at the weekly/monthly time scales. There were significant differences at the diurnal time scale when separating by both climate type and PFT for OWT NEE,MET and OWT GEP,MET , and significant differences at the seasonal and annual time scales among both climate type and PFT for all OWT NEE,MET and OWT GEP,MET except the PPFD cospectra. There were no PFT-related differences at hourly to seasonal time scales in OWT RE,T a ; significant climate and PFT-related differences for all OWT RE,MET emerged only at seasonal to annual time scales.

Multi-scale flux variability across global ecosystems
Spectral peaks at daily and seasonal/annual time scales are a general feature of the global NEE time series (Figs. 2  and 3). These peaks have been quantified previously in temperate forested ecosystems (Baldocchi et al., 2001b;Katul et al., 2001;Braswell et al., 2005;Stoy et al., 2005;Moffat et al., 2007; and further emphasize the dominant role of orbital motions at diurnal and seasonal/annual time scales in controlling the total variability in measured flux time series. Naturally, resolving such multi-scale variability alone does not necessarily translate to accurate estimation of long-term mean flux. These findings suggest that flux variability is concentrated in few time scales, i.e. the spectra are "unbalanced" . It has been previously argued that lowdimensional ecosystem models that capture these dominant modes of variability are a logical way forward for modeling flux at multiple time scales , but accurately modeling the seasonal and annual variability of fluxes has proven elusive for ecosystem models (Hanson et al., 2004;Siqueira et al., 2006), in part because of the shift from physical to biological control at longer time scales  as also evidenced by the PFT-related differences in OWT FLUX at longer time scales (Fig. 3). For example, a series of drought years can lead to a carry-over effect of reduced carbohydrate reserves for tree growth and a shift in carbon allocation in subsequent years, and incorporating these processes into ecological models is complex. Such dynamics were detailed at a semi-arid pine site in which NEE was much lower in a third year of drought than the previous years, suggesting a cumulative effect of ecosystem biology on biosphere-atmosphere flux .
It was found that the magnitude of the NEE and GEP spectra are not statistically different among climate zones or PFT until bi-weekly to monthly time scales are reached (Fig. 3). In other words, PFT is not a logical way to separate differences in the variability of flux time series at high frequencies in the FLUXNET database, at least when quantifying the total time series variance across scale rather than the changes in variance across scale and season. PFT does, however, separate the variability in observed NEE and GEP at the lower (monthly and longer) frequencies that are arguably more important for quantifying the role of terrestrial ecosystems in the global carbon cycle.
Climate and PFT-related differences in the RE spectra do not emerge until seasonal time scales (Fig. 3b). Given the ability of the eddy covariance system to resolve RE, PFT is not a logical way to separate RE variability across most of the time scales investigated here as expected. Because no statistically significant spectral differences emerge among PFT at higher frequencies, it may be argued that PFT is a scaledependent concept when considering the variability of flux in response to climate, and models may be simplified by ignoring vegetation-type distinctions at higher frequencies depending on the application and goals of the modeling analysis.
Regarding interannual and multi-annual flux variability, the 7.48 y variability in GEP and RE at Harvard Forest is greater than climatic variability (Fig. 4) and reflects the increase in flux magnitude over the course of measurements related to the species compositional shift toward red maple coupled with multi-annual recoveries from ecosystem disturbance (Urbanski et al., 2007). Such increases in variability at the low frequencies have been observed in other ecosystems P. C. Stoy et al.: CO 2 flux and climate spectra Mahecha et al., 2008), but a comprehensive synthesis of disturbance events on ecosystem C cycling requires more long term measurements from more ecosystems.

Ecosystem spectral transfer
Ecosystem CO 2 uptake strongly dampened the variability of most climatic drivers, on average (i.e. average EST NEE,MET and EST GEP,MET are less than zero) at time scales from weeks to years (Fig. 5), but ecosystem CO 2 loss was observed to be more variable than climatic drivers at these lower frequencies (i.e. average EST RE,MET >0). These results support a homeostatic mechanism in which ecosystem CO 2 uptake dampens instantaneous low frequency meteorological variability (see, for example, Odum, 1969), but low frequency amplification of climate by RE does not. Rather, climatic variability "excites" the RE spectra on average, which is consistent with the exponential transfer function commonly used to model RE variability at high frequencies and the strong spectral covariance between RE and T a at low frequencies (Figs. 5c and 6c). Also, significant differences in EST RE,MET emerged among different climate types, rather than PFT, across time scales (Fig. 5c). Together with findings from the previous subsection, these results agree with the notion that RE is determined by the effects of climate and moisture, more so than ecosystem type, on ecosystem C pools with different time scales of input and loss (Parton et al., 1987;Adair et al., 2008) despite the known coupling between GEP and RE (Högberg et al., 2001;Ryan and Law, 2005).

Cospectral relationships between flux and climate
Spectral peaks in NEE at longer (seasonal/annual) time scales are less-related to environmental drivers than the diurnal spectral peak (Fig. 6), which demonstrates a decreasing instantaneous influence between climate and NEE at longer time scales globally . Hence, it is no surprise that ecological models that use highly resolved PPFD variability (i.e. hourly), and to a lesser extent VPD variability, can explain the high-frequency energetic modes in NEE variability (Siqueira et al., 2006). At longer time scales, explaining flux variability via climatic variables alone is no longer effective (Siqueira et al., 2006;Stoy et al., 2008), in part because NEE and GEP tend to dampen climatic variability at low frequencies, supporting the homeostatic mechanism discussed in classic studies like Odum (1969). Results here agree with previous studies that have found a shift from exogenous to endogenous control of NEE in temperate deciduous broadleaf (DBF) and evergreen needleleaf (ENF) forests (Baldocchi et al., 2001b;Stoy et al., 2005;. Quantifying flux variability at longer time scales requires information on how ecosys-tems change in response to climatic variability, rather than how they merely respond to climatic variability.

Implications for ecological modelling
Significant climate and PFT differences in EST FLUX,MET appeared across more time scales than significant differences in OWT FLUX,MET , which are primarily restricted to diurnal and seasonal/annual time scales (Figs. 5 and 6). To the extent that changes in EST represent a change in the state variables or parameters that transfer climatic input to ecosystem output (Fig. 1), these state variables, parameters, and their variability are often different among PFT at bi-weekly to monthly time scales for NEE and GEP. Few PFT-related differences in EST RE,MET are significant, which again demonstrates that including PFT is less important for modeling multi-scale RE (see also Fig. 3). More ecosystem level ancillary data is required to determine how the shifts in transfer properties result from changes in ecosystem carbon stocks and/or the parameters of photosynthesis and respiration models Palmroth et al., 2005) to further investigate low frequency biological controls on flux . An emerging challenge is to model the multi-annual spectral energy of GEP and RE which, when combined via NEE, results in a low frequency dampening of carbon sequestration with respect to climatic variability.
The significant variation of EST across time scale indicates that the representation of ecosystems as simple instantaneous physical transformers, e.g. in simple soil-vegetationatmosphere models (SVATs) for photosynthesis, or time invariant regression equations for ecosystem respiration, is likely to fail at longer time scales (Carvalhais et al., 2008). This study confirms and generalizes former site-specific studies that have shown that biological variation in space and time is an important control of fluxes and may be empirically represented by changing biological ecosystem parameters that define the properties of the system Reichstein et al., 2003;Hibbard et al., 2005;Owen et al., 2007). In ecosystem process models, such temporal dynamics, which can change the response to climate variables across time-scales, may be represented by state variables (e.g. leaf area index, soil water content) and the appropriate representation of these state variables is pivotal for modeling NEE correctly across time scales, including pools of soil carbon with slow turnover times (Adair et al., 2008). Moreover, whether biophysical and ecophysiological parameters may be kept constant as is often the case in typical terrestrial biosphere models employed at global scale, or if parameters should be made temporally dynamic is a critical issue for modeling the land-surface C cycle (Curiel Yuste et al., 2004;Palmroth et al., 2005;Ryan and Law, 2005;Davidson and Janssens, 2006;Juang et al., 2008;Williams et al., 2009). Modeling or measuring rapid transients in endogenous variables (e.g. green-up of LAI in deciduous forests) may be critical for explaining flux variability over the stochastic weekly-P. C. Stoy et al.: CO 2 flux and climate spectra 2307 monthly spectral gap in some ecosystems . The importance of endogenous variables compared to climate drivers might well be one general reason why remote sensing based approaches that incorporate the biophysical state of the vegetation are comparably successful in quantifying terrestrial productivity, even after completely abandoning the use meteorological drivers (Jung et al., 2008).
A major motivation for the present study is to investigate the multi-scale coupling between climate, hydrology and carbon flux to identify potential improvements in ecosystem C models at the longer time scales at which they often fail (Hanson et al., 2004;Siqueira et al., 2006). The strong seasonal and annual OWT GEP,LE suggests that improving the representation of ecosystem hydrology and its coupling to CO 2 uptake is a logical step to improve models of CO 2 flux (Katul et al., 2003) in globally-distributed ecosystems. Interestingly, the mean seasonal/annual OWT GEP,VPD was of greater magnitude than OWT GEP,PPFD (Fig. 5b), further emphasizing the importance of hydrology in addition to radiation in determining the low frequency variability of photosynthesis.

Multi-scale coupling between GEP and RE
Quantifying the coupling between photosynthesis and ecosystem respiration has gained attention given that isotopic and ecosystem manipulation studies have consistently demonstrated a strong daily-to-monthly coupling between CO 2 input and output (Högberg et al., 2001;Barbour et al., 2005;Taneva et al., 2006). Despite the known limitation of the artificial correlation between GEP and RE in eddy covariance measurements (Vickers et al., 2009) [but see Lasslop et al. (2009)], some evidence for their coupling across time scales may be gleaned from a conservative investigation of their co-spectral properties.
The magnitude of the instantaneous GEP-RE covariance (OWT GEP,RE ) is large at seasonal to annual time scales; this relationship is on average as strong as the seasonal covariance between CO 2 uptake and water loss (Fig. 6). The GEP-RE covariance is also larger, on average, than the covariance between GEP and most meteorological drivers investigated, and is of the same average order of magnitude as OWT GEP,LE and OWT GEP,T a . GEP has been shown to explain 89% of the annual variance of RE across FLUXNET (Baldocchi, 2008) despite site-by-site differences in the strength of this relationship . Formally linking GEP and RE variability in models at lower frequencies is a logical step for improving model skill given the large low-frequency covariance between these fluxes, noting that the degree to which GEP and RE are coupled at low frequencies versus the degree to which they are controlled by similar environmental drivers across ecosystems remains unclear (Reichstein et al., 2007b). Future studies should investigate forthcoming FLUXNET data products that seek to minimize artificial correlations between GEP and RE (Lasslop et al., 2009). Interestingly, the average seasonal/annual covariance between RE and LE is as strong, or stronger, than the covariance between RE and T a . This may partly be explained by the tight coupling of photosynthesis, transpiration, and root respiration during the growing season (Drake et al., 2008). A strong correlation was observed between tree transpiration (sapflow) and root respiration measured on individual isolated trees at the semi-arid Metolius site, and soil respiration declined with soil water deficit as soil temperature continued to increase (Irvine et al., 2008) (see also Tang et al., 2005). The mechanisms for this relationship in the global database are less clear, and are likely linked through photosynthesis as well as soil moisture dynamics, which were not measured at a sufficient number of research sites to determine differences among climate zone and PFT. Measurements of soil moisture must be made to quantify the interaction between the terrestrial carbon and water cycles.
The covariance between GEP, RE and LE is consistently high at the seasonal/annual time scale, as is the covariance between GEP, RE and T a . This significant variability agrees with recent studies that demonstrate a statistically significant relationship between mean annual temperature and GPP and RE across European (Reichstein et al., 2007b) and Asian (Friedlingstein et al., 2006;Kato and Tang, 2008) ecosystems. In other words, despite the strong and wellcharacterized relationship between temperature and ecosystem respiration at short time scales (Lloyd and Taylor, 1994) (but see Janssens et al., 2001), warm years tend to be correlated with greater annual CO 2 uptake across FLUXNET (Delpierre et al., 2009;Kato and Tang, 2008;Reichstein et al., 2007b) except under instances of drought (Ciais et al., 2005;Reichstein et al., 2007a;Granier et al., 2007), and owing in part to the poor annual T a -RE relationships that are observed at selected sites Law et al., 2002;Stoy et al., 2008) (but see Reichstein et al., 2007b).

Conclusions and future studies
Eddy covariance measurement records are extending to time scales where they can be compared against independent measures of interannual variability and multi-annual variability in terrestrial ecosystem metabolism (Battle et al., 2000;Rocha et al., 2006;Rigozo et al., 2007). Multiannual eddy covariance measurement records are critical to understand the role played by ecosystems in controlling atmospheric CO 2 concentration (Houghton, 2000) and to contribute to global modeling analyses (Potter et al., 2005) for quantifying synchronous responses to disturbance events and climatic oscillations (Reichstein et al., 2007a;Qian et al., 2008) in an era of increasing climatic variability (Schär et al., 2004;Frich et al., 2002). EC measurements have nearly reached a point where predictions of classical ecological theories of long-term ecosystem function on the time scales of ecological succession (e.g. Odum, 1969;Luyssaert et al., 2008) can be quantified. Understanding these long-term dynamics will require additional information on long-term changes to ecosystem characteristics (e.g. species composition, Urbanski et al., 2007) and ecosystem C pools to rigorously test the predictions of both complicated ecosystem models (Williams et al., 2009) and simple ecological principles (Odum, 1969) (see, for example, Stoy et al., 2008).
The broader implications of this study suggest that, at shorter (hourly to weekly) time scales, the variability of C uptake across global ecosystems responds to climatic inputs in a similar matter regardless of PFT as revealed by the spectral analyses. PFT is a is a logical way to separate the variability in CO 2 uptake function at monthly to interannual time scales, but variability in ecosystem respiration was not clearly separated by PFT across most time scales (Figs. 3 and 5). The strong diurnal coupling between GPP and PPFD gave way to stronger covariance between the coupled C gain/water loss (via OWT GEP,LE ) and coupled C gain/C loss function of ecosystems at lower frequencies. These ecosystem-level mechanistic relationships between carbon and water cycling should be investigated in future studies of C flux variability, and modeling studies should seek to replicate these dynamics.

Modelling low frequency variability
To progress towards an understanding of the processes controlling multi-annual variability hinted at by the Harvard Forest time series, we investigated CANOAK model output for a DBF ecosystem driven by an 18-year meteorological record from Walker Branch, TN. CANOAK-modeled NEE, GEP and RE spectra demonstrated the expected seasonal and annual spectral peaks, and the GEP and RE spectra demonstrated an additional 7-11 year spectral peak (Fig. A1), recalling that GEP and RE were highly energetic at the 7.48 y time scale at Harvard Forest (Figs. 2 and 3) (note also Battle et al., 2000;Rocha et al., 2006;and Rigozo et al., 2007). CANOAK thus predicts multi-annual spectral peaks for a DBF ecosystem (Fig. A1) despite having a simple temperature-response function for a single soil C pool. Evidence for general multi-annual spectral peaks across biomes in the direct flux measurement record remains circumstantial given the paucity of direct observations to date. Until eddy covariance time series are extended, it remains unclear how vegetation and climate interact to create low-frequency flux oscillations, and the relative importance of low frequency oscillations to the carbon balance of global ecosystems (Rigozo et al., 2007;Qian et al., 2008). Uncovering the mechanisms for these multi-annual energetic events require longer time series coupled with careful process-level studies (Dengel et al., 2009).