Monitoring of carbon dioxide fluxes in a subalpine grassland ecosystem of the Italian Alps using a multispectral sensor

The study investigates the potential of a commercially available proximal sensing system – based on a 16band multispectral sensor – for monitoring mean midday gross ecosystem production (GEP m) in a subalpine grassland of the Italian Alps equipped with an eddy covariance flux tower. Reflectance observations were collected for 5 consecutive years, characterized by different climatic conditions, together with turbulent carbon dioxide fluxes and their meteorological drivers. Different models based on linear regression (vegetation indices approach) and on multiple regression (reflectance approach) were tested to estimate GEP m from optical data. The overall performance of this relatively low-cost system was positive. Chlorophyll-related indices including the red-edge part of the spectrum in their formulation (rededge normalized difference vegetation index, NDVI red-edge; chlorophyll index, CI red-edge) were the best predictors of GEPm, explaining most of its variability during the observation period. The use of the reflectance approach did not lead to considerably improved results in estimating GEP m: the adjustedR2 (adjR2) of the model based on linear regression – including all the 5 years – was 0.74, while the adj R2 for the multiple regression model was 0.79. Incorporating mean midday photosynthetically active radiation (PAR m) into the model resulted in a general decrease in the accuracy of estimates, highlighting the complexity of the GEP m response to incident radiation. In fact, significantly higher photosynthesis rates were observed under diffuse as regards direct radiation conditions. The models which were observed to perform best were then used to test the potential of optical data for GEPm gap filling. Artificial gaps of three different lengths (1, 3 and 5 observation days) were introduced in the GEP m time series. The values of adj R2 for the three gap-filling scenarios showed that the accuracy of the gap filling slightly decreased with gap length. However, on average, the GEP m gaps were filled with an accuracy of 73 % with the model fed with NDVIred-edge, and of 76 % with the model using reflectance at 681, 720 and 781 nm and PAR m data.


Introduction
In recent years, quantifying and understanding the dynamics and the main drivers of ecosystem carbon dioxide exchange, as well as upscaling the level of observations, have become critical challenges for the environmental scientific community (Canadell et al., 2000;Gamon et al., 2006;Running et al., 1999;Wohlfahrt et al., 2010).
The eddy covariance (EC) technique is a widely and commonly applied method to estimate carbon dioxide exchange between vegetation and the atmosphere at the ecosystem scale (Baldocchi, 2003;Burba, 2013;Geider et al., 2001).Although this method is able to provide direct, near-continuous and high-temporal resolution measurements of net gas exchange, it also has some limitations.
EC technique provides flux measurements of a relatively small area.The flux "footprint" varies from tens of meters to Published by Copernicus Publications on behalf of the European Geosciences Union.
several kilometers and depends on many parameters such as measurement height, wind velocity, surface roughness and atmospheric stability (Baldocchi, 2003;Kljun et al., 2001;Schmid, 1994).At the same time, the EC systems are relatively expensive -the typical cost for a complete EC system is on the order of USD 40-50 thousand, and the cost of site infrastructure is additional (Running et al., 1999).Considering all of these aspects, it is clear that, although EC measurements can be considered a solid basis for the ecosystem-scale CO 2 flux measurements, complementary methods are needed to extend the estimates to landscape and regional scales.
Important networks such as SpecNet, IMECC, and EU-ROSPEC have been investigating the potential of coupling spectral and EC observations (Balzarolo et al., 2011).In situ measurements can provide unique data sets with high spectral, spatial and temporal resolution, which represent a solid basis for validation of remote observations carried out at aircraft and satellite levels and further upscaling (Gamon et al., 2006;Gamon at al. 2010).As a result the number of sites where direct flux measurements are conducted simultaneously with in situ spectral measurements has increased significantly within the last decade.
The most commonly used approach to estimate the gross ecosystem production (GEP; µmol m −2 s −1 ) with proximal sensing is based on the light-use efficiency (LUE) model proposed by Monteith (Monteith and Moss, 1977;Monteith, 1972).This simple model assumes that GEP is driven by the absorbed photosynthetically active radiation (APAR; µmol m −2 s −1 ) and the photosynthetic radiation use efficiency expressing the carbon sequestration efficiency per amount of the absorbed solar energy (ε; µmol CO 2 µmol −1 APAR): where PAR is the incident photosynthetically active radiation (µmol m −2 s −1 ) and f APAR is the fraction of PAR absorbed by the vegetation canopy (%).
In non-stressed ecosystems characterized by strong seasonal dynamics such as some managed croplands, independent estimates of ε may be unnecessary due to its relation with the chlorophyll content (Gitelson et al., 2012;Peng and Gitelson, 2012;Peng et al., 2011;Rossini et al., 2012;Wu et al., 2009), and this is particularly true when integrating GEP over longer timescales, e.g., days (Gitelson et al., 2008).Therefore most of the variations in plant productivity in such ecosystems should be reflected by changes in APAR (Lobell et al., 2002).
Several studies modeled GEP as a function of VIs (Harris and Dash, 2010;Rossini et al., 2010;Sims et al., 2006;Sjöström et al., 2009;Xiao et al., 2004) and/or of VIs multiplied by PAR (Gitelson et al., 2006;Peng and Gitelson, 2012;Peng et al., 2011).Including PAR in the model should theoretically enhance the correlation with GEP, because the product of VI and PAR takes into account the seasonal changes in both biophysical parameters controlling the photosynthesis process (e.g., f APAR and chlorophyll content) and in the amount of radiation reaching the vegetation surface (Gitelson et al., 2012).
In the current study, 5 years of field multispectral data acquired with the CROPSCAN MSR16R system (CROP-SCAN Inc., Rochester, USA) deployed on the EC tower of the FLUXNET grassland site IT-MBo (Viote del Monte Bondone, Trento, Italy) are presented and analyzed.
In particular, the objectives of this paper are the following: (i) to investigate the potential of vegetation reflectance and narrow-band VIs for monitoring carbon dioxide fluxes exchanged between the dynamic grassland ecosystem and the atmosphere (ii) to analyze the relationships between spectral data and carbon dioxide fluxes during the 5 years of observations in order to determine how robust the relationships between vegetation spectral properties (reflectance and narrow-band VIs) and mean midday GEP (GEP m ) are (iii) to compare different approaches (correlation analysis and multiple regression) to estimate GEP m (iv) to evaluate the potential of spectral models to gap-fill GEP m data.

