Relationship between Hyperspectral Reflectance and Co 2 Exchange Printer-friendly Version Interactive Discussion on the Relationship between Ecosystem-scale Hyperspectral Reflectance and Co 2 Exchange in European Mountain Grasslands Bgd Relationship between Hyperspectral Reflectance and Co 2 Exchang

In this paper we explore the use of hyperspectral reflectance measurements and vegetation indices (VIs) derived therefrom in estimating carbon dioxide (CO 2) fluxes (net ecosystem exchange – NEE; gross primary production – GPP), and some key ecophys-iological variables related to NEE and GPP (light use efficiency – ε; initial quantum yield 5 – α; and GPP at saturating light – GPP max) for grasslands. Hyperspectral reflectance data (400–1000 nm), CO 2 fluxes and biophysical parameters were measured at three grassland sites located in European mountain regions. The relationships between CO 2 fluxes, ecophysiological variables and VIs derived using all two-band combinations of wavelengths available from the whole hyperspectral data space were analysed. We 10 found that hyperspectral VIs generally explained a large fraction of the variability in the investigated dependent variables and that they generally exhibited more skill in estimating midday and daily average GPP and NEE, as well as GPP max , than α and ε. Relationships between VIs and CO 2 fluxes and ecophysiological parameters were site-specific, likely due to differences in soils, vegetation parameters and environmen-15 tal conditions. Chlorophyll and water content related VIs (e.g. CI, NPCI, WI), reflecting seasonal changes in biophysical parameters controlling the photosynthesis process, explained the largest fraction of variability in most of the dependent variables. A limitation of the hyperspectral sensors is that their cost is still high and the use laborious. At the eddy covariance with a limited budget and without technical support, we suggest to 20 use at least dual or four channels low cost sensors in the the following spectral regions: 400–420 nm; 500–530 nm; 750–770 nm; 780–800 nm and 880–900 nm. In addition, our findings have major implications for up-scaling terrestrial CO 2 fluxes to larger regions and for remote and proximal sensing sampling and analysis strategies and call for more cross-site synthesis studies linking ground-based spectral reflectance with ecosystem-25 scale CO 2 fluxes.


