Remote Sensing of Grassland Biophysical Parameters in the Context of the Sentinel-2 Satellite Mission

This study investigates the potential of the Sentinel-2 satellite for monitoring the seasonal changes in grassland total canopy chlorophyll content (CCC), fraction of photosynthetically active radiation absorbed by the vegetation canopy (FAPAR), and fraction of photosynthetically active radiation absorbed only by its photosynthesizing components (GFAPAR). Reflectance observations were collected on a continuous basis during growing seasons by means of a newly developed ASD-WhiteRef system. Two models using Sentinel-2 simulated data (linear regression-vegetation indices (VIs) approach and multiple regression (MR) reflectance approach) were tested to estimate vegetation biophysical parameters. To assess whether the use of full solar spectrum reflectance data is able to provide an added value in CCC and GFAPAR estimation accuracy, a third model based on partial least squares regression (PLSR) and the ASD-WhiteRef reflectance data was tested.The results showed that FAPAR remained quite stable during the reproduction and senescence stages, and no significant relationships between FAPAR and VIs were found. On the other hand, GFAPAR showed clearer seasonal trends.The comparison of the threemodels revealed no significant differences in the accuracies of CCC and GFAPAR predictions and demonstrated a strong contribution of SWIR bands to the explained variability of investigated parameters.The promising results highlight the potential of the Sentinel-2 satellite for retrieving biophysical parameters from space.