Experimental site
The study site is a permanent alpine grassland located at 1550 m a.s.l. on the Viote del Monte Bondone plateau (46 • 00 N, 11 • 02 E; Italian Alps).The vegetation of the area is dominated by Festuca rubra (L.) (covering 25 % of the area), Nardus stricta (L.) (13 %) and Trifolium sp.(L.) (14.5 %), and is representative of a typical low-productive meadow of the Alps.The site is managed as an extensive meadow with low-mineral fertilization (applied in autumn) and is cut once a year, usually in mid-July (Gianelle et al., 2009).The maximum canopy height at the peak of the growing season (mid-June to early July) can reach approximately 30 cm.
The climate of this area is sub-continental (warm and wet summer) and is characterized by a mean annual temperature of 5.5 • C, with monthly averages ranging from −3.1 • C in February to 14.3 • C in July.The annual mean precipitation is 1244 mm, with maximum values in May (138 mm) and October (162 mm).The snow-free period lasts typically from early May to late October (Marcolla et al., 2011).
The site is characterized by a regular east-west wind circulation, showing along this direction an almost flat topography with a homogeneous vegetated fetch of more than 500 m.An experimental footprint analysis demonstrated that 30 % (in stable atmospheric conditions) to 80 % (in unstable conditions) of the total CO 2 flux originates within 30 m from the EC tower (Marcolla and Cescatti, 2005).

