Monitoring Grassland Seasonal Carbon Dynamics, by Integrating MODIS NDVI, Proximal Optical Sampling, and Eddy Covariance Measurements

This study evaluated the seasonal productivity of a prairie grassland (Mattheis Ranch, in Alberta, Canada) using a combination of remote sensing, eddy covariance, and field sampling collected in 2012–2013. A primary objective was to evaluate different ways of parameterizing the light-use efficiency (LUE) model for assessing net ecosystem fluxes at two sites with contrasting productivity. Three variations on the NDVI (Normalized Difference Vegetation Index), differing by formula and footprint, were derived: (1) a narrow-band NDVI (NDVI680,800, derived from mobile field spectrometer readings); (2) a broad-band proxy NDVI (derived from an automated optical phenology station consisting of broad-band radiometers); and (3) a satellite NDVI (derived from MODIS AQUA and TERRA sensors). Harvested biomass, net CO2 flux, and NDVI values were compared to provide a basis for assessing seasonal ecosystem productivity and gap filling of tower flux data. All three NDVIs provided good estimates of dry green biomass and were able to clearly show seasonal changes in vegetation growth and senescence, confirming their utility as metrics of productivity. When relating fluxes and optical measurements, temporal aggregation periods were considered to determine the impact of aggregation on model accuracy. NDVI values from the different methods were also calibrated against fAPARgreen (the fraction of photosynthetically active radiation absorbed by green vegetation) values to parameterize the APARgreen (absorbed PAR) term of the LUE (light use efficiency) model for comparison with measured fluxes. While efficiency was assumed to be constant in the model, this analysis revealed hysteresis in the seasonal relationships between fluxes and optical measurements, suggesting a slight change in efficiency between the first and second half of the growing season. Consequently, the best results were obtained by splitting the data into two stages, a greening phase and a senescence phase, and applying separate fits to these two periods. By incorporating the dynamic irradiance regime, the model based on APARgreen rather than NDVI best captured the high variability of the fluxes and provided a more realistic depiction of missing fluxes. The strong correlations between these optical measurements and independently measured fluxes demonstrate the utility of integrating optical with flux measurements for gap filling, and provide a foundation for using remote sensing to extrapolate from the flux tower to larger regions (upscaling) for regional analysis of net carbon uptake by grassland ecosystems.


Introduction
The need to better understand and predict future carbon-climate interactions makes the assessment of biosphere-atmosphere carbon exchange for the various biomes of the planet a critical topic for the scientific community and policy-makers [1].There is growing interest in assessing biospheric carbon storage with the hope that management regimes can optimize carbon sequestration.Also, soil processes play an important role in global climate change as soils have the potential to act as a net sink for CO 2 due to the large amount of carbon currently stored in soil organic matter [2,3].However, biosphere-atmosphere carbon fluxes are inherently dynamic, affected by short-term disturbance and weather events, making evaluation and prediction of carbon storage and climate feedbacks challenging.While much attention has been given to forests for their large feedback and carbon sequestration [4], the carbon storage potential of other biomes has been less thoroughly considered, even though many non-forest biomes cover large areas, contribute measurably to overall biospheric carbon uptake and storage, and provide economically important ecosystem goods and services.
Rangelands occupy close to 50% of the world's land area [5].They are defined as lands where indigenous or introduced vegetation is grazed or has the potential to be grazed and are generally considered and managed as natural ecosystems [6].Rangelands include natural grasslands, savannas, shrublands, deserts, tundra, alpine communities, marshes, and wet meadows [7].The grassland or prairie biome is particularly important in Alberta, as the Grassland Natural Region covers 14.4% of the province [8].The Dry Mixedgrass Natural Subregion accounts for 47.5% of the Grassland Natural Region area, and thus represents 7% of the area of Alberta province [9].Grassland provides a number of important services including carbon sequestration and livestock forage.The capacity of grassland ecosystems to sequester carbon is a leading factor in the consideration of policies and practices that maximize carbon storage without compromising managers' profits [10].The carbon sequestration of a grassland can vary considerably from year to year since it is influenced both by natural and anthropogenic factors such as temperature, rainfall, species composition, nutrient and water availability, light, grazing pressure and agricultural practices [11,12].Together, these factors make grassland carbon exchange responsive to climate change [13][14][15], and drive interest in monitoring grassland ecosystem responses for the purpose of developing optimal management regimes.
Historically, most attempts to assess ecosystem carbon storage have been based on field sampling of above-and below-ground biomass and carbon, which requires considerable investments in time and personnel [16], and may not capture the temporal dynamics in response to changing conditions.Ecosystem carbon dynamics have been modeled using different input parameters including climate [17,18] or soil and nutrient availability [19,20].However, models often lack spatially or temporally sufficient data inputs, so may not take into account all ecosystem dynamics and interactions.
Significant progress in understanding ecosystem dynamics has been accomplished using the eddy covariance (EC) technique, which quantifies carbon and water vapor fluxes between the biosphere and the atmosphere [21].This technique has been successfully applied to estimate carbon fluxes in grassland ecosystems [22,23].Nevertheless, the EC technique is expensive, subject to data gaps, and limited to relatively flat, large landscape regions, so cannot offer universal sampling of fluxes for all locations.
Remote sensing and proximal optical sampling methods can provide cost-effective, uniform sampling under conditions where eddy covariance techniques would be impractical.Satellite remote sensing combined with weather data have been used to model global carbon exchange for many years now, but these current models often disagree with local field measurements due to the coarse scale of their optical and meteorological inputs [24].For effective carbon policy, fine-grained measurements are needed to match the scale of local land management practices.Recent advances in automation and low-cost optical sensors are making a local-scale optical sensing approach increasingly attractive.Recent studies show that the combined application of remote and proximal optical sensors provide the opportunity to monitor productivity across a wide range of spatial and temporal scales and offer the possibility to relate remote sensing measurements with carbon fluxes as a foundation for modeling and upscaling [25].Networks such as SpecNet and EUROSPEC have been exploring the possibility of combining flux and spectral reflectance measurements for a more complete assessment of ecosystem function [26][27][28].These collaborations have increased the availability of both optical and flux data, allowing for improved models linking remote sensing measurements and carbon fluxes.A crucial first step is to quantitatively evaluate optical methods in comparison to flux methods to enable proper model parameterization.
The most widely applied model for combining carbon fluxes and optical measurements is the Light Use Efficiency (LUE) model.Originating with the work of Monteith [29,30], more recent parameterizations have been used to estimate gross primary productivity (GPP) based on the equation: GPP " ε ˆAPAR green (1) where APAR green (absorbed photosynthetically active radiation) represents the amount of radiation absorbed by green vegetation, and ε represents the efficiency with which that absorbed light is used by vegetation to fix carbon.To obtain annual GPP, the equation terms can be integrated (summed) over the full growing season [24,31].APAR green can be further defined or measured as the product of f APAR green (the fraction of photosynthetically active radiation absorbed by green vegetation) and PAR (photosynthetically active radiation, also called the photosynthetic photon flux density, or PPFD): Damm et al.
[32] reported that, especially in cropland and grassland, canopy structure and APAR covary in time during a seasonal cycle.Partitioning between photosynthetic and non-photosynthetic components of vegetation in APAR calculations significantly improves estimation of ecosystem productivity with respect to models driven by total f APAR [33,34].APAR green , a measure of green vegetation structure closely linked to green biomass or leaf area index, plays a dominant role in the model for most grassland ecosystems.The confusion of methods and LUE model definitions in the literature has led to uncertainty when comparing model results [34].To clarify the terminology and to emphasize the role of green (not total) vegetation in driving light absorption for photosynthesis and productivity, we adopted the terminology recommended by Gitelson and Gamon [34] and included the subscript "green" when referring to fAPAR and APAR.
Light use efficiency (ε) relates to vegetation physiology, and is presumably less important in determining seasonal GPP patterns for grasslands and other annual systems where canopy growth dominates the seasonal phenological signals detectable with remote sensing [35].Consequently, most of the variability in grassland productivity can be explained by changes in APAR green [36] when the ecosystem is not stressed [33,35,37].For this reason, we adopted the hypothesis of a constant ε over the season (May-October) in our study.
The strength of the LUE model is related to the possibility of deriving its input values from remote sensing observations [24,38,39].Several vegetation indices (VIs) have been used to estimate f APAR green , most notably the normalized difference vegetation index (NDVI) [40].When combined with PAR irradiance, either as a modeled or directly measured value (using PAR sensors), f APAR green readily provides a reliable measure of APAR green , the primary term in the LUE model (Equation ( 2)).
The rising interest in linking carbon flux measurements, proximal (in situ) spectral measurements and remote sensing data has motivated the establishment of a global network of experimental sites in various ecosystems that simultaneously collect optical and flux datasets [26,27].A key issue in combining optical and flux measurements is the evaluation of temporal scale, as optical and flux measurements vary on different time scales [26].To our knowledge, few studies explicitly consider the effect of temporal aggregation when comparing remote sensing and flux data, in part because many satellite sensors (e.g., the MODIS Aqua and Terra sensors) are in sun-synchronous orbit that limit the number of daily overpasses.In this study, we took advantage of the continuous nature of flux and optical data to enable an explicit consideration of temporal aggregation, as recommended by the SpecNet community [27].
The objectives of this study were to: ‚ assess the seasonal patterns of CO