Introduction
Spatial and temporal monitoring of biophysical parameters of vegetation canopies provides important information on their status and responses to changing environmental conditions [1].Canopy chlorophyll content (CCC) and fraction of absorbed photosynthetically active radiation (FAPAR) are widely applied in environmental studies concerning growth monitoring, stress detection, and yield estimation [2].Also, these variables are very important inputs for models of gas exchange between terrestrial ecosystems and the atmosphere and for models of vegetation productivity [3], which links atmospheric composition and vegetation processes [4].One of the most widely used approaches for gross ecosystem production (GEP) modelling using remote or proximal sensing data is based on the light use efficiency (LUE) concept developed by Monteith [5,6], according to which GEP is directly related to absorbed photosynthetically active radiation (APAR, mol m −2 s −1 ) and 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 FAPAR is the fraction of PAR absorbed by the vegetation canopy (-).From this general concept, APAR relates to vegetation structure and pigment pools, while  is linked to physiology through processes that influence the energy distribution within the photosynthetic system.However, the absorption and efficiency terms are linked in several ways, for example, through the chlorophyll and carotenoids or nitrogen content [7].
It is important to take into account the fact that both leaves and vegetation canopies include in their composition also nonphotosynthetic components (such as veins and walls present even in the green leaves or senescent leaves, branches, and stems at the canopy level) which will significantly affect the canopy FAPAR values [4].Therefore, partitioning of FAPAR into the fraction related only to photosynthesizing components ("green" FAPAR, GFAPAR) and estimating GEP with GFAPAR will most likely result in significant estimation improvement due to its higher consistency with the photosynthesis process [3,8].In fact, in the recent years, many studies tend to account for GFAPAR rather than FAPAR [1,9,10].
Furthermore, when no significant physiological stress is occurring, in ecosystems characterized by strong seasonal dynamics of green biomass such as croplands, grasslands, and deciduous forests [11][12][13][14][15], additional estimation of  may not be necessary due to its relation with the chlorophyll content.Also, when GEP is integrated over longer time scales, the importance of acquiring spectral data as a proxy of  is decreasing [16].
Due to spatial discontinuity of in situ measurements, effective monitoring of terrestrial biogeochemical cycles components at the global scale can be achieved only by satellite observations.However, ground-truth measurements are critical for calibrating and validating satellite data and atmospheric correction algorithms [31] and to ensure a clear and quantitative understanding of the satellite data [32,33].
Despite the importance of ground-truthing satellite observations, high spectral resolution measurements using standalone instruments are carried out only in a limited number of sites.The recent work within the COST Actions ES0903 "EUROSPEC" and ES1304 "OPTIMISE" highlighted the importance of enlarging the ground-truthing spectral networks and standardizing the observations, which have been so far carried out using system prototypes and nonstandardized setups and protocols [31,34].
ESA's recently launched (June 2015) Sentinel-2 satellite has begun to provide images of high spatial, spectral, and temporal resolution, offering high potential datasets for vegetation productivity and biophysical characteristics temporal and spatial monitoring which need to be validated at the eddy covariance (EC) towers.In particular, the MultiSpectral Instrument (MSI) of Sentinel-2 has made available a set of 13 spectral bands ranging from visible (VIS) and near infrared (NIR) to shortwave infrared (SWIR), featuring four bands at 10 m, six bands at 20 m, and three bands at 60 m spatial resolution (Table 1).In the full operational phase, a pair of twin Sentinel-2 satellites will deliver imagery every five days under cloud-free conditions (https://sentinel.esa.int/web/sentinel/missions/sentinel-2).Previous studies on evaluating the potential of the Sentinel-2 satellite for estimating vegetation biophysical variables relate mainly to agricultural ecosystems [35][36][37][38][39][40], with scarce examples of studies on forest [41] or grassland [38,42] ecosystems.Although further validation over crop ecosystems is yet required to reinforce the previous findings [36], extending the investigations to other ecosystem types (such as forests, grasslands, and peatlands) is a priority matter.
In this study, two years of CCC (2013 and 2014), FAPAR and GFAPAR (2014 and 2015), and simultaneous hyperspectral data acquired with an ASD-WhiteRef system [43] 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 evaluate the potential of the Sentinel-2 satellite (by simulating its spectral bands) for monitoring seasonal changes in grassland CCC, FAPAR, and GFAPAR.
(ii) To compare different approaches (correlation analysis and multiple regression) using Sentinel-2 simulated spectral vegetation indices (VIs) and reflectance data, respectively, to estimate CCC and GFAPAR.
(iii) To assess whether the use of full solar spectrum reflectance data provides an added value in CCC and GFAPAR estimation accuracy by means of partial least squares regression (PLSR) and the ASD-WhiteRef reflectance data.The vegetation of the investigated area is dominated by Festuca rubra (L.) (covering 25% of the area), Nardus stricta (L.) (13%), and Trifolium sp.(L.) (14.5%) and represents a typical low productive meadow of the Alps.The site is managed extensively by low mineral fertilization (applied in autumn) and one cut per year, occurring usually around mid-July [44].The maximum canopy height and the maximum plant area index (PAI) at the peak of the growth season (mid-June to early July) can reach approximately 30 cm [45] and 3 m 2 m −2 , respectively [43].

Material and Methods
The climate of this area is classified as subcontinental and is characterized by warm and wet summers, with mean annual temperature of 5.5 ∘ C (monthly averages ranging from −3.1 ∘ C in February to 14.3 ∘ C in July).The annual mean rainfall is 1244 mm, with maximum values in May (138 mm) and October (162 mm).The snow-free period lasts generally from early May to late October [46].