Introduction
Understanding the mechanisms that drive the carbon dioxide (CO 2 ) exchange of terrestrial ecosystems is one of the main challenges for ecologists working on climate change (Beer et al., 2010).Plant gross photosynthesis, also referred to as gross primary productivity (GPP), is one of the major components of the global carbon cycle.It interacts in complex ways with environmental factors such as radiation, nutrients, soil moisture, vapour pressure deficit, air temperature and soil temperature (Drolet et al., 2005).Plant biochemistry and structure determine many fundamental ecosystem patterns, processes and dynamics (Lambers et al., 1998;Waring and Running, 1998).The canopy nitrogen content regulates the canopy photosynthetic capacity and the canopy light use efficiency (ε) (Ollinger et al., 2008).In addition, the canopy chlorophyll content plays an important role in controlling ecosystem photosynthesis and carbon gain (Peng et al., 2011;Gitelson et al., 2006).
Optical remote sensing can help ecologists in qualitatively and quantitatively assessing plant and canopy properties (e.g.biomass - Vescovo et al., 2012;water content -Clevers et al., 2010;nitrogen content -Ollinger et al., 2008, Knyazikhin Published by Copernicus Publications on behalf of the European Geosciences Union.et al., 2012;chlorophylls -Gitelson et al., 2006;and photosynthetic rate -Inoue et al., 2008) that drive ecosystem processes related to the carbon cycle.
Empirically and physically based methods have been proposed by several authors to interpret optical plant and canopy properties.Empirical methods consist of, for example, linear regression analysis between plant or canopy properties and optical data.The most widely used empirical methods are hyperspectral index methods (Peñuelas et al., 1993;Sims and Gamon, 2002;Inoue et al., 2008) and multivariable statistical methods, e.g.stepwise linear regression, genetic algorithm and neural network (Grossman et al., 1996;Riaño et al., 2005a;Li et al., 2007).Physical methods are based on the use of radiative transfer models (RTMs) to simulate light absorption and scattering through the canopy as a function of canopy structure and leaf biochemical composition (Jacquemoud et al., 2000;Zarco-Tejada et al., 2003).Therefore, RTMs help in quantifying the contribution of canopy biophysical and biochemical variables to canopy reflectance.Such models can be used for identifying regions of the light spectrum that are of particular importance for specific biophysical properties of vegetation.For example, it was demonstrated that the red-edge region (between 680 to 730 nm) of the spectrum is sensitive to the leaf chlorophyll content and leaf area index (LAI) (Baret et al., 1992).It is also well-accepted that an increase in LAI includes a decrease in reflectance in the red and an increase in the nearinfrared (NIR) region (Jacquemoud, 1993).In the NIR region, effects of LAI and the leaf angle distribution equally contribute to the reflectance response (Bacour et al., 2002a).NIR reflectance between 800 and 850 nm is also related to canopy N content (Ollinger et al., 2008;Knyazikhin et al., 2012).In addition, the combination of the reflectance in NIR and in the short-wave infrared region (SWIR) is correlated with canopy water content (Colombo et al., 2008), but the reflectance between 1000 and 1400 nm is also highly sensitive to LAI.Therefore, some attention is needed when these spectral regions are used to retrieve water content, considering that the canopy properties in a given ecosystem often covary (Bacour et al., 2002b).
The drawback of such an approach consists in the fact that the process of building a model implies approximations and assumptions.For this reason we opted for a purely databased approach such as the hyperspectral index approach.This method consists of the use of spectral vegetation indices (VIs) defined as spectral band ratios or normalized band ratios between the reflectance in the visible (VIS)-vs.-NIR,VIS-vs.-VIS or NIR-vs.-NIRregion.
The typical optical sampling approach, which consists of linking spectral observations with CO 2 fluxes, is based on the Monteith (1972Monteith ( , 1977) ) equation: where ε is the light use efficiency and f APAR is the fraction of absorbed photosynthetically active radiation); both ε and f APAR can be retrieved by remote optical observations.A large number of VIs that can potentially be used to model the productivity of terrestrial ecosystems (as a proxy of ε and f APAR ) have been suggested (Inoue et al., 2008;Coops et al., 2010;Peñuelas et al., 2011;Rossini et al., 2012).The various VIs differ in their sensitivity to changes in photosynthetic status."Greenness indices" -such the widely used normalized difference vegetation index (NDVI) -have been demonstrated to be a good proxy for f APAR but are not sensitive to rapid changes in plant photosynthesis which are induced by common environmental and anthropogenic stressors (Gitelson et al., 2008;Hmimina et al., 2014;Soudani et al., 2014).However, in ecosystems characterized by strong dynamics (e.g.grasslands and crops with a strong greenup and senescence), other VIs are able to effectively monitor seasonal changes in biophysical parameters controlling canopy photosynthesis, such as f APAR and chlorophyll content, and, consequently, can be adopted to monitor the seasonal and spatial variability of carbon fluxes (Gitelson et al., 2012;Sakowska et al., 2014).Short-term changes in ε can be remotely detected through a spectral proxy of the xanthophyll cycle (photochemical reflectance index, PRI; Gamon et al., 1992).The PRI is one of the most promising VIs for a direct estimation of photosynthetic light use efficiency and of its seasonal and diurnal variations (Nichol et al., 2002).
Latest developments in the sun-induced fluorescence method may allow even more direct remote sensing of plant photosynthesis in the near future (Meroni et al., 2009;Rossini et al., 2010;Frankenberg et al., 2011).On the canopy scale, the relationship between PRI and ε was shown to be sitedependent (Garbulsky et al., 2011;Goerner et al., 2011) and strongly affected by environmental conditions (Soudani et al., 2014).
Whereas previous studies have demonstrated the ability of remote-sensing data to allow modelling ecosystem GPP on the ecosystem scale (e.g.Gianelle et al., 2009;Wohlfahrt et al., 2010;Rossini et al., 2012;Sakowska et al., 2014), a universal model for GPP estimation applicable across different ecosystems and a wide range of environmental conditions is still missing.Previous studies often focussed on single sites with specific characteristics (e.g.climate, vegetation composition, soil type; see Wohlfahrt et al., 2010).In addition, different studies often used different sensors, platforms and protocols (Balzarolo et al., 2011), making generalization difficult.Moreover, most of the studies have either relied on reflectance measurements in a few spectral wavebands (e.g.Wohlfahrt et al., 2010;Sakowska et al., 2014), missing potentially important information in under-sampled spectral regions.In order to overcome such heterogeneity in spectrometry measurements, SpecNet (http://specnet.info;Gamon et al., 2006), the European COST Action ES0903 (COST -Cooperation in Science and Technology; EUROSPEC; http: //cost-es0903.fem-environment.eu/)and the COST Action ES1309 (OPTIMISE; http://www.cost.eu/domains_actions/essem/Actions/ES1309) focused on the definition of a stan-  -170, 2006 (9) 122-303, 2006 (16) 129-201, 2005 (13)  124-192, 2006 (12)    * Range refers to days of the year (DOY); this is followed by the year; the number of hyperspectral measurement dates is given in brackets.
dardized protocol for making optical measurements at the eddy covariance CO 2 flux towers (Gamon et al., 2010).
The overarching objective of the present paper is thus to develop a common framework for predicting grassland carbon fluxes and ecophysiological parameters based on optical remote-sensing data across measurement sites exposed to diverse natural (climate) and anthropogenic (management) factors.To this end we combine eddy covariance CO 2 flux measurements with ground-based hyperspectral reflectance measurements for six different grasslands in Europe.In order to make the optical and fluxes measurements comparable, these were acquired at the six sites following a common protocol resulting in a unique standardized data set.We focused on European grasslands, which cover roughly 22 % (80 million ha) of the EU-25 land area and are thus among the dominating ecosystem types in Europe (EEA, 2005).Accordingly, their role in the European carbon balance has received a lot of scientific interest (Soussana et al., 2007;Gilmanov et al., 2007;Wohlfahrt et al., 2008;Ciais et al., 2010).Direct measurements of the grassland carbon exchange have been carried out and are still ongoing at a number of different grassland sites in Europe (e.g.Soussana et al., 2007;Cernusca et al., 2008).Scaling up these plot-level measurements to the continental scale requires a modelling approach typically based on or supported by remotely sensed data.Therefore, we believe that this study will improve the current knowledge on modelling the carbon dynamics of European grasslands.

Experimental site description
This study was carried out at six experimental mountain grassland sites in Europe covering different climatic and grassland management conditions existing in the mountain regions of Europe, which were already part of the preceding study by Vescovo et al. (2012).This data set combined in situ hyperspectral, biophysical and flux measurements based on common protocol (for more details, see Sects. 2.2,2.3 and 2.4).This data set is unique since no common protocol for hyperspectral measurements exists in the various eddy covariance networks (e.g.FLUXNET).In this study, three of these sites (Amplero, Neustift and Monte Bondone, see Tables 1 and S1 in the Supplement) composed the main data set used in the analysis, while the other three sites (Table S2) were used to independently validate the models obtained with the main data set. 1 and S1):

Amplero
The Amplero site is situated in the Mediterranean Apennine mountain region of Italy (41.90409 • N, 13.60516 • E) at 884 m a.s.l.This site is characterized by mild, rainy winters and by an intense drought in summer.Amplero is managed as a hay meadow, with one cut in late June and extensive grazing during summer and autumn (Balzarolo, 2008).

