Remote sensing-based estimation of gross primary production in a subalpine grassland

Abstract. This study investigates the performances in a terrestrial ecosystem of gross primary production (GPP) estimation of a suite of spectral vegetation indexes (VIs) that can be computed from currently orbiting platforms. Vegetation indexes were computed from near-surface field spectroscopy measurements collected using an automatic system designed for high temporal frequency acquisition of spectral measurements in the visible near-infrared region. Spectral observations were collected for two consecutive years in Italy in a subalpine grassland equipped with an eddy covariance (EC) flux tower that provides continuous measurements of net ecosystem carbon dioxide (CO2) exchange (NEE) and the derived GPP. Different VIs were calculated based on ESA-MERIS and NASA-MODIS spectral bands and correlated with biophysical (Leaf area index, LAI; fraction of photosynthetically active radiation intercepted by green vegetation, fIPARg), biochemical (chlorophyll concentration) and ecophysiological (green light-use efficiency, LUEg) canopy variables. In this study, the normalized difference vegetation index (NDVI) was the index best correlated with LAI and fIPARg (r = 0.90 and 0.95, respectively), the MERIS terrestrial chlorophyll index (MTCI) with leaf chlorophyll content (r = 0.91) and the photochemical reflectance index (PRI551), computed as (R531-R551)/(R531+R551) with LUEg (r = 0.64). Subsequently, these VIs were used to estimate GPP using different modelling solutions based on Monteith's light-use efficiency model describing the GPP as driven by the photosynthetically active radiation absorbed by green vegetation (APARg) and by the efficiency (e) with which plants use the absorbed radiation to fix carbon via photosynthesis. Results show that GPP can be successfully modelled with a combination of VIs and meteorological data or VIs only. Vegetation indexes designed to be more sensitive to chlorophyll content explained most of the variability in GPP in the ecosystem investigated, characterised by a strong seasonal dynamic of GPP. Accuracy in GPP estimation slightly improves when taking into account high frequency modulations of GPP driven by incident PAR or modelling LUEg with the PRI in model formulation. Similar results were obtained for both measured daily VIs and VIs obtained as 16-day composite time series and then downscaled from the compositing period to daily scale (resampled data). However, the use of resampled data rather than measured daily input data decreases the accuracy of the total GPP estimation on an annual basis.

Abstract.This study investigates the performances in a terrestrial ecosystem of gross primary production (GPP) estimation of a suite of spectral vegetation indexes (VIs) that can be computed from currently orbiting platforms.Vegetation indexes were computed from near-surface field spectroscopy measurements collected using an automatic system designed for high temporal frequency acquisition of spectral measurements in the visible near-infrared region.Spectral observations were collected for two consecutive years in Italy in a subalpine grassland equipped with an eddy covariance (EC) flux tower that provides continuous measurements of net ecosystem carbon dioxide (CO 2 ) exchange (NEE) and the derived GPP.
Different VIs were calculated based on ESA-MERIS and NASA-MODIS spectral bands and correlated with biophysical (Leaf area index, LAI; fraction of photosynthetically active radiation intercepted by green vegetation, f IPAR g ), biochemical (chlorophyll concentration) and ecophysiological (green light-use efficiency, LUE g ) canopy variables.In this study, the normalized difference vegetation index (NDVI) was the index best correlated with LAI and f IPAR g (r = 0.90 and 0.95, respectively), the MERIS terrestrial chlorophyll index (MTCI) with leaf chlorophyll content (r = 0.91) and the photochemical reflectance index (PRI 551 ), computed as (R 531 − R 551 )/(R 531 + R 551 ) with LUE g (r = 0.64).
Subsequently, these VIs were used to estimate GPP using different modelling solutions based on Monteith's lightuse efficiency model describing the GPP as driven by the photosynthetically active radiation absorbed by green vegetation (APAR g ) and by the efficiency (ε) with which plants use the absorbed radiation to fix carbon via photosynthesis.Results show that GPP can be successfully modelled with a combination of VIs and meteorological data or VIs only.Vegetation indexes designed to be more sensitive to chlorophyll content explained most of the variability in GPP in the ecosystem investigated, characterised by a strong seasonal dynamic of GPP.Accuracy in GPP estimation slightly improves when taking into account high frequency modulations of GPP driven by incident PAR or modelling LUE g with the PRI in model formulation.Similar results were obtained for both measured daily VIs and VIs obtained as 16-day composite time series and then downscaled from the compositing period to daily scale (resampled data).However, the use of resampled data rather than measured daily input data decreases the accuracy of the total GPP estimation on an annual basis.