Study Site and Experimental Design
The Mattheis Ranch is a 4977 ha working cattle ranch located near Duchess in southern Alberta, Canada, and is mostly characterized by native prairie vegetation.Recently, this site has become a research area managed by the Rangeland Research Institute (University of Alberta), offering the possibility of long-term monitoring of ecosystem dynamics under management regimes typical of this region.The climate of southern Alberta presents long winters and short, windy summers with extreme temperatures in both seasons, and peak precipitation periods occurring in early summer (June).The rest of the summer usually has lower precipitation, but summer rainfall regimes can vary markedly from year to year, with a similarly variable effect on productivity [12].The management regime consisted of light, rotational grazing, as further described in Wang et al. [41].
The two sites differed slightly in vegetation cover and microtopography; E5 appeared very flat and uniform, while the E3 landscape varied slightly more than E5 in microtopography and vegetation cover within the flux tower footprint [41], but was still a relatively uniform grassland-dominated site.The sites represent a typical dry mixed grass prairie community as described by Adams et al. [9].The E5 soil is a clay loam, and the plant community is dominated by needle-and-thread grass (Hesperostipa comata), Junegrass (Koeleria macrantha), blue grama grass (Bouteloua gracilis), and western wheatgrass (Pascopyrum smithii).Pasture sage (Artemisia frigida) and scarlet mallow (Sphaeralcea coccinea) are the most common forbs, and silver sagebrush (Artemisia cana) is the primary shrub species [42].In contrast, the soil at the E3 site is a sandy loam.Dominant grasses include sand grass (Calamovilfa longifolia), needle-and-thread grass (Hesperostipa comata), and low sedge (Carex stenophylla), while common forbs are scurf pea (Psoralidium lanceolatum), golden aster (Heterotheca villosa), and prickly-pear cactus (Opuntia polyacantha).Slight variations in topography allow forbs such as wild licorice (Glycyrrhiza lepidota) and golden bean (Thermopsis rhombifolia), and shrubs such as wild rose (Rosa woodsii) to flourish [42].Due to these differences in vegetation, soil type and associated hydrology (not shown), these two sites provided a natural experiment encompassing two contrasting landscapes typical of southern Albertan rangelands.The measurements at the two grassland sites (E3 and E5) included: -CO2 flux measurements from eddy covariance at the E3 and E5 sites, -continuous, proxy NDVI [27,43] from a set of 2-band radiometers, including PAR (photosynthetically active radiation) and PYR (pyranometer) sensors, comprising an "optical phenology station", -aboveground biomass samples from a 1 ha area around each flux tower site periodically collected following the sampling scheme reported in Figure 2, -NDVI680,800 measurements using a field spectrometer collected at regular intervals in the same area (1 ha) using a 10 m grid spacing (Figure 2), -incident and reflected incoming PAR measurements at each calibration point (numbered circles, Figure 2) -MODIS satellite NDVI data downloaded for the study areas.
Each phenology station was located approximately 10 meters to the southeast of its corresponding flux tower.These two sites (E3 and E5) were monitored for two consecutive growing seasons (May-September) during 2012 and 2013.The measurements at the two grassland sites (E3 and E5) included: -CO 2 flux measurements from eddy covariance at the E3 and E5 sites, -continuous, proxy NDVI [27,43] from a set of 2-band radiometers, including PAR (photosynthetically active radiation) and PYR (pyranometer) sensors, comprising an "optical phenology station," -aboveground biomass samples from a 1 ha area around each flux tower site periodically collected following the sampling scheme reported in Figure 2, -NDVI 680,800 measurements using a field spectrometer collected at regular intervals in the same area (1 ha) using a 10 m grid spacing (Figure 2), -incident and reflected incoming PAR measurements at each calibration point (numbered circles, Figure 2), and -MODIS satellite NDVI data downloaded for the study areas.
Each phenology station was located approximately 10 m to the southeast of its corresponding flux tower.These two sites (E3 and E5) were monitored for two consecutive growing seasons (May-September) during 2012 and 2013.

Biomass Estimation
From May to August 2012, vegetation samples were harvested from four biomass calibration points approximately every 20 days using a 30 cm diameter ring in both sites (E3 and E5) at a distance of 70 m from the tower in four directions (South-East, North-East, North-West, South-West) (Figure 2).The grass material within the ring was cut at ground level and placed into labeled paper bags.The following year, the sample size was increased, and samples from 12 points were collected from May to July about every 20 days (Figure 2).Note that each biomass sample was collected from a slightly different site within 2 m of each of the 12 locations indicated above (Figure 2) to avoid repeatedly harvesting the same locations.To develop calibrations between biomass, fAPARgreen, and NDVI680,800, biomass samples were collected at each of the harvest sampling points just after the collection of fAPAR and reflectance measurements used to calculated NDVI680,800.Each sample was manually sorted into green biomass and brown biomass, put into an oven at 60 °C for 24 h, and weighed.The selection was carried out considering the living tissue that was visibly green as green biomass and the visibly dead tissue as brown biomass.This selection of green and dead biomass allowed us to measure the current year's production, since the green tissue represented the current year's growth, and the dead tissue consisted of the previous year's growth.For each date, average green and total (green plus brown) above-ground biomass was calculated from the manually sorted and weighed biomass samples.For each sampling date, green biomass was calculated as the average of green biomass at all sampling locations (Figure 2), expressed as g•m −2 .By providing a direct metric of productivity, this measure of biomass provided an independent check on the validity of the optical (NDVI) measurements and the resulting LUE model.