Hyperspectral Radiometric Measurements.
Hyperspectral reflectance datawere acquired on a continuous basis during the growing seasons from 2013 to 2015 by means of the ASD-WhiteRef system designed for continuous and unattended acquisition of canopy reflectance data [31,43].
The core of the single-beam WhiteRef system is a commercially available ASD FieldSpec Pro spectroradiometer (Analytical Spectral Devices, Inc., Boulder, Colorado, USA; wavelength range between 350 and 2500 nm; spectral resolution, FWHM of a single emission line, approximately 3 nm around 700 nm; sampling interval 1.4 nm : 350-1000 nm, 2 nm : 1000-2500 nm; a cubic spline interpolation function is used to calculate spectra at 1 nm intervals), which was upgraded for standalone operations.The spectroradiometer consists of an array of 3 spectroradiometers.The first (VNIR) covers the 350-1050 nm range, while SWIR1 and SWIR2 ranges are 900-1850 and 1700-2500 nm, respectively.The splice between the VNIR and the SWIR1 spectrometers is at 964/965 nm, where the response of the SWIR1 spectrometer is superior to that of the VNIR photodiode array.The splice between SWIR1 and SWIR2 occurs at 1769/1770 nm [43].
The ASD-WhiteRef system was installed in May 2013 at the FLUXNET Monte Bondone site (IT-MBo).The ASD spectroradiometer equipped with the fiber optic with a field of view of 25 ∘ was kept in a waterproof ventilated box mounted on the EC tower.The reference target radiance was measured by automatically sliding a white reference panel (WR) under the fiber optic.The fiber optic installation height (6 m) resulted in an optical canopy footprint diameter of about 2.7 m.In order to protect the WR from adverse weather conditions, as well as from light, dust, and insects, it was housed inside a waterproof box and ejected only during the measurements.Each acquisition was preceded by a reading of a dedicated wetness sensor signal, and in case of rainfall or dew the measurements were not conducted.Additionally, during each ejection and insertion phase, the WR was sprayed with compressed air in order to remove eventual dust/insects from its surface [31,43].
The incident solar radiance and reflected radiance (average of 25 scans) were collected every 5 minutes with dedicated LabVIEW software.Vegetation target measurements were sandwiched between two WR measurements conducted on average 30 seconds apart and the radiance of the WR panel at the time of the target measurements was estimated by linear interpolation.
Further details regarding the ASD-WhiteRef system setup, hardware, software, data acquisition, and quality controls can be found in Sakowska et al. [43].
In order to limit the signal noise ("sawtooth behavior") in the SWIR2 region caused by low and/or variable light levels, only reflectance data acquired during clear and stable sky conditions were included in the analysis.Stable sky conditions were defined as those when the change of the illumination conditions (the change between the WR1 and WR2 radiance measurements) was not higher than 10%.Incoming radiation quality was assessed by computation of Diffusion Index (DI, the ratio between diffuse and total incident PAR, measured with BF3H Sunshine Sensor, Delta T Devices Ltd., Cambridge, UK, and recorded at 1 min intervals by CR3000 data logger, Campbell Scientific Inc., Logan, Utah, USA).DI of 0.25 was adopted as the clear sky conditions threshold.Furthermore, reflectance data in the water absorption bands (1350-1460 nm and 1790-1960 nm) [47] and wavelengths above 2300 nm affected by noise due to low solar irradiance [48] were removed from the dataset, which resulted in 1668 available bands.
In addition, data were excluded (i) when the site was covered by snow, (ii) when precipitation was recorded 2 h prior to the measurements, and (iii) during the annual hay cutting and drying operations, which can take a few days, especially when the rainy weather conditions are not allowing for the immediate removal of the dried cut hay from the footprint of the ASD-WhiteRef system (and the EC tower) after the cut event.
The ASD FieldSpec Pro hyperspectral spectroradiometer spectral resolution is much higher than the one of the MSI of Sentinel-2.To match the spectral response of the two instruments, the ASD signal was resampled to new band positions using Gaussian models defined by the MSI fullwidth half-maximum (FWHM) values.The simulation was done by means of the resample2 function of the prospect R package [49,50].The twelve Sentinel-2 bands simulated in this study are presented in Table 1.
Finally, both the ASD-WhiteRef narrow band and Sentinel-2 simulated reflectance spectra were averaged over 2 h close to the solar noon (11:00 a.m.-1:00 p.m. of local solar time) to minimize solar angle effects and then used for computing the VIs.Although many different VIs were investigated, only the most commonly used and the best performing in CCC, FAPAR, and GFAPAR estimation are presented in this study.The list of the nine presented VIs is reported in Table 2, and charts (Figure A) displaying the seasonal trends of these VIs were included in the Supplementary Materials available online at http://dx.doi.org/10.1155/2016/4612809.
The schematic view of the ASD-WhiteRef hyperspectral data preprocessing flow is shown in Figure 1.

