Modeling Gross Primary Production (GPP) of a Mediterranean Grassland in Central Spain Using Sentinel-2 NDVI and Meteorological Field Information

Abstract

Mediterranean grasslands provide different ecosystems and social and economic services to the Mediterranean basin. Specifically, in Spain, pastures occupy more than 55% of the Spanish surface. Farmers and policymakers need to estimate the Gross Primary Production (GPP) to make sustainable management of these ecosystems and to study the role of grasslands acting as sinks or sources of Carbon in the context of climate change. High-frequency satellites, such as Sentinel-2, have opened the door to study GPP with a higher spatial and lower revisit time (10 m and 5 days). Therefore, the overall objective of this research is to estimate an ecosystem light use efficiency (eLUE) GPP model for a Mediterranean grassland in central Spain using Sentinel-2 NDVI Normalized Difference Vegetation Index (NDVI), complemented with meteorological information at the field scale for a relatively long period (from January 2018 to July 2020). The GPP models studied in this research were the MODIS GPP product, as well as the four eLUE models built with MODIS or Sentinel-2 NDVI and complemented by the inclusion of minimum temperature (T_min) and soil water content (SWC). The models were validated through the GPP obtained from an eddy-covariance flux tower located in the study site (GPP_T). Results showed that the MODIS GPP product underestimated the GPP_T of the grassland ecosystem. Besides this, the approach of the eLUE concept was valid for estimating GPP in this Mediterranean grassland ecosystem. In addition, the models showed an improvement using Sentinel-2 NDVI compared to MODIS GPP product and compared to the models that used MODIS NDVI due to its higher spatial and temporal resolution. The inclusion of T_min and SWC was also a determinant in improving GPP models during winter and summer periods. This work also illustrates how the main wind directions of the study area must be considered to appropriately estimate the footprint of the eddy covariance flux tower. In conclusion, this study is the first step to efficiently estimating the GPP of Mediterranean grasslands using the Sentinel-2 NDVI with complementary meteorological field information to make the management of these ecosystems sustainable.

1. Introduction

Mediterranean grasslands provide different ecosystem, social, and economic services to the Mediterranean basin, occupying approximately 1.3 million km² of the terrestrial surface [1,2]. They help to maintain biodiversity, improve soil protection, and lead to important water and climate interactions [3,4]. They are also determinants for livestock production, especially sheep and cattle [5,6]. As an example, pastoral exploitation is the first economic activity in the Spanish mountains, constituting the economic base of the rural environment. As a whole, forest and agricultural pastures occupy more than 27 million hectares, i.e., more than 55% of the Spanish surface [7].

Mediterranean grassland dynamics are conditioned by the influence of the Mediterranean climate, which results in large inter-annual and seasonal variability. This leads to a complex pattern of CO₂ exchange between the atmosphere and the ecosystem. On the one hand, the assessment of CO₂ balance dynamics through photosynthesis (Gross Primary Production, GPP) and ecosystem respiration (Reco) has become a fundamental issue due to the potential impact that it may have on climate change [8,9,10]. Thus, it is important to estimate and model GPP to quantify whether Mediterranean grasslands act as sinks or sources of CO₂ over the years depending on the various factors that can interact with them, such as climate, disturbing agents (plagues, fires, natural disasters…), and human factors, among others [8,11]. Furthermore, in the Mediterranean environment, all the atmospheric general circulation models predict an increase in the frequency and intensity of drought periods as a consequence of climate change [12]. Hence, water stress in vegetation is expected to increase, leading to changes in biomass production as well as in the structure and spatial distribution of plant ecosystems [13,14,15]. On the other hand, estimating GPP is essential for farmers to make sustainable management of ecosystems, especially grasslands, combining social, economic, and environmental factors [16,17,18]. This motivates the study of grasslands from a dynamic point of view, that is, analyzing its temporal evolution at different spatial scales.

In this sense, remote sensing (RS) is accepted as the most powerful tool for studying ecosystems at different spatial and temporal scales [19,20]. In the last two decades, several research studies have assessed the use of RS to estimate GPP in different ecosystems [21,22,23]. Most RS models used to estimate GPP are based on the Light Use Efficiency concept proposed by Monteith [24] (Equation (1)):

$$GPP=\epsilon \times APAR=\epsilon \times fPAR\times PAR$$

(1)

where ε is the light use efficiency conversion coefficient, and the Absorbed Photosynthetically Active Radiation (APAR) is the product of the Fraction of Photosynthetically Active Radiation (fPAR) by the Photosynthetically Active Radiation (PAR).

