Estimation of carbon fluxes from eddy covariance data and satellite-derived vegetation indices in a karst grassland (Podgorski Kras, Slovenia)

.................................................................................................................................................. 3 Resumo alargado ..................................................................................................................................... 5 List of Figures ......................................................................................................................................... 8 List of Abbreviations ............................................................................................................................... 9


Acknowledgements
First of all, I would like to thank God for helping me throughout this important step in my life.
I would like to express my sincere gratitude to Sofia Cerasoli and Mitja Ferlan for their very good supervision and support all the way. I thank them for contributing greatly to reach this quality of the dissertation.
I am also very thankful to the president of jury Prof. Pedro Ochôa and the jury member João Silva for their assessment of the document, and their valuable comments.
I am really grateful to the MEDfOR commission and the Education, Audiovisual and Culture Executive Agency (EACEA) of the European Union for this opportunity to study in some top universities in Europe.
Finally, I would like to thank my family and friends for their unfailing support so that whenever I felt down, I had someone to turn to.

Introduction
Grasslands are one of the most widespread vegetation types worldwide, covering between 14 and 26% of the earth surface (Mason and Zanner, 2005;Scurlock and Hall, 1998). Moreover, they are increasing in area recently due to the abandonment of former agriculture lands in some parts of the world. In fact, the transformation of the agricultural context has released many areas from use in Europe, and such areas would naturally undergo succession (Benjamin et al., 2005) from grasslands to shrublands and forests.
In Slovenia, abandoned agriculture lands in the karst region are known to go through some successional stages such as grasslands and woody ecosystems (Ferlan et al., 2016;Kaligarič et al., 2006). Each succession stage has a different carbon balance. In fact, while forests are known to act as a sink accumulating carbon in their woody biomass (Ferlan et al., 2016;Pan et al., 2011;Post and Kwon, 2000), grasslands may act as a source of carbon (Ferlan et al., 2011) particularly in periods of drought when the lack of precipitation can decrease photosynthetic carbon uptake (Arnone III et al., 2008;Gilmanov et al., 2007;Meyers, 2001), or as a consequence of other disturbances (e.g. fire events). This led to a poor recognition of the role of grassland ecosystems in the global carbon cycle Hall and Scurlock, 1991). However, grasslands play a tremendous role as they stock an important amount of carbon in their soils, estimated at 23% of the global soil carbon (Buringh, 1984). Given this fact, it is important to study the extent of carbon sequestration in grasslands and understand how they contribute in mitigating the effects of climate change.
The Eddy Covariance (EC) method permits to measure, at ecosystem scale, the atmosphere-ecosystem exchange of water, energy and CO2 fluxes (Papale et al., 2006). The method consists of measurements of the net exchange of gases between the atmosphere and the ecosystem above the canopy where turbulence can be considered more or less constant (Ferlan, 2013). It has been used for the first time in the 1970s (Baldocchi et al., 1988;Desjardins, 1974) and since then has been widely employed in different ecosystem types around the world (Baldocchi, 2008;Ferlan et al., 2011Ferlan et al., , 2016Haszpra et al., 2005; Y.-L. Li et al., 2008;Peichl et al., 2012;Ruimy et al., 1995;Saigusa et al., 2002;Yan et al., 2015;Yao et al., 2018). The eddy covariance method provides a reliable direct measurement of different gas compounds together with meteorological variables, at high temporal detail, permitting to ascertain the influence of climate drivers or other disturbances on ecosystem fluxes (Burba and Anderson, 2010).
The eddy covariance measurements represent fluxes in an area around the tower (named footprint), the size and shape of which depend on the set-up of the equipment, the structure and height of the canopy and varies with prevalent wind direction and speed. Usually, the footprint extends over a distance ranging from tens of meters to more than 1km from the tower (Göckede et al., 2008). In the Slovenian karst grassland, footprint analyses showed that the mean distance from the tower is about 195 m (Ferlan, 2013). This spatial limitation raises the necessity to find a way to estimate carbon fluxes outside the footprint of the eddy covariance tower, since it would be costly and unfeasible to install towers to cover all areas of interest.
NEE is the result of the balance between two components, the respiration of the ecosystem (Reco) and the Gross Primary Productivity (GPP), which represents the amount of carbon fixed through photosynthesis. A possible approach to estimate GPP outside the limits of the footprint of an eddy covariance tower is the use of remotely sensed information. In fact, optical data such as satellite images can be applied to estimate carbon fluxes by establishing relationships between remotely sensed information and fluxes in order to extrapolate GPP outside the boundaries of a tower footprint.
Reflectance values obtained from remote platforms in specific wavelengths are generally employed to calculate normalized differences to obtain spectral vegetation indices (VIs).
Light-Use Efficiency models are the most frequently applied to estimate GPP from vegetation indices (VIs) remotely retrieved (Nestola et al., 2016). Besides LUE models, empirical models based on the relationship observed between VIs and GPP are also frequently applied (Gilmanov et al., 2005;Zhengquan Li et al., 2007;Nestola et al., 2016;Rossini et al., 2012). In addition, since photosynthesis (GPP) and respiration (Reco) are generally positively correlated (Baldocchi, 2008;Baldocchi et al., 2015;Ma et al., 2016), it is often possible to infer both NEE and GPP from remote sensing products.
The general objective of this study is to estimate NEE and GPP for a karst grassland in the Podgorski Kras plateau by combining both eddy covariance and satellite data in order to provide a basis for largescale monitoring of the carbon balance in the Podgorski karst grassland.
In order to reach that objective, this study aims 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 a grassland area in the Podgorski Kras Plateau.