fAPARgreen and APARgreen Calibration
At each of the 12 calibration points (Section 2.2 above; numbered circles, Figure 2), we also measured incident, transmitted and reflected PAR (photosynthetically active radiation) using a light bar (AccuPAR LP-80, Decagon, Pullman, WA, USA) as well as NDVI680,800 from spectrometer measurements (Figure 2).At each site, light bar measurements consisted of downwelling PAR above the canopy (S), upwelling PAR reflected from the canopy (R), downwelling PAR below the canopy

Biomass Estimation
From May to August 2012, vegetation samples were harvested from four biomass calibration points approximately every 20 days using a 30 cm diameter ring in both sites (E3 and E5) at a distance of 70 m from the tower in four directions (South-East, North-East, North-West, South-West) (Figure 2).The grass material within the ring was cut at ground level and placed into labeled paper bags.The following year, the sample size was increased, and samples from 12 points were collected from May to July about every 20 days (Figure 2).Note that each biomass sample was collected from a slightly different site within 2 m of each of the 12 locations indicated above (Figure 2) to avoid repeatedly harvesting the same locations.To develop calibrations between biomass, f APAR green , and NDVI 680,800 , biomass samples were collected at each of the harvest sampling points just after the collection of f APAR and reflectance measurements used to calculated NDVI 680,800 .Each sample was manually sorted into green biomass and brown biomass, put into an oven at 60 ˝C for 24 h, and weighed.The selection was carried out considering the living tissue that was visibly green as green biomass and the visibly dead tissue as brown biomass.This selection of green and dead biomass allowed us to measure the current year's production, since the green tissue represented the current year's growth, and the dead tissue consisted of the previous year's growth.For each date, average green and total (green plus brown) above-ground biomass was calculated from the manually sorted and weighed biomass samples.For each sampling date, green biomass was calculated as the average of green biomass at all sampling locations (Figure 2), expressed as g¨m ´2.By providing a direct metric of productivity, this measure of biomass provided an independent check on the validity of the optical (NDVI) measurements and the resulting LUE model.

fAPAR green and APAR green Calibration
At each of the 12 calibration points (Section 2.2 above; numbered circles, Figure 2), we also measured incident, transmitted and reflected PAR (photosynthetically active radiation) using a light bar (AccuPAR LP-80, Decagon, Pullman, WA, USA) as well as NDVI 680,800 from spectrometer measurements (Figure 2).At each site, light bar measurements consisted of downwelling PAR above the canopy (S), upwelling PAR reflected from the canopy (R), downwelling PAR below the canopy (T), and upwelling PAR below the canopy (U).Ten measurements were made at each calibration point and averaged, then f APAR was calculated, following the procedure described in Nassar et al. [44] and the equation: fAPAR where t is the fraction of radiation transmitted through the canopy (t = T/S), r is the fraction of radiation reflected by the canopy (r = R/S), and r s is the reflectance of the soil surface (r s = U/T).These f APAR measurements were multiplied by the green:total biomass fraction to calculate a f APAR green value (fraction of PAR absorbed by green canopy material) for comparison with the proxy NDVI values (phenology stations) and NDVI 680,800 values collected at each harvest site (Figure 2).To calculate the APAR green term of the LUE model (Equation ( 2)), we derived the following NDVI-f APAR green relationships from the f APAR calibration sites: fAPAR green " p0.8799 ˆNDVI 680,800 q ´0.1394 pR 2 " 0.79q (4) fAPAR green " p1.3916 ˆproxy NDVIq ´0.3603 pR 2 " 0.82q fAPAR green " p0.9879 ˆMODIS NDVIq ´0.2131 pR 2 " 0.79q In these equations, proxy NDVI and MODIS NDVI were simulated from narrowband reflectance (as described in Section 2.4 below).In the final model, we used Equation (5) (with proxy NDVI) since it yielded the best calibration (highest R 2 ).In accordance with Equation (2), continuous APAR green , was calculated as the product of PAR (from phenology station) and f APAR green values.