On the other hand, several researchers have evidenced that GPP is related to spectral vegetation indices. Ref. [25] proposed to multiply the chlorophyll content by PAR in the Greenness and Radiation model, obtaining better results in estimating crops’ GPP. After that, Peng et al. [26] applied this model using several spectral indices to estimate crops’ GPP, obtaining good results. Ma et al. [27] proposed the concept of ecosystem light use efficiency (eLUE), that is, the product of fPAR by ε, which is related to spectral vegetation indices (Equation (2)):

$$GPP=eLUE\times PAR$$

(2)

They proposed that eLUE was a function of the Enhanced Vegetation Index (EVI) and that multiplying EVI by the PAR at the Top of the Atmosphere (PAR_TOA) was appropriate to obtain a robust estimation of GPP in savannah ecosystems. In this research, we hypothesize that this approach, along with the information from high-frequency satellites, could be adequate for estimating GPP in Mediterranean grasslands.

High-frequency satellites, such as Sentinel-2, have opened the door to the study of ecosystems with a high spatial and temporal resolution [28,29]. Although extensive research exists about estimating GPP in grasslands with moderate-resolution sensors, such as MODIS [30,31,32], the estimation of GPP in this ecosystem using Sentinel-2 images is relatively new. In this sense, Sakowska et al. [33] highlighted the potential of the Sentinel-2 satellite to monitor biophysical parameters in a dynamic subalpine grassland ecosystem. More recently, Misra et al. [34] reviewed the status of phenological research using Sentinel-2 data in different ecosystems. They showed that Sentinel-2 data is useful for capturing the spatial and temporal variations in grasslands due to management practices and climate. Besides, Wang et al. [35] studied the aboveground biomass of two pastures in central Oklahoma using Sentinel-1, Sentinel-2, and Landsat images, and they showed that the combination of these data improved biomass estimation to improve management in this ecosystem. Similarly, Sharma [36] developed the Dominant Species-Physiognomy-Ecological (DSPE) classification system for differentiating plant ecological communities using high-spatial-resolution Sentinel-2 images, demonstrating the effectiveness of Sentinel-2 in ecological studies at a regional scale. However, only a few works have estimated GPP using Sentinel-2 data in Mediterranean environments. Chrysafis et al. [37] evaluated the use of Sentinel-2 data to estimate the growing stock volume of a Mediterranean forest, obtaining better results than those obtained by Landsat imagery. Cerasoli et al. [38] achieved good estimations of GPP of Mediterranean grasslands with hyperspectral measurements, which were resampled to reproduce the spectral bands of both Sentinel-2 and Landsat 8. In addition, Gómez-Giráldez et al. [39] used Sentinel-2 and MODIS data for estimating Net Primary Production (NPP) at different scales of a holm oak Mediterranean forest, obtaining promising results. Recently, Spinosa et al. [40] upscaled the GPP of a flux tower located in a Mediterranean wetland using Sentinel-2 data. They showed that their GPP model predictions improved those obtained from using the MODIS product. However, they used only one year, and further research is needed to estimate GPP through longer periods. In this sense, this study proposed using Sentinel-2 data to estimate the GPP of a Mediterranean grassland for three growing seasons.

GPP derived from remote sensing data must be evaluated at the field scale using eddy covariance (EC) measurements from flux towers. EC technique has allowed quantifying more accurately the flux of CO₂ between the atmosphere and ecosystems (Net Ecosystem exchange, NEE) and deriving GPP and Reco through different techniques at different spatial-temporal scales [41,42,43], as Equation (3) shows:

$$NEE=Reco-GPP$$

(3)

The overall objective of this research is to estimate an ecosystem light use efficiency GPP model for a Mediterranean grassland using Sentinel-2 NDVI and meteorological field information. The specific objectives are the following:

(1)

To validate the remote sensing model with measurements from an EC flux tower in terms of quantity and seasonal and interannual dynamics.

(2)

To compare the model at different spatial and temporal resolutions using Sentinel-2 and MODIS data.

2. Materials and Methods

2.1. Study Area

The Mediterranean grassland is located close to the Forest of La Herreria, in the village of El Escorial in the province of Madrid (central Spain) at 40°34′56″ N and 4°08′16″ W and at an altitude of 920 m.a.s.l. (Figure 1). The climate of this site corresponds to a Mediterranean climate Csa according to Köppen climate classification [44] with hot and dry summers and cold winters, with an average annual precipitation of around 500 mm and an average annual temperature of 12.5 °C. The soil of the area is acid soil classified as a Haplic Cambisol according to the World Reference Base of Soil Resources [45]. The vegetation is composed of supra-Mediterranean mesophytic grasslands with some trees of Fraxinus angustifolia Vahl. and some bushes of Rosa sp. and Rubus sp. [46]. The study area is subject to cattle grazing throughout the year.

2.2. Flux Tower Data

In the study site, there is an eddy covariance flux tower that belongs to the Guadarrama Meteorological Network [47,48,49].