Principle of the eddy covariance method
Airflow consists of numerous eddies. The general principle of eddy covariance measurements can be understood as the covariance between the concentration of interest and vertical wind speed in the eddies (Burba and Anderson, 2010). Put simpler, it consists of measuring how many particles of a component of interest are moving up and down over time and how fast they are (Ferlan, 2013). The horizontal airflow over an investigated area is composed of numerous rotating eddies that can be represented at a single point on the tower by the Figure 1. At a given moment (time 1), eddy 1 moves air parcel c1 downward with the speed w1. At the same point, the next moment (time 2), eddy 2 moves air parcel c2 upward with speed w2. Given that each air parcel has its own characteristics i.e. gas concentration, temperature, humidity, etc., if these characteristics and the speed of vertical air movement are measured, the vertical upward or downward fluxes can be known (Burba and Anderson, 2010). The eddy covariance method has the advantage of performing direct measurements for different types of gas (CO2, H2O, CH4, etc.), with high precision and detail (Burba and Anderson, 2010). However, there are some problems related to the eddy covariance method. The main problem is the occurrence of gaps in the data, which need to be filled by statistical regressions with different methods. Gaps occur due to power breaks (mostly when power system is based on solar panels), damages to instruments, for instance, due to animals or lightning (Aubinet et al., 2012). In addition, equipment malfunctioning such as the anemometers that might not work during heavy precipitation events would make the eddy covariance system glitch during rainy periods (Ferlan, 2013). Another limitation of the method is its restriction to flat areas (Burba and Anderson, 2010).
The Net ecosystem Exchange (NEE) of CO2 between the atmosphere and the biosphere measured by eddy covariance can be partitioned into the two components of carbon fluxes, Gross Primary Productivity (GPP) and ecosystem respiration (Reco) (Lasslop, Reichstein, Papale, et al., 2010;Reichstein et al., 2005). GPP refers to the total amount of carbon fixed in the process of photosynthesis by plants in an ecosystem; Reco is the amount of carbon lost by autotrophic and heterotrophic respiration. NEE refers to the balance between GPP and carbon losses due to ecosystem respiration (Aubinet et al., 2012;Kirschbaum et al., 2001) as in the equation (1) hereafter.

NEE = GPP + Reco (1)
Where NEE is the Net Ecosystem Exchange, GPP is the Gross Primary Productivity and Reco is the Ecosystem Respiration. We adopted in this study the atmospheric sign convention where a flux toward the surface (carbon uptake, i.e. GPP) is negative whereas a flux upward the atmosphere (carbon release, i.e. Reco) is positive (Baldocchi, 2008;Lasslop, Reichstein, Detto, et al., 2010). Consequently, a negative NEE indicates that the ecosystem is acting as a carbon sink while a positive NEE indicates that the ecosystem is acting as a source of carbon.

