VIS-NIR, Red-Edge and NIR-Shoulder Based Normalized Vegetation Indices Response to Co-Varying Leaf and Canopy Structural Traits in Heterogeneous Grasslands

Red-edge (RE) spectral vegetation indices (SVIs)—combining bands on the sharp change region between near infrared (NIR) and visible (VIS) bands—alongside with SVIs solely based on NIR-shoulder bands (wavelengths 750–900 nm) have been shown to perform well in estimating leaf area index (LAI) from proximal and remote sensors. In this work, we used RE and NIR-shoulder SVIs to assess the full potential of bands provided by Sentinel-2 (S-2) and Sentinel-3 (S-3) sensors at both temporal and spatial scales for grassland LAI estimations. Ground temporal and spatial observations of hyperspectral reflectance and LAI were carried out at two grassland sites (Monte Bondone, Italy, and Neustift, Austria). A strong correlation (R2 > 0.8) was observed between grassland LAI and both RE and NIR-shoulder SVIs on a temporal basis, but not on a spatial basis. Using the PROSAIL Radiative Transfer Model (RTM), we demonstrated that grassland structural heterogeneity strongly affects the ability to retrieve LAI, with high uncertainties due to structural and biochemical PTs co-variation. The RENDVI783.740 SVI was the least affected by traits co-variation, and more studies are needed to confirm its potential for heterogeneous grasslands LAI monitoring using S-2, S-3, or Gaofen-5 (GF-5) and PRISMA bands.


Introduction
Canopy structural organization describes the three-dimensional geometric distribution of the aboveground photosynthetic and non-photosynthetic vegetation components [1]. Canopy structure is described by plant traits (PTs) such as leaf area index (LAI), aboveground biomass (AGB) and other canopy and leaf structural traits such as leaf angle distribution (LAD), gap fraction, leaf clumping, the proportion of photosynthetic and non-photosynthetic elements [2][3][4], specific leaf area (SLA) and leaf dry matter, which can influence absorption and scattering light dynamics [5,6].
Remote sensing can provide fundamental spatial and temporal information, which can be used in monitoring PTs related to plant biochemistry, photosynthetic processes and canopy structure. During the last years, the proximal sensing approach was used to fill the scaling gap between leaf and satellite measurements, linking vegetation characteristics and spectral responses from the leaf level to increasing pixel sizes [7][8][9]. While the visible (VIS) and shortwave infrared (SWIR, 1100-2500 nm) parts of the reflectance spectrum are mainly determined by pigments and water content absorption, respectively, in the near infrared (NIR, 750-1400), reflectance is high compared to the VIS domain because individual leaves and whole plant canopies strongly scatter NIR, and the degree of NIR scattering is driven by the internal leaf structure alongside with canopy structure and the ratio between green and non-photosynthetic components [5].
In herbaceous plants, LAI is a spatially-and temporally-dynamic key trait related to ecosystem functions (e.g., productivity and evapotranspiration), and remote sensing data have been widely used to capture its variability at various scales [10,11]. However, the structural and biochemistry variability among leaves, plants and ecosystems-particularly in natural grasslands, characterized by extreme heterogeneity [9]-is strongly affecting our ability to link spectral variation and LAI. Simultaneously with LAI, factors such as leaf anatomy and LAD are also varying in space and time (e.g., across heterogeneous canopies or due to phenological changes), and this has a significant and often unpredictable impact on scattering across the spectrum. When more structural traits co-vary, LAI estimation based on spectral data may be challenging [5], as reflectance is sensitive to multiple leaf and canopy traits and disentangling LAI from structural and biochemical drivers is difficult [12]. The impact of vegetation structural heterogeneity on the ability of different optical-based models to retrieve LAI has not been sufficiently described in the literature, and new knowledge is needed to quantify the uncertainties of such models and disentangle the impact of structural and biochemical heterogeneity on LAI estimations.
One of the main remote sensing approaches to estimate PTs focuses on empirical models, which are used to quantify relationships between PTs and canopy reflectance or spectral vegetation indices (SVIs) [13,14]. Over the last decades, the SVIs-based methods traditionally used combinations between NIR and VIS bands [15][16][17] to estimate LAI. However, several authors demonstrated that SVIs based on the NIR and red-edge (RE) spectral domains can significantly improve LAI estimations [13,18,19]. The RE is defined as the spectral region between 680 and 750 nm where a sharp change in the vegetation reflectance can be observed [20][21][22]. Such spectral domain is on the transition between chlorophyll absorption in the red wavelengths and leaf/canopy scattering in the NIR wavelengths [22,23]. The use of narrow-band SVIs is becoming extremely important in the context of the recently-launched satellite missions such as Sentinel-2 (S-2) and Sentinel-3 (S-3), as well as within the context of new hyperspectral missions such as Gaofen-5 (GF-5) and PRISMA.
The Radiative Transfer Models (RTMs) based on physical principles [24] are capable of simulating the interaction of light with vegetation at leaf and canopy levels and provide an explicit method for estimating the vegetation biophysical variables from canopy reflectance [25]. Another approach to retrieve LAI is based on the inversion of RTM which simulate the interactions of radiation with vegetation elements and the soil [26]. Such inversion approaches demonstrated to be challenging when the model is not well-suited for the observed vegetation type [26] and when suitable ancillary data and regularization methods to optimize the inversions for an efficient parameterization are lacking [27]. In addition, RTMs are mainly focused on chlorophyll and not on other pigments such as brown pigment content (polyphenols; Cbrown) which play an important role in shaping the spectral response of grasslands at varying phenological stages. This trait is rather poorly analyzed in the literature and the RTM parameter itself lacks of a proper physically meaning, and it is thus not measurable with field observations [28].