The flux tower provides data from 2017 to the present at a frequency of 10 Hz of several micrometeorological variables, including those for calculating CO₂ flux (Net Ecosystem exchange, NEE). This study collected and analyzed nearly three years of variables from January 2018 to July 2020 (

Table 1. Variables obtained for the study from January 2018 to July 2020 in the Mediterranean grassland: Gross Primary Production obtained from the flux tower (GPP_T). Minimum air temperature (T_min), Soil water content at 4 cm depth (SWC), Photosynthetically Active Radiation from the Top of the Atmosphere (PAR_TOA), NDVI obtained from Sentinel-2 (Sentinel-2 NDVI), NDVI obtained from MODIS (MODIS NDVI), Gross Primary Production derived from the MODIS product (MODIS_GPP).

Variable	Source	Frequency	Daily Average	5-Day Average	8-Day Average
GPP_T	Flux tower	30 min	Yes	Yes	Yes
Tmin	Flux tower	30 min	Yes	Yes	Yes
SWC	Flux tower	30 min	Yes	Yes	Yes
PARTOA	Solar Radiation	daily	Yes	Yes	Yes
SENTINEL-2 NDVI	Sentinel-2	5-day	No	Yes	No
MODIS NDVI	MODIS	8-day	No	No	Yes
MODIS_GPP	MODIS	8-day	No	No	Yes

). Firstly, raw files from the IRGASON (Campbell Scientific, Inc., Logan, UT, USA) every three hours were obtained for the study period, and the raw files were divided into 30-min files through the CardConvert 1.9 software (Campbell Scientific, Logan, UT, USA; hereafter CSI). Secondly, these files were processed with EddyPro 7.0.9 software (Li-Cor, Inc., Lincoln, NE, USA) to obtain the CO₂, water vapor, and sensible heat fluxes every 30 min for the study period. Half-hourly averaged fluxes were excluded when more than 10% of the raw data records were missing for any of the three components of wind velocity, temperature, and/or CO₂ concentration. To ensure the quality of EC flux data, we utilized quality control flags 0 and 1 from the outputs of EddyPro [50]. Data gaps were filled using the marginal distribution sampling technique described by Reichstein et al. [51]. Flux partitioning into Ecosystem Respiration (Reco) and Gross Primary Production (GPP) was carried out according to the method by Lasslop et al. [52]. Finally, minimum air temperature (T_min), soil water content at 4 cm depth (SWC), and tower GPP (GPP_T) were averaged daily and then every 5 days to match the temporal resolution of remote sensing images. GPP_T was then smoothed with a Savitzky-Golay filter [53] to reduce noise using a second-order polynomial with a smoothing window of seven and using the same filter as the MODIS and Sentinel-2 NDVI.

2.3. Photosynthetically Active Radiation from the Top of the Atmosphere (PAR_TOA)

The extraterrestrial solar radiation (R_a, [MJ m⁻² day⁻¹]) was obtained daily as follows [54], Equation (4):

$$R_a=\frac{24\ast 60}{\pi } G_{sc}{d}_r\left[{\omega }_s sin\left(\phi \right)\sin \left(\delta \right)+cos\left(\phi \right)\cos \left(\delta \right) \mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_s\right)\right]$$

(4)

where:

G_sc is the solar constant = 0.0820 [MJ m⁻² min⁻¹], dr is the inverse relative distance Earth-Sun [rad], ${\omega }_s$ is the sunset hour angle [rad], $\phi$ is the latitude [rad] and $\delta$ is the solar declination [rad].

The Photosynthetically Active Radiation from the Top of the Atmosphere (PAR_TOA) was obtained as 40% of R_a [27,55]. Then, it was averaged every 5 days to match the temporal resolution of the remote-sensing images.

2.4. Remote Sensing Data

2.4.1. MODIS