NDVI 680,800
Narrow-band reflectance measurements were obtained with a dual channel spectrometer (UniSpec-DC, PP-Systems, Amesbury, MA, USA), which has a spectral range of 305-1130 nm and a «3 nm nominal bandwidth (10 nm full width at half maximum).The spectrometer was fitted with two optical fibers, one looking upward and one looking downward, enabling simultaneous sampling of downwelling irradiance and target radiance and correction for variable sky conditions [45].The upward-looking detector was fitted with a hemispherical cosine head (UNI435, PP Systems, Amesbury MA, USA) and sampled downwelling radiation; the downward-looking detector was fitted with a fiber optic (UNI684, PP Systems, Amesbury MA, USA) and field-of-view restrictor (hypotube, UNI688, PP Systems, Amesbury MA, USA) and sampled upwelling radiation with a nominal field-of-view of approx.20 degrees from a distance of approx. 2 m.
The spectrometer sampling procedure involved walking a series of transects across a 1-hectare grid centered on the flux tower, with each sample approximately 10 m apart (Figure 2).This grid sampling was completed approximately every 20 days during the growing season within one hour of solar noon.This allowed us to calculate the NDVI 680,800 of an area approximating the flux footprint and to compare the field spectral measurements with those derived from the MODIS satellite sensors.A separate study [41] demonstrated that, due to the uniformity of the landscape, the NDVI values of these different sampling regions were very similar, with an error of 4% or less.In addition to the grid sampling, the field spectrometer was also used to measure spectral reflectance at exactly the same location of each aboveground biomass sampling site for the purpose of NDVI-biomass calibration (numbered circles, Figure 2).Spectral processing software (Multispec, http://specnet.info)was used to calculate raw reflectance values, calculated as the ratio between upwelling surface radiance (R target ) and downwelling solar irradiance (I downwelling ): Raw reflectance " R target {I downwelling (7) The next step was to calculate a cross-calibration value comparing downwelling solar irradiance (I downwelling ) to a 99% reflective white standard panel (R panel ) (Spectralon, Labsphere Inc., North Sutton, NH, USA): Cross ´calibration " I downwelling {R panel (8) This cross-calibration was then used to calculate a corrected reflectance value as follows: R corrected " pR target {I downwelling q ˆpI downwelling {R panel q (9) This application of the cross-calibration was corrected for optical differences among sensors and for changing irradiance [45].Both target and panel readings were taken under similar illumination and sun angle conditions.White panel measurements were taken at the beginning and ending of every grid sampling procedure.
For NDVI 680,800 calculations, reflectance was interpolated to 1-nm wavebands.NDVI 680,800 was obtained from spectrometer grid data using the 680 nm waveband for the red region of the spectra and the 800 nm waveband for the infrared region, and was calculated using the following formula: NDVI 680,800 " pρ 800 ´ρ680 q{pρ 800 `ρ680 q (10) where ρ represents the reflectance at a given waveband, and subscripts indicate the wavelength values used.NDVI 680,800 was compared to green biomass and f APAR green (see Sections 2.2 and 2.3), but because of the low sampling frequency (about every 20 days) was not compared to daily flux values or used directly in the LUE model.For f APAR green calibration, reflectance spectra from the calibration sites were convolved against PAR, PYR (pyranometer), and MODIS NDVI band responses by multiplying each reflectance spectrum by the sensor band response for each wavelength, and dividing the summed product by the sum of the band responses across all wavelengths.These bands were then used to simulate proxy and MODIS NDVI values for calibration against f APAR green (see Section 2.3 above).

Proxy NDVI Measurements
Each site was monitored with an optical phenology station consisting of a data logger (H21-001, Onset Computer Corporation, Bourne, MA, USA) and two-band radiometer mounted on a boom and tripod, 3 m above the ground, placed approximately 10 meters southwest of each tower.The sensors detected incoming light (irradiance) and reflected light from the canopy.One band consisted of two PAR sensors (S-LIA, Onset Computer Corporation) that measured the photosynthetic photon flux density (PPFD) within the PAR band (400-700 nm), and the other band consisted of two PYR sensors (S-LIB, Onset Computer Corporation) measuring irradiance or radiance across a spectral range from 300 to 1100 nm.Both up-and down-looking sensors had cosine (nominally 180 degree) foreoptics, allowing us to sample a relatively large area («100 m 2 ) immediately adjacent to the flux tower, but restricted to a small portion of the total flux footprint (typically >1 ha, with variations depending on windspeed and wind direction) [41].This phenology station provided a continuous, high temporal resolution proxy NDVI by means of continuous broad-band reflectance in both PAR and PYR bands sampled every 1-min during the whole growing season and logged as 15-min averages, thus monitoring seasonal changes in the vegetation photosynthetic phenology [43].Additionally, temperature and relative humidity were monitored via appropriate sensors (RH; S-THB-M002, Onset Computer Corporation) allowing rainy-day optical data to be identified and filtered from the dataset (see Section 2.6, below).To derive proxy NDVI values, each sensor pair had been previously cross-calibrated against one another by comparing both upward-and downward-looking sensors in an upward-looking (irradiance) configuration, yielding a coefficient that corrected for any sensor differences (typically less than 5%).These coefficients were then applied to each sensor pair (PAR and PYR) prior to calculating proxy NDVI.Subsequent analyses used the 15-min averages (aggregated over different time periods from 1 to 5 h) to calculate midday average APAR green and proxy NDVI values for comparison with midday average net CO 2 flux from the eddy covariance tower (Figure 3).Since 5-h aggregation periods yielded the highest correlation with CO 2 flux (see Supplementary Materials), and fully encompassed the overpass times of the MODIS Aqua and Terra sensors, we chose that aggregation period for all subsequent analyses.The phenology station data were used to compute a proxy NDVI [27,43] following the formula: where ρPYR is the reflectance of the solar radiation (PYR band), calculated as the ratio of the reflected solar radiation to the incoming solar radiation and ρPAR is the reflectance of the photosynthetically active radiation (PAR band) calculated as the ratio between the reflected PAR and the incoming PAR.This provided a continuous proxy NDVI , which was then averaged over a 5-h midday period (Figures 2 and 3) for comparison with the daily flux measurements.Using the calibration derived from the harvest sites (Equation ( 5)), the proxy NDVI time series was used to derive fAPARgreen for the LUE model (Figure 3).

MODIS NDVI Measurements
The NASA Terra and Aqua satellites, which have orbited Earth since 1999 and 2002, respectively, each carry a MODIS (Moderate Resolution Imaging Spectroradiometer) sensor.These satellites pass daily over most of Earth's surface and provide a nominal 250 m spatial resolution dataset in 36 spectral bands.NDVI values derived from MODIS NDVI products (MOD13Q1 from Terra and MYD13Q1 from Aqua) were downloaded from the Oak Ridge National Laboratory Distributed Active Archive Center website (http://daacmodis.ornl.gov/cgi-bin/MODIS/The phenology station data were used to compute a proxy NDVI [27,43] following the formula: proxy NDVI " pρ PYR ´ρPAR q{pρ PYR `ρPAR q (11) where ρ PYR is the reflectance of the solar radiation (PYR band), calculated as the ratio of the reflected solar radiation to the incoming solar radiation and ρ PAR is the reflectance of the photosynthetically active radiation (PAR band) calculated as the ratio between the reflected PAR and the incoming PAR.This provided a continuous proxy NDVI , which was then averaged over a 5-h midday period (Figures 2 and 3) for comparison with the daily flux measurements.Using the calibration derived from the harvest sites (Equation ( 5)), the proxy NDVI time series was used to derive f APAR green for the LUE model (Figure 3).

MODIS NDVI Measurements
The NASA Terra and Aqua satellites, which have orbited Earth since 1999 and 2002, respectively, each carry a MODIS (Moderate Resolution Imaging Spectroradiometer) sensor.These satellites pass daily over most of Earth's surface and provide a nominal 250 m spatial resolution dataset in 36 spectral bands.NDVI values derived from MODIS NDVI products (MOD13Q1 from Terra and MYD13Q1 from Aqua) were downloaded from the Oak Ridge National Laboratory Distributed Active Archive Center website (http://daacmodis.ornl.gov/cgi-bin/MODIS/GLBVIZ_1_Glb/modis_subset_order_global_col5.pl)The data were derived from a single MODIS pixel (nominally 250 m ˆ250 m) overlapping each eddy covariance tower locations for each site (E3 and E5).MODIS bands 1 (red, 620-670 nm) and 2 (infrared, 841-876 nm) were used for MODIS NDVI calculation.For each MODIS sensor (Aqua and Terra), the 16-day aggregation period was based on the best observations during the composite period, and the actual collection date for each optimal observation was used to provide more accurate time series.The date-corrected time series was then interpolated to produce a daily MODIS value for comparison with daily flux values (Figure 3).

CO 2 Flux Measurements
Net CO 2 fluxes were measured using the eddy covariance (EC) technique [21].Identical EC flux towers were deployed at our E3 (50.8672 ˝, ´111.9045˝) and E5 (50.9057 ˝, ´111.8823˝) grassland sites.Sites exhibited uniform vegetation cover (described above), providing measureable fluxes from all wind directions (i.e., the flux footprint) except those passing through the tower structures.Each tower was equipped with an open-path infrared gas analyzer (IRGA; LI-7500, LI-COR, Lincoln NE, USA) and a three-dimensional sonic anemometer (CSAT3; Campbell Scientific, Logan UT, USA) to quantify vertical CO 2 fluxes.Each EC sensor was affixed at 2.9 m (E3 site) and 3.0 m (E5 site) above ground level and each IRGA was horizontally separated 15-17 cm from each sonic anemometer.All raw EC fluxes were measured at a sampling rate of 10 Hz.Other sensors to quantify environmental conditions and weather were also fixed to each tower to measure half-hour averages of, for example, air temperature, relative humidity (RH) and soil conditions.All collected data were stored by a datalogger (CR5000, Campbell Scientific).
EC data were analyzed using the software package EddyPro (LI-COR, v. 5.1) to quality check raw data, remove outliers and apply standard corrections to calculate corrected vertical fluxes of CO 2 .Raw 10Hz CO 2 fluxes were preconditioned by eliminating data spikes greater than 3.5 standard deviations (σ), temporary drop-outs (10% per bin), and heavily-skewed data (+2 > x < ´2).Anemometer tilt (mean pitch: E3 = ´1.6 ˝; E5: +1.5 ˝) was corrected using a double rotation by forcing lateral and vertical wind components to zero.Data were removed when the time lag between measurements of the anemometer and IRGA exceeded ˘0.4 s.Flux data were further detrended using a block-averaging approach over half-hour time increments.Calculated fluxes were corrected for density fluctuations using the Webb et al. [46] procedure.High-pass spectral corrections were implemented after Moncrieff et al. [47], while low-pass spectral corrections, integrating in-situ conditions to determine system cut-off frequencies, were used after Ibrom et al. [48].Fluxes were further rejected when EC sensors malfunctioned or were affected by moisture, when wind passed through the tower before contacting EC sensors, or when friction velocities fell below 0.1 m¨s ´1 (after Wille et al. [49]).Turbulence tests using the approach of Mauder and Foken [50] were used to remove the poorest-quality fluxes (level 2) when they did occur.Half-hour fluxes below ´3σ or above +3σ were considered outliers and removed.Corrections applied to daytime data resulted in the removal of 23.33% of all calculated fluxes across all EC measurements at both sites.For final comparison with optical data and derived biophysical parameters (proxy NDVI and APAR green ), a midday average net CO 2 flux was calculated based on a 5-h average, matching the period of averaging used for optical data (Figure 4).
were considered outliers and removed.Corrections applied to daytime data resulted in the removal of 23.33% of all calculated fluxes across all EC measurements at both sites.For final comparison with optical data and derived biophysical parameters (proxy NDVI and APARgreen), a midday average net CO2 flux was calculated based on a 5-h average, matching the period of averaging used for optical data (Figure 4).

Filtering and Averaging Process
Preliminary analysis of the optical data indicated invalid proxy NDVI values during periods of rainfall.These periods were easily detectible by evaluating relative humidity (RH) values.Proxy NDVI values associated with RH values greater than 76% tended to be outliers (data not shown) and were discarded.Similarly, automatic gain control (AGC), an IRGA diagnostic value, identified obstructions blocking the sensor's optical path.Water droplets, rime, dew, dust or pollen on the optical path impede sensor function and reduce data quality.Proxy NDVI data with associated AGC values greater than 64% were filtered from the dataset.To achieve data comparable to the MODIS values, a midday 5-h averaging (aggregation) was also applied to both proxy NDVI and flux data.In preliminary analyses, temporal averaging of proxy NDVI over 5 h (Figure 4) reduced short-term variability and clarified seasonal patterns.The impact of temporal aggregation was quantified by calculating the R 2 for the three time intervals considered (1, 3, 5 h midday periods).R 2 values increased from 1-h to 5-h averaging periods considering the whole season; similar results were found when considering green up and senescence periods separately (Table S1, Supplementary Materials).Since the relationship between proxy NDVI and net CO 2 flux showed the best fit using the 5-h aggregation period (Table S1), and encompassed the MODIS overpass times, this time period was used in all subsequent analyses of proxy NDVI.

Gap Filling
One objective was to test the utility of optical data as a proxy of EC flux measurements, thus allowing gap-filling of the flux dataset.Non-gap-filled flux data were augmented by both filtered proxy NDVI data (providing the most continuous time-series of all three NDVIs used in this study, but not necessarily the most spatially representative NDVI) and MODIS data (providing a less continuous, but more spatially representative NDVI).The relationships between these different NDVIs and fluxes were then used to estimate missing flux data by applying several parameterizations of the LUE model (Equation ( 1)), each assuming an invariant efficiency (ε) but using variable NDVI inputs to obtain APAR green .According to the LUE model, we used and compared different model terms: parameterizations ranged from simple NDVI (proxy NDVI and MODIS NDVI) to APAR green .
When comparing proxy NDVI to CO 2 fluxes, both datasets were averaged for a 5-h interval around solar noon (11 am-4 pm; Figure 4) and these values were used as the daily value for flux and proxy NDVI, respectively.These midday average values were then used to model and fill data gaps in the flux data from optical measurements.
When using MODIS NDVI values, NDVI for each day over the entire season was calculated through a linear interpolation of actual date values from MODIS NDVI products (MOD13Q1 from Terra and MYD13Q1 from Aqua).The relationships between daily interpolated MODIS NDVI and net CO 2 fluxes were then used to estimate flux data throughout 2012 for comparison with the proxy NDVI method.
Lastly, the midday average APAR green was calibrated against midday average net CO 2 flux (Figure 3).The relationships between APAR green and net CO 2 flux were used to estimate midday flux data for 2012, and to gap-fill missing flux data.Finally, measured fluxes were compared with the gap-filled, modeled fluxes based on proxy and MODIS NDVI (Figure 3).

Seasonal Patterns of CO 2 Exchange, Above-Ground Biomass, and NDVI Measurements
At both sites, NDVI values obtained from the three different methods all showed similar patterns over the summer growing season, reaching peaks within the first half of July, while green biomass showed maximum values later in the month (Figure 5).An exception was a single MODIS Terra NDVI value in each year (8 July 2012 and 11 July 2013) that yielded one anomalously high NDVI value for each site and date (Table 1, and Figure 5).In most cases (particularly in 2012), the E3 site showed slightly higher productivity than the E5 site, demonstrated by the generally higher NDVI values for E3 (Table 1).In most cases, higher values of NDVIs occurred in 2013 as compared with 2012 (Figure 5, Table 1).A particularly low peak proxy NDVI value in E5 during 2012 (asterisk, Table 1) was the result of a mid-season data gap (spanning from 26 June to 27 August 2012).In general, MODIS Terra values were higher than the corresponding Aqua values, and MODIS values were higher than the ground-based NDVI values, in agreement with previous studies (e.g., [51]).

Green Biomass Estimation Using Different NDVI Formulations
Strong correlations emerged between each of the NDVI variants and green biomass (Figure 6).No significant differences occur between the slopes and the intercepts of the relationships for E3 and E5, thus allowing us to fit a single regression to the NDVI-green biomass data from both sites.Of the three NDVIs, the NDVI 680,800 yielded a slightly stronger correlation (R 2 = 0.82) than the proxy (R 2 = 0.73) or MODIS (R 2 = 0.73) NDVIs (Figure 6).All the relationships were statistically significant, showing p values <0.01.
Remote Sens. 2016, 8, 260 13 of 24 0.73) or MODIS (R 2 = 0.73) NDVIs (Figure 6).All the relationships were statistically significant, showing p values <0.01.c).For this purpose, only results from the E3 site during 2012 were selected because this site and year had the fewest gaps in the optical measurements.To explore the possible effect of hysteresis in our datasets, relationships between flux and NDVI (or APAR green ) are presented for the whole season (single fit) and for the green up phase and the maturity phase separately (separate fits).In this analysis, we used the peak flux value of the season to separate the season in the two parts.Better model estimates (higher R 2 values) were obtained when separate equations were used for the two parts of the season (Figure 7), indicating seasonal hysteresis in the flux-NDVI and flux-APAR relationships.Flux model estimates derived from MODIS NDVI or APAR green (based on proxy NDVI) yielded slightly better agreement with measured fluxes than the estimates based on proxy NDVI (Figure 7).