SVIs and LAI Empirical Models: Does Trait-Covariation Matter?
SVIs combining NIR and RE bands have been used in the last years to estimate both structural and biochemistry-related traits. The impact of structural traits on reflectance in RE spectral regions and on RE SVIs, however, is arousing controversy. According to some authors, RE SVIs are able to reduce structure-related artifacts in retrieved biochemical-related traits, as the RE is thought to be sensitive to chlorophyll content and largely unaffected by structural properties [12]. More recently, Peng et al. [41] showed that SVIs based on NIR band and RE band of 740 nm (RE740) demonstrated to be good predictors of canopy chlorophyll in crop types with contrasting canopy structure. The authors concluded that using NIR and RE740 band combinations provided good chlorophyll estimation due to the "reduced sensitivity of the RE to hysteresis driven by different canopy and leaf structures". Conversely, many authors agree on the fact that both biochemistry and structure contribute in determining the spectral response also within the RE domain [12,[42][43][44]. At the same time, Ollinger [5] pointed out that the variation in SVIs involving VIS and NIR bands is often driven to a greater extent by the variation in NIR reflectance than by variation in the VIS reflectance. In this context, the impact of canopy structure on SVIs still needs to be clarified.
The aforementioned controversies highlight that the response of structural traits, in combination with biochemical ones, on reflectance and on the SVIs models accuracy has not been fully characterized for different canopy types. To this regard, more efforts are needed to characterize the spectral response in the RE reflectance domain-in heterogeneous canopies-at varying biochemistry and structure, and more specifically to analyze the impact of canopy structural and leaf traits co-variation on SVIs-LAI relationships [45]. Such characterization is particularly important for multi-species natural grassland canopies, characterized by high spatial heterogeneity and temporal phenological changes, as presented in Section 2.3.4. Darvishzadeh et al. [46] concluded that LAI estimation in grasslands with mixed species and heterogeneous architecture is challenging, and that detailed investigations are needed to assess the suitability of different remote sensing models when many combinations of several species are observed.
In this framework, the objectives of the present study were: • To compare the ability of different SVIs including information from the RE and the NIR-shoulder spectral regions to estimate LAI at both temporal and spatial scales using ground hyperspectral data • To analyze the potential of Sentinel band combinations across the RE and the NIR-shoulder spectral regions using S-2 and S-3 simulated bands to estimate LAI in two grassland ecosystems of the Alps with contrasting structures • To determine the impact of grassland structural and biochemical heterogeneity on LAI estimations by analyzing the spectral reflectance response to co-varying biochemical and structural leaf and canopy traits across the RE and NIR-shoulder spectral domain using an RTM approach • To identify the best performing S-2 and S-3 SVIs for monitoring grasslands with heterogeneous structure by describing the impact of co-varying leaf and canopy structural traits on the relationships between LAI and SVIs calculated from S-2 and S-3 bands, as well as comparing RTM and empirical approaches.

Study Sites
This study was conducted in two different grassland sites of the Italian and Austrian Alps (Figure 1), situated in the subalpine and montane vegetation belts, respectively, and characterized by contrasting management types.
The first site is a permanent meadow at the Viote del Monte Bondone plateau (46.0147N-11.0458E, Italian Alps). The plateau meadow is managed mostly extensively, with low mineral fertilization and one cut annually around mid-July. In its central part, the meadow hosts a Fluxnet Eddy Covariance (EC) tower (IT-MBo, Italy). Due to its heterogeneous management and orography, the plateau is characterized by the presence of different grassland types characterized by extremely varying LAI and biomass [9]). Such heterogeneous meadows include several different vegetation types (two most abundant associations, one of them including two variants [9]). The Sieverso-Nardetum strictae association (average aboveground biomass: 236 g·m −2 [9]) coves a high portion of the plateau, including the EC footprint, and is characterized by short canopies and not intensively managed. The Scorzonero Aristatae-Agrostidetum tenuis (average aboveground biomass: 384 g·m −2 ) association is also very common on the plateau; it grows on calcareous soils and includes very productive species, which, although typical of much lower altitudes, can be frequent in some of the most fertile and well-exposed areas of the plateau (e.g., Arrenatherum elatius and Dactylis glomerata). The plateau, in its Eastern part, consists of small peatland associations of Caricion fuscae and Caricion davallianae characterized by very high productivity [9].
At the IT-MBo, several different ecosystems with extremely contrasting structures and productivity can be found within a few hundreds of meters distance ( Figure S1). The vegetation in sampling Plots 47 was very tall (maximum height reached 120 cm) and their structure was more erectophile, which is representative of the Arrenatherion alliance considered one of the most common vegetation types of Central Europe, at medium-lower altitudes. Vegetation in Plots 1-3 was very short with small dense Nardus tussocks and a limited number of scattered Festuca spikes (maximum height 50 cm). These plots represent, in their physiognomy, a typical grassland (meadow or pasture) on lime-deficient or acidified soils from the lower mountains of Europe, up to the lower alpine belt above the timberline. Vegetation in Plots 8-10 showed an intermediate height (maximum 90 cm) and was characterized by the presence of medium-productive grasses (such as Agrostis Tenuis and Trisetum flavescens) and forbs. This vegetation type is representative of a typical Centro-European species-rich mesophile grasslands at intermediate altitudes (between the Arrenatherion and Nardus types) of the montane and sub-alpine levels. The IT-MBo meadow soil can be classified as a Typic Hapludalfs, lyme loamy, mixed, mesic with the following characteristics in the 0-30 cm horizon: total soil organic content (SOC) = 9.4 ± 0.4 kg C m −2 ; total N = 0.29 ± 0.02 kg N m −2 ; and soil bulk density = 0.79 ± 0.29 g cm −3 [47].
The second site (AT-Neu, Fluxnet site) is a meadow located in Neustift (47.1162 N, 11.3204 E, Tyrol, Austria) classified as a Pastinaco-Arrhenatheretum [48]. The meadow shows very high productivity (with aboveground biomass values of up to 700 g·m −2 [36]) and is intensively managed with three cuts in mid-June, at the beginning of August and at the end of September. The vegetation type includes a few dominant graminoids (Dactylis glomerata, Festuca pratensis, Phleum pratense and Trisetum flavescens) and forbs which are abundant in terms of biomass and are characterized by wider leaves such as Ranunculus acris, Taraxacum officinale, Trifolium repens, Trifolium pratense and Carum carvi ( Figure S2 [48]). The 1-m deep soil profile in AT-Neu meadow soil has been classified as a Fluvisol (FAO classification) with a very thin organic layer (up to 2 cm) and beyond which it is described as a sandy loam. The vegetation roots reach down to 50 cm, but 80% of them are concentrated in the upper 13 cm of the soil [48].

