Improving the Performance of 3-D Radiative Transfer Model FLIGHT to Simulate Optical Properties of a Tree-Grass Ecosystem

Abstract

The 3-D Radiative Transfer Model (RTM) FLIGHT can represent scattering in open forest or savannas featuring underlying bare soils. However, FLIGHT might not be suitable for multilayered tree-grass ecosystems (TGE), where a grass understory can dominate the reflectance factor (RF) dynamics due to strong seasonal variability and low tree fractional cover. To address this issue, we coupled FLIGHT with the 1-D RTM PROSAIL. The model is evaluated against spectral observations of proximal and remote sensing sensors: the ASD Fieldspec^® 3 spectroradiometer, the Airborne Spectrographic Imager (CASI) and the MultiSpectral Instrument (MSI) onboard Sentinel-2. We tested the capability of both PROSAIL and PROSAIL+FLIGHT to reproduce the variability of different phenological stages determined by 16-year time series analysis of Moderate Resolution Imaging Spectroradiometer-Normalized Difference Vegetation Index (MODIS-NDVI). Then, we combined concomitant observations of biophysical variables and RF to test the capability of the models to reproduce observed RF. PROSAIL achieved a Relative Root Mean Square Error (RRMSE) between 6% to 32% at proximal sensing scale. PROSAIL+FLIGHT RRMSE ranged between 7% to 31% at remote sensing scales. RRMSE increased in periods when large fractions of standing dead material mixed with emergent green grasses —especially in autumn—; suggesting that the model cannot represent the spectral features of this material. PROSAIL+FLIGHT improves RF simulation especially in summer and at mid-high view angles.

1. Introduction

Structure and biochemical constituents govern the optical responses of the vegetation [1]. This fact allows remote sensing to provide spatial and temporal information on vegetation properties, which is necessary to understand key processes and variables such as photosynthesis or health status of the plants [2,3,4]. Remote sensors provide spectro-directional reflectance factors (RF) at the top of the canopy that can inform on both vegetation structure and leaf properties. Linking these variables is not trivial since the problem is ill-posed and other factors, such as woody elements, background scattering, illumination conditions, or the viewing geometry confound their relationship [5,6]. Radiative Transfer Models (RTM), spectral unmixing, and data fusion are the main ways to address these issues in remote sensing [7,8]. RTM are physical models capable of simulating the interaction of light with vegetation at leaf and canopy levels. These models are based on the radiative transfer equation that studies the propagation of radiation through an absorptive, dispersive, and/or emissive medium. The most frequently used RTM describes the vegetation canopy as a turbid medium where infinitesimally small leaves and gaps are distributed randomly [9]. These models assume that leaves are bi-Lambertian [10] and that Leaf Area Index (LAI) and Leaf Angle Distribution (LAD) mainly control the architecture of the canopy [11]. This 1-D approach can acceptably represent homogeneous vegetation canopies like most of the crops, grasslands, and closed forests; but might not hold for heterogeneous and discontinuous vegetation canopies [12,13]. More complex RTM provide a 3-D representation of the vegetation that is more appropriate for these ecosystems [7]. 3-D RTM overcome traditional approaches in radiative transfer based on the Kubelka–Munk theory [14,15,16] and use more suitable methods to perform the spatial and spectral scaling [7]. Among the most common approaches is the Monte Carlo Ray Tracing method that simulates a chain of dispersion events experienced by a photon in its trajectory from the source to the receiver or until its absorption [17,18]. Monte Carlo models allow accurate representation of the canopy, especially, when a large number of photons—usually over 10⁶ photons—are processed, which results in large computational costs and processing times [19]. In recent years, increasing computing capabilities have expanded the use of Monte-Carlo-based RTM to a number of applications [20,21,22,23]. For instance, Kötz et al. 2004 [13] demonstrated the possibility of estimating live fuel moisture and green fuel loading in a coniferous forest with 3-D RTM FLIGHT [17]. Simulated spectra compared to images of the airborne Reflective Optics System Imaging Spectrometer (ROSIS) and Hyperion Satellite Imaging Spectrometer obtained a relative deviation of around 5% in all wavelengths. Errors found were related to uncertainties in the variables describing the canopy structure such as LAI or the fractional Vegetation Cover (F_Cover).

3-D RTM usually involve a large number of parameters to describe, in detail, the vertical and horizontal distribution of the canopy, which might require field data to accurately simulate the directional reflectance of a given pixel. For example, Schneider et al. [24] collected substantial field data in order to accurately simulate a temperate mixed forest with a Discrete Anisotropic Radiative Transfer (DART) model [25,26]. Such simulations were successfully compared with images acquired with the Airborne Prism Experiment (APEX) sensor. The study analyzed the sources of error and concluded that these were partly related to an inadequate representation of small-scale structures and the variability of leaf optical properties between species and individual trees. They also carried out an analysis of the background influence on simulated RF by replacing a background based on the spectral properties of the main species in the study area by a black 100% absorbing background. Results showed an impact that was higher than the authors expected for a relatively dense mixed forest. The impact should be much higher in ecosystems with low tree fractional cover and high clumping, such as tree-grass ecosystems (TGE).

Mixed tree-grass and shrub-grass vegetation associations are one of the most spatially extensive and widely distributed forms of terrestrial vegetation on Earth and are often linked to livestock production [27]. These ecosystems experience great pressure from changes in land use, at the same time that they are expected to increase in surface area and shift their geographical distributions as a consequence of climate change [28]. The expected up-coming changes and increased role of TGE demand the improvement of our understanding of their functioning and our capability of monitoring their status from space.