Introduction
The availability of simultaneous acquisition of near-surface spectral observations and gas flux measurements quantified with the eddy covariance (EC) technique (Baldocchi et al., 1996) has notably increased in recent years (Sims et al., 2006a;Nakaji et al., 2007Nakaji et al., , 2008;;Hilker et al., 2008a;Cheng et al., 2009;Middleton et al., 2009;Rossini et al., 2010) due to its potential to identify effective links between optical signals and photosynthesis at canopy level (Gamon et al., 2006(Gamon et al., , 2010)).Currently, several research groups have developed different automatic devices to collect canopy spectral properties (Leuning et al., 2006;Hilker et al., 2007;Nakaji et al., 2007Nakaji et al., , 2008;;Daumard et al., 2010;Hilker et al., 2010;Ide et al., 2010;Balzarolo et al., 2011;Meroni et al., 2011) for the purpose of gaining new insights into the quantification and monitoring of plant photosynthesis on a temporal scale.Such devices are generally operated automatically for long periods in the sampling area of flux towers.The increased availability of coupled spectral and flux measurements acquired with comparable temporal and spatial scales has encouraged the revision of existing approaches to modelling photosynthesis and the assessment of the potential for using remotely sensed inputs to spatially extrapolate at landscape level predictions of carbon exchange from information acquired at tower sites.
One of the most widely applied approaches to modelling gross primary production (GPP) based on remote sensing (RS) data is the light-use efficiency (LUE) model proposed by Monteith (1972Monteith ( , 1977)), in which GPP is modelled as a function of the incident photosynthetically active radiation absorbed by vegetation (APAR), determined as the product of the fraction of photosynthetically active radiation absorbed by vegetation (f APAR) and the incident photosynthetically active radiation (PAR), and the conversion efficiency of absorbed energy to fixed carbon (light-use efficiency, ε).
f APAR is usually modelled as a function of VIs.Besides the normalized difference vegetation index (NDVI, Rouse et al., 1974), several recent satellite products or indexes (e.g.enhanced vegetation index (EVI), Huete et al., 2002) have been explored to estimate f APAR.With the advent of hyperspectral RS and the availability of commercial sensors and field instruments, the exploration of a number of different wavelengths and VIs has been promoted to estimate f APAR (Inoue et al., 2008).
A more challenging component of the Monteith model to be inferred from RS is ε.In most LUE models, ε is expressed as a biome-specific constant at its potential maximum, adjusted for unfavourable environmental conditions (e.g.limitations of temperature, humidity, soil moisture, etc.) (Nouvellon et al., 2000;Veroustraete et al., 2002;Heinsch et al., 2006).Some attempts have recently been made to directly infer ε from RS data by exploiting variations in vegetation spectral properties resulting from photoprotection, a process closely linked to photosynthesis.For this purpose, Gamon et al. (1990) originally proposed to exploit changes in reflectance in a narrow-waveband interval centered at 531 nm to track the xanthophyll de-epoxidation state and formulated the photochemical reflectance index (PRI, Gamon et al., 1992): where R ref is a xanthophyll-insensitive reference band.Several studies have demonstrated that ε can be successfully estimated with PRI at leaf (Meroni et al., 2008a), canopy (Evain et al., 2004;Meroni et al., 2008b) and ecosystem (Drolet et al., 2005(Drolet et al., , 2008;;Middleton et al., 2009Middleton et al., , 2011) ) scale.
An alternative approach recently proposed to directly infer ε from RS exploits the link between carbon fixation and sun-induced chlorophyll fluorescence, derived from the oxygen absorption band located at 760 nm (Meroni et al., 2009).Tests of this method have been limited to few studies (Damm et al., 2010;Rossini et al., 2010;Frankenberg et al., 2011;Joiner et al., 2011), and, consequently, the potential of this approach has not yet been fully evaluated.
Another approach to estimating GPP proposed in recent years builds on a simplified version of Monteith's model, which does not need independent estimates of the f APAR and the ε terms.Based on the assumption that chlorophyll is related to the presence of photosynthetic biomass, which is essential for primary production and thus conceptually related to GPP (Sellers et al., 1992), recent studies (Gitelson et al., 2008;Harris and Dash, 2010) suggest that GPP can be estimated through direct correlation with chlorophyll-related indexes.Successful results have been obtained in agricultural crops (Gitelson et al., 2008).In this study, the investigated crops did not suffer short-term environmental stresses.In such conditions, an independent estimate of ε can be unnecessary due to its correlation with chlorophyll content, allowing the use of chlorophyll-related VIs as a proxy of photosynthesis or primary productivity (Sims et al., 2006a).However, these models are unable to model high frequency GPP variations due to changing illumination conditions.To take into account these variations, several studies modelled GPP as the product of VIs and the incident PAR (Gitelson et al., 2006;Wu et al., 2009;Peng et al., 2011).
From this overview of RS approaches currently adopted to estimate GPP, it is evident that, although the ability to model GPP has increased considerably in recent years, a unique model for GPP estimation valid across different ecosystems and a wide range of environmental conditions has not yet been identified (Hilker et al., 2008b;Coops et al., 2010;Penuelas et al., 2011).The proposed research strives to improve our understanding of the links between optical and flux measurements to help developing models suitable for determination of global productivity from space.
In this study, two years of field spectroscopy measurements acquired with an automatic spectral system (Meroni et al., 2011) on a subalpine grassland equipped with an EC tower have been analysed to: (1) evaluate the potential of automatic continuous spectral measurements to monitor the seasonal development of a grassland ecosystem and (2) test the performances of different LUE model formulations driven by RS indexes and meteorological data to estimate GPP.While several studies have evaluated the possibility of modelling grassland GPP based on RS indexes derived from satellite data (Sims et al., 2006b;Li et al., 2007;Harris and Dash, 2010), we are aware of only one study, by Wohlfahrt et al. (2010), that investigated the relationship between EC-derived carbon fluxes and ground measurement of NDVI collected at similar temporal (i.e.daily) and spatial scale in a mountain grassland.In this study, nearsurface spectral measurements were resampled at the same spectral and temporal resolution as the National Aeronautics and Space Administration's (NASA) Moderate Resolution Imaging Spectrometer (MODIS) and the European Space Agency's (ESA) Medium Resolution Imaging Spectrometer (MERIS) onboard Envisat to evaluate the usefulness of currently available global satellite mission observations for modelling GPP by means of the LUE approach.Thus, the research presented in this paper is expected to advance our current ability to monitor and model grassland photosynthesis and it should be useful for the future application of these models to better quantify CO 2 fluxes in different terrestrial ecosystems.

Experimental site
The study site is an unmanaged grassland of the subalpine belt located in the northwestern Italian Alps (45 • 50 40 N, 7 • 34 41 E, Torgnon, Aosta Valley) at 2160 m a.s.l.(Migliavacca et al., 2011).The vegetation of the site is composed mainly of matgrass, and the dominant species are Nardus stricta, Arnica montana, Trifolium alpinum and Carex sempervirens.The area is classified as an intra-alpine region with semi-continental climate with an annual mean temperature of 3.1 • C and mean annual precipitation of about 920 mm (Mercalli and Berro, 2003).The snow-free period lasts generally from late May to early November.

