MODIS vegetation products as proxies of photosynthetic potential along a gradient of meteorologically and biologically driven ecosystem productivity

A direct relationship between gross ecosystem productivity (GEP) estimated by the eddy covariance (EC) method and Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation indices (VIs) has been observed in many temperate and tropical ecosystems. However, in Australian evergreen forests, and particularly sclerophyll and temperate woodlands, MODIS VIs do not capture seasonality of GEP. In this study, we re-evaluate the connection between satellite and flux tower data at four contrasting Australian ecosystems, through comparisons of GEP and four measures of photosynthetic potential, derived via parameterization of the light response curve: ecosystem light use efficiency (LUE), photosynthetic capacity (Pc), GEP at saturation (GEPsat), and quantum yield (α), with MODIS vegetation satellite products, including VIs, gross primary productivity (GPPMOD), leaf area index (LAIMOD), and fraction of photosynthetic active radiation (fPARMOD). We found that satellite-derived biophysical products constitute a measurement of ecosystem structure (e.g. leaf area index – quantity of leaves) and function (e.g. leaf level photosynthetic assimilation capacity – quality of leaves), rather than GEP. Our results show that in primarily meteorological-driven (e.g. photosynthetic active radiation, air temperature, and/or precipitation) and relatively aseasonal ecosystems (e.g. evergreen wet sclerophyll forests), there were no statistically significant relationships between GEP and satellite-derived measures of greenness. In contrast, for phenology-driven ecosystems (e.g. tropical savannas), changes in the vegetation status drove GEP, and tower-based measurements of photosynthetic activity were best represented by VIs. We observed the highest correlations between MODIS products and GEP in locations where key meteorological variables and vegetation phenology were synchronous (e.g. semi-arid Acacia woodlands) and low correlation at locations where they were asynchronous (e.g. Mediterranean ecosystems). However, we found a statistical significant relationship between the seasonal measures of photosynthetic potential (Pc and LUE) and VIs, where each ecosystem aligns along a continuum; we emphasize here that knowledge of the conditions in which flux tower measurements and VIs or other remote sensing products converge greatly advances our understanding of the mechanisms driving the carbon cycle (phenology and climate drivers) and provides an ecological basis for interpretation of satellite-derived measures of greenness. Published by Copernicus Publications on behalf of the European Geosciences Union. 5588 N. Restrepo-Coupe et al.: MODIS vegetation products as proxies of photosynthetic potential