Monte Bondone
The Monte Bondone site is situated in the Italian Alps (46.01468 • N, 11.04583 • E) at 1550 m a.s.l.This site is characterized by a typical subcontinental climate, with mild summers and precipitation peaks in spring and autumn.Monte Bondone is managed as an extensive meadow, with one cut in mid-July.

Neustift
The Neustift grassland site is located in the Austrian Alps (47.11620 • N, 11.32034 • E) at 970 m a.s.l.The climate of this area is continental and Alpine, with precipitation peaks during the summer (July).This site is intensively managed as a hay meadow, with three cuts in mid-June, at the beginning of August and at the end of September.

Längenfeld
The site Längenfeld is located in the Austrian Alps (47.0612 • N, 10.9634 • E) at 1180 m a.s.l.The climate of this area is continental/Alpine; however, compared to the other Alpine sites in this study, the site receives comparably less precipitation due to rain shadowing effects from both the north and south.The site is intensively managed as a hay meadow, with three cuts in mid-June, mid-August and mid-October.

Leutasch
The site Leutasch is located in the Austrian Alps (47.3780 • N, 11.1627 • E) at 1115 m a.s.l.The climate of this area is Alpine, with substantial precipitation due to its position on the north range of the Alps.The site is extensively managed as a hay meadow, with two cuts at the end of June and beginning of September.

Scharnitz
The site Scharnitz is located very close to Leutasch (47.3873 • N, 11.2479 • E) at 964 m a.s.l., and the climate is thus very similar to Leutasch.The site is extensively managed as a hay meadow, with two cuts at the beginning of July and beginning of September.

Hyperspectral reflectance measurements
The canopy hyperspectral reflectance measurements were collected at each site under clear-sky conditions close to solar noon (between 11:00 to 14:00 Central European Time) using the same model of a portable spectroradiometer (ASD Field-Spec HandHeld, Inc., Boulder, CO, USA) at all sites.The spectroradiometer acquires reflectance values between 350 and 1075 nm, with a full-width half maximum (FWHM) of 3.5 nm and a spectral resolution of 1 nm.In order to achieve a better match between the eddy covariance flux footprint and optical measurements, a cosine diffuser foreoptic (ASD Remote Cosine Receptor, Inc., Boulder, CO, USA), calibrated by the manufacturer, was used for nadir and zenith measurements (Gianelle et al., 2009;Fava et al., 2009;Meroni et al., 2011).The ASD's cosine receptor is designed with a geometry and material that provides a hemispherical field of view (FOV) of 180 • and optimizes the cosine response.To reduce the nadir FOV contamination (i.e.sky irradiance and canopy irradiance) due to the hemispherical view of the sensor the instrument was placed on a 1.5 m long horizontal arm at a height of 1.5 m above the ground.To avoid the zenithal FOV contamination, the measurements were taken at least at a 15 m distance from the eddy covariance tower (maximum height of the tower was 6 m).The vegetation irradiance (sensor-pointing nadir) and sky irradiance (sensor-pointing zenith) were measured by rotating the spectroradiometer alternately to acquire spectra from the vegetation and from the sky.Hemispherical reflectance was derived as the ratio of reflected to incident radiance.Each reflectance spectrum was automatically calculated and stored by the spectroradiometer as an average of 20 readings.Before starting each spectral sampling, a dark current measurement was done.For more details on the experimental set-up, see Vescovo et al. (2012).Spectral measurements were collected from spring until the cutting date at Amplero and Monte Bondone, while at the site in Neustift, which is cut three times during the season, spectral measurements were taken about once per week throughout the growing season of 2006.

Biophysical and biochemical canopy properties
Samples for dry phytomass, nitrogen and water content measurements were collected at the time of the hyperspectral measurements in the field of view of the hyperspectral sensor (see Vescovo et al., 2012 for more details).A similar data set was collected in 2013 at Monte Bondone by combining hyperspectral data with chlorophyll measurements.Chlorophyll samples were collected in the field of view of the hyperspectral sensor and chlorophyll content was detected by UV-VI spectroscopy.First, the samples were ground with liquid nitrogen and then immersed in 80 % acetone solution (0.1 g per 10 mL), shaken for 10 min in an automatic shaker at 250 rpm (Universal Table Shaker 709) and centrifuged at 4000 rpm for 10 min (Eppendorf 5810 R) in order to remove particles from the solution.The absorbance of extracted solutions was measured at 470, 646.8 and 663.2 nm by a UV-VIS spectrophotometer (Shimadzu UV-1601), and the concentrations of chlorophyll a (C a ), chlorophyll b (C b ) and carotenoids (C x+c ) were calculated as proposed by Lichtenthaler (1987).The weight of sampled sediment was used to calculate pigment concentrations per unit leaf mass (mg g −1 ), and the weight of green biomass per ground area was used to obtain the total chlorophyll content (mg m −2 ).