Biochemical and structural field data
Leaf area index (LAI) was determined destructively every two weeks during the two growing seasons (2009 and 2010) at 12 plots of 30 × 30 cm.Collected phytomass was kept on ice and transported to the laboratory.Sample leaves were run through an area meter (Model LI-3100, Li-Cor, Inc., Lincoln NE) and the LAI was determined.Total LAI for the 12 plots was then averaged to obtain a site-level value.Further-more, in correspondence to the 12 plots identified for LAI estimation, a nadir picture of an area of 50 × 50 cm of the canopy (identified by a square positioned on the ground) was acquired every week.The collected images were then analysed with the WinCAM software (Regent Instruments Inc., Quebec, Canada) to classify the percentage of photosynthetic (green) and non-photosynthetic (yellow and dry leaves) components of the canopy during the growing season.The percentage of green components of the canopy derived from image classification was fitted with a forth-order polynomial to obtain the seasonal courses of greenness at daily time-step.
The fraction of photosynthetically active radiation intercepted by vegetation (f IPAR) was computed using measurements from four LI-190 PAR sensors (Li-Cor, Inc., Lincoln NE): one sensor was installed above the canopy at a height of 2.20 m, while three sensors were positioned on ground below the canopy at a distance of approximately 2 m from one another.fIPAR was then computed as where PAR i , PAR t and PAR r are the incident, transmitted and reflected PAR, respectively.In the grassland studied, yellow and/or dead biomass represented a significant fraction of the above-ground biomass during much of the growing season.To adjust the interception of this non-photosynthetic biomass in the calculation of f IPAR, the fraction of standing green vegetation derived from the analysis of nadir pictures was multiplied by f IPAR to give an estimation of "green" or photosynthetic f IPAR (f IPAR g , Hall et al., 1992).
Furthermore, in 2010, leaf samples were collected every ten days at the 12 plots used for LAI estimation.Leaf samples were immediately stored in sealed plastic bags, kept fresh in an ice chest until transported to the laboratory and stored at −80 • C. Leaf pigments were extracted in the following days with N, N-dimethylformamide (DMF) from 100 mg of fresh biomass.The tissue samples were crushed by adding liquid nitrogen, ground in 10 ml DMF for 2 h and then centrifuged (Thermo Electron Corporation Mod.PK110) at 4000 rpm for 25 min to remove particulates.The absorbance of the extracted solutions was measured at 663.8 and 646.8 nm by a Varian UV-Visible Cary100 spectrophotometer.Chlorophyll a and chlorophyll b concentrations per unit leaf mass (µg g −1 ) were then calculated using the extinction coefficients derived by Porra et al. (1989).