Ground Biophysical Measurements
In the growing season of 2013, the fraction of absorbed photosynthetically active radiation (fAPAR) in IT-MBo was quantified at different vegetation development stages by periodic measurements (8 measurements in the period between May and July 2013) of incoming, reflected and transmitted (2 repetitions) PAR using the SunScan Canopy Analysis System (Delta T Devices Ltd., Cambridge, UK). In 2014, the fAPAR in IT-MBo was estimated using continuous measurements in the period between June and July 2014 by means of the Li-COR PAR sensors (Li-COR Inc., Lincoln, Nebraska, USA). Two  Quantum sensors were installed above the canopy level, measuring both incoming and reflected PAR, while the third sensor  Line Quantum) was placed at the ground level, measuring PAR transmitted through the vegetation canopy.
The temporal patterns of fAPAR in AT-Neu in the growing season of 2018 were computed using continuous measurements in the period between April and May 2018 carried out with BF2H (Delta T Devices Ltd., Cambridge, UK) and QSO-Sun (Apogee Instruments, Inc., Logan, UT, USA) sensors measuring incoming and reflected PAR, respectively, and two SQ-316 Line Quantum sensors (Apogee Instruments, Inc., Logan, UT, USA) measuring transmitted PAR. In both ecosystems (IT-MBo, AT-Neu), the temporal scale measurements of transmitted PAR with line sensors were performed within the footprint of the ASD-WhiteRef hyperspectral system. All the PAR data were recorded by a data logger (CR3000 in IT-MBo 2014, CR1000 in AT-Neu 2018; Campbell Scientific Inc., Logan, Utah, USA) at 1-min intervals and averaged over solar noon (11:00-13:00 local solar time) to match the time period used for vegetation spectral properties calculations. The spatial patterns of fAPAR in the biomass peak season of 2017 in IT-MBo were quantified at 10 different grassland plots using the SunScan Canopy Analysis System (Delta T Devices Ltd., Cambridge, UK). Two fAPAR measurements were conducted along the transect of the sampling plots with the SunScan instruments, which consist of 64 PAR sensors embedded in a 1-m-long portable probe positioned underneath the grass canopy. In each plot, the probe was positioned along 2 diameter axes (one orthogonal to the other), thus providing 64 × 2 = 128 individual measurements, which were averaged prior to fAPAR calculation. Additionally, two measurements of reflected and incoming light were performed with the probe right above the canopy. The summary of biophysical measurements acquired in the study is presented in Table 1.  The fAPAR was calculated as: where PARi, PARr and PARt are incident, reflected and transmitted PAR, respectively. LAI was estimated non-destructively by an indirect method based on canopy PAR transmission using PAR sensors data and a physical model of radiative transfer "INVERSION" [49]. The model of radiative transfer considers the canopy as a turbid medium in which multiple scattering occurs due to elements of turbidity (phytoelements). The model uses four adjustable parameters: (1) phytoelement dispersion factor; (2) phytoelement reflection; (3) transmission coefficients; and (4) soil reflection. The LAD function for the investigated vegetation type was assumed as erectophile [50], and it was used in "INVERSION" parameter settings. None of the other parameters (dispersion, reflectivity, transmissivity and soil reflectivity) were determined for the investigated canopies, thus default values [49] were used in the model to estimate the LAI.

Hyperspectral Reflectance Measurements
Hyperspectral reflectance data at both sites (in IT-MBo in the growing season of 2013 and 2014 and in AT-Neu in the growing season of 2018) were acquired (Table 1) on a continuous basis by means of the ASD-WhiteRef system [35], allowing measurements in the wavelength range between 350 and 2500 nm. The installation height (6 m in IT-MBo and 2.6 m in AT-Neu) and the system FOV (25 deg) resulted in an optical canopy footprint diameter of about 2.7 and 1.1 m in IT-MBo and AT-Neu, respectively. Both the ASD-WhiteRef narrow band reflectance spectra were averaged over 2 h close to the solar noon (11:00-13:00 local solar time) to minimize the solar angle effects and then used for computing the SVIs.
In addition, in 2017, spectral observations were performed in the spatial domain in IT-MBo plateau by deploying the ASD FieldSpec Pro spectroradiometer (Analytical Spectral Devices, Inc., Boulder, CO, USA) equipped with a fiber optic with the field of view of 25 deg and a white reference panel on a custom-made aluminum portable system (of a height of 2 m and the resulting FOV of ca. 0.9 m, Figure S3) designed for periodic nadir observations. The portable system consisted of a vertical pole equipped with two horizontal arms (installed one above another; the upper one was fixed and served the fiber optic assembly; the second arm-placed 20 cm below-was rotatable and allowed installation of the white reference panel) enabling alternating observations of the reference and the vegetation target and a ground structure ensuring system stability.
The spatial observations were carried out in July 2017 in 10 different grassland plots covering the aforementioned vegetation types and characterized by quite diverse canopy structure ( Figure S1) and productivity [9]. The canopy structure of the investigated grasslands was ranging from very short and dense canopies with Nardus Stricta (typical of high altitudes: V1-V3; Figure S1) to very productive grasslands with tall species such as Dactylis Glomerata and Arrenatherum elatius (V4-V7; Figure S1). For the spatial observations, the field campaign was conducted at the biomass peak to explore as many plots as possible within a short period (3 days) to avoid grass phenological changes (browning of the vegetation) and to ensure clear-sky conditions. The number of spatial scale observations was in general lower than that of temporal observations, as the latter were carried out on a more continuous basis.

Best Band Combination and Hyperspectral NDIs
To identify the band combinations most sensitive to LAI, the normalized difference index (NDI = (B1 − B2)/(B1 + B2), where B1 and B2 refer to reflectance values at specific ASD-WhiteRef bands) was calculated with all possible two-band combinations based on ASD ground hyperspectral data within the VIS and NIR spectral range (350-1000 nm) at 1-nm resolution. The definition of the two-band NDIs used in this study is presented in Figure 2. The 2D correlograms for B1 (350-1000 nm) versus B2 (350-1000 nm) highlights the performance of SVIs in the different spectral regions: the VIS region (400-680 nm), the RE region (680-750 nm) and the NIR-shoulder region (750-900 nm). The correlograms were generated using the "Spectral Indices (SI) assessment toolbox" from a modular software package ARTMO (Automated Radiative Transfer Models Operator) [51]. The SI ARTMO toolbox facilitates the assessment of spectral-domain prediction efficiency based on the adopted SI formulation and generates the correlation (R 2 ) matrices with all possible two-band combinations between measured and estimated values for any biophysical parameter [52]. The 2D correlograms illustrate the performance of all two-bands normalized difference combinations in retrieving LAI.

