A geostatistical synthesis study of factors affecting gross primary productivity in various ecosystems of North America

Abstract. A coupled Bayesian model selection and geostatistical regression modeling approach is adopted for empirical analysis of gross primary productivity (GPP) at six AmeriFlux sites, including the Kennedy Space Center Scrub Oak, Vaira Ranch, Tonzi Ranch, Blodgett Forest, Morgan Monroe State Forest, and Harvard Forest sites. The analysis is performed at a continuum of temporal scales ranging from daily to monthly, for a period of seven years. A total of 10 covariates representing environmental stimuli and indices of plant physiology are considered in explaining variations in GPP. Similarly to other statistical methods, the presented approach estimates regression coefficients and uncertainties associated with the covariates in a selected regression model. Unlike traditional regression methods, however, the approach also estimates the uncertainty associated with the selection of a single "best" model of GPP. In addition, the approach provides an enhanced understanding of how the importance of specific covariates changes with the examined timescale (i.e. temporal resolution). An examination of changes in the importance of specific covariates across timescales reveals thresholds above or below which covariates become important in explaining GPP. Results indicate that most sites (especially those with a stronger seasonal cycle) exhibit at least one prominent scaling threshold between the daily and 20-day temporal scales. This demonstrates that environmental variables that explain GPP at synoptic scales are different from those that capture its seasonality. At shorter time scales, radiation, temperature, and vapor pressure deficit exert the most significant influence on GPP at most examined sites. At coarser time scales, however, the importance of these covariates in explaining GPP declines. Overall, unique best models are identified at most sites at the daily scale, whereas multiple competing models are identified at longer time scales.


Introduction
Vegetation ecosystem dynamics are governed by processes operating over a wide range of spatio-temporal scales (Levin, 2003;Osmond, 1989).
A key challenge in understanding these dynamics lies in determining factors that govern gross primary productivity (GPP) and ecosystem respiration (R E ) at specific scales of interest.Over the years, understanding of vegetation dynamics that influence GPP and R E has improved considerably at small scales (Jarvis, 1995).However, the transfer of this knowledge to larger spatial and longer temporal scales remains difficult and fraught with potential errors, due to the spatio-temporal heterogeneity of vegetation dynamics (Gardner et al., 2001;Bradford and Reynolds, 2006).
These scaling issues also complicate modeling of responses to climate change, because the range of scales over which specific governing processes can be assumed to remain consistent is poorly understood (Jarvis, 1993).Consequently, biospheric model predictions have large uncertainties, in part because the importance of process-based functional relationships is not known at different spatio-temporal scales.To reduce this uncertainty, improved understanding of the varying influence of eco-climatic drivers on GPP and R E is required across a broad range of spatio-temporal scales.With respect to temporal scale, the sensitivities of GPP and/or R E to these drivers can be established by examining these cycles at flux tower sites (Ma et al., 2007;Pielke, 2000).
V. Yadav et al.: A geostatistical synthesis study of factors affecting gross primary productivity In general, longer records (typically ≥ 5 years) of continuous net ecosystem exchange (NEE) measurements are needed to assess the importance of eco-climatic drivers influencing GPP and R E at various time scales (Baldocchi et al., 2001a;Ito et al., 2005).These long term measurements provide a holistic perspective of an ecosystem, by measuring its response to a range of climatic conditions and extremes (Barr et al., 2007;Bradford and Reynolds, 2006;Gilmanov et al., 2006;Houghton, 2000).Presently, however, there are few flux sites where continuous NEE measurements are available for long time periods (Baldocchi, 2008).
Previous studies at longer time scales (>5 years) for these flux sites have focused primarily on examining GPP/R E associations to eco-climatic drivers at a single temporal scale, ranging from daily to inter-annual temporal variability (Barr et al., 2007;Dunn et al., 2007;Powell et al., 2006;Urbanski et al., 2007).For flux sites where less data is available, these relationships have mostly been assessed at daily or sub-daily scales (Schmid et al., 2000).The results from these studies indicate that different combinations of eco-climatic drivers explain variability of GPP/R E in different biomes and climate regimes (Ma et al., 2007).However, as most of these studies have been performed using data at a single temporal scale, the relative influence of these drivers on GPP/R E at both shorter (diurnal) and longer (monthly, inter-annual) temporal resolutions is poorly understood for all sites.As a consequence, the temporal scales at which given eco-climatic drivers are (or are not) important, and which drivers dominate at which scales, is not understood.
Although a few studies have examined the relationship of NEE, GPP, and R E to eco-climatic drivers at multiple temporal scales, their focus has mostly been on assessing correlations between fluxes and groups of eco-climatic drivers (e.g.climate versus vegetation) using either Fourier (Baldocchi et al., 2001b) or wavelet analysis (Stoy et al., 2005(Stoy et al., , 2009)).Though these methods suggest that relationships between NEE, GPP/R E , and these groups of drivers change as a function of scale, it is difficult to distinguish between spurious and causative implications of these correlations (See, for example, Granger Causality, Granger, 1969;Hacker and Hatemi, 2006) without applying more rigorous statistical testing (e.g.assessing the specific relationship, and its associated uncertainty, between a covariate and GPP/R E ).This creates a need for the application of more robust statistical inference methods for examining the relative importance of various drivers in explaining variations in GPP and/or R E at different time scales.One such statistical regression method, applied in Mueller et al. (2010), is extended in this study to understand the relationship between GPP and its eco-climatic drivers in different ecosystems across North America.
In addition to the studies described above, the increasing availability of eddy-covariance estimates of NEE from different regions and biomes has made it possible to develop canopy-level GPP models using empirical statistical methods (Makela et al., 2008).The most widely used statistical model for GPP is based on the concept of light use efficiency (LUE).LUE assumes a linear dependence of GPP on absorbed photosynthetically active radiation, modified by parabolic effects of temperature, and reductions in carbon uptake caused by increasing vapor pressure deficit and/or decreasing soil water content (Landsberg and Waring, 1997;McMurtrie et al., 1994).
The majority of past studies of factors affecting the temporal variability in GPP have relied on either a pre-specified parametric model, or on linear and/or nonlinear regression against a single factor or covariate (e.g.Mahadevan et al., 2008;Yuan et al., 2007).These relationships, though useful, do not give insight into the uncertainty associated with the form of the regression model itself and with the choice of covariates included in it.Moreover, in some cases, application of linear regression is not appropriate, because GPP observations (and/or regression residuals) are temporally autocorrelated, especially at shorter temporal resolutions (Mueller et al., 2010).To address these problems, it is desirable to include multiple covariates simultaneously in a regression framework, and to account for the correlation structure of the GPP observations in the statistical model.Three issues have to be taken into account while building such a model: the model needs to (1) account for any autocorrelation in the portion of the observations (i.e., GPP in this study) not explained by the ancillary variables, (2) identify a valid set of regressors/covariates to include in the regression model, and (3) account for multicollinearity among covariates.The first of these three issues requires the selection of a regression modeling framework that takes autocorrelation into account, whereas the last two can be accommodated by implementing an appropriate model selection scheme.Rigorous application of the model selection scheme in these situations also has the capacity to provide a metric to compare the importance of a given covariate across different regressions (or temporal scales, in this study).(Burnham and Anderson, 1998), and also to determine the relative importance of several covariates in explaining GPP.
In addition to the above, the approach presented and applied in this work provides a method for quantifying the uncertainty associated with the choice of the best set of covariates, when comparing competing regression models.As such, the analysis accounts for the possibility of there being multiple sets of covariates that provide comparable fits to the available GPP observations, which has not been considered in past studies at the AmeriFlux sites (Hui et al., 2003;Law et al., 2002).Finally, the presented analysis is performed across a continuum of temporal scales, ranging from daily to monthly, in order to identify specific timescales at which the importance of specific covariates changes.
The main implications of this work are, first, that it makes both methodological and theoretical advancements in the study of GPP and its explanatory covariates across temporal Second, an enhanced understanding of the environmental variables explaining GPP will help in better constraining the spatial and temporal trends of carbon flux in inverse models (e.g.Gourdji et al., 2008;Mueller et al., 2008).Lastly, improved knowledge of the dominant covariates, and changes in their importance as a function of temporal scale, provides an empirical foundation for improving mechanistic models of carbon flux.