Eddy covariance and meteorological data
The turbulent vertical fluxes of CO 2 and latent and sensible heat were measured using the EC technique (Baldocchi et al., 1996).According to EUROFLUX methodology (Aubinet et al., 2000), only half-hourly data in which the theoretical requirements of the EC technique are fulfilled were retained for the following analyses and gap filling techniques were www.biogeosciences.net/9/2565/2012/Biogeosciences, 9, 2565-2584, 2012 used to re-create continuous net ecosystem exchange (NEE) time series.To evaluate temporal variations of CO 2 fluxes and compare these data with spectral measurements, halfhourly measurements of NEE were partitioned into ecosystem respiration and GPP.For the gap-filling and partitioning, the marginal distribution sampling (MDS) method and the partitioning method described in Reichstein et al. (2005), implemented in the online tool (http://www.bgc-jena.mpg.de/bgc-mdi/html/eddyproc/), were used.Different CO 2 flux metrics were used in the analyses: daily midday average GPP (GPP m , µmol CO 2 m −2 s −1 ) computed for the same time period used for spectral properties (11:00-13:00 local solar time) and daily cumulated GPP (GPP d , g C m −2 d −1 ).A detailed description of the EC flux measurements and flux footprint is reported in Migliavacca et al. (2011).Since only PAR absorbed by photosynthetic pigments (approximated with IPAR g in this study) enables photosynthesis processes, to provide more realistic LUE estimates, a "green" LUE (LUE g , Zhang et al., 2009) was computed: Along with EC fluxes, the main meteorological variables were measured with a time step of 30 min; among these, the incident PAR and air temperature were measured above the grassland by means of a quantum sensor (LI-190s, LI-COR Inc.) and a shielded thermo-hygrometer (HMP45C, Vaisala Inc., Woburn MA, USA), respectively.Precipitation was measured using a tipping bucket rain gauge (CS700, Campbell Scientific, Logan, Utah, USA); soil water content (SWC) was measured with water content reflectometers (CS-616, Campbell Scientific, Logan, Utah, USA) installed at two different depths (5-30 cm).

Radiometric measurements and spectral index computation
Canopy radiance spectra were collected using the hyperspectral irradiometer (HSI, Meroni et al., 2011).This instrument is designed for unattended high temporal frequency acquisition of high spectral resolution radiometric measurements.HSI employs a rotating arm equipped with a cosine-response optic to observe alternately the sky and the target surface, thus allowing the computation of the bi-hemispherical reflectance factor (BHR, Schaepman-Strub et al., 2006).HSI uses two HR4000 (Ocean Optics, USA) spectrometers sharing the same optical signal: one covering the visible and near-infrared range (400-1000 nm) with a full width at half maximum (FWHM) of 1 nm; and the other providing higher spectral resolution (0.1 nm FWHM) within a narrower spectral interval (700-800 nm) in the near-infrared.In this study, only the visible and near-infrared spectrometer was used.The spectrometer was spectrally calibrated with a source of known characteristics (CAL-2000 mercury ar-gon lamp, Ocean Optics, USA), while the radiometric calibration was inferred from cross-calibration measurements performed with a calibrated FieldSpec FR Pro spectrometer (ASD, USA).This spectrometer is calibrated by the manufacturer with yearly frequency.Furthermore, the stability of the spectral calibration is regularly assessed during the season using field measured data and the SpecCal algorithm (Meroni et al., 2010;Busetto et al., 2011).The instrument was installed in the proximity of the EC tower at a height of 3.5 m above the investigated surface using a dedicated tower, thus allowing the measurement of the BHR with a nadir viewing geometry.With this configuration, 97 % of the total signal comes from a circular ground area with a radius of about 20 m.Unattended operations were carried out during the snowfree season in 2009 and 2010.During 2009, the instrument was operated between 9 June and 17 October and in 2010 from 20 May to 15 October.Spectral measurements were acquired every 5 min during daylight hours.Only data collected close to solar noon (between 11:00 and 13:00 local solar time) were used for the analyses to minimize changes in solar angle.The spectral system was operated automatically through dedicated software (Meroni and Colombo, 2009).For each acquisition session, the following spectra were collected: spectrometer dark current, incident irradiance, upwelling irradiance and finally incident irradiance again.The target measurement was "sandwiched" between two downwelling irradiance measurements collected some seconds apart.The incident irradiance at the time of target measurement was then computed by linear interpolation.For every acquisition, ten scans were averaged and stored as a single file.
Collected data were processed with a specifically developed IDL (ITTVIS IDL 7.1.1)application.This application allowed the basic processing steps of raw data necessary for the computation of BHR and the application of a set of quality criteria for automatic data selection, described in Meroni et al. (2011).These criteria are intended to identify poorquality data due to unfavourable meteorological conditions (e.g.clouds, rain or fog) or instrumental causes (e.g.problems in the optimization procedure).Whenever one of the quality criteria was not satisfied, the measurement was rejected and excluded from further analyses.
For each retained measurement, canopy reflectance spectra were used to simulate MERIS and MODIS spectral bands, on the basis of the spectral bandwidths and spectral response functions of the two sensors.The different positions of MODIS and MERIS bands are shown in Fig. 1 on the grassland spectra.The list of spectral indexes investigated in this study is reported in Table 1.The NDVI, EVI and PRI spectral indexes were computed from MODIS-simulated data, while the MTCI index was computed from MERISsimulated data.In particular, the MODIS PRI was calculated using the MODIS band 11, centered at 531 nm, which is affected by the xanthophyll de-epoxidation state, and the spectral bands 1 (620-670 nm) (PRI 645 ), 4 (545-565 nm) (PRI 555 ), 12 (546-556 nm) (PRI 551 ), and 13 (662-672 nm) (PRI 667 ) as potential reference bands, in accordance with recent studies (Drolet et al., 2005(Drolet et al., , 2008;;Goerner et al., 2011).
Daily time series of solar-noon spectral indexes were then computed by daily averaging the index values collected between 11:00 and 13:00 local solar time.Sixteen-day composite time series of the different indexes were finally derived from the daily data using the maximum value composite technique to simulate the 16-day dataset routinely produced from MODIS NDVI and EVI and MERIS MTCI acquisitions.The 16-day composite VI time series were then smoothed with a cubic smoothing spline to downscale from the compositing period to daily VI values (Bradley et al., 2007).We will refer to this VI time series as resampled VIs hereafter in this paper.

Testing of different RS models to estimate GPP
Four groups of models with an increasing data requirement and complexity (i.e.number of model parameters) were tested to estimate GPP: Furthermore, the widely used LUE model MOD17 (Heinsch et al., 2006), which is the algorithm used for the MODIS GPP product, was also included in the model comparison.
MOD17 is driven by meteorological variables (minimum air temperature (T min ) and vapour pressure deficit (VPD)), PAR and f APAR.In this study, f APAR was estimated as a linear function of VI g , so the resulting model formulation is where ε max was the maximum radiation-use efficiency (g C MJ −1 ); f (VPD) and f (T min ) varied linearly between 0 and 1 as a consequence of suboptimal temperatures and water availability for photosynthesis.We used site measurements of PAR, VPD and T min to feed the MOD17 algorithm.
To take into account the nonlinear relationship between GPP and the incident PAR (Gilmanov et al., 2007), the inclusion of the logarithm of PAR (ln(PAR)) instead of PAR in model formulations was also tested.Models 1 to 4 were tested using both the measured and resampled VI time series and midday average or daily value of the measured meteorological variables.The performances of MOD17 were evaluated using measured and resampled VI time series and daily meteorological variables.

Statistical analysis
Pearson's correlation analysis was used to test the significance of relationships of VIs and biochemical and structural field data.Model coefficients were derived by fitting each model against both GPP m and GPP d for each day where HSI data were available and for the resampled VI time series.Model coefficients and their relative standard errors were estimated using the Gauss-Newton nonlinear least square optimization method (Bates and Watts, 1988), implemented in the R standard package (R, version 2.6.2,R Development Core Team, 2011).The main fitting (determination coefficient r 2 and root mean square error RMSE) and cross-validated statistics (r 2 cv and RMSE cv ) obtained with the k-fold cross-validation procedure were computed to compare performances of different groups of models.The k-fold cross-validation approach (Hastie et al., 2001) divides the data into k subsets, then the model is fitted using (k − 1) subsets as the training set and the validation is conducted using the omitted subset.In this study, the k subsets were defined by partitioning the dataset into 10 ordinal subsets of equal length (each subset corresponds to 1-month data).This approach is more restrictive than the random definition of the k subsets and was chosen to assess the model performances when large gaps occurred in the data time series (Richardson et al., 2006).Finally, the Akaike information criterion (AIC, Akaike, 1973) was adopted to compare performances of the various model formulations.
For the same period, the average midday air temperature was 19.2 and 19.7 • C in 2009 and 2010, respectively (Fig. 2b).The total amount of precipitation (Fig. 2c) during the snow-free period markedly differed in the two years: 172 mm in 2009 and 362 mm in 2010.The precipitation amount recorded in 2010 was similar to the longterm average (400 mm, 1927-2001) in the same area (Mercalli and Berro, 2003), while 2009 was particularly dry.SWC was strongly related to precipitation inputs during the growing season.As a consequence, the average seasonal SWC in 2010 was higher (24.6 mm 3 mm −3 ) than in 2009 (15.0 mm 3 mm −3 ).In particular, in 2010 there was a precipitation event exceeding 65 mm (DOY 226) that considerably affected SWC.
LAI increased from May and reached its annual maximum in mid-July in both years (Fig. 3a).The maximum LAI in 2009 was 2.7 m 2 m −2 (DOY 194), slightly lower than the 3 m 2 m −2 in 2010 (DOY 201).LAI decreased earlier and steeper in 2010 than in autumn 2009.The variation of leaf chlorophyll content during the growing season was measured only in 2010 (Fig. 3b).The first available sampling was on DOY 176, when the chlorophyll content was already high.It peaked at around DOY 201, as did LAI, and after that it started to decrease.The seasonal pattern of IPAR m (Fig. 3c) showed an increasing trend at the beginning of the growing season: it reached a maximum in about early August and then remained quite stable, with a slightly decreasing course.Therefore, IPAR m failed to detect the reduction of PAR absorbed by the canopy and thus used for CO 2 fixation at the end of the growing season when the canopy was dominated by yellow and dead material.This trend was instead captured by (IPAR g ) m (Fig. 3d).Both IPAR m and (IPAR g ) m were also characterised by considerable day-to-day oscillations due to variations in the ratio of direct to diffuse radiation.

Seasonal variability of spectral data
HSI was operated for 130 days in 2009 (9 June-17 October) and 148 days in 2010 (20 May-15 October).A total of 7331 spectra were collected in the time window used in the present study.Of these, 32.5 % was not considered in the following analyses since they did not fulfil the data quality criteria.Most data were rejected due to instable meteorological conditions, typical of the study site, while only a small percentage was rejected due to instrument failures.
Figure 4 shows the daily time series of midday VIs computed from HSI data.
In both years, measurements started about two weeks after snow melting when the grassland was already greening.As the growing season proceeded, all the VIs except PRI 555 and PRI 551 increased as a result of green biomass accumulation, reaching maximum values in July (around DOY 190) at the same time as maximum LAI and (IPAR g ) m .Then, in the senescent phase of the grassland (from August on), indexes decreased due to plant yellowing and wilting.The patterns of MTCI resembled that of NDVI but, due to the higher sensitivity of MTCI to chlorophyll content with respect to NDVI (Dash and Curran, 2004), it started to decrease earlier and showed year-to-year variability.The EVI dynamics, as compared with other VIs, showed a higher scatter, in particular for high EVI values.This result confirmed a previous study by Miura et al. (2000) that demonstrated that EVI uncertainties tended to increase with increasing VI values and attributed this uncertainty to the inclusion of the blue band in VI formulation for EVI values above 0.4 (between DOY 180 and 225 in our study).PRI 645 and PRI 667 exhibited a pattern similar to other VIs, while PRI 555 and PRI 551 showed an opposite trend characterised by a progressive decrease at the beginning of the growing season up to maximum canopy development and a slower increase in the senescent phase.The most notable differences between the two years analysed were observed in the seasonal dynamics of MTCI and PRI 645/677 between DOY 220 and 250.

Retrieval of biochemical, biophysical and ecophysiological variables from HSI data
The higher sensitivity of MTCI to chlorophyll content was confirmed by the correlation analysis.Chlorophyll was best correlated to MTCI (r = 0.91, p < 0.001) and two PRI indexes using red reference bands (PRI 645 , PRI 667 ) (r = 0.86 and 0.84 , respectively, p < 0.001).NDVI provided a lower correlation (r = 0.80, p < 0.01) whereas the relationship for EVI was not significant (Table 2).NDVI was the VI that related best to LAI (r = 0.90, p < 0.001) and f IPAR g (r = 0.95, p < 0.001).LUE g was best explained by PRI 551 obtained with MODIS band 4 (r = 0.64, p < 0.001); similar results were obtained for PRI 555 with MODIS band 12. Therefore, LUE g was best correlated to PRI indexes based on green reference bands (551, 555 nm), providing results about 20 % better than those obtained using the PRI indexes based on red reference bands (645, 667 nm).

Comparison of VIs and micrometeorological measurements
The comparison between seasonal variations of VIs and variables derived from EC measurements (Fig. 5) showed that the temporal behaviour of VIs related to canopy chlorophyll content (i.e.MTCI, Fig. 4) tracked GPP m quite well (Fig. 5a).
In 2010, both MTCI and GPP m tracked a rebound around DOY 240, probably caused by a rain pulse that occurred on DOY 226 (Fig. 2c).Seasonal courses of PRI 667 and PRI 645 www.biogeosciences.net/9/2565/2012/Biogeosciences, 9, 2565-2584, 2012 Table 2. Coefficients of correlation (r) between the HSI VIs and ancillary and eddy data (LUE g ) measured at the study site.n is the number of samples for each correlation analysis.The VI best-correlated with each variable is in bold print; the second most correlated is in italic.were more similar to those of MTCI than PRI 555 or PRI 551 .
Although PRI 667 and PRI 645 were supposed to be a proxy for (LUE g ) m , the correlation analysis confirmed that in this study they were instead correlated most with leaf chlorophyll content and fIPAR g .The different concavity of PRI 555 or PRI 551 compared to PRI 667 or PRI 645 can be explained by the different position of the reference bands on the grassland spectra (Fig. 1).Bands 4 and 12 (551 and 555 nm band center wavelength, respectively) fell on the peak of vegetation reflectance in the green region, while bands 1 and 13 (645 and 667 nm, respectively) were in the chlorophyll absorption well in the red region of the spectrum.Thus, bands 4 and 12 always had a higher value than band 11 (531 nm) during the growing season.On the contrary, bands 1 and 13 were much higher than band 11 at the beginning and end of the growing season, while they had similar values during maximum canopy development (July and August).(LUE g ) m (Fig. 5b) showed high values at the beginning of the growing season, with a maximum around DOY 165, corresponding to a LAI of 1.5 m 2 m −2 in both 2009 and 2010, and then it suddenly started to decrease.In both years, (LUE g ) m exhibited indistinct seasonality from around DOY 190 to 250, but day-to-day fluctuation, especially on cloudy and partly cloudy days (Fig. 2b) and in correspondence of sharp meteorological events.As an example, the 68 mm precipitation event that occurred on DOY 226 in 2010 caused a sudden increase in SWC and appeared to stimulate (LUE g ) m , which started to increase and reached a value of 0.048 µmol CO 2 µmol −1 photon on DOY 236.From DOY 250 on, (LUE g ) m started to increase probably due to the reduction of the incoming PAR (Fig. 2b).Finally, (LUE g ) m dropped after DOY 280, when the canopy was composed almost entirely of yellow and dead material.A similar seasonal course was observed for both PRI 555 and PRI 551 , thus suggesting that these PRI formulations were the best suited to track (LUE g ) m in this ecosystem.Better performances of these reference bands confirmed previous studies by Middleton et al. (2009) on a Douglas-fir forest and Goerner et al. (2011) on non-boreal/savanna sites.This result supports the hypothesis that the suitability of different reference wavelengths may depend on species composition and stand structure (Gamon et al., 1992;Goerner et al., 2011).

Measured time series
The summary statistics in fitting and cross-validation of the different models tested for GPP estimation are shown in Tables 3 and 4.
Results of model 1 (simple regression analysis) showed that midday VIs explained most of the variability in both GPP m and GPP d : MTCI was the best predictor with a RMSE cv of 1.50 µmol CO 2 m −2 s −1 and 0.74 g C m −2 d −1 , respectively, followed by NDVI and EVI.The inclusion of incident PAR as a multiplicative term of VI g in model formulation (model 2) decreased model performances in GPP m estimation up to a RMSE cv of almost double relative to the corresponding model 1.As an example, RMSE cv of the model using MTCI increased from 1.50 up to 3.30 µmol CO 2 m −2 s −1 and AIC increased from 157 to 417.Similar results were obtained on including the PAR in the form of model 3 for both GPP m and GPP d estimation.Thus, in the majority of cases, the direct use of PAR did not appear to be a useful model component in estimating GPP.On the contrary, results obtained with model 2 and 3 including the logarithm of the incident PAR in the model showed an improvement of the performances in both GPP m and GPP d estimation.The extent of the improvement changed with the different indexes considered.The inclusion of PRI to estimate ε generally increased model performances, in particular when it was used in combination with MTCI and ln(PAR).Model 4 using MTCI, PRI 555 and ln(PAR m ) showed the best performances in estimating GPP m with a RMSE cv of 1.42 µmol CO 2 m −2 s −1 .It is interesting to note that this model also showed the lowest AIC, despite the increase in the number of model variables with respect to model 1.The best-performing model in estimating GPP d was instead model 2 with ln(PAR d ) and MTCI.MOD17, in which ε was expressed as constant ε at its potential maximum adjusted for unfavourable T min and VPD, showed a RMSE cv between 0.78 g C m −2 d −1 for the model driven by MTCI and ln(PAR d ) and 1.57 g C m −2 d −1 for the model driven by EVI and ln(PAR d ).These results were slightly poorer than those obtained on estimating ε as a function of PRI, and, due to the higher complexity of this model, it had a higher AIC.

Resampled time series
Results obtained on performing the same analysis on data aggregated at the 16-day time scale and then downscaled to a daily time step (Tables 5 and 6) confirmed overall those obtained by feeding models with data measured at a daily step.MTCI was the best estimator of f APAR in models 1, 2 and 3 for both GPP m (Table 5) and GPP d (Table 6) estimation.As before, ln(PAR) performed better than linear PAR in models 2, 3 and 4, and the improvement was higher for GPP m estimation.Regarding MOD17, the use www.biogeosciences.net/9/2565/2012/Biogeosciences, 9, 2565-2584, 2012 To obtain an overall view of the capability of different models to represent the seasonal time courses of GPP, we compared EC daily observations (EC-GPP d ) and daily model outputs obtained with the best-performing model for each class fed with both measured and resampled daily inputs (RS-based estimation of GPP, RS-GPP) (Fig. 6).Both measured and resampled daily RS-GPP values agreed quite well with EC-GPP d concerning both amplitude and seasonal phase and successfully described the seasonal dynamics captured by tower fluxes.As noticeable in Fig. 6, the limitation behind the use of resampled rather than measured daily inputs to model seasonal GPP trends was their inability to model GPP day-to-day variations.Even though the statistics in fitting (r 2 and RMSE) and cross-validation (r 2 cv , RMSE cv and AIC) of different models fed with resampled VIs were in most cases better than their daily counterpart (Tables 3  and 4), these models had poorer performances in predicting the sums of daily GPP related to the two growing seasons.Figure 7 shows the sums of daily GPP estimated using the best-performing model for each class.On days for total GPP made it possible to obtain good estimates with both measured and resampled daily inputs.However, absolute average errors in GPP estimation using daily inputs ranged from 0.5 to 1.2 % with models 2 and 1, respectively, and from 1.8 to 2.8 % with models 3 and 1 respectively using resampled data inputs.MOD17 fed by both RS and meteorological inputs produced an average error of 2.4 % in GPP estimation using daily inputs and 1.8 % using resampled data inputs.In general, in 2009 RS-GPP res tended to underestimate the EC-GPP d .This was caused by the inability of RS-GPP res to track the peak of EC-GPP d occurring between DOY 180 and 210 and the recovery of EC-GPP d at the end of the growing season (DOY 260-290).On the contrary, in 2010, RS-GPP res tended to overestimate the EC-GPP d , especially at the end of the growing season.

Discussion
Unattended high temporal and spectral resolution canopy spectra coupled with EC data were acquired for two consecutive years on a subalpine grassland to exploit different strategies for evaluating the potential of RS in estimating carbon uptake.Collected data were processed using automatic procedures which took into account a series of quality criteria related to the illumination conditions during the acquisition and the system performances, and reliable time series of VIs, providing useful information on the time course of different grassland variables, have been obtained.In particular, MTCI was the index most related to chlorophyll content and NDVI to f IPAR g and LAI, confirming previous studies on different ecosystems (Dash and Curran, 2004;Huemmrich et al.,  2).To our knowledge, this is the first study showing the potential of PRI to estimate ε expressed in terms of LUE g , representing a more physiologically realistic way of quantifying the PAR effectively used for photosynthesis compared to ε more widely computed as GPP/APAR or GPP/incident PAR (see the recent review by Garbulsky et al. (2011)).It is worth noting that, as opposed to PRI 555/551 , PRI computed using a reference band positioned in proximity of the chlorophyll absorption well (645 and 667 nm) was more related to leaf chlorophyll concentration than LUE g (Table 2).Therefore, the choice of the reference band used to compute PRI appears to play a key role in the determination of the sensitivity of this index to photosynthetic efficiency.This result confirmed recent studies by Middleton et al. (2009) and Goerner et al. (2011), although we believe that further studies are needed to explore the best reference band for estimating PRI across vegetation types and temporal scales.Furthermore, the translation of these findings to more complex ecosystems (e.g.forests) is not trivial due to the effects of canopy structure on the relationship between PRI and LUE (Barton and North, 2001;Hilker et al., 2008a;Cheng et al., 2010Cheng et al., , 2011)).Most VIs peaked in the first half of July, in correspondence to maximum canopy development, attested by maximum values of LAI and GPP (Figs. 3,4 and 5).However, due to the different sensitivity of VIs to grassland variables, their minimum and maximum values occurred at different DOYs and their slope changed in time.For example, PRI 555 and PRI 551 had a less distinct seasonal course and they reached minimum values about 10-20 days after full canopy development.This time-lag observed between the peak of PRI 555/551 and indexes using red bands can be explained by considering selective light absorption by photosynthetic pigments.Chlorophyll controls the energy flux that can be transferred to the dark reaction of photosynthesis and, because of the lower chlorophyll absorption of green light (Terashima et al., 2009), indexes based on green wavebands may therefore reach their peak later in the season compared to indexes involving a strong chlorophyll absorption band in the red spectral region.
The analysis conducted with LUE models indicated that GPP can be successfully modelled using RS indexes or combining RS indexes with meteorological data.Results of model 1 confirmed that VIs related to canopy greenness, and specifically to chlorophyll content, explained most of the variability in GPP in an ecosystem characterised by a strong seasonality in green-up and senescence such as grasslands and crops (Gitelson et al., 2006;Wu et al., 2009;Peng et al., 2011).MTCI was the best predictor for both GPP m and GPP d , confirming its better performances with respect to EVI in estimating GPP in grassland ecosystems (Harris and Dash, 2010).However, as highlighted by Gitelson et al. (2008), this kind of model is not able to describe variations in GPP due to short-term (hours to days) variations of illumination or environmental stresses (such as temperature and water availability).This limitation was overcome by exploiting models 2 and 3, which take into account variations related to changing incident irradiance.Somewhat surprisingly, the inclusion of incident PAR in model formulation did not result in improved estimation of GPP.However, using ln(PAR) instead of PAR in model parameterization, the accuracy of GPP estimation improved.This means that the grassland increases its efficiency at low values of incident PAR, while, given its moderate LAI and erectophile leaf angle distribution, it is not able to fully exploit high radiation loads.This higher efficiency at low PAR can probably result from more diffuse light scattered within the canopy and less photoinhibition on the top of the canopy, which lead to a reduced tendency toward saturation (Chen et al., 2009).Furthermore, in our case, low PAR conditions can probably be associated with precipitation events, associated with high SWC and low temperatures, which are known to stimulate photosynthetic efficiency in alpine plants (Billings and Mooney, 1968;Korner and Diemer, 1987;Polley et al., 2011).
To account for stress-induced changes in photosynthetic efficiency, the PRI was also tested to directly infer ε from RS data.The inclusion of PRI in model formulation showed slight improvement in GPP estimation, in particular for GPP m .Physiologically, this means that in our ecosystem, APAR g is coupled with ε, and the inclusion of the ε term in the model slightly improves its ability to track seasonal variations.Similar results were obtained by Rossini et al. (2010) and Gitelson et al. (2006) in other ecosystems characterised by strong seasonal variability (crops).Modelling ε as a function of meteorological conditions generally results in lower accuracy in GPP estimation (Table 4).
To evaluate the effect of the temporal resolution of VI time series on GPP estimation, 16-day composite time series of MODIS-(i.e.NDVI, EVI and PRI) and MERIS-derived (MTCI) products were then simulated and downscaled to daily frequency and results were compared.Short-term variability (hours to days) in both VIs and flux data is dampened out by averaging data over two weeks, thus leading to good performances when fitting GPP against resampled VIs (Tables 5 and 6).However, when these models are used to simulate annual GPP, they inevitably provide a decrease in the accuracy of total GPP estimation.The results from models driven only by RS and PAR variables were as good as, and in many cases better than, the more complex MOD17 GPP model, which requires meteorological and vegetation type data inputs in addition to RS indexes.As with several previous studies on VIs, since the estimation of model coefficients is based on a semiempirical regression technique and is conducted only for a single site, further verification studies should be conducted under other vegetation and climatic conditions and at different sites potentially characterised by a more complex structure to fully explore the efficacy of this method and make general inferences.
This study provides a conceptual background for GPP estimation using real satellite data and a better understanding of the spatio-temporal variations of productivity.The choice of the index depends on the spectral characteristics of the satellite sensor being used.In particular, MTCI can be derived from satellite systems with spectral bands in the red edge region (MERIS in this study), EVI and NDVI from satellites having blue, red and near-infrared bands (MODIS in this study) and PRI from satellites with a narrow green band centered at 531 nm (MODIS in this study).Our results show that red edge indexes like MTCI can be used both as single variables or in combination with PRI and meteorological variables to obtain accurate estimations of GPP in a grassland ecosystem.Unfortunately, the computation of MTCI and PRI from a single satellite is currently only feasible from the NASA Earth Exploring One (EO-1) Hyperion sensor, which is near the end of its lifetime with 12 yr in orbit (launched November 2000).The launching of new image spectrometers, such as the NASA HyspIRI or the DLR EnMAP, will allow the calculation of a greater number of indexes, including MTCI and PRI, thus offering significant potential to enhance the accuracy of the assessment of CO 2 uptake in terrestrial ecosystems from space.Finally, we remark that NDVI and EVI showed poorer performances when used as single variables to predict GPP and it is preferable to use these indexes in combination with PRI and meteorological variables to improve accuracy in GPP modelling.

Conclusions
This study investigated the potential of automatic continuous near-surface spectral measurements to monitor the seasonal development of a grassland ecosystem and to evaluate different strategies for terrestrial ecosystem GPP estimation.The main outcomes of this research can be summarized as follows: -continuous field spectroscopy measurements provided reliable information on the seasonal variations of vegetation biophysical and ecophysiological variables with daily temporal resolution.The correlation analysis between VIs and different canopy variables suggested the possibility of using NDVI as an indicator for LAI and f IPAR g (r = 0.90 and 0.95, respectively), the MTCI for leaf chlorophyll content (r = 0.91) and the PRI 551 for LUE g (r = 0.64); -the spectral vegetation index MTCI, designed to be more sensitive to chlorophyll content, explained most of the variability in GPP in the ecosystem investigated, which was characterised by a strong seasonal dynamic of green-up and senescence; -accuracy in GPP estimation improved when taking into account high frequency modulations of GPP driven by incident PAR (in the form of ln(PAR)) or modelling LUE g with the PRI in model formulation; the model formulation that gave the best results in GPP estimation was based on f APAR g (estimated as a function of MTCI) and ε (as a function of PRI 551 ); -results from models driven only by PAR and RS indexes were in many cases better than those obtained from the MOD17 model, which requires meteorological and vegetation type data inputs in addition to RS indexes; -the use of VIs obtained as 16-day composite time series, simulating the 16-day dataset produced from satellite acquisitions, and then downscaled from the compositing period to daily scale rather than measured daily input data, decreased the accuracy of the total GPP estimation on the annual basis.
The approach proposed in this study finds application within the framework of the established SpecNet (Gamon et al., 2010) and recent activities related to the EuroSpec COST action which propose to collect spectral data continuously, regularly and from a worldwide network in connection with the well-established network of flux towers (FLUXNET).Furthermore, improvements in operational LUE algorithms for monitoring global GPP are desirable in the context of efforts to understand trends in global carbon uptake.

Fig. 1 .
Fig. 1.Temporal changes of monthly grassland reflectance spectra collected at midday during 2009.Grey shaded areas represent the position and bandwidth of the MODIS spectral bands: B1 centered at 645 nm, B2 at 858.5 nm, B3 at 469 nm, B4 at 555 nm, B11 at 531 nm, B12 at 551 nm and B13 at 667 nm.White areas represent those of the MERIS sensor: b8 centered at 681.25 nm, b9 at 708.75 nm and b10 at 753.75 nm.
i. model 1, direct linear relationship between GPP and a VI related to canopy greenness (VI g ) GPP = aVI g + b; (4) ii.model 2, direct linear relationship between GPP and the product of a VI g and PAR GPP = a(VI g × PAR) + b; (5) iii.model 3, LUE model assuming constant ε and f APAR estimated as a linear function of a VI g GPP = ε × (aVI g + b) × PAR.(6) Finally, to overcome the limitation of a constant ε, a fourth set of models in which ε is estimated as a linear function of PRI was tested: iv.model 4, assuming ε and f APAR estimated as a linear function of PRI and VI g , respectively

Fig. 6 .
Fig. 6.Time courses of GPP d (g C m −2 d −1 ) estimated from EC measurements (EC-GPP d ) (filled circles), GPP d modelled (open circles) with models fed with measured daily inputs (RS-GPP d ) and GPP d modelled (filled triangles) with models fed with resampled daily inputs (RS-GPP res ) in 2009 (left panels) and 2010 (right panels) for the best performing formulation of each class of models: (a and b) model 1 parameterized with MTCI; (c and d) model 2 parameterized with MTCI and ln(PAR); (e and f) model 3 parameterized with MTCI and ln(PAR); (g and h) model 4 parameterized with MTCI, PRI 555 and ln(PAR); and (i and j) MOD17 parameterized with MTCI and ln(PAR).

Table 1 .
Spectral vegetation indexes investigated in this study: normalized difference vegetation index, NDVI; enhanced vegetation index, EVI; MERIS terrestrial chlorophyll index, MTCI; photochemical reflectance index, PRI.R is the reflectance at the specified wavelength (nm).FWHM is full width at half maximum in nm.

Table 3 .
Summary of statistics in fitting (r 2 and RMSE) and cross-validation (r 2 cv , RMSE cv and AIC) of different models tested in this study using average GPP m and PAR m data.The best-performing models in each class are in bold print.The most successful of all models is grey highlighted.

Table 4 .
Summary of statistics in fitting (r 2 and RMSE) and cross-validation (r 2 cv , RMSE cv and AIC) of different models tested in this study using GPP d and PAR d data.The best-performing models in each class are in bold print.The most successful of all models is grey highlighted.

Table 5 .
Summary of statistics in fitting (r 2 and RMSE) and cross-validation (r 2 cv , RMSE cv and AIC) of different models tested in this study using average GPP m and PAR m data and resampled VI time series.The best-performing models in each class are in bold print.The most successful of all models is grey highlighted.

Table 6 .
Summary of statistics in fitting (r 2 and RMSE) and cross-validation (r 2 cv , RMSE cv and AIC) of different models tested in this study using GPP d and PAR d data and resampled VI time series.The best-performing models in each class are in bold print.The most successful of all models is grey highlighted.