Index
Formulation Reference  ( The weight of sampled sediment was used to calculate the pigments concentrations per unit leaf mass (mg⋅g −1 ) and the weight of biomass per m −2 was used to obtain the total canopy chlorophyll content (CCC, mg⋅m −2 ).
where PARi, PARr, and PARt are incident, reflected, and transmitted PAR, respectively.In the studied grassland, the nonphotosynthetic component represents a considerable fraction of the aboveground phytomass during a substantial part of the growing season, and furthermore its contribution increases significantly towards the end of the growing season [1,10].Therefore, in order to estimate the fraction of PAR absorbed only by living and photosynthesizing components of the canopy ("green" FAPAR, GFAPAR), FAPAR was multiplied by a normalized (by scaling between 0 and 1) greenness index (GI, calculated as a ratio between the digital number values of green and the sum of red, green, and blue digital number values) derived from the analysis of time series of Phenocam (NetCam SC H.264, StarDot Technologies, Buena Park, CA, US) digital images acquired two times per day around solar noon, using the Phenopix R package [52].

Models for Vegetation Biophysical Parameters Estimation.
In this study, three groups of models were tested to estimate grassland biophysical parameters: (i) linear regression (model 1), (ii) multiple regression (MR, model 2), and (iii) partial least squares regression (PLSR, model 3).The first two approaches utilized Sentinel-2 simulated data, while the third model was making use of the full set of hyperspectral ASD-WhiteRef spectra.The first model assumed direct linear relationship between CCC, FAPAR, GFAPAR, and VIs.In the second approach, the interaction effects between different variables were explored by running a stepwise bidirectional multiple regression model, in which CCC or GFAPAR were set as dependent variables and reflectance values at twelve Sentinel-2 simulated bands as explanatory variables.In the third approach, the full set of ASD-WhiteRef spectra (1668 bands between 350 and 2300 nm) was used simultaneously to predict CCC and GFAPAR.The PLSR method iteratively transforms predictor (the ASD-WhiteRef reflectance data) and response (either CCC or GFAPAR) variables into latent vectors and generates band-wise calibration factors used to create a predictive linear model.PLSR aims at maximizing the covariance between independent and dependent variables, maintaining at the same time orthogonality in the factors derived from spectra [53].
The abovementioned models were tested for all the available years together (CCC: 2013 and 2014, 32 observations; GFAPAR and FAPAR: 2014 and 2015, 94 observations) in order to obtain general models for estimation of vegetation biophysical parameters.

Statistical Analysis.
Pearson's correlation analysis was used to test the significance of the relationships between VIs and biophysical data (CCC, FAPAR, and GFAPAR).
In the multiple stepwise bidirectional linear regression model, the set of Sentinel-2 simulated reflectance data that best fits the CCC and GFAPAR was chosen according to Akaike's information criterion (AIC) [54].In the following step, the variance inflation factor (VIF) [55] was computed to measure the degree of (multi)collinearity of the th independent variable with the other independent variables in the regression model.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 [20,56].The procedure was repeated until none of the VIF factors exceeded an 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  2 (Adj. 2 )-after adding a subsequent predictor variable to the multiple regression model-was significantly different from zero (at significance level  = 0.001) by means of the Fisher test [57].
The PLSR analyses were carried out by means of the PLS package of R software environment [58].The predictor PLSR matrix consisted of 1668 spectra and 32 and 94 observations for CCC and GFAPAR estimates, respectively.The optimal number of components in the PLSR analysis was determined by minimizing the leave-one-out cross-validated root mean square error of predictions (RMSEP CV ) [59].The variable importance projection (VIP) statistic was computed in order to evaluate the relative importance of reflectance at different wavelengths in the PLSR model [60].VIP values >1 indicated high importance to the PLSR model [61].
Each of the three model's coefficients was obtained by fitting each model against the measured variables (CCC, FAPAR, and GFAPAR).The main goodness-of-fit statistics (adjusted coefficient of determination-Adj. 2 , root mean square error-RMSE, RMSE normalized to the range of data-NRMSE, and probability value-) and crossvalidated statistics (Adj. 2  CV , RMSE CV , and NRMSE CV ) obtained with leave-one-out procedure were computed to compare the performance of the different models.
All the statistical analyses were performed by means of the R software (version 3.0.3,https://www.r-project.org/).

Seasonal Variation of Vegetation Biophysical Parameters.
Seasonal patterns of vegetation biophysical parameters were driven by both environmental variables (such as air temperature and precipitation) and grassland management.The grassland cut (occurring on 18 July (DOY 199), 25 July (DOY 206), and 7 July (DOY 188) in 2013, 2014, and 2015, resp.)split the growing seasons into two subperiods.
Figure 2 presents the variations of CCC (mg m −2 ) measured in the growing seasons of 2013 and 2014.The first sampling took place in the second half of May (DOY 143 and DOY 142 in 2013 and 2014, resp.) when the CCC was still very low, and the grassland had just begun to green up.CCC peaked around DOY 182 in 2013 and after that it started to decrease until the grassland was cut.In 2014, the highest value of CCC also occurred in the first part of July (DOY 192), but, in this case, the following gradual decline was not observed.After the cut event, the Monte Bondone canopy regrowth usually reaches the peak in the first part of September [57], and so it was in 2014.However, year 2013 was characterized by very low productivity of the second part of the growing season, caused by unusual drought affecting the grass growth after the cut event, which can be considered as an exception in this regard.The precipitation to temperature ratio for a 2-week period after the cut was 4 times smaller in 2013 compared to 2014.Moreover, the precipitation amount recorded during the entire month of July was equal to 53 mm, while the precipitation sum for the same month in 2014 was 6 times higher.
The seasonal variations of FAPAR and GFAPAR showed different dynamics in 2014 and 2015 (Figure 3).The second investigated year was characterized by earlier snowmelt; thus, new green plant shoots appeared approximately one month earlier than in 2014.In both measurement years, FAPAR showed a gradual increase during the vegetative state reaching the maximum around DOY 194 and DOY 165 in 2014 and 2015, respectively, and then remained practically invariant until the grassland cut event.The FAPAR dynamics showed the same pattern also in the second part of the growing seasons: after a progressive increase related to the canopy regrowth, this parameter reached the maximum (around DOY 258 and DOY 250 in 2014 and 2015, resp.) and generally remained at a constant level.On the contrary, the GFAPAR (characterized by lower absolute values) showed a decreasing course during both ripening stage (at the beginning of July in 2014 and in the second part of June in 2015) and senescence stage (starting in mid-September in 2014 and at the beginning of September in 2015).

Retrieval of Vegetation Biophysical Parameters with
Remotely Sensed Data.The linear regression analysis (model 1, Table 3) showed that the presented VIs explained at least 69% (Adj. 2 CV ) of the variability of CCC in the growing seasons of 2013 and 2014 considered simultaneously.
On the other hand, very low accuracy of the same model was reported for FAPAR estimates when the whole growing seasons of 2014 and 2015 were considered together (max.Adj. 2 CV = 0.15, Table 4).Converting FAPAR into GFAPAR resulted in a significant improvement of the general model 1 performance (Table 4).Adj. 2 CV was on average 94% higher and NRMSE CV 22% lower for GFAPAR than FAPAR estimates.For example, Adj.  2  CV of model 1 fed with RECI increased from 0.02 to 0.77 (-) and NRMSE CV decreased from 0.24 to 0.12 (-).
Taking into account the poor results obtained in the FAPAR estimations (discussed in Section 3.2) considering the general model 1, this variable was omitted from further analysis presented in this study.
The stepwise bidirectional procedure (MR, model 2) selected Sentinel-2 simulated reflectance () at 1610 and 945 nm as significant drivers of CCC and 865, 705, and 1610 for GFAPAR (Table 5).It is interesting to note that NIR and SWIR1 bands were included as important predictors of both CCC and GFAPAR.However, MR did not yield any improvement in the explained variance of GFAPAR (Adj. 2 CV = 0.78, RMSE CV = 0.09, and NRMSE CV = 0.12 for model 2; Table 5).Also, MR resulted only in a very slight improvement in the accuracy of CCC estimation compared to the best linear regression model based on VI.Adj. 2 CV increased from 0.88 (NDWI model) to 0.90, while RMSE CV decreased from 163.42 to 147.95 mg m −2 (Tables 3 and 5).
PLSR modelling indicated that 89% of the variation in CCC and 82% of the variation in GFAPAR were explained by algorithms derived from reflectance at 1668 ASD-WhiteRef narrow bands (Adj. 2 CV , Table 6).RMSE CV of the PLSR models was equal to 148.10 mg m −2 for CCC and 0.08 (-) for GFAPAR (Table 6), requiring 3 and 5 latent components, respectively, to achieve this predictive capacity (Figure 4).The VIP values displayed fairly consistent patterns across the spectrum, as well as similar ranges for both models (varying from 0.18 to 1.88 for CCC PLSR model and from 0.11 to 1.86 for GFAPAR model) (Figure 5).For example, wavelengths in the region between 735 and 1130 nm were almost uniformly distributed and particularly important for both variables.A close agreement was found also for the SWIR region from 1590 to 1780 nm; however, the VIP threshold of 1 was exceeded only for the GFAPAR model.A key difference was observed solely in the first section of the red-edge region (690-735 nm), which appeared to be important only in the GFAPAR model with a local VIP peak occurring at around 710 nm.In the CCC model, this peak occurred at a slightly shorter red-edge wavelength (∼700 nm) and was characterized by 2-fold lower absolute value.
Scatter plots of observed versus predicted CCC and GFAPAR using Sentinel-2 simulated (model 2) and full ASD-WhiteRef (model 3) spectra are presented in Figure 6.

Discussion.
The ASD-WhiteRef was demonstrated to be a reliable system for unattended, continuous, and long-term hyperspectral reflectance measurements, as no significant failures during the three years of observations (2013-2015) were recorded.There are only few hyperspectral systems available worldwide that succeeded in continuous proximal sensing from EC towers (such as UNIEDI operating in Hyytiälä, Finland [62]; Hyperspectral Irradiometer installed in Torgnon, Italy [63]; tram system from Skye Oaks, USA [64]; AMSPEC-MED operating at Las Majadas, Spain [65]).Furthermore, to our knowledge, only one system for continuous hyperspectral measurements including the SWIR region (however not covering the SWIR bands above 1800 nm) was tested before [66,67].For this reason, the ASD-WhiteRef stands out for its suitability for in situ spectral validation of various satellite products incorporating also SWIR bands, such as, for example, Sentinel-2, Landsat 8, or WorldView-3.
From the presented data, it appeared to be clear that FAPAR failed in capturing the vegetation greenness decrease related to both ripening and senescence stages, as the main driver of FAPAR is the total phytomass.During these stages, the contribution of necromass increased significantly; therefore, although the canopy was still intercepting PAR, it was used for photosynthesis to a much lower extent.As observed in other studies, no significant relationship was observed between this parameter and VIs (Table 4) [1].On the contrary, GFAPAR showed to be able to represent well the Monte Bondone grassland greenness seasonal trends.Model 1 showed high performance and good agreement for both CCC and GFAPAR, which is a consequence of the strong linear   relationship between these biophysical parameters (Figure 7).In fact, according to Peng et al. [9], CCC is the main driver of GFAPAR.
The results of model 1 confirmed the findings of previous studies indicating that (i) the most commonly used VI, NDVI, tends to saturate under conditions of moderate-to-high aboveground biomass due to the saturation of reflectance in the visible bands and due to the limitation of the normalized difference approach [28,68,69]; (ii) linearizing the relationship between VIs and CCC or GFAPAR can be obtained by parameterizing the model either with VIs based on the red-edge part of the spectrum or with modified variants of NDVI, such as, for example, WDRVI [1,28,70].
Regarding the highest accuracy of CCC estimation obtained with NDWI, the good performance of VIs based on NIR and SWIR bands might be attributed to the strong sensitivity of reflectance at SWIR bands to the foliage water content [71], which in turn correlates well with leaf chlorophyll content [72].Consequently, at full canopy cover, SWIR-based VIs indirectly relate to the quantity of aboveground biomass and other vegetation biophysical parameters.Although incorporating SWIR bands into VIs formulation was proposed already many years ago and good relationships with biophysical variables were obtained [67,[73][74][75][76][77], the potential of using these bands in both proximal and remote sensing of vegetation is still not fully explored and strongly underutilized.Further studies are needed to highlight the potential of using SWIR bands for improving the estimation of biophysical parameters and, in turn, GEP.The Sentinel-2 simulated NIR and SWIR reflectance were selected as the main drivers of both CCC and GFAPAR also within the MR method (model 2).
Using MR reflectance-based approach instead of the bestperforming VIs did not lead to significantly improved results in investigated parameters estimation.In more detail, MR resulted only in a 2% increase in explained variance of CCC and no improvement in GFAPAR predictions, respectively (Tables 3, 4, and 5).Despite the fact that the red-edge bands were not included in the final MR CCC model (705 was selected as the third predictor, and its partial  2 of 0.003 did not result in statistically significant change in the cumulative Adj. 2 ), the successful performance of VIs including 705 in their formulation confirms the important role of this part of the spectrum in monitoring the dynamics of CCC.The VIP statistics associated with each wavelength were analyzed in order to identify the wavelengths that contribute the most to predicting the variable of interest (in this case, either CCC or GFAPAR).These analyses revealed that CCC and GFAPAR models did not utilize all regions of the solar spectrum.Moreover, a high degree of coincidence was observed between the PLSR selected narrow bands and both Sentinel-2 simulated bands selected within the MR and bands used in the best-performing VIs models.This means that, in general, the same spectral regions showed to be important in all the methods.A clear example of this degree of coincidence was CCC, where the best fit in model 1 was obtained using NDWI based on bands centered at 865 and 1610 nm.Both NIR and SWIR bands were chosen also within the MR model and such wavelengths are included in the area where PLSR VIP achieved values higher than 1.The reason why the relevant bands for CCC estimation were located in the NIR and SWIR regions is related to the fact that, in the Monte Bondone grassland, most of the variation in CCC was driven by the variation in biomass, as the seasonal changes of leaf chlorophyll content (data not shown) were too low to influence the CCC significantly.This is in agreement with previous studies demonstrating a strong contribution of SWIR bands to the explained variability of vegetation biophysical parameters by means of MR (CCC, LAI [20]) and PLSR models (nitrogen, biomass, and canopy water content [78]).
It was observed that, for both investigated variables, predictive models based on PLSR did not produce considerably higher accuracies compared to VIs or MR models based on Sentinel-2 bands (CCC Adj. 2 CV = 0.88, 0.90, and 0.89 for VIs-based, MR, and PLSR models, resp., and GFAPAR Adj. 2 CV = 0.77, 0.78, and 0.82 for VIs-based, MR, and PLSR models, resp.).Although this is in contrast with some other studies [79,80], our results indicate that VIs or few spectral bands alone are able to provide sufficient and accurate information on vegetation biophysical parameters seasonal changes in dynamic grassland ecosystems (Figure 6).
The promising results obtained with Sentinel-2 simulated bands underline the potential of this satellite mission for retrieving the investigated parameters from space.However, it has to be stressed out that similar studies, conducted in different types of ecosystems, and further testing of the above models for GEP estimations are required in order to fully evaluate the quality of Sentinel-2 future products.

Conclusions
This study aimed at assessing the performance of the newly developed automatic ASD-WhiteRef system in acquiring hyperspectral reflectance data on a continuous basis and in using this ground-truth database to investigate the potential of the Sentinel-2 satellite to monitor biophysical parameters (such as CCC, FAPAR, and GFAPAR) in a dynamic subalpine grassland ecosystem.For this purpose, the modelling results obtained only on the basis of the Sentinel-2 simulated data (correlation with VIs and multiple regression with reflectance at 12 bands as predictor variables) were compared with predictive fit achieved with the PLSR modelling using narrow-band solar reflectance.The lack of failures during the 3-year observation period proved the ASD-WhiteRef system reliability in continuous and long-term field spectroscopy measurements.The use of full spectrum instead of Sentinel-2 simulated data approach did not lead to considerably improved results in CCC and GFAPAR estimation.Although further research in other ecosystems is required, such results, considering the high spatial resolution and short revisit time of the Sentinel-2 satellite, indicate a strong potential of this mission for retrieving vegetation biophysical parameters in dynamic ecosystems, which needs to be confirmed by satellite observations when appropriate time-series will be available.

Figure 1 :
Figure 1: Schematic view of the ASD-WhiteRef hyperspectral data preprocessing flow.

Figure 2 :
Figure 2: Seasonal courses of total canopy chlorophyll content (CCC, mg m −2 ) in the growing seasons of 2013 and 2014.

Figure 3 :
Figure 3: Seasonal courses of fraction of photosynthetically active radiation absorbed by vegetation canopy (FAPAR, -) and fraction of PAR absorbed only by photosynthesizing components of the vegetation canopy (GFAPAR, -) in the growing seasons of 2014 and 2015.

Figure 4 :
Figure 4: Cross-validated root mean square error of predictions (RMSEP CV ) curves: RMSEP CV as a function of the number of components in the partial least squares regression (PLSR) models for total canopy chlorophyll content (CCC, mg m −2 ) and fraction of photosynthetically active radiation absorbed by photosynthetic vegetation components (GFAPAR, -).

Figure 7 :
Figure 7: Relationship between the fraction of photosynthetically active radiation absorbed by photosynthetic vegetation canopy components (GFAPAR, -) and the total canopy chlorophyll content (CCC, mg m −2 ) in the growing season of 2014.Adj. 2 : adjusted coefficient determination (-).

Table 1 :
Specification of the MultiSpectral Instrument (MSI) of the Sentinel-2 satellite system.The twelve bands simulated in this study are shown in bold.
At the Monte Bondone site, vegetation sampling was undertaken at different vegetation development stages in two campaigns carried out during the growing seasons of 2013 and 2014, with a total of 12 observation dates in 2013 and 20 observation dates in 2014.During each date, three 10 m long and 0.1 m wide stripes (representing 1 m 2 ) were cut within the footprint of the EC tower (and nearby the footprint of the automatic ASD-WhiteRef hyperspectral system) in order to collect phytomass samples (50 g per each stripe).Collected phytomass was stored in sealed plastic bags and kept fresh in an ice chest while being transported to the laboratory, where it was divided manually into biomass and necromass.Both biomass and necromass were then dried at 80 ∘ C for 48 hours and green herbage ratio (GR = biomass/(biomass + necromass)) was calculated.Furthermore, randomly chosen additional green leaf samples (>30 g) were used for chlorophyll measurement and characterization by UV-VIS spectroscopy.In the first step, the green tissue samples were ground in the presence of liquid nitrogen.

Table 3 :
Summary of the statistics (: number of observations; Adj. 2 : adjusted coefficient of determination; Adj. 2 CV : cross-validated Adj. 2 ; RMSE: root mean square error; RMSE CV : cross-validated RMSE; NMRSE: normalized RMSE; NRMSE CV : cross-validated NRMSE) of the linear regression (model 1) between measured canopy chlorophyll content (CCC, mg m −2 ) and spectral vegetation indices (VIs) considering both years of observation together (2013 and 2014).The three best-fitting models are printed in bold.The best performing model is additionally highlighted in italic.All the regressions were statistically significant ( < 0.001).

Table 5 :
Summary of the multiple regression models: partial adjusted coefficient of determination (partial Adj. 2

Table 6 :
Summary of partial least squares regression modeling and cross-validation: number of observations (), number of PLSR model components (PLSR comp.), adjusted  2