MODIS images for the product MOD09Q1 Version 6.1 were downloaded from the NASA Land Processes Distributed Active Archive Center (LP DAAC) webpage (https://lpdaac.usgs.gov/ accessed on 10 September 2023). This product provides an estimate of the surface spectral reflectance for the Near Infrared (NIR) and Red (R) bands at a spatial resolution of 250 m and a temporal resolution of eight days. These bands were obtained for the calculation of the NDVI for the eddy covariance flux tower pixel from January 2018 to July 2020 (117 images), following Equation (5):

$$NDVI=\frac{(NIR-R)}{(NIR+R)}$$

(5)

Then, a Savitzky-Golay filter [53] was applied to reduce noise and smooth the time series using a second-order polynomial with a smoothing window of seven.

The product MOD17A2HGF Version 6.1 Gross Primary Productivity (GPP) product is a cumulative 8-day composite of values with a 500-m pixel size [56]. This product was downloaded for the pixel of the flux tower with the AppEEARS Tool [57].

2.4.2. Sentinel-2

Remote sensing images were obtained from the Sentinel-2 MultiSpectral Instrument for the tile 30TVK from January 2018 to July 2020 (186 images). Level-2A images (atmospherically corrected) were downloaded from (https://scihub.copernicus.eu/dhus/#/home accessed on 10 September 2023). Sentinel-2 images with the highest spatial resolution (10 m) and a temporal resolution of 5 days were used. The bands Red (R, b4: 646–685 nm) and Near InfraRed (NIR, b8: 767–908 nm) were used for the calculation of the NDVI following Equation (5).

Then, a Savitzky-Golay filter [53] was applied to reduce noise and smooth the time series using a second-order polynomial with a smoothing window of 7. Sáenz et al. [58] evidenced that this type of filter was adequate for smoothing Sentinel-2 NDVI in a variety of croplands and natural ecosystems. Thus, this filter was also applied to the MODIS and Sentinel-2 NDVI of the grassland ecosystem.

2.5. Model Construction

Models were built according to the ecosystem light use efficiency concept (eLUE) (Equation (2)). The first model (Model 1) was similar to that proposed by Ma et al. [27], but using the NDVI obtained from MODIS or Sentinel-2 instead of the MODIS EVI, as Equation (6) shows:

$$GPP=f\left(NDVI\right)\times {PAR}_{TOA}$$

(6)

The second one (Model 2) was a new proposal based on a multiple linear regression model with three parameters: (1) NDVI × PAR_TOA, (2) T_min, and (3) SWC, as Equation (7) shows:

$$GPP=a_0+b\times \left(NDVI\times {PAR}_{TOA}\right)+c\times T_{min}+d\times SWC$$

(7)

2.6. Statistical Methods

The two models were built with data from the years 2018 and 2019 and validated with the year 2020 through the comparison of the observed GPP from the flux tower data versus the estimated GPP by the models using the coefficient of determination (R²).

The first type of model (Equation (5)) was built based on the relationship between the ecosystem light use efficiency (eLUE) and the MODIS or Sentinel-2 NDVI (Equation (2)).

The multiple linear regression model (Model 2) was adjusted according to the least square method and an analysis of variance (ANOVA) was performed to evaluate the significance of the model using the F-statistic. The significance of the model’s coefficients was evaluated by the T-statistic.

The coefficient of determination (R²) was used to measure the proportion of the variance explained by the models. The Root Mean Square error (RMSE) is calculated for n different observations as the square root of the mean of the squares of the deviations, as Equation (8) shows:

$$RMSE=\sqrt{\frac{\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{n}}$$

(8)

where $y_i$ is the observed variable, ${\widehat{y}}_i$ is the estimated variable, and n is the number of observations.

The statistical analyses were performed through Statgraphics Centurion XIX (StatPoint Technologies Inc., Warrenton, VA, USA) and MATLAB (R2022b, MathWorks, Natick, MA, USA).

3. Results

3.1. Dynamics of GPP_T in the Mediterranean Grassland

Figure 2 shows the dynamics of tower GPP (GPP_T) and of the soil water content at 4 cm depth (SWC). As it is observed, GPP_T had a minimum during winter (around 2 g C m⁻² day⁻¹) due to minimum temperatures. Then, it increased during spring, having a clear maximum in May–June (around 10–12 g C m⁻² day⁻¹). After that, during summer, GPP_T shows another minimum (around 2 g C m⁻² day⁻¹) due to the lack of SWC. During autumn, when SWC recovers, GPP_T shows a secondary growing season. According to the classification of the Spanish State Meteorological Agency (AEMET), the year 2018 was classified as humid, while the year 2019 was classified as very dry from January to October [59].

3.2. MOD17A2HGF v 6.1 Gross Primary Productivity (GPP) Product (MODIS_GPP)

Figure 3 shows the MODIS_GPP with a spatial resolution of 500 m and a temporal resolution of 8 days for the pixel of the tower versus the GPP_T. As observed, MODIS_GPP underestimated GPP throughout the study period except for August and September in 2018 and 2019. However, MODIS_GPP was able to capture the seasonal and interannual dynamics of GPP well. The R² for all the time series was 0.80, and the RMSE was 2.24 g C m⁻² day⁻¹.

3.3. eLUE Models Using MODIS

The MODIS models were built for the pixel of the tower with a spatial resolution of 250 m and a temporal resolution of 8 days.

The first model (GPP_M1) was based on the relationship between eLUE (GPP_T/PAR_TOA) and MODIS NDVI, as Figure 4 shows:

The linear regression model (n = 92, R² = 0.347, ANOVA p-value = 0.00) obtained is shown in Equation (9):

$$eLUE=1.152\times NDVI-0.162$$

(9)

Therefore, the GPP_M1 model for obtaining GPP is shown in Equation (10):

$$GPP\_M1=\left(1.152\times \left(NDVI-{Soil}_{bg}\right)-0.162\right)\times {PAR}_{TOA}$$

(10)

Being ${Soil}_{bg}$ the soil background coefficient (${Soil}_{bg}=0.161)$ obtained from the linear regression model between MODIS NDVI and GPP. The soil background coefficient was the intercept of the line, that is, the value of NDVI when GPP was zero. Although the linear regression model showed a poor adjustment, we believe that this coefficient was important to compare the model proposed by Ma et al. [27].

The time series obtained for this model (GPP_M1) from January 2018 to July 2020 versus the tower GPP (GPP_T) is shown in Figure 5. As it is observed, the model underestimated the GPP during the study period. The RMSE for the study period was 2.78 g C m⁻² day⁻¹.

The validation of the model was performed using the year 2020, obtaining an R² = 0.78 and a RMSE of 2.66 g C m⁻² day⁻¹. Figure 6 shows the values of GPP_M1 versus the values of GPP_T:

The second model (GPP_M2) was a multiple linear regression model with three parameters: (1) NDVI × PAR_TOA, (2) T_min, and (3) SWC, as Equation (7) shows. The result of the analysis of variance (ANOVA) is shown in

Table 2. Analysis of variance (ANOVA) of the multiple regression model: G$PP\_M2=a_0+b\times \left(NDVI\times {PAR}_{TOA}\right)+c\times T_{min}+d\times SWC$ from January 2018 to December 2019 in the Mediterranean grassland.

Source	Sum of Squares	Degrees of Freedom	Mean Square	F-Statistic	p-Value
Model	560.466	3	186.822	325.84	0
Residual	50.4555	88	0.573
TOTAL	610.921	91

.

Table 3. Estimated coefficients of the multiple regression model: G$PP\_M2=a_0+b\times \left(NDVI\times {PAR}_{TOA}\right)+c\times T_{min}+d\times SWC$, standard error of the coefficients and T-statistic from January 2018 to December 2019 in the Mediterranean grassland.

Coefficient	Estimation	Standard Error	T-Statistic	p-Value
Constant (a0)	−5.148	0.541	−9.522	0
NDVI × PARTOA (b)	1.214	0.047	25.951	0
Tmin (c)	0.091	0.025	3.609	0.001
SWC (d)	12.859	1.617	7.951	0

shows the estimated coefficients of the model and their significance according to the T-statistic:

Therefore, the GPP_M2 model for obtaining GPP with a significative coefficient of determination (R²) of 0.91 is shown in Equation (11):

$$GPP\_M2=-5.148+1.214\times \left(NDVI\times {PAR}_{TOA}\right)+0.091\times T_{min}+12.859\times SWC$$

(11)

The time series obtained for this model (GPP_M2) from January 2018 to August 2020 versus the tower GPP (GPP_T) is shown in Figure 7. As it is observed, the model slightly underestimated the spring GPP maxima in 2018 but it slightly overestimated the maxima in 2019 and 2020. In addition, the model overestimated the GPP minima during summer and during winter showed similar values to the GPP_T. The RMSE for all the time series is 0.94 g C m⁻² day⁻¹.

The validation of the model was obtained using the year 2020 with an R² = 0.75 and a RMSE of 1.42 g C m⁻² day⁻¹. Figure 8 shows the values of GPP_M2 versus the values of GPP_T:

3.4. Footprint Estimation for Sentinel-2

Sentinel-2 overpassed our study area at 10:56 UTC, therefore, we have selected the average of the footprint (x_peak according to Eddypro software) of all days at 11.00 UTC (110 m) calculated by the EddyPro 7.0.9 software (Li-Cor, Inc., Lincoln, NE, USA) from January 2018 to July 2020. A circle of 110 m was built centered in the tower. In addition, the histogram of wind directions at 11.00 UTC was built from January 2018 to August 2020 to consider the main directions of the wind at this time (Figure 9). The wind rose during the study period at 11.00 UTC as shown in Figure 10. As can be observed, the main wind direction is from the southeast. Finally, the pixels inside the green area, which was delimited by the main wind directions, were selected for calculating the spatial average of NDVI (Figure 11).

3.5. eLUE Models Using Sentinel-2

The models were built as the MODIS ones for the spatial average of the pixels mentioned above, which had a temporal resolution of 5 days. The first model was based on the relationship between eLUE (GPP/PAR_TOA) and NDVI, as Figure 12 shows:

The linear regression model showed a higher R² value (R² = 0.51) than the MODIS one, and it is shown in Equation (12):

$$eLUE=0.886\times NDVI+0.0395$$

(12)

Therefore, the GPP_S1 model for Sentinel-2 for obtaining GPP is shown in Equation (13):

$$GPP\_S1=\left(0.886\times \left(NDVI-{Soil}_{bg}\right)+0.0395\right)\times {PAR}_{TOA}$$

(13)

Being ${Soil}_{bg}$ the soil background coefficient (${Soil}_{bg}=-0.05)$ obtained from the linear regression model between Sentinel-2 NDVI and GPP. The soil background coefficient was the intercept of the line, that is, the value of NDVI when GPP was zero.

The time series obtained for this model (GPP_S1) from January 2018 to July 2020 versus the tower GPP (GPP_T) is shown in Figure 13. As it is observed, this model adjusted better than GPP_M, underestimating the spring GPP maxima and overestimating the GPP minima during summer and winter. The RMSE for the study period was quite lower than the MODIS model (GPP_M1), and it was 0.93 g C m⁻² day⁻¹.

The validation of the model was obtained using the year 2020, obtaining an R² = 0.87 and a RMSE of 0.99 g C m⁻² day⁻¹. Figure 14 shows the values of GPP_S1 versus the values of GPP_T:

The second model (GPP_S2) was a multiple linear regression model with three parameters as the second MODIS model (GPP_M2): (1) NDVI × PAR_TOA, (2) T_min, and (3) SWC, as Equation (14) shows. The result of the analysis of variance (ANOVA) is shown in

Table 4. Analysis of variance (ANOVA) of the multiple regression model: G$PP\_S2=a_0+b\times \left(NDVI\times {PAR}_{TOA}\right)+c\times T_{min}+d\times SWC$ from January 2018 to December 2019 in the Mediterranean grassland.

Source	Sum of Squares	Degrees of Freedom	Mean Square	F-Statistic	p-Value
Model	901.512	3	300.504	476.88	0
Residual	89.4805	142	0.630145
TOTAL	990.992	145

.

Table 5. Estimated coefficients of the multiple regression model: G$PP=a_0+b\times \left(NDVI\times {PAR}_{TOA}\right)+c\times T_{min}+d\times SWC$, standard error of the coefficients and T-statistic from January 2018 to December 2019 in the Mediterranean grassland.

Coefficient	Estimation	Standard Error	T-Statistic	p-Value
Constant (a0)	−2.160	0.413	−5.236	0
NDVI × PARTOA (b)	1.100	0.035	31.542	0
Tmin (c)	0.097	0.019	4.980	0
SWC (d)	3.846	1.314	2.926	0.004

shows the estimated coefficients of the model and their significance according to the T-statistic:

Therefore, the second model for obtaining GPP with a coefficient of determination (R²) of 0.91 is shown in Equation (14):

$$GP{P}_{S2}=-2.160+1.100\times \left(NDVI\times {PAR}_{TOA}\right)+0.097\times T_{min}+3.846\times SWC$$

(14)

The time series obtained for this model (GPP_S2) from January 2018 to August 2020 versus the tower GPP (GPP_T) is shown in Figure 15. As observed, the model underestimated the spring GPP maxima in 2018 and 2019. In addition, the GPP minima during summer and winter are similar to the GPP_T. The RMSE for all the time series is 0.85 g C m⁻² day⁻¹.

The validation of the model was obtained using the year 2020 with an R² = 0.87 and a RMSE of 1.06 g C m⁻² day⁻¹. Figure 16 shows the values of GPP_S2 versus the values of GPP_T:

3.6. Summary of the Models

Table 6. Coefficient of determination (R²) and Root Mean Square Error (RMSE) of the five models studied in the Mediterranean grassland for the year 2020 and for all the time series (TS).

Model	R2 (2020)	RMSE (2020)	R2 (TS)	RMSE (TS)
MODIS_GPP	0.81	2.49	0.80	2.24
GPP_M1	0.78	2.66	0.74	2.78
GPP_M2	0.75	1.42	0.87	0.94
GPP_S1	0.87	0.99	0.89	0.93
GPP_S2	0.87	1.06	0.90	0.85

summarizes the statistics of the five models studied in this research. The R² of the models were relatively high (above 0.7) both for the year 2020 and for all the time series. In both cases, the highest RMSE occurred with the MODIS_GPP product and with the GPP_M1, while the lowest RMSE occurred with the Sentinel-2 models (GPP_S1 and GPP_S2).

4. Discussion

4.1. Overall Discussion

In this research, the assessment of the dynamics of GPP and the estimation of a GPP model through remote sensing data for a Mediterranean grassland of Central Spain has been carried out. The dynamics of the GPP from an eddy covariance flux tower have been studied during three growing seasons (2018–2020). As it is observed (Figure 2), GPP_T depends strongly on SWC having two growing seasons during the year, the main growing season during spring and a secondary growing season during autumn. On one hand, the length and the maximum of the spring season depend on the SWC [60,61]. In fact, during 2018, the growing season was the maximum, and it was extended until mid-June due to the precipitation in May, which entailed an increase in SWC. On the other hand, the maximum of autumn also depends on the recovery of SWC, and the extension depends on the arrival of the minimum temperatures at the end of autumn [31,62].

Results showed that MODIS_GPP underestimated GPP throughout the study period except for the minima of summer periods in 2018 and 2019. However, MODIS_GPP was able to capture the seasonal and interannual dynamics of GPP well. Although the MODIS GPP product has been improved in recent years, the sources of error of this product could be the spatial resolution that is not enough to assess finer ecosystem processes [63,64] or the parameters used in MODIS GPP [65].

The proposed models were based on the eLUE approach proposed by Ma et al. [27] for different savannah ecosystems in Australia, using NDVI instead of EVI. Huete et al. [66] concluded that NDVI showed problems of saturation and a lower dynamic range over the more humid forested sites. On the other hand, they showed that NDVI could present a higher range in values over the semiarid sites, as in the case of our study area. In fact, we obtained the Sentinel-2 EVI time series for our pixels, and it presented more noise than the time series of NDVI. Therefore, the use of NDVI instead of EVI was justified in our case. However, EVI should be considered in other ecosystems due to its higher sensitivity to canopy structural variations and improved resistance to atmospheric conditions [66,67]. Further research is needed to evaluate the use of other indices in the models, such as the Normalized Difference Water Index (NDWI) [68] or Spectral Shape Indices (SSI) [69].

PAR_TOA was used to compare the same model proposed by Ma et al. [27] in Australian savannahs in Mediterranean grasslands. This research has already discussed the use of PAR_TOA instead of the Photosynthetic Active Radiation from the top of the canopy (PAR_TOC), and they have obtained better results with PAR_TOA. Therefore, the use of PAR_TOA in this research was justified mainly for two reasons: (1) To apply and compare the same model as Ma et al. [27] in our Mediterranean grassland, (2) we do not have the measure of PAR_TOC in our site.

In this work, the eLUE approach is valid for a Mediterranean grassland ecosystem. Firstly, in Mediterranean ecosystems, the light use efficiency conversion coefficient (ε) shows a strong seasonality [31,70] as in savannahs. Therefore, it seems that, in Mediterranean grasslands, the hypothesis that FAPAR and ε are synchronized and related both with spectral indices, in our case NDVI, is also fulfilled. This fact is due to the patterns of Mediterranean grasslands that are conditioned by the dynamics of the Leaf Area Index, which depends on environmental variables, such as SWC [71,72]. However, in this research, the correlation between eLUE and MODIS NDVI was much lower (R² = 0.35) than the one obtained by Ma et al. [27] between eLUE and MODIS EVI in savannahs (R² = 0.81). The reason was probably due to the different ecosystems that have been compared. In fact, in one type of semi-arid Savanna (Acacia woodland) studied by Ma et al. [27], this relationship was also lower. The Acacia semi-arid woodland is more similar to our Mediterranean grassland ecosystem with small bushes, trees, and herb species than the other types of studied savannahs. In this type of savannah, the eLUE model failed to capture the maxima and the minima of the GPP, but in our eLUE MODIS model 1, the model infra-estimated the values during all the study period (GPP_M1, Figure 4). The model improved clearly (with a lower RMSE) when SWC and T_min were considered (GPP_M2, Figure 6). Pasture phenology was linked to the dynamics of soil water availability, water being the main limiting factor during the growing cycle. Thus, during summer, annual herb species died, and only green trees and bushes could contribute to GPP due to their longer roots that could absorb water from deeper layers of the soil [73,74]. This fact was the reason for a value greater than 0 g C m⁻² day⁻¹ during this season. Therefore, the inclusion of SWC in the model helped to capture the minima during summer. This factor was also relevant in the savannahs of Australia, especially in the semi-arid ones [27,75]. In the Mediterranean climate, it is known that SWC is an essential parameter of GPP (Figure 2) [72,76]. On the other hand, temperature is only a limiting factor during winter [31,77], although some species could be active during this season [61], but the inclusion of this factor in the model contributed to capturing the GPP minima during winter. In addition, both factors affected the capture of the maximum GPP during spring.

The ability of Sentinel-2 indices to model the GPP of grasslands has been demonstrated. Lin et al. [78] also concluded that several Sentinel-2 vegetation indices, especially those based on red-edge data, followed the phenological variations of two grasslands in Australia. The higher temporal (5 days) and spatial resolution (10 m) of Sentinel-2 have served to improve the model results. In fact, only by increasing these resolutions did the GPP_S1 adjust better than the MODIS GPP_M1. In this sense, Wang et al. [35] showed that the combination of high spatial and temporal resolution sensors, such as Sentinel-2, improved biomass estimation in grassland ecosystems. Besides, Misra et al. [34] also showed the potential of Sentinel-2 compared with other remote sensing data to assess the phenology of ecosystems. In addition, Phiri et al. [79] showed that Sentinel-2 has a positive impact on land cover/use monitoring due to its higher temporal and spatial resolution.

In this work, the main wind directions of the satellite’s overpassed time were taken into account when selecting the pixels for the study. This footprint estimation approach was done due to the location of the flux tower, which is influenced by mesoscale thermally driven flows (mountain breezes) due to the presence of the Guadarrama Mountains, located quite close to the station [48,49]. According to these studies, the dynamics of the planetary boundary layer and these breezes influenced the CO₂ mixing ratio of the study area. Therefore, it was important to consider the main direction of wind when selecting the Sentinel-2 pixels that would influence the observed and estimated GPP. In this sense, Aguirre-García et al. [80] followed a similar strategy; they also considered the wind direction and the satellite overpassed time to calculate the footprint of the eddy covariance tower for modeling evapotranspiration in a Mediterranean olive orchard.

The validation of the Sentinel-2 models for the year 2020 showed that both models showed very similar R² and RMSE values. This fact is due to the model having only been validated until mid-July 2020, and the inclusion of the meteorological parameters would have more weight during all the summer and winter periods.

4.2. Study Limitations and Further Research

The GPP obtained from the eddy covariance method has been considered a ground-truth measurement. This method might have some limitations due to assumptions, physical phenomena, instrument problems, and specifics of terrain and setup [43]. Anyway, eddy covariance estimations are generally considered in the literature as the best data to validate models [81].

Although Sentinel-2 models showed good performance over almost three years, future research is needed to obtain a longer time series of GPP to make these site-specific models more robust due to the inter-annual variations of the Mediterranean climate. This fact is especially important to obtain accurate models under different climate change scenarios. In addition, adjustments in the model are necessary to better capture maxima values in spring and minima values in summer and winter, which could be crucial for ecological and economic grassland management. In this sense, it could be very useful to obtain a GPP model based only on remote sensing data to avoid dependence on site field variables. Further research is needed to try new complex models, such as machine learning techniques, to achieve a better adjustment. In addition, with longer time series, statistical methods such as time series analysis could be applied for also assessing the temporal evolution of variables through years, i.e., inter-annual variation and their dynamics [82,83], and this fact could enhance GPP forecasting.

The models studied here are valid for a single study area. In order to test the models, they must be applied to other grasslands with similar climatic conditions. Secondly, to upscale the model, they must be applied to other Mediterranean grasslands with different ecological and climatic characteristics. In this sense, the main limitation will be the existence of numerous flux towers to make a robust validation. Therefore, the Fluxnet [84] and ICOS [85] communities offer an excellent opportunity to try this model in other grassland sites with an Eddy Covariance flux tower.

Finally, it is essential to transfer these models to farmers and policymakers. It would be important to make an application for mobile devices that could indicate the status of the grassland in real time by remote sensing images and the development of robust models that could estimate and forecast the Net Primary Production and the yield of the pasture. However, this fact should require training in aspects such as remote sensing techniques and information management by these stakeholders.

5. Conclusions

This research has studied the dynamics of Gross Primary Production (GPP) obtained from an eddy covariance flux tower of a Mediterranean grassland in central Spain for almost three years. While minimum temperature was limiting for growth during winter, the soil water content was the limiting factor during the rest of the year, especially in summer. It must be noted that the presence of a secondary growing season during autumn in the Mediterranean grasslands is associated with the recovery of soil water content after summer and its extension limited by the arrival of low temperatures.

In this work, the estimation of the GPP of this grassland through remote sensing data and meteorological information has been carried out through the approach of the eLUE proposed by Ma et al. [27]. This approach, which was originally applied to Australian savannahs, has been demonstrated to be valid for this Mediterranean grassland ecosystem. This model was initially constructed with MODIS data. In this research, MODIS data were used to compare results with those of other studies. However, the approach has been extended to the higher spatial and temporal resolution of Sentinel-2, demonstrating the ability of this satellite to improve the GPP site-specific models over Mediterranean grasslands.

The NDVI has been used instead of the EVI due to it could present a higher range in values over this Mediterranean site. The inclusion of site meteorological variables (minimum temperature and soil water content) has improved the model, especially in the summer and winter periods when these factors were constraining GPP. In addition, the main wind directions must be considered when estimating the footprint of the eddy covariance tower and selecting the appropriate Sentinel-2 pixels to build the model.

This study is the first step toward a model for the GPP of Mediterranean grasslands with Sentinel-2 data to make the management of these ecosystems sustainable. It would be essential to transfer these models to farmers and policymakers of these ecosystems by the calculation of the Net Primary Production and the yield of the pasture.