Multispectral Sentinel 2 and 3 SVIs
To obtain the reflectance values in the S-2 and S-3 bands, a simulation approach considering the average reflectance over the bandwidth of the respective S-2 and S-3 bands was adopted following Peng et al. [41]. SVIs based on S-2 and S-3 simulated bands (Table S1) in the VIS-NIR, RE and NIR-shoulder spectral regions were calculated, and their potential for estimation of grasslands LAI was tested. Two SVIs (NDVI 865.665 and MTCI) were based on VIS and NIR reflectance (referred to as VIS-NIR SVIs), three SVIs (RENDVI 783.740 , RENDVI 783.705 and RENDVI 865.740 ) were based on RE and NIR reflectance (referred to as RE SVIs) and the other two SVIs (NSDI 779.754 and NSDI 865.783 ) were based on the NIR-shoulder reflectance (≥ 750 nm; referred to as NIR-shoulder SVIs). We chose to compare the performance of two VIS-NIR SVIs (as a reference), alongside the performance of three RE making use of a RE band and two NIR-shoulder indices ( Table 2).

Global Sensitivity Analysis
Different factors (e.g., background soil and Solar-object-sensor geometry parameters) besides the biochemical and structural traits affect the canopy reflectance. To investigate the impact of structural PTs on reflectance, the ARTMO Global Sensitivity Analysis (GSA) toolbox was used [45]. The Variance-based GSA enables evaluating the relative importance of the spectral bands by identifying relative contribution (SI) of the key input variables that drive RTM spectral outputs [49]. The PROSAIL model input parameters used in this study for the GSA analysis are summarized in Table 3. The biochemical and structural input variables of the PROSAIL RTM model co-varied between the minimum and maximum value and for the geometrical parameters used the fixed averaged value. Table 3. Input parameters used in the global sensitivity analysis (GSA). LAI was estimated from the canopy fAPAR measurements at IT-MBo. The ranges of other input parameters (N, Cab, Car, Cw, Cm, Cbrown, soil and LAD) were selected based on existing literature [54][55][56][57][58].

SVIs Performance Using Simulated Spectra Under Different Temporal and Spatial Scenarios
RTMs are used to understand light interception by plant canopies and for the interpretation of vegetation reflectance in terms of biophysical characteristics [54]. In this study, the PROSAIL RTM [52] was used to assess the influence of leaf and canopy PTs on the SVIs calculated based on VIS, RE and NIR-shoulder spectral regions. PROSAIL couples two separate models: (a) the PROSPECT leaf optical model; and (b) the SAIL canopy reflectance model [28]. PROSAIL can simulate the canopy bidirectional reflectance in the spectral range between 400 and 2500 nm as a function of up to 16 input parameters. In this study, the input parameters were constrained according to the ranges of biophysical parameters measured in the field and literature-based values (Table 3). A look-up table (LUT) generated using the Latin hypercube sampling (LHS) method [59] was adopted to achieve a uniform distribution of the model inputs within the given boundaries for each scenario. The LHS method divides the cumulative density function into n bins of the same size from which data are randomly selected.
To investigate the impact of structural (LAI, LAD, leaf structural parameter (N) and leaf dry matter (Cm)) and biochemical (leaf chlorophyll content (Cab), carotenoid content (Car), brown pigment content (Cbrown) and leaf water content (Cw)) input parameters of PROSAIL on canopy spectral reflectance, various scenarios were assumed.
Using co-variation between input variables to constrain the model, we adopted parameters ranges observed in the literature (Table 3). With this approach, we aimed at exploring the impact of traits co-variation at both the spatial scale (due to vegetation heterogeneity) and the temporal scale (due to grassland phenological changes). Hence, running PROSAIL with a 100-iteration step, we analyzed four scenarios (1t-4t) based on the LAI range observed at the temporal scale at IT-MBo (during the growing period), as follows: (1t) LAI was set to vary between the minimum (at the beginning of the season) and the maximum value (at the biomass peak) and other PROSAIL parameters were kept constant (average value). This corresponds to a theoretical scenario where only grassland canopy LAI changes during the growth period, and no change is assumed in both other canopy structural parameters and biochemistry.
(2t) LAI and LAD values were co-varying between the minimum and the maximum value (Table 4) during the growth period and the rest of the PROSAIL parameters were kept constant (average value). This scenario can be considered representative of grasslands with slight canopy architectural and leaf structural dynamics at different growth stages.
(3t) LAI, LAD, N and Cm values were co-varying between the minimum and the maximum value (Table 4) and the rest of the PROSAIL parameters (Table 3) were kept constant at the average value. This corresponds to a scenario where all canopy and leaf structural traits change during the growing season (due to, e.g., ecological factors such as temperature, soil moisture and light competition), but biochemical traits remain constant. This scenario can be associated with grasslands where phenology, species composition and management practices (e.g., architectural changes due to variation of species composition after mowing and leaf biochemical changes due to fertilization) are determining more evident both canopy architectural and leaf structural traits dynamics at the temporal scale.
(4t) All the PROSAIL input variables (Table 4) were allowed to co-vary from the minimum to the maximum value. This corresponds to a scenario where there is a relevant temporal co-variation of both canopy structural and biochemical traits. Such scenario can be considered representative of grasslands with more extreme phenology dynamics, due to: (i) stronger changes of ecological factors related to climate dynamics and water availability (e.g., in Mediterranean grasslands with strong seasonal variations of leaf water content, leaf chlorophyll content and proportion of brown dead leaves); or (ii) management practices (e.g., architectural changes due to variation of species composition after mowing and leaf biochemical changes due to fertilization).
The PROSAIL model was also run with a 100-iteration step for four scenarios (1s-4s) using the range of LAI measured during spatial observations (minimum and maximum values measured in the field at the biomass peak) as an input parameter (Table 4) as follows: (1s) LAI was set to vary between the minimum (in the less productive grassland) and the maximum value (in the most productive grassland) and other PROSAIL parameters were kept constant at their average value. This corresponds to a theoretical scenario where grassland is spatially homogeneous in terms of both structure and biochemistry.
(2s) LAI and LAD values were set to vary between the minimum and the maximum value and other PROSAIL parameters (Table 4) were kept constant (average value). This scenario can be considered representative of grasslands with slightly spatially-heterogeneous species composition and functional types.
(3s) LAI, LAD, N and Cm values were co-varying between the minimum and the maximum value and other biochemical parameters were kept constant. This scenario corresponds to a grassland with strong structural spatial heterogeneity, while biochemistry is homogeneous. It can be associated with grasslands with heterogeneous species composition and functional types, associated with spatial variation of ecological factors (e.g., soil moisture and pH, topographic aspect, etc.).
(4s) All the PROSAIL parameters were co-varying. This scenario corresponds to a grassland with strong structural and biochemistry spatial heterogeneity and can be considered representative of grasslands with extreme variations of species composition, plant functional types, ecological factors and management regimes.
All PROSAIL simulations were performed using the MATLAB environment (The MathWorks Inc. 2019a).