CO 2 flux measurements
Continuous measurements of the net ecosystem CO 2 exchange (NEE) were made by the eddy covariance (EC) technique (Baldocchi et al., 1996;Aubinet et at., 2012) at the six study sites using identical instrumentation.The three wind components and the speed of sound were measured using ultrasonic anemometers, and CO 2 molar densities were measured using open-path infrared gas analysers (IRGAs), as detailed in Tables S1 and S2.Raw data were acquired at 20 Hz and averaged over 30 min time windows in postprocessing.Turbulent fluxes were obtained from raw data by applying block averaging (Monte Bondone, Neustift, validation sites) or linear de-trending (Amplero) methods with a time window of 30 min.A 3-D coordinate correction was performed according to Wilczak et al. (2001).The CO 2 fluxes were corrected for, the effect of air density fluctuations as proposed by Webb et al. (1980).Low-and high-pass filtering was corrected for following Aubinet et al. (2000) (Amplero, Monte Bondone) or Moore (1986) (Neustift, validation sites).Data gaps due to sensors malfunctioning or violation of the assumptions underlying the EC method were removed and filled using the gap-filling and flux-partitioning techniques as proposed in Wohlfahrt et al. (2008).Ecosystem respiration (R eco ) was calculated from the y intercept of the light response model (see Eq. 4).Gross primary productivity (GPP) was calculated as the difference between NEE and R eco .Half-hourly NEE and GPP values were averaged between 11:00 to 14:00 solar local time (at the time window of optical measurements) to allow for direct comparison with the hyperspectral data, and daily sums were also computed.At each site the following supporting environmental measurements were acquired: photosynthetically active radiation (PAR; quantum sensors), air temperature (T a ; PT100, thermistor and thermoelement sensors) and humidity (RH; capacitance sensors) at some reference height above the canopy and soil temperature (T s ; PT100, thermistor and thermoelement sensors) and volumetric water content (SWC; dielectric and time-domain reflectometry sensors) in the main rooting zone.In this study we used CO 2 flux and meteorological data of the years 2005 and 2006 for Monte Bondone and of 2006 for the other sites.

Estimation of grassland ecophysiological parameters
Canopy light use efficiency (ε) was derived from photosynthetically active radiation (PAR) absorbed by the canopy (APAR) as and was estimated both at midday and daily time resolutions.We estimated the fraction of PAR absorbed by the canopy (f APAR ) from measured values of the leaf area index (LAI) using the Lambert-Beer law: where k is the canopy extinction coefficient (fixed at k = 0.4 as defined for a southern mixed-grass prairie in Texas; Kiniry et al., 2007) and 0.95 is the proportion of intercepted PAR that is absorbed by plants (Schwalm et al., 2006).LAI was quantified non-destructively by an indirect method based on canopy PAR transmission using line PAR sensors (SunScan, Delta-T, UK) and an inversion of an RTM (Wohlfahrt et al., 2001).These measurements were done within the footprint area of the spectroradiometer simultaneously with the hyperspectral measurements.Three additional key parameters of the response of NEE to PAR were extracted by fitting measured NEE and PAR to a simple Michaelis-Menten-type model: where α represents the apparent quantum yield (µmol CO 2 µmol −1 photons), F sat the asymptotic value of GPP (µmol CO 2 m −2 s −1 ), PAR the photosynthetically active radiation (µmol photons m −2 s −1 ) and R eco the ecosystem respiration (µmol CO 2 m −2 s −1 ).For all sites, using the Levenberg-Marquardt (Marquardt, 1963) algorithm the parameters of Eq. ( 4) were estimated by fitting Eq. ( 4) to both day and nighttime data, which were pooled into 3-day blocks centred on the date of the hyperspectral data acquisition.
For each acquisition date, we then used Eq. ( 3) to derive GPP at an incident PAR of 1500 µmol m −2 s −1 , referred to as GPP max in the following.

Hyperspectral data analysis
In order to explore the information content of the hyperspectral data for estimating CO 2 fluxes (i.e.midday and daily average of NEE and GPP) and ecophysiological parameters (i.e.α, ε and GPP max ), we performed a correlation analysis between spectral reflectance indices (independent variables) and these (dependent) variables.To this end, we derived spectral ratio (SR; Eq. 5), spectral difference (SD; Eq. 6) and normalized spectral difference (NSD; Eq. 6) indices using all possible two-band (i and j ) reflectance (ρ) combinations between 400 and 1000 nm (180 600 combinations).These three formulations were selected since they represent the most common equations used to compute vegetation indices (see Table 2).
Linear regression analysis was performed among all possible wavelength combinations for all three index types (SR, NSD and SD) and the investigated dependent variables.The performance of linear models in predicting dependent variables (i.e.carbon fluxes and ecophysiological parameters) was evaluated by the coefficient of determination (R 2 ) and root mean square error (RMSE).The coefficients of determination (R 2 ) resulting from the linear models were visualized in correlograms, as in the example in Fig. 1.
We also calculated four SR and seven NSD indices which are commonly used in relation to vegetation activity and CO 2 fluxes (Table 2).Figure 1 shows the location of these indices in the waveband space of the correlograms.In this analysis, we also considered the enhanced vegetation index (EVI), which is one of the most frequently used vegetation indices to predict CO 2 fluxes.In Fig. 1 the location of the EVI is not shown since this index is computed by the combination of three spectral bands as shown in Table 2.
The robustness of the model selected on the basis of the best band combinations for all ecophysiological parameters for each site and all sites pooled was tested by the leave-oneout cross-validation technique.The predictive performance was expressed as the cross-validated coefficient of determination (R 2 CV ) and the cross-validated root mean square error (RMSE CV ).In addition, the capability of the selected models in predicting different ecophysiological parameters was tested by applying the selected models to the validation data set (Table S2) composed of three different grasslands not used in the previous analysis.This data set was selected because the hyperspectral and flux data were collected by using exactly the same protocol applied for the main data set (see Sect. 2.1).
In order to explore the basis of the correlation between the selected band combinations and ecophysiological variables (e.g.α, GPP max , GPP, ε), the relationship between the selected bands and biophysical parameters such as dry phytomass, nitrogen and water content collected during the field campaign in the same footprint of the hyperspectral measurements was examined.

Band selection based on the combination of random forests and genetic algorithm (GA-rF)
In order to complement the more conventional analysis described in the previous section, we also explored the use of a hybrid feature selection strategy based on a genetic algorithm and random forests (GA-rF).The first method was used for the feature selection and the second one as regression for predicting the target variables.First of all, the original data set was aggregated to 10 nm bands in order to reduce the effects of autocorrelation in frequency space.The algorithm generates a number of possible model solutions (chromosomes) and uses these to evolve towards an approximation of the best solution of the model.In our case the genes of each chromosome correspond to the wavebands.We made use of five genes for each chromosome in order to overcome overfitting.Each population of 1000 chromosomes evolved for 200 generations.The mutation chance was set to the inverse of a population size increase of 1.The fitness of each chromosome was measured by applying the random forest algorithm (Breiman, 2001).This was used as an ensemble method for regression that is based on the uncontrolled development of decision trees (n = 100).We opted for this method because of its demonstrated efficiency with large data sets.In combining the two methods, we choose the mean squared error as the target variable to be minimized.