Carbon fluxes estimates integrating remote sensing data
With the increasing availability of spatial data thanks to technological progress (Lillesand et al., 2004), their integration into GPP models raised considerably (Mäkelä et al., 2007). LUE models are the most widely used models for integrating carbon fluxes and optical measurements (Lees et al., 2018;Nestola et al., 2016), thanks to the ease of applicability, with the possibility to retrieve indirectly all variables (Bagnara et al., 2018).
LUE models can be expressed in a general form as follows (Yuan et al., 2014): Where PAR is the incident photosynthetically active radiation (MJ m -2 ); fAPAR is the fraction of absorbed photosynthetically active radiation; LUEmax is the potential LUE (g C m -2 MJ -1 APAR) under ideal environmental conditions; and M is a modifier that depends on environmental conditions, and constrains LUEmax to its actual value.
One main limitation of LUE models is related to the LUEmax term in the models. It is expressed as a biome-specific constant in most of the models (Goerner et al., 2011;Rossini et al., 2012).
In grasslands, APAR can be assumed to explain most of the variability in GPP (Lobell et al., 2003) and therefore LUE (LUEmax*M in equation 1) is generally considered constant (Nestola et al., 2016).
However, considering LUE as a constant leads often to errors in the estimate of GPP. In their study, Nestola et al. (2016) obtained better results splitting the analysis in two parts (greening and senescence phases) along the growing season. By doing so, and considering a different LUE for each of the 2 periods, carbon fluxes were more accurately estimated.
Ground cover and leaf area are significant variables that determine absorption of PAR by the canopy.
Thanks to the empirical relation that exists between fAPAR and vegetation indices (Running et al., 2004;Yuan et al., 2014), the latter are used in many studies as proxy of carbon fluxes (Nestola et al., 2016;Yan et al., 2015;Y. Zhou et al., 2014). In fact, given the robust relationship between fAPAR and Leaf Area Index (X. Zhou et al., 2002), the fAPAR can be determined based on vegetation indices derived from remote observations of surface spectral reflectance (Myneni and Williams, 1994).
Despite the increasingly available panoply of vegetation indices, their use for estimating fAPAR has been limited to one or two (Y. Zhou et al., 2014). Some of the most widely used indices include NDVI and EVI (Yan et al., 2015;Y. Zhou et al., 2014). In fact, NDVI and fAPAR increase with ground cover and plant leaf area, and their good relationship made it possible to estimate fAPAR from NDVI in many studies (Myneni and Williams, 1994). Nestola et al. (2016) confirmed the effectiveness of NDVI as a metrics of green biomass, making it a useful parameter in a simple expression of the LUE model for a grassland ecosystem.
Despite NDVI being a good metrics of green vegetation, it is quite sensitive to background reflectance and tends to saturate at high leaf area. In such conditions, EVI could be employed as an alternative to NDVI because it is less sensitive to these limitations (Rocha and Shaver, 2009).
Other vegetation indices similar to NDVI and EVI could also be explored. For instance, by using the Green instead of the Red band for the calculation of NDVI, the Green NDVI (GNDVI) is more sensitive to chlorophyll content in plants than NDVI (Gitelson et al., 1996). The Soil Adjusted Vegetation index (SAVI) is modified from NDVI by including a soil-adjustment factor (Jovanović et al., 2016).
To overcome the main limitation of LUE models (difficulties in the estimation of the LUE term of the models), some authors attempted to estimate NEE or GPP through linear regression, considering vegetation indices as independent variables (Nestola et al., 2016;Rossini et al., 2012).
In these empirical models, other satellite-derived indices not only related with fAPAR can be explored, for example vegetation indices known to be able to depict surface water content. These indices include the Normalized Difference Senescent Vegetation Index (NDSVI), the Land Surface Water Index (LSWI) and the Modified Normalized Difference Water Index (MNDWI) (Hill, 2013;John et al., 2008;Yan et al., 2015). Water-related vegetation indices can be very good indicators of plant activity in summer when other vegetation indices could be stationary due to some greenness of the plants despite the very low photosynthetic activity. In fact, they are more sensitive to drought than greenness related vegetation indices (Bajgain et al., 2015).
While vegetation indices, and NDVI in particular, proved in many cases to be effective in approximating GPP through the fAPAR component in a LUE model (Myneni and Williams, 1994;Nestola et al., 2016), they are less correlated with ecosystem respiration making Reco usually the most important source of uncertainty in NEE estimation through remote sensing (Yan et al., 2015). Moreover, ecosystem respiration is the sum of heterotrophic (microbes, soil fauna) and autotrophic (plant roots) respiration (Bond-Lamberty et al., 2004;Hanson et al., 2000), which would make its estimation more difficult in some complex ecosystems involving important contribution of heterotrophic respiration. However in some cases it was possible to estimate NEE adopting models integrating remote sensing products (e.g. Nestola et al. 2016).

