1. Introduction
Grasslands are among the most widespread vegetation types worldwide, covering between 14% and 26% of the earth’s surface [
1,
2]. They play an important role in the carbon (C) cycle, as they store about 23% of the global soil C [
3] and, given the contribution of grasslands in mitigating climate change, it is important to monitor their C sequestration potential.
The Eddy Covariance (EC) method permits to measure, at an ecosystem scale, the atmosphere-ecosystem exchange of water, energy and CO
2 fluxes [
4]. It provides a reliable direct measurement of different gas compounds together with meteorological variables, at high temporal detail, making it possible to ascertain the influence of climate drivers or other disturbances/stress factors on fluxes at ecosystem scale [
5]. The EC measurements represent fluxes in an area around the tower (named footprint) whose size and shape depend on the set-up of the equipment, the structure and height of the canopy and wind direction and speed. This spatial limitation raises the necessity to estimate C fluxes outside the footprint of an EC tower, since it would be costly and unfeasible to install towers to cover all areas of interest.
The net ecosystem exchange (NEE) is the balance between two components, the respiration of the ecosystem (R
eco) and the gross primary production (GPP), which represents the amount of carbon fixed through photosynthesis [
6,
7]. A possible approach to estimate GPP outside the limits of the footprint of an EC tower is the use of remotely sensed information [
8]. Light-use efficiency (LUE) models are the most frequently applied to estimate GPP from spectral vegetation indices (VIs) remotely retrieved [
9].
The LUE concept states that GPP is directly related to the absorbed photosynthetically active radiation (APAR) and photosynthetic efficiency by Photosynthetically Active Radiation (PAR) unit, which can be reduced under environmental stresses such as low temperatures or water stress [
10,
11]. Thanks to the empirical relationship that exists between the fraction of absorbed photosynthetically-active radiation (fAPAR) and VIs [
10,
11,
12], the latter are used in many studies as proxies of GPP. The most widely used VIs to that end are the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI), which are empirically known to be good metrics of green biomass [
8,
9,
13,
14]. NDVI is more sensitive to background reflectance and tends to saturate at high leaf area, which makes EVI a good alternative, as it is less sensitive to these limitations [
14,
15]. Some VIs similar to NDVI and EVI have also been used, such as the Green Normalized Difference Vegetation Index (GNDVI) which is more sensitive to chlorophyll content in plants than NDVI [
16,
17], and the Soil-Adjusted Vegetation Index (SAVI) which is derived from NDVI by including a soil-adjustment factor [
18]. In grasslands, several greenness-related VIs demonstrated their suitability as proxies of GPP. The MERIS Terrestrial Chlorophyll Index (MTCI) was found to provide the best estimates of GPP in subalpine grasslands [
19,
20], EVI and SAVI were the best predictors of GPP in a temperate-alpine grassland [
14], MTCI, PSRI (Plant Senescence Reflectance Index) and GNDVI revealed their suitability in a Mediterranean grassland [
17]. Other studies revealed the usefulness of water-related vegetation indices. In this regard, other VIs used to depict the water content in plants such as Normalized Difference Surface Water Index (NDSWI), Land Surface Water Index (LSWI) and Modified Normalized Difference Water Index (MNDWI) [
13,
21,
22] could also be useful, since they are more sensitive to drought than greenness-related VIs, which makes them likely the best VIs in the estimation of C fluxes in ecosystems undergoing seasonal drought [
23]. In most studies, the selection of VIs [
14,
17,
19] is limited by the availability of remote information from satellite platforms.
An important input of LUE models is the LUE term. While some studies recently attempted to estimate LUE based on some indices such as the photochemical reflectance index (PRI) or the sun induced chlorophyll fluorescence (SIF) [
19,
24], other studies considered the LUE term as a biome-specific constant [
25], one which is usually dependent on seasonal changes [
19,
26]. Besides LUE models, empirical models based on the relationship observed between VIs or the product VIs *PAR and GPP were also successfully frequently applied [
9,
19,
27,
28]. In addition, since photosynthesis (GPP) and respiration (R
eco) are generally positively correlated [
29,
30,
31], reliable estimates of both NEE and GPP have been obtained from remote sensing products.
Grasslands in the Slovenian Karst represent an important resource, as they are one of the first colonizers of shallow soils found in the area. In addition, they are important biodiversity hotspots, are used as extensive pastures by local populations, and contribute to sequestering C in soils [
32]. With climate change, more uncertainty could arise regarding the provision of the aforementioned ecosystem goods and services. Thus, thematic digital maps reporting GPP and NEE estimates in the region might provide short-term information on the ecosystem’s response to climate in an easy, cost-efficient way, with advantages for stakeholders engaged in grasslands productivity and sustainability.
The objective of this study was to estimate NEE and GPP for a karst grassland in the Podgorski Kras plateau (Slovenia) by combining both EC and satellite data in order to provide a basis for large-scale monitoring of the C balance in the grasslands of this region. To that end, we aimed to:
- (i)
Evaluate the ability of different VIs retrieved from remote sensing platforms to represent GPP and NEE trends in a karst grassland;
- (ii)
Compare the performance of different models, integrating VIs in the estimation of GPP and NEE;
- (iii)
Apply obtained results to map NEE and GPP for the grassland area in the Podgorski Kras Plateau comparing the suitability of different remote platforms.
3. Results
3.1. Carbon Fluxes and Environmental Variables Measured by the EC Tower
Fluxes (NEE, GPP, R
eco) and main meteorological variables (VPD, Tair, Rg and P) are reported for the whole study period (2014–2017) in
Figure 2. Large gaps due to power problems are noticeable for flux data. Tair and Rg show no gaps because they were gap-filled using data from a nearby meteorological station in Koper.
3.2. Vegetation Indices
The temporal trends of the considered spectral vegetation indices are reported in
Figure 3. All the vegetation indices that represent vegetation greenness (NDVI, EVI, SAVI and GNDVI) have quite similar trends. In addition, NDVI from the two different sources (Landsat and Proba-V) match quite well, with the difference that NDVI
PV has slightly higher values than NDVI. All the vegetation indices increase from the beginning of the growing season, to reach a peak around end of May or beginning of June. After that, they start decreasing slowly for about 3 months. Around September, there is again a slight increase in the vegetation index. This last increase is better seen with NDVI
PV especially in years 2016 and 2017.
Similar patterns in fluxes were visible every year: GPP shows some seasonality and has two peaks during the growing season, a high peak around end of May–beginning of June and a low peak in October. In between the two peaks, there is a period of low C uptake (lower GPP) values in July and August. Similar trends were visible in NEE and Reco. In fact, a strong correlation was observed between NEE and GPP during the growing season, both for midday and daily averages (r = 0.96 and 0.89 respectively).
The maximum values of VPD and Tair match the low GPP period. Rg, however, reaches its yearly peak earlier than VPD and Tair, somewhere between May and June. Precipitation data shows no real pattern of seasonality, as the distribution over the year seems quite random. However, less precipitation was observed in July and August in years 2016 and 2017.
The year 2016 shows an unusual trend, with a very low value recorded for some vegetation indices in August, due to a fire that occurred at that period. This unexpected disturbance was visible in GNDVI, NDVI, EVI, SAVI, LSWI and MNDWI, but less visible in NDVIPV and NDSVI.
LSWI is a water-related vegetation index, but it has a trend that is similar to the greenness-related vegetation indices. NDSVI, which is also a water-related vegetation index, does not decrease during June, July and August as MNDWI.
3.3. Correlation Charts of Fluxes and Vegetation Indices
A correlation matrix of midday and daily aggregates of fluxes with all vegetation indices is shown in
Figure 4. Correlations were evaluated for the wet phase (
Figure 4, blue line), dry phase (red line) and the entire growing season. The good relationship among most vegetation indices, as previously mentioned, is confirmed by the high correlation coefficients. Only MNDWI and NDSVI were poorly correlated with the other vegetation indices.
On the correlation charts, only separate fit lines are shown (wet and dry phases) in order not to make the graphs overloaded (left bottom of the matrix). However, correlation coefficients were computed also for the whole season (values in black in the matrix). Low correlation coefficients were obtained for the whole season whereas the consideration of the two separate phases in the growing season resulted in higher correlation coefficients.
All the greenness-related VIs (NDVIPV, NDVI, EVI, SAVI and GNDVI) and LSWI gave higher correlation coefficients with GPP and midday NEE in the wet phase compared to the dry phase. In contrast, for daily NEE aggregates, the VIs, except GNDVI, gave a better correlation coefficient during the dry phase. MNDWI gave high correlation coefficients during the dry phase for GPP and NEE, both for midday and daily aggregates. NDSVI gave better correlation coefficients during the wet phase than the dry phase for GPP and NEE with both aggregates, but generally lower than correlation coefficients observed with other VIs.
3.4. Comparison of the Different Models
Values of accuracy metrics for the three models and for midday average fluxes are reported in
Table 3. Similarly,
Table 4 shows accuracy metrics for daily average fluxes. The values of these two tables directly reflect trends and correlation coefficients previously presented in
Figure 4. No exceptions were found in our study, as high values of R
2 corresponded always to low values of RMSE and AIC.
The best model for the entire growing season considering Landsat VIs was always obtained with LSWI for NEE and EVI for GPP in a direct correlation (model 1), either with midday or daily average fluxes. In contrast, NDVIPV gave lower R2 values in all models and for both aggregations of flux data (midday and daily) when the whole growing season was considered.
Overall, separate fits gave higher R2 and lower RMSE for all models than considering the whole growing season. For the wet phase, NDVI was the best predictor of midday aggregates for both GPP and NEE with all three model types. For daily aggregates, NDVI was the best vegetation index only in models 2 and 3 whereas in model 1, LSWI and GNDVI were the best vegetation indices for NEE and GPP respectively. For the dry phase, MNDWI was the best predictor of midday aggregates of GPP and NEE with all three model types. For daily aggregates, MNDWI was the best VI except in model 2 and model 1 with NEE, where LSWI gave the highest R2.
Midday fluxes gave higher R2 than daily fluxes during the wet phase for the model 1 (direct correlation carbon fluxes-VIs). However, during the dry phase (all models) and also the wet phase with models including PAR or Rg (2 and 3), daily fluxes gave higher R2 than midday fluxes.
The bootstrapped R2boot values were very similar to the R2 values, with very narrow 95% confidence intervals. The only time where the best results differ between R2 and R2boot was for midday estimate of NEE, which seems to be the best with MNDWI when considering R2 and LSWI if R2boot were considered instead.
Based on these accuracy metrics, the best regression models are presented in
Table 5, according to the two aggregates (midday and daily), sources of vegetation indices (Landsat and Proba-V) and fluxes (GPP and NEE) for the two phases of the growing season. Direct correlation (model type 1) was found to be the best option for midday average fluxes, except for the estimation of GPP with NDVI
PV during the wet phase, where the model type 2 was the best. For daily averages, however, the model type 1 performed better with NEE whereas model types 2 and 3 were the best for GPP.
During the wet phase, the best models obtained for GPP gave higher coefficients of determination (R2) compared to those obtained with NEE. The opposite was observed for the dry phase.
Based on the best models reported in
Table 5, the
Table 6 gives midday and daily estimates of GPP and NEE from NDVI
PV data for the growing season (March to October) of the years 2014 to 2017. Average flux estimates of the four years are −711 g C m
−2 growing season
−1 and −199 g C m
−2 growing season
−1 for GPP and NEE respectively.
3.5. Flux Maps Using the Best Models
Figure 5 shows a map of GPP and NEE for the whole study area, computed based on the equations in
Table 5. These fluxes represent estimated averages for the periods 15 May 2017 to 19 May 2017 for Landsat and 11 May 2017 to 20 May 2017 for Proba-V. Similarly,
Figure 6 shows estimated average of NEE and GPP for the periods 2 July 2017 to 6 July 2017 for Landsat and 1 July 2017 to 10 July 2017 for Proba-V. The dates were chosen to show an example for both wet and dry phases.
During the wet phase (
Figure 5), GPP or NEE estimated from the two different sources of vegetation indices have similar spatial distributions in spite of their different spatial resolution. In the case of midday average fluxes, NEE and GPP seem to have the same distribution with differences only in their values. This is a consequence of the very good linear relationship that was found between midday GPP and NEE in our study. Daily fluxes, however, show more difference between NEE and GPP, as the correlation observed was not as good as that of midday fluxes.
During the dry phase (
Figure 6), the spatial distribution of fluxes is quite different between flux estimates from the two sources of VI (Landsat and Proba-V), with the exception of estimates of daily aggregates of NEE, where the estimate from NDVI
PV appears like a generalization of the estimate from LSWI.
A difference can be noticed between the wet and the dry phases: GPP estimates are very low during the dry phase, and the carbon balance (NEE estimates) became positive in some parts of the study area.
4. Discussion
In our study, the low C uptake observed during the period June to August with high VPD and Tair, and relatively lower precipitation confirmed that summer water stress is a major constraint in C sequestration for grasslands in the karst region characterized by shallow soils [
36]. Regression models obtained demonstrated the ability of spectral vegetation indices to represent the impact of climate constraints on both NEE and GPP.
4.1. Performance of the Different Vegetation Indices
Overall, NDVI was the best predictor for the estimation of NEE and GPP during the wet phase of the growing season, reinforcing the conclusion of previous studies that NDVI is a suitable index for depicting photosynthetic C assimilation [
9,
64,
65]. In addition, the best results obtained with NDVI instead of EVI or GNDVI suggest that the frequently reported issue of NDVI saturation was not observed in our study. Indeed, NDVI values observed are below the threshold of 0.8 reported for grasslands above which saturation has been observed [
66,
67]. On the other hand, the fact that greenness related vegetation indices were less performant when the entire growing season was considered is probably due to the influence of increased VPD on GPP and tissue water content in the karst grassland addressed by our study, during the dry phase. In fact, LSWI and MNDWI proved to be ideal for estimating fluxes during the dry period, confirming the fact that water-related vegetation indices are the best indicators of photosynthetic activity during dry periods, as they are more sensitive to water availability than the other vegetation indices [
23]. NDSVI, which is also a water-related vegetation index, did not give as good results for the dry phase as LSWI and MNDWI. Other studies reported the low performance of NDSVI in depicting non-photosynthetic vegetation in a semi-arid grassland and a cropland [
68,
69], which is confirmed by the low performance observed in our study.
The two sources of vegetation indices (Landsat and Proba-V) could not really be compared, given that Landsat images were single date images whereas NDVIPV are 10 days composite data. In addition, the difference in the time step of aggregation of fluxes (10 days for NDVIPV and 5 days for Landsat vegetation indices) makes a comparison inappropriate. However, it should be noted that the highest correlation coefficient obtained with NDVIPV is generally lower than the one obtained with Landsat vegetation indices. This is probably due to the greater variability in fluxes and NDVIPV within the aggregation timeframe of the NDVIPV data (10 days).
The similarity in the correlation observed for most of the VIs with GPP and NEE is probably the consequence of the good relationship observed between the two variables, which often occurs when the analysis is focused on the growing season [
29,
30,
31].
Despite the fire event that occurred in August 2016, we did not notice any outliers corresponding to that period in our regressions. This might be explained by the occurrence of the fire event towards the end of the dry phase, and the grasses recovered during the second part of the wet phase in September. Indeed, grasslands are known to recover quickly following a fire event [
70].
4.2. Performance of the Different Models
In agreement with recent evidence [
9,
17], we observed an increase of model fitness when data were split into a dry and a wet phase, in contrast with the general agreement that LUE can be considered constant in grasslands along the whole growing season. Recent studies addressed this problem by attempting independent LUE estimations from the PRI and SIF [
19,
24]. Other models constrain LUE by decreasing a maximum value (LUE
max) with meteorological variables (scalars) [
11]. In the present study, good results have been achieved by splitting the growing season into a wet and dry phases, reducing model input and further calculations.
The higher correlation coefficients observed between VIs and midday average fluxes compared to daily average fluxes during the wet phase could be due to the fact that most C uptake occur around midday during that phase of the growing season. During the dry phase however, depending on daily environmental conditions, there can be more variation in the midday average fluxes compared to daily averages. Hence during dry periods, daily averages showed more robust and higher correlations with VIs than midday averages. Overall, despite the differences observed, daily estimates were robust enough in our study to suggest their use instead of midday average flux estimates, being more representative of the total productivity of the ecosystem.
The validation of the models using the bootstrapping resampling method confirmed the assessment made with the regular accuracy metrics. The narrow 95% confidence intervals obtained for the R2boot reflects the low variability of the coefficient of determination obtained from the different bootstrapping samples.
From our carbon flux estimates, we could deduce that during the growing season, midday GPP represents over half of the whole day fluxes. NEE values are negative, which means that at least during the growing season, the grassland of our study area acts as a carbon sink. To our knowledge, there are no studies that estimated carbon fluxes in a karst grassland. However, the average GPP and NEE of −711 and −199 g C m
−2 per growing season respectively are in range with estimates, of −867 and −132 g C m
−2 per season respectively for GPP and NEE, obtained by [
71] in a Mediterranean grassland.
4.3. GPP and NEE Maps
The use of regression models in our study for mapping GPP or NEE helped to appreciate the spatial variations in C uptake and C sequestration in a homogeneous area and estimate the contribution of grasslands. However, the spatial variability in C fluxes depended on the spatial variability of the vegetation indices. The obtained models should therefore be used only in homogenous environmental conditions.
Carbon fluxes maps obtained from NDVIPV appeared as a generalization of that obtained from Landsat derived VIs, which reflects the different spatial resolution of data from the two platforms. This suggests the possibility to rely on NDVIPV data for its availability, especially when a continuous estimation of carbon fluxes is needed over a certain period of time.
During drought periods, an ecosystem can quickly shift from carbon sink to source [
72]. The difference observed in the maps of NEE between the wet and the dry phases illustrates this fact, since NEE values became positive during the dry phase, demonstrating the importance of comparing maps of the same area over time.