Seasonal variation of meteorological variables, LAI and CO 2 fluxes
Environmental conditions and the seasonal development of LAI, NEE, GPP α, ε and GPP max during the study period are shown in Fig. 2. A strong influence of the typical climatic conditions at the three study sites was evident: Amplero was characterized by a Mediterranean climate, with highest incoming radiation and temperatures and the lowest amount of precipitation, which translated into a substantial seasonal drawdown of soil moisture.Monte Bondone and Neustift, more influenced by a continental Alpine climate, experienced comparably lower temperatures with higher precipitation and soil moisture with respect to Amplero (Fig. 2).Maximum LAI values were similar at Monte Bondone and Amplero (2.8-3.4 m 2 m −2 ), while twice as much leaf area developed at the more intensively managed study site Neustift.The latter was also characterized by higher NEE and GPP (i.e. more photosynthesis and net uptake of CO 2 ).The reductions in leaf area related to the cuts of the grasslands were associated, as expected, with marked increases and reductions in NEE and GPP, respectively.The canopy light use efficiency, ε, was inversely related to GPP and LAI, peaking at the beginning of the season at Amplero and Monte Bondone (0.01-0.10 µmol photons µmol −1 CO 2 ),  Vegetation index improved for soil and atmospheric effects Huete et al. (1997) while for Neustift ε showed the highest values after the cuts (0.01-0.20 µmol photons µmol −1 CO 2 ).At Amplero, α and GPP max peaked in spring and then decreased during the summer drought period, while at Neustift and Monte Bondone, temporal patterns of α and GPP max were more strongly affected by management.A number of interesting insights may be gained from Figs. 3-5 and S2-S4, which we summarize in the following:

Hyperspectral data and their relation to CO 2 fluxes and ecophysiological parameters
1.The correlograms exhibited quite different patternssome correlograms showed that a wide range of band combinations was able to explain the simulated quanti-ties (e.g.GPP at Amplero; Figs. 3 and S2), while some correlograms exhibited very pronounced patterns, with the R 2 value changing greatly with subtle changes in band combinations (e.g.ε at Neustift; Figs. 3 and S2).
2. Maximum R 2 values were often clearly higher than the surrounding areas of high predictive power (e.g.ε at Amplero; Fig. 3).
3. The different types of indices (compare Figs. 3-5) yielded similarly high correlations with the same dependent variable at the same site in similar spectral regions.This indicates that band selection is more important for explanatory power than the mathematical formulation of the VI (i.e.ratio vs. difference, with/without normalization).SR and NSD indices (Figs. 3 and 4) yielded similar results compared to SD indices (Fig. 5).    the visible and NIR (Fig. S1), resulted in comparably lower correlations.

The highest correlations
5. For midday and daily time resolutions different band combinations were selected (e.g.NEE at Amplero; compare Figs. 4 and S3).For similar selected regions, daily averages were characterized by higher explanatory power compared to midday averages (e.g.ε at Neustift).
Figure 6 shows the performance of linear regression models for the best NSD-type indices for midday ecophysiological parameters for each site and all sites pooled (Fig. S7 shows the results of the same analysis for daily averages).Large differences existed between the study sites regarding the explanatory power of the same index for the same dependent variable.The highest R 2 cv and values were generally obtained for Amplero, followed by Neustift and then Monte Bondone, and the lowest R 2 cv values resulted when data from all three sites were pooled, confirming the difficulties in finding a general relation valid among sites.
For Amplero and Neustift the NIR-vs.-NIRcombinations showed a positive linear regression model with α, GPP max and GPP, while for Monte Bondone a negative linear correlation was observed.For Amplero the VIS-vs.-VIScombination showed a good performance in predicting ε; the NIR-vs.-NIRcombinations showed good performance for Neustift, as did the VIS-vs.-NIRcombination for Monte Bondone.The linear models for NEE were site-specific.In fact, Amplero and Monte Bondone showed a positive linear regres-

Correlation between conventional VIs, ecophysiological variables and CO 2 fluxes
The correlation analysis between the conventional VIs, the midday CO 2 fluxes (Table 3) and ecophysiological parameters (Table 4), generally confirmed the results obtained with the hyperspectral data.
For the same dependent variable (α, GPP max , GPP, ε and NEE), the performance of the various VIs showed large differences between sites.For example, for GPP max all of the investigated indices except NPQI resulted in significant linear correlations at Amplero, explaining 41-89 % of the variability in GPP max .In contrast, only NDVI, PRI, NPCI and SRPI showed a slightly significant linear performance (17-26 %) for GPP max at Neustift.
The different VIs performed differently in predicting the same dependent variable at the different study sites.For all dependent variables (Tables 3, 4 and S3), the VI resulting in the highest R 2 values was never the same at all sites.Often the best-fitting VI at one site resulted in a non-significant correlation at another site.Therefore, none of the dependent variables clearly emerged as the one best predicted (Tables 3,  4 and S3).
When data from all sites were pooled, models showed the same performance for the same VI and dependent variable except for GPP and NEE.The best-performing VI for GPP and NEE was SIPI; NPCI performed best for α, GRI for ε and SIPI for GPP max .
The choice of the averaging period (midday vs. daily) applied to ε, NEE and GPP generally did not modify the ranking of the VIs, but the R 2 values tended to be similar or somewhat higher on the daily timescale (compare Tables 3 and 4 with Table S3).Overall, the results of the validation were mixed.Good performance was observed mainly for the Neustift, Monte Bondone and pooled M.Bondone&Neustift models (see Table S4).In particular, the best performance values were obtained for (i) α for the pooled Monte Bondone and Neustift model for SR-type indices; (ii) GPP max for NSD-, SR-and SD-type indices and models except for Amplero; and (iii) for midday GPP for all NSD-, SR-and SD-type indices and models.It is interesting to note that lower performances were generally found for the models based on the Amplero parameterization.This is understandable as Monte Bondone and in particular Neustift were structurally and functionally much more similar to the validation sites compared to Amplero (Tables 1 and S2).Considerably poorer performance was observed for ε and NEE across all model-index type combinations (Table S4).The validation on a daily timescale always resulted in a poorer performance compared to the midday average timescale (Fig. S10 and Table S5).