Parameterization of the LUE Model Using Proxy NDVI, MODIS NDVI and APARgreen
Remote Sens. 2016, 8, 260 14 of 24 datasets, relationships between flux and NDVI (or APARgreen) are presented for the whole season (single fit) and for the green up phase and the maturity phase separately (separate fits).In this analysis, we used the peak flux value of the season to separate the season in the two parts.Better model estimates (higher R 2 values) were obtained when separate equations were used for the two parts of the season (Figure 7), indicating seasonal hysteresis in the flux-NDVI and flux-APARgreen relationships.Flux model estimates derived from MODIS NDVI or APARgreen (based on proxy NDVI) yielded slightly better agreement with measured fluxes than the estimates based on proxy NDVI (Figure 7).The resulting time trends of modeled and measured fluxes using these two approaches (whole-season vs. two halves) are shown in Figure 8. Regardless of the model formulation, a single fit for the whole season tended to underestimate fluxes in the first half of the season and overestimate them in the second, and separate fits improved the agreement (Figure 8).In our final exercise, we used APAR green as independent variable and net CO 2 flux as a dependent variable to derive a modeled flux (Figures 7c and 8c).Unlike NDVI, APAR green also explicitly considers the dynamic variation in PAR irradiance (i.e., PPFD), and thus incorporated hourly and day-to-day variation in illumination.Relative to the proxy NDVI (from which it was derived), the APAR green -based model yielded a better fit with measured fluxes.The APAR green -based model also yielded more dynamic patterns of flux measurements relative to NDVI-based models (cf. Figure 8c to Figure 8a,b).
Remote Sens. 2016, 8, 260 15 of 24 The resulting time trends of modeled and measured fluxes using these two approaches (whole-season vs. two halves) are shown in Figure 8. Regardless of the model formulation, a single fit for the whole season tended to underestimate fluxes in the first half of the season and overestimate them in the second, and separate fits improved the agreement (Figure 8).In our final exercise, we used APARgreen as independent variable and net CO2 flux as a dependent variable to derive a modeled flux (Figures 7c and 8c).Unlike NDVI, APARgreen also explicitly considers the dynamic variation in PAR irradiance (i.e., PPFD), and thus incorporated hourly and day-to-day variation in illumination.Relative to the proxy NDVI (from which it was derived), the APARgreen-based model yielded a better fit with measured fluxes.The APARgreen-based model also yielded more dynamic patterns of flux measurements relative to NDVI-based models (cf. Figure 8c to Figure 8a,b)