Study sites
The six sites selected for the presented analysis are described in Table 1  The selection of this time period and these sites was based on two major criteria: (1) the sites should have more than five years of continuous GPP data, supplemented by other auxiliary environmental variables, and (2) for the same time period, the sites should have continuous gap filled remote sensing (RS) time series of MODIS-derived variables such as Leaf Area Index (LAI), Fraction of Photosynthetically Active Radiation (FPAR), Enhanced Vegetation Index (EVI), and Normalized Difference Vegetation Index (NDVI).These gap-filled data products were created as part of the MODIS for North American Carbon Program (NACP) project, and are available from http://ladsweb.nascom.nasa.gov/(Level 1 and Atmosphere Archive and Distribution System).Details about the algorithms used for creating these gap-filled products are presented in Gao et al. (2008) and Jonsson and Eklundh (2004).These two criteria ensure that the number of GPP observations is sufficiently large across the range of examined temporal scales to provide insights into the controls of GPP at different scales.Inclusion of RS data is important because, if RS indices can be used to explain the pattern of GPP observed at tower sites, then these indices can also be used to reduce the uncertainty associated with carbon flux estimates at unsampled locations.

GPP and auxiliary environmental data
AmeriFlux level four gap-filled GPP, available from the Carbon Dioxide Information Analysis Center (see, ftp: //cdiac.esd.ornl.gov/pub/AmeriFlux/data/)was analyzed in this study.For all six study sites, GPP obtained from the marginal distribution sampling (MDS; Reichstein et al., 2005) method was used.In this method, the Lloyd and Taylor (1994) regression model with constraints on soil or air temperature range is used to determine R E , which is then subtracted from NEE to obtain GPP.Hence, MDS-based GPP is not directly based on any prescribed functional form between environmental drivers and NEE.Moreover, it shows low annual bias, and is a recommended method for flux data synthesis studies (e.g.Moffat et al., 2007;Richardson et   2008).To further ascertain that no bias is introduced due to the choice of a particular GPP product for analysis, the results of the statistical modeling methodology described in Sect. 3 were repeated at all study sites using GPP derived from artificial neural networks (Papale and Valentini, 2003).No significant differences were found, confirming that the presented results are not dependent on the choice of MDS for estimating GPP.Supported by this observation and findings of Desai et al. (2008), who found that most methods of GPP estimation only differ by up to approximately 10%, GPP obtained from the MDS method was used in the results presented here, because it is available for all six study sites.
AmeriFlux level four gap filled data products also include data for six common covariates collected across flux sites.These six covariates, in conjunction with four MODIS measures of vegetation properties (phenology, density and light absorption), were used to model GPP.Additional data are available at some sites, including soil heat flux, wind speed, and others.However, to maintain consistency for cross-comparison of sites, we avoided inclusion of data that are not available across all sites.The ten examined covariates are therefore: (1) Composed LAI, (2) Composed FPAR, (3) Composed EVI, (4) Composed NDVI, (5) Vapor pressure deficit, (6) Global radiation, (7) Air temperature, (8) Precipitation, (9) Soil temperature, and (10) Soil water content.The first four of the above were obtained from MODIS (available from; http://ladsweb.nascom.nasa.gov/),and are documented in Huete et al. (1999), whereas details on the last six can be obtained from the Carbon Dioxide Information Analysis Center ftp site (ftp://cdiac.esd.ornl.gov/pub/AmeriFlux/data/). For algorithmic details on MODIS-based covariates, see Huete et al. (1994Huete et al. ( , 1999) ) and Myneni et al. (1997).Linear interpolation was used to obtain daily values from 8-day MODIS data products.All ten covariates were employed in the statistical modeling of GPP at five of the six sites; however, soil water content data were not available for Harvard Forest, and therefore only the nine remaining covariates were used for this site, as well as for jointly modeling the GPP of all study sites at the daily scale.In addition to these ten covariates, interactions among some of the covariates (e.g.Air Temperature × Vapor pressure deficit, EVI × Soil Moisture × Vapor pressure deficit) were also initially considered for modeling GPP.However, these combined variables were never found to be significant in explaining GPP using the methods outlined in Sect.3, and were therefore excluded from the final analysis.
The majority of the analysis was conducted at daily, 8-day, and monthly time scales, as these are the most commonly utilized scales for examining trends in GPP.However, for understanding changes in the importance of covariates in explaining GPP across a continuum of temporal scales, some results are also presented for temporal resolutions ranging from one day to 30 days, in one-day increments.

Statistical methodology
Regression modeling has been adopted as the primary method for understanding the relationship between GPP and its eco-climatic drivers at AmeriFlux sites (see for e.g., Hui et al., 2003;Powell et al., 2006;Urbanski et al., 2007).However, the traditional multiple linear regression (MLR) approach does not account for the temporal correlation observed in regression residuals in flux tower analyses (Law et al., 2002), which can lead to an underestimation of the uncertainty associated with regression coefficients (Chatfield, 2003), which, in turn, can make certain ancillary variables erroneously appear to be significant.In addition, if there are gaps within the time series (e.g.examining seasonal GPP observations over multiple years), then regression coefficients may also be biased (Cressie, 1993;Hoeting et al., 2006).Alternately, the geostatistical regression approach presented in Mueller et al. (2010), and extended here, quantifies and accounts for the temporal correlation of regression residuals of GPP time series.Software implementing the approach described here is available at http://puorg.engin.umich.edu.

Geostatistical regression
Geostatistical regression (GR, e.g.Erickson et al., 2005) is a parametric method used to model environmental phenomena, such as GPP, that are correlated in space and/or time.In GR, a variable such as GPP is expressed as the sum of a deterministic component that represents the portion of the signal that can be explained by a set of covariates, and a stochastic component modeled as a spatially and/or temporally autocorrelated random function.Mathematically, a GR model can be expressed as: where z (n×1) are observations of the parameter of interest (GPP in this study) at specific locations and/or times, X (n×p) is a pre-specified design matrix of covariates, β(p×1) are the coefficients relating individual covariates in X to the dependent variable z, and ε are zero-mean intrinsically stationary residuals (Cressie, 1993).If in Eq. ( 1) the residuals are uncorrelated in space and time, then inference about the relationship between X and z can be derived from an MLR model.The residuals ε in this case would be independent, and their covariance would be a scaled identity matrix.In this work, z are GPP observations at a given flux tower site as a function of time, and X contains a subset of the environmental variables described in Sect.2.2, defined at the same times as the GPP observations.The approach for selecting the subset(s) of environmental variables is described in Sect.3.2.A covariance function is used in GR to model the correlation structure of the stochastic component ε.In this study, an exponential covariance function combined with a nugget effect is used, as such a model has been shown to represent the variability of NEE and GPP in the past studies (Gourdji et al., 2008;Michalak et al., 2004;Mueller et al., 2010).This function is defined as: where Q is the (n×n) covariance matrix of the regression residuals, σ 2 N is the variance of the variability that is uncorrelated in time (this can include measurement error and/or the variability at time scales below the averaging time used for GPP observations), h is the (n×n) matrix of time lags between GPP observations, σ 2 S is the variance of the variability that is temporally correlated, l is the correlation range parameter, such that 3l represents the time after which the correlation between GPP residuals becomes negligible, and δ h=0 is the Kronecker delta function, equal to one when the time lag is equal to zero (i.e. on the diagonal of the matrix Q) and zero otherwise.Equation (2) assumes that the GPP residuals (ε) are homoscedastic.Although a more complicated model could be implemented to include heteroscedasticity, we do not expect the variance to change significantly within the growing season for a particular site.
Restricted maximum likelihood (e.g., Kitanidis, 1995;Michalak et al., 2004) is used to optimize the parameters σ 2 N , σ 2 S and l.For GR, the restricted maximum likelihood objective function can be written as: where L(σ 2 N ,σ 2 S ,l;z) is minimized with respect to σ 2 N , σ 2 S and l in Q, and || denotes the matrix determinant.The reader is referred to Michalak et al. ( 2004) and Mueller et al. (2010) for a more detailed description of this approach in the context of NEE and GPP/R E modeling.After obtaining Q, best estimates of the GR coefficients and their associated uncertainties are obtained as: where β are the (p×1) best estimates of the GR regression coefficients β, and V β is the (p × p) covariance matrix representing their uncertainties and cross-covariances.The V β as given in Eq. ( 5) is further used to compute uncertainties of the contribution of the covariates X β for explaining the observed GPP, which is defined as: where the diagonal elements of V X β (n × n) represent the uncertainties of X β. the relationship of a particular covariate with GPP within the GR model, the correlation coefficients of the estimated regression coefficients are quantified as: where W(p ×p) is a diagonal matrix of the square root of the diagonal entries in V β .The performance of GR is evaluated by computing the coefficient of determination (R 2 ) between X β and z, termed "variance reduction" in this study.

Model selection using the Bayesian Information Criterion
Model selection in GR is the process of selecting covariates with the goal of creating a regression model of optimal complexity.Although scientific understanding should form the primary basis in selecting covariates, existing theory provides limited guidance with regard to which eco-climatic covariates should be included in explaining GPP at different sites and time scales of interest (for example see different formulations of LUE models (Yuan et al., 2007)).Hence, the need arises for a statistical method that can identify regression models of optimal complexity.The approach implemented here, as in Mueller et al. (2010), is based on the Bayes Information Criterion (BIC, Schwarz (1978)).
Relative to other methods of regression model selection, such as Akaike's information criterion (AIC, Akaike (1974)), BIC is more informative in an inferential framework.In addition, BIC can also be used to determine the relative likelihood of several alternative regression models, and the relative importance of specific covariates, in terms of their posterior probabilities.For the GR models in this study, this is achieved by comparing BIC values of all possible linear combinations of the ten or nine covariates at a particular flux site.This is in contrast to more traditional hypothesis-based model selection techniques that are based on F-tests, and that can only be used to compare nested statistical models, and that may therefore overlook covariates that are jointly, but not individually, significant in explaining GPP.
The process of model selection also helps in eliminating any highly-correlated covariates that provide redundant information.The presence of correlated covariates in a regression model increases the size and complexity of the model, but does not significantly improve its fit.BIC accounts for this correlation, and therefore does not select highly correlated covariates, thereby reducing the problem of multicollinearity.
Bayesian model selection is based on the idea that candidate models should be compared in terms of prior and posterior evidence for a model over an alternative model.This idea is expressed in terms of Bayes factors that can be used to compare several nested and non-nested models simultaneously, where a set of nested models is one where one model is a subset of the other.Mathematically the standard form of the Bayes factor B(z) can be written as: where p(M 1 |z)and p(M 2 |z) represent posterior probabilities of models M 1 and M 2 given the available measurements z and, p(M 1 ) and p(M 2 ) represent the prior probabilities of the two models, which are assumed even in this study (i.e., p(M 1 ) = p(M 2 )= 1/2), because no prior information is available in favor of any model.B(z) quantifies the relative support for various models, and does not lead to simple accept and reject decisions as in hypothesis-based methods.Jeffreys (1961) proposed a scale for interpreting Bayes factors, suggesting that if B(z) > 1 then the measurements favor M 1 over M 2 , and when B(z) < 1 they favor M 2 .He also suggested the following grades of evidence for B(z): when 1 ≤ B(z) ≤ 3 there is "very weak evidence" for M 1 over M 2 , when 3 ≤ B(z) ≤ 10 the evidence is "positive", when 10 ≤ B(z) ≤ 100 it is "strong", and when B(z) > 100 it is "decisive."Unfortunately, finding Bayes factors involves evaluating the integrals of the likelihoods of different models, which is difficult to do analytically.An approximation to this approach, based on a maximum log-likelihood estimator, is given by the BIC (Schwarz, 1978): where n is the number of observations, k is the number of covariates and L * σ 2 N ,σ 2 S ,l;z is the log-likelihood of the model under consideration.For GR, assuming that the residuals ε in Eq. (1) follow a Gaussian distribution, L * σ 2 N ,σ 2 S ,l;z can be expressed as (Mueller et al., 2010): In this study, all possible sets of covariates for GPP, at all sites and timescales of interest, are compared by computing their BIC, with the set of covariates corresponding to the lowest BIC being identified as the "best" GR model.

Quantifying and accounting for model uncertainty
Although the model with the lowest BIC can be thought of as the "best" GR model, multiple models could potentially explain the observed GPP to a similar degree, especially given the similarity among some of the candidate covariates (Sect.2.2).Hence, in these circumstances, it is more suitable to select a set of candidate models from all possible GR models.A Bayesian solution for selecting these candidate models was proposed by Leamer (1978), and implemented in a BIC framework by Raftery (1995).In this study, we implement a modification of this method as described below.
When comparing a group of regression models, M = {M 1 ,M 2 , M 3 ,..., M N }, the posterior probability of a particular model M i from the group of N possible models in M can be given as (Kass and Raftery, 1995): where the a priori probability of all the models in the group is equal in this study, i.e., p(M 1 )=p(M N )=1/N , and p(z|M i )∝ exp − 1 2 BIC M i .For further discussion see Hoeting et al. (2000), Kass and Raftery (1995) and Wasserman (2000).Hence: However, when the total number of possible regression models N is large, the posterior probability as given in Eq. ( 12) would be small even for the best regression model chosen from BIC.The number N is particularly large when the number of candidate covariates potentially explaining the response variable is large, and M has regression models formulated from all possible subsets of these individual covariates.In this situation, M includes many models that have no basis of support.Hence, it is desirable to reduce the size of M before computing p(M i |z) in Eq. ( 12).In this study, a two step approach was adopted to reduce the number of candidate models.First, GR models that are poorly supported by the available data were eliminated if they were at least 20 times less likely (analogous to a 0.05 significance level) than the "best" model as determined using the BIC (Sect.3.1).This is equivalent to removing all GR models that have a BIC value of at least 6 above than that of the "best" model (i.e., BIC M i − BIC bestmodel > 6), and corresponds to "strong" evidence in favor of "best" model on the Jeffreys (1961) scale described in Sect.3.2 (i.e., 2 log (Bayes factor); 2 log (20) ∼ = 6).Second, the number of candidate models in M was further reduced by comparing the log-likelihoods of any remaining nested models as given in Eq. ( 10) (Burnham and Anderson, 1998), and eliminating any models that are not significantly better than a model consisting of a subset of their covariates, at a 0.05 significance level based on a chi-squared test.The subset of models remaining after these two preliminary steps was considered to be that that could adequately explain GPP for a given site and time scale.In the analyses performed here, the number of remaining models N ranged from 1 to 6 (Table 2) across the examined sites and time scales (see Sect. 4.2).The posterior probability i.e., p(M i |z) of these models was then computed from Eq. ( 12) to assess their relative probabilities.Although the "best" model as identified using the BIC criterion alone (Sect.3.2) always has the highest probability of being the "correct" model, this probability will not be 100% unless only one model remains in the selected subset (N = 1), and may in some cases will be below 50% (see Sect. 4.2)

Quantifying importance of individual covariates
In addition to identifying the set of GR models that explain GPP and the posterior probability associated with each of these models, we were also interested in finding the relative importance of including a specific covariate in the full set of possible GR models.In simple terms, this importance represents the posterior probability of a covariate being included in a particular set of models M. Mathematically, given the superset of covariates X={X 1 , X 2 , X 3 ...X p }, where p = 9 or 10 for the sites examined in this study (Sect.2.2), the importance of a covariate can be expressed as (Raftery et al., 1997): where p(X j |z) is the posterior probability of the j th covariate, and p(M i |z) is calculated as in Eq. ( 12).This quantity indicates the "importance" of a particular covariate, and is equal to the sum of p(M i |z) for all possible subsets of covariates that include X j .The posterior probability of the covariates in Eq. ( 13) is computed without reducing the set of possible models using the procedure described in Sect.3.3.This metric provides an indication of the relative importance of specific covariates for a particular time scale and site.According to Raftery (1995), p(X j |z) > 0.50, 0.75, 0.95 and 0.99 indicates "weak", "positive", "strong", and "very strong" importance, respectively.Similar probabilityevidence strength relationships were also suggested with regard to model probabilities, p(M i |z) in Eq. ( 12).In Sect.4, we use these probability thresholds to discuss both model probabilities, i.e., p(M i |z) and the importance of individual covariates, i.e., p(X j |z), in explaining GPP.
With regard to GPP, we posit that covariates that explain GPP vary in their importance across temporal scales.Therefore, the importance of all covariates was calculated for all time scales ranging from one-day to 30-day intervals.This is achieved by averaging daily GPP and auxiliary eco-climatic data time series to these 30 different time scales.The variability of p(X j |z) from 1-day to 30-day time scales illustrates how the importance of specific covariates in explaining GPP evolves as the examined time scale increases, and provides a key foundation for understanding how processes governing GPP vary with temporal scale.

Results and discussion
Results are presented primarily at monthly, 8-day and daily temporal scales.Some of the detailed supporting findings not www.biogeosciences.net/7/2655/2010/Biogeosciences, 7, 2655-2671, 2010 Table 2. Final number N of geostatistical regression (GR) models considered, and the posterior probability of the "best" model (see Eq. 11 and subsequent text) at the daily, 8-day and monthly scales.The posterior probability of the covariates (see Eq. 12 and Sect.3.4) is shown for all covariates.The covariates that are included in the "best" GR model are shaded, with green indicating a positive correlation with gross primary productivity (i.e.βi >0) and red indicating a negative correlation (i.e.βi <0).See Tables A1-A3 for further details.
directly related to the main thesis of this work are presented in the supplementary material.

Relationship between covariates and observed GPP
Examination of the covariates selected for the "best" model for each site and time scale, based on BIC analysis, as well as of the covariates selected to model GPP when pooling data from all sites at the daily scale (Table 2), provides an indication of the factors that explain GPP at different time scales across ecosystems.At the daily scale, Global radiation, Air temperature, Vapor pressure deficit, and Precipitation are the most common variables selected across sites.The selection of these covariates is not surprising, and reconfirms the light use efficiency approach adopted in earlier studies.At the monthly scale, however, the covariates used in light use efficiency models do not perform well, because several competing covariates can explain the seasonality of GPP, and the variables that best explain this seasonality vary by site.Among MODIS-based variables, EVI or NDVI perform better in explaining GPP at shorter temporal scales, whereas at the monthly scale this demarcation is not clear, with LAI, FPAR and EVI each being selected for at least one site in the "best" GR model.Overall, however, there is no mechanistic reason for which GPP at one site would be better explained by EVI and another by NDVI or LAI, as all three are indicative of changes in phenology.The observed differences are simply indicative of better "statistical performance" of a particular vegetation index in explaining the phenology at a particular site.For the examined sites, BIC does not select more than one of these four variables for a given site at most scales (Table 2), confirming that these variables are providing similar (and therefore redundant) information.
To understand the relationship between selected covariates and GPP, the sign of the regression coefficients β was examined for the "best" GR model identified using BIC at the daily scale.A positive sign on βi indicates a positive correlation with GPP (i.e., carbon uptake), whereas a negative sign indicates a negative correlation (i.e., a reduction in carbon uptake).The average individual and total contribution of these regressors to GPP for the "best" GR model was also analyzed at the daily scale, and is shown in Fig. 2. Since the correlations ρ β (Eq.7) between the uncertainties of the regression coefficients for Air temperature and Vapor pressure deficit were greater than 0.70 at the Tonzi Ranch and Blodgett Forest sites, and for Global radiation and Air temperature at the Harvard Forest site, the combined (rather than individual) contributions of these covariates are shown in Fig. 2. Standardized regression coefficients for the "best" GR model at daily, 8-day and monthly temporal scales are presented in Tables A-1, A-2, and A-3 in the supplementary material.
In general, the signs of the regression coefficients are consistent across sites and ecosystems.At all sites, LAI, EVI, NDVI, Global radiation, and Air temperature are positively correlated with GPP, whereas Precipitation and FPAR are negatively correlated with GPP (Table 2).This indicates that the overall impact of these variables on GPP is consistent across ecosystems and examined sites.The influence of these variables on photosynthesis has been well recognized in earlier studies (Gourdji et al., 2008;Iio et al., 2004;Powell et al., 2006;Schmid et al., 2000;Sims et al., 2008;Davidson et al., 2000;Day, 2000).The sign of the coefficient for Vapor pressure deficit, on the other hand, varied across sites.Vapor pressure deficit can be an indicator of water stress, and, as Vapor pressure deficit rises, stomatal conductance declines, which reduces the photosynthetic rate and the efficiency of plants to use light to fix carbon.This yields a negative regression coefficient for VPD at most sites.The Vapor pressure deficit at which photosynthetic rate gets suppressed varies by plant species, and depends on the availability of water.At the Kennedy Space Center site, moisture stress is insignificant, and Vapor pressure deficit is mostly below 1.5 kilopascals.As Vapor pressure deficit is not an inhibitor of photosynthesis at Kennedy Space Center, it is positively associated with GPP, and serves primarily to capture the seasonal cycle of GPP at the monthly scale.
Among the ten candidate covariates, the largest contribution to carbon uptake is associated with EVI or NDVI at the Kennedy Space Center, Vaira Ranch, Tonzi Ranch, Morgan Monroe and Harvard Forest sites, with Global radiation accounting for the second largest contribution at all sites (Fig. 2) except Morgan Monroe.Conversely, Vapor pressure deficit and Precipitation are associated with minor reductions in uptake at most sites.For the remaining covariates, the strength of contribution varies across sites.From a broader perspective, the analysis of the size of the contribution brings to the forefront the fact that light and indicators of vegetation phenology/density are most strongly associated with carbon uptake, whereas moisture stress, indicated by Vapor pressure deficit, plays a minor role in reducing GPP.
This conclusion is further confirmed when GPP from all six study sites are jointly regressed against the covariates described in Sect.2.2 at the daily scale (Table 2, first line).As part of this joint regression, EVI, Global radiation, Air temperature, Soil temperature, Vapor pressure deficit and Precipitation were selected as explaining GPP across all examined sites.Out of these six covariates, the first three were positively correlated with GPP, and the remaining three were negatively correlated with GPP.As in the site-level analysis, the largest contribution to carbon uptake is associated with EVI, followed by Global radiation (Fig. 2).

Model performance
Overall, the ability of the selected variables to explain observed variability in GPP improves at coarser temporal scales (Table 3).Even at the shortest examined temporal resolution (daily scale), however, a substantial amount of variability in GPP can be explained by the selected covariates.For example, the GR of daily GPP across all six sites had an R 2 of 0.72.For GR at individual sites, Fig. 2. Total average contribution of each covariate included in the best GR model to GPP estimated at the daily scale for the examined sites.The total contribution is calculated by averaging the portion of GPP predicted from each covariate (X βi ) over the entire examined period.Numbers in bold for each site represent the average daily predicted GPP over the examined period (See Table 1 for the length of the annual growing season and the examined time period for each site).
the variance reduced is highest and the root mean squared error of the "best" GR model is lowest at the monthly scale (Table 3).This is to be expected, as GPP time series at finer resolutions have more variability and noise.Due to this higher variability, more covariates are required to explain GPP at the daily scale (Table 2) relative to the monthly and 8-day scales.GPP observations and GPP predicted (X β) by the "best" GR model are also presented at the 8-day scale for each site in Fig. 3.The one standard deviation uncertainty bounds in Fig. 3 have been estimated as the square root of the sum of the uncertainty variance of X β (Eq.6), and the estimated variance of the residuals (σ 2 N + σ 2 S , Eqs. 1 and 2).Results for other time scales are qualitatively similar.
Comparatively, across sites and scales, the selected "best" GR models of GPP are more effective in explaining GPP variability in ecosystems that have a definite seasonal cycle, including the Harvard Forest, Morgan Monroe, Tonzi Ranch and Vaira Ranch sites.Conversely, the variance reduced is lower at the Blodgett Forest and Kennedy Space Center sites, where the growing season is longer in comparison to the other study sites (Table 1).Ecosystems exhibiting little seasonal variability, such as the Kennedy Space Center site, are also found to have shorter temporal correlation lengths (3 l) in their regression residuals, which imply that GPP observation residuals at shorter temporal lags are independent of one another.Absence of strong seasonality at these sites also means that phenological factors do not exert dominant controls, and GPP is primarily governed by processes operating at shorter time scales, some of which www.biogeosciences.net/7/2655/2010/Biogeosciences, 7, 2655-2671, 2010  are identified through proxy covariates (see Sect. 4.1).In addition, some of the explanatory environmental variables omitted due to the cross-site data consistency restrictions imposed in this study could also play an important role in improving the GR model performance at both these sites.Some examples of these additional covariates include aerosols and albedo, which impact radiation and friction velocity, and, in turn, GPP (Chambers et al., 2004;Oliveira et Fig. 4. Posterior probability of covariates at 1-day to 30-day temporal scales for examined AmeriFlux study sites (see Sect. 3.4 and Eq.12).The probabilities across temporal resolutions have been smoothed using a three day running average, and only covariates selected in the best GR model at the daily, 8-day and/or monthly scale are presented for each site.al., 2007).Deficiencies in the candidate variables available to the GR model are also possible for the other examined sites, but any such missing variables had a lesser impact on the ability of the model to reproduce the observed variability.
At the monthly and 8-day scales, several competing candidate models that can explain GPP are identified at most sites (Table 2), because more than one subset of covariates explains GPP to a similar extent.Although the numbers of plausible models is higher at longer time scales, several common covariates are observed across these plausible models.These are generally variables that have "positive" to "very strong", importance, as defined in Sect.3.4.Hence, although the number of candidate models selected is larger at coarser temporal scales, the dominant covariates can still be identified through their posterior probabilities, as shown in Table 2.
The variance reduced and the posterior probability (Eq.12) of the GR models at the daily, 8-day and monthly temporal scales highlight the problems that can be encountered in the traditional statistical modeling of GPP that relies on studying relationships based on a single "best" model.At shorter temporal scales, the "best" GR models have higher posterior probability (Table 2) but explain less variance in the GPP signal (Table 3).At coarser temporal scales, the variance reduced for the "best" GR model is higher, but several competing GR models can explain the seasonality of GPP, which reduces the posterior probability of the single "best" GR model (Table 2).This clearly shows that using a single "best" statistical model of GPP without accounting for its uncertainty relative to other competing models can yield misleading conclusions about the importance of specific covariates in explaining GPP.The application of the method presented here provides a more objective basis for evaluating the uniqueness of the selected GR model of GPP, by computing its uncertainty within a Bayesian model selection framework.

Importance of covariates across continuum of temporal scales
Examining the relative importance of specific covariates across multiple temporal scales provides an indication of their role in explaining GPP as a function of temporal resolution.This is illustrated in Table 3, which shows the posterior probability (Eq.13) of specific covariates at the examined sites at daily, 8-day and monthly time scales.
For example, at Harvard Forest, the importance of Global radiation is "very strong" at the daily scale, but this variable is not significant at the monthly scale.Similarly, at the Tonzi Ranch site, Precipitation is very important at the daily scale, but is not significant at the 8-day scale.Analogous examples of varying importance were found at all sites.These differences in the posterior probability of covariates indicate that the importance of environmental variables in explaining GPP at one temporal scale does not imply similar importance at another scale.As a result, it becomes necessary to examine these relationships at a continuum of temporal scales, in order to highlight (1) covariates that are important at all scales, and consequently are dominant explanatory variables governing GPP in a particular ecosystem, and (2) as suggested by Wu and Li (2006), time scale thresholds below or above which certain variables become important in explaining observed variability in GPP, and thus reflect fundamental shifts in controlling factors or processes across scales.
To study the change in importance of covariates in explaining GPP at a continuum of time scales, the posterior probability (Eq.13) of covariates was computed at 1-day to 30-day time scales.Results reveal that most covariates vary in their importance as a function of temporal resolution, and that few cross-scale or scale invariant explanatory covariates of GPP are present at all examined sites.Fig. 4 depicts changes in the posterior probability, or importance, of individual covariates from the daily to the monthly scale, and includes all covariates that were selected in the daily, 8-day, and/or monthly analyses discussed previously for each site.Any covariate that has greater than "positive" importance (posterior probability greater than 0.75) at most time scales in Fig. 4 is considered to have a scale-invariant importance.
Among the ten variables examined in this work, the only covariates exhibiting scale-invariant importance are LAI and Vapor pressure deficit for the Vaira Ranch site, Global radiation for the Tonzi Ranch and Blodgett Forest sites, and EVI for the Kennedy Space Center and Morgan Monroe sites.Therefore, the fact that a given model or covariate is an important predictor of GPP at a particular time scale does not, in general, imply that this is also true at other scales.In earlier studies, some of the covariates included in this study have been found to significantly explain variability of GPP at a particular temporal scale.But in the absence of a multi-temporal scale analysis, such as the one presented here, it is difficult to judge whether the strength of modeled relationships holds for other temporal scales, or is specific to a particular examined scale (Baldocchi et al., 2004(Baldocchi et al., , 2005;;Goldstein et al., 2000;Lee et al., 2002;Ma et al., 2007;Urbanski et al., 2007;Wang et al., 2008;Xu and Baldocchi, 2004;Xu et al., 2001;Powell et al., 2006;McMillan et al., 2008;Saito et al., 2009).
Most sites exhibit at least one temporal scale threshold, above or below which several covariates become significant in explaining GPP.All six study sites exhibit at least one prominent scaling threshold below the 20-day scale, indicating that environmental variables and processes that are important in explaining GPP variability at synoptic scales are different from those operating at larger temporal resolutions (Fig. 4).For instance, at the Vaira Ranch site, Air temperature, Global radiation and Precipitation are important at the synoptic scale, but their influence declines gradually beyond the 15-day scale.As expected, these thresholds vary both with sites and covariates.For example, at the Harvard Forest site, the scaling thresholds for different covariates lie between the 1-day and 15-day temporal resolutions.Multiple temporal scale thresholds are also observed in some cases, such as for Blodgett Forest site, where shifts in the importance of Soil temperature and Precipitation are observed between the 1-day and 5-day scales, whereas the importance of Global radiation and Vapor pressure deficit changes between the 20-day and 25-day scales.For some sites, such as Kennedy Space Center, no clear overall temporal scale thresholds are observed, although the importance of individual variables still changed more gradually across temporal scales.
Beyond well defined time scale thresholds, the importance of individual covariates can also be assessed, which helps in identifying different types of scaling behaviors, and therefore may reveal relationships between different time scales in an ecosystem (for further discussion see, Ehleringer and Field, 1993;Jarvis, 1995).The importance of individual variables across temporal scales can generally be classified into four groups, which are defined here as: (1) Scale invariant dominance (2) Declining dominance from daily to coarser scales, (3) Emergent dominance from daily to coarser scales, and (4) Varying dominance across scales.
Scale invariant dominance was discussed in the previous paragraphs.Declining dominance is common across sites, indicating that some covariates are only important in explaining GPP at finer temporal resolutions.For example, the importance of Precipitation and Global radiation decreases gradually from the daily to the monthly scale at the Kennedy Space Center site.This kind of scaling behavior is commonly observed in ecology, and represents processes that are only important at finer scales, and do not lead to any significant understanding of the pattern (GPP in this case) observed at coarser scales (cf.Levin, 1992).Emergent dominance is less common, but is observed, for example, in the case of LAI at the Tonzi Ranch site, and LAI and NDVI at the Harvard Forest site.Minor increases in importance are also observed for some covariates at other sites.This behavior can be caused by physiological processes that decouple systems from their primary controls at finer scales by creating a web of indirect interactions (Greig-Smith, 1979;Woodward, 1987).However, at coarser scales, these interactions no longer remain important, and other processes that govern a system at these scales become dominant (Levin, 1989;Wiens, 1989).The last scaling behavior, that of varying dominance across scales, also occurs at many sites.For example, at the Morgan Monroe site, the importance of Global radiation is high up to the 5-day scale, after which it declines, only to rise again at scales greater than a 15-day resolution.Similarly, at the Tonzi Ranch site, the importance of Air temperature is high at the daily scale, declines until the weekly scale, and rises again slowly at coarser scales.As comparative scaling behaviors at a continuum of scales have not been studied previously, they are not documented, and these last results cannot be related to findings of earlier studies.However, given the results obtained, we conclude that the primary controls of GPP or any ecological phenomenon are strongly dependent on temporal resolution.For some covariates, the importance follows previously established general scaling behavior, whereas completely different scale-dependent patterns emerge in other cases.
Environmental forcing and biotic responses vary across temporal (and spatial) scales.The results shown here demonstrate that factors that explain GPP differ according to the scale of assessment.Hence, evaluating the importance of covariates across a continuum of temporal scales can assist in understanding the varying importance of specific environmental variables, and thereby lead to insights into processes that best explain GPP at a particular temporal scale.Furthermore, scrutiny of cross-scale controls of GPP can also help in identifying factors that are important irrespective of scale.

Conclusions
Mechanistic modeling of complex ecosystems for assessing the terrestrial carbon cycle requires understanding of numerous interlinked cause and effect relationships, which, in most open systems, can only be achieved on the foundations laid by empirical models.These empirical models, though not a substitute for detailed process-based understanding of ecosystems, go a long way in bringing to the forefront the primary factors that explain GPP and/or R E , and provide a framework for improving mechanistic models.Within this context, we propose and implement improved methods for empirical modeling of GPP observations.The statistical technique presented in this study builds on the work presented in Mueller et al. (2010), and utilizes Bayesian model selection in a geostatistical regression framework for assessing associations between GPP and environmental variables representing plant function and external forcing.The presented approach is shown to be particularly useful for discerning the scale-dependent importance of specific environmental covariates, and the uncertainty associated with the selection of covariates to be included in empirical models.As a result, the presented methods can play a valuable role in identifying candidate process-based mechanisms for explaining variability in GPP for diagnostic analysis of ecosystems.Furthermore, though not presented here, the presented approach also has the capacity to comparatively evaluate alternative representations of process-based mechanisms (e.g.different nonlinear temperature functions) in terms of their ability to explain GPP.We believe this type of future application can provide a testbed for building nonlinear or mixed linear-nonlinear empirical models of GPP and/or R E .
The results from this research confirm that Global radiation, Air temperature, and Vapor pressure deficit are the key variables that explain GPP across ecosystems, but the contribution of this work does not lie in this expected result.Instead, this work demonstrates and quantifies how the importance of environmental variables included in the empirical model of GPP varies across temporal scales.With regard to carbon cycling, the results presented in this work clearly show scale-dependence in the importance of specific covariates in explaining GPP at the examined AmeriFlux sites.Hence, having established that scale is of such primary significance, predictive or explanatory relationships of GPP based on any single scale can lead to erroneous conclusions regarding the importance of specific covariates.To avoid these pitfalls, scale-specific explanatory environmental variables should be identified and used in models formulated for a particular scale.In situations where appropriate scaling laws are unknown (which is mostly the case!), covariates with cross-scale importance identified in this research provide the best foundation for building mechanistic biospheric models.The varying importance of covariates explaining variability in GPP at different spatial and temporal scales does not imply the existence of different physiological processes at different scales, but only indicates the presence of different dominant environmental stimuli that can explain the observed pattern.However, once the scale-dependence of the importance of environmental variables is understood, the effect of these variables on plant physiology under changing environmental conditions can also be understood.Such an examination would enhance our ability to understand the response of ecosystems to climate change, and could lead to a reduction in the uncertainty surrounding future changes in the global carbon cycle.

Fig. 3 .
Fig. 3. AmeriFlux level 4 GPP and GPP estimated from the best GR model (X β) at the 8-day scale, with one standard deviation uncertainty bounds estimated as the square root of the sum of the uncertainty variance of X β (Eq.6), and the estimated variance of the residuals 2 N +σ 2 S , Eqs. 1 and 2)

Table 1 .
Description of AmeriFlux sites used in the analysis.Note that the Kennedy Space Center Scrub Oak site is referred to as "Kennedy Space Center" in subsequent tables and figures, and the Morgan Monroe State Forest site is referred to as "Morgan Monroe".

Table 3 .
Variance reduced (R 2 ) and Root mean squared error (RMSE; MgC/day/km 2 ) of the best GR model data at the daily, 8-day and monthly scales.The column "Daily (Joint Fitting)" represents the statistics based on deriving a single daily GR model and set of regression coefficients using data from all sites combined, and the associated R 2 and RMSE of this model for the individual sites.