Statistical Analysis
To compare the performance of the investigated NDIs in LAI estimation, the following linear and second-order polynomial regression statistics were computed: R 2 , coefficient of determination; Adj. R 2 , adjusted coefficient of determination; and RMSE, root mean square error. All statistical analyses were performed by means of the R software (version 3.6.0, https://www.r-project.org/).

Relationship between the Measured Spectra and LAI
Temporal trends of LAI in IT-MBo (2013 and 2014) and in AT-Neu (2018) are shown in Figure 3A-C, respectively. Figure 3D presents LAI corresponding to spatial observations. For temporal scale observations, the values of LAI at IT-MBo ranged from 0.3 to 3.7, while, at AT-Neu, LAI ranged from 0.2 to 8.8 and showed a smoother LAI increase compared to IT-MBo. At the spatial scale in IT-MBo, as the measurements were carried out at the biomass peak, the variability was lower and the LAI ranged from 2.0 to 4.5.
The canopy reflectance obtained from field measurements performed at different vegetation growth stages is plotted in Figure 3A Figure 3D. In the case of temporal observations, canopy reflectance showed a typical response to LAI increase, mostly characterized by a gradually decreasing reflectance in the VIS region, a gradually increasing reflectance in the NIR-shoulder region and increasing NIR-shoulder slope as observed by Vescovo et al. [36]. On the contrary, the reflectance values of spatial scale observations at IT-MBo 2017 showed a much more complex pattern ( Figure 3D), where increasing values of LAI did not always result in higher NIR and lower visible reflectance. As an example, the grassland plot highlighted in yellow (with a relatively high LAI value of 3.4) had a very low NIR reflectance, while the grassland plot highlighted in red (LAI 4.5, corresponding to the highest LAI value) showed the highest reflectance in the visible wavelengths.
The reflectance values at 740 nm appeared to be relatively closer to the ones of the NIR shoulder plateau at AT-Neu 2018 compared to IT-MBo 2013, IT-MBo 2014 and IT-MBo 2017, suggesting a lower chlorophyll absorption at this wavelength in this ecosystem. The NIR-shoulder slope, between 760 and 900 nm, observed in AT-Neu 2018 was generally slightly less steep compared to the slope observed in IT-MBo in both temporal and spatial scale observations.   Figure 4A,B), the R 2 values displayed very consistent patterns across VIS-NIR band combinations, but with slightly different R 2 ranges for both investigated years: the maximum R 2 was around 0.8 for IT-MBo 2013 and 0.9 for IT-MBo 2014 and the minimum R 2 was slightly higher for IT-MBo 2013 compared to IT-MBo 2014. At the AT-Neu study site, generally lower R 2 values were observed compared to the IT-MBo temporal observations. In the VIS spectral region, a slightly different pattern of the R 2 values was observed for both study sites ( Figure 4C). The R 2 patterns were more different across the RE and NIR-shoulder regions. In particular, an evident shift of the well-correlated areas towards the lower wavelengths (from around 750 nm to 740 nm) was observed for AT-Neu ( Figure 4C) compared to IT-MBo temporal observations ( Figure 4A,B).

Best Band Combination and Hyperspectral NDIs
All the two bands NDIs (  Figure 4A,B), except for the NSDI 865.783 index (R 2 < 0.3, Figure 4A,B). For AT-Neu, a slightly different pattern of R 2 values was observed, where RENDVI 865.740 and NSDI 779.754 showed very low correlations compared to IT-MBo (R 2 > 0.2, Figure 4C). The best correlations at the AT-Neu site were observed for the RENDVI NDIs (RENDVI 783.740 and RENDVI 783.705 ) with R 2 values exceeding 0.6 ( Figure 4A-C). Therefore, the RENDVI 783.740 , which is still not very commonly used in the literature, performed well at both sites. Although the NIR-shoulder NSDI 779.754 showed a significant correlation with high R 2 values (R 2 > 0.85 Figure 4A,B) at the IT-MBo site, much lower R 2 values were observed for the same band combinations at At-Neu (R 2 < 0.2, Figure 4C). Other than the aforementioned SVIs, there is a wide range of band combinations which showed a strong correlation with LAI. For both sites, an area of high R 2 values (R 2 > 0.7) was observed for combinations within the NIR part of the spectrum (B2 around 950-970 nm and B1 of 900-950 nm: Figure 4A-C). Such band combinations are commonly used to calculate water band SVIs which are indicators of water status [60][61][62][63][64]. On the other hand, water band SVIs demonstrated to be good proxies of structure-related parameters such as LAI and phytomass [10,65]. A strong correlation between Normalized Water Index (NWI; calculated with PRISMA bands B2 of 962 nm and B1 of 897 nm [62,63,66]) and LAI was observed (R 2 value for IT-MBo 2013 = 0.77, for IT-MBo 2014 = 0.83 and for AT-Neu 2018 = 0.85).
The correlogram based on S-2 ( Figure 4D-F) and S-3 ( Figure 4G-I) bands provides an overview of the performance of band combinations obtained with the average reflectance over the bandwidth of the respective S-2 and S-3 bands. For both study sites, the S-2 graph indicated that the RENDVI based on wavelength 740 and 783 nm is among the best correlated with LAI (R 2 > 0.65, Figure 4D-F). For S-3 simulated data, NSDI 779.754 showed very high correlations for IT-MBo (R 2 > 0.75, Figure 4G,H) but very low correlations at AT-Neu. Moreover, we also analyzed the correlograms of RMSE based on temporal observations ( Figure S4), showing that the RMSE patterns were mostly similar but inverse, as high R 2 values corresponded with low RMSE values and vice versa.
Results from the observations at the spatial scale at IT-MBo 2017 showed poorer correlations ( Figure 5A) than those obtained with the multi temporal data. Only band combinations from the spectral range within 350-400 nm showed slightly higher R 2 values compared to the rest of the band combinations. The correlograms of RMSE based on spatial observations presented in Figure S5 showed overall high RMSE for most of the two band combinations.