Study area
The present study was conducted in the Podgorski Kras plateau located in the sub-mediterranean region of south-west Slovenia. The area underwent major human influences, due to its position at the transition between the Mediterranean and central Europe. In fact, agricultural practices such as overgrazing in the past centuries led to a pronounced destruction of the vegetation cover, causing severe soil erosion and resulting into a stony and bare landscape. However, thanks to the economic development causing the near-abandonment of agriculture practices in the area, a succession is taking place and we can observe in the plateau different vegetation types ranging from grasslands to secondary oak forests (Ferlan, 2013).
The bedrock is mainly composed of limestone from Paleocene and Eocene (Knez et al., 2015). The chemical weathering known as karst phenomena led to the formation of Leptosols and Cambisols, which represent insoluble fractions of carbonates. As a result, the soil is superficial, with depths ranging from 0 cm to several decimeters in soil pockets between rocks. The organic matter represents about 12-15% of the topsoil (Ferlan, 2013).
The climate of the area is transient between mediterranean and continental. It is more humid than a typical Mediterranean climate, with less pronounced dry period in summer and colder winter. It is often referred to as sub-mediterranean climate, with a mean annual temperature of 10.5°C, a mean daily  Air temperature and global radiation data were gapfilled based on data from a meteorological station located in Koper (at a distance of 15 Km from the tower). NEE data were partitioned into GPP and Reco according to  using daytime data-based estimates, considering temperature sensitivity of respiration and VPD limitation of GPP.

Spectral vegetation indices
In this study, two sources of data were considered. The 300m resolution NDVI data, consisting of 10 Subtraction in order to remove atmospheric components such as scattering and absorption of solar energy in the atmosphere and obtain Top of Canopy reflectance (Chavez, 1996). The image processing was conducted with the ENVI 5.1 software.  (Table 2), an average value was computed for the pixels in the footprint (Figure 2) using the Raster package in R software.  (2006) L is a constant dependent on the vegetation cover and takes values from 0 (for very green vegetation) to 1 (areas with no green vegetation), assumed 0.5 here.
Despite the availability of flux data since 2008, only a timeframe of four years was considered in order to match the timeframe of the remote sensing information used in this study.

Data analysis
In this study, we considered two types of aggregation of flux data. The first type of aggregation consisted

GPP = (a*VI+b)*PAR (5)
All models were tested for the entire growing season (single) or splitting the growing season in two phases (green and dry). The separation of the growing season was based on preliminary tests during which we plotted GPP or NEE as a function of VIs and tried a separation based on months. The months of June, July and August allowed a visual identification of a different group. However, the separation was not perfect and different values of midday aggregates of Tair and VPD (averaged over 10 days for SPOT-Vegetation NDVIs and 5 days for Landsat VIs) were tried instead of months. VPD proved to be the best for separating the growing season into 2 phases, with a threshold of 1500 Pa. We defined the greening phase as the period of the growing season with midday average VPD less than or equal to 1500 Pa, and the dry phase as the period of the growing season with midday average VPD greater than 1500 Pa.
In order to compare the performance of the models obtained from the different regressions, three accuracy metrics were computed, namely the coefficient of determination (R²), the Root Mean Square Error (RMSE) and the Akaike Information Criterion (AIC). The best models are the ones with a high value of R² and a low value of RMSE and AIC.
All analyses were done using the R software, version 3.4.4.
The best models selected were used to create illustrative GPP and NEE maps of the study area for two

Carbon fluxes and environmental variables
The Figure 3 below presents, for the period 2014-2017, fluxes (NEE, GPP, Reco) and some main environmental variables including vapor pressure deficit (VPD), air temperature (Tair), global radiation (Rg) and precipitation (P). Large gaps are noticeable for flux data and VPD, whereas Tair and Rg show no gaps because they were gap-filled from another meteorological station.
GPP shows some seasonality and has two peaks during the growing season, a high peak around end of May or beginning of June and a low peak in October. In between the two peaks, there is a period of low carbon uptake translating into low 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, both for midday and daily averages ( Figure 5).
The maximum values of VPD and Tair match the low GPP period. Global radiation 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, there is generally less precipitation in July and August than in any other months, in some years.