Eddy covariance and meteorological data
Continuous EC measurements of CO 2 , water vapor and sensible heat fluxes were performed at the Monte Bondone FLUXNET site from the beginning of August 2002.In the present study, data from 2008 to 2012 were used, to match the available spectral data set.
The eddy covariance (EC) system consisted of a Licor Li-7500 open-path infrared gas analyzer (Li-COR Inc., Lincoln, Nebraska, USA) and a Gill R3 3-D ultrasonic anemometer (Gill Instruments Ltd., Lymington, UK), mounted at a height of 2.5 m.Raw data were recorded at a frequency of 20 Hz and stored by means of the EDISOL software package (Moncrieff et al., 1997).The EdiRE software (version 1.4.3.1021,R. Clement, University of Edinburgh) was used to compute turbulent CO 2 fluxes from the raw data.
Along with EC flux measurements, the main meteorological and soil physical variables were measured, including the following: short-and long-wave radiation components (Kipp & Zonen CNR1, Delft, the Netherlands), incoming total and diffuse PAR (LI-COR LI-190SA, Lincoln, USA; and Delta-T BF3H, Cambridge, UK), precipitation (Young 52202H, Traverse City, Michigan, USA), air humidity and temperature (Rotronic MP103A, Crawley, UK), soil temperature profile at depths of 2, 5, 10, 20 and 50 cm (STP01, Hukseflux, Delft, the Netherlands), and volumetric soil water content at depths of 10 and 20 cm (CS615 reflectometers, Campbell Scientific Inc., Logan, Utah, USA).All meteorological variables were recorded at 1 min intervals and averaged over 30 min; both 1 min data and half-hourly averages were stored on a CR23X data logger (Campbell Scientific Inc., Logan, Utah, USA).
Half-hourly measurements of net ecosystem exchange (NEE) were gap-filled and partitioned into ecosystem respiration (Reco) and gross ecosystem production (GEP) by means of the online tool developed by Reichstein et al. (2005) (http://www.bgcjena.mpg.de/bgcmdi/html/eddyproc).However, only non-gap-filled data were analyzed in this study.To maintain consistency between the time window used for calculating vegetation reflectance and narrow-band VIs, the mean midday gross ecosystem production (GEP m , µmol m −2 s −1 ) and mean midday incoming photosynthetically active radiation (PAR m , µmol m −2 s −1 ) were calculated for the same time period used for vegetation spectral properties (11:00 a.m.-1:00 p.m. of local solar time).
Further details regarding the EC instrumentation, data elaboration and quality control can be found in Marcolla et al. (2011).

Multispectral reflectance and narrow-band vegetation indices
Multispectral data were acquired on a continuous basis from 2008 to 2012 by means of the CROPSCAN multispectral radiometer system MSR16R (CROPSCAN Inc., Rochester, USA).The system consists of a 16-band radiometer (simultaneously measuring reflected and incoming radiation in narrow spectral bands) and a data logger controller (DLC) storing the acquired data (Table 1).For each band, the incoming solar irradiance is measured through a cosine diffuser, while reflected radiance is measured through a 28 • field of view foreoptic.The system was installed on the existing EC tower at a height of 6 m, which allowed the observation of a 3.0 m diameter vegetation surface.Before the beginning of each growing season, the system was calibrated using the method recommended by the manufacturer, based on the use of a white reference panel with known reflectance (http://www.cropscan.com/wsupdn.html).Additionally, CROPSCAN, Inc. provided cosine response calibration data with each upward-facing MSR16 module and temperature sensitivity calibration data.Both cosine and temperature corrections were included in the postprocessing software (POSTPROC program) provided with the MSR system.
Incident irradiance and reflected radiance were collected every 10 min, and reflectance at given wavelengths was calculated.In order to minimize solar angle effects, reflectance data were finally averaged over 2 h close to a solar noon (11:00 a.m.-1:00 p.m. of local solar time).
Due to the noisy and unreliable optical signal beyond 1000 nm (bands no. 15 and 16; Table 1), only the data of the first 14 bands were included in the analyses.In addition, data were excluded when (1) the site was covered by snow, (2) precipitation was recorded 2 h prior or during the midday averaging period, and (3) the weather conditions did not allow for the removal of the cut biomass from the footprint of CROPSCAN system (and EC tower) straight after the cut event.According to these quality criteria, 24 % of the data were discarded, mainly due to the meteorological conditions.
Canopy reflectance spectra were then used for computing the VIs.Although many different VIs were investigated (Table A1), only the most commonly used and the best performing in GEP m estimation -considering all the 5 years of observations -are presented in the study.The list of the five presented VIs is reported in Table 2.

Models for GEP m estimation
In order to estimate GEP m we used two approaches: one based on linear regression (using the concept of the LUE model, i.e., Eq. 1) and the other on multiple regression.The first approach assumed a direct linear relationship between GEP m and VIs (model 1) and between GEP m and the product of VIs and PAR m (model 2).In the second approach, the interaction effects between different variables were explored by running two stepwise bidirectional multiple regression models, in which GEP m was set as a dependent variable and reflectance (model 3), or reflectance and PAR m (model 4), as explanatory variables.The abovementioned models (Table 3) were tested both for each year on a separate basis, and for all the years together in order to obtain the general models for the estimation of GEP m .

Statistical analysis
Pearson's correlation analysis was used to test the significance of the relationships between GEP m and VIs • PAR m .
In order to evaluate how robust the relationships between GEP m and VIs were, the slopes of the linear regressions be-tween the best performing VI against GEP m were analyzed.In particular, the slopes of the regressions obtained for each year and obtained in the general model 1 (including all 5 years) were compared by means of a t test to check whether the regression coefficients were statistically different.
Besides, a multiple stepwise bidirectional linear regression was used to explore the interaction effects between variables (considering GEP m as a dependent variable and reflectance at 14 analyzed wavelengths (model 3), or reflectance values and PAR m (model 4), as explanatory variables) to find the model that best fits the data according to Akaike's information criterion (AIC; Akaike, 1973).The variance inflation factor (VIF; Mason et al., 2003) was used to measure the degree of (multi)collinearity of the ith independent variable with the other independent variables in the regression models.
When VIF for any of the predictors reached the threshold value of 10, the (multi)collinearity was reduced by eliminating one independent variable (the last one selected by the automatic stepwise bidirectional regression) from the analysis (O'Brien, 2007).The procedure was repeated until none of the VIF factors exceeded the acceptable threshold value; thus the subset of explanatory variables was free of significant (multi)collinearity issues.
The final subset of the predictor variables was selected by testing whether the increase of the adjusted R 2 (adjR 2 ) after adding a subsequent predictor variable to the multiple regression model was significantly different from zero (at significance level α = 0.001).Multiple regression models were compared by means of the Fisher test.
Each of the four model's coefficients was obtained by fitting each model against GEP m .The main goodness of fit statistics (adjusted coefficient of determination -adjR 2 , root mean square error -RMSE, percentage root mean square error -PRMSE and probability valuep) were computed to compare the performance of the different models.
Additionally, a validation of the best performing general models using training/validation splitting approach, in which 1 year at a time was excluded from the data set, was conducted.The remaining 4-year subset was used as a training set and the excluded year as a validation set.The model was fitted (calibrated) against each training set and the resulting parameterization was used to predict the GEP m of the excluded year.Validation accuracy was evaluated in terms of RMSE.
All the statistical analyses were performed by means of the R software (version 2.15.2, http://www.r-project.org).

The gap scenarios
In order to evaluate the ability of spectral models to gap-fill CO 2 flux data, secondary data sets were generated by flagging ∼16 % of the five growing seasons data as unavailable (artificial gaps constituting 90 observation days out of 573 available observation days).The percentage of artificial gaps  Table 3.The four models for GEP m estimation tested in the presented study.
Model Model formulation was chosen due to the fact that during the observation period of the study (∼May to November, 2008-2012) the EC data set had an average of 16 % of missing or rejected values of NEE data collected during midday hours.Following Moffat et al. (2007) these artificial gaps were superimposed on the already incomplete data, without regard for the distribution of real gaps in the time series.Three gap length scenarios were considered: gaps of 1, 3 and 5 observation days.The artificial gaps were distributed randomly, and each of the three scenarios was permuted 10 times and results were averaged (Moffat et al., 2007).Secondary data sets with artificial gaps were used to calibrate the models that were applied for filling GEP m data.The gap-filling statistical metrics (adjR 2 , RMSE, PRMSE) were calculated using the EC derived GEP m in these artificial gaps to validate the predictions of filling technique.

Results
Figure variability was observed in total precipitation recorded from May to November (Fig. 2).The differences in precipitation sums between the investigated years reached up to 50 %.The precipitation amount in 2011 (1008 mm) was similar to the 20-year period average (990 mm, 1993-2012).The growing season of 2010 (1473 mm) was particularly wet, with a precipitation sum 49 % higher than the long term average, while 2009 (744 mm) was fairly dry, with a total precipitation Seasonal patterns of GEP m were driven by both environmental variables (such as incoming PAR and air temperature) and grassland management (Marcolla et al., 2011).The grassland cut occurred around mid-July, and split the growing season into two sub-periods.The maximum gross CO 2 flux rates were recorded in the early summer (end of Junemid July).After the cut event, the canopy regrowth generally reached a peak at the beginning of September.
The VIs showed a similar behavior to GEP m , and the peaks of these time series were almost synchronous.Starting from the early part of September, VIs began decreasing gradually in all the investigated years due to the senescence phase (characterized by a progressive canopy yellowing and wilting), but at varied rates.
Examples of seasonal courses of investigated VIs and GEP m measured in 2012 are shown in Fig. 3.For better visualization and easier comparison, both GEP m and VIs were normalized by scaling between 0 and 1.The graphs which refer to other years of observations can be found in Appendix Fig. B1.
The linear regression analysis (Table 4) showed that the presented VIs explained at least 50 % of the variability of GEP m .
The highest accuracy of model 1 was obtained in 2009 and 2012 (adjR 2 up to 0.81).On the other hand, the lowest accuracy of the same model was reported in 2011 (max adjR 2 = 0.64).This low value of adjR 2 could be explained by the fact that during this year the CROPSCAN sensor was not operated during the autumn period, and thus the range of VIs and GEP m was smaller as the senescence phase was missed (Table 4).The estimation accuracy was also dependent on the VIs used for the parameterization of model 1 (Table 4).VIs, including the red-edge band in their formulation, turned out to be the best candidates for GEP m estimations considering both the general model and the 5 different years on a separate basis.The MSR, although it is based on the NIR and red bands, also showed reliable performance.Taking into account the models for the single years, MSR, DR, and CI red-edge were included in the group of the three best fitting models 3, 2 and 4 times, respectively.NDVI red-edge was in the group of the three best performing models in each investigated year.On the contrary, NDVI was never included among the best predictors of GEP m (Table 4).
The inclusion of incoming PAR m into the model resulted in a general decrease of its performance.The PRMSE was on average 14.64 % higher in model 2 than in model 1 considering all of the 5 years of observations.As an example, the adjR 2 of the general model (2008-2012) fed with NDVI red-edge decreased from 0.74 to 0.61, RMSE increased from 3.41 to 4.19 µmol m −2 s −1 and PRMSE increased from 16.40 to 20.18 %.A similar pattern was observed in each of the investigated years (Table 4).
In order to investigate the impact of radiation quality on these results, the light response of half-hourly GEP (data collected between 11:00 a.m. and 1:00 p.m; during the snowfree period of 2012) considering different levels of diffuse radiation was investigated.Two different relationships between GEP and incoming PAR were found: one for cloudy conditions (when diffusion index -DI, which is the ratio between diffuse and total incident PAR, exceeded 0.7) and one for sunny conditions (DI < 0.3) (Fig. 5).The data when the abovementioned populations were overlapping (PAR from  800 to 1350 µmol m −2 s −1 ) indicated that, in the Monte Bondone grassland site, photosynthesis rates were significantly higher under diffuse compared to direct radiation.A stepwise bidirectional procedure selected reflectance (R) at 681, 781 and 720 nm (model 3) and R681, R781, PAR m and R720 (model 4) as significant drivers of GEP m , considering each of the 5 years of observations simultaneously (Table 5).
It is interesting to note that in both model 3 and 4, referring to each observation year on a separate basis (data not shown), the red-edge bands were included as important predictors in all of the 5 investigated years.In model 3 the red region was chosen as a highly predictive variable in 40 % of cases, while the NIR region in three out of five investigated growing seasons.In model 4, red and NIR bands contributed to the stepwise regression model in 3 and 2 out of 5 observation years, respectively.PAR m , as an additional variable of model 4, was included in the model three out of five times.The range of adjR 2 values for different years considered on a  A stepwise bidirectional multiple regression with reflectance at 681, 781 and 720 nm as predictors did not yield any improvement in the explained variance of GEP m when the entire data set was considered (adjR 2 =0.74 -general model 1; adjR 2 =0.73 -general model 3; Table 4 and 5, respectively).Also, adding PAR m as an independent variable of the model resulted only in a slight improvement in the accuracy of the GEP m estimation compared to the general linear regression model 1 based on NDVI red-edge .In fact, the adjR 2 increased from 0.74 to 0.79, while the PRMSE decreased from 16.40 to 14.75 % (Tables 4 and 5).
Validation of model 1 based on NDVI red-edge showed that there was no relevant difference in prediction accuracy among validation years (RMSE was varying between 3.12 and 3.85 µmol m −2 s −1 , Fig. 6).Validation results of general model 4 showed that considering all the 5 validated years RMSE was on average 3.26 µmol m −2 s −1 .
The differences in the adjR 2 performance of the gapfilling scenarios showed that the accuracy of gap filling decreased slightly with gap length, while the range of the goodness of fit statistics (adjR 2 , RMSE, PRMSE) generally increased with gap size (Table 6).However, on average, GEP m gaps were filled with an accuracy of 73 % with model 1 fed with NDVI red-edge (RMSE = 3.40 µmol m −2 s −1 , PRMSE = 16.48 %), and with an accuracy of 76 % (RMSE = 3.14 µmol m −2 s −1 , PRMSE = 15.25 %) with model 4 using reflectance at 681, 720 and 781 nm and PAR m data.

Discussion
Continuous and simultaneous measurements of narrow-band canopy reflectance and EC carbon dioxide fluxes have been successfully performed for 5 consecutive years in a subalpine grassland ecosystem.The multispectral CROPSCAN MSR16R system demonstrated to be a reliable instrument for monitoring carbon dioxide fluxes.The results of this study provided important information on how consistent and robust the relationships between VIs and GEP m are in such a dynamic ecosystem.Additionally, they allowed the comparison of different approaches (correlation analysis and multiple regression) for predicting GEP m .
Although several studies have already compared VIs obtained from in situ observations against EC CO 2 fluxes (Gitelson et al., 2003b;Inoue et al., 2008;Peng and Gitelson, 2012;Peng et al., 2011;Rossini et al., 2010;Sims et al., 2006), and a few studies have focused on very similar canopies (Gianelle et al., 2009;Rossini et al., 2012;Wohlfahrt et al., 2010), we are not aware of any study based on such a long time series, acquired on a continuous basis during the growing seasons.
From the data presented, it follows that MSR and DR indices, which are modified and improved variants of the most commonly used VIs, showed generally a slightly stronger linear relationship with GEP m when compared to NDVI.Nevertheless, considering all of the observation years, the most robust estimates of GEP m were obtained when NDVI red-edge and CI red-edge were used to parameterize the model (Table 4).These results confirmed the findings of previous studies on both similar (Rossini et al., 2012) and different ecosystems (Gitelson et al., 2003b;Peng and Gitelson, 2012;Peng et al., 2011;Rossini et al., 2010), indicating that VIs based on the red-edge part of the spectrum are the most sensitive to the seasonal GEP dynamics due to their better linearity with chlorophyll content (Gitelson et al., 2003a;Sims and Gamon, 2002;Wu et al., 2008), and with green leaf area index -green LAI (Gitelson et al., 2003c;Viña et al., 2011).In general, VIs (such as NDVI) calculated as a normalized difference between NIR bands -characterized by a high reflectance due to leaf and canopy scattering, and visible bands (e.g., red), where absorption by the chlorophyll pigments is predominant (Jackson and Huete, 1991) -tend to lose their sensitivity to moderate-high aboveground biomass due to the saturation of reflectance in the visible bands and due to the limitation of the normalized difference approach www.biogeosciences.net/11/4695/2014/Biogeosciences, 11, 4695-4712, 2014  ( Fava et al., 2007;Gao et al., 2000;Mutanga and Skidmore, 2004).Better performances of NDVI red-edge and CI red-edge stem from the fact that even though the red-edge part of the spectrum is characterized by lower absorption by chlorophyll, it still remains sensitive to changes in its content, reducing the saturation effect and enhancing the sensitivity of these VIs to moderate-high vegetation densities (Clevers and Gitelson, 2013;Wu et al., 2008).
Incorporating PAR m into the model resulted in a general decrease in the goodness of fit of the linear regression.One reason for this is that sunlight is used by plants more efficiently under cloudy than clear sky conditions due to a more uniform illumination of the canopy, and thus a smaller fraction of the canopy likely to be light saturated (Baldocchi and Amthor, 2001;Chen et al., 2009;Mercado et al., 2009).Accordingly, significantly higher photosynthesis rates under diffuse as regards direct radiation conditions (with similar values of PAR) were noted in the Monte Bondone site (Fig. 5).Similar results have been reported by Rossini et al. (2012), who also pointed out that, in a similar subalpine grassland ecosystem, the inclusion of incident PAR in a model formulation did not result in an improved estimation of GEP.However, in several other studies referring to other dynamic ecosystems, GEP was successfully estimated as a product of VIs and PAR (Peng and Gitelson, 2012;Rossini et al., 2010;Wu et al., 2009).A recent study of Peng et al. (2013) confirmed that the use of PAR in the model can introduce noise and unpredictable uncertainties in GEP estimations.As suggested by these authors, the response of productivity to changes in PAR is quite complex and is influenced by many variables such as vegetation physiological status, canopy structure and light distribution in the canopy.Some other authors also brought to light some important aspects related to the use of PAR.Sims et al. (2008) showed that the variation in PAR is a more relevant determinant of GEP over very short timescales, and appears to be important for diurnal trends.Gitelson et al. (2012) demonstrated that seasonal variation of PAR potential (defined as the maximal value of incident PAR that may occur when the concentrations of atmospheric gasses and aerosols are minimal) can be used to improve the performance of the models.Therefore, further analyses of the response of different vegetation types to various levels of diffuse radiation are required, and the hypothesis that the DI and PAR potential can improve the performance of the models including radiation as an input parameter needs to be verified.
The use of the reflectance approach instead of the VI approach did not lead to considerably improved results in estimating GEP m .Including additional predictors in multiple stepwise regression resulted in only a 6 % improvement of the explained variance, considering all of the 5 years of observations collectively.We believe this was partly due to the limited number of available bands of the CROPSCAN system, and that further studies are needed to explore the benefits of using hyperspectral data for predicting CO 2 uptake across different terrestrial ecosystem types.
A detailed analysis of the full vegetation spectrum and of the various spectral absorption features appears to be particularly meaningful for providing a solid basis for upscaling of GEP estimations using airborne and satellite platforms.
In this study the reflectance value at 720 nm, which was used in the multiple regression models, did not bring a relevant increase in the adjR 2 values (partial adjR 2 was 0.04 and 0.03 for model 3 and 4, respectively).On the other hand, the successful performance of VIs using this band confirms the important role of this part of the spectrum in monitoring the dynamics of ecosystem carbon dioxide fluxes.
Validation results of general model 1 fed with NDVI red-edge showed that RMSE increased on average from 3.41 to 3.48 µmol m −2 s −1 , compared to non-validated general model 1 (averaging the values obtained from the 5 different validation years).Validation results of general model 4 showed that RMSE increased on average from 3.06 to 3.26 µmol m −2 s −1 , compared to non-validated general model 4. The highest decrease of the GEP m estimation accuracy was noted in the growing season of 2012 (Table 4, Fig. 6), which was presumably caused by the unusual drought which occurred just after the cut event.The precipitation to temperature ratio for a 15-day period after the cut in the growing season of 2012 was more than 10 times lower than in the other years, and this fact could have affected GEP m to a higher extent than VIs related to canopy "greenness".As a consequence, models calibrated with the first 4 years of the data set overestimated the GEP m measured in the second part of the growing season of 2012.
During the observation period, the study site experienced a high variability in both precipitation and air temperature (covering approximately 88 % and 54 % of the variability observed in a 20-year period for precipitation and temperature, respectively) (Fig. 2).General model 1 parameterized with NDVI red-edge (adjR 2 = 0.74), and general model 3 (adjR 2 = 0.73) and 4 (adjR 2 = 0.79) based on the reflectance data were successful in capturing the inter-annual variability of GEP m among the 5 years characterized by different climatic conditions.Therefore, these results support the use of ground spectral measurements for monitoring GEP m in a long-term framework.We must however emphasize that the possible limitation of the approach based on VIs related to "canopy greenness" is that variations of GEP due to the short term environmental stresses cannot be monitored by these VIs, unless these stresses affect chlorophyll content (Gitelson et al., 2008).
Combining proximal sensing with EC observations may be relevant also for the EC data gap filling.In fact, the accuracy and reliability of the EC measurements depend on certain theoretical assumptions (e.g., requirement for turbulent and non-advective atmospheric conditions, stationarity of the measured fluxes) which often cannot be fulfilled in real field conditions (Foken et al., 2004;Göckede et al., 2004;Papale et al., 2006).The need for rejecting data acquired during periods when the abovementioned micrometeorological conditions were not met or due to other reasons such as nonoptimal wind directions, equipment failures etc. results in data set gaps constituting from 20 % to 60 % of annual data (Falge et al., 2001;Hui et al., 2004;Moffat et al., 2007).One of the most widely used gap-filling routines is based on the modeling of flux data with available environmental variables by means of nonlinear regression (Aubinet et al., 2000;Falge et al. 2001).This technique uses two equations -one for the response of ecosystem respiration (R eco ) to temperature and one for the light response of GEP (Moffat et al., 2007) allowing their predictions during gaps.The implementation of VIs into the light response model might help to improve the gap-filling results, especially in very dynamic ecosystems such as croplands, grasslands or deciduous forests.This could be particularly useful in case of long gaps in the EC data, which are inherently associated with a large degree of uncertainty (Moffat et al., 2007;Richardson and Hollinger, 2007;Wohlfahrt et al., 2010) and in case of managed ecosystems, where carbon dioxide uptake depends not only on the incoming radiation seasonality but also on cutting and grazing events.The results of a simple gap-filling approach presented in this study (based on creating and superimposing randomly distributed artificial gaps of three different lengths on the real data set and comparing GEP m values derived from EC with GEP m values filled with the best performing spectral models) encourage the use of spectral data in the gap-filling procedures of EC flux time series.The spectral based models were able to predict GEP m values with a performance comparable with others methods (Moffat et al., 2007) with adjR 2 ranging from 0.70 (5-day-long gap, general model 1 parameterized with NDVI red-edge ) to 0.78 (1-day-long gap, general model 4 based on reflectance at 681, 720 and 781 nm and PAR m data) (Table 6).

Conclusions
This study investigated the potential of a commercially available system -based on a 16-band multispectral sensor -for monitoring mean midday gross ecosystem production (GEP m ) in a dynamic subalpine grassland ecosystem of the Italian Alps.Chlorophyll-related indices including the red-edge part of the spectrum in their formulation (such as NDVI red-edge and CI red-edge ) were the best predictors of GEP m , and were able to explain most of its variability (adjR 2 = 0.74 for NDVI red-edge , adjR 2 = 0.73 for CI red-edge ) during the 5 consecutive years of observations, characterized by different climatic conditions.Our results confirm the findings of the literature regarding the complexity of the response of ecosystem productivity to change in PAR (Peng et al., 2013).This response is influenced by many variables.In fact, in our study, the accuracy of GEP m estimation decreased after including incident PAR m into the linear regression model, and the photosynthesis process was shown to be more efficient under diffuse compared to direct radiation.Further investigations are planned in order to explore the utility of including DI and PAR potential in the models to improve their performances.Also, the use of the reflectance approach instead of the VI approach did not lead to considerably improved results in estimating GEP m .Although a more detailed analysis of the full vegetation spectrum is desirable (for providing best performing algorithms and a solid basis for in situ R681)/(R750 + R681) Rouse et al. (1973) MSR (R750/R681 − 1)/(R750/R681 + 1) 1/2 Haboudane et al. (2004) DR (R750 − R720)/(R750 − R681) Datt (1999) NDVI red-edge (R750 − R720)/(R750 + R720) Gitelson and Merzlyak (1994) CI red-edge (R750/R720) − 1 Gitelson et al. (2003a)

Figure 2 .
Figure 2. Cumulative precipitation versus average daily air temperature for the period May-November.

Figure 3 .
Figure 3. Seasonal courses of normalized spectral vegetation indices -nVIs (-) and normalized mean midday gross ecosystem production -nGEP m (-) in the growing season of 2012; adjR 2 between GEP m estimated from EC measurements and GEP m obtained with model 1 fed with the various VIs.

Figure 4 .
Figure 4. Relationship between the red-edge normalized difference vegetation index (NDVI red-edge ) and mean midday gross ecosystem production (GEP m ), considering both the 5 years of observation together and annual observations.Dashed and solid trend lines refer to the general model 1 (considering each year of observation) and model 1 based on the observations from a specific year, respectively.

Figure 5 .
Figure 5.Light response of half-hourly gross ecosystem production (GEP; from 11:00 a.m. to 1:00 p.m.) to the incident photosynthetically active radiation (PAR) in the snow-free period of 2012 (May-November).Diffusion index (DI) is the ratio between diffuse and total incident PAR.It ranges from 0 to 1.

Figure 6 .
Figure 6.Root mean square error (RMSE) of the validated models based on the red-edge normalized difference vegetation index (NDVI red-edge ).

Table 2 .
Spectral vegetation indices presented in this study: normalized difference vegetation index (NDVI), modified simple ratio (MSR), difference ratio (DR), red-edge normalized difference vegetation index (NDVI red-edge ), and chlorophyll index (CI red-edge ).R refers to reflectance at a specific band (nm).

Table 4 .
Summary of the statistics (n -number of observations, adjR 2 -adjusted coefficient of determination, RMSE -root mean square error, PRMSE -percentage root mean square error) of the two linear regression models tested in this study both annually and considering all of the 5 observation years together.The three best-fitting models in each group are printed in bold.The best performing model is additionally highlighted in italic.All the regressions were statistically significant (p < 0.01).

Table 5 .
Summary of the general multiple regressions: partial adjusted R 2 , variance inflation factor (VIF), significance levels of the predictor variables (p), number of observations (n), cumulative adjusted R 2 , root mean square error (RMSE) and percentage root mean square error (PRMSE).R refers to reflectance at a given waveband (e.g., R720 -reflectance at 720 nm).

Table 6 .
Summary of the statistical metrics of gap-filling procedure: adjusted R 2 (adjR 2 ), root mean square error (RMSE) and percentage root mean square error (PRMSE).