The Performance of Multispectral Sentinel 2 and 3 SVIs
Scatterplots between VIS-NIR SVIs (calculated with data resampled to S-2 bands) and LAI for both sites are presented in Figure 6. Both investigated VIS-NIR SVIs (NDVI 865.665 and MTCI) showed a high correlation with an R 2 > 0.75 and RMSE < 0.5 m 2 ·m −2 ( Figure 6A,B and Table 5) for observations on a temporal basis at IT-MBo. For the AT-Neu study site (which has very high productivity and LAI values up to 8.8), NDVI 865.665 showed a strong saturation effect ( Figure 6A) compared to MTCI, resulting in a lower R 2 value (R 2 = 0.55 and R 2 = 0.81 for NDVI 865.665 and MTCI, respectively) and RMSE (RMSE = 1.71 m 2 ·m −2 and RMSE = 1.11 m 2 ·m −2 , respectively). Conversely, when considering the observations at the spatial scale performed at IT-MBo, no significant correlation was found between the investigated VIS-NIR SVIs and LAI ( Figure 6C,D and Table 5).
Scatterplots between RE-based SVIs (calculated with data resampled to S-2 bands) and LAI for both sites are presented in Figure 7. The R 2 values of the linear relationship between RE SVIs and LAI in IT-MBo were high with R 2 > 0.8 (up to 0.96 and RMSE < 0.45 m 2 ·m −2 ) for the temporal observations at IT-MBo, but at the AT-Neu site the RE SVIs showed slightly lower R 2 (< 0.8 with RMSE > 1 m 2 ·m −2 ) values, as a result of the saturation effect ( Table 5). The saturation effect was particularly strong for NDVI 865.665 and for AT-Neu site, where a polynomial regression model showed to increase the R 2 from 0.55 to 0.71 and decrease RMSE from 1.71 to 1.37 m 2 ·m −2 . For all other SVIs, polynomial relationships showed only a slight increase of the performance of the model (Table S2). The RE SVIs including RENDVI 783.740 and RENDVI 783.705 from temporal scale observations showed a strong correlation with an R 2 > 0.65 for both study sites, but RENDVI 865.740 lost its predictive power at AT-Neu (R 2 = 0.2 and RMSE = 2.28 m 2 ·m −2 ). As for VIS-NIR SVIs, the observations at the spatial scale showed very low R 2 (< 0.06) values.  Concerning the NIR shoulder SVIs, the NSDI 779.754 showed a strong positive correlation with an R 2 > 0.8 (up to 0.95 and RMSE < 0.40 m 2 ·m −2 ) for the temporal observations at IT-MBo, but at the AT-Neu study site the correlation was much weaker (R 2 < 0.1, Figure 8A; RMSE = 2.44 m 2 ·m −2 , Table 5). The index based on 865 and 783 nm showed an inverse relationship with LAI for both ecosystems, with weaker correlation at IT-MBo (R 2 < 0.3) compared to AT-Neu (R 2 of 0.58, Figure 8B). No significant correlations were observed on the spatial basis ( Figure 8C,D).
Considering all the VIS-NIR, RE and NIR shoulder SVIs, RENDVI 783.740 was always among the three best-performing SVIs, for observations carried out at the temporal scale at the two investigated sites ( Table 5). The other SVIs with good performance at the temporal scale were RENDVI 865.740 , MTCI and NSDI 779.75 . In all panels, solid lines represent a linear fit to the data. Table 5. Summary of the statistics (N, Number of observations; R 2 , coefficient of determination; Adj. R 2 , adjusted coefficient of determination; RMSE, root mean square error) of the linear regression between leaf area index (LAI, m 2 ·m −2 ) estimated from fraction of absorbed photosynthetically active radiation (fAPAR) and the spectral vegetation indices (SVIs) calculated from measured spectra for IT-MBo 2013, IT-MBo 2014 and AT-Neu 2018 at temporal scale observations and for IT-MBo 2017 at spatial scale observations. The three best-fitting models are highlighted in bold. Asterisk indicates significance of correlation: *** p < 0.001; ** p < 0.01; * p < 0.05. n.s., not significant (Pearson's correlation test).

Global Sensitivity Analysis of the Spectral Bands
The results concerning the impact of leaf and canopy parameters on the different spectral regions through the GSA are illustrated in Figure 9. Background (soil) demonstrated to have a rather homogeneous impact on canopy reflectance across all wavelengths. LAI and Cab demonstrated a major influence on reflectance in the VIS part of the spectrum and Cab showed two SI peaks (around 560 and 705 nm) and an impact on the reflectance response was observed up to 760 nm. Within the RE (680-750 nm) part of the spectrum, the influence of LAD was significantly increasing at longer wavelengths, with SI reaching about 40% at 740 nm, while at 705 nm the SI was less than 10%. For spectral bands in the RE region, a slight effect of Cbrown and Cm was also observed, while the impact of LAI was significantly increasing from RE to NIR-Shoulder spectral bands. LAI and LAD were the main drivers of reflectance also in the far RE region (740-750 nm) of the spectrum. Within the NIR-shoulder region (750-900 nm), reflectance showed to be driven mostly by LAI and LAD, while Cm showed an average SI value of less than 5%. In the spectral range 740-820 nm, Cbrown also showed SI values up to 5%. Leaf water content (Cw) response only started beyond 930 nm (Figure 9).  Table 3.

SVIs Calculated from Modeled PROSAIL Reflectance Simulations vs. LAI
To analyze the impacts of co-variation of structural and biochemical traits for LAI monitoring using the PROSAIL modeled reflectance output, a series of scatterplots between the SVIs and LAI is presented in this section. The modeled reflectance was obtained from PROSAIL simulations constraining the input parameters according to the eight scenarios presented in Section 2.3.4.
For the investigated grasslands at IT-MBo, the LAI range observed at the temporal scale also included low LAI values (0.3), while the minimum LAI value for spatial scale observations was 2.0. In general, when LAI values were restricted to the spatial range, the lack of low LAI had a strong effect on the predictive power of VIS-NIR SVIs with a noticeable decrease in correlation coefficients when two or more traits were co-varying ( Figure 10B-D F for NDVI 865.665 and Figure 10F-H for MTCI) and increase in RSME (Tables 6 and 7).

Discussion
The grassland spectral response across different spectral regions (VIS, RE and NIR-shoulder) showed to be both site-specific and scale-dependent. The NIR-shoulder slope showed different trends at two sites confirming the results of Vescovo et al. [36] on the site-specificity of NIR-shoulder indices. Moreover, the NIR-shoulder slope response at the spatial scale appeared to be more complex than the one at the temporal scale and did not strictly follow the typical temporal response at increasing LAI, characterized by an increase of NIR-shoulder slope corresponding to an increase of LAI.
The importance of using hyperspectral and superspectral sensors adopting SVIs with bands in the RE and NIR shoulder spectral region [19,36] to estimate canopy structure was highlighted in this study for montane temperate grassland ecosystems, and novel band combinations were explored. RTM observations demonstrated that chlorophyll absorption can determine the spectral response well beyond 700 nm, and thus the optimal band combinations to estimate LAI often included far red bands (as shown also in [32,33,67]). Ollinger's paradox [5] stating that (as chlorophyll absorption range was assumed by this author to be 400-700 nm) "the physiological activity of vegetation is often more strongly related to reflectance at wavelengths that are not used in photosynthesis than to those that are" has been-at least partially-explained during the last years as many studies have been highlighting that photosynthetic pigments' absorption is active in the far red domain. In our study, the RTM models indicate that chlorophyll absorption can be observed up to 760 nm (Figure 9), while in vivo observations previously indicated 740 nm as an upper limit [67]. For both sites, high R 2 values were detected for band combinations used in water indices. Such result is opening interesting perspectives as these indices can be calculated using both PRISMA and GF-5 satellite data.
The GSA showed that the VIS part of the spectrum is mainly influenced by LAI and Cab, in the RE spectral region reflectance is determined by LAI, Cab, LAD, Cm and Cbrown. In previous studies the RE spectral region showed optimal performances in retrieving LAI and canopy chlorophyll content. The RE part of the spectrum is characterized by lower absorption by chlorophyll, but remains sensitive to changes in its content, reducing the saturation effect and enhancing the sensitivity of these SVIs to moderate-high vegetation densities reducing the typical SVI saturation effect [32,41]. The evident shift of the well-correlated areas towards the lower wavelengths at the AT-Neu site suggests that the performance of some of the RE SVIs is site-specific, probably partly due to different absorption thresholds [29,36]. It is interesting to notice that, despite such shift, the area corresponding to RENDVI 783.740 was well correlated at both sites, and for both the two years of observations at IT-MBo.
Very contrasting results were achieved at the temporal and spatial scales. At IT-MBo temporal scale strong correlations (R 2 > 0.8) were observed between LAI and both traditional RE and NIR-shoulder (NSDI 779.754 ) SVIs. Differently from previous studies [20,36,37], the performance of such indices at the spatial scale was particularly poor with an R 2 < 0.1. Such poor performance could be partially explained by the different LAI ranges at the temporal and spatial scales. At the spatial scale, saturation of some SVIs may be observed above certain LAI values, and this can constrain the ability to retrieve LAI when only full-canopy cover ecosystems are observed [10]. The poor performance of the SVIs in retrieving LAI at the spatial scale observations is noteworthy and confirming the observations of Dong et al. [68] on the strong response of canopy reflectance to canopy structural traits. Darvishzadeh et al. [46] showed that LAI, in heterogeneous grasslands, could be estimated at the spatial scale using the SVIs approach with intermediate accuracy (R 2 cv values from 0.49 to 0.69). Atzberger et al. [25] evaluated the PROSAIL RTM suitability for grasslands and demonstrated that PROSAIL is well suited for LAI estimations. However, in their study, the PROSAIL-generated and the in-situ hyperspectral-derived correlation plots, across the RE and NIR-shoulder ranges, demonstrated different R 2 patterns between observed and modeled LAI, indicating that RTM parameterization for LAI retrieval is challenging in these spectral domains. The selection of leaf structural and canopy architectural settings is key to achieve an accurate LAI retrieval, and can be complex when heterogeneous canopies are modeled. Following the indications of Darvishzadeh et al. [46] on the impact of grassland species composition and canopy architecture on remote sensing models, we hypothesized that the structural heterogeneity of the investigated montane temperate grasslands at the spatial scale might play a crucial role on LAI retrieval. Our results confirm that: (i) the co-variation of all structural traits (such as LAI, LAD, Cm and N, at the spatial scale) could explain the poor performance of most SVIs; and (ii) due the co-variation of both structural and biochemistry traits, no SVI is able to provide reliable LAI spatial estimations. In our study, we used the PROSAIL RTM in forward mode to study the impact of co-variation of PTs on LAI estimation. The PROSAIL results confirmed that trait-covariation resulting from extreme grassland heterogeneity has a strong impact on LAI estimation accuracy. Such findings agree with the observations of Ollinger [5] on "the difficulty of assessing the relative importance of individual traits that co-vary with a suite of plant properties". For the investigated montane temperate grasslands types-characterized by extreme spatial heterogeneity-structural and biochemical intraspecific drivers linked to heterogeneous canopy species composition seem to have a stronger impact on traits estimation than interspecific drivers related to phenology. In other ecosystem types (e.g., arid and Mediterranean grasslands characterized by strong seasonality and extreme phenology dynamics due to changes of ecological factors related to climate and water availability) observations at the temporal scale are more challenging and significant limitations for remote sensing analysis are posed. In such grasslands, senescence can take place at varying rates and periods, increasing the variability of surface biophysical and optical properties [69].
RENDVI 783.740 showed a good performance and it is not very commonly used in the literature, as it was introduced relatively recently. Peng et al. [41] demonstrated that RENDVI 783.740 was accurate in estimating canopy chlorophyll in crops with contrasting architectures (maize and soybean). Our results show how RENDVI 783.740 can be used to monitor grassland LAI. RENDVI 783.740 demonstrated to be the most insensitive to grassland structural traits co-variation. For this reason, considering the increasing availability of hyperspectral and superspectral sensors such as S-2, GF-5 and PRISMA, more studies are needed to investigate its full potential for monitoring grasslands and other spatially-heterogeneous ecosystems.
The impact of structural PTs on the relationships between SVIs and LAI, although is well known in the literature, should be taken more carefully into account [70]. Structural traits directly determine the interactions between light and both leaf and canopy elements. In particular, LAD (which is a key canopy trait whose effect has usually not been considered when applying common vegetation SVIs for mapping LAI) showed to have a strong influence for wavelengths > 705 nm. The impact of LAD has been determined for crops [71], but limited research was carried out for grasslands [36].
Grassland canopy structure heterogeneity may impact the applicability of algorithms to detect vegetation changes due to phenology also at the temporal scale [72,73], and this could be the case of temperate grasslands where species composition and coverage is varying due to, e.g., light competition dynamics. As the co-variation of biochemical traits was included in the PROSAIL RTM simulations, a significant impact on the models accuracy was observed, indicating that both structural and biochemical factors play a major role in models' performance. Biochemical heterogeneity can be determined not only by different vegetation types at the spatial scale, but also by phenological changes. In Mediterranean grasslands, where severe droughts and grassland curing are taking place during the summer months, both biochemical and structural changes are expected [74,75].
Using a combination of two NIR bands (both beyond 750 nm), the NIR-shoulder SVI NSDI 779.75 (calculated from S-3 simulated bands) performed very well at IT-MBo, but not at the AT-Neu. This is probably due to the fact that chlorophyll absorption threshold is different at the two grassland ecosystems, and absorption is still present at 754 nm only in the IT-MBo grassland. This result is defining the green-dependency of NIR shoulder SVIs investigated in Vescovo et al. [36], which were thought to be related to scattering mechanisms and not chlorophyll absorption. Our study highlighted that this spectral region, which has been poorly investigated in the literature, is largely affected, in addition to chlorophyll until 760 nm, by LAI, LAD, Cbrown and Cw.
The suitability of well-known and widely adopted SVIs for retrieving LAI in grasslands with heterogeneous structure was also questioned in this paper. Many widely-adopted SVIs, e.g., NDVI 865.665 and MTCI, exhibited a strong correlation with LAI when only a few traits were co-varying, while a much weaker correlation was observed when more traits were co-varied. This agrees with the findings of Peng et al. [41] and Horler et al. [22] and advises to carefully evaluate potential uncertainties of satellite-based vegetation products such as LAI, fAPAR in spatially-heterogeneous canopies. In our work, we demonstrated that SVIs such as RENDVI 783.740 seem to be less influenced by canopy architecture, leaf structure and biochemical traits co-variation, and need further testing.
In this work, some constraints were highlighted on the use of statistical approaches based on SVIs. However, strong limitations of RTM inversion are also implied, and a reliable LAI spatialization in heterogeneous canopies needs to be based, in the future, on more detailed parameterization of the traits which are co-varying with LAI. RTM models show evident intrinsic limitations in their capacity to simulate heterogeneous canopies [76,77], and do not take fully into account some parameters, e.g., the presence of non-photosynthetic material in the canopy [58]. Furthermore, the GSA showed the impact of many different input parameters on LAI retrieval, and only very preliminary information on the spectral impact of some RTM parameters (e.g., canopy brown pigments and anthocyanins) is available in the literature [76,78].

Conclusions
The potential of Sentinel bands combinations across the RE and the NIR-shoulder spectral region such as RENDVI 783.740 and NSDI 779.754 (which are novel or not commonly used in the literature) was highlighted. Such SVIs are worth more attention to ascertain their performance on other canopy types. Moreover, the hyperspectral analysis highlighted the suitability of the spectral regions related to water absorption features for LAI estimations [10,65].
The impact of grassland structural and biochemical heterogeneity on LAI estimations was demonstrated to be strong and, for this reason, no reliable field LAI estimation was possible at the spatial scale with any SVI. The results of the empirical approach were confirmed by the simulations performed with the RTM PROSAIL, when both structural and biochemical traits were co-varied. In this context, the uncertainties of satellite-based LAI products (in grassland canopies with either spatially or temporally-heterogeneous structure) need to be carefully taken into account adopting a modeling approach which is minimizing the impact of canopy structural heterogeneity. Despite the fact that the sensitivity analysis demonstrated that LAD impact is quite strong starting from 705 nm, RENDVI 783.740 proved to be the best performing S-2-based SVIs for monitoring grasslands with heterogeneous structure. Given the fact that our study was carried out in two sites in the Alps, and that spatial observations were carried out only in a limited number of plots, more studies are needed in other grassland ecosystems and/or in other geographic areas to confirm the potential of SVIs using this spectral domain (alongside with the water absorption features) for vegetation monitoring, in the context of the Sentinel, GF-5 and PRISMA missions.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/14/2254/s1, Figure S1: RGB images of the different plots at the IT-MBo study site used for spatial analysis. Figure S2: RGB images of the plot at the AT-Neu study site used for temporal analysis. Figure S3: The portable system used for spectral measurements at the IT-MBo for spatial scale observations. Table S1: Specifications of the multispectral instrument (MSI) and ocean and land color instrument (OLCI) on the S-2 and S-3 satellite system, respectively. The NIR-shoulder bands investigated in this study are shown in bold. Table S2: Summary of the statistics (N, Number of observations; R 2 , Coefficient of determination; Adj. R 2 , adjusted coefficient of determination; RMSE, Root mean square error) of the second-order polynomial regression between leaf area index (LAI, m 2 ·m −2 ) estimated from fraction of absorbed photosynthetically active radiation (fAPAR) and the spectral vegetation indices (SVIs) calculated from measured spectra for IT-MBo 2013, IT-MBo 2014 and AT-Neu 2018 at temporal scale observations and for IT-MBo 2017 at spatial scale observations. The three best-fitting models are highlighted in bold. Asterisk indicates significance of correlation: *** p < 0.001; ** p < 0.01; * p < 0.05. n.s., not significant (Pearson's correlation test).