Vegetation indices
The Figure 4 represents all the spectral vegetation indices used in this study. 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 SPOT-Vegetation) match quite well, with the difference that NDVIs has slightly higher values than NDVI. All these vegetation indices increase from the beginning of the growing season, to reach a peak around end of May or beginning of June.
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 NDVIs.  Table 2). Different symbols represent different years.
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 NDVIs and NDSVI.
LSWI is a water-related vegetation index, but it has a trend that is similar to the greenness-related vegetation indices. This can be explained by the fact that vegetation greenness reflects the water content in plant leaves, which is detected through LSWI. MNDWI on the contrary has an overall different trend since it is really meant to detect water bodies (Xu, 2006), but it also shows a decrease during the months of June, July and August. NDSVI, which is also a water-related vegetation index, does not decrease during June, July and August as MNDWI.
These observations about the vegetation indices along with the trend of fluxes and environmental variables suggest that the growing season can be subdivided in two phases, a greening phase (March to May/June and September to October) and a dry phase (June to August) when the higher Tair and VPD, together with a decrease in precipitation lead to a decrease in GPP. However, these two phases could not be clearly defined by months, since climate is quite variable from one year to another. Our separation, based on VPD instead of months, proved to better split the two phases of the growing season with the threshold of 1500 Pa.

Correlation charts of fluxes and vegetation indices
The Figure  On the correlation charts, only separate fit lines are showed in order not to make the graphs overloaded (left bottom of the matrix). However, correlation coefficients were computed also for the single fit (values in black in the matrix). Low correlation coefficients were obtained for the single fit whereas the consideration of the two phases of the growing season gave higher correlation coefficients.
All the greenness-related VIs (NDVIs, NDVI, EVI, SAVI and GNDVI) and LSWI gave higher correlation coefficients in the greening phase compared to the dry phase for GPP (both midday and daily aggregates) and midday NEE. For daily NEE aggregates, those VIs except GNDVI gave a better correlation coefficient during the dry phase instead. 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 greening phase than the dry phase for GPP and NEE with both aggregates, but generally lower than correlation coefficients observed with other VIs. Values in black, blue and red represent Pearson correlation coefficients for a single fit, the greening phase (VPD < 1500 Pa) and the dry phase (VPD > 1500 Pa) respectively. The remark "_24" is used to distinguish daily average fluxes from midday average fluxes. *, **, *** indicates correlations significant at 95%, 99% and 99.9% Midday average fluxes generally show a better correlation with VIs than daily averages during the greening phase, while the opposite was observed during the dry phase. This observation was the same both for GPP and NEE. The best correlation coefficients observed for midday averages during the greening phase could be due to less fluctuation in the midday average fluxes due to the reduced time scale considered and the fact that most carbon uptake occur around midday during that phase of the growing season. During the dry phase however, there is less carbon uptake during the day, and midday averages may fluctuate more depending on particular environmental conditions of every single day (higher temperatures of some days may lead to more respiration). Daily averages would hide those fluctuations, explaining why we obtained better correlations with VIs than midday averages during the dry phase.

Comparison of the different models
The Table 3 shows values of accuracy metrics for the three types of model considered in this study, for midday average fluxes. Similarly, Table 4 shows accuracy metrics considering daily average fluxes. The values of these tables directly reflect trends and correlation coefficients previously presented in Figure   5. No exceptions were found in our study, as high values of R² corresponded always to low values of RMSE and AIC. The three metrics are all however useful, since slight differences can be noticed in the RMSE and AIC when R² is the same, and this helped in the choice of the best model.
The VIs performed differently from one model to another, according to the phase of the growing season being considered and the flux aggregate being estimated.
What was mentioned in the previous section, about a poor correlation between fluxes and vegetation indices with a single fit, can be confirmed by generally low values of R² or higher values of RMSE and AIC. The best model with a single fit considering Landsat VIs was always obtained with LSWI for NEE and EVI for GPP in a direct correlation (model 1), be it with midday average or daily average fluxes.
NDVIs however gave lower R² values in all models and for both aggregations of flux data (midday and daily) when a single fit for the whole season was considered, and the highest R² observed was 0.43.
Overall, separate fits gave higher R² and lower RMSE and AIC in all models compared to a single fit.
In fact, the highest R² obtained in a single fit was 0.59 for midday fluxes and 0.62 for daily fluxes, but it was possible to obtain 0.85 and 0.91 during the greening phase for midday and daily aggregates respectively, 0.69 and 0.86 during the dry phase for midday and daily aggregates respectively.
For the greening phase, NDVI was the best predictor of midday aggregates of 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 R².
Midday fluxes gave higher R² than daily fluxes during the greening phase for the model type 1 (direct correlation). However, during the dry phase (all model types) and also the greening phase with models including PAR (2 and 3), daily fluxes gave higher R² than midday fluxes.  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 SPOT-Vegetation) 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 NDVIs during the greening 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 greening phase, the best models obtained for GPP gave higher coefficients of determination (R²) compared to those obtained with NEE. The opposite was observed for the dry phase.