TGE combine an overstory of evergreen or deciduous trees, sometimes a middle layer of shrubs, and a grass understory which—in some bioclimatic regions as in the Mediterranean—, can present high biodiversity as well as large spatial variability and strong temporal dynamics. Several approaches have been applied to simulate the optical properties of inhomogeneous ecosystems using physical or empirical methods [25,29,30]. A number of 3-D RTM have been developed to model complex bi-layered ecosystems such as the Forest Radiative Transfer(FRT) model [31] or the 5-scale model [32]. These complex models are the result of an ensemble of RTM operating at different scales (leaf, canopy, different layers…). In addition, vegetation RTM are sometimes coupled to atmospheric RTM in order to simulate realistic illumination conditions [31] or retrieve vegetation parameters from top of the atmosphere radiances observed from remote sensors [33]. Several of these models have been tested in the third phase of the RAdiation Transfer Model Intercomparison (RAMI) project (http://rami-benchmark.jrc.ec.europa.eu) [34].

FLIGHT is one of the 3-D models tested in the third phase of the RAMI experiment [34]. It presents some features that make it suitable to model vegetation in Mediterranean environments. On one side, FLIGHT is not directly coupled to a leaf RTM. The model demands leaf reflectance and transmittance factors, which can be provided externally by models or databases. This feature is convenient in those cases where RTM cannot accurately represent leaf optical properties of certain species since measurements can be directly provided. For example, in Mediterranean ecosystems, sclerophyllous species such as Quercus ilex L. can develop different protective structures on leaf surfaces. Pacheco-Labrador et al., [35] found that trichomes can strongly modify leaf optical properties, hampering the inversion of leaf RTM PROSPECT to estimate leaf biophysical properties. On the other hand, the FLIGHT model allows for the representation of two different types of leaves. This is suitable to represent canopies combining leaves of different ages or in different stages of development. For example, in Mediterranean ecosystems, evergreen species can present annual regrowth of new leaves which are very different from the older leaves remaining in the crown from previous years. A number of studies have used FLIGHT to simulate canopy reflectance in different ecosystems [20,21,23]. Most of these works mainly focused on the representation of RF as a function of tree crowns’ properties standing on bare soil. However, in TGE the background is also composed to vegetation and can represent high fractions of the surface. Therefore, an adequate representation of the radiative transfer of this layer is necessary to mechanistically link vegetation properties to the RF in these ecosystems. This study aims to improve the capability of FLIGHT to simulate the spectral RF of TGE with a grass understory by coupling this model with PROSAIL [36] (hereafter PROSAIL+FLIGHT). In this model, the 3-D FLIGHT model represents tree crowns whereas the 1-D RTM PROSAIL features the underlying grassland.

PROSAIL+FLIGHT is tested in a Mediterranean TGE against RF measurements at proximal and remote sensing scales. First, we tested the capability of the model to represent the ranges and the seasonal variability of spectral RF measured in the study area. To do so, the main phenological phases (hereafter, phenophases) were characterized using Moderate Resolution Imaging Spectroradiometer (MODIS) time series analysis in the study site (Section 2.1). This allowed us to limit the ranges and define the covariance between the biophysical parameters of vegetation in each phenophase according to field measurements. We used this information to generate Look Up Tables (LUT) of RF characteristics of each of these periods with PROSAIL and PROSAIL+FLIGHT; the RTM represented the herbaceous background and the ecosystem, respectively. Secondly, field biophysical data were used to simulate RF with PROSAIL and PROSAIL+FLIGHT, which were compared with spectral measurements acquired in the same spots at different spatial scales (ground, airborne, and satellite). Finally, we compared RF predicted by PROSAIL, FLIGHT, and PROSAIL+FLIGHT at different angles in the principal and the cross-principal plane.

2. Materials and Methods

2.1. Study Site

The study site is a well-monitored experimental station included in the FLUXNET network (http://fluxnet.ornl.gov/site/440) located at Majadas de Tiétar, in the West of Spain (39°56′29″N, 5°46′24″W), see Figure 1. It is a Mediterranean TGE—locally known as “dehesa”—with an extensive livestock use (<0.3 cows ha⁻¹). The climate is Continental Mediterranean with a mean annual temperature of 16 °C and a mean annual precipitation of 650 mm. Rain events are concentrated especially in the autumn and spring seasons [37].

Vegetation is structured in two layers; a tree stratum and a grass understory, see Figure 1b–c. The tree layer covers ~20% of the surface and is primarily composed of Holm oak (Quercus ilex, subsp. ballota L.), an evergreen species that undergoes an annual process of partial leaf renewal every spring. As a consequence, leaves of three vegetative periods may coexist in the tree crown. During the leaf expansion process in spring, the newly sprouted leaves present biochemical and spectral features that differ greatly from the leaves sprouted in previous years [38]. In spring, the outer canopy layer is composed of a dynamic changing mixture of current and the previous year’s leaves. The biochemical and spectral behavior of the canopy achieve a homogeneous distribution pattern in autumn and winter time. The ratio between old and new leaves is driven by the availability of water and nutrients, as well as the presence of diseases and plagues [39]; this changes throughout the year and can be also quite dynamic in space. The grass layer covers almost the totality of the soil, even under the tree canopy, and combines different functional types: grasses, forbs, and legumes [40]. This layer is highly dynamic [41] and biodiverse in space and time, depending on the availability of resources [42,43,44]. Its strong temporal dynamics are characterized by two seasonal biomass peaks with a maximum in spring and a second maximum after the summer decay, following the soil re-wetting phase. Besides, the grass cover, this Mediterranean TGE features large fractions of dry biomass accumulated during the grass regrowth and decay processes [41]. The grassland layer dominates the optical signals and carbon fluxes of the ecosystem [45].

2.2. PROSAIL and FLIGHT Input Parameters

This work combines 3-D RTM FLIGHT with 1-D RTM PROSAIL to simulate directional reflectance factors in a Mediterranean TGE. PROSAIL combines the leaf model PROSPECT [46] and the canopy model SAIL [16]. PROSPECT is a leaf RTM, which simulates optical directional-hemispherical reflectance (ρ_leaf) and transmittance (τ_leaf) factors as a function of the leaf internal structure and biochemical parameters in a spectral range from 400 to 2500 nm [46]. The model is based on Allen’s plate theory [14]. This study runs PROSPECT-5B [36] where Chlorophyll a and b (C_ab) and Carotenoids (C_ar) are separately included as model inputs, see

Table 1. Inputs required by the PROSPECT, SAIL (for the grass layer), and FLIGHT (for the tree layer) models. SAIL and FLIGHT models share the underlined variables.

Variable	Symbol	Units	Source (In Tis Work)
PROSPECT-5B
Leaf structural parameter	N	Unitless	Literature
Leaf water content	Cw	mg/cm2	Destructive sampling
Dry matter content	Cm	mg/cm2	Destructive sampling
Chlorophyll content	Cab	µg/cm2	Destructive sampling
Carotenoids content	Car	µg/cm2	Destructive sampling
Brown pigments	Cbrown	Unitless	Literature
SAIL
Leaf Area Index	LAIgrass	m2/m2	Destructive sampling
Leaf Inclination Distribution Function or Leaf Angle Distribution 1	LIDF/LAD	(°)	Literature and field observation
Hot Spot	q	m/m	Literature
Solar Zenith Angle	θs	(°)	Function of azimuth angle and date
View Zenith Angle	θv	(°)	Nadiral observation
Relative azimuth angle between sun and sensor	φsv	(°)	Solar calculator and satellite overpass geometry
Soil Reflectance	ρsoil	%	Field spectroscopy
FLIGHT
Leaf Area Index	LAItree	m2/m2	Tree sampling
Leaf Size	LFsize	cm	Destructive sampling
Leaf Area Distribution	LAD	(°)	Not measured
Fraction of Vegetation Cover	FCover	%	airborne image
Fraction of Current-year leaves	Fcurrent	%	Destructive sampling
Fraction of Previous-year leaves	Fprevious	%	Destructive sampling
Fraction of shoots	FS	%	Destructive sampling
Solar Zenith Angle	θs	(°)	Solar calculator
Solar Azimuth Angle	φs	(°)	Solar calculator
View Zenith Angle	θv	(°)	Nadiral observation
View Azimuth Angle	φv	(°)	Satellite overpass
Aerosol Optical Depth	AOT	Unitless	Field measurements and Sentinel-2 product
Crown Shape	Cs	Unitless	Airborne Laser Scanning
Crown radious	Exy	m	Airborne Laser Scanning
Centre to top distance	Ez	m	Airborne Laser Scanning
Height to first branch	Ht	m	Airborne Laser Scanning
Trunk Diameter at Breast Height	DBH	m	Field measurements
Tree leaf reflectance	ρleaf	%	Field spectroscopy
Tree leaf transmittance	τleaf	%	Field spectroscopy
Soil Reflectance	ρsoil	%	Field spectroscopy
Bark reflectance	ρbark	%	Field spectroscopy

¹ Note that the LAD represents the distribution of leaf zenith angles within the canopy [11], which can be typified in several distributions: Planophile, Erectophile, Spherical, Pagliophile, Extremophile and Uniform; while the LIDF is characterized by two parameters: (a) LIDFa, which controls the average leaf slope (from −1 to 1) and (b) LIDFb which controls the distribution’s bimodality (from −1 to 1).

. SAIL [16] is a 1-D canopy model based on the Kubelka–Munk theory and on the model developed by Suits [15]. It simulates directional RF of a canopy turbid medium by using a relatively small amount of inputs, see Table 1.

The FLIGHT model uses the Monte Carlo ray tracing method to predict the directional RF of vegetation canopies by simulating photon transport within a pixel or a scene. It represents tree crowns as geometric volumes randomly or deterministically distributed and a simple growth model in order to limit the overlap between neighboring crowns [22]. In the model, probability density functions determine the photon paths and fate when interacting with the crowns and the underlying background. These functions define the photon probability of scattering, absorption, or transmission [22,47]. The forest structure is simulated using a set of geometric primitives that allow for modeling of the main characteristic of the forest with reasonable computational resources.

2.2.1. Vegetation Biophysical Variables

PROSAIL and FLIGHT model inputs, shown in Table 1, were measured in the study site in 21 sampling campaigns between 2013 and 2017 in the context of the research projects FLUXPEC (http://www.lineas.cchs.csic.es/fluxpec/) and SynerTGE (http://www.lineas.cchs.csic.es/synertge/). PROSAIL variables were obtained by destructive sampling of the grass layer in 11 plots of 25 m × 25 m located in an area of around 1 km², see Figure 1b. Depending on the campaign, either two or three grass samples of 25 cm × 25 cm were collected in each plot in areas visually identified as representative of the sub-plot variability. In those plots containing trees, one of the sampling quadrants was located underneath one of the tree crowns in order to take into account the potential variability induced by tree–grass interactions. Spectral information was simultaneously acquired in the sampling plots (Section 2.2.2). All the grass samples were stored in sealed plastic bags, weighed in the field, and kept in a portable cooler with ice packs until they were stored in a freezer in the laboratory (see details in [48]).

Chlorophyll and carotenoids where measured only in samples collected during 2017. Pigment samples required a specific protocol for conservation and laboratory processing. Therefore, these samples were acquired in 25 cm × 25 cm quadrants contiguous to the ones were the samples for the rest of the biophysical parameters were collected. For pigment analysis, only the green fraction of the grass within the quadrant was collected. Samples were weighed in the field and immediately frozen in dry ice. Frozen samples were preserved at −20 °C at the laboratory until analysis was performed (see details in [49]).

Seasonal measurements of the tree canopy LAI were performed in 28 trees and nine field campaigns during 2015 and 2016. Direct LAI estimations were based on normalized surface canopy litterfall measurements of six vessels hanging in the tree crown and leaf removal rate measurements on eight external shoots per tree. In addition, seasonal measurements of the proportion of current and previous-year leaves and the shoots of the trees were performed from manual counts of branch samples located on the external part of the crown in north and south orientations. At the same time, measurements of the size of the leaves were made on the sampled trees. Airborne Laser Scanning (ALS) provided by the Spanish Program of Aerial Orthophotography (PNOA) was used to estimate the geometric structural parameters of the trees, shown in Table 1 (see details in [30]). Finally, trunk diameter at breast height (DBH) was obtained from field measurements made with a caliper located at a height of 1.30 m.

2.2.2. Ground Spectral Measurements

Simultaneous to the grass samplings, hemispherical conical reflectance factors (HCRF) were measured using an ASD FieldSpec^® 3 spectroradiometer (Analytical Spectral Devices Inc., Boulder, CO, USA) with a spectral range between 350 and 2500 nm and a nominal Field of View (FOV) of 25°. Around 30 spectra were acquired at regular intervals and at an approximate height of 1.20 m in two diagonal transects (NW-SE and NE-SW) in each of the 25 m × 25 m plots where the vegetation samples were collected. Before each transect, the dark current was recorded and subtracted. Down-welling radiance optimized the instrument parameters and was recorded at the beginning and the end of each transect on a calibrated 99% reflective Spectralon® panel (Labsphere Inc., North Sutton, NH, USA). All measurements were taken under clear sky conditions in a period of +/−2 h of solar noon. In 2017, spectral measurements were additionally acquired over 1 × 1 m sub-plots encompassing destructive grass sampling quadrants. Around 20 spectra per sub-plot were acquired. Dark current and down-welling radiance reference measurements followed the same protocol as that previously described for plot transects.

ρ_leaf and τ_leaf required in the FLIGHT model, see Table 1, were measured using a Li-Cor-1800-12 integrating sphere (Li-Cor, Lincoln, NE, USA) attached to the same ASD Fieldspec^® 3 spectroradiometer. Sampling took place at relevant phases of leaf development when significant differences in their biophysical and spectral properties occur. These periods were set at the beginning of spring when the new leaves sprout and at the end of the summer when the new leaves have almost achieved maturity and exhibit spectral behavior quite similar to previous year’s leaves. The integrating sphere measurement protocol is described in [50]. The operational spectral range of the sphere was limited to 300–1000 nm. In order to extend the spectral range to longer wavelengths, we characterized the RF of the internal coating using five calibrated Teflon® sheets of different widths up to 1600 nm. Above this wavelength, the signal was too noisy and the data was unusable [35]. The C_ab, C_ar, C_m, and C_w of leaves measured in the sphere were later determined from SPAD measurements (C_ab, C_ar) and gravimetric methods (C_m, C_w). ρ_leaf and τ_leaf in the interval 1600 to 2500 nm were predicted from inverting the PROSPECT model against the leaf RF within the operational range of the sphere [35].

Additionally, tree-bark RF required by FLIGHT was measured using the ASD Fieldspec^® 3 attached to an ASD Plant Probe device provided with a low-intensity bulb specially designed for collecting non-destructive data from vegetation and other heat-sensitive targets. The contact area of the probe and the tree trunk was sealed with black foam to prevent spurious entrance of ambient light due to a heterogeneous bark surface. Measurements were carried out in different areas of the trunk with different characteristics. Three trees were assessed in total and 25 spectra per tree were collected.

2.2.3. Airborne Hyperspectral and Satellite Images

Four Airborne campaigns were carried out on 08-April-2014, 23-May-2015, 03-July-2015, and 03-May-2016 using the CASI hyperspectral sensor operated by the Spanish National Institute of Aerospace Technology (INTA). Airborne overpasses were completed during approximately +/− 1 h of solar noon under clear sky conditions. The CASI-1500i sensor collects information in 144 bands in the visible (VIS) and near-infrared (NIR) regions with a full width at half maximum (FWHM) of approximately 7.5 nm and a spatial resolution of 1.5 m. Geometric correction was performed using onboard data and ground reference points while atmospheric correction was carried out using ATCOR4 (https://www.rese-apps.com/software/atcor-4-airborne/index.html) as described in [51]. Radiometric correction was further refined using the empirical line correction method [52,53] based on reference surfaces measured simultaneously to airborne overpasses using an ASD Fieldspec^® 3 spectroradiometer. A supervised classification was applied to the CASI images using the Mahalanobis distance method [54] in order to estimate F_Cover. For validation purposes, reflectance values of the pixels included in the sampling plots were extracted for all the bands of each image and average values per plot were calculated and further used to compare with simulated data.

Three MultiSpectral Instrument (MSI) L1C images acquired from Sentinel-2 on three different dates (21-February-2017, 22-April-2017, and 22-May-2017) were downloaded from the Sentinel Scientific Data Hub (https://scihub.copernicus.eu). Atmospheric correction was performed using the Sen2cor atmospheric correction algorithm, based on the Atmospheric Pre-corrected Differential Absorption (APDA) algorithm [55]. The Sen2cor tool was provided by the European Space Agency (ESA).

2.3. Generation of Seasonal Look-Up Tables of Reflectance Factors

This section describes the use of PROSAIL+FLIGHT to generate LUT of RF representative of the different phenophases identified in the study area. Phenophases determined ranges and covariance of RTM input parameters that were used to produce their respective LUT. A 16-year time series of the NBAR-MODIS product (MCD43A4) was used to analyze the spectral variability of the ecosystem in order to determine four different phenophases, mainly driven by the seasonality of the grass understory. Phenological transition dates in the study area were defined using the Normalized Difference Vegetation Index (NDVI) calculated from the MCD43A4 pixel centered in the study site. TIMESAT open software [56] was used to derive seasonal information from MODIS-NDVI data; start and end of the growing season as well as the seasonal amplitude. In addition, the maximum and minimum NDVI peaks were calculated each year. In order to remove extreme values and reduce noise, a Gaussian asymmetric model was applied to the whole time series [57]. The phenological year was classified into four phenophases: Phenophase (1) Summer drought, phenophase (2) Autumn regrowth, phenophase (3) Biomass peak; and phenophase (4) beginning of the grass decay season.

2.3.1. PROSAIL Simulation of the Grass Layer

For each of the phenophases, we produced a LUT containing PROSAIL parameters distributions, as well as bi-directional reflectance factors (BRF) and hemispherical-directional reflectance factors (HDRF) representing black and white sky illumination conditions. These RF were later used as background inputs of the FLIGHT model (Section 2.3.2). PROSAIL inputs were constrained according to the ranges and covariance of biophysical parameters observed in the field in each of the phenophases, see

Table 2. Minimum and maximum values of the PROSAIL input variables (units in Table 1).

Variable	Phenophase 1		Phenophase 2		Phenophase 3		Phenophase 4
	Min.	Max.	Min.	Max.	Min.	Max.	Min.	Max.
PROSPECT-5B
Mesophyll structure	1	1	1	1	1	1	1	1
Leaf water content	0.0	0.019	0.006	0.030	0.010	0.032	0.013	0.048
Dry matter content	0.0095	0.0182	0.003	0.009	0.003	0.009	0.005	0.013
Chlorophyll a+b content	1	10	12.6	63.4	12.6	63.4	12.6	63.4
Carotenoids content	0	1	3.7	12.7	3.7	12.7	3.7	12.7
Brown pigments	2	3	0	1	0	1	0	1
SAIL
Leaf Area Index	0.1	2.2	0.3	3.2	0.8	3.2	0.4	3.0
Hot Spot	0.001	0.001	0.001	0.001	0.001	0.001	0.001	0.001
Leaf Inclination Distribution Function (LIDFa)	−1	1	−1	1	−1	1	−1	1
Leaf Inclination Distribution Function (LIDFb)	−1	1	−1	1	−1	1	−1	1
Solar Zenith Angle	16.0	63.3	26.9	70.0	16.7	70.0	16.5	50.7
Observer Zenith Angle	0	0	0	0	0	0	0	0
Solar Azimuth Angle	114.1	225.2	125.3	212.2	114.2	232.0	111.5	229.9
Observer Azimuth Angle	105	105	105	105	105	105	105	105

. LUT generated using the Latin Hypercube Sampling (LHS) scheme [58] were used in order to achieve a uniform distribution of the model inputs within the boundaries found in each phenophase. The LSH technique divides the cumulative density function into n bins of the same size from which data is randomly selected. LHS enables a significant reduction in the number of simulations needed to completely map the desired space with respect to standard randomized sampling. We ran 1000 simulations per phenophase. Minimum and maximum values for all PROSAIL input parameters are shown in Table 2.

C_ab, C_ar, C_m, and C_w distributions were constrained by the observed covariance with LAI. This was achieved by fitting linear models where LAI was the independent variable and adding random noise sampled from a normal distribution of 0 mean and the standard deviation equals the Root Mean Squared Error (RMSE) [59] of the linear model fit. The structural parameter N was set to 1 assuming that grass species in the study site are predominantly monocotyledons [60]. The hot spot parameter was set to a fixed value used in previous works to model Mediterranean grasslands with PROSAIL [61]. The grass layer shows strong spatial and temporal biodiversity, with a varying combination of three functional types (grass, forbs, and legumes) which results in a variety of LIDF [40]. Therefore, LHS distributions of LIDFa and LIDFb were limited to the range [−1, 1] and forced to limit the sum of absolutes to 1. For a yearly period, a LUT containing day of the year (DoY), solar zenith (θ_s), and azimuth angles (φ_s) between 10 and 14 Coordinated Universal Time (UTC) with a 15-minute interval was generated. DoY and θ_s were obtained via LHS, and the corresponding φ_s were calculated via interpolation. Regarding the viewing angles, nadir observation and azimuth angle at the time of the Sentinel-2 overpass were set for the corresponding simulation.

The PROSAIL model requires information regarding the spectral properties of the soil. In this study, soil optical properties were estimated by constraining the brightness-shape-moisture (BSM) model included in the Soil Canopy Observation Photosynthesis Energy (SCOPE) model [62] against field spectral measurements of bare soil as described in [63]. This RTM simulates soil Lambertian reflectance as a function of a brightness factor, color coordinates, volumetric soil content, and capacity and water film optical thickness. Several soil measurements with different moisture conditions acquired in the study site with an ASD FieldSpec^® 3 determined the soil parameters via BSM model inversion. Soil Water Content (SWC) was parameterized as a function of vegetation C_w using a model fit with available observations from a network of Temporary Domain Reflectometry (TDR) sensors at a depth of 5 cm installed in the study site and the C_w estimated through destructive sampling. Then, SWC values estimated from the C_w values in the LHS were used to generate a range of soil spectra from very dry to wet conditions.

2.3.2. PROSAIL+FLIGHT Simulation of the Ecosystem

Table 3. Minimum and maximum values of the FLIGHT input parameters (units in Table 1).

Variable	Phenophase 1		Phenophase 2		Phenophase 3		Phenophase 4
	Min.	Max.	Min.	Max.	Min.	Max.	Min.	Max.
Solar Zenith Angle	16.0	63.3	26.9	70.0	16.7	70.0	16.5	50.7
Solar Azimuth Angle	114.1	225.2	125.3	212.2	114.2	232.0	111.5	229.9
View Zenith Angle	0	0	0	0	0	0	0	0
View Azimuth Angle	105	105	105	105	105	105	105	105
Aerosol Optical Depth	0.11	0.36	0.11	0.36	0.11	0.36	0.11	0.36
Leaf Size	2.16	3.36	2.16	3.36	2.16	3.36	2.16	3.36
Tree-Vegetation Cover	5.0	35.0	5.0	35.0	5.0	35.0	5.0	35.0
Leaf Area Index	2.00	2.26	1.31	2.00	1.17	1.67	1.67	2.26
Fraction of current leaves	0.67	0.67	0	0.33	0	0.61	0.66	0.77
Fraction of previous leaves	0.23	0.23	0.57	0.9	0.39	1	0.13	0.23
Fraction of shoots	0.1	0.1	0.1	0.1	0.1	0.1	0.1	0.1
Leaf Area Distribution	Spherical		Spherical		Spherical		Spherical
Optical Measurements
Leaf reflectance (ρleaf)	Li-Cor 1800-12 + ASD Fieldspec® 3 measurements
Leaf transmittance (τleaf)	Li-Cor 1800-12 + ASD Fieldspec® 3 measurements
Bark reflectance (ρbark)	Plant Probe + ASD Fieldspec® 3 measurements
Background reflectance(ρs)	PROSAIL LUT

shows the ranges of FLIGHT input parameters used in each phenophase.

As in PROSAIL, some of the variables used in the FLIGHT model were limited according to field data measured in each phenophase. The tree parameters with the larger seasonal variability were LAI_tree, F_current, and F_previous. For the sake of simplicity, other variables such as leaf size or F_Cover were considered seasonally independent. Current and previous ρ_leaf and τ_leaf, measured with the integration sphere (Section 2.2.2), were selected randomly as model inputs for each phenophase. Even if they are not model inputs, the corresponding leaf biophysical properties of these leaves were known and, therefore, could be related to RF in the LUT. Tree canopy structural parameters estimated from ALS [30] or field measurements were used as fixed inputs, see

Table 4. Mean tree canopy structure values obtained for the study area.

Crown Shape (m)	Crown Radius (m)	Centre to Top Distance (m)	Height to First Branch (m)		Trunk DBH (m)
			Min.	Max.
Ellipsoid	4.20	2.70	1.50	3.50	0.45

. As, at this stage, we did not aim to realistically represent the ecosystem at a very high spatial resolution, but rather simulate reflectance factors on an ecosystem scale, averaged values and arbitrary locations of tree crowns were used. However, the FLIGHT model can describe the position as well as specific structural parameters of each tree individually, if required.

Figure 2 describes the methodology followed to couple the PROSAIL and FLIGHT models. As shown in Table 1, both models share several variables related to the illumination and observation conditions.

The diffuse fraction of the down-welling irradiance was estimated following the approach implemented in FLIGHT (Equation (1)):

$$diff=ryfrac·\left(1.0-\exp \left(-1.0·T_{\mathrm{Rayleigh}}\right)\right)+aerfrac·\left(1.0-\exp \left(-1.0·T_{\mathrm{aero}}\right)\right),$$

(1)

where ryfrac (0.5) and aerfrac (0.75) are the fraction of Rayleigh- and aerosol-scattered light to reach the ground surface. $T$_Rayleigh is the Rayleight Optical Depth and is given by the following expression [64]:

$$T_{\mathrm{Rayleigh}}=\frac{\frac{P}{P_0}·\left(0.008569·W{b}^{-4}·\left(1.0+0.0113\ast W{b}^{-2}+0.00013·W{b}^{-4}\right)\right)}{\cos \left({\theta }_{\mathrm{s}}\right)},$$

(2)

where P is the atmospheric pressure (1013.25), P₀ is the atmospheric pressure at an altitude of 0 m (1013.25), Wb is the Wavelength, and θ_s is the Solar Zenith Angle. Then, Τ_aero is the Aerosol Optical Depth calculated as:

$$T_{\mathrm{aero}}=\frac{dmeff·a·W{b}^k}{\cos \left({\theta }_{\mathrm{s}}\right)},$$

(3)

where dmeff is the effective diffuse fraction (0.55) used in FLIGHT and k (−1.25) is an approximate Angstrom coefficient for a continental aerosol [65], and a = Aerosol Optical Depth (AOT)/(0.55^k). Once the diffuse radiation was computed, the fractions of diffuse (diff_frac) and direct (dir_frac) radiation were calculated from the following expressions:

$$dif{f}_{\mathrm{frac}}=\frac{diff}{dir+diff},$$

(4)

$$di{r}_{\mathrm{frac}}=1-dif{f}_{\mathrm{frac}},$$

(5)

where dir is the direct radiation calculated as:

$$dir=\exp \left(-1.0·T_{\mathrm{tot}}\right),$$

(6)

where Τ_tot is the sum of Rayleigh Optical Depth and AOT. Finally, the directional reflectance factor was calculated through the weighting of the two fractions:

$$RF=di{r}_{\mathrm{frac}}·BRF+dif{f}_{\mathrm{frac}}·HDRF,$$

(7)

Predicted RF is a linear combination of the black and white sky simulations of FLIGHT using the modeled grass BRF and HDRF, respectively, as a background. Both runs of the model are weighted according to the fractions of diffuse and direct down-welling radiance in each band. The FLIGHT model includes a bidirectional reflectance distribution function (BRDF) simulation based on Hapke equations [66] that is used to simulate background directional reflectance. The computation of the BRDF grid (5° resolution) for each of the spectral bands simulated would have been too expensive. For this reason, we simplified the problem by using the HDRF as a proxy of the integrated background RF for diffuse illumination. As a consequence, the contribution of BRF is slightly overestimated since this factor is also used as a response for photons scattered from and through the tree canopy; therefore, changing their orientation. However, this effect is expected to be low in remote sensing applications where high solar elevation and observation angles are predominant. In this case, most of the observed background is directly illuminated from the sun and reflected photons show lower orders of magnitude. On the other hand, this approximation achieves a better description of the grass background BRDF than the original soil model.

2.4. Quantitative Comparison of Simulated and Measured Reflectance Factors

In this section, the performance of the PROSAIL and PROSAIL+FLIGHT models have been evaluated against simultaneous or quasi-simultaneous observations of RF and biophysical parameters on proximal (ground) and remote sensing (airborne and satellite) scales, respectively.

2.4.1. PROSAIL Simulated Versus Ground-Measured Reflectance Factors of the Grass Layer

In order to investigate the ability of PROSAIL to mimic RF of the grass cover, we compared PROSAIL’s predictions with field observations. Model inputs were sampled in three campaigns carried out in 2017 (16-Mar-2017; 21-April-2017; 19-May-2017) over five 25 cm × 25 cm quadrants. Grass destructive sampling provided the model parameters C_ab, C_ar, C_m, C_w, and grass LAI (LAI_grass). In PROSAIL, N parameter was set to 1 assuming predominant monocotyledons species in the grass layer. Since no observations of LAD were available, a RF was simulated using each of the six different predefined LAD types [67]: (a) Planophile (b) Erectophile (c) Plagiophile (d) Extremophile (e) Spherical, and (f) Uniform. Then, we selected the LAD achieving the best fit with each observation.

AOT values were randomly set between the minimum (0.06) and maximum (0.16) according to values provided by Sentinel-2 L2A-AOT products in clear days during 2017. Regarding C_brown, a preliminary test (not shown) of the variable was performed by running the model with minimum and maximum values set to 0 and 3, respectively, and steps of 0.25. Due to the strong uncertainties both in the definition and the realistic boundaries of this parameter in our ecosystem we decided to use a conservative fixed value of 0.1.

The similarity between simulated and measured spectra was quantified using the Relative Root Mean Square Error (RRMSE) [59] and the Spectral Angle Mapper (SAM) technique [68]. SAM calculates the spectral angle (SA) between the observed and the predicted spectra by treating them as vectors in an n-dimensional space where n is the number of spectral bands [68], see Equation (8).

$$SA\left(V,W\right)={\cos }^{-1}\left(\frac{{{\displaystyle \sum }}_{i=1}^n{V}_i{W}_i}{{\left({{\displaystyle \sum }}_{i=1}^n{V}_i{}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\ast {\left({{\displaystyle \sum }}_{i=1}^n{W}_i{}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\right),$$

(8)

where V is the observed spectrum and W is the modeled spectrum.

The SAM technique was selected because it offers complementary information to the RRMSE. It enables the quantification of the differences in the shape of the compared spectra. This can be especially relevant if they affect regions/bands used for the calculation of vegetation indices as the changes in the shape can alter the relative difference between those bands and, therefore, the index value. This is particularly important when the simulated dataset is used to explore the performance of spectral indices to estimate biochemical or structural vegetation parameters [69,70,71].

2.4.2. PROSAIL+FLIGHT Simulated Versus Airborne and Satellite Measured Reflectance Factors at the Ecosystem Level

The coupled PROSAIL+FLIGHT model simulated the RF at the ecosystem level of pre-defined scenes for two remote sensors: the airborne CASI, see

Table A1. Input variables for the generation of simulated scenes of the CASI sensor. LAI was specified for each of the trees in the plot (units in Table 1).

Variable	08/April/2014		23/May/2015		03/July/2015		03/May/2016
	PlotA	PlotB	PlotA	PlotB	PlotA	PlotB	PlotA	PlotB
PROSPECT-5B
Leaf structural parameter	1	1	1	1	1	1	1	1
Equivalent Water Thickness	0.024	0.019	0.030	0.017	0.009	0.013	0.032	0.023
Dry matter content	0.007	0.007	0.009	0.007	0.014	0.013	0.008	0.006
Chlorophyll a+b content	37.48	45.86	35.24	20.89	0.1	0.1	45.41	18.67
Carotenoids content	10.13	8.22	10.54	12.39	0.1	0.1	6.22	6.80
Brown pigments	0.1	0.1	0.1	0.1	0.1	0.1	0.1	0.1
SAIL
Leaf Area Index	2.24	1.65	1.54	1.65	0.39	0.64	3.02	2.16
Hot Spot	0.001	0.001	0.001	0.001	0.001	0.001	0.001	0.001
Grass Leaf Area Distribution	---	---	---	---	---	---	---	---
FLIGHT
Leaf Size	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0
Tree-Vegetation Cover	22.1	12.9	22.1	12.9	22.1	12.9	22.1	12.9
Leaf Area Index	---	---	---	---	---	---	---	---
Fraction of green leaves	0.3	0.3	0.3	0.3	0.15	0.15	0.3	0.3
Fraction of senescent/shoot	0.7	0.7	0.7	0.7	0.85	0.85	0.7	0.7
Tree Leaf Area Distribution	Spherical		Spherical		Spherical		Spherical
Shared Variables
Solar Zenith Angle	33.1	33.1	27.8	27.8	17.9	17.9	24.4	24.4
Solar Azimuth Angle	168.6	168.6	168.7	168.7	159.2	159.2	168.4	168.4
View Zenith Angle	0	0	0	0	0	0	0	0
View Azimuth Angle	105	105	105	105	105	105	105	105
Aerosol Optical Depth	0.06	0.06	0.11	0.11	0.16	0.16	0.16	0.16

, and the MSI onboard Sentinel-2 see

Table A2. Input parameters for the generation of simulated scenes of Sentinel-2 (units in Table 1).

Variable	21/February/2017		22/April/2017		22/May/2017
	PlotA	PlotB	PlotA	PlotB	PlotA	PlotB
PROSPECT-5B
Leaf estructural parameter	1	1	1	1	1	1
Equivalent Water Thickness	0.033	0.020	0.020	0.017	0.034	0.023
Dry matter content	0.008	0.007	0.009	0.008	0.009	0.009
Chlorophyll a+b content	38.79	31.58	25.34	25.06	33.01	54.19
Carotenoids content	9.67	8.26	6.28	6.36	8.62	13.86
Brown pigments	0.1	0.1	0.1	0.1	0.1	0.1
SAIL
Leaf Area Index	1.46	0.88	0.75	0.97	0.98	1.39
Hot Spot	0.001	0.001	0.001	0.001	0.001	0.001
Grass Leaf Area Distribution	---	---	---	---	---	---
FLIGHT
Leaf Size	3.0	3.0	3.0	3.0	3.0	3.0
Tree-Vegetation Cover	22.1	12.9	22.1	12.9	22.1	12.9
Leaf Area Index	---	---	---	---	---	---
Fraction of current leaves	0.3	0.3	0.3	0.3	0.15	0.15
Fraction of previous leaves	0.7	0.7	0.7	0.7	0.85	0.85
Tree Leaf Area Distribution	Spherical		Spherical		Spherical
Shared Variables
Solar Zenith Angle	42.1	42.1	28.3	28.3	20.7	20.7
Solar Azimuth Angle	168.2	168.2	168.7	168.7	166.9	166.9
View Zenith Angle	0	0	0	0	0	0
View Azimuth Angle	105	105	105	105	105	105
Aerosol Optical Depth	0.12	0.12	0.41	0.41	0.29	0.29

. The results obtained with the model were resampled to the bands of both sensors using the corresponding spectral response functions. At the same time, simulations with the original PROSAIL and FLIGHT models were performed separately for comparison. The simulated scenes corresponded to two field plots of 75 m × 75 m with a tree F_Cover of 22.1% (PlotA) and 12.9% (PlotB), see Figure 3.

For these simulations, FLIGHT and PROSAIL+FLIGHT accounted for the actual positions of the trees and the individual LAI_tree measured in the field, and PROSAIL used the value of the biophysical variables measured in the field to predict grass RF. In addition, ρ_leaf and τ_leaf of current and previous leaves measured with an integrating sphere Li-Cor-1800-12 (see Section 2.2.2) in different phenological stages of the trees were used as model inputs. For the bark spectrum, an average reflectance value from field measurements described in Section 2.2.2 was used, assuming that this tree element does not undergo variations throughout the year. Table A1 and Table A2 show the input variables used for each plot and date for the sensors CASI and MSI, respectively. For the CASI sensor, the AOT data used were obtained from the CE312 IR radiometer CLIMAT (https://www.cimel.fr/?instrument=radiometer-ir-climat-benchmark&lang=en) operated by INTA during the airborne campaigns. For the Sentinel-2 images, the L2A-AOT product (https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/processing-levels/level-2) provided by this sensor was used.

Field measurements of vegetation biophysical parameters feed the different models. Observations were simultaneous to the airborne campaigns, with a maximum difference of one day. In the case of two of the Sentinel-2 images selected for this comparison, field campaigns occurred less than three days from the dates of the overpasses (22-April-2017 and 22-May-2017). However, in the case of the 21-February-2017 overpass, fieldworks could only take place 22 days later. Despite the mismatch, we decided to keep this date for the analysis as we could confirm that, in 2017, grasses remained in a stable situation in the study site between February and mid-March (data not shown).

2.5. Angular Evaluation of PROSAIL+FLIGHT

In this section, we compared the RF generated at different observation angles PROSAIL, FLIGHT, and PROSAIL+FLIGHT. A comparison was carried out for two spectral bands—Red (670 nm) and NIR (800 nm)—in the principal and the cross-principal planes. Input parameters for this simulation are summarized in

Table A3. Input parameters for the comparison between the PROSAIL and FLIGHT models and the PROSAIL+FLIGHT coupling.

Variable	Grass	Tree	Eco
Leaf estructural parameter	1.0	2.1	1.2
Equivalent Water Thickness	0.024	0.014	0.023
Dry matter content	0.007	0.015	0.008
Chlorophyll a+b content	55.70	70.40	57.72
Carotenoids content	10.56	18.76	11.69
Brown pigments	0	0	0
Leaf Area Index	2.24	1.79	2.60
Hot Spot	0.001	0.001	0.001
Leaf Area Distribution	Spherical	Spherical	Spherical
Leaf Size	---	3.0	3.0
Tree Fcover	20.0	20.0	20.0
Fraction of current year leaves	---	0.25	0.25
Fraction of previous year leaves	---	0.75	0.75
Solar Zenith Angle	20.0	20.0	20.0
Solar Azimuth Angle	180.0	180.0	180.0
View Zenith Angle	−75 to 75	−75 to 75	−75 to 75
View Azimuth Angle	Principal an cross plane	Principal an cross plane	Principal an cross plane
Aerosol Optical Depth	0.06	0.06	0.06

. For PROSAIL and PROSAIL+FLIGHT, BRF, HDRF, and RF were generated since the models provide BRF, HDRF as outputs that are then combined into the RF as a function of the spectral fractions of dir_frac and diff_frac respectively. FLIGHT directly provided the RF since the code internally accounts for the interaction of direct and diffuse irradiance with the pixel elements, and the contributions of the direct and the diffuse flux are not provided as an output.

While PROSAIL+FLIGHT grass parameters (P_grass) and tree parameters (P_tree) are provided separately, for PROSAIL and FLIGHT we scaled them to ecosystem level (P_eco). To do so, we assumed Spherical LAD for grasses and trees. LAI was the only structural parameter scaled to ecosystem level (LAI_eco), fot that we assumed that the grasses covered the whole surface and the trees only a fraction of it, the F_cover, see Equation (9):

$$LA{I}_{\mathrm{eco}}=LA{I}_{\mathrm{grass}}+LA{I}_{\mathrm{tree}}·\left(F_{\mathrm{cover}}/100\right),$$

(9)

where LAI_tree is the LAI of tree crowns. In the case of leaf parameters, these were scaled according to the LAI of each of the vegetation types and trees F_cover, see Equation (10):

$$P_{\mathrm{eco}}=\frac{P_{\mathrm{grass}}·LA{I}_{\mathrm{grass}}+P_{\mathrm{tree}}·LA{I}_{\mathrm{tree}}·\left(F_{\mathrm{cover}}/100\right)}{LA{I}_{\mathrm{eco}}},$$

(10)

3. Results

3.1. PROSAIL Simulation of the Grass Layer

Time series of MODIS-NDVI showed the strong seasonal variability of vegetation dynamics, see Figure 4, which is mainly driven by the phenology of the grass cover. NDVI presented values from 0.55 to 0.80 between October/November to May. Values were lower than 0.40 during the summer period (July to September) when the grass layer is completely dry. Inter-annual variability was also strong; exceptionally dry years in 2005 or 2009 contrast with very humid ones, as in 2012 and 2015. In the dry years, the growth season and green peak are significantly reduced with peak NDVI below 0.60.

This strong inter-annual variability justified the definition of transition dates with overlap between two phenophases, so that the accuracy of the RTM is not conditioned to LUT selection during transition periods. Therefore, we defined (a) constant transition dates with overlap, taking into account the maximum amplitude of each phenophase observed in the MODIS-NDVI time series, see

Table 5. Phenophases and constant transition dates established in the study area for the period 2000–2016.

Phenophase	Start	End
Phenophase 1: Summer drought	02-July	11-November
Phenophase 2: Autumn regrowth	06-October	23-April
Phenophase 3: Biomass peak	23-December	14-June
Phenophase 4: Beginning of the grass decay	12-April	29-September

; and (b) annual boundaries that were used to classify each of the individual field campaigns (2013–2017) in one of these phenophases for model parameterization and validation. Therefore, the input values in PROSAIL simulation were defined accordingly using field data collected in each period, see

Table 6. Start date of each phenophase between 2011 and 2016. The sampling campaigns carried out in this period were assigned to each of the phenophases following this classification.

Year	Phenophase 1: Summer Drought	Phenophase 2: Autumn Regrowth	Phenophase 3: Biomass Peak	Phenophase 4: Beginning of the Grass Decay
2011–2012	24-July-2012	05-November-2011	02-February-2012	30-April-2012
2012–2013	07-August-2013	16-October-2012	05-February-2013	25-May-2013
2013–2014	25-July-2014	20-October-2013	31-January-2014	13-May-2014
2014–2015	03-July-2015	06-October-2014	23-January-2015	13-May-2015
2015–2016	31-December-2016 (End of the MODIS time-series)	24-October-2015	12-February-2016	02-June-2016

.

The longest phenophases are number 2 and 3, corresponding with the beginning and end of the growing season for the grass layer. The interannual variability of these two phenophases is higher because they depend, to a large extent, on the timing and amount of rainfall collected during the year and, especially, at the beginning of the spring [72,73].

Results of the PROSAIL simulations for each phenophase are shown in Figure 5 (grey lines). The variability derives from the configuration of the LUT, see parameters and bounds in Table 2, and the subsequent LHS. For qualitative assessment, we have overlaid spectra measured in the study site in different field campaigns using the ASD FieldSpec^® 3 spectroradiometer.

In most of the cases, field observations fell within the distributions simulated from the LUT. Observations presented a lower variability as they correspond to specific campaigns and thus, they do not cover the full range of grass conditions considered in each phenophase. However, there are some significant discrepancies between the measured and simulated spectra in phenophase 2, corresponding with autumn regrowth, see Figure 5b. In this case, NIR RF are systematically overestimated by the model. A similar effect is observed in phenophase 1, corresponding with the summer drought period, see Figure 5a, where RF are slightly overestimated between 900 and 1300 nm. The spectral behavior of the grass layer in the autumn regrowth and in the biomass peak is quite similar as fresh green plants are present in both periods. Although, in the case of autumn regrowth, we can see the influence of standing-dead material remnant from the senescent period, see Figure 5b. This effect is further analyzed in the discussion section.

3.2. PROSAIL+FLIGHT Simulation of the Ecosystem

The combination of the grass BRF and HDRF LUT generated by coupling PROSAIL and FLIGHT models provides the ecosystem RF for the different phenophases, see Figure 6. Maximum (red line) and minimum (blue line) boundaries of the simulated directional reflectance of the grass layer using the PROSAIL model have been overlaid to PROSAIL+FLIGHT simulations (grey lines) so we can have a first qualitative insight into the differences between the two models. Seasonal differences are still remarkable between the summer drought and the rest of the phenophases. In Mediterranean TGEs, the tree cover is less dynamic than the grass layer and tree F_Cover is low; therefore, for the simulated conditions (at nadir observations and relatively low θ_s), the effect of the trees is reduced. In general, a small reduction of RF is observed in PROSAIL+FLIGHT compared with PROSAIL simulations, especially in the short-wave infrared (SWIR) region. These changes are more noticeable in phenophase 1 (summer drought), where there is a larger difference between the RF of the grass and the tree canopies, as the grass layer is completely dry but the trees are still green.

Figure 7 compares the PROSAIL+FLIGHT RF simulated in each phenophase with the CASI imagery. CASI RF have been extracted from the eleven 25 m × 25 m sampling plots in the four images used in this study). For phenophase 1 (summer drought) a CASI image acquired on 03-July-2015 was used to compare with simulated spectra, see Figure 7a; for phenophase 3 (biomass peak) the images of 08-April-2014 and 03-May-2016 were selected, see Figure 7b; and finally, the image of 23-May-2015 was used to compare with phenophase 4 (beginning of the grass decay) simulation, see Figure 7c. For phenophase 2 (autumn regrowth) no CASI image was available. In general, a good correspondence between simulated and measured spectra is observed.

As with the grass LUT generated with PROSAIL, shown in Figure 5, the range of PROSAIL+FLIGHT simulated spectra is greater than the one collected in the field during specific campaigns since the LUT represents the variability of the whole phenophase. A slight tendency to overestimate RF in the NIR region of the simulated spectra is observed in phenophase 1 (summer drought).

3.3. Quantitative Comparison of Simulated and Measured Reflectance Factors

3.3.1. PROSAIL Simulated versus Ground Measured Reflectance Factors of the Grass Layer

The RF simulated with PROSAIL were compared with those acquired using the ASD Fieldspec^® 3 on 1 m × 1 m grass sub-plots measured in three different campaigns (16-March-2017; 21-April-2017; 19-May-2017) and differences were quantified using RRMSE and SA metrics.

Figure 8 shows data from 16-March-2017. PROSAIL simulations show good agreement with measured RF with a RRMSE between 9.64% and 19.04%, except for plot P5, see Figure 8e, where RRMSE is above 30 %. The angular distances between simulated and measured spectra never exceed 5.0°. In the graphs, larger differences seem to be located in the NIR region, especially in the NIR plateau, and in the SWIR region, see Figure 8b,d,e. However, if we analyze the RRMSE obtained in each of the three main spectral regions observed: VIS, NIR, and SWIR, see

Table A4. RRMSE (%) between measured (ASD FieldSpec^® 3) and simulated (PROSAIL) spectra in the VIS (400–700 nm), NIR (701–1200 nm) and SWIR (1201–2500 nm) spectral regions and for the whole spectrum (400–2500 nm). In bold regions with highest RRMSE per plot.

Plot	16-March-2017				21-April-2017				19-May-2017
	VIS	NIR	SWIR	Total	VIS	NIR	SWIR	Total	VIS	NIR	SWIR	Total
P1	5.17	7.70	10.79	9.64	30.59	19.48	8.49	9.88	5.19	4.27	12.35	8.10
P2	22.29	19.47	17.30	19.04	27.51	4.15	5.89	6.42	21.29	17.01	19.41	17.28
P3	18.84	12.81	6.74	10.84	8.12	9.87	10.01	10.13	10.86	4.91	8.16	11.37
P4	32.69	9.83	25.91	18.79	19.33	7.20	11.35	11.39	18.72	12.52	18.33	14.64
P5	24.19	17.37	40.37	31.32	14.71	18.61	8.56	13.72	21.35	17.13	12.38	14.77

, we find the largest RRMSE in the VIS region in plots 2 to 4 and in the SWIR in plots 1 and 5. This field campaign occurred in the overlap period between phenophases 2 (autumn regrowth) and 3 (biomass peak), when the grass cover in some of the sampling plots still presents a high mixture between standing-dead material and fresh green grass. This mixture produces a systematical overestimation of PROSAIL simulations in the NIR plateau, as can be observed in Figure 8b,c. On this date, Planophile LAD most often achieved the best match between the simulated and measured spectra.

Figure 9 shows the data corresponding to 21-April-2017. On this date, a strong mixture of green and senescent fractions was observed. Senescent material was predominant in some of the sampled sub-plots. This mixture is typical of phenophase 4 (beginning of the grass decay), that occurred quite early during 2017 due to limited rainfall during the first months of the spring. Precipitation slightly increased at the end of April, producing a re-greening. In this period, the agreement between simulated and measured spectra is better than on 16-March-2017, see Figure 8, with a RRMSE below 14% in all the plots, and a SA between 2.7° and 5.2°. The relative differences for spectral regions show larger differences in the VIS for most of the plots, see Table A4. On this date, the selected LAD that best matched the simulated and measured spectra was diverse: Planophile, Spherical, Plagiophile, and Uniform.

Figure 10 shows data from a re-greening episode on 19-May-2017 after some late spring rain events. On this date, the grass status varied between phenophase 2 (biomass peak) and 3 (beginning of the grass decay). In some plots, Figure 10a,c, a good correspondence between the measured and simulated spectra is observed. RRMSE are below 12% and SA are close to 3°, and the largest absolute differences are found in the SWIR region. In other cases, Figure 10b,d, PROSAIL overestimates field observations. RRMSE reaches up to 17.28%, SA is still ~ 3.0°. In P5, Figure 10e, a mismatch between simulated and measured spectra is observed in the VIS and NIR regions with RRMSE are 21.35% and 17.31%, respectively, see Table A4. This is quite similar to what occurs in the P2, see Figure 10b, but with a higher SA of 7°. The Erectophile LAD is the most frequently selected on this date; it is typical of grass species which are predominant at this phenophase as they endure longer than the forbs (Spherical) and legumes (Planophiles).

3.3.2. PROSAIL+FLIGHT Simulation

Figure 11 compares simulated PROSAIL+FLIGHT versus measured CASI RF for two 75 m × 75 m plots located in the study site, see Figure 3. Due to the low tree F_Cover (<22%), the ecosystem reflectance factor was quite similar to the background spectrum simulated with PROSAIL for a quasi-nadir observation, see the black dashed lines in Figure 11; however, both spectra are more dissimilar in the summer period, see Figure 11d, when the grass cover is completely dry but the tree cover still provides a “green vegetation” signal. Differences between simulated and measured spectra were lower for all dates in PlotB that also has lower tree F_Cover (~13%) than in PlotA (F_Cover ~ 22%). The lowest RRMSE (8.09%) was found in PlotB for the April campaign, while the highest (31.07%) was found in PlotA for the July campaign. In spite of this, the shape of the simulated and measured spectra are quite similar for this summer date, as confirmed by a low SA (3.3°) calculated using SAM. The largest SA (5.1°) was recorded in PlotB during the 03-May 2016 field campaign. In this case, the PROSAIL+ FLIGHT simulation increases the differences with the CASI measurement when compared with the PROSAIL simulation of the grass layer.

Comparison with CASI images for the two main spectral regions, see

Table A5. RRMSE (%) between measured (CASI) and simulated (PROSAIL+FLIGHT) spectra in the VIS (400–700 nm) and NIR (701–1200 nm) spectral regions and for the whole spectrum (400–1200 nm).

Plot	08-April-2014		03-May-2016		23-May-2015		03-July-2015
	VIS	NIR	VIS	NIR	VIS	NIR	VIS	NIR
PlotA	25.83	15.89	43.03	24.87	35.44	20.57	46.35	25.43
PlotB	25.57	5.55	50.69	10.87	32.40	9.00	35.49	13.26
Total RRMSE
PlotA	19.81		30.75		25.37		31.07
PlotB	8.09		16.33		12.75		18.86

, shows that, as observed in the PROSAIL simulation, see Table A4, differences between simulated and measured spectra at the ecosystem level are significantly lower in the NIR than in the VIS region; where RRMSE reached 50% in the worst of the cases (plotB in May 2016). On 08-April-2014, RRMSE were the lowest, especially in the NIR (~5.00%). On the contrary, the date disclosing the highest differences was 03-July-2015, especially in plotA where the RRMSE in the visible and in the NIR were 46.35% and 25.43%, respectively. These are similar to the RRMSE obtained for the same plot in 03-May-2016.

Figure 12 compares PROSAIL+FLIGHT simulations with MSI Sentinel-2 RF for the same 75 m × 75 m plots (PlotA and PlotB) used in the previous comparison with CASI images and for three dates: 21-February-2017, 22-April-2017, and 22-May-2017. PlotB in the Sentinel-2 image on 22-May-2017 was discarded due to the presence of clouds. Figure 12 also compares the RF generated by PROSAIL analogously to Figure 11. As can be observed, simulated and measured spectra were very similar on 22-May-2017 for PlotA, see Figure 11c, with RRMSE close to 7% and SA ~ 3.0°. A similar pattern was found for this plot in April but only in the NIR region, while MSI bands 2 and 4 (blue and green spectral regions) showed larger errors. PlotB presents a strong mismatch on both dates with RRMSE greater than 25% in all cases, unless SA values are much lower than in PlotA, as for 21-February-2017.

Table A6. RRMSE (%) between measured and simulated RF for MSI Sentinel-2 bands.

Band	21-February-2017		22-April-2017		22-May-2017
	PlotA	PlotB	PlotA	PlotB	PlotA	PlotB
B2 (490 nm)	42.00	9.78	52.70	57.83	25.43	---
B3 (560 nm)	22.22	19.63	24.95	39.25	11.10	---
B4 (665 nm)	42.72	5.86	41.14	45.62	18.95	---
B5 (705 nm)	20.10	9.63	20.83	36.76	12.77	---
B6 (740 nm)	25.19	26.41	3.27	24.80	1.89	---
B7 (783 nm)	30.29	29.55	2.03	23.15	0.30	---
B8a (865 nm)	23.44	26.50	6.16	23.55	1.41	---
B11 (1610 nm)	22.16	30.50	0.44	11.06	6.57	---
B12 (2190 nm)	35.60	39.82	5.19	2.61	1.44	---
Total	30.22	31.77	13.96	26.64	7.07	---

shows the RRMSE values between simulated and measured RF for all MSI Sentinel-2 bands simulated with PROSAIL+FLIGHT. As for CASI images, larger differences were found in the visible bands but especially for bands 2 (Blue) and 4 (Red) with maximum RRMSE over 50% on 22-April-2017. The best agreement between simulated and measured RF was found in the Red Edge NIR and SWIR bands (bands 5 to 12), with minimum RRMSE below 0.5% in April and May, both in PlotA.

3.4. Angular Evaluation of PROSAIL+FLIGHT

Figure 13a–d compares the BRF predicted by PROSAIL+FLIGHT and PROSAIL in the two main planes and two spectral bands (650 and 800 nm) at different view zenith angles (θ_v). In the principal plane, PROSAIL predicts brighter BRF than PROSAIL+FLIGHT, however, near the hotspot in the NIR, PROSAIL+FLIGHT show larger values, see Figure 13a,c. PROSAIL BRF are also brighter in the cross-principal plane, see Figure 13b,d. NIR BRF show different behaviors at large θ_v for each of the models, increasing with θ_v for PROSAIL and decreasing in the case of PROSAIL+FLIGHT due to the stronger influence of the volumetric or the geometric scattering in each case. In the Red region, both models behave similarly due to the stronger absorption of light. PROSAIL+FLIGHT HDRF show similar angular responses to PROSAIL HDRF, but always with lower values, see Figure 13e–h. PROSAIL+FLIGHT HDRF presents a slight instability, especially in the Red region due to Monte Carlo noise of the FLIGHT model [34]. Figure 13i–l presents the RF resulting from the combination of direct and diffuse illumination in the models. In this case, we also present the RF predicted by FLIGHT, using a Lambertian background reflectance which equals the RF (combining BRF and HDRF) of the grass layer predicted by PROSAIL. In the Red region, FLIGHT presents the brighter RF except near the hotspot region in the principal plane, where the other models present more acute hotspots, see Figure 13i,j. In the NIR region, PROSAIL and FLIGHT simulate the brighter and the darker RF, respectively, see Figure 13k,l. PROSAIL+FLIGHT present intermediate values and as FLIGHT, RF decreases with θ_v, denoting a geometric scatter. PROSAIL predicts RF that decrease and increase with θ_v in the Red and the NIR regions, respectively. In all the cases, FLIGHT presents a reduced hotspot.

4. Discussion

The combination of the 3-D RTM FLIGHT and the 1-D RTM PROSAIL allowed the simulation of RF representative of Mediterranean TGE as a function of the parameters of the two vegetation layers present in the ecosystem. We tested the performance of the PROSAIL and PROSAIL+FLIGHT models at grass and ecosystem scales respectively. Lower RRMSE found at close range can be, in part, explained by spatial scaling uncertainties. The comparison of both models on an ecosystem scale and at different angular configurations suggest that PROSAIL could mimic RF of the ecosystem when optical and/or biophysical properties of the grass and the tree layers are similar and observation and illumination angles are high. However, the predictive capability of the 1-D model would decrease when differences between layer properties or sun/view angles increase and the contribution of tree crowns and the 3-D structure of the ecosystem play a stronger roll on the signal.

During the model evaluation, larger uncertainties were found at the ecosystem level than at the grass plot scale. These differences are, in part, explained by the larger spatial mismatch existing between the remote sensing spectral footprint and the small quadrants where biophysical-sampling-provided model inputs were taken within the 75 m × 75 m plots. The limited biophysical sampling might not be able to represent the properties of a large pixel in grassland characterized by a high spatial variability. However, in the simulations performed with Sentinel-2, see Figure 12, the coupled PROSAIL+FLIGHT represented pixel RF better than the PROSAIL simulation of the background, see Figure 12b,c. Other sources of uncertainty that explain the differences between simulations and observations could be a) the uncertainties in the laboratory estimation of the biophysical variables, b) the selection of AOT and other parameters not measured such as LAD, C_brown, or the structural parameter N, and c) inaccuracy of the models to represent vegetation radiative transfer of specific ecosystem components. In fact, our results suggest that current RTM do not properly represent the optical properties of standing death material, which induces large differences between predictions and observations when this is present. This jeopardizes the use of physical models to simulate RF and estimate vegetation parameters in Mediterranean TGE and, in general, in grasslands which accumulate large amounts of grass litter over time. Also, PROSAIL+FLIGHT does not account for the change of angular distribution of photons falling onto the grass due to scattering and diffuse transmission within the tree crowns when BRF is represented. Uncertainties induced by this model error would be larger in the NIR than in the VIS region, where crown transmittance and scattering are larger. At high observation and illumination angles, tree crowns will tend to occlude the same areas that are occluded to direct illumination, minimizing, in general, the observation of areas where the background BRF is more overweighed by the model with respect to HDRF.

PROSAIL+FLIGHT simplifies the representation of the background BRDF without mapping a BRDF for each spectral band but rather separating the model into BRF and HDRF components. This approach speeds up computation but is still more accurate than assuming a Lambertian vegetation background. Errors induced by this simplification are expected to be comparable to other errors induced by model assumptions such as Lambertian symmetric leaves, omission of clumping, or the homogeneous distribution of leaves within the canopy, as well as the problems of the RTM regarding the accurate representation of vegetation optical properties ([45], this article). Also, the simulation of the illumination conditions adopted by FLIGHT may become an oversimplification. The combination of vegetation and atmospheric radiative transfer models could improve the atmospheric correction and retrieval of vegetation parameters [74]. Inaccuracies in the simulation of illumination conditions may account for part of the errors found in the different comparisons performed here with remote sensing data. The use of complex atmospheric RTM such as MODTRAN [75,76], or 6S [77] could help to improve the quality of the simulations if their input parameters can be accurately quantified.

Field, airborne, and satellite data assessed the capability of the model to reproduce RF measured at different scales. Due to the lack of observations of some of the input variables, see Table 1, the comparison between simulated and measured spectra cannot be strictly considered a validation. However, it allowed us to understand which parts of the model might produce most of the uncertainties. PROSAIL and PROSAIL+FLIGHT simulations were separately compared with ground and airborne/satellite observations. This exercise revealed that most of the differences between simulated and observed RF at the ecosystem scale were induced by the background RTM (PROSAIL), which failed to adequately represent the optical properties of death stand material. Another factor explains part of the uncertainties found in the comparison of models with remote sensing data; the high spatial variability of the grass layer hampers the use of a few samples to represent the average biophysical state of a remote sensing plot. This issue is inherent to grasslands with a high biodiversity and heterogeneous distribution of species in space and time, such as the ecosystem under study. Migliavacca et al., [40] found that differences in LAD related to nutrient availability could explain strong differences in the optical properties of the grassland layer in the same study site of this article. Former works accounted for 145 different species throughout the year in the study area [44]. Such biodiversity makes unfeasible scaling through species-based metrics such as in [78], and makes valuable use of modeling approaches that can be validated at a proximal sensing scale to later produce remote estimates of vegetation properties, as in this work.

The comparison of field observations with the simulated RF at a proximal sensing scale minimized spatial mismatch between spectral and biophysical footprints. This allowed us to detect the systematic overestimation of NIR RF in phenophase 2: Autumn regrowth, as well as in phenophase 1: Summer drought, see Figure 5a,b. We hypothesize that such an overestimation could be related to the presence of, not only senescent, but decomposing grass material, whose optical properties are not well represented by the PROSPECT model, see Figure 14. This standing death material would feature a low reflectance in the NIR caused by a strong deterioration of the cell walls [79] and the degradation of spongy mesophyll. This would translate into a change in the leaf refractive index in PROSPECT that has not been characterized. Such a model error would increase uncertainties in the estimation of biophysical variables both from RTM (this work) and empirical approaches [80]. Future studies could focus on characterizing the optical properties of this decomposing litter, which would probably improve the simulation of grass RF in phenophase 2 (autumn regrowth) at the beginning of the growth season. Despite these imbalances, it is noteworthy that the angular differences measured with SAM between the simulated and measured spectra of the grass layer do not generally exceed 5°; therefore, the shape of the spectra is similar but only differences in albedo are observed.

Uncertainties, in the SWIR and NIR regions, might also be related to errors in the representation of soil RF. In this work, soil RF was predicted as a function of SWC measured at a depth of 5 cm at a few points located in the center of the study area. SWC can greatly vary with meso- and micro-topography; therefore, a spatial heterogeneity in the ecosystem might not have been accurately simulated [81,82]. The soil influence on the observed signal depends on grass LAD and LAI [83], in combination with observation and illumination angles. Figure 15 shows simulated grass, see Figure 15a and soil, see Figure 15b, RF observed at nadir for different SWC values. As can be seen, the largest changes occur when SWC is low (less than 20%), which can concur with LAI_grass values between 0.1 and 2.2 during phenophase 1.

Regarding the temporal variability of the ecosystem and the effect on the modeling, phenophases that were defined according to the time series analysis over MODIS-NDVI data were found to be appropriate to represent the phenology of the TGE under investigation. However, large overlaps between phenophases, especially between phenophase 2 and 3, reveals the strong inter-annual variability of this ecosystem. LUT design has been oriented so that each LUT is capable of representing RF in each of the surrounding transition periods with no loss of accuracy. In addition, other sources of information could improve the phenological characterization of the ecosystem, such as phenocams [41,84] or spectral indices adapted to phenological analyses [85]. Enhanced information would contribute to improving the definition of the phenophases and the classification of field remote sensing data in each of them.

5. Conclusions

In this study, the FLIGHT model has been adapted through a simplified coupling strategy with PROSAIL to improve the simulation of the optical properties of a Mediterranean TGE. The evaluation of the PROSAIL model against proximal sensing observations of the grass layer revealed a decrease of model performance under the presence of standing death material. The evaluation of PROSAIL and PROSAIL+FLIGHT at the ecosystem scale revealed that the uncertainties found in the modeling of the grass background were a major contributor to the errors in the simulation of ecosystem RF. In the growing season, PROSAIL can predict similar RF as those of the coupled model for high view and zenith angles. However, PROSAIL+FLIGHT performs better when properties of grass and trees differ and represents ecosystem optical responses at larger view angles more realistically. RTM evaluation against airborne and satellite imagery revealed the challenge of representing the average state of ecosystems with a high spatial variability from a reduced number of vegetation samples. Scaling uncertainties explain, in part, the better performance of PROSAIL at a proximal sensing scale than PROSAIL+FLIGHT at a remote sensing level. This work proves the suitability of proximal sensing for validating RTM instead. PROSAIL+FLIGHT and the seasonally-based modeling framework developed in this work could improve the estimation vegetation biophysical variables at different scales and angular configurations in TGE either via RTM inversion [78,86] or empirical models [69,70] in the future.