Effects of canopy structure on selected band combinations
Tables 5 and S6 show the results of the correlation analysis between the selected NSD-, SR-and SD-type indices for ecophysiological variables and fluxes and biophysical properties of vegetation, such as dry phytomass, nitrogen and water content.Overall, the spectral response in the selected band combinations for NSD-, SR-and SD-type indices was strongly related to vegetation properties of the three grasslands (e.g.nitrogen and dry phytomass), which impacted on their spectral response in the NIR and VIS regions.For the Mediterranean site (Amplero) and for all ecophysiological parameters (e.g.α, GPP max ), dry phytomass was the main driving factor of the spectral response in the selected bands, while nitrogen content drove the spectral response in the NIR region for Neustift.For Monte Bondone, both dry phytomass and nitrogen content affected the spectral response of the grassland.Similar results were obtained for SR-and SD-type indices.
Figures 8 and S11 show the correlation analysis between the selected NSD-, SR-and SD-type indices for all ecophysiological parameters (e.g.α, GPPmax) and chlorophyll content for Monte Bondone in 2013.The chlorophyll content showed a very good correlation for all selected models and for all indices.The values of R 2 were always higher than the values of R 2 obtained for the other biophysical variables (Tables 5 and S6).In Fig. 9, it is possible to see that NSD-and SR-type indices for the selected bands for estimating GPP (i.e.996 and 710 nm) are strongly correlated with canopytotal chlorophyll content (R 2 > 0.80).

Band selection using the GA-rF method
Figure 9 shows the results of the band selection based on the GA-rF method.In particular, each plot represents the frequency of the occurrence of each band in the genetic algorithm.
Overall, using the GA-rF method it was possible to identify portions of the spectrum that were of particular significance for estimating specific properties of the different ecosystems.For example, for predicting midday GPP (Fig. 9b) for all sites pooled together, the bands at 430, 630, 660 and 710 nm showed the best results.The bands at 505 nm, 660 nm and 710 nm played an important role in predicting midday GPP for Amplero, Neustift and Monte Bondone, respectively.Some differences were found for the different time resolutions (compare Fig. 9b and c).For example, the bands at 580 and 800 nm showed the best results for Amplero, and bands at 530 nm showed the best results for Neustift.Figure 10 shows the results for the band selection by GA-rF methods for biophysical variables (i.e.dry phytomass, nitrogen and water content).For the variables related to slow processes, the GA-rF method highlighted different bands for different sites; a much higher between-site variability for the variables related to ecophysiological processes (e.g.ε, α and GPP max ) was detected, and we weren't able to identify common "hot spots".

Discussion
This study aimed at evaluating the potential of hyperspectral reflectance measurements to simulate CO 2 fluxes and ecophysiological variables of European mountain grasslands over a range of climatic conditions and management practices (grazing, harvest).To this end, we combined eddy covariance CO 2 flux measurements with ground-based hyper-spectral measurements at six mountain grassland sites in Europe.