Flux maps using the best models
The Figure 6 shows a map of fluxes for the study area, computed based on the equations of Table 5  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 7), the spatial distribution of fluxes is quite different between flux estimates from the two sources of VI (Landsat and SPOT-Vegetation), with the exception of estimates of daily aggregates of NEE, where the estimate from NDVIs appears like a generalization of the estimate from LSWI.

Figure 7:
Maps of average midday and daily GPP and NEE estimates for the periods 02/07/2017 to 06/07/2017 (for Landsat) and 01/07/2017 to 10/07/2017 (for SPOT), using the best models obtained for the dry phase. The remark "_24" is used to distinguish daily from midday average fluxes.
A difference can be noticed between the greening 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.

Discussion
Several climatic factors, including temperature, precipitation, solar radiation and water deficit interact to influence the variability observed in carbon fluxes (Griffis et al., 2000), making it difficult to separate their individual effects clearly (Hui et al., 2003). In our study, the low carbon uptake observed during a period of the year (June to August) with high VPD and Tair, and relatively lower precipitation confirmed that climate is a key driver of carbon fluxes. The impact of climatic drivers could be stronger in this karst grassland, characterized by shallow soils. In fact, in ecosystems facing drought periods (due to low precipitation rates and shallow soils), GPP and consequently NEE are primarily controlled by climatic factors such as precipitation (Ferlan et al., 2011). In addition, the fact that grasses use water intensively (Rodriguez-Iturbe et al., 2001) makes precipitation an important variable that affects the variability in carbon fluxes.
The separation of the growing season into two phases as in Nestola et al. (2016) proved to be useful for a better correlation between fluxes and vegetation indices. In our study, the threshold of 1500 Pa for VPD showed to be a threshold across years in spite of differences observed in precipitation.
The good relationship between GPP and NEE (Baldocchi, 2008;Baldocchi et al., 2015;Ma et al., 2016) was confirmed in our study, allowing to infer both NEE and GPP from vegetation indices.
Overall, NDVI was the best predictor for the estimation of NEE and GPP during the greening phase of the growing season, reinforcing the conclusion of previous studies that NDVI is an index of choice for depicting photosynthetic activity (Gamon et al., 1995;Myneni et al., 1995;Nestola et al., 2016) which translates into CO2 assimilation (GPP component of carbon fluxes). On the other hand, LSWI and MNDWI proved to be ideal for estimating fluxes during the dry period, confirming the fact that waterrelated vegetation indices are the best indicators of photosynthetic activity during the dry season, as they are more sensitive to water stress than the other vegetation indices (Bajgain et al., 2015). NDSVI, which is also a water-related vegetation index, did not give good results for the dry phase as for LSWI and MNDWI. Other studies reported the low performance of NDSVI in depicting non-photosynthetic vegetation in a semi-arid grassland and a cropland (Zhaoqin Li and Guo, 2018;Sonmez and Slater, 2016), which relates to the low performance observed in our study.
Strong direct correlations were observed between fluxes and vegetation indices (model type 1) as in Rossini et al. (2012). This confirms that APAR (which is highly dependent on the vegetation greenness) explains most of the variability of GPP in an ecosystem characterized by a seasonality in greening and senescence such as grasslands and croplands (Gamon, 2015;Gitelson et al., 2006;Lobell et al., 2003).
However, the fact that GNDVI did not perform well generally suggests that chlorophyll content is not the only factor that comes into play in this particular karst grassland, which seemed to have an important water stress influence during the dry phase.
Despite the fire event that occurred in August 2016, which was easy to depict in the VIs profiles, 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 greening phase in September. Indeed, grasslands are known to recover quickly following a fire event (Keeley and Keeley, 1984).
An interesting point about this study is that reasonable results were obtained also for daily fluxes. We consider daily average flux estimates more useful than midday average flux estimates in terms of monitoring the total productivity of an ecosystem. Even though a daily (24 hours) average flux was considered instead of the sum due to missing data, we considered that the average without the missing data represents a good proxy of the total daily flux.
The two sources of vegetation indices (Landsat and SPOT) could not really be compared, given that Landsat images were single date images whereas NDVIs is a 10 days composite data. In addition, the difference in the time step of aggregation of fluxes (10 days for NDVIs and 5 days for Landsat vegetation indices) rendered a comparison inappropriate. However, it should be noted that the highest correlation coefficient obtained with NDVIs is generally lower than the one obtained with Landsat vegetation indices. This is probably due to the greater variability in fluxes and NDVIs within the timeframe considered for the NDVIs data (10 days).
The use of the regression models in our study for mapping GPP or NEE seems interesting as it helped to appreciate the spatial distribution of carbon fluxes. However, the spatial variability in carbon fluxes only depended on the spatial variability of the vegetation indices. The obtained models should therefore be used only in homogenous environmental conditions.
During drought periods, an ecosystem can quickly shift from carbon sink to source (Lei et al., 2016).
The difference observed in the maps of NEE between the greening 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.
In the approach adopted in this study, VPD values from the eddy covariance tower are still needed to identify the greening and dry phases for a given date. This need of VPD in order to choose which model to apply could hinder the applicability over large areas where VPD values may change significantly.
The possibility to retrieve VPD from meteorological stations could be explored in that case.
Another limitation of our study is the fact that there was no validation of the regression models. This is due the limited number of Landsat 8 images available along with missing flux data, not permitting to spare some data for validation. However, we consider that the model selection based on the good correlations observed and the regression accuracy metrics was robust enough and indicative of the utility of empirical regressions.