Discussion
Eddy covariance measurements, while they provide the "gold standard" for net CO 2 flux measurements, are often affected by a number of problems, including instrumentation malfunctioning or low friction velocity, which can cause long data gaps.For this reason, several methods have been implemented to overcome periods of missing data [52][53][54].The novelty of our approach lies in the particular combination of sensors used and in the comparison of different methods of obtaining NDVI as an input to the LUE model.This demonstrates the possibility of modeling and filling the gaps in a carbon flux dataset using empirical methods based on irradiance and spectral reflectance measurements.We address the data aggregation in the time domain and take advantage of flux station meteorological variables to filter optical measurements, illustrating how one dataset can be used to inform the processing of another.
Over the years, NDVI has emerged as a widely used reflectance index for detecting seasonal changes related to photosynthetic activity [37,[55][56][57][58].The effectiveness of NDVI as a metric of green biomass and as a useful parameter in a simple expression of the LUE model for a grassland ecosystem is confirmed in our study.The strong correlations between NDVI, APAR green , and midday CO 2 flux for this prairie grassland provide a strong foundation for using optical measurements to model fluxes and gap-fill eddy covariance data, as proposed by Wohlfahrt et al. [54].Our findings are in agreement with past studies that have shown good correlations between NDVI and green biomass, APAR green , or total canopy chlorophyll content for many ecosystems [37,[59][60][61].Strong agreements between NDVI or NDVI-based modeled fluxes and eddy covariance measurements, particularly for ecosystems dominated by annual grasslands, are reported by similar studies [37,60,62].These observations illustrate the utility of combining optical and flux measurements to characterize seasonal photosynthetic patterns.
The best relationship with NDVI and green biomass was established with NDVI 680,800 , but similar, strong patterns were found for all NDVI formulations (Figure 6).More work would be needed to understand the slight differences, which are presumably related to the different band formulations, temporal sampling periods, solar angle influence, or footprints of the three NDVI methods.Despite its broad-band formulation, proxy NDVI was able to depict most of the seasonal variability in net CO 2 flux, and had the benefit of continuous sampling (similar to eddy covariance sampling).This outcome highlights the utility of broadband proxy NDVI for monitoring seasonal photosynthetic activity ("photosynthetic phenology") using a set of simple, inexpensive optical measurements [43].These findings confirm that, for the purpose of measuring the seasonal course of photosynthetic activity in grassland ecosystems, broadband dual channel low-cost sensors can replace more complex and expensive hyperspectral sensors at eddy covariance sites when the budget is limited [27,43,63,64].Of the three NDVI methods used, the proxy NDVI had the most limited sampling footprint (estimated at «100 m 2 ), which might explain why it sometimes deviated from the other methods in the depiction of seasonal trends (Figure 5), and why it yielded slightly lower agreement with measured fluxes relative to MODIS NDVI (Figure 7).Further study with mobile sensors or fixed arrays of several proxy NDVI sensors located within the flux tower footprint would possibly improve the agreement with fluxes and clarify this question.
Several versions of the LUE model driven by satellite remote sensing of NDVI have been in wide use for many years [24,31,65,66].While these models can depict broad global patterns, their agreement with local field measurements varies considerably across ecosystems or with different model formulations [67,68].In fact, the selection of input to key parameters (e.g., PAR, fAPAR, APAR, ε) or their precise definitions (e.g., whether APAR and f APAR are based on total or green canopy material) can broadly affect the results [34,69].Site-specific studies incorporating proximal optical measurements can offer insights into these causes of variation.Several recent studies have focused on estimation of carbon dioxide fluxes using proximal optical measurements as proxies for CO 2 fluxes [54,60,63].However, a universal LUE model parameterization to estimate fluxes from proximal optical measurements for all ecosystems, while often discussed, does not yet exist [34,35].Despite these limitations, it is clear that simple models that use linear relationships between fluxes and greenness indices such as those used in our study can successfully be used to estimate seasonal patterns of CO 2 flux for certain ecosystems [70].
In our study, we compared the estimates of carbon fluxes using different predictor variables (proxy NDVI, MODIS NDVI and APAR green derived from proxy NDVI), which reproduced the seasonal net CO 2 flux patterns with varying degrees of fidelity.This experiment allowed us to show that the APAR green -based method improved the estimates over NDVI alone.The proxy and MODIS NDVI methods could not explain all flux variability since they could not capture high-frequency flux changes due to irradiance dynamics.Unlike NDVI, APAR green is not only influenced by canopy structure, but also incorporates PAR irradiance (PPFD) [34].Consequently, the incorporation of PAR irradiance (PPFD) in the relationship better captures the light dynamics affecting fluxes.This is probably one of the main reasons why the APAR green -based model provided a more accurate representation of day-to-day variation in net CO 2 fluxes than NDVI alone (Figure 8).While similar findings have been reported in some cases [71], other studies [64] have sometimes reported the opposite, finding that incorporation of PAR to yield APAR degrades agreement with CO 2 fluxes relative to NDVI or f APAR alone.The reasons for these discrepancies deserve further consideration, and could partly result from methodological differences, or the degree of stress (and variation in light-use efficiency), which can be assessed by the shape of the APAR green -flux curve.Additionally, the period of data aggregation could also affect these relationships, as the effects of including PAR are often more readily visible at fine scales (Figure 8).In general, while validation of modeled productivity using in situ measurements provides essential information on model performance, sampling strategies and methodological differences in data collection complicate comparability [72].One benefit of the proxy NDVI method used here is that it can provide an objective and repeatable sampling method that removes much of the subjective choices often used in field sampling.On the other hand, a limitation of this proxy NDVI method is its relatively limited footprint, which could be addressed by networking multiple sensors.
The simple regression approach used in our study appears to work particularly well in ecosystems characterized by strong seasonal dynamics in above-ground green biomass such as grasslands and croplands [37,59,60,62,63].In our grassland ecosystem, where seasonal CO 2 fluxes largely scaled with canopy development (e.g., as measured by green, above-ground biomass, Figures 5 and 6), an APAR green -based model assuming constant efficiency yielded good agreement with fluxes (Figures 7  and 8), despite our ignoring the efficiency (ε) term in the LUE model.In the literature, there is not a universal agreement on implementing ε, and this lack of a standard approach or definition for ε has recently been addressed by Gitelson and Gamon [34].In past studies, ε has been assumed constant for all vegetation [56], fixed for each vegetation type (varying among different plant communities) [69,73] or variable depending both on vegetation type and stress linked to environmental factors [24].As has been discussed elsewhere [74], we would not normally expect very large variations in efficiency (ε) over the season in grassland ecosystems.For this largely deciduous ecosystem, CO 2 fluxes can be reasonably approximated by f APAR green (estimated by NDVI alone) or by APAR green , allowing us to assume a constant efficiency (ε) over much of the growing season.This supports the hypothesis that, in these annual ecosystems, controls on ecosystem fluxes are largely related to canopy structure, with relatively little physiological control evident [35].These findings are in agreement with a broad range of studies that note strong relationships between CO 2 fluxes and green canopy structure, expressed as APAR green) , LAI, or total canopy chlorophyll [37,59,62,75].In such cases, modeling efforts should focus primarily on accurate estimation of the APAR green term in the LUE model.
We note that not all ecosystems lend themselves to such simple analyses.In many cases, a simple APAR green -based model based may be insufficient, particularly for evergreen ecosystems where dynamic light-use efficiency can be an important determinant of CO 2 fluxes and may have to be included in the model [35,74,76].Additionally, the exact definition and parameterization of the LUE model may itself affect observations regarding light-use efficiency.Many studies define LUE as the flux rate divided by PAR irradiance (incident PPFD) [71,77], causing the efficiency term to be confounded with green canopy structure.In deciduous canopies, this definition of the model terms can lead to large seasonal variations in apparent efficiency that may actually be caused by variations in green canopy structure (the APAR green term of our model), not efficiency as defined in our analysis.As introduced, we defined APAR based on radiation absorbed by green canopies (APAR green ), not incident radiation or total absorbed radiation.This distinction and its consequences has been discussed in more detail elsewhere [34].These differences in model definition and formulation make it difficult to directly compare results across sites and ecosystems, hindering meta-analyses, and confounding larger conclusions beyond the seemingly disparate results of different studies.Consequently, we recommend further attention be given to model formulation with the goal of a standardized definition and parameterization of model terms.
Despite the generally good fit between a simple APAR green -based model and measured CO 2 fluxes, some unexplained features remain.In our study, we found that the fit between CO 2 fluxes and both NDVI and APAR green presented some hysteresis between the first and second halves of the growing season, a finding that agrees with other recent studies at similar annual-dominated sites [61,78].Using different equations for these two growth phases considerably improved the flux estimates over the use of a single equation.Such hysteresis could have several causes.It is likely that ecosystem carbon balance between photosynthesis and respiration shifts over the course of the season, contributing to the hysteresis in APAR green and measured fluxes, and this could also cause a seasonal shift in the APAR green -flux patterns that could be interpreted as a change in light-use efficiency [79].In grasslands, the dead biomass that partially hides new green biomass in the early stage, which is the carry over effect of the previous season's growth, along with the gradually declining chlorophyll and nitrogen content of the leaves during the senescence phase, can also affect the optical-flux relationship.It is possible that our method of partitioning vegetation into "green" and "brown" failed to capture the full dynamics of leaf pigmentation over the season, leading to an apparent change in LUE.Variation in canopy structure (e.g., leaf rolling or leaf inclination angle over the season) could also contribute to this hysteresis.These possibilities warrant further investigation and should be considered in future efforts to model fluxes from optical sensors, particularly if we are to apply the LUE model to all phases of the season.In a similar study, Flanagan et al. [78], documented that ε could decline in the senescence phase due to low soil moisture and high vapour pressure difference (VPD).Such stress in the late season can explain why a unique equation might not be accurate for modeling seasonal CO 2 fluxes.Other studies [80,81] reported that grassland plants quickly respond to changes in moisture conditions and that seasonal patterns of ε are regulated by water availability [82].Thus, ε could decrease from the green-up phase to the maturity phase owing to increasing water stress (and decreasing temperatures) [79].This finding revealed that environmental stresses that are not related to chlorophyll content or change in structure, cannot always be fully detected using NDVI alone [64].We conclude that the hypothesis of a constant ε for our ecosystem seems violated over the entire season (May-October), but holds for half season periods.
Obtaining an adequate optical sampling within the eddy covariance footprint can be a challenging issue [35,45,83].Due to the high variability of the eddy covariance footprint, which varies with wind speed and direction [84], it is not always feasible to directly couple optical and flux footprints, so approximations (e.g., Figure 2) are often necessary.Preliminary results of the flux footprint analysis for our site suggest a typical flux footprint area slightly larger than 1 ha [41].In our study, the 1 ha grid NDVI 680,800 measurements were intended to approximate the footprint of the eddy covariance tower in size, but not necessarily in shape due to the temporal dynamics of the flux footprint.Because NDVI 680,800 involved intensive manual sampling, it was limited to periodic measurement and could not provide a sufficiently dense time series to allow close comparison with the continuous flux measurements.Similarly, the spatial coverage of the MODIS NDVI (pixel size nominally 250 m ˆ250 m) provided a reasonable sample of the eddy covariance footprint but could not capture its full temporal dynamics.MODIS NDVI provided a sufficiently frequent time series that (once interpolated to daily data) allowed comparison with seasonal flux measurements.Of the three NDVI methods, the proxy NDVI from the automated phenology station provided the most detailed time series, matching the sampling frequency to that of the flux measurements.However, the proxy NDVI footprint was limited to a relatively small area (estimated at a few tens of meters squared), which may have reduced its effectiveness in modeling the CO 2 uptake form the much larger and more dynamic flux footprint.These differences in "optical footprint" among NDVI methods undoubtedly influenced our results, presumably contributing to the slight contrasts in seasonal patterns (Figure 5) and the different strengths of correlation between NDVIs and biomass (Figure 6) or net CO 2 flux (e.g., Figures 7 and 8).However, a separate analysis of flux footprint assumptions revealed that the errors due to sampling footprint size and shape were much smaller than differences due to site (E3 vs. E5) or season due to the uniform landscape of this site [41].For other sites with less uniform vegetation stands, the selection of NDVI sampling method and footprint can presumably introduce significant error into the LUE model.
Temporal data aggregation is another important consideration, analogous to the challenge of footprint characterization and matching.Averaging flux and proxy NDVI measurements reduced fine-scale variability of both datasets suggesting that proper temporal aggregation can improve model accuracy when comparing optical and flux data (Table S1, Supplementary Materials).Our findings show that it is crucial to optimize the temporal aggregation when fluxes are related to optical measurements since aggregation can deeply influence the accuracy of the resulting model.For our study, a 5-h averaging period (Table S1) worked well, but we would expect the optimal aggregation period to vary slightly with conditions, e.g., depending upon the stability of the fluxes and optical sampling, which would vary with latitude and daylength across sites.Thus, the topic of optimal aggregation period deserves further study and should be considered as an important issue when predicting fluxes from optical data.
Clearly, one of the challenges of integrating remote and proximal sensing with LUE model in the context of the flux footprint is the need for defining an appropriate temporal and spatial sampling scale.More efforts are needed in the future to further clarify this issue.A detailed investigation of the flux footprint and its dynamics was beyond the scope of this study but was a topic of a separate investigation at this site [41].While sampling footprint error was low for our site, a possible improvement might derive from integrating a high number of wireless proxy NDVI sensors in order to expand the footprint of optical sampling while maintaining high frequency to better match the footprint and time scale of eddy covariance tower.At our site, we are examining this method of improving the spatial representation of proxy NDVI sampling while preserving the rich time series of automated measurements.More attention should also be given to the potentially confounding effects of spatial or temporal autocorrelation, which were not explicitly considered in our study.In the future, UAVs offer a way to improve the footprint of field optical sampling, but obtaining sufficiently rich time series data will remain a challenge.
The integration of optical and flux measurements can add insights and capabilities beyond those of either measurement alone, in part by using one dataset to inform the other, and by addressing the common problem of missing data.Ongoing integration of optical and flux measurements in the framework of the LUE model offers rich potential for improved evaluation of carbon uptake and sequestration for rangeland sites.Comparisons of optical and flux measurements can provide a basis for extrapolating fluxes beyond the scale of the flux tower footprint to larger areas ("upscaling"), a topic of separate research at our site [41].Further attention should be given to standardizing these methods across sites and research groups if we are to reach the full potential of integrating optical with flux measurements to deepen our understanding of ecosystem carbon fluxes.

Conclusions
The main take-home messages of this study exploring the potential of blending different datasets together for monitoring grassland ecosystems are:

Figure 1 .
Figure 1.Map of study location (top left), with detail of study area (bottom left), showing the E3 and E5 calibration sites.Further details of each site are provided on the right.Image credits: ESRI, DigitalGlobe, GeoEye, i-cubed, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

Figure 1 .
Figure 1.Map of study location (top left), with detail of study area (bottom left), showing the E3 and E5 calibration sites.Further details of each site are provided on the right.Image credits: ESRI, DigitalGlobe, GeoEye, i-cubed, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

Figure 2 .
Figure 2. Optical and biomass sampling design around each eddy covariance (EC) and phenology station (PS) location.Canopy-level NDVI680,800 was sampled at regular intervals within a one hectare region surrounding each flux tower using 10 m grid spacing (approx.100 samples; sampling locations represented by black dots).Biomass sampling and fAPARgreen calibration occurred at the locations indicated by numbered circles (1-12), with the sampling years indicated by each circle.

Figure 2 .
Figure 2. Optical and biomass sampling design around each eddy covariance (EC) and phenology station (PS) location.Canopy-level NDVI 680,800 was sampled at regular intervals within a one hectare region surrounding each flux tower using 10 m grid spacing (approx.100 samples; sampling locations represented by black dots).Biomass sampling and f APAR green calibration occurred at the locations indicated by numbered circles (1-12), with the sampling years indicated by each circle.

Figure 3 .
Figure 3. Experimental design, summarizing steps used in derivation of NDVI and APARgreen (arrows), for comparisons with net CO2 fluxes (double lines).As indicated, APARgreen was estimated from the different NDVI metrics by regression of measured fAPARgreen against proxy NDVI, and multiplying fAPARgreen with PAR (photosynthetic photon flux density, measured with PAR sensors).Because the proxy NDVI provided continuous measurements, it was used to derive APARgreen for comparison against Net CO2 flux.Further explanation can be found in the text.

Figure 3 .
Figure 3. Experimental design, summarizing steps used in derivation of NDVI and APAR green (arrows), for comparisons with net CO 2 fluxes (double lines).As indicated, APAR green was estimated from the different NDVI metrics by regression of measured f APAR green against proxy NDVI, and multiplying f APAR green with PAR (photosynthetic photon flux density, measured with PAR sensors).Because the proxy NDVI provided continuous measurements, it was used to derive APAR green for comparison against Net CO 2 flux.Further explanation can be found in the text.

Figure 4 .
Figure 4. Sample diurnal course of proxy NDVI (from phenology station, circles) and net CO2 flux (from eddy covariance, triangles), showing average values (open symbols) calculated for the 5-h midday period (arrow between two vertical lines).See Methods and Supplemental Materials for further details on averaging (aggregation).Data from site E3, 7 June 2012.

Figure 4 .
Figure 4. Sample diurnal course of proxy NDVI (from phenology station, circles) and net CO 2 flux (from eddy covariance, triangles), showing average values (open symbols) calculated for the 5-h midday period (arrow between two vertical lines).See Methods and Supplemental Materials for further details on averaging (aggregation).Data from site E3, 7 June 2012.

Figure 6 .
Figure 6.Correlations between green biomass and (a) NDVI680,800; (b) proxy NDVI (5 h averaged) and (c) MODIS NDVI for both E3 (black dots) and E5 sites (white dots) in 2012 and 2013.Error bars for green biomass are expressed as standard error of the mean.Resulting equations and correlations (R 2 values) are indicated for each fit above.

Figures 7 and 8
Figures 7 and 8 present the different approaches in flux-optical comparison, LUE model testing, and subsequent gap-filling of the flux dataset using optical datasets, including results based on proxy NDVI (panel a), MODIS NDVI (panel b), or APARgreen derived from proxy NDVI (panel c).For this purpose, only results from the E3 site during 2012 were selected because this site and year had the fewest gaps in the optical measurements.To explore the possible effect of hysteresis in our

Figure 6 .
Figure 6.Correlations between green biomass and (a) NDVI 680,800 ; (b) proxy NDVI (5 h averaged) and (c) MODIS NDVI for both E3 (black dots) and E5 sites (white dots) in 2012 and 2013.Error bars for green biomass are expressed as standard error of the mean.Resulting equations and correlations (R 2 values) are indicated for each fit above.

3. 3 .
Parameterization of the LUE Model Using Proxy NDVI, MODIS NDVI and APAR green

Figures 7
Figures 7 and 8 present the different approaches in flux-optical comparison, LUE model testing, and subsequent gap-filling of the flux dataset using optical datasets, including results based on proxy NDVI (panel a), MODIS NDVI (panel b), or APAR green derived from proxy NDVI (panelc).For this purpose, only results from the E3 site during 2012 were selected because this site and year had the fewest gaps in the optical measurements.To explore the possible effect of hysteresis in our datasets,

Figure 7 .
Figure 7. Relationship between proxy NDVI (5 h average, panel a); MODIS NDVI (panel b); APARgreen (panel c) and filtered net CO2 fluxes (5 h average) for the E3 site in 2012.Fits are shown for the whole season (dashed line), for the green-up phase (black dots, solid line), and for the senescence phase (white dots, solid line).R 2 values are reported for each panel considering the whole season (dashed line), the green-up phase (black dots) and the senescence phase (white dots).p values are <0.001 for all the relationships.

Figure 7 .
Figure 7. Relationship between proxy NDVI (5 h average, panel a); MODIS NDVI (panel b); APAR green (panel c) and filtered net CO 2 fluxes (5 h average) for the E3 site in 2012.Fits are shown for the whole season (dashed line), for the green-up phase (black dots, solid line), and for the senescence phase (white dots, solid line).R 2 values are reported for each panel considering the whole season (dashed line), the green-up phase (black dots) and the senescence phase (white dots).p values are <0.001 for all the relationships.

Figure 8 .
Figure 8.Time series of observed and modeled net CO2 fluxes based on proxy NDVI (panel a); MODIS NDVI (panel b); and APARgreen (panel c).The original fluxes are shown as solid lines.Fluxes modeled using a single fit are shown as dashed lines, and fluxes modeled using two separate fits are shown as dotted lines.All results are for the E3 site in 2012.

Figure 8 .
Figure 8.Time series of observed and modeled net CO 2 fluxes based on proxy NDVI (panel a); MODIS NDVI (panel b); and APAR green (panel c).The original fluxes are shown as solid lines.Fluxes modeled using a single fit are shown as dashed lines, and fluxes modeled using two separate fits are shown as dotted lines.All results are for the E3 site in 2012.