Upscaling of in situ relationships between VI indices and CO 2 fluxes and ecophysiological parameters
Despite the fact that we focused on a single type of ecosystem, our results showed that large differences existed among the investigated sites in the relationships between hyperspectral reflectance data and CO 2 fluxes and ecophysiological parameters.For all study sites pooled, hyperspectral reflectance data explained 40-68 % of the variability in the dependent variables (Figs.3-5).The conventional VIs yielded a maximum of 47 % of explained variability in the data (Tables 3-4).This is the first study comparing different grasslands characterized by different plant species and environmental con- ditions.The use of simple models based on a linear relationship between GPP and VIs, related to canopy greenness, has proven to be a good proxy for the GPP of ecosystems with strong green-up and senescence (Peng et al., 2011;Rossini et al., 2012).The loss of this relationship may be related to low ε variability due to abiotic and biotic stressors, the dependency of PRI on LAI, leaf and canopy biochemical structure (e.g.leaf orientation), and xanthophyll cycle inhibition or saturation, and zeaxanthin-independent quenching (Gamon et al., 2001;Filella et al., 2004;Rahimzadeh-Bajgiran et al., 2012;Hmimina et al., 2014).For alpine grasslands, a key meteorological variable that played a relevant role in stimulating ε was high soil water content associated with low temperatures (Polley et al., 2011).Low soil water contents also triggered a decrease in leaf conductance as well as in ε and in α for two oak and beech ecosystems (Hmimina et al., 2014).However, no significant differences in leaf biochemical and structural properties of the canopy at the lowest and highest water content were found.In addition, in this special issue, Sakowska et al. (2014) showed that ε is also strongly affected by the directional distribution of incident PAR, i.e. the ratio of direct to diffuse PAR.
Considering all sites pooled together (Figs. 3 and S2), NSD-type indices showed a very poor correlation in the VISvs.-NIRband combinations (i.e.traditional greenness indices; see Table 2) with GPP.It is well-known in the literature (Rossini et al., 2010(Rossini et al., , 2012;;Peng et al., 2010;Sakowska et al., 2014) that greenness indices, for grasslands and crops, are often good proxies of f APARgreen (and thus carbon fluxes).Interestingly, in our study their performance was considerably poorer than expected.The NSD-type index showed a better performance in VIS-vs.-VISband combinations than VISvs.-NIRones.VIS-vs.-VISband combinations for NSD-type indices (e.g.green vs. blue or red and green vs. green wavelengths; see, e.g., Inoue et al., 2008) are defined as greenness indices (Fig. 1), although their performance is generally much poorer than NSD VIS-vs.-NIRindices.These results are likely due to the confounding effects of the different canopy structures and consequently of the different NIR response of the investigated grasslands (see Fig. S1).In fact, the different grassland structures (spatial distribution of photosynthetic and also non-photosynthetic material, leaf angles, etc.) affect our ability to use traditional indices to estimate f APARgreen (and fluxes) when we consider different grasslands together because the structural effects on scattering are very complex in the NIR (Jacquemoud et al., 2009;Knyazikhin et al., 2012).These results are of importance for the community, which still relies on these relationships a lot, also favoured by the availability of affordable narrow-band sensors that allow continuous monitoring of, e.g., NDVI.These results suggest that waveband combinations not exploited by presently used (conventional) VIs may offer considerable potential for predicting grassland CO 2 fluxes; this has implications for the design and capabilities of future space, airborne or ground-based low-cost sensors.In particular, these results also have a strong impact on our ability to upscale grassland f APARgreen and carbon fluxes using future sensors (e.g.Sentinel 2).
The evaluation of the models for the main data set against three new sites showed that at least some of these models can be transferred to predict carbon fluxes and ecophysiological parameters for similar grasslands (Fig. 7).However, for some parameters (e.g.ε), the independent validation indicated a poor performance; it challenges the current practice in upscaling to larger regions by grouping all grasslands into a single plant functional type (PFT).We advocate more studies to be conducted merging CO 2 flux with hyperspectral data by means of models which use a more process-oriented and coupled approach to simulating canopy CO 2 exchange and reflectance in order to explore the causes underlying the observed differences between seemingly closely related study sites.

Grassland structural characteristics and their spectral response
Although we considered similar ecosystems (belonging to the same vegetation type), the investigated canopies were very different and included Mediterranean, extensive alpine NSD-type Amplero Dry phytomass (g m −2 ) 0.66 For Amplero and Neustift, NSD-type indices performed well for NIR-vs.-NIRband combinations for all investigated parameters, while Monte Bondone showed best performances in the VIS-vs.-NIRband combinations for GPP max and ε (Fig. 6).The dry phytomass was the main driving factor in the spectral response in NIR-vs.-NIRband combinations for Amplero, while nitrogen content drove the spectral response in NIR-vs.-NIRband combinations of Neustift for all parameters except for α (Tables 5 and S6).Interestingly, for Monte Bondone both dry phytomass and nitrogen content explained the spectral response of the grassland in VIS-vs.-NIRband combinations for GPP max and ε, while no significant relationships with biophysical variables were found for α, GPP or NEE.These results partially confirm the findings of Vescovo et al. (2012), who highlighted a strong relationship, for several grassland types, between an NSD-type index and phytomass.
For Monte Bondone, NSD-and SR-type indices for the selected bands for estimating all variables except α were strongly correlated with canopy-total chlorophyll content (R 2 > 0.85).
The chlorophyll indices (e.g.red-edge NDVI and CI; see Tables 3 and 4) -which are considered the best indices for estimating carbon fluxes in grasslands and crops -showed a good performance for Amplero and Monte Bondone in our data set but performed poorly for Neustift.
It was demonstrated by many authors that the red-edge domain, where reflectance changes from very low in the absorption region to high in the NIR, is one of the best descriptors of chlorophyll concentration.On the other hand, it is wellknown that the canopy structure can be a very strong confounding factor.Our results confirm that this topic needs to be further investigated, as this finding has a relevant impact concerning the use of Sentinel 2 to upscale f APAR and carbon flux observations.
It is interesting to see that the NSD-type indices in the NIR-vs.-NIRband combinations appeared to be the best proxy for GPP fluxes when all the grasslands were pooled together.These results can be linked to the controversial paper focused on the strong impact of structure on the ability to estimate canopy nitrogen content (Knyazikhin et al., 2012) and confirm the need for more studies in this direction.Good relationships were found between the NIR-vs.-NIRband combinations (> 750 nm wavelengths) and fluxes; the physical basis of these relationships needs to be further investigated.In fact, it is important to highlight that the literature indicates that the wavelengths in the NIR (> 750 nm) are not sensitive to chlorophyll content, but they are related to leaf, canopy structure and -around the 970 nm area -to water.
As confirmed by comparing the correlation matrix approach with the GA-rF approach, we could not find a universal relationship between reflectance in specific wavelengths of the light spectrum and biophysical properties of vegetation.We think that this is strongly linked to vegetation structure effects.For this reason we believe that further research is necessary to disentangle the impact of factors such as bidirectional reflectance distribution function and scaling effects.

Conclusions
The present study focused on understanding the potential of hyperspectral VIs in predicting grassland CO 2 exchange and ecophysiological parameters (α, ε and GPP max ) for different European mountain grasslands.
The major finding of this study is that the relationship between ground-based hyperspectral reflectance and the ecosystem-scale CO 2 exchange of mountain grasslands is much more variable than what might be supposed given this closely related group of structurally and functionally similar ecosystems.As a consequence, the unique models of mountain grassland CO 2 exchange, i.e. the best-fitting models for all sites pooled, explained 47 and 68 % of the variability in the independent variables when established VIs and optimized hyperspectral VIs, respectively, were used.Interestingly, VIs based on reflectance either in the visible or NIR part of the electromagnetic spectrum were superior in predicting mountain grassland CO 2 exchange and ecophysiological parameters compared to commonly used VIs, which are based on a combination of these two wavebands.The band selection based on the GA-rF algorithm confirmed that it is difficult to define a universal band range able to describe ecophysiological parameters, carbon fluxes and biophysical variables, even for a closely related group of ecosystems.
The take-home message from this study thus is that continuing efforts are required to better understand differences in the relationship between ecosystem-scale reflectance and CO 2 exchange and to improve models of this relationship which can be employed to upscale the land-CO 2 exchange to larger spatial scales based on optical remote-sensing data.Initiatives such as SpecNet (http://specnet.info;Gamon et al., 2006), the COST Action ES0903 (EUROSPEC; http: //cost-es0903.fem-environment.eu/)and the COST Action ES1309 (OPTIMISE; http://www.cost.eu/domains_actions/essem/Actions/ES1309) are instrumental in this as they provide the scale-consistent combination of hyperspectral reflectance and CO 2 exchange data.
The Supplement related to this article is available online at doi:10.5194/bg-12-3089-2015-supplement.

Figure 1 .
Figure 1.A selected example of a correlogram between NSD-type indices and midday average GPP for all sites pooled.The correlogram shows all R 2 values, the white dot indicates the two-band combination with the highest R 2 value and the red dots indicate the location of the reference VIs reported in Table 2 (SR: simple ratio; GRI: green ratio index; WI: water index; SRPI: simple ratio pigment index; NDVI: normalized difference vegetation index; NPQI: normalized phaeophytinization index; NPCI: normalized pigment chlorophyll index; CI: chlorophyll index; red-edge NDVI; SIPI: structural independent pigment index; PRI: photochemical reflectance index).

Figures 3 -
Figures 3-5 show correlograms between NSD-, SR-and SD-type indices, respectively, and the investigated dependent midday ecophysiological parameters and fluxes.The correlograms for daily data can be found in the Supplement (Figs.S2-S4 in Supplement).Selected examples of key spectral signatures of the investigated grasslands are shown in Fig. S1 in the Supplement.A number of interesting insights may be gained from Figs. 3-5 and S2-S4, which we summarize in the following:

Figure 3 .
Figure 3. Correlograms of R 2 values for α, GPP max , midday averaged GPP, ε and NEE, on the one hand, and NSD-type indices for Amplero, Neustift, Monte Bondone (both study years pooled) and all sites pooled on the other.The white dots indicate the position of paired band combinations corresponding to the maximum R 2 .

Figure 4 .
Figure 4. Correlograms of R 2 values for α, GPP max , midday averaged GPP, ε and NEE, on the one hand, and SR-type indices for Amplero, Neustift, Monte Bondone (both study years pooled) and all sites pooled on the other.The white dots indicate the position of paired band combinations corresponding to the maximum R 2 .

Figure 5 .Figure 6 .
Figure 5. Correlograms of R 2 values for α, GPP max , midday averaged GPP, ε and NEE, on the one hand, and SD-type indices for Amplero, Neustift, Monte Bondone (both study years pooled) and all sites pooled on the other.The white dots indicate the position of paired band combinations corresponding to the maximum R 2 .
sion model for NEE, but the VIS-vs.-VISband combination was selected for Amplero and the NIR-vs.-NIRcombination for Monte Bondone.Neustift performed well with NEE for NIR-vs.-NIRcombinations, but with an inverse relationship.The different types of indices (compare Figs. 6, S5 and S6) resulted in similar models.The different time resolutions gave different models (e.g.GPP, ε and NEE at Monte Bondone; compare Figs. 6 and S7, Figs.S5 and S8, or Figs.S6 and S9).

Figure 7
Figure7shows the results of the validation for each ecophysiological parameter and midday averaged fluxes and NSD-, SR-and SD-type indices against data from the validation sites.The models used in the validation are based on the best models determined for each site (i.e.Amplero, Neustift and Monte Bondone) and on pooling together the two alpine

Figure 8 .
Figure 8. Correlation between selected (a) NSD-, (b) SR-and (c) SD-type indices for the α, GPP max , midday GPP, midday ε and midday NEE (plots in the columns) and the total chlorophyll content for Monte Bondone in 2013.R 2 -coefficient of correlation; RMSE -root mean square error; R 2 cv -cross-validated coefficient of correlation; RMSE cv -cross-validated root mean square error.The solid red lines indicate the fitted models, and the dotted red lines represent the 95 % upper and lower confidence bounds.The selected bands for computing NSD-, SR-and SD-type indices are reported in brackets.

Figure 9 .
Figure 9. Results of the GA-rF method for band selection for Amplero, Neustift, Monte Bondone and all sites pooled for (a) α and GPP max , (b) midday average ε and CO 2 fluxes (NEE and GPP); (c) daily average ε and CO 2 fluxes (NEE and GPP).

Figure 10 .
Figure 10.Results of the GA-rF method for band selection for Amplero, Neustift, Monte Bondone and all sites pooled for dry phytomass, water and nitrogen content.

Table 1 .
Description of the study sites and period.

Table 2 .
Summary of the vegetation indices characteristics used in this study.

Table 3 .
Results of statistic of linear regression models between VIs and ecophysiological parameters: α, ε

Table 4 .
Results of statistic of linear regression models between VIs and midday average CO2 fluxes: NEE and GPP.R 2 -coefficient of determination; RMSE -root mean square error.

Table 5 .
Results of the correlation (R 2 -coefficient of determination) between the best NDS-, SR-and SD-type indices selected for α, GPP max , midday GPP, midday ε and midday NEE and dry phytomass, nitrogen and water content for Amplero, Neustift, Monte Bondone and all sites pooled.The selected bands for computing NSD-, SR-and SD-type indices are reported in brackets.Statistical significance is indicated as * (p < 0.05), * * (p < 0.01) and * * * (p < 0.001).