Conclusion
This study was particularly interesting, as it investigated the ability of several vegetation indices from two different remote sources in estimating carbon fluxes. It was equally compelling because it was conducted in a karst grassland, which is a particular ecosystem due to the rocky nature of the bedrock, and its relatively shallow soils.
The eddy covariance method represents a suitable method for measuring gas fluxes, which is widely used to assess and monitor the atmosphere-ecosystem carbon exchange. The fact that its measurements are limited to a relatively small footprint was addressed in this study, by attempting to estimate GPP and NEE for a homogenous area larger than the tower footprint. GPP showed some patterns of seasonality with two peaks, one at the end of May and another in October. In between the two peaks, a period of low GPP was identified, which we later called dry phase while the rest of the growing season was referred to as greening phase. VPD was used to distinguish between the two phases, with the threshold of 1500 Pa above which begins the dry phase. Spectral vegetation indices from two sources (Landsat and SPOT-Vegetation) were explored, to test their ability as predictors in the estimation of NEE and GPP.
The main results of our study are the followings:  A strong relationship was found between most of the vegetation indices and also between NEE and GPP, explaining the comparable performance of some of the regression models.
 NDVI proved to be the best predictor of GPP and NEE during the greening phase as the carbon uptake during that phase is mostly explained by changes in the green biomass in the plants. For the dry phase however, water-related vegetation indices such as LSWI and MNDWI were the best predictors of GPP and NEE, since they are more sensitive to drought.
 Good results were obtained not only for midday but also for daily (24 hours) fluxes, which is interesting since they could be more useful in terms of monitoring the carbon balance of an ecosystem, instead of considering only midday fluxes.
 Model type 1 (direct correlation between fluxes and vegetation indices) was the best for estimating midday average fluxes (GPP and NEE) and daily average of NEE whereas the models types 2 (direct correlation between fluxes and the product of PAR and vegetation indices) and 3 (simplified LUE model) were the best for daily average GPP.
 The maps of GPP and NEE could help appreciate spatial and temporal variability of fluxes. The two sources of NDVI data gave similar maps, in most of the cases. In those cases, the maps derived from NDVIs appeared as a generalization of maps derived from Landsat NDVI due to a lower spatial resolution of SPOT-Vegetation data. Following the spatial detail, the range of GPP and NEE values in the maps derived from Landsat is wider than the range of GPP and NEE derived from SPOT-Vegetation.
The results of this study seemed quite interesting, and the best regression models could be used for estimating carbon fluxes outside of the eddy covariance tower. Suggestions for further studies would be the validation of the models once more data would be available. The availability of Sentinel 2 images could also be considered. Alternatively, Landsat 7 images could be explored in order to use all the available eddy covariance data prior to 2014.