## 2 Materials and Methods
### Data Analysis and Modelling  
The general approach consisted of flux and environmental parameters aggregation, followed by regressions between aggregated fluxes and VIs. Two types of aggregation of flux data were considered. For the first data aggregation type, midday (between 11 a.m. and 4 p.m.) averages of half-hourly fluxes were calculated as in [9]. For the second type of aggregation, we calculated daily averages of half-hourly fluxes. For further analysis, midday and daily aggregated flux data were subsequently averaged over different time steps as follows: for NDVI\({}_{\text{PV}}\) data, aggregated flux data were averaged over 10 days period to match temporal aggregation provided by Proba-V. For VIs derived from Landsat 8 images, midday and daily flux data were averaged over 5 days (acquisition date plus 4 days prior to this date), since a preliminary test showed a better correlation if aggregated fluxes were considered instead of the overpass day only. Three types of model linking VIs and EC data were tested in this study after [19]:  
(i) Model 1 assuming a direct linear relationship between GPP or NEE and a vegetation index  
\[\text{NEE or GPP}=\text{a*VI}+\text{b}, \tag{1}\]  
where NEE is the net ecosystem exchange, GPP is the gross primary production, VI is a vegetation index, a and b are regression constants. (ii) Model 2 assuming a direct linear relationship between GPP and the product of a VI and incoming PAR or Rg in this study  
\[\text{GPP}=\text{a*(VI*Rg)}+\text{b}, \tag{2}\]  
where GPP is the gross primary production, VI is a vegetation index, Rg is the short-wave radiation, a and b are regression constants. (iii) Model 3, a simplified LUE model, in which the LUE term is considered a constant (integrated in the model) and APAR is obtained by multiplying the fraction of absorbed photosynthetically active radiation (fAPAR), estimated as a linear function of a VI, by PAR (replaced by Rg in this study). By this approach, LUE and fAPAR estimates are conceptually independent. \[\text{GPP}=\text{LUE}^{*}\text{(a*VI}+\text{b)*Rg}, \tag{3}\]  
where GPP is the gross primary production, LUE is the light use efficiency, VI is a vegetation index, Rg is the short-wave radiation, a and b are regression constants. All three models were implemented for the entire growing season (single) or splitting the growing season in two phases (wet and dry). When stressed by low soil water content (SWC) or high vapor pressure deficit (VPD), plants close stomata, reducing transpiration and photosynthesis [57; 10; 58]. Therefore, SWC and VPD could be used as proxies for water stress. Based on our measurements, evapotranspiration was more sensitive to an increase in VPD than a decrease in SWC and a threshold value of 1500 Pa was considered to separate the wet phase from the dry phase of the growing season. A VPD value of 1500 Pa was also reported to be the optimal VPD for stomatal conductance in grasslands [59; 60]. Above this value, plant activity would reduce progressively. We therefore defined the wet 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 (June-August). The equation 3 was applied to each phase, assuming a constant LUE within the phase. In order to compare the performance of the different vegetation indices as proxies of carbon fluxes in the different regressions, three accuracy metrics were computed, namely the coefficient of determination (R\({}^{2}\)), the Root Mean Square Error (RMSE) and the Akaike Information Criterion (AIC). Additionally, in order to validate our models, we applied a bootstrap simulation with 1000 iterations to obtain the mean coefficient of determination (R\({}^{2}\)\({}_{\text{boot}}\)) with a 95% confidence interval for each regression [61; 62; 63]. The best models are the ones with a high value of R\({}^{2}\) and a low value of RMSE and AIC.