Abstract.A direct relationship between gross ecosystem productivity (GEP) estimated by the eddy covariance (EC) method and Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation indices (VIs) has been observed in many temperate and tropical ecosystems.However, in Australian evergreen forests, and particularly sclerophyll and temperate woodlands, MODIS VIs do not capture seasonality of GEP.In this study, we re-evaluate the connection between satellite and flux tower data at four contrasting Australian ecosystems, through comparisons of GEP and four measures of photosynthetic potential, derived via parameterization of the light response curve: ecosystem light use efficiency (LUE), photosynthetic capacity (Pc), GEP at saturation (GEP sat ), and quantum yield (α), with MODIS vegetation satellite products, including VIs, gross primary productivity (GPP MOD ), leaf area index (LAI MOD ), and fraction of photosynthetic active radiation (fPAR MOD ).We found that satellite-derived biophysical products constitute a measurement of ecosystem structure (e.g.leaf area index -quantity of leaves) and function (e.g.leaf level photosynthetic assimilation capacity -quality of leaves), rather than GEP.Our results show that in primarily meteorological-driven (e.g.photosynthetic active radiation, air temperature, and/or precip-itation) and relatively aseasonal ecosystems (e.g.evergreen wet sclerophyll forests), there were no statistically significant relationships between GEP and satellite-derived measures of greenness.In contrast, for phenology-driven ecosystems (e.g.tropical savannas), changes in the vegetation status drove GEP, and tower-based measurements of photosynthetic activity were best represented by VIs.We observed the highest correlations between MODIS products and GEP in locations where key meteorological variables and vegetation phenology were synchronous (e.g.semi-arid Acacia woodlands) and low correlation at locations where they were asynchronous (e.g.Mediterranean ecosystems).However, we found a statistical significant relationship between the seasonal measures of photosynthetic potential (Pc and LUE) and VIs, where each ecosystem aligns along a continuum; we emphasize here that knowledge of the conditions in which flux tower measurements and VIs or other remote sensing products converge greatly advances our understanding of the mechanisms driving the carbon cycle (phenology and climate drivers) and provides an ecological basis for interpretation of satellite-derived measures of greenness.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Introduction
Eddy flux towers constitute a powerful tool to measure and study carbon, energy, and water fluxes.Even though the number of eddy covariance (EC) sites has been steadily increasing (Baldocchi, 2014;Baldocchi et al., 2001), instrumentation, personnel costs, and equipment maintenance limit the establishment of new sites.This is demonstrated by the distribution of flux towers around the world and in particular the under-representation of tropical and semi-arid locations in the Southern Hemisphere (Australia, Africa, and South America) (http://fluxnet.ornl.gov/maps-graphicsand Beringer et al., 2007).The first EC tower was established in 1990 at Harvard Forest (Wofsy et al., 1993) followed by five other sites in 1993 (Baldocchi, 2003).In Australia, only two locations, Howard Springs (AU-How; Hutley et al., 2000) and Tumbarumba (AU-Tum; Leuning et al., 2005), have a record that extends more than 10 years.
Many applications rely on large-scale, remotely sensed (RS) representations of vegetation dynamics (greenness) to (1) upscale water and carbon fluxes from the limited tower footprint (radius < 10 km) representative of eddy covariance measurements, (2) scale fluxes in time and extend a longer time series from limited tower data, (3) fill gaps due to quality control in the flux measurements, (4) study continental phenology to be validated at flux tower sites, and (5) parameterize land surface (LSMs) and agricultural models to be tested at EC locations.Past studies have focused on the relationship between the Moderate Resolution Imaging Spectroradiometer (MODIS) VIs, such as the enhanced vegetation index (EVI), and tower-based measurements of gross ecosystem productivity (GEP) (Gamon et al., 2013;Huete et al., 2006Huete et al., , 2008;;Maeda et al., 2014;Sims et al., 2006;Wang et al., 2004).In these studies, satellite-derived vegetation indices (VIs) represented a community property of chlorophyll content, leaf area index (LAI), and fractional vegetation cover.A simple linear regression between seasonal (monthly or 16-day) EVI and GEP has previously provided a good coefficient of determination (R 2 ) for different ecosystems: where b 0 and b 1 are the fitted coefficients.Huete et al. (2006) reported an R 2 of 0.5 for Eq. ( 1) in tropical forests and converted pastures over the Amazon basin, and an R 2 of 0.74 in dry to humid tropical forest sites in Southeast Asia (Huete et al., 2008).Over the North Australian mesic and xeric tropical savannas, R 2 ranged from 0.52 at a wooded grassland (Alice Springs, AU-ASM) to 0.89 in woodlands (Howard Springs, AU-How) (Ma et al., 2013).Similar relationships to Eq. ( 1) have been explored using monthly maximal net ecosystem exchange (NEE max ): (2) This regression showed an improved fit in forests (R 2 = 0.83 for deciduous and R 2 = 0.72 for coniferous forests) com-pared to the GEP-EVI model (R 2 = 0.81 for deciduous and R 2 = 0.69 for evergreen forests) (Olofsson et al., 2008).
Other approaches to link carbon fluxes to RS products include radiation-greenness (GR) models, where both a meteorological driver, represented by the photosynthetic active radiation (PAR), and a vegetation phenology driver, represented by EVI or by the normalized difference vegetation index (NDVI), are implicitly included in the model (Ma et al., 2014;Peng and Gitelson, 2012).By definition, the GEP / PAR ratio is commonly referred to as ecosystem light use efficiency (LUE), where (3) However, the EVI vs. LUE relationship has shown lower R 2 values (0.76) compared to the EVI vs. GEP regression (0.92) for a group of North American ecosystems that included evergreen needleleaf and deciduous forests, grasslands, and savannas (Sims et al., 2006).Hill et al. (2006) also reported an R 2 of ∼ 0.2 for the NDVI vs. LUE relationship for the Australian sclerophyll forest of Tumbarumba (AU-Tum); however, this result was not statistically significant (p > 0.05).To better represent GEP in rainfall-driven semi-arid ecosystems, Sjöström et al. (2011) increased the level of complexity of the GR model by scaling down observations of PAR using the evaporative fraction (EF) term from EC measurements (a proxy for water availability); thus GEP was calculated as where EF is the ratio between latent heat flux (LE) and the surface turbulent fluxes (H + LE), and H is defined as the sensible heat flux, EF = LE/(H +LE).The model increased the predictive power of the GR model in some ecosystems; however, it was not applicable at regional scales due to its reliance upon supporting tower measurements.Temperature-greenness (GT) models use the MODIS land surface temperature product (LST) and VIs to calculate GEP as in Sims et al. (2008).The GT GEP model for nine North American temperate EC sites was calculated as where m is a function of mean annual LST and plant functional type (different formulation provided for evergreen and deciduous vegetation), LST scaled is the minimum of two equations (LST/30) and (2.5 − (0.05 × LST)), and EVI scaled is EVI − 0.10.A similar GT model, used by Wu et al. (2011), showed high correlation in deciduous forests (R 2 = ∼ 0.90) and lower R 2 values in non-forest areas (R 2 = 0.27 to 0.91) and evergreen forests (R 2 = 0.28 to 0.91).Other more complex derivations, including the C-Fix model (Veroustraete et al., 2002) and the MODIS gross primary productivity product (GPP MOD ), rely on biomespecific relationships that include (1) vegetation phenology represented by MODIS-derived fraction of absorbed PAR that a plant canopy absorbs for photosynthesis and growth (fPAR MOD ); and (2) air temperature (T air ), water vapour pressure deficit (VPD), and PAR as climate drivers (Running et al., 2000).When applied to Australian ecosystems, the GPP MOD (collection 4) was able to estimate the amplitude of the GEP annual cycle in an temperate evergreen wet sclerophyll forest (Eucalyptus-dominated); however, it was out-of-phase (Leuning et al., 2005).For a tropical savanna (AU-How), GPP MOD (collection 5) overestimated dry-season GEP (Kanniah et al., 2009).Even though GPP MOD (collection 4.8) at AU-How accurately represented seasonality in productivity, low estimates of PAR and other model input variables were compensated by abnormally high fPAR MOD values (Kanniah et al., 2009) -a clear indication of obtaining a good result for the wrong reasons.
Besides the difficulties inherent in determining GEP in diverse ecosystems, all of the complex models (e.g.GPP MOD and GT model) require in situ measurements of water fluxes, PAR, and/or biome classification information to calibrate or derive some variables, and consequently, regression coefficients do not necessarily extend to ecosystem types other than those for which the derivation was obtained.Our first objective was to revisit the GEP vs. EVI and GEP vs. GPP MOD regressions at different sites to gain a better understanding of ecosystem behaviour, rather than simply to determine the "best performing model".We look at particularly challenging land cover classes: seasonal wet-dry and xeric tropical savannas, Mediterranean environments characterized by hot and dry summers (mallee), and temperate evergreen sclerophyll forests.The selected locations are part of the OzFlux eddy covariance network and represent sites where previous studies have shown satellite-derived GEP models to be unable to replicate in situ measurements.
Our second objective was to derive, using the light response curve, different ground-based measures of vegetation photosynthetic potential: quantum yield (α), photosynthetic capacity (Pc), GEP at saturation light (GEP sat ), and ecosystem light use efficiency (LUE), in an attempt to separate the vegetation structure and function (phenology) from the climatic drivers of productivity.We explored the seasonality of the four measures of photosynthetic potential (α , Pc, LUE, GEP sat ), and aimed to determine whether EVI was able to replicate absolute value and their annual cycle rather than photosynthetic activity (GEP), based on linear regressions.Similarly, we included other MODIS biophysical datasets (NDVI, LAI MOD , and fPAR MOD ) in our analysis in an effort to understand how to interpret different satellite measures of greenness and how these products can inform modellers and ecologists about vegetation phenology.In contrast to biomespecific classification approaches, we treated the relationship between greenness and photosynthetic potential to be a continuum, and therefore, we explored multiple site regressions.
Our third objective was to combine satellite-derived meteorology (radiation, precipitation, and temperature) and biological drivers (vegetation phenology) to determine site-  Kottek et al. (2006) and Rubel and Kottek (2010), where Aw is equatorial winter dry climate, BSh is arid steppe, BWh is hot arid desert, BWk is cold arid desert, Cfb is warm temperate fully humid warm summer, Cfa is warm temperate fully humid hot summer, and Cwa is warm temperate winter dry hot summer.Other climate classes are equatorial fully humid (Af) and monsoonal climate (Am), arid summer dry and cold desert (Bsk), and warm temperate hot summer (Csa) and warm summer (Csb) steppes.specific and multi-biome GEP values using multiple regression models.In this study, we evaluated the advantages of introducing both types of variables; we investigated whether the regressions hold across biomes, and whether productivity processes are driven by phenology, light, water availability, and/or temperature; and we infer which of these variables govern the GEP seasonal cycle for each particular ecosystem.These results advance our understanding of driving mechanisms of the carbon cycle (climate, biological adaptation, or a combination of both) and temporal and spatial scaling, and they provide an ecological basis for the interpretation of satellite-derived measures of greenness and phenology products.

Study sites
The OzFlux infrastructure network is operated by a collaborative research group and was set up to provide the Australian and global ecosystem modelling communities with CO 2 and H 2 O flux and meteorological data (Beringer et al., 2016).We selected four contrasting long-term eddy flux (EC) sites from the OzFlux network (Fig. 1  Fluxes at all towers were measured by the EC method with an open-path system.Simultaneously, an array of different sensors measured meteorological data including air temperature (T air ), relative humidity (RH), incoming and reflected shortwave radiation (SW down and SW up ), and incoming and reflected longwave radiation (LW down and LW up ).Refer to each site reference for complete information regarding ecosystem and measurement techniques.

Eddy covariance data
We used Level 3 OzFlux data that include an initial OzFlux standard quality control (QC) (Isaac et al., 2016).All data were subject to the same quality assurance (QA) procedures and calculations, providing methodological consistency among sites and reducing the uncertainty of the calculated fluxes.We performed additional quality checks and removal of outliers, and data were corrected for low turbulence periods (see Sect. 2.2.1).Ecosystem respiration (R eco ) and GEP were calculated from EC measurements of net ecosystem exchange (NEE) as presented in Sect.2.2.2.Finally, we derived different measures of ecosystem vegetation photosynthetic potential (Sect.2.2.3).

Eddy covariance and meteorological measurements
Incoming and outgoing radiation, both shortwave (SW down and SW up ) and longwave (LW down and LW up ), were measured using a CNR1 net radiometer instrument (Campbell Scientific).All sensors were placed above the canopy at the same height or higher than the EC system.As there were no measurements of PAR radiation available at AU-ASM, AU-Tum, and AU-Cpr, we assumed PAR = 2 × SW (Papaioannou et al., 1993;Szeicz, 1974), where PAR is measured as the flux of photons (µmol m −2 s −1 ) and SW down as the heat flux density (W m −2 ).We understood this as an approximation because PAR radiation (0.4-0.7 nm) is a spectral subset of SW down (0.3-3 nm).
At AU-Tum, the NEE is calculated as the sum of the turbulent flux measured by eddy covariance (F C ) plus changes in the amount of CO 2 in the canopy air space (storage flux, S CO 2 ), where NEE = F C + S CO 2 .At all other sites, given the sparse vegetation cover and the smaller control volume over the vegetation that is lower in height (< 15 m), F C is assumed to be representative of NEE.
Hourly fluxes measured during rainy periods, when the sonic anemometer and the open-path infrared gas analyser (IRGA) do not function correctly, were identified and removed from the time series.We also removed isolated observations (between missing values).We identified any residual spikes from the hourly NEE data using the method proposed by Papale et al. (2006) and modified by Barr et al. (2009).
For each hour (i), the measure of change in NEE (d i ) from the previous (i − 1) and next (i + 1) time step is calculated as A spike is identified if the change is outside a given range: where M d is the median of the differences (d i ), ±0.6745 are the quartiles for a standard normal distribution, and the constant z was conservatively set to 5 (Restrepo-Coupe et al., 2013).Night-time hourly NEE values were corrected for periods of low turbulent mixing by removing them from the time series data.Low turbulent missing periods were determined when friction velocity (u * in m s −1 ) was below a threshold value (u * thresh ) as described in Restrepo-Coupe et al. (2013).Table 1 presents site-specific u * thresh values and the corresponding upper and lower confidence bounds.Night-time NEE was assumed to be representative of ecosystem respiration (R eco ), and it was calculated by fitting R eco to a second-order Fourier regression based on the day of the year (DOY) as in Richardson and Hollinger (2005): where fo, e, s 1 , c 1 , s 2 , and c 2 are the fitted coefficients and Dpi = DOY × 360/365 in radians.This method calculates R eco with minimal use of environmental covariates.In order to determine the consistency of the Fourier regression method and the low friction velocity (u * ) filter on the modelled R eco (directly dependent on night-time NEE values), we compared the results presented here to R eco values based on the intercept of the relation (rectangular hyperbola) between NEE and SW down (for no incoming radiation, SW down = 0) (Suyker and Verma, 2001) (Supplement Fig. S1).
Gross ecosystem exchange (GEE) was calculated as the difference between NEE and R eco (GEE = NEE+R eco ).We defined gross ecosystem productivity (GEP) as negative GEE (positive values of GEP flux indicate carbon uptake).For a 16-day moving window, we fitted two rectangular hyperbolas on the relationship between incoming PAR and GEP observations (separating morning and afternoon values) as in Johnson and Goody (2011) and based on the Michaelis and Menten formulation (1913): where α is the ecosystem apparent quantum yield for CO 2 uptake (the initial slope), and GEP sat is GEP at saturating light (the asymptote of the regression) (Falge et al., 2001) (Fig. 2).Our intention was to compare 16-day MODIS data to observations rather than to model a complete time series.We therefore filled infrequent GEP missing values only if there were 30 h of measurements in a 16-day period.
We obtained similar seasonal patterns and good agreement using different methods for calculating GEP and R eco (Fig. S1).We observed no statistically significant seasonal differences between calculating R eco as the intercept of the light response curve (Falge et al., 2001), and NEE not subject to u * thresh correction (R eco LRC ), and calculating R eco using the Fourier regression method (slope ∼ 0.87 and R 2 = 0.94 linear regression between R eco LRC and R eco ).This comparison increased our confidence in using either method to derive From the rectangular hyperbola: quantum yield (α, µmolCO 2 µmol −1 ) (blue dashed line) and GEP at saturation (GEP sat , µmolCO 2 m −2 s −1 ) (blue dotted line).Photosynthetic capacity (Pc, µmolCO 2 m −2 s −1 ) (black dashed line) was calculated as the 16-day mean GEP at mean annual daytime PAR (PAR) ±100 µmol m −2 s −1 (grey area) and mean annual VPD (VPD) ±2 standard deviations.Light use efficiency (LUE, µmolCO 2 µmol −1 ) was defined as the ratio between daily GEP and PAR, the slope of the linear regression (blue line).
the GEP and R eco fluxes from the EC data, their absolute values, and the seasonality presented here.
GEP and GPP (true photosynthesis minus photorespiration; Wohlfahrt and Gu, 2015) have been used interchangeably in the literature.However, GPP in this study was distinguished from GEP; thus as GEP does not include CO 2 recycling at leaf level (i.e.re-assimilation of dark respiration) or below the plane of the EC system (i.e.within canopy volume) (Stoy et al., 2006), differences may be important when comparing flux-tower observations of GEP to the MODIS GPP product (see next section).

Four measures of ecosystem photosynthetic potential: α, LUE, GEP sat , and Pc
Measures of photosynthetic potential constitute an attempt to separate the inherent vegetation properties that contribute to photosynthetic activity (GEP) from the effects of the meteorological influences on productivity using the parameterization of the 16-day light response equation.The variables α, LUE, GEP sat , and Pc were intended to represent an ecosystem property, a descriptor of the vegetation phenology similar to leaf area index (LAI) or above ground biomass (AGB).We calculated 16-day mean α and GEP sat , which are the two coefficients that define the GEP vs. PAR rectangular hyperbola (Eq.5) as a measure of the vegetation structure and function (Fig. 2).Both α (µmol CO 2 mmol −1 ) and GEP sat (µmol CO 2 m −2 s −1 ) values are known to vary with vegetation type, temperature, water availability, and CO 2 concentration.The GEP sat represents the ecosystem response at saturating levels of PAR, usually constrained by high vapour pressure deficit (VPD), air temperature (T air ), water availability, and foliar N, among other variables (Collatz et al., 1991;Ehleringer et al., 1997;Tezara et al., 1999).By contrast, α is measured at low light levels, when diffuse radiation is high (cloudy periods, sunset, and sunrise).Ecosystem light use efficiency (LUE) was defined as the mean daily GEP / PAR ratio.Therefore, LUE includes the effect of day length, the radiation environment (diffuse vs. direct), water availability, and other physical factors.We used the relationships between tower-measured GEP, PAR, and VPD to characterize the photosynthetic capacity of the ecosystem (Pc), where Pc was defined as the average GEP for incoming radiation at light levels that are non-saturating -values between the annual daytime mean PAR ± 100 µmol m −2 s −1 (940, 1045, 788, and 843 µmol m −2 s −1 at AU-How, AU-ASM, AU-Tum, and AU-Cpr, respectively), and VPD ranges between annual daytime mean ± 2 standard deviations (Fig. 2) (Hutyra et al., 2007;Restrepo-Coupe et al., 2013).Pc was interpreted as a measure of the built capacity without taking into account the dayto-day changes in available light, photoperiod, and extreme VPD and PAR values.The derivation of Pc did not take into account other variables such as T air or soil water content.

Moderate Resolution Imaging Spectroradiometer (MODIS)
We retrieved MODIS reflectances, VIs, and other products from the USGS repository covering the four eddy flux locations.Data were subject to quality assurance (QA) filtering, and pixels sampled during cloudy conditions and pixels adjacent to cloudy pixels were rejected (for a complete list of QA rules see Supplement Table S1).Other QA datasets and/or fields related to the above products that were not included in the original metadata were not examined as part of the quality filtering process.
At each site we extracted either a 1 km window (or a 1.25 km window depending on MODIS product resolution -see Table 2) centred on the location of the flux tower.The mean and standard deviation of all pixels were assumed to be representative of the ecosystem.The derivative data collection included the following MODIS data (also see Table 2).

MCD43A1
The 8-day 500 m (collection 5) nadir bidirectional reflectance distribution function (BRDF) adjusted reflectance (NBAR) product was used to derive the enhanced vegetation index (EVI SZA30 ) and the normalized vegetation index (NDVI SZA30 ) at fixed solar zenith an-gle of 30 • (available for 2003 to 2013): where R SZA30 , NIR SZA30 , and B SZA30 are the red, nearinfrared, and blue band BRDF-corrected reflectances, and coefficients G = 2.5, C1 = 6, C2 = 7.5, and L = 1 (Huete et al., 1994).Both VIs are measures of greenness and have been designed to monitor vegetation, in particular photosynthetic potential and phenology (Huete et al., 1994;Running et al., 1994).However, the EVI has been optimized to minimize the effects of soil background, and to reduce the impact of residual atmospheric effects.
We labelled the NBAR VIs as EVI SZA30 and NDVI SZA30 to differentiate them from the MOD13 VI product (EVI and NDVI), and emphasize that the values presented here include a BRDF correction that aims to remove the influence of sun-sensor geometry on the reflectance signal (Schaaf et al., 2002).

MOD15A2
The leaf area index (LAI MOD ), and fraction of photosynthetically active radiation (fPAR MOD ) absorbed by vegetation from atmospherically corrected surface reflectance products (Knyazikhin et al., 1999) were recorded.Data were filtered to remove outliers present in the fPAR MOD and LAI MOD time series using Eq.(3).A threshold value of 6 for the z coefficient was calibrated to remove 8-day variations of ±50 % on fPAR MOD , and ±3-4 units in LAI MOD .
MOD17A2 The 8-day gross primary production (GPP MOD ) and net photosynthesis (PsnNet) (collection 5.1) are calculated.The GPP MOD is calculated using the formulation proposed by Running et al. (2000) and relies on satellite-derived shortwave downward solar radiation (SW down ), fPAR MOD , maximum light-use-efficiency (ε max ) obtained from a biome-properties look-up table, and maximum daily VPD (VPD max ) and minimum daily air temperature (T min ) from forcing meteorology: where only the highest quality data were selected for the analysis.
MOD11A2 Daytime land surface temperature (LST day ) 8day time series was included in the analysis in order to study the effect of T air , another important ecosystem carbon flux driver.Thus, LST or skin temperature (temperature at the interface between the surface and the atmosphere) has been proven to be highly correlated to T air (Shen and Leptoukh, 2011).(Huffman et al., 2007).We used 1.0 • resolution monthly surface shortwave flux down (all sky) in watts per square metre from the Clouds and the Earth's Radiant Energy System (CERES) experiment (Gesch et al., 1999).The CERES Energy Balanced And Filled top of the atmosphere (EBAF) Surface_Ed2.8product provided fluxes at surface, consistent with top of the atmosphere fluxes (CERES EBAF-TOA) (Kato et al., 2012).No quality control was performed on the rain (Precip TRMM ) or shortwave (SW CERES ) satellite-derived time series.We used satellite-derived meteorological variables instead of in situ measurements as the independent variable in GEP models (see Sect. 2.5); thus, our findings (e.g.regressions) can be extrapolated to regional and continental scales.

Mean values
All analyses were done on 16-day data; therefore, 8-day MODIS products were resampled to the match the selected temporal resolution.We interpolated lower frequency satellite remote sensing time series (e.g.CERES and TRMM), using a linear regression from the original dataset to 16-days, where the original value corresponds to the centre of the month defined as day 15, and the newly interpolated value will be representative of the middle of the 16-day period.
Mean fluxes and variables from the eddy covariance are reported on a 30 min or hourly basis.Daily averages were calculated if at least 45 out of 48, or 21 out of 24 data points were available for the day.Bi-weekly values were calculated if at least 4 out of the 16 days were available.For analysis and presentation purposes, we averaged all existing 16day values of EC and RS data to produce a single-year, seasonal cycle.We understand measures of photosynthetic potential to be dependent on the selection of the aggregation period.However, the 16-day interval has been shown to be representative of important ecological processes, in particular, leaf appearance to full expansion (Jurik, 1986;Varone and Gratani, 2009), greenup of soil biological crusts in response to precipitation events (Cleverly et al., 2016a), and reported ecosystem-level changes in ecosystem water use efficiency (Shi et al., 2014).

Evaluation of synchronicity between remote sensing and flux tower data
We fitted type II (orthogonal) linear regressions that account for uncertainty in both variables (satellite and EC).
We obtained an array of very simple models of productivity and photosynthetic potential.For example, GEP RS , where GEP RS = b 0 + b 1 × RS, b 0 and b 1 were site-specific coefficients, and RS are satellite-derived products (EVI, fPAR, etc.).We compared the different models to the observations (GEP vs. GEP EVI , GEP vs. GEP NDVI , etc.) using Taylor single diagrams (Taylor, 2001), where the radial distances from the origin are the normalized standard deviation, and the azimuthal position is the correlation coefficient between the GEP RS and GEP or any other measure of ecosystem photosynthetic potential (Fig. S2).We determined at each site which combination of carbon flux and MODIS index showed good agreement based on statistical descriptors: coefficient of determination, p value, root-mean-square error (RMSE), standard deviation (SD) of the observation and model, and the Akaike's information criterion (AIC).Thus, we analysed site-specific and cross-site multiple regression models to compare different biological (greenness) and environmental controls (precipitation, temperature, and radiation) on productivity.In each ecosystem, GEP was modelled as a linear regression using a single independent variable, two variables, and bivariate models that included an interaction term.For example, (1  compared to the site-specific results.We inferred ecosystem adaptation responses to climate (e.g.light harvest adaptation, water limitation, among other phenological responses) from the bivariate models.This analysis is useful for the interpretation of satellite-derived phenology metrics and understanding the biophysical significance of different measures of greenness when incorporated into land surface models as representative of vegetation status (Case et al., 2014).

Seasonality of in situ measurements
In this section we describe the seasonality of in situ meteorological measurements to better understand ecosystem carbon fluxes, and to contextualize the differences in vegetation responses to climate.In particular, we contrast seasonal patterns of air temperature (T air ), precipitation, and VPD across sites, and compare observations of the annual cycle of photosynthetic activity (productivity) and potential (biophysical drivers of productivity) for each ecosystem.
With the exception of AU-How, all sites showed strong seasonality in T air (Fig. 3).However, the timing of mean daily T air minimum and maximum, and the amplitude of the annual values, varied according to site.The smallest range in T air (5 • C) occurred at the northern tropical savanna (AU-How), and the largest amplitude (15 • C) occurred at the southern temperate locations (AU-Cpr and AU-Tum).The annual cycle of VPD followed T air at all locations except AU-How where summer and autumn rains (February-March) led to a increase in atmospheric water content (Fig. 3).Precipitation at AU-How was higher and more seasonal than at any other site, with a mean monthly rainfall of 152 mm (1824 mm yr −1 ) and ranging from 1 to 396 mm month −1 .Incoming radiation at the tropical savanna site (AU-How) did not show clear seasonality (Fig. 3).In this tropical savanna (latitude 12.49 • S) the summer solstice, where top of the at-mosphere (TOA) radiation is highest, coincides with monsoonal cloudiness, resulting in reduced surface radiation.By contrast, at temperate sites like AU-Cpr and AU-Tum, the difference in mean daily PAR between summer and winter was ∼ 460 µmol m −2 s −1 .Rainfall was aseasonal at AU-Tum (∼ 78 mm month −1 ) and was very low at the semi-arid sites of AU-Cpr and AU-ASM, with mean precipitation values of 34 and 37 mm month −1 respectively.
Productivity in the four ecosystems ranged from a high at AU-How and AU-Tum (Fig. 4) (peak 16-day multi-year average GEP of 8.4 and 7.7 gC m −2 d −1 respectively) to a low at AU-Cpr and AU-ASM (peak 16-day annual average GEP average of 2.4 and 3.4 gC m −2 d −1 respectively) (Fig. 4).There was a clear seasonal cycle in photosynthetic activity with maxima in the summer at AU-How and AU-Tum (November-March) and in the autumn (March-April) at AU-ASM and AU-Cpr.The peaks were broader at AU-Tum than at AU-How and at AU-ASM (Fig. 4).An additional shortlived increase in GEP was apparent at AU-ASM in the spring (October) before the summer wet period (Fig. 4a).Figures S3  and S4 show the diel cycles of VPD, GEP, and other meteorological and flux variables in example summer (January) and winter months (July).
Vegetation phenology, as indicated by the seasonal cycle of photosynthetic potential (Pc, LUE, α, and GEP sat ), diverged from photosynthetic activity (GEP) at the southern locations of AU-Tum and AU-Cpr as shown by the differences in the timing of maximum and minimum GEP compared to vegetation phenology (Figs. 4 and S5).At the tropical savanna site (AU-How), ecosystem quantum yield (α) increased gradually in the spring (September), reaching a maximum during the summer month of January in synchrony with GEP.In the sclerophyll forest (AU-Tum), α remained at a constant value of ∼ 1.4 gC MJ −1 until the middle of the autumn (April-May) when it reached a value of 1.76 gC MJ −1 .Maximum GEP sat occurred during the summer at this site (∼ 36 gC m −2 d −1 ), and gradually decreased by the start of the autumn with a winter minimum (20 gC m −2 d −1 ).At AU- Tum, the GEP sat and α were out-of-phase (Fig. 4) and although seasonality was limited in GEP sat and α, neither of them matched seasonal fluctuations in VPD (cf.Figs. 3 and  4).Similar to GEP sat , LUE decreased during the summer months and experienced a winter maximum opposite to the annual cycle of GEP.Given the high degree of seasonality of GEP at AU-Tum, it is interesting that the photosynthetic potential was comparatively less seasonal and asynchronous to productivity.Figure S5 shows the relationships between the different measures of ecosystem performance, indicating that they are not always linear.

Seasonality of satellite products
In the tropical savanna (AU-How) the annual cycles of RS products synchronously reached an early summer maximum in January, and high values extended throughout the autumn (Fig. 4d and e).By contrast at AU-Cpr, both NDVI SZA30 and EVI SZA30 peaked in autumn-winter, coinciding with the lowest GEP values (Fig. 4p and s).EVI SZA30 and NDVI SZA30 at AU-ASM captured the autumn peak in GEP with a maximum in March; however, a spring VI minimum (November) was not observable in GEP.At the two semi-arid sites (AU-ASM and AU-Cpr), fPAR MOD was relatively aseasonal, and the amplitude of the annual cycle was ∼ 0.09, with a 0.25-0.34range at AU-Cpr and lower values between 0.17 and 0.26 at AU-ASM (Fig. 4o).LAI MOD at AU-Cpr reached a maximum of 0.50 during the autumn (March) and a spring minimum (September) of 0.39.At AU-ASM, the LAI MOD product ranged from 0.17 (December) to 0.27 (April) (Fig. 4t).Most RS products (e.g.EVI SZA30 and LAI MOD ) showed no clear seasonality at AU-Tum (Fig. 5i and j).

Relationship between MODIS EVI and GPP and in situ measures of ecosystem photosynthetic activity (GEP)
In this study we used a simple linear model to predict GEP from EVI SZA30 and GPP MOD .We observed three patterns.First, in the tropical savanna site (AU-How) there was a highly significant correlation between photosynthetic activity and EVI SZA30 , where EVI SZA30 explained 82 % of GEP www.biogeosciences.net/13/5587/2016/Biogeosciences, 13, 5587-5608, 2016 (Fig. 5a).Similarly at AU-ASM, productivity was statistically related to EVI SZA30 (R 2 = 0.86, p < 0.01).However, GPP MOD only explained 49 % of GEP at AU-How and 48 % at AU-ASM (Fig. 5e and g).
A second pattern was observed in the sclerophyll forest site (AU-Tum), where the relationship between GEP and EVI SZA30 was not statistically significant (R 2 < 0.01 and p = 0.93, Fig. 5b).At AU-Tum there was a clear seasonal cycle in GEP (low in winter and high during the summer) that was not captured by the small amplitude of the satellitederived data (Fig. 3).Of the four ecosystems examined, AU-Tum was the only site where GPP MOD showed an improvement (higher predictive value of GEP) compared to EVI SZA30 .However, as reported in previous works (Leuning et al., 2005), the GPP MOD product was unable to capture the seasonality of the sclerophyll forest as it underestimated the observed summer peak in GEP which corresponded to a second minimum in GPP MOD .
Finally, at the semi-arid site (AU-Cpr), we observed R 2 values significantly different from 0 but a small R 2 value of 0.34 and 0.24 (p < 0.01) for GEP vs. EVI SZA30 and GEP vs. GPP MOD , respectively.This demonstrated the low predictive power of both satellite products to determine seasonal GEP values in this Mediterranean ecosystem.In particular the GEP EVI and GPP MOD models tended to underestimate productivity at low levels (Fig. 5d and h).
The relationship between productivity and EVI SZA30 was complex across the different Australian ecosystems (Fig. 5).
The semi-arid site of AU-Cpr and the sclerophyll forest of AU-Tum are particularly interesting because of the inability of EVI SZA30 to seasonally replicate GEP (Fig. 5).An additional analysis that considers the amplitude and phase of the annual cycle (based on all available 16-day observations) was conducted using Taylor plots (Fig. S7).This analysis showed that EVI SZA30 was in-phase and able to predict the range of productivity values at AU-How and AU-ASM, while at the AU-Cpr site the EVI SZA30 captured the amplitude of seasonal GEP; however, the linear model was out-ofphase.At AU-Tum, the EVI SZA30 -based model consistently preceded in situ observations (asynchronous) and exaggerated GEP seasonality (ratio between the standard deviation of the model and observations was 4.98).

Relationship between EVI SZA30 and measures of photosynthetic potential (α, LUE, GEP sat , and Pc)
In this section we reconsider our understanding of EVI SZA30 by relating it to different measures of photosynthetic potential (α, LUE, GEP sat , and Pc) across the four sites (Fig. 6).Similar to Sect.3.3, we used a very simple linear model in which EVI SZA30 was expected to predict α, LUE, GEP sat , and Pc.In the regression models for photosynthetic potential the R 2 values were similar to the GEP models for AU-How and AU-ASM (cf.Fig. 6c and g).However, EVI SZA30 vs. α at AU-How R 2 was relatively low (R 2 < 0.4, p < 0.01).At the AU-Cpr site, the EVI SZA30 -based model was able to improve the timing and amplitude of the annual cycle when used to calculate LUE, Pc, and GEP sat instead of GEP (Fig. 6 and  S7).
At the sclerophyll forest site (AU-Tum) the EVI SZA30 was able to predict vegetation phenology rather than productivity.For example we observed that Pc (but not α) was significantly related to EVI SZA30 (R 2 = 0.16, p < 0.01; Fig. 6 and Table S4).Even though the regressions between LUE, GEP sat , and Pc against EVI SZA30 showed higher correlation (R 2 ∼ 0.13, p < 0.01) than the GEP vs. EVI SZA30 relationship (R 2 = 0.04, p = 0.25) at AU-Tum, R 2 values were still low.The low R 2 can be explained by the small dynamic range of both seasonal measures of photosynthetic potential and EVI SZA30 (cf.Figs. 4 and 6).

Satellite products compared to flux-tower-based measures of ecosystem potential
In this section we explore other MODIS products (LAI MOD , fPAR MOD , and NDVI SZA30 ) to determine whether the predictive power of EVI SZA30 as a measure of photosynthetic potential (e.g.Pc) can be generalized across other satellitederived biophysical parameters.We aimed to determine, for each location, which of the MODIS products capture the seasonality and phenology of vegetation, thereby gaining some insight into the significance of the different VIs and other satellite-derived ecosystem drivers.At AU-How and AU-ASM the MODIS LAI MOD , fPAR MOD , and VIs showed a larger or similar correlations to LUE and Pc in comparison to GEP (Table S4, Fig. 7a and b and i and j, respectively).At AU-How, AU-ASM, and AU-Cpr, based on our analysis using Taylor plots, most RS products were in-phase with the various measures' phenology (R 2 > 0.8 and low RMSE) (Figs. 7, S2 and Table 4).However, there was a tendency for most RS indices to underestimate the seasonality of the LUE annual cycle at all sites (i.e., standard deviation was smaller for LUE RS than the observed, Fig. 7).With the exception of AU-Tum, all products were able to capture seasonal changes in Pc (Figs. 6 and 7).Similar to EVI SZA30 , most of the MODIS indices, and in particular fPAR MOD and LAI MOD , showed strong linear relationships with LUE and Pc at the Mediterranean ecosystem AU-Cpr, where the introduction of phenology represented an important improvement over the RS-derived models (Figs. 6  and 7).Similarly, comparable to EVI SZA30 , other MODIS www.biogeosciences.net/13/5587/2016/Biogeosciences, 13, 5587-5608, 2016 products were unable to replicate GEP at AU-Tum (Fig. 7).However, the small amplitude of seasonality in LUE and Pc were well characterized by LUE RS and Pc RS , including a winter maximum similar to that in LUE (Fig. 4), despite underestimating the annual seasonal cycle in the sclerophyll forest (Figs. 4 and 7e-h).

Multi-biome-derived linear relationships between
VIs and photosynthetic potential (phenology) and activity (productivity) Our objective was to investigate whether one relation fits all flux sites, and which RS products and equations would enable us to extend our analysis from these four key Australian ecosystems to a continental scale.The all-site relationship for MODIS EVI SZA30 , NDVI SZA30 , LAI MOD , and fPAR MOD products (in that order) shows the best agreement (phase and amplitude) to seasonality of LUE and Pc (Fig. 7).Correlations increased for relationships built using data for all the ecosystems instead of the site-specific equations, with the exception of the AU-ASM site (Table 3 and Figs.7 and 8).
Improvements in how satellite products can model biological drivers (photosynthetic potential) instead of productivity per se, are clearly seen at the evergreen temperate forest of AU-Tum.At AU-Tum the relationship between GEP and any of the satellite products was not statistically significant (R 2 < 0.1), with the exception of LST day (Figs. 5 and  7).However, skin temperature (LST day ) is a meteorological driver rather than a direct measure of productivity, and the low all-site LST day vs.GEP correlation was an indication of this (R 2 = 0.66, p = 0.03; Fig. 8).The wet sclerophyll forest introduced the greatest uncertainties to the linear models across all sites (Fig. 8).For example, regressions involving EVI SZA30 were exponential; therefore, significantly increasing GEP and LUE translated into slightly higher EVI SZA30 values, a behaviour mostly driven by the observations at AU-Tum.In particular, the relationships between LUE and fPAR MOD and LUE and NDVI SZA30 at AU-Tum were problematic as fPAR MOD and NDVI SZA30 appeared to "saturate" at 0.9 and 0.8, respectively (Fig. 8).
EVI SZA30 explained 81 % of Pc seasonality based on an all-site regression (Table S4).Similarly, NDVI SZA30 showed a high coefficient of determination (0.70 for GEP NDVI , 0.75 for LUE NDVI , and 0.79 for Pc NDVI ) (Table S4).The null hypothesis of no correlation was rejected (p < 0.01) for all regressions between MODIS VIs, LAI MOD , and fPAR MOD vs. photosynthetic potential (phenology) and activity (productivity) (Table S4).However, statistical significance of GEP vs. GEP RS was driven by the AU-ASM and AU-How ecosystems.
Multiple linear regression models used to predict GEP by combining satellite-derived meteorology and biologic parameters (Table 3) showed large correlations when both drivers were introduced (weather and vegetation phenology), with the exception of the AU-Tum site where SW CERES and LST day explained 60 and 58 % of GEP, respectively, and the AU-ASM and AU-How sites where EVI SZA30 and NDVI SZA30 explained ∼ 84 and ∼ 80 % of the variations in GEP, respectively.In particular, at the AU-How site, no significant improvement to the GEP model was obtained when combining MODIS VIs with any meteorological variable (R 2 remains similarly high: R 2 ∼ 0.82).By contrast, at the AU-ASM site, EVI SZA30 , satellite-derived incoming shortwave (SW CERES ), and the interaction of both significantly increased model correlation with an R 2 of 0.88 and a lower AIC (Akaike's information criterion as a measure of model quality) when compared to models relying only on EVI SZA30 (R 2 = 0.85, AIC = 64) or SW CERES (R 2 = 0.02, AIC = 209) (Table 3).Similar results were obtained for those regressions driven by EVI SZA30 and precipitation at this rainfall-pulse-driven site (R 2 = 0.88, AIC = 42).At the AU-Cpr site, temperature-greenness models were highly correlated to GEP (R 2 > 0.64); however, the best results (higher R 2 and lower AIC) were obtained for radiationgreenness models, explaining 71 % (EVI SZA30 − SW CERES and NDVI SZA30 − SW CERES ) of GEP.For a complete version of Table 3 that includes all available variable combinations, see Table S3.

Discussion
4.1 Derivation of measures of photosynthetic potential at tropical savannas, sclerophyll forests, and semi-arid ecosystems In this study we were able to separate the biological (vegetation phenological signal) from the climatic drivers of productivity using eddy covariance carbon exchange data.Using the parameterization of the light response curve, we derived different measures of vegetation photosynthetic potential (α, LUE, GEP sat , and Pc) (Balzarolo et al., 2015;Wohlfahrt et al., 2010).At seasonal timescales (e.g.16-days, monthly), our analysis looks at the biotic drivers of productivity; whereas at shorter timescales (e.g.hourly, daily), photosynthetic potential can be limited or enhanced by meteorological controls, thus linked to resource scarcity (i.e.high VPD or water constraints), or availability (e.g. increase radiation or access to soil water), and correspondent ecosystem responses (e.g.stomatal closure, CO 2 fertilization) will determine GEP (Ainsworth and Long, 2005;Doughty et al., 2014;Fatichi et al., 2014).The variables α, LUE, GEP sat , and Pc have different biophysical meanings; therefore, we were able to establish physiological explanations for describing why and which RS products and environmental variables relate to them in each ecosystem.For example, GEP sat measured at high levels of PAR is prone to be influenced by various environmental factors (VPD, T air , and soil water availability), and therefore may be a good indicator of canopy stress.As observed at AU-How, GEP sat was highly and negatively correlated to periods of low precipitation and negatively correlated with VPD (Table S4).Seasonal values of GEP sat at the semi-arid sites (AU-Cpr and AU-ASM) did not show a direct relationship with VPD or precipitation.This does not mean that there is no effect of atmospheric demand or soil moisture content on carbon fluxes at shorter timescales (hourly or daily) (Cleverly et al., 2016b;Eamus et al., 2013).Compared to GEP sat , we expected α to be less dependent on VPD and to better reflect vegetation phenology, as α represents the canopy photosynthetic response at low levels of PAR characteristic of cloud cover (diffuse light) during early morning or late afternoon periods (Kanniah et al., 2012(Kanniah et al., , 2013)).However, among all measures of phenology, α showed one of the lowest site-specific correlations when compared to any of the RS products presented in this study.Our results show that LUE and Pc showed the best correlations to VIs; this is confirmation that this research deals less with the instantaneous responses (GEP sat and α) and rather www.biogeosciences.net/13/5587/2016/Biogeosciences, 13, 5587-5608, 2016   focuses on the mid-term, 16-day seasonal descriptors of vegetation phenology (Pc and LUE).
The influence of other environmental factors apart from PAR and VPD, such as soil water content and T air , is difficult to isolate from the derivation of vegetation descriptors as there may be a high degree of cross-correlation between the different variables (e.g.VPD vs. T air ).Moreover, to what degree it is feasible to untangle the relations between climate and vegetation is complex and not well understood, as the feedback processes are essential in ecosystem function (leaf flush, wood allocation, among other vegetation strategies respond to available resources), species competition, and herbivory cycles (Delpierre et al., 2015).Our results show that VIs were highly related to Pc, which is interpreted as a phenology descriptor that does not consider the day-to-day changes in available light or photoperiod or the vegetation response to high and low VPD and PAR values.By contrast, implicit in the derivation of LUE were the day length and anomalous climatic conditions.This finding has important implications when using EC data for the validation of satellite-derived phenology (Restrepo Coupe et al., 2015).

Seasonality and comparisons between satellite
products and flux-tower-based measurements of carbon flux: photosynthetic activity (productivity) and potential (phenology) Previous satellite-derived models of productivity usually apply to locations where the seasonality of GEP is synchronous with climatic and vegetation phenology drivers (Mahadevan et al., 2008;Sims et al., 2008;Wu et al., 2010;Xiao et al., 2004), such as in temperate deciduous forests, where temperature and incoming radiation coincide with changes in ecosystem structure and function (e.g.autumn sub-zero temperatures may initiate leaf abscission; Vitasse et al., 2014).
In our analysis, productivity was only synchronous with all measures of photosynthetic potential at the savanna site (AU-How), where clouds and heavy rainfall in the summer wet season resulted in low VPD, reduced TOA (aseasonal PAR), and minimal fluctuations in T air .At AU-How, we observed a consistently large correlation between MODIS VIs and productivity and no improvement in GEP when accounting for meteorology.Moreover, the highly significant EVI SZA30 vs. GEP relationship at AU-How could be generalized to other satellite-derived biophysical products.Arid and semi-arid vegetation dominate ∼ 75 % of the Australian continent, and in these ecosystems a characteristic mix of grasses (understorey) and woody plants (overstorey) contributes to total annual GEP at different times of the year.More importantly, the phenology of grasses and trees is driven by, or responds differently to, various climatic drivers (e.g.trees greening up after spring rainfalls while grasses remain dormant (Cleverly et al., 2016a;Ma et al., 2013;Shi et al., 2014).The changing seasonal contributions to the reflectance signal and to GEP are generally related to soil water content thresholds.Our study presents two semiarid Acacia and Eucalyptus woodlands, where we found that models relating VIs with photosynthetic potential (phenology), rather than activity (productivity), improved the predictive power of RS greenness indices (AU-Cpr) or showed sim-ilar statistical descriptors (AU-ASM).At the woodland Acacia site, LAI MOD and fPAR MOD overestimated the periods of low capacity (associated with browndown phases) (Ma et al., 2013).This can be better understood if we account for small but non-negligible photosynthetic activity in Acacia after the summer rains have ended (Cleverly et al., 2013;Eamus et al., 2013).At this particular site (AU-ASM), the high LAI MOD and VIs observed during dormancy may not be interpreted as high photosynthetic potential.Satellite data, and even some ground-based measurements of LAI MOD , cannot differentiate between the different fractional components: photosynthetic active vegetation (fPAV), and non-photosynthetic vegetation (fNPV).Future work requires phenocams or biomass studies in which fPV and fNPV may be spectrally or mechanically separated.
In low productivity ecosystems (AU-ASM and AU-Cpr), satellite and EC data/noise ratio may have a considerable effect on the site-specific regressions (e.g.sun geometry influence on VIs seasonal values, and EC uncertainties).However, differences between AU-ASM and AU-Cpr regressions (e.g.EVI SZA30 was highly correlated to GEP only at AU-ASM), and the fact that the VI product has been corrected for BRDF effects, increases our confidence in the analysis presented here.Moreover, the lower VIs vs. GEP correlation values obtained at AU-Cpr compared to AU-ASM could be attributed to mallee site productivity being more dependent on meteorological drivers than photosynthetic potential, or GEP being driven by climate (e.g.autumn precipitationwhen Pc remains constant) or by vegetation phenology (e.g.summer LAI and canopy chlorophyll content, among others) at different times of the year.
Similar to Mediterranean ecosystems (AU-Cpr), in wet sclerophyll forests (AU-Tum) without signs of water limitation, the VIs were unable to replicate seasonality in GEP.In particular, the dominant species of sclerophyll forests, Eucalyptus, Acacia, and Banksia, show very little seasonal variation in canopy structure as seen in aseasonal LAI observations (Zolfaghar, 2013), and leaf longevity (Eamus et al., 2006).Leaf quantity (e.g.LAI) and quality (e.g.leaf level photosynthetic assimilation capacity) are two key parameters in driving photosynthetic potential; when these are aseasonal, asynchronous, or lagged, they may confound the interpretation of seasonal measures of greening.Thus, the observed increasing predictive power of VIs as a measure of photosynthetic potential (e.g.EVI SZA30 vs. Pc, R 2 = 0.16 at AU-Tum) may not be comparable to similar relationships at sites where vegetation phenology showed a larger dynamic range (e.g.EVI SZA30 vs. Pc, R 2 = 0.79 at AU-How).

Considerations for the selection of RS data to be used on GEP models and phenology validation studies
This study reports high correlations for Pc vs. EVI SZA30 (R 2 = 0.81) and Pc vs. NDVI SZA30 (R 2 = 0.80).The fact that a brighter soil background results in lower NDVI values than with a dark soil background for the same quantity of partial vegetation cover (Huete, 1988;Huete and Tucker, 1991) may have a positive effect on the all-site Pc vs. NDVI SZA30 regressions (increase R 2 ).However, darkened soils following precipitation also raise NDVI values for incomplete canopies (Gao et al., 2000), and may similarly suggest higher vegetation or soil biological crust activity.On the other hand, soil brightness and moisture may have a negative effect on the confidence interval of the x intercept for the proposed relationships (e.g.Pc vs. NDVI SZA30 , for NDVI SZA30 ∼ 0).Moreover, at certain times the AU-ASM and AU-Cpr sites were at the low end of the vegetation activity range, and the observed RS signal may have been dominated by soil water content rather than by photosynthetic potential.Furthermore, caution is needed when using fPAR MOD and other products as we observed a threshold value above which in situ changes were undetectable (e.g.MODIS fPAR > 0.9, NDVI SZA30 > 0.8).This might have been due to the NDVI saturating at high biomass (Huete et al., 2002;Santin-Janin et al., 2009).Temperature-greenness models of GEP (Sims et al., 2008;Xiao et al., 2004) take into account the meteorological and biophysical drivers that determine productivity.Nevertheless, correlations between photosynthetic characteristics and LST day were weaker than for VIs.Moreover, if the seasonality of GEP is driven by local climatology, as in the case of AU-Tum where GEP was statistically correlated to LST day , our intent is to understand the relation between vegetation characteristics and RS products rather than indiscriminately use any satellite-derived index to describe phenology or photosynthetic potential.Our study demonstrates that multiple linear regression models that combine satellite-derived meteorology and biological parameters to describe GEP fit better when both drivers are introduced rather than when only one factor drives the relation (a single meteorology or greenness variable).However, two exceptions to this rule were observed: (1) at AU-Tum where SW CERES was able to explain 60 % of GEP, and (2) in the tropical savanna at AU-How where EVI SZA30 was able to explain ∼ 82 % of the variation in GEP, and where we did not obtain any significant improvement to the GEP model when combining MODIS VIs and any meteorological variable (R 2 remains similarly high: R 2 > 0.82).In summary, in evergreen sclerophyll forests, even when GEP is highly seasonal, GEP is driven by meteorology as seen by the fact that most of the measures of photosynthetic potential showed small seasonal changes, similar to different MODIS products.By contrast, at sites where most of the GEP seasonality was driven by vegetation status (Pc as a proxy) rather than the meteorological inputs (PAR, air temperature, and precipitation), or where meteorology and phenology were synchronous, VIs were strongly correlated to both GEP and Pc (e.g.tropical savanna).This was in agreement with the expectation than RS products constitute a measurement of ecosystem photosynthetic potential rather than productivity per se.
Our analysis shows how MODIS greenness indices were able to estimate different measures of ecosystem photosynthetic potential across biomes.At only one site (AU-Tum) was there very little seasonal variation in EVI SZA30 , compared to other evergreen ecosystems.Both the strong correlations among VIs and Pc from in situ EC carbon flux measurements at the remaining sites (AU-How, AU-ASM, and AU-Cpr), and the positioning of each ecosystem along a continuum of MODIS-derived variables representing vegetation phenology, confirm the usefulness of satellite products as being representative of vegetation structure and function.This research confirms the viability of remote-sensingderived phenology to be validated and, more importantly, understood, using eddy-flux measurements of Pc.However, an increase in effort in determining seasonal patterns of carbon allocation (partition between leaves and wood), understorey and overstorey responses, and leaf carbon assimilation and chlorophyll content over time may be required to obtain a more meaningful understanding of RS indices and their biophysical significance.Moreover, the reader should be aware that rapid changes in vegetation phenology (e.g.α and GEP sat ) caused by short-term environmental stresses (e.g.T air , humidity, soil water deficit, or waterlogging) may not be accurately estimated by RS products, and may require the employment of in situ high-frequency optical measurements (e.g.phenocams), land surface vegetation models, or direct EC measurements.
For this study we included all available 16-day data corresponding individually to more than 10 years at AU-How and AU-Tum, and 2 to 3 years at AU-Cpr and AU-ASM.The long-term sampling implies that we were likely to be capturing a large range in mean ecosystem behaviour.RS products may over-or under-represent the canopy response to periods of extreme temperature and precipitation, although the time series in this study included years that were warmer than normal and heatwaves, e.g.2012-2013(BOM, 2012, 2013;;van Gorsel et al., 2016) and years that were wetter than normal, e.g.2011 (Fasullo et al., 2013;Poulter et al., 2014), which led to GEP that is larger than normal at AU-ASM and AU-Cpr (Cleverly et al., 2013;Eamus et al., 2013;Koerber et al., 2016).It is beyond the scope of this work to evaluate the inter-annual variability of the vegetation responses to disturbance (e.g.insect infestation or fire) or extreme climatic events (e.g.flooding or long-term drought).Improvements to satellite-derived phenology can be related to an increasing number of EC sites and samples, thereby emphasizing the importance of long-term time measurements and sampling of diverse ecosystems.

Conclusions
Satellite vegetation products have been widely used to scale carbon fluxes from eddy covariance (EC) towers to regions and continents.However, in some key Australian ecosystems, MODIS gross primary productivity (GPP) product and vegetation indices (VIs) do not track seasonality of gross ecosystem productivity (GEP).In particular, we found that EVI SZA30 was unable to represent GEP at the temperate evergreen sclerophyll forest of Tumbarumba (AU-Tum) and at the Mediterranean ecosystem (mallee) of Calperum-Chowilla (AU-Cpr).This result extends across satellite products overall: MODIS GPP MOD , LAI MOD , fPAR MOD , and other VIs.
We aimed for a greater understanding of the mechanistic controls on seasonal GEP, and proposed the parameterization of the light response curve from EC fluxes as a novel tool to obtain ground-based seasonal estimates of ecosystem photosynthetic potential (light use efficiency (LUE), photosynthetic capacity (Pc), GEP at saturation (GEP sat ), and quantum yield (α)). Photosynthetic potential refers to the presence of photosynthetic infrastructure in the form of ecosystem structure (e.g.leaf area index -quantity of leaves) and function (e.g.leaf level photosynthetic assimilation capacity -quality of leaves) independent of the meteorological and environmental conditions that drive GEP.Based on basic linear regressions, we demonstrated that MODIS-derived biophysical products (e.g.VIs) were a proxy for ecosystem photosynthetic potential rather than GEP.We reported statistically significant regressions between VIs (e.g.NDVI SZA30 and EVI SZA30 ) and long-term measures of phenology (e.g.LUE and Pc), in contrast to ecosystem descriptors subject to short-term responses to environmental conditions (e.g.GEP sat and α).Our results should extend to other methods and measures of greenness, including VIs and chromatic indices from phenocams and in situ spectrometers.
We found that the linear regressions between MODIS biophysical products and photosynthetic potential converged on a single function across very diverse biome types, which implies that these relationships may persist over very large areas, thus improving our ability to extrapolate in situ phenology and seasonality to continental scales, across longer temporal scales, and enabling us to identify rapid changes due to extreme events or spatial variations at ecotones.We further found that saturation of fPAR MOD and NDVI SZA30 restricted their usefulness, except in comparatively low biomass ecosystems (savannas and arid and semi-arid savannas and woodlands).
We quantified how much GEP seasonality could be explained by different variables: radiation (SW down ), temperature (T air ), precipitation (Precip), or phenology (VIs as proxy).Our analysis showed that the relationship between RS products and GEP was only clear when productivity was driven by either (1) ecosystem phenology and climate, synchronously driving GEP, as was observed at Alice Springs www.biogeosciences.net/13/5587/2016/Biogeosciences, 13, 5587-5608, 2016 mulga woodland (AU-ASM), and similar to many temperate deciduous locations, or (2) solely the vegetation photosynthetic potential, as observed at the tropical savanna site of Howard Springs (AU-How).At AU-How, radiation and temperature were constant throughout the year, although ecosystem photosynthetic activity (GEP) and potential (e.g.Pc and LUE) fluctuated with the highly seasonal understorey.However, RS products do not follow GEP when (3) phenology is asynchronous with key meteorological drivers such that GEP is driven by one or the other at different times of the year, as we observed at AU-Cpr; or when (4) GEP is driven by meteorology (SW down , T air , soil water availability, VPD, or different combinations), and photosynthetic potential is aseasonal, as observed at AU-Tum.At AU-Tum, changes in productivity were driven by SW down , while the ecosystem biophysical properties remained relatively constant throughout the year, represented by the small amplitude of the annual cycles in Pc and LUE (true evergreen forest).An understanding of why satellite vs. flux tower estimates of GEP relationships hold, or do not hold, greatly contributes to our comprehension of carbon cycle mechanisms and scaling factors that are at play (e.g.climate and phenology, among others).
The Supplement related to this article is available online at doi:10.5194/bg-13-5587-2016-supplement.

Figure 1 .
Figure 1.Location of four OzFlux eddy covariance tower sites included in this analysis: AU-How: Howard Springs (at Aw), AU-ASM: Alice Springs mulga (at BSh and BWh boundary), AU-Cpr: Calperum-Chowilla (at Bwk), and AU-Tum: Tumbarumba (at Cfa and Cfb boundary).Köppen-Geiger climate classification as published by Kottek et al. (2006) andRubel and Kottek (2010), where Aw is equatorial winter dry climate, BSh is arid steppe, BWh is hot arid desert, BWk is cold arid desert, Cfb is warm temperate fully humid warm summer, Cfa is warm temperate fully humid hot summer, and Cwa is warm temperate winter dry hot summer.Other climate classes are equatorial fully humid (Af) and monsoonal climate (Am), arid summer dry and cold desert (Bsk), and warm temperate hot summer (Csa) and warm summer (Csb) steppes.

Figure 4 .
Figure 4. Savanna (AU-How), wet sclerophyll (AU-Tum), mulga (AU-ASM), and mallee (AU-Cpr) ecosystems, OzFlux sites annual cycle (16-day composites) of eddy-flux-derived (a) gross ecosystem productivity (GEP; gC m −2 d −1 ) (black line) and MODIS gross primary productivity (GPP MOD ) product (light blue line); (b) GEP at saturation light (GEP sat ; gC m −2 d −1 ) (black line) and ecosystem quantum yield (α; gC MJ −1 ) (light blue line); (c) photosynthetic capacity (Pc; gC m −2 d −1 ) (black line) and the ratio of GEP over PAR (black line), the light use efficiency (LUE; gC MJ −1 ) (light blue line).At the bottom two panels, satellite-derived data of (d) MODIS enhanced vegetation index at fixed solar zenith angle of 30 • (EVI SZA30 ) (black line) and the normalized difference vegetation index (NDVI SZA30 ) (light blue line); (e) MODIS leaf area index (LAI MOD ) (black line) and MODIS fraction of the absorbed photosynthetic active radiation (fPAR MOD ) (light blue line).Grey boxes indicate Southern Hemisphere spring and summer September to March.The black dashed vertical line indicates the timing of maximum GEP.

Figure 7 .
Figure 7. Taylor diagrams showing model results for Howard Springs (AU-How), Tumbarumba (AU-Tum), Alice Springs (AU-ASM), and Calperum-Chowilla (AU-Cpr) based on site-specific and all-sites linear regressions between gross ecosystem productivity (GEP), light use efficiency (LUE), photosynthetic capacity (Pc), and ecosystem quantum yield (α) and different remote sensing products MODIS fixed solar zenith angle of 30 • enhanced vegetation index (EVI) and normalized difference vegetation index (NDVI), gross primary productivity product (GPP), daytime land surface temperature (LST), leaf area index (LAI), and fraction of the absorbed photosynthetic active radiation (fPAR).All-site relationships are labelled with an asterisk (e.g.EVI * ).EVI and NDVI labels are used instead of EVI SZA30 and NDVI SZA30 for displaying purposes.Missing sites indicate that the model overestimates the seasonality of observations -the model normalized standard deviation is > 2.
, • C), fixed solar zenith angle of 30 • normalized difference vegetation index (NDVI SZA30 ), precipitation from the Tropical Rainfall Measuring Mission (Precip TRMM , mm month −1 ) data product from 1998 to 2013 (NASA, 2014b), and surface shortwave incident radiation from the Clouds and the Earth's Radiant Energy System (SW CERES , W m −2 ) data product from 2000 to 2013 (NASA, 2014a).Model runs for AU-How: Howard Springs, AU-ASM: Alice Springs mulga, AU-Cpr: Calperum-Chowilla, and AU-Tum: Tumbarumba, and all available data (data include all sites).Bold fonts highlight values mentioned in the text.

Table 1 .
and Table 1) for this study.OzFlux sites presented in this study -location and additional information.

Table 2 .
Remote sensing data sources, cell size, sample size (eddy covariance tower site at the centre pixel), and time interval.

Table 3 .
Linear regressions obtained by a non-linear mixed-effects regression model for gross ecosystem productivity (GEP, gC m −2 d −1 ) vs. combinations of 16-day average MODIS products: fixed solar zenith angle of 30 • enhanced vegetation index (EVI SZA30 ), daytime and land surface temperature(LST day