fields from the farm and therefore present high uncertainties (see Béziat et al., 2009). These
validations are presented in section I-4.5.1.1.

To evaluate the performance of the model in terms of CO₂ flux simulations we used the data
from the two ICOS sites (FR-Lam and FR-Aur). The use of very contrasted climatic years and
management regimes at FR-Lam and FR-Aur (see section I.4.3.2.2) makes our validation approach
more robust. Those results are presented in section I-4.5.1.2.

Among the objectives of the present study, assessing annual carbon budget is the most challenging. Indeed, the performances of the model to simulate both the NEP and the amount of biomass exported at harvest must be quite high to compute NECB as those two terms are usually of opposite sign and partly compensate each other. A small absolute error on one of those two terms will end up in a large absolute error on the NECB. Section I-4.5.2.1 presents the ability of the model to reproduce the cumulated NEE, which leads to the NEP. Then sections I-4.5.2.2 and I.4.6 present and discuss, respectively, the capability of the model to assess the NECB.

**4.5. Results**

**4.5.1. Evaluation of the overall model performances**

**4.5.1.1. GAI, DAM and yield estimates**

Figure I-7 shows an overview of the statistical performances of the model in estimating GAI and DAM by comparing the estimates with satellite-derived GAI and destructive biomass measurements for the FR-Aur (even-numbered years) and FR-Lam (odd-numbered years) sites in the 8 cultivation years. In addition, Figures I-8a and I-8b present scatter plots of these two variables, distinguishing observations performed during the entire studied period (2006-2014) at the two flux sites and measurements collected during the 2011 field campaign. At the flux sites, the GAI relative root mean square errors (RRMSEs, i.e., RMSE normalized by the mean observed value) range between 6 and 24% (for the years 2014 and 2007, respectively), while the determination coefficients (R²) are between 0.919 and 0.998 (for the years 2009 and 2014, respectively, Figure I-7). When also considering fields from the 2011 campaign (Figure I-8a), the GAI is reproduced by the model with a very high R² (0.97) and with almost no bias (0.004 m².m⁻²). The magnitude of this performance confirms that the model allows the correct interpolation of the GAI derived from remote sensing, as has already been demonstrated in previous studies (Duchemin et al., 2015a; Fieuzal et al., 2011a; Hadria et al., 2010a).

On average, the model also reproduces the DAM with good precision. When considering all
investigated fields across the studied years (Figure I-8b), the RMSE, RRMSE and R² were 201 g.m⁻²,
26.6% and 0.90, respectively. These performances are consistent with those achieved using the SAFY
model (maize and sunflower had RMSEs of 252 and 145 g.m⁻², respectively, and RRMSEs of 24.7 and
39.1%, respectively, in Claverie et al., 2012), as well as the performances of other models such as the
APSIM-wheat model (Asseng et al., 1998a) and STICS (Brisson et al., 2002a), which estimated winter
wheat biomass with R² of 0.90 and 0.78, respectively, and RMSEs of 0.101 and 0.266 kg.m⁻²,
respectively. Nevertheless, the performances are more scattered when analysing annual statistics
(Figure I-7). Indeed, the DAM RRMSE values are between 11 and 39.2% (for the years 2007 and
2013, respectively), while the R² values are between 0.945 and 0.998 (for the years 2012 and 2006,
respectively). In general, the modelling approach tends to underestimate the highest biomass values at
the end of the season (observed at LAM2011, AUR2012, LAM2013 and AUR2014). The highest
levels of errors, for the DAM in 2012 and 2013, correspond to the strongest vegetation developments
of 1960 g.m⁻² and 2298 g.m⁻², respectively, reached at the end of the crop season. Those DAM underestimations are caused by the underestimation of the effective GAI produced by BV-NNET for the highest values since the clumping effect is not accounted for (Claverie et al., 2012a). This issue will be further discussed in Section I.4.6.2.

<figure><img src="image_20.png" /></figure>

Figure I-7: Statistical performances (relative root mean square errors and coefficients of determination correspond to the bars and dots, respectively) associated with the estimation of the GAI (green) and DAM (red) for winter wheat at the Auradé (even-numbered years) and Lamasquère (odd-numbered years) sites.

<figure><img src="image_21.png" /></figure>

Figure I-8 : Comparison between the observed and the estimated crop variables (GAI, DAM and yield) over the crop period. Sites refer to FR-Lam and FR-Aur.
On average, yields (only data from the 2011 field campaign are analysed here, see I.4.5.1.1) are estimated with good precision (RMSE=1.02 t.ha⁻¹, RRMSE=21.5%, R²=0.78), but underestimations are observed for the highest observed values. These underestimations are directly related to the DAM underestimations described above, as yield is estimated as a fraction of the final biomass. Nevertheless, the yield-estimation performances of SAFY-CO₂ are similar to those of SAFY (Duchemin et al., 2008a) and STICS (Brisson et al., 2002a) for wheat, with R² values of 0.64 and 0.65, respectively, and RMSEs of 0.5 and 1.6 t.ha⁻¹, respectively.

### 4.5.1.2. Components of the net CO₂ fluxes: daily GPP, Reco and NEE

In this section, the components of the net CO₂ fluxes simulated by SAFY-CO₂ are compared to the measured NEE at the FR-Aur and FR-Lam flux sites and with the GPP and Reco estimated following the partitioning of NEE. For the 8 investigated cropping years, the model performances are evaluated in terms of the error (RMSE) and correlation (R²) (see Figure I-9). These statistical parameters are also calculated for the following periods of each cropping year (see Figure I-10):

1. From the beginning of the cropping year (October 1st) until the emergence of the crop

2. From the emergence of the crop until the maximum vegetation

3. From the maximum vegetation until harvest

4. From harvest until the end of the cropping year (September 30th)

#### 4.5.1.2.1. GPP estimates

When considering the whole vegetative periods (i.e., crop development and post-harvest vegetative events) of the 8 site-years, the simulated GPP dynamics agree well with the measurements, showing R² values between 0.82 (LAM2011) and 0.94 (LAM2007) and RMSEs between 1.34 (LAM2009) and 2.39 (LAM2011) gC.m⁻².d⁻¹.

<figure><img src="image_22.png" /></figure>

**Figure I-9** : Statistical performances (root mean square errors and coefficients of determination are bars and dots, respectively) associated with the estimations of the daily GPP, Reco and NEE for the 8 winter wheat cropping years (October 1st to September 30th) at FR-Aur (even-numbered years) and FR-Lam (odd-numbered years).
<figure><img src="image_23.png" /></figure>

**Figure I-10 :** Boxplot of the pluriannual determination coefficients (R²) and the pluriannual root mean square errors (RMSEs) of the simulated GPP, R<sub>eco</sub> and NEE in the before-crop, crop-growth, crop-senescence, and post-harvest periods for the 8 cropping years of winter wheat at AUR and LAM. Only the years with post-harvest vegetative events (2006-2009-2011-2013) are considered for the boxplot of the post-harvest GPP.

In 2011, however, the simulated GPP during crop development showed poor statistics. A strong development of vegetation occurred in early 2011, while no flux data were recorded between January 1st and March 31th. Consequently, the CO2 fluxes were gap filled, and thus, the simulated GPP is underestimated over this period, leading to low R2 and high RMSE values.

During crop development (growth and senescence, excluding post-harvest vegetative events), the simulated GPP shows R² values between 0.86 and 0.96 and RMSE values between 0.90 and 2.79 gC.m⁻².d⁻¹. These periods are well reproduced by the model, the growth period shows slightly better performances than the senescence period. This trend can be explained by two phenomena: i) the underestimation of the simulated GPP during senescence due to the abovementioned limits of remote-sensed GAI (underestimation due to saturation effects) coupled with the potential lack of remote sensing observations, which prevents the correct reproduction of GPP during the senescence phase and ii) the presence of weeds growing inside the senescent crop (the understorey vegetation often observed at FR-Lam), which are undetectable by the satellite but still impact the flux measurements (see *Temporal Evolution* of LAM2013).

Considering the post-harvest period, only site-years with significant vegetative events are considered (i.e., weeds for AUR2006 and LAM2009, spontaneous crop re-growth for LAM2011 and cover crops for LAM2013). The values of the R² between the simulated and observed GPP are lower during this period than during the crop development period for several reasons. First, at this stage, the same parameterization as the winter wheat parameterization is used for simulating the different types of post-harvest vegetative events (except concerning the parameters calibrated based on GAI satellite estimates – see section I.4.4.2). Second the post-harvest vegetation is characterized be higher spatial heterogeneity than winter wheat, leading to likely divergence between observed and simulated CO2 fluxes. Indeed, the footprint of the EC system is probably not fully representative of the entire field simulated by the model. Nevertheless, the overall performances of the model always increase when post-harvest vegetative event are considered.

#### 4.5.1.2.2. Reco estimates

Overall, the ecosystem respiration follows the same dynamics as the GPP but in an attenuated way (i.e., the differences between bare-soil and vegetative periods are smaller). The simulated R<sub>eco</sub> dynamics are in good agreement with the observations, as shown in Figure I-9. The annual RMSE
values range from 0.82 (LAM2009) to 1.67 (LAM2013) gC.m⁻².d⁻¹, and the R² values range from 0.58 (LAM2013) to 0.84 (LAM2014).

The mean annual performances of the model in estimating R<sub>eco</sub> are lower than those in estimating GPP but are consistent with other studies that aim to estimate ecosystem respiration in agricultural fields. Zhan et al. (2019) found R<sup>2</sup> values between 0.84 and 0.87 and RMSEs between 1.52 and 1.65 gC.m<sup>-2</sup>.d<sup>-1</sup> for a maize-soybean system, while Lohila et al., (2003) estimated the total ecosystem respiration of barley with a simple soil temperature-dependent model with good precision (R<sup>2</sup> between 0.71 and 0.79). A more detailed analysis of our results shows that the R<sub>eco</sub> is always well estimated during the crop growing periods, showing R<sup>2</sup> values between 0.88 and 0.95, while the statistical performances are lower in other periods than in the crop growing period (Figure I-10). One explanation for this difference is that the formulation for estimating heterotrophic respiration is too simplistic to reproduce some of the processes occurring in the soil, especially after harvest. This issue will be further discussed in Section I.4.6.2.

Additionally, for GPP, the presence of weeds growing inside the senescent crop, which is hardly detectable by the considered satellite but still impacts flux measurements (see the temporal evolution of LAM2013 in section temporal evolution), cannot be considered in our modelling approach. On the other hand, the good performances during the crop growing periods indicate that the model accurately reproduces autotrophic respiration, which can represent 80% of the R<sub>eco</sub> during the crop season (Béziat, 2009). The modelled R<sub>eco</sub> estimates are thus satisfactory considering the simplicity of its representation as well as the limited number of inputs in this crop modelling approach.

#### 4.5.1.2.3. NEE estimates

The daily dynamics of the NEE are well reproduced by the SAFY-CO₂ model even when GPP and R<sub>eco</sub> are underestimated, since the errors of the two components compensate for each other (either because the partitioning process overestimates the in situ data or because the model underestimates the crop development, which affects both GPP and R<sub>eco</sub>). The model that shows good performances over the cropping year in terms of errors and correlations, with RMSEs between 1.09 (AUR2006) and 1.59 (LAM2011) gC.m<sup>-2</sup>.d<sup>-1</sup> and R<sup>2</sup> values between 0.78 (LAM2011) and 0.90 (AUR2012). The model achieves better performances during the vegetative stages than during fallow periods (as is the case for GPP and R<sub>eco</sub> due to the reasons mentioned above).

In order to compare the performances of SAFY-CO₂ to simulate NEE with other agronomical or land surface models, our results were confronted to those presented by Wattenbach et al., (2010). In their study, they compared the performances of DNDC (Li et al., 2005, 1994, 1992), ORCHIDEE-STICS (de Noblet-Ducoudré et al., 2004; Gervois et al., 2008b), SPA (Williams et al., 1996) and CERES-EGC (Gabrielle et al., 2006; Lehuger et al., 2009) in reproducing the GPP, R<sub>eco</sub> and NEE for several site-years, including the AUR2006 crop season. SPA and SAFY-CO₂ outperformed the other approaches by reproducing the dynamics of the cumulated NEE, with R<sup>2</sup> values of 0.993 and 0.995, respectively. In terms of errors, SAFY-CO₂ and CERES-EGC showed the lowest RMSEs, 33.6 gC.m<sup>-2</sup> and 44.16 gC.m<sup>-2</sup>, respectively.

#### 4.5.1.2.4. Temporal evolution

We decided to show the performances of SAFY-CO₂ in simulating the time courses of the GPP, R<sub>eco</sub> and NEE for three contrasting site-years (AUR2006, LAM2009 and LAM2013) among the eight cultivation years. The objectives here are i) to identify potential sources of errors in the GPP and R<sub>eco</sub> estimates that can affect the NEE and the net annual CO₂ fluxes (NEP) and therefore the NECB estimates, ii) to verify whether the proposed approach is robust for varying soils, management
practices and climatic years and iii) to analyse the potential of this approach to simulate contrasting post-harvest vegetative events.

<figure><img src="image_24.png" /></figure>

**Figure I-11** : Temporal evolution of the measured (in blue) and estimated (in red or red dashed lines) GPP, Reco and NEE for 2 site-years (AUR2006 and LAM2009). The red/red dashed lines represent the simulations that do/do not account for re-growth and weed events. The yellow envelopes represent the daily standard deviation of the 10 (/30) best simulations (i.e. smaller RMSE GAI error). The vertical dashed lines define the cropping year.
<figure><img src="image_25.png" /></figure>

**Figure I-12** : Temporal evolution of the measured (in blue) and estimated (in red or red dashed lines) cumulated GPP, R<sub>eco</sub> and NEE for LAM2013. The red/red dashed lines represent the simulations that do/do not account for the cover crop. The yellow envelopes represent the daily standard deviation of the 10 (/30) best simulations (i.e. smaller RMSE GAI error). The vertical dashed lines define the cropping year.

* Auradé 2006

The 2006 cropping year at FR-Aur (grain farm) is characterized by very clear sky conditions and strong radiation, little precipitation in spring (23.4, 29.7 and 32.8 mm of rain during April, May and June, respectively, in 2006; in contrast, the monthly means in these months over the 8 years of the study and both sites are 64.1, 61.1 and 95.4 mm, respectively) and several re-growth/weed events occurring after harvest. For this site-year (Figure I-11), the model correctly reproduces the GPP and the R<sub>eco</sub> in terms of errors (RMSEs of 1.38 and 0.87 gC.m<sup>-2</sup>.d<sup>-1</sup>, respectively) and dynamics (R<sup>2</sup> of 0.92 and 0.78, respectively). Consequently, the NEE for this year is accurately estimated (RMSE and R<sup>2</sup> of 1.09 and 0.88, respectively). Nevertheless, the modelled GPP is slightly overestimated after maximum development and at the beginning of senescence (Figure I-11). This overestimation could be related to the water stress conditions observed after the maximum GAI was reached (Béziat et al., 2013) that cannot be fully considered by the model. Moreover, spontaneous re-growth and weeds developed twice after harvest. The first event led to increases in the observed GPP and R<sub>eco</sub> just after harvest and was interrupted by soil work on July 31. The second vegetative event occurred from mid-August until late September 2006 and was interrupted by soil work on September 29. In a first attempt, we simulate only one vegetation cycle after harvest, so the two events are simulated as one (see red dashed line in Figure I-11); as a consequence, GPP and R<sub>eco</sub> are overestimated. NEE estimates during this period are improved (RMSE from 1.22 to 0.93 gC.m<sup>-2</sup>.d<sup>-1</sup>) by accounting for weed/re-growth development.

* Lamasquère 2009

The 2008-2009 cropping season at FR-Lam (milk and chicken production farm) was characterized by strong rains in November and December that saturated the soil, causing poor emergence and late winter wheat development. Additionally, weeds developed before harvest during
the winter wheat senescence. The GPP dynamics for LAM2009 are well reproduced, with an R² of 0.89. Additionally, R eco and NEE present R² values of 0.78 and 0.81, respectively, over the cultivation year. The errors of the GPP, R eco and NEE are also low, showing RMSE values of 1.34, 0.82 and 1.14 gC.m⁻².d⁻¹, respectively. However, during May 4, peaks are observed in measurements (13 days in total) that are not reproduced by the model. This divergence between simulated and observed GPP comes from the underestimation of the fELUE for days with high radiation and very clear sky conditions. After harvest (from mid-July to mid-September), the presence of weeds in the field is highlighted by the measured GPP dynamics. As for 2006, the dynamics and the range of the simulated post-harvest GPP and R eco are not correctly reproduced by the model since the parameterization is the same as that for winter wheat, and many phenomena are not considered (i.e., priming effect, nitrogen or water stress, etc.). Nevertheless, accounting for weeds in the model allows a better estimation of the NEE (RMSEs decreasing from 1.7 to 1.14 gC.m⁻².d⁻¹).

* Lamasquère 2013

We present the results of FR-Lam over the LAM2013 cropping year because it is the only site-year during which a cover crop was grown. Indeed, in 2013 at FR-Lam, after the harvest of the winter wheat, white mustard was sown on the 21st of August and incorporated in the soil on the 4th of December. Unlike re-growth and weed development, the development of the cover crop is rather homogeneous in the field and follows a growing cycle that is correctly reproduced by the model, in addition to the CO₂ fluxes (RMSE = 0.68 gC.m⁻².d⁻¹ and R² = 0.88 for GPP; RMSE = 1.14 gC.m⁻².d⁻¹ and R² = 0.62 for NEE).

LAM2013 is also marked by an early winter wheat development that benefited from good climatic and soil conditions for emergence. Therefore, the final biomass is high (2298 g.m⁻²) compared to that of other years (mean and standard deviation over the 8 studied years: 1566 +/- 453 g.m⁻²). This site-year is also characterized by weeds that developed during the senescence of winter wheat. These weeds could not be observed by the considered satellites, and the model was not able to simulate their effects on the CO₂ fluxes (particularly on the GPP and the R eco) observed just before and after harvest. The difference between the observed and simulated NEE dynamics is small, either because the “observed” GPP and R eco partly balanced each other or because the increases in the “observed” R eco and GPP are caused by errors in the NEE partitioning process.

<figure><img src="image_26.png" /></figure>

**Figure I-13** : Temporal evolution of the measured (in blue) and estimated (in red and red dashed lines) daily values of the GPP, R eco and NEE for 2 site-years: a) AUR2006 and b) LAM2009. The red/red dashed lines represent the simulations that do/do not account for re-growth and weed events. The yellow envelope represents the daily standard deviation of the 10 (/30) best simulations (i.e. smaller GAI RMSE).
4.5.2. From the cumulated NEE to yearly carbon budget

4.5.2.1. Analysis of the cumulated NEE dynamics

For the sake of conciseness, the analysis of the temporal behaviour of the cumulated NEE focuses on the same site-years ( i.e. , AUR2006 and LAM2009). These values are presented from October 1st until September 30th of the following year because this period corresponds to the agricultural cropping year in our area. The analysis of the cumulated NEE dynamics measured by the flux towers and modelled by SAFY-CO₂ allowed for the identification of the CO₂ net assimilation and release phases (Figure I-13). A negative slope in the cumulated NEE curve ( i.e. , corresponding to net assimilation, with GPP>R eco ) is observed during the growing season and during crop re-growth or weed or cover crop development. A positive slope ( i.e. , corresponding to CO₂ release, with GPP<R eco ) is observed during the bare-soil periods and the senescence stages.

The modelled cumulated NEE for AUR2006 slightly diverges from the observations before the growing season because of the underestimation of the simulated R eco during this period. The effect of the overestimation of the GPP at the end of the growing stage and beginning of the senescence stage (as discussed above) leads to an underestimation of the maximum simulated cumulated NEE value.

Finally, the slope becomes positive after harvest (when there is no more GPP). The modelled cumulated NEE is well estimated when re-growth and weeds considered (red dashed line) and diverges when they are not considered in the model (red solid line).

In 2009, the simulated cumulated NEE matches well with the observations before and during the vegetative period. Indeed, over this period, the model reproduces the GPP and the R eco well, leading to a good estimation of the NEE. After harvest, not considering weeds causes the model to diverge from the observations, as in 2006, while modelling post-harvest vegetation development attenuates this bias (even if the simulated weed growth starts later here than on the field).

When the post-harvest vegetative events are simulated, the difference between the observed and the simulated NEP is improved from 95 and 140 gC.m⁻² to 33 and 16 gC.m⁻² for AUR2006 and LAM2009, respectively.

These results emphasize the need to include the effects of re-growth events, weeds and cover
crops in the model because they represent non-negligible contributions to the NEP, which in turn
could have an important impact on the final annual NECB values.

4.5.2.1.1. Carbon budget over 8 agricultural seasons of winter wheat

The NECB model estimates and its components (NEP, C<sub>exp</sub>) are compared with those of the eight site-year measurements (Table I-7). In this table, the effects of post-harvest vegetative events on NEE and NECB are considered. For FR-Lam, the C<sub>inp</sub> is prescribed for the calculation of the modelled NECB, as the amount of C input as organic manure cannot be estimated by remote sensing. Also, C<sub>inp</sub> corresponding to the amount of seeds brought to the plots are prescribed for both sites (even if very small compared to the other terms). The NEP values estimated by SAFY-CO<sub>2</sub>, which vary from -191 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2009) to -486 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2011), are generally close to the measured values, which vary from -208 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2009) to -410 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2011). Table I-7 shows that for all site-years, winter wheat is a CO<sub>2</sub> sink, and the model is able to reproduce the inter-annual variability in this sink activity. The years showing the largest difference in terms of NEP are AUR2014, AUR2012 and LAM2013, with differences of 30%, 24% and 22%, respectively, compared to observations. For AUR2014, this difference can be explained in different ways. First, the senescent phase is not well reproduced by the model because only one satellite image was available during this
period, leading to an overestimation of the GPP. Second, after harvest, an increase in R<sub>h</sub> is observed and is not reproduced by the model. This is due to the incorporation of straw into the soil, which leads to an increase in the soil microbial activity. The difference observed for AUR2012 is partly due to two post-harvest increases in R<sub>h</sub> (priming effect) that could not be reproduced by the model. The first event, occurring at the beginning of July, was induced by rainfall. The second event followed ploughing that occurred at the beginning of August. Moreover, the NEP values of two of the three site-years where significant re-growth vegetative events occurred (during the cultural year) are better estimated once the re-growth is considered. Indeed, for AUR2006 and LAM2009, the differences between the simulated and observed NEP values are -30 and -68%, respectively, before taking re-growth into account and +10 and -7.8%, respectively, after taking re-growth into account. For LAM2011, NEP is first underestimated (-21%) and then overestimated (+19%).

<table><caption>**Table I-7 : Annual net ecosystem carbon budgets (NECB) and their components (NEP, Cinp, Cexp) derived from the in-situ (OBS) and modelled (SIM) data for 8 site-years. Uncertainties on observations are also shown (for more details see Béziat et al., 2009).**</caption><thead><tr><th rowspan="2">Site-Year</th><th colspan="4">NEP</th><th colspan="4">Cinp</th><th colspan="4">Cexp</th><th colspan="4">NECB</th></tr><tr><th>SIM</th><th>OBS<br>[gC/m²/yr]</th><th>Diff</th><th>Diff<br>[%]</th><th>OBS<br>[gC/m²/yr]</th><th>SIM</th><th>OBS<br>[gC/m²/yr]</th><th>Diff</th><th>Diff<br>[%]</th><th>SIM</th><th>OBS<br>[gC/m²/yr]</th><th>Diff</th><th>Diff<br>[%]</th><th>SIM</th><th>OBS<br>[gC/m²/yr]</th><th>Diff</th><th>Diff<br>[%]</th></tr></thead><tbody><tr><td>AUR2006</td><td>-355</td><td>-322 ± 20</td><td>-33</td><td>10</td><td>-6.25 ± 1.88</td><td>270</td><td>216 ± 56</td><td>54</td><td>25</td><td>-94</td><td>-113 ± 60</td><td>19</td><td>-17</td><td>-302</td><td>-371 ± 33</td><td>69</td><td>-19</td></tr><tr><td>LAM2007</td><td>-302</td><td>-371 ± 33</td><td>69</td><td>-19</td><td>-389 ± 95.4</td><td>387</td><td>322 ± 29</td><td>65</td><td>20</td><td>-304</td><td>-439 ± 105</td><td>135</td><td>-31</td><td>-191</td><td>-208 ± 19</td><td>17</td><td>-8</td></tr><tr><td>LAM2009</td><td>-253</td><td>-301 ± 47</td><td>48</td><td>-16</td><td>-6.25 ± 1.88</td><td>293</td><td>279 ± 25</td><td>14</td><td>5</td><td>-48</td><td>-78 ± 55</td><td>30</td><td>-38</td><td>-486</td><td>-410 ± 45</td><td>-76</td><td>19</td></tr><tr><td>AUR2010</td><td>-253</td><td>-301 ± 47</td><td>48</td><td>-16</td><td>-6.25 ± 1.88</td><td>253</td><td>204 ± 53</td><td>49</td><td>24</td><td>-6</td><td>-102 ± 71</td><td>96</td><td>-94</td><td>-362</td><td>-293 ± 34</td><td>-69</td><td>24</td></tr><tr><td>LAM2011</td><td>-486</td><td>-410 ± 45</td><td>-76</td><td>19</td><td>-166 ± 49.8</td><td>400</td><td>355 ± 32</td><td>45</td><td>13</td><td>-252</td><td>-221 ± 74</td><td>-31</td><td>14</td><td>-421</td><td>-345 ± 32</td><td>-76</td><td>22</td></tr><tr><td>AUR2012</td><td>-362</td><td>-293 ± 34</td><td>-69</td><td>24</td><td>-6.25 ± 1.88</td><td>326</td><td>223 ± 20</td><td>103</td><td>46</td><td>-41</td><td>-76 ± 40</td><td>35</td><td>-46</td><td>-316</td><td>-243 ± 34</td><td>-73</td><td>30</td></tr><tr><td>LAM2013</td><td>-421</td><td>-345 ± 32</td><td>-76</td><td>22</td><td>-178 ± 53.4</td><td>436</td><td>488 ± 40</td><td>-52</td><td>-11</td><td>-163</td><td>-36 ± 74</td><td>-127</td><td>355</td><td>-316</td><td>-243 ± 34</td><td>-73</td><td>30</td></tr><tr><td>AUR2014</td><td>-316</td><td>-243 ± 34</td><td>-73</td><td>30</td><td>-6.25 ± 1.88</td><td>285</td><td>214 ± 67</td><td>71</td><td>33</td><td>-31</td><td>-29 ± 75</td><td>-2</td><td>6</td><td>-316</td><td>-243 ± 34</td><td>-73</td><td>30</td></tr><tr><td>SIM vs. OBS</td><td></td><td>RMSE [gC/m²/yr]</td><td>61</td><td></td><td></td><td>RMSE [gC/m²/yr]</td><td>61</td><td></td><td></td><td>RMSE [gC/m²/yr]</td><td>77</td><td></td><td></td><td>RMSE [%]</td><td>-19.7</td><td></td><td></td></tr><tr><td>All Sites</td><td></td><td>RMSE [%]</td><td></td><td></td><td></td><td>RMSE [%]</td><td></td><td></td><td></td><td>RMSE [%]</td><td>-56.1</td><td></td><td></td><td>R²</td><td>0.58</td><td></td><td></td></tr><tr><td></td><td></td><td>R²</td><td></td><td></td><td></td><td>R²</td><td></td><td></td><td></td><td>R²</td><td>0.66</td><td></td><td></td><td>R²</td><td>0.83</td><td></td><td></td></tr></tbody></table>

The analysis of the amount of carbon exported from the ecosystem at harvest (C<sub>exp</sub>) shows that this amount varies considerably from one site-year to another. The simulated C<sub>exp</sub> varies from 253 gC.m<sup>-2</sup>.yr<sup>-1</sup> (AUR2010) to 436 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2013), while the observations range from 204 gC.m<sup>-2</sup>.yr<sup>-1</sup> (AUR2010) to 488 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2013). C<sub>exp</sub> values are often larger at FR-Lam than at FR-Aur due to the export of grain and straw and the model tends to overestimate the C<sub>exp</sub>. The differences between the observations and model estimates of C<sub>exp</sub> vary between 5% and 46%. However, the comparison of the modelled C<sub>exp</sub> and the observed C<sub>exp</sub> should be performed with caution since the precision of the yield provided by the farmer is questionable (especially at the Lamasquère site where most of the straws are exported) since those values are averaged over several fields of the farm (see section I.4.5.1.1). For this reason, Béziat et al. (2009) concluded that the uncertainty of in situ C<sub>exp</sub> is often higher than the uncertainty of the NEP at our sites. For AUR2012, for instance, the C<sub>exp</sub> estimated from the farmer's data is 223 gC.m<sup>-2</sup>.yr<sup>-1</sup>, while our destructive measurements encompassing over 30 subplots of 3.75*10<sup>-2</sup> m<sup>2</sup> in the field suggest a C<sub>exp</sub> of 406+/-53 gC.m<sup>-2</sup>.yr<sup>-1</sup>. When analysing the performance of the model against our own destructive field samples, the model generally performs much better (see Figure I-8b and I-8c). This indicates that our modelling approach may perform better for estimating C<sub>exp</sub> (and therefore NECB) than what is presented in Table I-7.

The NECB estimated from in situ data or from SAFY-CO<sub>2</sub> characterize all of the site-years as carbon sinks. The modelled NECB estimates vary from -304 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2007) to -6 gC.m<sup>-2</sup>.yr<sup>-1</sup> (AUR2010), while the in situ NECB estimates vary from -439 gC.m<sup>-2</sup>.yr<sup>-1</sup> (LAM2007) to -29 gC.m<sup>-2</sup>.yr<sup>-1</sup> (AUR2014). The FR-Lam site-years present the greatest variations between years, and the carbon inputs (NEP + C<sub>inp</sub>) are stronger than those of FR-Aur, partly because of organic fertilization.
As shown in Table I-7, the relative differences between the modelled and in-situ NECB vary from 6% (AUR2014) to 355% (LAM2013). In absolute terms, these differences vary from 2 gC.m⁻².yr⁻¹ (AUR2014) to 135 gC.m⁻².yr⁻¹ (LAM2007) but the sign of NECB is always similar between both approaches. Regarding all simulated years, the model shows RMSE of 77 gC.m⁻².yr⁻¹. Note however that, the modelled NECB match to that observed for five simulated years out of eight if the uncertainties are considered. Indeed, the uncertainties on the in-situ NECBs (derived from uncertainties on the observed NEP, C<sub>inp</sub>, C<sub>exp</sub>) range from 40 to 105 gC.m⁻².d⁻¹. These results highlight the importance of precisely estimating each of the terms that compose the NECB to obtain accurate estimations of the annual crop carbon budgets (both with the in-situ and the modelling approaches).

## 4.6. Discussion

### 4.6.1. Performances and benefits of our approach

In this study, our objective is to evaluate the potential of high resolution GAI products assimilation into a simple crop model for simulating the biomass, the yield, the net CO₂ fluxes components and the annual C budget of winter wheat crops at plot scale. The main advantage of this approach is that it requires few input data and little or no external information about management practices. Also, even with a limited number of equations and parameters compared to more complex crop models that require data on management practices, SAFY-CO₂ achieves equivalent or better performances regarding estimates of key components of the C budget: CO₂ fluxes, biomass and yield (see sections I-4.5.1.1 and I-4.5.1.2).

Next, we demonstrated the ability of the model to reproduce winter wheat dynamics, production and CO₂ fluxes under contrasted climatic and management conditions with the same parametrization. It shows the ability of such a remote sensing driven diagnostic approach (e.g. for calculating GPP) to account implicitly for the main stresses (N, drought, temperature...) and the main crop development limiting factors.

Finally, our methodology allows accounting for the effect of post-harvest spontaneous re-growth, weeds and cover crops on the CO₂ fluxes. As showed by Ceschia et al., (2010) this is essential for estimating accurately cropland C budgets and only remote sensing based approaches allow characterizing the dynamics and the spatial heterogeneity of the various post-harvest vegetative events. Therefore, in spite of a generic parameterization of the SAFY-CO₂ model for those post-harvest vegetative events, the overall performances of the model for simulating CO₂ fluxes and C budgets always increase when they are accounted for. The performances should improve thanks to the higher temporal resolution of the Sentinel missions. Next step could be to apply a specific parametrization, depending on the nature of those events or on their species composition (e.g. for cover crops), provided that the information is given by the farmer or can be retrieved by remote sensing (e.g. through cover crop classification).

### 4.6.2. Potential limitations of this approach and drawbacks for large scale application

The first main limitation of this approach based on optical remote sensing is that gaps in optical remote sensing observations during crucial periods of the crop development could lead to wrong estimates of the GAI dynamics, biomass and CO₂ fluxes or could even make our approach inoperative (e.g. in 2007-2008). Fortunately, recent Sentinel 2 satellite missions provide observations at high spatial resolution (10m) every 5 days all over the globe which could partially solve this problem. Also,
it was shown that the combined use of optical and Synthetic Aperture Radar (SAR) satellite data, like Sentinel 1, can overcome this issue (see Ameline et al., 2018; Baup et al., 2019; Betbeder et al., 2016; Remy Fieuzal et al., 2017; Revill et al., 2013). Indeed, the signal of the SAR satellites is not affected by clouds and they can even observe the surface at night. Another limitation of this optical remote sensing approach is that it cannot detect understorey vegetation (e.g. weeds) and their effect on the CO2 fluxes. Here again SAR data may overcome this issue as microwave signal are associated to deeper penetration capabilities (compared to optical reflectance), depending on the considered wavelength, and providing a valuable information on vegetation structure and water content (Brown et al., 2003).

The second main limitation of this approach concerns the availability of plot scale information regarding straw management and organic fertilization. Those practices cannot be detected or quantified by remote sensing at this stage and therefore the uncertainty on the C budgets estimates in areas where animal farming occurs is high. This issue may be overcome in the future if data from the Farm Management Information Systems (FMIS) become more easily and more widely accessible. Another limitation for applying our approach concerns the size and the shape of the agricultural plots. As mentioned above, the contours of the plots must be eroded so that the signal is not influenced by surrounding landscape elements. Thus, we consider that for plots below 0.5-1 ha, GAI products may not be of good enough quality to apply our approach. Also, high resolution GAI and crop maps are needed in our approach. In Europe, the later can be obtained via the Land Parcel Identification System (LPIS) and both data inputs should be available in a near future via the High Resolution Layers Copernicus Land Monitoring Service (https://land.copernicus.eu/pan-european/high-resolution-layers). The last limitation of our approach is that it cannot be used for forecasting, since it is based on satellite observations, although it is possible to test the effect of some scenarios on the C budgets (e.g. accounting or not for the effects of the post-harvest vegetative events or for the impact of exporting or not the straw from the plot).

Of course, the question of the transposability of our approach and of its domains (spatial and temporal) of validity should be considered with caution. Indeed, the current parametrization of the SAFY-CO2 model is adapted to the pedoclimatic conditions where it has been set and it should be adapted to other crop species. Also, the boundaries of the calibrated winter wheat phenological parameters are set for our pedoclimatic conditions. Applying this approach to areas where winter wheat has different periods of emergence and senescence would require to redefine those boundaries (e.g. in Northern countries). Note however that 1) the future High Resolution Phenology Copernicus Land Monitoring Service should provide, all over Europe since 2017, the dates of emergence and end of the growing for the crops and cover crops at plot scale for each cropping year by the end of this year and that 2) the transposability of the original SAFY model has already been tested in contrasted pedoclimatic conditions (France, Mexico and Morocco; see (Claverie et al., 2012; Duchemin et al., 2015, 2008; respectively) and for different crop species (corn, soybean, sunflower and wheat).

Also, the parametrisation of the Q₁₀ based approach for estimating heterotrophic respiration is well adapted to the type of soils and climates similar to the ones found in our area of study but it should probably be adapted to other soil types (e.g. organic or sandy soils) and climatic conditions. Still, with a similar approach, Delogu et al., (2017) obtained good R<sub>h</sub> estimates (RMSE comprised between 0.15 and 0.73 gC.m⁻².d⁻¹ and R² between 0.42 and 0.92 depending on the site) over contrasted pedoclimatic conditions. Another potential issue concerning our approach relates to the simplistic method for estimating R<sub>h</sub> which should be considered as a first step for estimating R<sub>h</sub>. In the future, this method could be improved, with little changes in the formalisms, by considering a R<sub>10</sub> parameter that depends on top soil slow carbon content as it is proposed in Delogu et al., (2017). Another step of
improvement could be to account for the priming effect following the incorporation of fresh organic
matter into the soil (Kuzyakov et al., 2000a). Indeed, with our modelling approach we could already
estimate the amount of crop residues, cover crop, weeds and spontaneous re-growth incorporated in
the soil. However, accounting for the effect of organic fertilization is not an option at this stage, since
this kind of information is not yet available at plot scale over large areas. For similar reasons, we did
not account for the effect of soil work on soil respiration and also because it was shown that it has no
significant effect on soil respiration (Eugster et al., 2010b) and no clear effect on SOC mineralisation
(e.g. in Dimassi et al., 2014; Powlson et al., 2016; Virto et al., 2012). Also, a potential drawback of
our approach for estimating R<sub>h</sub> is that the effect of the soil water content is not accounted for. As for
the SAFY model (Battude et al., 2017a; Duchemin et al., 2015a), we have already tested the coupling
of a soil water module (FAO56; Allen et al., 1998) to the SAFY-CO<sub>2</sub> model (see Veloso, 2014). Such
a coupling allows accounting for soil water content effects on photosynthesis and on heterotrophic
respiration, and requires that accurate data concerning soil properties (e.g. texture, depth) at the plot
scale are available. Unfortunately, the current products mapping soil properties (e.g. GlobalSoilMap,
SoilGrids) either have a too coarse resolution for our area of study (250 m, for SoilGrids) or have too
little accuracy to meet the needs of a SAFY-CO<sub>2</sub> – soil water module coupled approach at plot scale
(e.g. the performances of GlobalSoilMap in France are R<sup>2</sup> = 0.27 and RMSE = 128 g.kg<sup>-1</sup> for clay
content). Therefore, in the perspective of up-scaling our current approach for estimating annual
cropland C budgets, we chose to rely only on currently available and sufficiently accurate data at plot
scale. This is the reason why we chose a simple Q<sub>10</sub> approach for estimating heterotrophic respiration
at this stage. Of course, if this choice is likely suitable for plot scale annual C budget estimates, at crop
rotation scale or over longer periods of study, the coupling of the SAFY-CO<sub>2</sub> model with a soil
organic matter model (e.g. RothC, Coleman and Jenkinson, 1996; AMG, Saffih-Hdadi and Mary,
2008) should be considered. Such a step would benefit from 1) the improvement of the current soil
products (e.g. GlobalSoilMap) which could be achieved by developing the current methods of
inversion based on high resolution multi-spectral or hyperspectral remote sensing data (see Castaldi et
al., 2019; Vaudour et al., 2019) and 2) from an easier and more systematic access to the FMIS data.

Another limitation of our modelling approach is that it tends to underestimate the highest
biomass values at the end of the season and therefore also the yield, impacting the C budgets
estimates. These underestimations may be partly due to the underestimation of satellite-derived GAI
(Claverie et al., 2012a). In such a context, the assimilation of both GAI derived from optical images
and dry biomass estimated from SAR images into the agro-meteorological model overcomes the
limitation, as presented by Betbeder et al., (2016) in the specific case of soybean. Furthermore, taking
into account the clumping effect in radiative transfer model would make it possible to limit the
underestimation of GAI values (derived from optical images) when vegetation becomes dense.

Still, we show that, within the limit of its domain of application, our approach was able to
reproduce correctly the GAI, biomass and CO₂ flux dynamics and it was able to estimate the NEP with
a satisfactory level of accuracy. The relatively large error of prediction on the C budgets (mean RMSE
of 77 gC.m⁻².yr⁻¹ and rRMSE of 56 %) has to be tempered considering the uncertainties on the NECB
calculated from the in-situ data. Indeed, besides the inherent uncertainties on the NEP associated to
measurements errors and data processing, there is a strong uncertainty on the in-situ Cexp term for the
two ICOS sites as mentioned previously. Considering the resulting uncertainties on the in-situ NECB,
we conclude that the modelled NECB match the observations for five years out of eight. Also, section
4.1.2.4 showed that the model was able to estimate the yield with a rather good precision, suggesting
that the error on the NECB could be reduced when comparing our estimates with in-situ NECB
calculated with more accurate yield data. Still, even if SAFY-CO₂ provided accurate estimations of the
annual components of the NECB we cannot claim at this stage that the model can reproduce accurately carbon budgets, especially over the long term. More accurate in-situ data, a larger dataset of validation and/or simulations on longer periods evaluated against estimates of soil C stock changes based on soil analysis would be needed to conclude.

Of course, other approaches allowing to estimate carbon budget exist, such as soil organic matter models which are designed to simulate the evolution of soil C stocks. The two most widely used and validated SOM are Roth-C (Coleman and Jenkinson, 1996) and CENTURY (Parton et al., 1987). Those models estimate soil C stock changes have been evaluated against long term experiments. Contrary to the proposed approach, they need information about soil texture, management practices or residue quality. Their relative error in estimating soil C stock changes is comprised between 2-30% (Falloon and Smith, 2006; Guo et al., 2007; Smith et al., 1997) for Roth-C model and between 1.8-16.4% for CENTURY (Cong et al., 2014; Falloon and Smith, 2003). These results, which are more accurate than those achieved with SAFY-CO₂ model should be tempered by the fact that they represent two different approaches, requiring different input and designed for different purposes and time-scales.

Also, as in other studies (Ceschia et al., 2010b; Schmidt et al., 2012) our results showed that for all cropping years the plots behave as a net CO₂ sinks and our results concerning the potential C storage of winter wheat crop are consistent with other studies (e.g. Aubinet et al., 2009; Ceschia et al., 2010). In addition, our results show that, in soils with low SOC content, post-harvest vegetative events (e.g. cover crops) increase soil organic carbon storage which is consistent with other studies (Kaye and Quemada, 2017; Pellerin et al., 2019; Poeplau and Don, 2015; Tribouillois et al., 2018).

Finally, in spite of the limitations and potential drawbacks of this approach, it seems to be a good compromise for estimating the components of the annual C budgets over large areas at this stage and we think that it offers great perspectives of development and applications at large scale thanks to the new satellite missions and Copernicus services.

## 4.7. Conclusion

In this work, we demonstrate the potential of high-resolution remote sensing data assimilation in a semi-physical crop model (SAFY-CO₂) to successfully provide estimates of some of the main components of cropland annual carbon budgets (i.e., net CO₂ flux components and yield). While this modelling approach is promising because it requires few input parameters and no management data for estimating crop production and net CO₂ fluxes, this approach should be considered a first step for filling the gap in obtaining spatially explicit representation of the main components of cropland carbon budgets at the regional scale for a crop rotation or longer. Indeed, the main limitation of this approach is that, in areas concerned with animal farming, the calculation of the carbon budget requires data on i) organic amendments and ii) the fraction of straw exported at harvest, which presently cannot be retrieved by remote sensing at this stage. The second main limitation relates to the fact that the simple Q₁₀ based approach for estimating heterotrophic respiration does not allow us to estimate accurately the C budget for periods longer than the cropping year. For longer periods of study, the benefit of coupling our model with a soil module should be investigated. Another limitation concerns the availability of satellite observations, since our approach is data driven. However, because of recent HTRS satellite missions (Sentinel 2 and Landsat-7&8), this type of approach could be generalized and more accurate and robust. Synthetic aperture radar satellites (e.g., Sentinel 1) could also be used to overcome cloudy conditions (Veloso et al., 2017). In addition, our results show that the performance of the model in estimating net CO₂ fluxes and thus C budgets are significantly improved by
considering the development of weeds and crop re-growth after harvest. These events, as well as the
presence of cover crops in crop rotations, are rarely or never accounted for in regional or global
modelling of CO2 fluxes, although they significantly impact cropland carbon budgets.

In the perspective of future global-scale applications, our approach could be strengthened (validated for a wider range of climates and management regimes) and extended to other crops by using data from international flux networks (e.g., ICOS and FLUXNET) and from recent HTRS satellite missions.
## 4.8. Supplementary material

### 4.8.1. Effect of diffuse radiation on the ELUE

The fraction of diffuse radiation affects photosynthesis (Béziat et al., 2009a) and should be accounted for when simulating GPP. Furthermore, because the photosynthetic rate of leaves is usually saturated under high incoming radiation, leaves with lower irradiance will be more efficient than those with higher irradiance, and a reduction in the volume of shade leaves within the canopy should result in an increase in the efficiency of the canopy in the presence of low and diffuse radiation (Roderick et al., 2001). Thus, the photosynthetic efficiency is expected to increase as the diffuse solar radiation increases. To quantify this effect, we used data from 5 European flux sites, including the Lamasquère and Auradé sites. The 3 other sites were Lonzee (LON) in Belgium, Grignon (GRI, located near Paris) in France and Oensingen (OEN) in Switzerland. For all sites, the ELUE increased non-linearly with the ratio of diffuse over total global radiation (Fig. I-14). As a consequence, the relationship between the effective light-use efficiency and the ratio between diffuse and direct radiation at ground level was defined as an exponential function, with parameter "b" fixed to 1.34 and parameter "a" calibrated based on the assimilation of GAI derived from satellite observations. The relationship seemed relatively generic; the correlation coefficient of the regression was 0.63.

<figure><img src="image_27.png" /><figcaption>Figure I-14 : Relationship between the ELUE and the ratio between diffuse (Rdf) and global (Rg) radiation measured at 5 European flux sites (Lonzée, Grignon, Oensingen, Lamasquère & Auradé).</figcaption></figure>

### 4.8.2. Temporal evolutions of the in-situ and simulated net CO₂ flux components

The purpose of this section is to compare the temporal evolutions of in-situ and simulated CO₂ fluxes (GPP, R<sub>eco</sub>, NEE and cumulated NEE) at FR-Lam and FR-Aur for the site years discussed in the text but not shown.
<figure><img src="image_28.png" /><figcaption>Figure I-15 : Temporal evolution of the measured (in blue) and estimated (in red) GPP, R<sub>eco</sub> and NEE, top, and cumulated NEE, bottom, for LAM2007. The yellow envelopes represent the daily standard deviation of the 10 (/30) best simulations (i.e. smaller RMSE GAI error). The vertical dashed lines define the cropping year.</figcaption></figure>
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

Figure I-16 : Temporal evolution of the measured (in blue) and estimated (in red) GPP, R<sub>eco</sub> and NEE, top, and cumulated NEE, bottom, for AUR2010. The yellow envelopes represent the daily standard deviation of the 10 (/30) best simulations (i.e. smaller RMSE GAI error). The vertical dashed lines define the cropping year.
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

Figure I-17 : Temporal evolution of the measured (in blue) and estimated (in red or red dashed lines) GPP, $R_{eco}$ and NEE, top, and cumulated NEE, bottom, for LAM2011. The red/red dashed lines represent the simulations that do/do not account for re-growth and weed events. The yellow envelopes represent the daily standard deviation of the 10 (/30) best simulations (i.e. smaller RMSE GAI error).
<figure><img src="image_33.png" /></figure>

Figure I-18 : Temporal evolution of the measured (in blue) and estimated (in red lines) GPP, R<sub>eco</sub> and NEE, top, and cumulated NEE, bottom, for AUR2012. The yellow envelopes represent the daily standard deviation of the 10 (30) best simulations (i.e. smaller RMSE GAI error). The vertical dashed lines define the cropping year.
<figure><img src="image_34.png" /></figure>

Figure I-19 : Temporal evolution of the measured (in blue) and estimated (in red or red dashed lines) cumulated NEE for LAM2013. The red/red dashed lines represent the simulations that do/do not account for re-growth and weed events. The yellow envelopes represent the daily standard deviation of the 10 (/30) best simulations (i.e. smaller RMSE GAI error).
<figure><img src="image_35.png" /></figure>

Figure I-20 : Temporal evolution of the measured (in blue) and estimated (in red lines) GPP, R<sub>eco</sub> and NEE, top, and cumulated NEE, bottom, for AUR2014. The yellow envelopes represent the daily standard deviation of the 10 (/30) best simulations (i.e. smaller RMSE GAI error). The vertical dashed lines define the cropping year.
4.8.3. Comparison of the performance of several
cropland ecosystem models

We used the study from Wattenbach et al. (2010) to compare the statistical performances of SAFY-CO₂ (with and without accounting for post-harvest vegetative events) with those of other cropland ecosystem models.

The following models were compared:

1 – CERES-EGC (Gabrielle et al., 2006; Lehuger et al., 2009)

2 – DNDC (Li et al., 2005, 1994, 1992)

3 – SPA (Williams et al., 1996)

4 – ORCHIDEE-STICS (de Noblet-Ducoudré et al., 2004; Gervois et al., 2008b)

5 – SAFY-CO₂ without taking post-harvest vegetation events into account

6 – SAFY-CO₂ including post-harvest vegetation events simulation

In their study, Wattenbach et al. (2010) compared the ability of those 4 models (1 to 4) to reproduce CO₂ fluxes (GPP, R<sub>eco</sub> and NEE) and evapotranspiration (ETa) for winter wheat and maize in 5 European sites, including Auradé. The models were evaluated for the Auradé 2006 in situ data, which we also simulated in our study (Fig. 1-21).

Overall, the version of SAFY-CO₂ accounting for post-harvest vegetative events achieved better performances than the other models (Table I-8). The coefficient of determination was between 0.862 (ORCHIDEE) and 0.993 (SPA & SAFY-CO₂ re-growth), while the RMSE [gC.m⁻²] was between 123.22 (ORCHIDEE) and 22.71 (SAFY-CO₂ re-growth) (Table. The 4 models (other than SAFY-CO₂) showed difficulties in correctly reproducing phenology. ORCHIDEE-STICS, SPA and CERES-EGC overestimated the amplitude of crop development, leading to an overestimation of the maximum biomass. The authors also observed a phase shift in the DNDC simulations that resulted in crop development that occurred too early in the season. These results show the power of assimilating GAI derived from remote sensing to constrain crop modelling. By reproducing satellite observations, the model can reproduce well the crop phenology, biomass, yield and CO₂ fluxes well without management data.

<figure><img src="image_36.png" /></figure>

**Figure I-21** : Cumulated NEE simulated by 5 different models (CERES-EGC, DNDC, SPA, ORCHIDEE-STICS and SAFY-CO2 with and without accounting for post-harvest vegetation events) at the Auradé site during the 2006 cultural year (winter wheat).
<table><caption>Table I-8 : Statistic performances (R², RMSE and RRMSE) of 5 cropland ecosystem models : CERES-EGC, DNDC, SPA, ORCHIDEE-STICS and SAFY-CO2 (with and without accounting for post-harvest vegetation events)</caption><thead><tr><th></th><th>CERES</th><th>DNDC</th><th>SPA</th><th>ORCHIDEE</th><th>SAFY-CO2</th><th>SAFY-CO2 re-growth</th></tr></thead><tbody><tr><th>R<sup>2</sup></th><td>0,970</td><td>0,971</td><td>0,993</td><td>0,862</td><td>0,968</td><td>0,995</td></tr><tr><th>RMSE [gC/m<sup>2</sup>]</th><td>45,16</td><td>49,04</td><td>46,49</td><td>123,22</td><td>53,14</td><td>33,62</td></tr><tr><th>RRMSE [%]</th><td>-15,6</td><td>-16,9</td><td>-16,0</td><td>-42,5</td><td>-18,3</td><td>-11,6</td></tr></tbody></table>

Note that the best estimation of cumulated NEE was achieved by the version of SAFY-CO₂ that accounts for post-harvest vegetative events. This was the only model able to reproduce post-harvest vegetative events thanks to the use of remotely sensed GAI.

♠ Ma contribution

Cette étude avait été entamée au cours de la thèse d'Amanda Veloso (2014) mais elle n'avait pas été publiée. Le modèle avait été mis en place et la culture de blé avait été simulée. En reprenant ce travail, j'ai recodé le modèle pour m'assurer qu'il n'y avait pas d'erreur dans les formalismes ou dans l'architecture du modèle et cela m'a aussi permis de me l'approprier. J'ai ainsi corrigé quelques erreurs mineures dans les formalismes et j'ai également récupéré toutes les données nécessaires à l'étude (GAI, satellite, sites flux Européen, campagne, etc.) pour les mettre en forme. J'ai ensuite réalisé une étude de sensibilité me permettant de fixer avec une meilleure précision les bornes des paramètres à optimiser. J'ai ensuite obtenu tous les résultats présents dans la publication à l'exception des sorties des autres modèles présentées dans la section I.4.8.3, directement extraites d'une autre étude. Les objectifs, analyses et enjeux de l'étude ont été discuté avec mes encadrants. J'ai écrit l'article, je l'ai soumis et j'ai répondu aux relecteurs.

✓ Ce qu'il faut retenir

♠ Il n'existe pas d'estimations précises du bilan de C des terres arables à large échelle.

♠ Les modèles de culture existants sont mal adaptés à l'estimation des bilans de C des parcelles agricoles sur de grands territoires car ils nécessitent de nombreuses données d'itinéraires techniques.

♠ Le modèle de culture SAFY-CO₂, piloté par des produits satellites, a été créé à partir de SAFY et permet l'estimation du flux net de CO₂ et de ses composantes.

♠ SAFY-CO₂ a été validé pour la culture de blé et présente des performances similaires à d'autres modèles de culture largement utilisés.

♠ Le modèle permet également de simuler la végétation pouvant se développer pendant l'interculture (repousses, adventices, CI)
♠ Si les données concernant l'application d'amendements organiques et l'export des pailles étaient disponibles à la parcelle sur de larges territoires, il serait possible d'estimer leurs bilans de C.

♠ Il est nécessaire de calibrer et valider d'autres cultures avec le modèle SAFY-CO₂ afin de pouvoir simuler des rotations de cultures et estimer les bilans de carbone associés à ces rotations ou à telle ou telle pratique.

## 5. Conclusion de l'article en français

Dans cette étude, nous avons démontré le potentiel de l'assimilation de données satellites à hautes résolutions spatiale et temporelle dans un modèle de culture simple (SAFY-CO₂) pour estimer avec précision les principales composantes du bilan annuel de carbone des parcelles agricoles (i.e. le flux net annuel de CO₂ et le rendement). Cette approche de modélisation est prometteuse puisqu'elle nécessite très peu de paramètres et aucune donnée de pratique pour estimer les productions agricoles et les flux nets de CO₂. Elle devrait être considérée comme un premier pas pour combler le vide qui existe pour obtenir une représentation spatiale explicite des composantes principales du bilan de carbone des parcelles agricoles à l'échelle régionale et pour des rotations de cultures. En effet, la limitation principale concerne les zones de polycultures élevage puisque dans ces cas, l'estimation du bilan de carbone nécessite des données sur i) l'application d'engrais organique et ii) la fraction de paille exportée lors de la récolte, qui ne peuvent pas être estimés par télédétection actuellement. La seconde limitation est celle de la représentation simplifiée de la respiration hétérotrophe estimée par une approche Q₁₀ qui ne permet pas une estimation précise du bilan de carbone à une échelle pluriannuelle. Pour ce faire, le couplage de SAFY-CO₂ avec un modèle de sol devra être réalisé. Une autre limitation de l'approche concerne la disponibilité des images satellites qui pilotent le modèle. Cependant, et grâce aux récentes missions satellites HSTR (Sentinel 2, Landsat 7/8), ce genre d'approche pourrait être généralisée, plus robuste et précise. Les satellites radar à synthèse d'ouverture (e.g. Sentinel 1) pourraient aussi être utilisés en cas de ciel nuageux (Ameline et al., 2018; Revill et al., 2013a). De plus, nos résultats montrent que les performances du modèle pour estimer les flux nets de CO₂ et donc les bilans de carbone sont nettement améliorées une fois les évènements de repousses ou d'adventices post-récoltes pris en compte. Ces évènements, tout comme la présence de culture intermédiaires dans les rotations de culture, sont rarement, voire jamais pris en compte dans les simulations de flux de CO₂ à l'échelle régionale ou globale, bien qu'ils aient un impact significatif sur le bilan carbone des parcelles agricoles.

Dans un objectif d'application à grande échelle, notre approche pourrait être renforcée (validée dans des conditions pédoclimatiques plus variées) et étendue à d'autres cultures en utilisant les données de réseaux de mesures internationaux (e.g. ICOS, FLUXNET) et des missions satellites HSTR récentes.
# 6. Extension à d'autres cultures : cas du maïs

Le choix de la culture de maïs comme culture à calibrer après celle du blé a été motivé en premier lieu par la superficie qu'elle représente et par le fait que cette culture était déjà calibrée dans SAFY (Claverie et al., 2012a). La culture de maïs représente 16.5% des surfaces agricoles en France (d'après statistiques Agreste 'Grandes cultures') ce qui en fait la 3ème culture la plus cultivée après le blé et l'orge en termes de surface.

La simulation d'une nouvelle culture dans SAFY-CO₂, comme dans la plupart des modèles de cultures nécessite un jeu de paramètre spécifique. Les paramètres propres au maïs ont été, pour la majorité, extraits de la littérature. Notons ici qu'il n'existe pas de valeur unique pour ces paramètres et que souvent ils sont assez variables d'une étude à l'autre. On cherchera dans ces cas, les études les plus génériques possibles afin que les paramètres soient représentatifs du plus de configurations possibles (conditions pédoclimatiques, rotation, etc.). Pour les paramètres ne pouvant pas être déduits de la littérature, il a fallu les estimer à partir de mesures in situ. Les sections suivantes détaillent ces travaux.

## 6.1. Calibration du paramètre ELUEb

<figure><img src="image_37.png" /></figure>

**Figure I-22 :** ELUE estimée à partir des données météorologiques et de GPP acquises à Lamasquère en fonction du ratio rayonnement diffus sur rayonnement direct. Les paramètres de la régression sont indiqués dans l'encadré de gauche.

Comme présenté précédemment, le paramètre ELUEb permet de prendre en compte l'effet du rayonnement diffus sur l'efficience d'utilisation de la lumière pour la photosynthèse (voir section I.4.8.1). Ce paramètre n'existe pas dans la littérature puisque l'équation dans laquelle il intervient a été spécialement conçue pour ces travaux. Pour la culture de blé, ce paramètre a été déterminé grâce aux données de plusieurs sites Européen de mesure de flux de CO₂ (réseau ICOS). En revanche pour la culture de maïs seules les données de Lamasquère, où du maïs est cultivé à peu près une année sur
deux, ont été utilisées. Des données d'autres sites ICOS ont toutefois été collectées mais toutes les données nécessaires à l'estimation du paramètre ELUEb (LAI, rayonnement diffus/direct, température, etc.) n'ont pas encore été obtenues. Dès que ces données seront disponibles, ce paramètre sera ajusté en fonction.

A partir de la formule de la GPP [eq. I-1] on peut en déduire la ELUE :

$$ ELUE = \frac{GPP}{R_g * \varepsilon_c * f_{APAR} * f_T(T_a) * sR10} \qquad I-12 $$

Une fois la ELUE estimée à partir des données terrain, le paramètre ELUEb peut être estimé grâce à une relation exponentielle. Un filtre a été appliqué sur les données correspondantes à la période de sénescence car au cours de cette période la plante n'est pas à 100% de ses capacités photosynthétiques. Pour exprimer la ELUE en fonction du ratio rayonnement diffus sur rayonnement direct, nous avons utilisé le PPFD (pour 'photosynthetic photon flux density') et le PPFDd (diffus), tous deux mesurés in situ avec un capteur BF4 (DeltaT, Cambridge, UK).

La Figure I-22 présente les valeurs de ELUE estimées à partir des mesures in situ ainsi que l'équation exponentielle permettant de déduire ELUEb. La valeur déterminée pour le maïs est de 1.23.

## 6.2. Définition du ratio 'root-to-shoot'

<figure><img src="image_38.png" /><figcaption>Figure I-23 : Evolution temporelle 1) du GAI – simulé par SAFY-CO₂ (trait vert), issu des données satellites (point vert plein) et des données in situ destructives (cercle vert) – et 2) du ratio biomasse aérienne sur biomasse racinaire – dérivé de Amos and Walter (2006) (trait marron), dérivé de Yu et al., 2015 (trait pointillé marron) et observé (données in situ, point marron).</figcaption></figure>

Le formalisme utilisé jusqu'ici pour séparer la biomasse totale en parties aériennes et racinaires a été définit pour la culture de blé à partir de Baret et al., (1992). Il a donc fallu trouver une autre équation pour estimer ce ratio, noté RtS, pour le maïs. Il nous est impossible de déterminer une équation à partir de mesures terrains puisque, pour cela, il faudrait un jeu de données conséquent de biomasses aériennes et racinaires sur maïs, que nous n'avons pas. Des recherches bibliographiques ont
permis de trouver 2 équations de RtS pour maïs (Amos et Walters, 2006; Yu et al., 2015). La Figure I-23 présente la dynamique de végétation du maïs en 2014 ainsi que le RtS déduits des 2 études suscitées et de quelques mesures in situ de biomasses aériennes et racinaires issues de la thèse d'Enguerrand Burel (2018) sur le site de Lamasquère.

La formulation de Amos et Walters (2006) a été retenue car elle offrait une meilleure corrélation avec les observations (issues des méta-analyses). Aussi, la comparaison des RtS estimés à partir de Amos et Walter (2006) et des données in-situ offre des performances satisfaisantes (R² = 0.94). Il a donc été décidé de conserver cette dernière dans le modèle [eq. I-13].

$$RtS = 0.15 + 0.53 * \exp^{-0.03*DAS} \qquad I-13$$

Où DAS représente le nombre de jour après le semis ('day after sowing').

## 6.3. Indentification d'un problème concernant les données de flux sur maïs

Malgré les inversions forcées sur la biomasse et les flux de CO₂, le modèle n'est pas parvenu à représenter correctement ces variables. Quand la biomasse était bien estimée par le modèle, les flux de CO₂ étaient surestimés. Inversement, quand les flux étaient correctement estimés, la biomasse était sous-estimée. Plusieurs tests ont été réalisés afin de corriger ces divergences. Le paramètre de surface spécifique foliaire (SLA pour 'specific leaf area') a été modifié comme cela avait été fait dans Battude et al., (2016). Cette modification permet de représenter l'épaississement des feuilles et correspond au moment où l'assimilation surpasse la translocation. Cette modification n'a pas suffi à corriger les divergences entre modèle et observations. Ensuite le RtS et les paramètres de respiration ont été modifiés afin que le modèle produise plus de biomasse pour des flux de CO₂ identiques. Aucun de ces tests n'a permis d'obtenir des sorties de modèle correspondantes aux observations, tout en restant réaliste d'un point de vue agronomique et physiologique.

Suite à cette batterie de tests, il a été conclu qu'un problème pouvait exister au niveau des mesures de flux de CO₂. Cette hypothèse a été renforcée par l'analyse du dispositif de mesure. Ainsi, au moment de l'installation de la tour à flux en 2004, le constructeur de l'analyseur de gaz (LICOR-7500, LiCor, Lincoln, NE, USA, permettant l'estimation des flux de CO₂ et d'eau) indiquait qu'une distance minimale de 1,5 m entre le haut de la canopée et le capteur devait être respectée. Ultérieurement (i.e. en 2013), une recommandation préconisait une distance de 2m. Or cet instrument a été initialement installé à Lamasquère à une hauteur de 3.6 m, et la hauteur maximale de la végétation sur toutes les années culturales en maïs jusqu'en 2015 était comprise entre 2.42 et 2.99 m. La distance minimale de 2 m au développement maximal du maïs n'a donc jamais été respectée. Il est donc possible que les flux soit correctement acquis en début de saison (lorsque la végétation est encore basse), et que la mesure soit détériorée le reste de la saison. La méthode des fluctuations turbulentes, comme son nom l'indique, mesure des flux turbulents souvent représentés par des tourbillons ('eddies' en anglais). Lorsque le capteur est trop près de la canopée, les plus grands tourbillons ne sont plus mesurés et une partie de l'information est perdue.

Afin de vérifier ces hypothèses, les valeurs de biomasse du blé et du maïs ont été estimées à partir des flux de CO₂ mesurés au niveau de la tour. En soustrayant la respiration hétérotrophe (estimé à partir de la température du sol, [eq. I-9]) à la NEE, on obtient directement la NPP. Puis, en divisant la NPP par le contenu en carbone de la plante, on obtient la biomasse totale sèche. Par ailleurs, on dispose également de données de biomasses aériennes destructives. Théoriquement, la biomasse
(totale) dérivée des flux de CO₂ devrait être supérieure à la biomasse (aérienne) mesurée. Sur la Figure I-24, les 2 graphes de gauche correspondent aux données acquises sur blé, tandis que les graphes de droite correspondent aux données acquises sur maïs. On constate effectivement que pour la culture de maïs, les premières mesures de biomasses destructives (pendant les 50 jours suivant la levée) sont cohérentes avec les valeurs de biomasse dérivées des flux (i.e. légèrement inférieures car les données destructives correspondent aux parties aériennes seulement). En revanche, une fois le maïs développé, on remarque un décrochage qui indique un problème sur la mesure des flux, et renforce l'hypothèse d'une distance canopée/capteur trop faible.

<figure><img src="image_39.png" /></figure>

Figure I-24 : Haut : Scatter plots de la biomasse totale estimée à partir des flux de CO₂ en abscisses et de la biomasse aérienne obtenue par mesures destructives en ordonnées. Bas : Biomasse totale dérivée des flux de CO₂ en trait bleu et biomasses aériennes destructives avec incertitudes. Les 2 graphes de gauche concernent la culture de blé et les 2 graphes de droite concernent la culture de maïs.

Ces résultats indiquent que la validation des simulations pour la culture de maïs ne pourra pas être réalisée à partir des mesures de flux. En revanche, il a tout de même été décidé de valider une première étape de paramétrisation pour la culture de maïs (en attendant de pouvoir disposer de données de flux plus récentes ou sur d'autres sites) mais seulement à partir des données de GAI, de rendement et de biomasse destructive. La section suivante présente ces travaux.

## 6.4. Validation de la culture de maïs
Il a été remarqué lors de la calibration de la culture de blé, une forte sensibilité du modèle au paramètre ELUEa. À cause des compensations entre paramètres du modèle, il est possible d'obtenir des dynamiques de GAI comparables pour des valeurs de ELUEa très différentes. Ce problème d'équifinalité avait déjà été constaté dans la version initiale du modèle (Duchemin et al., 2008a). Or le paramètre ELUEa (qui définit la ELUE) impacte la GPP au premier ordre. Une variation de ce paramètre entraîne une modification significative de la GPP et donc de la biomasse. Il est donc primordial de borner ce paramètre correctement, afin de réduire la dispersion dans les sorties du modèle.

Afin de déterminer les bornes de ELUEa, le modèle a été inversé sur le GAI ainsi que sur la biomasse, tout en laissant de la liberté au paramètre ELUEa (borné sur l'intervalle $[1 – 2.5]$). La méthode utilisée (simplex contraint (Lagarias et al., 1998)) cherche à minimiser une fonction coût, par modification des valeurs des paramètres à optimiser. Cette méthode itérative offre l'avantage de pouvoir résoudre rapidement des problèmes d'optimisation linéaire à plus de trois variables. L'inconvénient est qu'elle peut aboutir à des minimums locaux, et par conséquent échouer à trouver le minimum global. Pour pallier à ce problème, le simplex est lancé une trentaine de fois, en choisissant des points de départs de manière aléatoire (contraints par les valeurs maximum et minium de chaque paramètre à optimiser) à chaque optimisation.

Trente inversions ont donc été lancées sur le GAI satellite et la biomasse, pour les 5 années culturales de maïs. Les résultats de ces inversions sont présentés sur la Figure I-25. Sans surprise le GAI et la biomasse sont très bien représentés par le modèle (puisque ces deux variables sont prises en comptes dans la fonction coût). La distribution de ELUEa, issue des inversions sur les 5 années culturales de maïs est aussi représentée Figure I-25, en considérant la meilleure optimisation (de chaque année) ou les 10 meilleures optimisations (sur 30). Cette distribution nous permet de définir les bornes de la ELUEa pour le maïs à $[1.6 – 1.9] gC.MJ⁻¹$. Ces bornes restreintes devraient permettre de limiter la dispersion dans les sorties et ainsi de pallier aux problèmes d'équifinalité du modèle.

<figure><img src="image_40.png" /><figcaption>Figure I-25 : Résultats d'inversion de SAFY-CO₂ sur le GAI satellite et sur les biomasses in-situ pour les années en maïs à FR-Lam. Scatterplot du GAI (gauche), de la biomasse (centre) et distribution du paramètre optimisé ELUEa issu de la meilleure optimisation et des 10 premières optimisations (sur 30).</figcaption></figure>

L'étape suivante a consisté à valider ces nouvelles bornes en s'appuyant sur des données de biomasses et de rendement spatialisées. Voir section 4.3 pour le détail de ces mesures.

Le modèle a donc été appliqué sur 14 ESU (6 en 2014 et 8 en 2015), situées sur différentes parcelles dans la zone de l'OSR, avec les nouvelles bornes de ELUEa. La Figure I-26 présente les comparaisons entre les biomasses aériennes et les rendements estimés par le modèle face aux
observations. Seules les mesures faites en 2014 disposent d'une incertitude. Le modèle a reproduit les
biomasses et les rendements avec précision (rRMSE = 21 et 22% pour les biomasses et les
rendements, respectivement). On notera ici que le modèle n'est pas parvenu à reproduire les plus
faibles biomasses (et rendement) qu'il suretime, et les plus fortes biomasses, qu'il sous-estime. La
sous-estimation des fortes valeurs de biomasse a déjà été observée sur la culture de blé. En effet, les
plus fortes valeurs de biomasses finales de blé étaient mal reproduites par le modèle, ce biais
s'expliquait principalement par la saturation du signal satellite pour les fortes valeurs de GAI. Ainsi,
au-delà d'un GAI > 3.5 m².m⁻², les couverts très développés sont mal représentés par le modèle.
Concernant les plus faibles valeurs de biomasse, la surestimation du modèle est attribuée à son
comportement général qui ne prend pas en compte certains facteurs. Les développements faibles sont
souvent dus à des problèmes de croissance, liés à un ravageur, à des évènements météorologiques
extrêmes ou encore à des maladies. Ces facteurs ne sont pas pris en compte dans notre approche et ne
peuvent donc pas être représentés. Ceci peut expliquer la surestimation des productions par le modèle
pour de faibles développements bien qu'il soit ici difficile d'affirmer que c'est le cas.

<figure><img src="image_41.png" /><figcaption>Figure I-26 : Biomasse aérienne (gauche) et rendement (droite) de maïs estimés par SAFY-CO₂ face aux observations acquises sur 14 ESUs (6 en 2014 et 8 en 2015).</figcaption></figure>

Cette validation limitée à des données de biomasse et de rendement permet toutefois de
proposer un premier jeu de paramètres pour la culture de maïs (voir Tableau I-9). L'objectif premier
de SAFY-CO₂ étant d'estimer les bilans de C sur des parcelles agricoles ainsi que leurs composantes,
des travaux de validation sur flux de CO₂ (idéalement en multisites) seront donc nécessaires. Ce travail
pourra s'appuyer sur les données du réseau d'observation ICOS dont certaines ont commencé à être
collectées pendant la thèse.

Une fois les données de flux de CO₂ exploitables, il faudra en premier lieu redéfinir le
paramètre ELUEb. En effet, il a été calibré à l'aide de données de flux considérées comme
partiellement erronées, et uniquement sur le site de Lamasquère. Il faudra ensuite valider la nouvelle
paramétrisation sur un large jeu de donnée de biomasse, de rendement et de flux de CO₂. Ce travail
sera réalisé en 2021 dans le cadre des projets Naturellement Popcorn
(https://www.popcorn.fr/actualites/consortium-sur-le-projet-naturellement-popcorn) et Soil Carbon
Farming.
<table><caption>Tableau I-9 : Paramètres du modèle SAFY-CO₂ spécifiques à la culture de maïs. Les bornes des paramètres à optimiser sont indiquées ainsi que les sources des paramètres fixes.</caption><thead><tr><th>Description</th><th>Notation</th><th>Unité</th><th>Valeur / Bornes</th><th>Méthode</th><th>Source</th></tr></thead><tbody><tr><td>Climatic efficiency</td><td>e<sub>c</sub></td><td>-</td><td>0.48</td><td>Literature</td><td>Varlet-Grancher (1982)</td></tr><tr><td>Light-interception coefficient</td><td>K<sub>EXT</sub></td><td>-</td><td>0.63</td><td>Literature</td><td>Claverie (2012)</td></tr><tr><td>Minimal temperature for growth</td><td>T<sub>min</sub></td><td>°C</td><td>8</td><td>Literature</td><td>Ritchie (1991)</td></tr><tr><td>Maximal temperature for growth</td><td>T<sub>max</sub></td><td>°C</td><td>45</td><td>Literature</td><td>Ritchie (1991)</td></tr><tr><td>Optimal temperature for growth</td><td>T<sub>opt</sub></td><td>°C</td><td>30</td><td>Literature</td><td>Ritchie (1991)</td></tr><tr><td>Polynomial degree</td><td>b</td><td>-</td><td>2</td><td>Literature</td><td>Duchemin (2008)</td></tr><tr><td>Corrective factor over GPP during senescence</td><td>C<sub>s</sub></td><td>-</td><td>1</td><td>In situ data</td><td>-</td></tr><tr><td>Effective light-use efficiency parameter A</td><td>ELUE<sub>a</sub></td><td>gC.MJ<sup>-1</sup></td><td>[1.6-1.9]</td><td>Calibration</td><td>-</td></tr><tr><td>Effective light-use efficiency parameter B</td><td>ELUE<sub>b</sub></td><td>-</td><td>1.23</td><td>In situ data</td><td>In situ data - FR-Lam</td></tr><tr><td>Maintenance respiration parameter: Q<sub>10</sub></td><td>Q<sub>10</sub></td><td>-</td><td>2</td><td>Literature</td><td>Amthor (2000)</td></tr><tr><td>Maintenance respiration parameter: R<sub>10</sub></td><td>R<sub>10</sub></td><td>gC/gDM</td><td>0.0015</td><td>literature</td><td>Yang (2004)</td></tr><tr><td>Growth respiration conversion efficiency parameter</td><td>Y<sub>G</sub></td><td>-</td><td>0.75</td><td>Literature</td><td>Johnson (1983)</td></tr><tr><td>Root fraction parameters</td><td>f<sub>1</sub> - f<sub>2</sub></td><td>-</td><td>0.15 / 0.53</td><td>Literature</td><td>Amos (2006)</td></tr><tr><td>Carbon content coefficient</td><td>C<sub>veg</sub></td><td>gC/gveg</td><td>0.46</td><td>Literature</td><td>Béziat (2009)</td></tr><tr><td>Harvest index</td><td>HI</td><td>-</td><td>0.5</td><td>Literature</td><td>Cabelguenne (1999)</td></tr><tr><td>Heterotrophic respiration parameter: R<sub>href</sub></td><td>R<sub>href</sub></td><td>gC.m<sup>-2</sup>.d<sup>-1</sup></td><td>0.34</td><td>Literature</td><td>Suleau (2011)</td></tr><tr><td>Heterotrophic respiration parameter: Q<sub>10</sub></td><td>Q<sub>10</sub></td><td>-</td><td>2.3</td><td>Literature</td><td>Suleau (2011)</td></tr><tr><td>Conversion factor of Ta into Ts</td><td>t</td><td>-</td><td>1.07</td><td>In situ data</td><td>-</td></tr><tr><td>Specific leaf area</td><td>SLA</td><td>m<sup>2</sup>.g<sup>-1</sup></td><td>[0.001 - 0.04]</td><td>Calibration</td><td>-</td></tr><tr><td>Partition-to-leaf function parameter A</td><td>Pl<sub>A</sub></td><td>-</td><td>[0.01-0.5]</td><td>Calibration</td><td>-</td></tr><tr><td>Partition-to-leaf function parameter B</td><td>Pl<sub>B</sub></td><td>-</td><td>[0.0001-0.01]</td><td>Calibration</td><td>-</td></tr><tr><td>Sum of temperature for senescence</td><td>Sen<sub>A</sub></td><td>°C</td><td>[800-1600]</td><td>Calibration</td><td>-</td></tr><tr><td>Rate of senescence</td><td>Sen<sub>B</sub></td><td>°C.day<sup>-1</sup></td><td>[10<sup>3</sup>-3.10<sup>4</sup>]</td><td>Calibration</td><td>-</td></tr><tr><td>Day of plant emergence</td><td>D<sub>0</sub></td><td>day</td><td>[1er Mai-20 Juil]</td><td>Calibration</td><td>-</td></tr></tbody></table>

✓ Ce qu'il faut retenir

▲ Un problème d'acquisition concernant les flux de CO₂ a empêché une validation complète de la culture de maïs.

▲ Une paramétrisation basée sur des données de la littérature et des mesures effectuées sur le site de Lamasquère est tout de même proposée. Elle a été partiellement validée grâce à des données de biomasses et de rendements en multisites.

▲ Une validation reste à réaliser sur des flux de CO₂ pour pouvoir introduire le maïs dans les rotations de culture et pouvoir simuler les bilans associés.
# 7. Extension à d'autres cultures : cas du tournesol

De la même manière que pour la culture maïs, le choix de paramétrer/valider la culture de tournesol a été réalisé de manière pragmatique. Au niveau national, c'est la culture de colza qui est cultivée en majorité après le blé et le maïs, avec 16.1% des surfaces occupées en moyenne sur la période 2014-2018 (d'après les données Agreste grandes cultures. Surfaces 2019) suivi de la culture de tournesol avec 6.4% des surfaces. En revanche, le tournesol est une culture très présente dans notre zone d'étude. Le tournesol représentait en effet 34 et 31% des surfaces en grandes cultures pour la Haute-Garonne et le Gers, respectivement, en 2019 tandis que le colza représentait seulement 7 et 5% respectivement. De plus, la rotation blé/tournesol (ou tournesol/blé) est la rotation la plus observée dans la zone d'étude. Enfin, le tournesol fait aussi partie des cultures ayant déjà été simulées avec SAFY (Claverie et al., 2012a).

Les paramètres propres au tournesol ont été dérivés de la littérature, et certains formalismes ont dû être modifiés. Les prochaines parties détaillent ces adaptations à la culture, et une validation est ensuite proposée. Ce travail a donné lieu à la rédaction d'un proceeding, ainsi qu'une présentation au M2GARSS (https://2020.m2garss.org).

## 7.1. Calibration du paramètre ELUEb

Aucune donnée autre que celles acquises sur le site d'Auradé n'a pu être utilisée pour l'estimation du paramètre ELUEb. En effet, des données de GPP, de rayonnement diffus, direct, de GAI, de température etc. sont nécessaires à l'estimation, mais difficilement accessibles sur d'autres sites. En outre, parmi les 2 années culturales de tournesol à Auradé, une seule a pu être utilisée pour l'estimation de ce paramètre, un problème étant survenu sur le capteur de rayonnement diffus en 2016.

<figure><img src="image_42.png" /></figure>

Figure I-27 : ELUE estimée à partir des données météorologiques et de GPP acquises à Auradé en 2007 en fonction du ratio rayonnement diffus sur rayonnement direct. Les paramètres de la régression sont indiqués dans l'encadré de gauche.

La Figure I-27 présente la ELUE déduite des mesures in situ réalisées à Auradé au cours de
l'année 2007. La valeur retenue pour le paramètre ELUEb est de 1.14. Un filtre a été appliqué sur les valeurs de GAI inférieures à 0.5 m².m⁻² car des valeurs aberrantes de ELUE (>8 gC.m⁻²) étaient observées pour les faibles GAI. Comme le nombre de points utilisé pour réaliser cette régression est limité et que cette analyse est cantonnée à un seul site, il sera nécessaire de réaliser une étude plus poussée en multisites et sur plusieurs années de mesures.

## 7.2. Définition du ratio ‘root-to-shoot’

La recherche d'un formalisme pour définir le RtS du tournesol n'a pas été facile. En effet, les recherches bibliographiques n'ont pas été fructueuses et aucune autre formulation n'a pu être directement introduite dans le modèle. Il a donc été décidé de repartir de la formulation de Baret et al., (1992), initialement utilisée pour le blé, et de l'adapter au cas du tournesol. N'ayant aucune donnée spécifique à nos sites pour adapter cette formulation, les données de l'études de Ma et al., (2017) ont été utilisées. Dans cette étude les auteurs ont mesuré la biomasse aérienne et racinaire à plusieurs moments de la croissance du tournesol, permettant de calculer le RtS à trois dates au cours de la saison et sous différentes concentrations d'azote et de salinité. La Figure I-28 présente la dynamique de végétation du tournesol à Auradé en 2007 simulée par SAFY-CO₂. Les données de l'étude de Ma et al., (2017) y sont également représentées, ainsi que le RtS adapté de Baret et al., (1992). Comme pour le paramètre ELUEb, le RtS pourra être définit de manière plus précise une fois des données de biomasse aérienne et racinaire disponibles, ou encore en s'appuyant sur les sorties de simulation d'autres modèles agronomiques (i.e. STICS).

<figure><img src="image_43.png" /><figcaption>Figure I-28 : Evolution temporelle 1) du GAI – simulé par SAFY-CO₂ (trait vert), issu des données satellite (point vert plein) et des données destructives (cercle vert) – et 2) du ratio biomasse aérienne sur biomasse racinaire – issus de Ma et al., 2017 (ombrage marron) et simulé (trait marron) –.</figcaption></figure>
## 7.3. Validation de la culture de tournesol

Afin de valider la paramétrisation définie pour la culture de tournesol, le modèle a été appliqué sur les deux années culturales de tournesol disponibles à Auradé. Il a ensuite été validé sur le GAI, la biomasse ainsi que les flux de CO2. Les données utilisées ici sont détaillées section 4. Les prochaines sections se concentrent sur la validation du modèle.

### 7.3.1. Estimations du GAI et de la biomasse

La Figure II-29 présente la comparaison entre le GAI et la biomasse simulés par SAFY-CO2 et les observations. Les valeurs de GAI observées sont dérivées des satellites FORMOSAT-2 et SPOT2/4/5. Les observations de biomasse sont des mesures destructives réalisées sur la parcelle d'Auradé (voir section 4.1). Dans l'ensemble, les GAI et les biomasses simulés concordent bien avec les observations. Sur l'année 2007, le GAI est estimé avec peu d'erreur (rRMSE = 8.6%), une très bonne corrélation (R2 = 0.97) et presque pas de biais (0.003 m².m⁻²). Le GAI est également très bien estimé sur l'année 2016 (rRMSE = 2.4%, R2 = 0.9997, biais = 0.013 m².m⁻²). Ces statistiques montrent que le modèle interpole le GAI satellite avec précision, comme cela a pu être montré dans de précédentes études (Duchemin et al., 2015a; Fieuzal et al., 2011a). Le modèle estime aussi la biomasse avec une grande précision, avec des coefficients de détermination de 0.99 et 0.97, et des rRMSE 12.3% et 20.0%, pour les années 2007 et 2016 respectivement. Ces performances sont similaires à une assimilation combinée de signaux satellites optique et radar dans le modèle (Fieuzal et al., 2017).

<figure><img src="image_44.png" /><figcaption>Figure I-29 : Scatter plot entre GAI (gauche) et biomasse (droite) observé et simulé avec SAFY-CO2. Les observations de GAI sont dérivées des mesures satellites, tandis que les biomasses observées sont obtenues par mesures destructives sur le site d'Auradé.</figcaption></figure>

Les deux années culturales de tournesol ont été marquées par de faibles biomasses qui s'expliquent par un été froid en 2007, accompagné de faibles précipitations en juillet. En 2016, le mois de juillet a été peu lumineux, tandis que le mois suivant a été très sec par rapport aux normales. Ces conditions climatiques ne sont pas favorables à une forte photosynthèse. La biomasse simulée en 2016 est surestimée de 9.4%, erreur qui pourrait provenir de la surestimation du GAI satellite. En effet, cette année-là, le GAI destructif maximum était de 1.6 m².m⁻² tandis que le GAI satellite maximum était de 3.6 m².m⁻². Cet écart peut aussi être expliqué par la forte hétérogénéité de développement spécifique à la culture de tournesol, qui rend la mesure destructive peu représentative de la parcelle entière.
### 7.3.2. Validation des flux de CO₂

Les performances obtenues sur le flux net de CO₂ et ses composantes sont résumées dans le Tableau I-10. Comme il n'y a pas de photosynthèse simulée avant l'émergence et après la récolte, seule la période de végétation est considérée pour le calcul des performances sur la GPP. Pour la $R_{ECO}$ et la NEE, toute la période de simulation est considérée.

<table><caption>Tableau I-10 : Performances statistiques associées à la comparaison entre composantes du flux net de CO₂ simulées et observées pour les années culturales 2007 et 2016 de tournesol à Auradé.</caption><thead><tr><th></th><th>GPP</th><th>R<sub>eco</sub></th><th>NEE</th></tr></thead><tbody><tr><td>Nb de données</td><td>214</td><td>672</td><td>731</td></tr><tr><td>R<sup>2</sup></td><td>0.90</td><td>0.86</td><td>0.87</td></tr><tr><td>RMSE (gC.m<sup>-2</sup>.d<sup>-1</sup>)</td><td>1.61</td><td>0.85</td><td>0.81</td></tr></tbody></table>

En plus des statistiques, les évolutions temporelles de flux net de CO₂ et de ses composantes sont représentées (Fig. I-30), de manière à analyser les périodes de divergences entre simulation et observation.

Le modèle a reproduit la dynamique de GPP avec une bonne précision ($R^2$ = 0.90) et peu d'erreur (RMSE = 1.61 gC.m⁻².d⁻¹) sur les deux années de simulation. La GPP simulée est légèrement sous-estimée, en particulier pendant la période de sénescence. À la fin de la saison culturale 2007, on observe deux pics de GPP (Fig. I-30) qui sont dus à des périodes de forts rayonnements et à de fortes températures qui n'ont pas été correctement reproduites par le modèle. En 2016, le GAI satellite indique une sénescence précoce par rapport au GAI in situ, ce qui entraîne une sous-estimation de la GPP sur cette période. Néanmoins, la dynamique de GPP est correctement reproduite par le modèle malgré les différences observées entre GAI satellite et in situ (dues aux hétérogénéités de développement comme évoqué précédemment).

La simulation de la $R_{ECO}$ montre des performances légèrement moins bonnes en termes de dynamique ($R^2$ = 0.86) qui s'expliquent par la représentation très simple de la respiration hétérotrophe dans le modèle (formalisme seulement basé sur la température du sol). En effet, bon nombre de phénomènes ne sont pas pris en compte dans notre approche, comme l'effet du contenu en eau du sol ou du contenu en matière organique. De la même manière que pour la GPP, la $R_{ECO}$ est sous-estimée par le modèle à la fin de la période de végétation (Fig. I-30) à cause de la sous-estimation de la respiration autotrophe (dépendante de la GPP) sur cette période. L'erreur associée à l'estimation de la $R_{ECO}$ reste néanmoins faible (RMSE = 0.85 gC.m⁻².d⁻¹) malgré la représentation simple de la respiration hétérotrophe.

La dynamique journalière de la NEE est correctement reproduite par le modèle en termes de corrélation ($R^2$ = 0.87) et d'erreur (RMSE = 0.81 gC.m⁻².d⁻¹). Les principales sources d'erreur dans l'estimation de la NEE proviennent de i) la représentation simple de la respiration hétérotrophe qui ne permet pas une estimation précise des flux pendant les phases de sol nu et de ii) la sous-estimation de la GPP à la fin de la période de végétation en 2007 qui entraîne une sous-estimation de la NEE. Néanmoins, bien qu'aucun stress hydrique ne soit simulé par cette version du modèle, les performances du modèle dans l'estimation des flux de CO₂ sont très bonnes.
<figure><img src="image_45.png" /><figcaption>Figure I-30 : Evolutions temporelles des observations (rouge) et des simulations (bleu) de GPP (haut), de R<sub>ECO</sub> (milieu) et de NEE (bas) pour les années culturales 2007 et 2016 de tournesol à Auradé.</figcaption></figure>

### 7.3.3. Paramétrisation du tournesol

La validation précédente a permis de déterminer un jeu de paramètre pour la culture de tournesol. Ce jeu de paramètres est présenté dans le Tableau I-11.
<table><caption>Tableau I-11 : Paramètres du modèle SAFY-CO₂ spécifiques à la culture de tournesol. Les bornes des paramètres à optimiser sont indiquées ainsi que les sources des paramètres fixes.</caption><thead><tr><th>Description</th><th>Notation</th><th>Unité</th><th>Valeur / Bornes</th><th>Méthode</th><th>Source</th></tr></thead><tbody><tr><td>Climatic efficiency</td><td>e<sub>c</sub></td><td>-</td><td>0.48</td><td>Literature</td><td>Varlet-Grancher (1982)</td></tr><tr><td>Light-interception coefficient</td><td>K<sub>EXT</sub></td><td>-</td><td>0.8</td><td>Literature</td><td>Flénet (1996)<br>Steer (1996)</td></tr><tr><td>Minimal temperature for growth</td><td>T<sub>min</sub></td><td>°C</td><td>8</td><td>Literature</td><td>STICS</td></tr><tr><td>Maximal temperature for growth</td><td>T<sub>max</sub></td><td>°C</td><td>42</td><td>Literature</td><td>STICS</td></tr><tr><td>Optimal temperature for growth</td><td>T<sub>opt</sub></td><td>°C</td><td>25</td><td>Literature</td><td>STICS</td></tr><tr><td>Polynomial degree</td><td>b</td><td>-</td><td>2</td><td>Literature</td><td>Duchemin (2008)</td></tr><tr><td>Corrective factor over GPP during senescence</td><td>C<sub>s</sub></td><td>-</td><td>1</td><td>In situ data</td><td>-</td></tr><tr><td>Effective light-use efficiency parameter A</td><td>ELUE<sub>a</sub></td><td>gC.MJ<sup>-1</sup></td><td>[0.75-0.95]</td><td>Calibration</td><td>-</td></tr><tr><td>Effective light-use efficiency parameter B</td><td>ELUE<sub>b</sub></td><td>-</td><td>1.14</td><td>In situ data</td><td>In situ data - FR-Lam</td></tr><tr><td>Maintenance respiration parameter: Q<sub>10</sub></td><td>Q<sub>10</sub></td><td>-</td><td>2</td><td>Literature</td><td>Amthor (2000)</td></tr><tr><td>Maintenance respiration parameter: R<sub>10</sub></td><td>R<sub>10</sub></td><td>gC/gDM</td><td>0.0036</td><td>literature</td><td>Yang (2004)</td></tr><tr><td>Growth respiration conversion efficiency parameter</td><td>Y<sub>G</sub></td><td>-</td><td>0.76</td><td>Literature</td><td>Johnson (1983)</td></tr><tr><td>Root fraction parameters</td><td>f<sub>0</sub> - f<sub>inf</sub> - c</td><td>-</td><td>0.3 / 0.1 / 1.5</td><td>Literature</td><td>Baret (1992) - Ma (2017)</td></tr><tr><td>Carbon content coefficient</td><td>C<sub>veg</sub></td><td>gC/g<sub>veg</sub></td><td>0.46</td><td>Literature</td><td>Béziat (2009)</td></tr><tr><td>Harvest index</td><td>HI</td><td>-</td><td>0.29</td><td>Literature</td><td>Cabelguenne (1999)</td></tr><tr><td>Heterotrophic respiration parameter: R<sub>href</sub></td><td>R<sub>ref</sub></td><td>gC.m<sup>-2</sup>.d<sup>-1</sup></td><td>0.34</td><td>Literature</td><td>Suleau (2011)</td></tr><tr><td>Heterotrophic respiration parameter: Q<sub>10</sub></td><td>Q<sub>10</sub></td><td>-</td><td>2.3</td><td>Literature</td><td>Suleau (2011)</td></tr><tr><td>Conversion factor of Ta into Ts</td><td>t</td><td>-</td><td>1.07</td><td>In situ data</td><td>-</td></tr><tr><td>Specific leaf area</td><td>SLA</td><td>m<sup>2</sup>.g<sup>-1</sup></td><td>[0.001 - 0.04]</td><td>Calibration</td><td>-</td></tr><tr><td>Partition-to-leaf function parameter A</td><td>Pl<sub>A</sub></td><td>-</td><td>[0.01-0.5]</td><td>Calibration</td><td>-</td></tr><tr><td>Partition-to-leaf function parameter B</td><td>Pl<sub>B</sub></td><td>-</td><td>[0.0001-0.01]</td><td>Calibration</td><td>-</td></tr><tr><td>Sum of temperature for senescence</td><td>Sen<sub>A</sub></td><td>°C</td><td>[200-1800]</td><td>Calibration</td><td>-</td></tr><tr><td>Rate of senescence</td><td>Sen<sub>B</sub></td><td>°C.day<sup>-1</sup></td><td>[10<sup>3</sup>-3.10<sup>4</sup>]</td><td>Calibration</td><td>-</td></tr><tr><td>Day of plant emergence</td><td>D<sub>0</sub></td><td>day</td><td>[4 Avril-20 Juil]</td><td>Calibration</td><td>-</td></tr></tbody></table>
✓ Ce qu'il faut retenir

▲ La culture de tournesol a été introduite dans le modèle SAFY-CO₂ et validée sur des données de biomasse et de flux de CO₂ acquises sur le site d'Auradé.

▲ Certains paramètres nécessiteraient une calibration plus robuste (sur la base de données acquises durant plusieurs années culturales, et sur différents sites) car peu de données étaient disponibles jusqu'à présent.

▲ Comme le tournesol est une culture fréquemment soumise au stress hydrique, le couplage de SAFY-CO₂ avec un module de bilan hydrique devra être testé.

▲ Avec trois grandes espèces paramétrées, le modèle pourra être utilisé pour simuler des rotations de culture et ainsi pouvoir quantifier le bilan de C de certaines rotations

8. Conclusion

Dans cette section du manuscrit, le modèle SAFY-CO₂ a été appliqué avec succès à la culture de blé. Le modèle a été capable de reproduire avec précision les biomasses, rendements et flux de CO₂ (GPP, RECO et NEE). Comme des données concernant l'export des pailles ainsi que l'application d'engrais organique (à Lamasquère) étaient disponibles, des bilans de carbones ont pu être estimés par cette approche. Il a notamment été montré que le modèle était capable d'estimer les bilans annuels de C sur parcelles agricoles en blé pour des années climatiques contrastées avec une précision de 77 gC.m⁻².an⁻¹.

Concernant la culture de tournesol, elle a également été simulée avec succès par SAFY-CO₂. Le modèle a montré de bonnes performances dans l'estimation des biomasses et des flux de CO₂. Toutefois, des analyses supplémentaires devront être réalisées afin de rendre la paramétrisation plus robuste/générique. En effet, le paramètre ELUEb a été calibré avec un nombre très limité d'observation, sur le même site où les validations ont été réalisées. Enfin, le formalisme RtS pourrait également être définit avec plus de données ou avec un formalisme différent et validé pour cette culture. Cette paramétrisation a pu être appliquée pour une analyse des performances du modèle à simuler des biomasses et rendements spatialisés. Cette étude sera présentée dans le chapitre suivant.

Les mesures de flux de CO₂ acquises à Lamasquère sur maïs n'ont pas pu être exploitées pour valider le modèle sur cette culture. En revanche, des données de biomasses et de GAI ont été utilisées et ont permis d'aboutir à une paramétrisation qui a été validée sur des données spatialisées de biomasses et de rendements. La culture de maïs a néanmoins besoin de validations supplémentaires, moyennant des flux de CO₂ mesurés sur d'autres sites et d'autres années. Le paramètre ELUEb aura aussi besoin d'être recalibré car sa valeur actuelle a été obtenue avec les données de CO₂ qui sont considérées comme partiellement erronées.

La paramétrisation du modèle pour la culture de colza faisait partie des objectifs initiaux de la
thèse. En effet, avec les trois cultures suscitées, elles représentent les quatres les plus cultivées en France. Les problèmes observés sur les flux de CO₂ acquis sur maïs ont également été constatés sur une des années de mesure sur colza à Auradé. La paramétrisation de cette culture a donc été repoussée. Néanmoins, un grand jeu de donnée est désormais disponible (grâce au projet Colza Digital et via d'autres sites de flux en Europe) et permettra une validation rapide de cette culture.

Enfin, l'étude réalisée sur la culture de blé a permis de mettre en exergue l'importance de la
prise en compte des repousses/adventices pour une estimation précise des bilans de carbone.
Actuellement les repousses sont simulées avec une paramétrisation générique, sans distinguer la nature
de la végétation présente (i.e. repousse, CI, adventices). Un travail devra être mené pour affiner cette
paramétrisation en fonction de leur nature.
# Chapitre II
# Chapitre II : Intégration d'un module de bilan hydrique à SAFY-CO₂

## 1. Pourquoi simuler les composantes du bilan d'eau

L'estimation des flux d'eau sur les parcelles agricoles et sur de grands territoires a un intérêt à la fois sociétal et environnemental. En 2013, la France a prélevé environ 33.5 milliards de m³ d'eau (Aillery et al., 2018), essentiellement en surface au niveau des lacs et des rivières (83%) mais également sous terre. L'essentiel de cette eau est destiné aux centrales nucléaires et sert à les refroidir (Fig. II-1). La quasi-totalité de cette eau 'de refroidissement', comme celle utilisée pour le maintien des canaux navigables, est restituée au milieu naturel. Les prélèvements d'eau potable, ainsi que ceux liés à l'usage agricole et industriel représentent 11 milliards de m³ en 2013. Le secteur agricole consomme à lui seul 2.7 milliards de m³ dont 37% proviennent des nappes phréatiques. Notons ici que les secteurs agricoles et industriels prélèvent chacun la même quantité d'eau.

<figure><img src="image_46.png" /></figure>

Champ : France métropolitaine.

Source : Agence française de la biodiversité, Banque nationale des prélèvements quantitatifs

en eau (BNPE). Traitements : SDES, 2016

<caption>**Figure II-1 : Prélèvement d'eau en France en 2013 par secteur. Exprimé en milliard de m³**</caption>

Sur les 2.7 milliards de m³, 80% sont utilisés pour l'irrigation, le reste servant à nettoyer les machines agricoles, les bâtiments et à abreuver le bétail. Pourtant, seulement 6% de la surface agricole utile est irriguée (Le service de l'observation et des statistiques (SOeS), 2017). En 2013, l'irrigation représentait donc en moyenne 1700 m³.ha⁻¹ sur les surfaces irriguées. Néanmoins, l'irrigation est très
contrastée entre régions ; certaines n'irriguent pas du tout (région de la Lorraine ou département du Doubs), tandis que onze départements consomment plus de la moitié de l'eau douce utilisée par l'agriculture (Gers, Bouches-Du-Rhône, Landes, Gironde, Loiret, etc.). Le Sud-Ouest de la France est une zone fortement consommatrice d'eau (Fig. II-2).

Concernant la disponibilité de cette ressource, l'augmentation globale des températures, entraînant des épisodes de sécheresse estivale répétés, tend à la raréfier. Ces épisodes de fortes chaleurs et de faibles précipitations, qui sont amenés à s'accentuer à l'avenir selon les modèles climatiques (Brisson and Levrault, 2010) imposent au monde agricole de s'adapter rapidement. L'été 2020, notamment, a été marqué par de fortes chaleurs et de faibles précipitations entraînant d'importantes restrictions sur la consommation en eau. Au 15 Septembre 2020, ce n'était pas moins de 77 départements qui étaient soumis à une restriction (au-delà de la vigilance) et 179 arrêtés étaient en cours (http://propluvia.developpement-durable.gouv.fr).

Dans le même temps, des tensions concernant la gestion de la ressource en eau sont déjà observables sur le territoire, notamment concernant les nouvelles retenues d'eau opposant agriculteurs et associations environnementales. Ces tensions risquent de s'exacerber dans l'avenir et vont forcer le milieu agricole à s'adapter à cette nouvelle situation (Delgoulet, 2014).

<figure><img src="image_47.png" /></figure>

Note : * "la méthodologie et la carte indiquant les noms des sous-bassins figurent en fin de document.
Source : Onema, Banque nationale des prélèvements quantitatifs en eau (BINPE). Traitements : SOeS, 2016

**Figure II-2 : Prélèvements d'eau en France en 2013 par sous-bassins hydrographiques.** Exprimé en million de m³.

Dans ce contexte, l'estimation spatialisée des composantes du bilan d'eau en agriculture via des modèles agronomiques, et donc de la consommation en eau des cultures, permettrait de mieux gérer cette ressource notamment via le calcul d'indicateurs d'efficience de l'utilisation de l'eau (Tallec et al., 2013). En outre, un autre intérêt à intégrer la simulation des composantes du bilan d'eau dans l'approche de simulation présentée dans le chapitre précédent, est que cela pourrait nous permettre d'améliorer les estimations de biomasse, de rendement, de flux de CO₂ et donc les bilans de C. En effet, le stress hydrique, qui peut être plus ou moins marqué en fonction des cultures et des conditions météorologiques, impacte directement les flux de CO₂ (fermeture stomatique). Pour prendre en compte ces effets, il est communément admis qu'il est nécessaire de simuler l'évolution des stocks d'eau dans le sol, eux même pilotés par les flux d'entrée et de sortie (précipitation, drainage, évaporation,
transpiration et flux diffusifs). Cela nécessite donc de représenter les formalismes régissant ces flux et d'avoir des informations sur les capacités de rétention d'eau du sol.

L'ajout d'un module de bilan hydrique au sein de SAFY-CO₂ a donné naissance au modèle SAFYE-CO₂. Il peut être vu comme le couplage entre SAFY-CO₂ (Pique et al., 2020a) et la composante bilan hydrique du modèle SAFY-WB (Duchemin et al., 2015a) qui est adaptée de la méthode FAO-56 (Allen et al., 1998). Le module de bilan hydrique est tout d'abord détaillé dans la section suivante, puis le modèle est ensuite appliqué sur la parcelle d'Auradé et au niveau régional.

## 2. La simulation des flux d'eau

### 2.1. Généralités concernant la méthode FAO-56

Le module de bilan hydrique qui a été introduit dans SAFY-CO₂ est donc adapté du modèle FAO-56 proposé par la FAO (Allen et al., 1998). Ce modèle prend en compte la demande évaporative de l'atmosphère via l'évapotranspiration de référence (ET₀). ET₀ correspond, selon la FAO, au « taux d'évapotranspiration d'un tapis gazonné vert de hauteur uniforme comprise entre 8 et 15 cm, en pleine croissance, ombrageant totalement le sol et ne manquant pas d'eau ». À partir de ET₀ et d'un coefficient cultural Kc (Fig. II-3), il est alors possible d'estimer l'évapotranspiration d'une culture sous conditions standards (ETc), c'est-à-dire sans aucun stress hydrique. Comme ces conditions sont rarement rencontrées dans la nature, la FAO propose le calcul d'une évapotranspiration ajustée (ETcadj) via un coefficient de stress hydrique Ks, permettant de prendre en compte le stress hydrique subit par la plante. Le module hydrique de SAFYE-CO₂ est basé sur cette approche de calcul de ETcadj.

<figure><img src="image_48.png" /></figure>

**Figure II-3** : Evapotranspiration de référence ET₀, évapotranspiration des cultures en condition optimale de ressource en eau du sol (ETc) et en condition de stress hydrique (ETcadj). Extrait de Allen et al., (1998).
## 2.2. Description du module de bilan hydrique

Le module de bilan hydrique simule au pas de temps journalier le contenu en eau du sol, ainsi que les flux d'évaporation et de transpiration. Dans cette approche, le sol est divisé en deux ou trois couches horizontales. Une couche de surface d'épaisseur constante, dite « couche évaporative » qui est concernée par l'évaporation et le prélèvement d'eau par les plantes pour leur transpiration. Une couche intermédiaire, qui existe seulement en présence de végétation, d'épaisseur variable en fonction de la croissance racinaire, dite « couche racinaire » dans laquelle les plantes puisent de l'eau pour la transpiration. Une couche profonde d'épaisseur constante qui ne contribue ni à l'évaporation ni à la transpiration.

Pour chacune de ces trois couches, il est possible de définir son contenu en eau maximal (accessible par la plante) à l'aide des paramètres de capacité au champ ($\theta_{FC}$, pour 'field capacity', exprimé en m$^3.m^{-3}$) et de point de flétrissement ($\theta_{WP}$, pour 'wetting point', exprimé en m$^3.m^{-3}$). $\theta_{FC}$, définit comme la teneur maximale en eau disponible, correspond à l'humidité du sol après drainage par gravité. $\theta_{WP}$, est définit comme l'humidité en dessous de laquelle la plante ne peut plus extraire l'eau du sol, autrement dit, quand les forces de rétention de l'eau par le sol deviennent supérieures aux forces de succion de la plante. La différence entre ces deux paramètres ($\theta_{FC} - \theta_{WP}$) donne la réserve utile (RU) qui est donc l'eau disponible pour la transpiration de la plante.

À partir de ces paramètres de sol, on peut définir le contenu en eau maximum (TWC pour 'total water content') pour les trois couches (évaporative : 'evap', racinaire : 'root' et profonde : 'deep'). Les équations [Eq. II-1 à II-3] sont établies à partir des $\theta_{FC}$, $\theta_{WP}$ et des épaisseurs Z respectives de ces trois couches.

$$TWC_{evap} = (\theta_{FC\_evap} - \theta_{WP\_evap}) * Z_{evap} \qquad \text{II-1}$$

$$TWC_{root} = (\theta_{FC\_root} - \theta_{WP\_root}) * Z_{root} \qquad \text{II-2}$$

$$TWC_{deep} = (\theta_{FC\_deep} - \theta_{WP\_deep}) * Z_{deep} \qquad \text{II-3}$$

Le contenu en eau de la couche évaporative (AWC$_{evap}$ pour 'actual water content') est mis à jour quotidiennement [Eq. II-4 et II-5], en fonction du contenu en eau de la veille, de l'eau apportée par les précipitations et/ou par l'irrigation (PE), de l'évaporation (E) et de l'eau prélevée par la plante dans cet horizon pour la transpiration (T). Comme 2 couches sont concernées par la transpiration, le partage de prélèvement entre ces dernières se fait grâce à une fonction (ftr) dépendante du contenu relatif de la couche évaporative par rapport à celui de ces deux couches.

$$AWC_{evap}(i) = AWC_{evap}(i-1) + PE(i) - E(i) - T(i) * ftr(i) \qquad \text{II-4}$$

$$ftr(i) = \frac{AWC_{evap}(i-1)}{AWC_{evap}(i-1) + AWC_{root}(i-1)} \qquad \text{II-5}$$

Les flux diffusifs ne sont pas pris en compte dans le calcul du contenu en eau du sol car il a été observé dans de précédentes études que ce terme n'apportait pas plus de précision aux estimations de l'ETR.
La quantité d'eau arrivant au sol (PE) est calculée [Eq. II-6] en fonction de l'eau apportée par les précipitations (P) et/ou l'irrigation (I), ainsi qu'en fonction de l'eau interceptée (Interc) par la végétation. En effet, avant d'atteindre le sol, l'eau peut s'évaperer. L'eau interceptée va alors dépendre [Eq. II-7] de la couverture de végétation au sol ($F_{COVER}$).

$$PE(i) = \max\{P(i) + I(i) - Interc(i), 0\} \qquad \text{II-6}$$

$$Interc(i) = \min\{(P(i) + I(i)) * F_{COVER}, ET_0\} \qquad \text{II-7}$$

La variable $F_{COVER}$ est calculée quotidiennement [Eq. II-8] en fonction du GAI et de deux coefficients spécifiques à la culture considérée, $K_{COV}$ et $E_{COV}$.

$$F_{COVER}(i) = K_{COV} * (1 - \exp(-E_{COV} * GAI(i))) \qquad \text{II-8}$$

La fonction d'estimation de l'évaporation du sol a été modifiée par rapport aux anciennes versions du modèle (SAFY-WB ; Battude et al., 2017, Duchemin et al., 2015a). Elle est décrite plus en détail dans la section suivante.

La transpiration est définie selon la méthode FAO-56 [Eq. II-9] et dépend de $ET_0$, du coefficient cultural $K_{CB}$ (CB pour 'crop basal') et d'un coefficient de stress ($K_s$).

$$T(i) = K_{CB} * K_s(i) * ET_0(i) \qquad \text{II-9}$$

Le coefficient $K_{CB}$ est définit [Eq. II-10] à partir du GAI et de deux paramètres, $K_{CBmax}$ et $E_{TRP}$.

$$K_{CB}(i) = K_{CBmax} * (1 - \exp(-E_{TRP} * GAI(i))) \qquad \text{II-10}$$

$K_s$ est définit [Eq. II-11] en fonction de l'humidité relative ($H_R$) des couches évaporative et racinaire, et d'un paramètre $H_{CRIT}$, spécifique à chaque culture et définissant l'humidité à partir de laquelle un stress apparaît.

$$K_s(i) = \min\left\{1 - \frac{HRS(i - 1) - H_{CRIT}}{H_{CRIT}}\right\} \qquad \text{II-11}$$

$$HRS(i - 1) = \max\{H_{R_{evap}}(i - 1), H_{R_{root}}(i - 1)\} \qquad \text{II-12}$$

$$H_{R_x} = \frac{AWC_x}{TWC_x} \qquad \text{II-13}$$

Le contenu de la couche racinaire ($AWC_{root}$) est calculé [Eq. II-14] en fonction du contenu en eau de la veille, de l'eau excédentaire de la couche supérieure (Pe pour 'percolation', Eq. II-15), et de l'eau extraite par transpiration.

$$AWC_{root}(i) = AWC_{root}(i - 1) + Pe_{evap}(i) - T(i) * (1 - ftr(i)) \qquad \text{II-14}$$

$$Pe_x(i) = \max\{TWC_x(i) - AWC_x(i), 0\} \qquad \text{II-15}$$

Enfin le contenu en eau de la couche profonde ($AWC_{deep}$) est simplement estimé [Eq. II-16] à partir du contenu de la veille, et de l'eau excédentaire de la couche supérieure (i.e. racinaire ou évaporative en fonction de la présence ou non de végétation).
$$AWC_{deep}(i) = AWC_{deep}(i-1) + Pe_{root}(i) \qquad \text{II-16}$$

## 2.3. Nouvelle formulation de l'évaporation

D'après la méthode FAO-56, l'évaporation est définie comme un produit de l'ET₀ par une fonction de stress liée à l'eau disponible dans le sol [Eq. II-17]. Cette fonction de stress, qui approxime la fonction bilinéaire de Allen et al., (1998) dépend de l'humidité relative de la couche évaporative et d'un coefficient de stress β.

$$Evap = (1 - F_{COVER}) * ET_0 * (1 - (1 - H_{Revap})^{\beta}) \qquad \text{II-17}$$

Il a été remarqué que cette fonction ne représentait pas correctement l'évaporation pendant les périodes de sol nu (où FCOVER = 0). Ce constat qui avait déjà été établi par Soares et al. (1988a) s'explique aisément par la définition de l'ET₀ (voir section II-2.1), qui ne représente pas l'évaporation d'un sol nu mais l'évapotranspiration d'une surface couverte de gazon. Afin de corriger cette équation pour les périodes de sol nu (ex. entre deux cultures), le rapport de l'évapotranspiration mesurée (ETRMES), sur les deux sites flux (Auradé et Lamasquère), sur l'ET₀ a été calculé en fonction de l'humidité du sol pour différentes valeurs d'ET₀. Il a été remarqué que i) l'ETRMES était parfois supérieure à l'ET₀ et que ii) l'ETRMES était dépendante de l'ET₀. Le premier point s'explique par la différence d'évaporation entre une surface en gazon et un sol nu, qui a été corrigée par un facteur KS équivalent à 1.2. La seconde observation peut s'expliquer par le fait que lorsque la demande évaporative est forte, les premiers centimètres de sol vont sécher rapidement et vont créer une 'barrière' aux futures remontées capillaires. Le sol aura donc plus de mal à répondre à la demande évaporative de l'atmosphère quand cette dernière sera forte. Pour prendre en compte ce phénomène, la pente α de la fonction linéaire interpolant le rapport ETRMES/ET0 en fonction de l'humidité relative de la couche évaporative a été calculée pour différentes valeurs de ET0 et pour chaque site instrumenté (Auradé et Lamasquère) en ayant corrigé l'ETR par le facteur KS. Cette même pente, qui théoriquement devrait être égale à 1, a été recalculée pour les deux sites en pondérant la valeur de la pente α par le nombre d'observations. La Figure II-4 résume ces résultats. Ce travail m'a permis de définir une nouvelle équation de l'évaporation [Eq. II-18].

$$Evap = (1 - F_{COVER}) * K_S * (Ev_A * ET_0 + Ev_B) * ET_0 * \left( \frac{SM - \theta_{WP\_evap}}{\theta_{FC\_evap} - \theta_{WP\_evap}} \right) \qquad \text{II-18}$$

Afin de s'assurer de l'apport de cette nouvelle équation dans l'estimation de l'évaporation, cette dernière a été appliquée aux sites instrumentés. La Figure II-5 présente les cumuls d'évaporation dérivés des observations, ceux estimés à partir de l'équation de la FAO [Eq. II-17], et avec la nouvelle équation [Eq. II-18], ainsi que la comparaison entre observations journalières d'ETR et les valeurs obtenues avec les deux formalismes. La nouvelle équation offre des résultats légèrement meilleurs que l'équation de la FAO. L'erreur obtenue sur les données de Lamasquère passe de 0.54 mm.j⁻¹ à 0.49 mm.j⁻¹, et de 0.52 mm.j⁻¹ à 0.44 mm.j⁻¹ à Auradé. Notons également qu'une légère amélioration de la corrélation est constatée. La nouvelle équation de l'évaporation est donc conservée pour la suite de ce travail. L'équation de la FAO est toutefois conservée dans le modèle et peut être choisie par l'utilisateur.
<figure><img src="image_49.png" /></figure>

Figure II-4 : Valeur de la pente α en fonction de ET₀. α correspond au coefficient directeur de la fonction affine interpolant le rapport ETR/ET₀ en fonction de l'humidité relative du sol.

<figure><img src="image_50.png" /></figure>

Figure II-5 : Evapotranspiration cumulée estimée à partir des mesures (en bleu), de l'équation de la FAO (Eq. II-17, en rouge) et de la nouvelle équation d'évaporation (New – Eq. II-18, en orange) sur les périodes de sol nu et scatter plots entre ETR observée et estimée à partir des deux équations suscitées pour les sites de Lamasquère (à gauche) et Auradé (à droite).

Suite à l'ajout de ce module de bilan hydrique dans SAFY-CO₂, une étude a été réalisée sur la culture de tournesol. Le nouveau modèle, SAFYE-CO₂, a été évalué sur le site flux d'Auradé et de manière spatialisée grâce à un jeu de données de biomasses et de rendements acquis sur 2 années dans la zone de l'OSR. Elle a donné lieu à une publication en 2020 dans la revue Remote Sensing. Cet article est présenté ci-après.

## 3. Résumé de l'article en Français

L'augmentation des besoins alimentaires mondiaux dans le contexte de changements climatiques requiert une compréhension précise du fonctionnement des terres agricoles et de leur impact sur les cycles biogéochimiques. Or, bien que les échanges gazeux entre les surfaces agricoles et l'atmosphère soient mesurables sur le terrain, ils sont difficilement quantifiables à l'échelle de la
parcelle et sur de grands territoires à cause de l'hétérogénéité des paysages et des pratiques culturales. Pourtant, évaluer précisément les bilans de carbone et d'eau des parcelles agricoles est essentiel pour promouvoir des pratiques agronomiques durables, et réduire la demande en eau et l'impact climatique des cultures, tout en maintenant des rendements suffisants. À partir de ce constat nous avons développé un modèle de culture, SAFYE-CO₂, qui assimile des produits satellites à hautes résolutions spatiale et temporelle (HRST) pour estimer au pas de temps journalier la biomasse, les flux d'eau et de CO₂, le rendement et le bilan de carbone des cultures à l'échelle parcellaire sur de grands territoires. Cette approche de modélisation a été évaluée pour la culture de tournesol grâce à deux jeux de données de terrain. Dans un premier temps, les sorties du modèle ont été comparées à deux années culturales de mesures acquises sur le site instrumenté ICOS ('integrated carbon observation system') d'Auradé dans le sud-ouest de la France. Le modèle a simulé précisément les flux net de CO₂ (RMSE = 0.97 gC.m⁻².d⁻¹ et R² = 0.83) et d'eau (RMSE = 0.68 mm.d⁻¹ et R² = 0.79). Les performances du modèle ont ensuite été évaluées grâce à des données de biomasse et de rendement collectées sur 80 parcelles situées dans le sud-ouest de la France. Le modèle a été capable d'estimer les dynamiques de biomasse ainsi que les rendements avec précision (RMSE respectivement de 66 et 54 g.m⁻²). Afin d'étudier le potentiel d'application de cette méthode à grande échelle, et sachant que les propriétés de sol sont un facteur important en entrée du modèle, une analyse de sensibilité à deux produits sols existants (GlobalSoilMap et SoilGrid) a été réalisée. Nos résultats montrent que ces produits ne sont pas suffisamment précis pour pouvoir être utilisés à une échelle parcellaire en entrée du modèle, qui nécessite des informations plus précises sur les capacités de rétention d'eau du sol nécessaire à l'estimation des flux d'eau. De plus, nous affirmons que dans notre approche, aucun stress hydrique n'a besoin d'être considéré pour le calcul de la photosynthèse, puisque ce stress est déjà pris en compte implicitement via l'information satellite qui permet de décrire la dynamique et l'intensité de développement de la végétation via le GAI. Cette étude devrait être considérée comme une première étape pour combler le vide qui existe concernant l'estimation des composantes des bilans de carbone à l'échelle parcellaire sur de grands territoire, et pour estimer précisément l'effet de certaines pratiques agricoles, comme les cultures intermédiaires, sur les bilans de carbone et d'eau.
# 4. Combining high-resolution remote sensing products with a crop model to estimate carbon and water budget components: Application to sunflower (published in Remote Sensing)

Gaétan Pique <fnref n="1" />,<fnref n="2" />*, Rémy Fieuzal <fnref n="1" />, Philippe Debaeke <fnref n="3" />, Ahmad Al Bitar <fnref n="1" /> and Tiphaine Tallec <fnref n="1" />,
Eric Ceschia <fnref n="1" />,<fnref n="4" />

## 4.1. Abstract

The global increase in food demand in the context of climate change requires a clear understanding of cropland function and of its impact on biogeochemical cycles. However, although gas exchange between croplands and the atmosphere is measurable in the field, it is difficult to quantify at the plot scale over relatively large areas because of the heterogeneous character of landscapes and differences in crop management. However, assessing accurate carbon and water budgets over croplands is essential to promote sustainable agronomic practices and reduce the water demand and the climatic impacts of croplands while maintaining sufficient yields.

From this perspective, we developed a crop model, SAFYE-CO₂, that assimilates high-spatial- and-temporal-resolution (HSTR) remote sensing products to estimate daily crop biomass, water and CO₂ fluxes, annual yields and carbon budgets at the parcel level over large areas. This modeling approach was evaluated for sunflower against two in situ datasets. First, the model's output was compared to data acquired during two cropping seasons at the Auradé ICOS instrumented site in southwestern France. The model accurately simulated the daily net CO₂ flux (RMSE = 0.97 gC.m⁻<fnref n="2" />.d⁻<fnref n="1" /> and R<fnref n="2" /> = 0.83) and water flux (RMSE = 0.68 mm.d⁻<fnref n="1" /> and R<fnref n="2" /> = 0.79). The model's performance was then evaluated against biomass and yield data collected from 80 plots located in southwestern France. The model was able to satisfactorily estimate biomass dynamics and yield (RMSE = 66 and 54 g.m⁻<fnref n="2" />, respectively). To investigate the potential application of the proposed approach at a large-scale, given that soil properties are important factors affecting the model, a sensitivity analysis of two existing soil products (GlobalSoilMap and SoilGrids) was carried out. Our results show that these products are not sufficiently accurate for inclusion as inputs to the model, which requires more accurate information on soil water retention capacity to assess water fluxes. Additionally, we argue that no water stress should be considered in the crop growth computation since this stress is already present because of remote sensing information in the proposed approach

This study should be considered a first step to fulfill the existing gap in quantifying carbon budgets at the plot scale over large areas and to accurately estimate the effects of management practices such as the use of cover crops or specific crop rotations on cropland C and water budgets.

<footnote n="1">¹ CESBIO, Université de Toulouse, CNES/CNRS/INRAE/IRD/UPS, Toulouse, France; fieuzal@cesbio.cnes.fr (R.F.), albitara@cesbio.cnes.fr (A.A.B.), tallect@cesbio.cnes.fr (T.T.)</footnote>

<footnote n="2">² Agence De l'Environnement et de Maîtrise de l'Energie (ADEME), Angers Cedex 1, France;</footnote>

<footnote n="3">³ INRAE, UMR 1248 AGIR, Castanet-Tolosan, France; philippe.debaeke@inra.fr</footnote>

<footnote n="4">⁴ INRAE, USC 1439 CESBIO, Toulouse, France; ceschiae@cesbio.cnes.fr</footnote>

<footnote>* Correspondence: piqueg@cesbio.cnes.fr</footnote>
**Keywords:** crop modeling; carbon budget; water budget; remote sensing; soil properties

## 4.2. Introduction

Sunflower (*Helianthus annuus* L.) was introduced by Spaniards in Europe during the sixteenth century as an ornamental plant species before its oil started to be used for food during the nineteenth century. Currently, sunflower is cultivated on more than the 5 continents, with Ukraine and Russia being the largest producers, followed by the EU-28. In France, sunflower crops cover 590,000 ha, with an annual production of 1.36 M tonnes, resulting in an average yield of 2.3 t.ha⁻¹ (source: FAO, http://www.fao.org/faostat/en, 10-year average between 2007 and 2016). In France, this oilseed crop species is grown mainly in the central-western (Nouvelle Aquitaine) and southwestern (Occitanie) regions. Sunflower is a summer crop; it is mainly grown under rainfed conditions and has a low water-use (Tallec et al., 2013a). In southwestern France, the temperate, Mediterranean-like climate (dry and warm summers) results in sunflower crops being subjected to frequent and severe water stress conditions that impact crop productivity at both field and landscape/regional scales. These variations should be quantified to better manage croplands. Moreover, in this context, accurate estimation of yield, carbon (C) and water budgets at the field level over large areas remains challenging. The heterogeneous character of landscape in terms of climatic and soil conditions as well as the diversity of agricultural practices make the upscaling of agronomic models difficult. Thus, a precise understanding of the cropland C budget is needed, as well as agronomic indicators (e.g. water-use efficiency) useful for monitoring crop species such as sunflower.

Crop models were first used for research purposes (de Witt and Brouwer, 1970; G. F. Arkin et al., 1976) in the early 1970s. During the 1980s, they became more diversified and complex to meet management needs and to support field-scale decision making ((G. G. Wilkerson et al., 1983b), (D. P. Swaney et al., 1983) and farming applications (McKinion et al., 1989)). Today, a considerable number of crop models exist. They can be classified into 3 groups:

C-driven models, or what we call the “de Wit school” crop growth simulation models (Bouman et al., 1996a), are based on the assimilation of C by photosynthetic parts of the plant, such as the SUCROS ((van Laar, 1992)) or GECROS ((Yin, 2005)) models. These models, on the basis of their hierarchical architecture, are complex, require numerous parameters, and are difficult to calibrate; thus, they are difficult to upscale to large areas.

Water-driven models are based on a predetermined relation between crop growth rate and transpiration according to water productivity parameters. Compared to that of the previous class, the approach of this class results in models whose structure is less complex, leading to fewer parameters. Among the most representative ones are the AquaCrop (Steduto et al., 2009b) or AqYield (Constantin et al., 2015) models.

Radiation-driven models, which have been constructed on the basis of (Monteith and Moss, 1977b) approach, convert intercepted or absorbed solar radiation into aboveground or total biomass via a radiation-use efficiency parameter, significantly reducing the number of photosynthetic parameters needed. Most crop models are radiation driven, including the CERES (Jones et al., 1986b), EPIC (J. R. Williams et al., 1989b), STICS (Brisson et al., 2003a), SAFY (Duchemin et al., 2008b) and SAFY-WB (Duchemin et al., 2015b) models.

The model presented in the current study belongs to both the first and third groups. Indeed, it evolved from i) the SAFY-WB model (Duchemin et al., 2015b), which belongs to the third group and
allows the simulation of biomass production, yield and water fluxes and budget, and from ii) the SAFY-CO₂ model (Pique et al., 2020a), which calculates biomass production as the difference between plant assimilation and respiration.

Regarding sunflower, crop growth models are either generic, with a set of parameters specific to sunflower (AquaCrop, CropSyst (Stöckle et al., 2003b), EPIC, STICS or WOFOST (Diepen et al., 1989b)) or specific to sunflower (OilCrop-SUN, (Villalobos et al., 1996); SUNFLO, (Casadebaig et al., 2011b)).

To reduce the number of model parameters and to increase the robustness of the models upon upscaling, several remote sensing data assimilation approaches for crop models have been developed. The benefits of using HSTR satellite products in crop models have already been demonstrated in several studies (Battude et al., 2016b; Pique et al., 2020a; Sibley et al., 2014; Sus et al., 2013; Wu et al., 2012). These approaches are either based on establishing relations between remote sensing products and crop productivity (Running et al., 2000) or on calibrating the model's parameters (Cicek et al., 2010; Panda et al., 2010; Shanahan et al., 2001). Such approaches have already been applied to sunflower crops for biomass estimation (Claverie et al., 2012b; Ferrant et al., 2014). Those studies showed that accounting for the spatial heterogeneity in sunflower development is crucial for the accurate estimation of crop dynamics; nevertheless, these models are not able to estimate cropland C budgets. In this context, we propose a newly-developed agro-meteorological model, the Simple Algorithm For Yield, Evapotranspiration and CO₂ fluxes estimates (SAFYE-CO₂), that assimilates HSTR products, and we validate the model for use in sunflower. Adapted from models developed by (Pique et al., 2020a) and [17], this model allows us to simulate cropland growth, yield and CO₂ and water fluxes at a large scale because it requires little input data (e.g. little to no crop management data).

Initially, the Simple Algorithm For Yield estimates (SAFY) model was developed to estimate crop biomass and production and was first validated for use in winter wheat (Duchemin et al., 2008a), followed by sunflower and maize (Claverie et al., 2012a; Battude et al., 2016a). Those studies demonstrated the benefits of assimilating the HSTR green area index (GAI) derived from optical remote sensing satellite data for crop monitoring by calibrating several phenological parameters as well as light-use efficiency. Later, this model was coupled with a water module based on the FAO-56 approach (Allen et al., 1998) to estimate cropland evaporation (Ev) and transpiration (TR) and thus water stress and water needs. This resulted in the SAFY-WB (water balance) model, which was validated for use in wheat (Duchemin et al., 2015a) and maize (Battude et al., 2017a). More recently, Pique *et al.* (Pique et al., 2020a) showed the benefit of using HSTR GAI products for assessing the components of the C budget over winter wheat crops by adapting the SAFY model to a C module, which led to a new model called SAFY-CO₂. In line with these studies, a new model, coupling the SAFY-CO₂ model with the same water module as that in the SAFY-WB model, has been constructed and is evaluated here. The new model is called the SAFYE-CO₂ model [E stands for evapotranspiration (ETR)], which aims to assess daily biomass production and CO₂ and water fluxes, as well as annual yields and the C and water budget over croplands. The development and validation of the SAFYE-CO₂ model for sunflower was driven by the following perspectives: i) simulation of the water flux and of the soil water provides an opportunity to estimate the effects of water stress on crop growth, which could improve biomass and CO₂ flux estimates for summer crops (e.g. sunflower), and ii) SAFY-CO₂ has already been validated for use in winter wheat (Pique et al., 2020a), and validation for use in sunflower would allow us to simulate C and water budgets for sunflower-winter wheat crop rotations, which are very common in southwestern France.
In terms of upscaling the model over large areas and simulating the effects of water stress on
crop production and C budgets, the objective of this study is twofold: i) to analyze whether a coupled
remote sensing/modeling approach such as that involving SAFYE-CO₂ is able to accurately reproduce
sunflower function, i.e. biomass, yield and CO₂ and water fluxes, at the ground-atmosphere interface
and ii) to evaluate the impact of current soil map products on the simulation accuracy as well as the
additive value of coupling SAFYE-CO₂ with a water module to improve sunflower simulation.

In the following sections, the materials and methods are described. Afterward, the results
section starts with a validation of the SAFYE-CO₂ model against biomass, yield, and CO₂ and water
flux measurements at the Auradé integrated carbon observation system (ICOS) instrumented site. A
sensitivity analysis of the model to the soil parameters is then performed. Finally, SAFYE-CO₂ is
validated via biomass and yield data collected from more than 80 fields during five years, covering a
wide range of soil and climatic conditions. Maps of the components of the C budget over a cropping
season are then presented. The study ends with a discussion on the benefits and limitations of our
approach, raising issues concerning the assessment of accurate C and water budgets at the plot scale
over large areas and discussing potential improvements for simulating those variables.

**4.3. Materials**

The proposed study used a large amount of data from several sources to calibrate, validate and
upscale our model. In addition to meteorological inputs and soil property data needed for the model,
high-resolution remote sensing data were required for calibration. The validation of the model is based
on CO2 and water flux measurements, biomass production and yield. The last two originated either
from destructive measurements carried out at the FR-Aur site (see section II-4.3.2.2) or from surveys
carried out in 2013, 2014 and 2015 (Champolivier et al., 2019). The following sections present these
data as well as the study area.

**4.3.1. Study area**

The study area is located in southwestern France and encompasses 4 administrative departments (Haute-Garonne, Gers, Aude and Ariège, see Figure II-6). It is part of the Regional Spatial Observatory – Southwest [OSR-SO, https://osr.cesbio.cnrs.fr/] located near Toulouse and includes the instrumented ICOS agricultural site of Auradé (FR-Aur). The region is characterized by a temperate and Mediterranean-like climate, with dry and warm summers (see Table II-1). The Garonne River flows across the study area from south to north and has shaped the topography into terraces and hill landscapes that resulted from the erosion of old terraces. This has led to heterogeneous soils called "boulbènes" (silty soil on terraces) and "terrefort" (clayey-limestone soils on hills), characterized by contrasting textures (especially with respect to the clay content) and with variations in depth within each type of soil.

According to the 2015 French Land Parcel Identification System (LPIS, or the Registre Parcellaire Graphique in French), agricultural activity represents 66% of the study area (focus in Figure 1), which includes forests, grasslands, orchards, etc. Sunflower crops constituted 19.1% of the agricultural area.

In the region, sunflower is usually sown between April and May and is harvested between the
end of August and September. The climatic conditions during the crop cycle of sunflower are
characterized by high radiation and strong evaporative demand, while low and erratic rainfall occurs in
the summer (Table 1). The selected study area is very seldom irrigated, like elsewhere in France (<5%
of the sunflower area; (Champolivier, 2011)). Owing to its deep root system, its thick and hairy epidermis and its adaptive capability (through a reduction of the photosynthetic rates), sunflower is moderately sensitive to water deficit. However, hydric stress occurring before pre-flowering can lead to early senescence, ultimately causing losses in yield and oil content.

### 4.3.2. Flux site and experimental instrumentation

Since 2005, the FR-Aur site (Figure II-6) has been highly instrumented for the continuous collection of micrometeorological, meteorological, soil and vegetation data (see (Béziat et al., 2009a) for more details) in the context of the OSR-SO. FR-Aur is also part of the European Research Infrastructure Consortium ICOS [https://www.icos-cp.eu/], which is a European network documenting long-term greenhouse gas concentrations in the atmosphere and greenhouse gas exchange between the surface and the atmosphere in several ecosystems under different pedoclimatic conditions (crops, grasslands, forests, etc.). The plot area is 23.5 ha and is highly heterogeneous because of the topography (17.3 m height difference), which results in different soil depths (from 100 cm to more than 200 cm) and thus different crop yield potentials. The FR-Aur field belongs to a conventional grain farm where only mineral fertilization is applied. The plot is characterized by a rapeseed/winter wheat/sunflower/winter wheat rainfed rotation. Sunflower was grown during the 2006-2007 and 2015-2016 cropping seasons (see Table II-1 for sowing and harvest dates). Hereafter, the crop years will be abbreviated by “CY” followed by the harvest year (e.g. CY-07 for the 2006-2007 crop season).

<figure><img src="image_51.png" /></figure>

Figure II-6: Top left: Fr-Aur field and the location of the measurement station; Top-right: Location of the study area in France; Bottom: Location of the fields in which biomass and yield data were collected in 2013, 2014 and 2015, as well as the selected area for the model's upscaling (see section 4.3.3).

#### 4.3.2.1. Meteorological data

The SAFYE-CO₂ model requires meteorological data as the main input, while vegetation remote sensing was used to constrain the modeled dynamics. Specifically, air temperature (T AIR) at 2 m, incoming global radiation (R G), rain and reference ETR (ET 0) were needed. These variables are
provided either from *in situ* measurements acquired at the FR-Aur site for local validation or from SAFRAN reanalysis data (Durand et al., 1993a) for the spatial analysis.

At FR-Aur, meteorological data are recorded every 30 minutes and averaged on a daily basis (see (Béziat et al., 2009a) for more details). These data were used to calculate the ET₀ using the Penman-Monteith formula (Allen et al., 1998).

<table><caption>**Table II-1 :** Sowing and harvest dates of sunflower cultivated at FR-Aur and climatic variables averaged from the 15th of April to the 15th of September for all crop years.</caption><thead><tr><th>Cropping year</th><th>Sowing</th><th>Harvest</th><th>Mean Temperature [°C]</th><th>Mean Radiation [MJ.m<sup>-2</sup>]</th><th>Cumulative precipitation [mm]</th></tr></thead><tbody><tr><td>2007</td><td>10th of April</td><td>20th of September</td><td>18.0</td><td>19.3</td><td>339</td></tr><tr><td>2013</td><td>-</td><td>-</td><td>17.3</td><td>20.1</td><td>364</td></tr><tr><td>2014</td><td>-</td><td>-</td><td>17.6</td><td>19.8</td><td>178</td></tr><tr><td>2015</td><td>-</td><td>-</td><td>19.0</td><td>20.2</td><td>275</td></tr><tr><td>2016</td><td>20th of April</td><td>23th of September</td><td>18.1</td><td>20.8</td><td>294</td></tr><tr><td>Average 2007-2016</td><td>-</td><td>-</td><td>18.2</td><td>20.3</td><td>267</td></tr></tbody></table>

The SAFRAN meteorological data were used as inputs to simulate the sunflower experiments of 2013, 2014 and 2015. SAFRAN reanalysis data are provided over a grid of 8 x 8 km² over France. In this study, the climatic data were estimated for each field using a bilinear interpolation of the 4 closest grid points.

The climatic conditions of all simulated years are summarized in Table II-1. Each variable is averaged from the 15th of April to the 15th of September, corresponding to the average sunflower growing season in the region. The crop cycles of CY-07 and CY-16 were quite similar, although 2016 was slightly hotter and sunnier than 2007 (Table II-1), leading to slightly higher biomass production in 2016. In 2013 and 2014, the growing seasons were colder compared to the season during 2015. The 2013 growing season (CY-13) was the wettest, and the CY-14 growing season was the driest. Note for the sake of clarity, that only sowing and harvest dates of CY-07 and CY-16 are indicated in Table II-1. Indeed, more than one field was simulated the others CY (respectively 9, 21 and 49 for CY-13, CY-14 and CY-15), thus the sowing and harvest dates vary by plot which explained the empty cases of Table 1.

### 4.3.2.2. Biomass and yield data

Destructive measurements were made at FR-Aur to characterize the sunflower dry aboveground mass (DAM) dynamics during CY-07 and CY-16 (see Table II-2). In CY-07, 30 plants were collected monthly inside the flux tower footprint area distributed on either side of the eddy covariance tower towards the prevailing winds. The DAM was obtained by drying samples in an oven at 65°C for 48 h. The total DAM produced (before harvest) was measured by cutting plants from a 24 m² area and measuring them. The grain yield (YLD) was provided by the farmers for CY-07, and the grains were separated from the heads and weighed to obtain yield for CY-16.

Extensive field campaigns were conducted in 2013, 2014 and 2015 to collect biomass and yield data across the study area (Figure II-6). In 2013, 9 ESUs (20 x 20 m²) corresponding to 8 different fields were monitored during the growing season. For each ESU and at 4 times during crop
development (the last sampling was performed just before harvest), 5 samples of aboveground vegetation consisting of 10 adjacent plants were collected (except for the first date, when only 3 plants were collected), oven dried at 50°C for 48 h and then weighed. DAM was calculated by dividing the total biomass by the area sampled, calculated as the product of the interrow width and the length of the sampling area. On the last sampling date, the grains were separated from the heads and weighed to calculate the yield.

<table><caption>Table II-2: Dataset of DAM and yield as well as their measurement scale.</caption><thead><tr><th>Year</th><th>Number of fields monitored</th><th>Dry Aboveground Mass</th><th>Final Yield</th><th>Workable final yield</th><th>CO<sub>2</sub> and H<sub>2</sub>O fluxes</th><th>Scale</th></tr></thead><tbody><tr><td>2007</td><td>FR-Aur</td><td>7</td><td>1</td><td>1</td><td>Yes</td><td>Field/ESU*</td></tr><tr><td>2013</td><td>9</td><td>36</td><td>9</td><td>9</td><td>No</td><td>ESU</td></tr><tr><td>2014</td><td>70</td><td>-</td><td>36</td><td>21</td><td>No</td><td>Field</td></tr><tr><td>2015</td><td>116</td><td>-</td><td>79</td><td>49</td><td>No</td><td>Field</td></tr><tr><td>2016</td><td>FR-Aur</td><td>3</td><td>1</td><td>1</td><td>Yes</td><td>ESU</td></tr></tbody></table>

*Elementary Sampling Units

In 2014 and 2015, yield data were collected for more than 70 and 116 fields, respectively. They were provided by the farmers, and because some of the yield values were averaged over several fields, only yield data corresponding to one field were considered (see Table II-2). Two or three visits were made to each field in the 2014 and 2015 campaigns to characterize the crop growth status (e.g. canopy heterogeneity, presence of weeds, pest and disease pressure). From those observations, we defined quality scores, and low-quality data were discarded. These data corresponded to fields where i) vegetation-free areas larger than 2 m² were observed and ii) weeds covered more than 50% of the field. After discarding all these data (115), we kept data for 70 fields for the model evaluation (Table II-2). None of the fields monitored in 2013, 2014 and 2015 were irrigated.

### 4.3.2.3. Flux data

Turbulent fluxes of CO₂, water vapor (ETR and latent heat), sensible heat, and momentum were measured continuously using the eddy covariance method (Aubinet et al., 2012a; Baldocchi, 2003a; Moncrieff et al., 1997a) by the combination of an open-path fast infrared gas analyzer (LI-7500, LI-COR) and a 3D sonic anemometer (CSAT3, Campbell Scientific). The flux tower was installed near the center of the field (see Figure II-6). EdiRe software (Clement, 1999) was used to calculate the ETR and net CO₂ fluxes (see (Béziat et al., 2009a) and (Pique et al., 2020a) for details). Flux filtering, quality controls, and gap filling procedures were performed in accordance with the CarboEurope-IP recommendations (http://www.carboeurope.org/).

The net ecosystem exchange (NEE) of CO₂ was partitioned into gross primary productivity (GPP) and ecosystem respiration (R<sub>ECO</sub>) according to the methods proposed by (Reichstein et al., 2005) and adapted for croplands by (Béziat et al., 2009a). During the period when bare soil was exposed (deduced from field observations), the GPP was set to 0, while the R<sub>ECO</sub> was set equal to the NEE.

### 4.3.2.4. Soil data

To run the SAFYE-CO₂ model, data on the soil moisture at field capacity (θ<sub>FC</sub>) and the permanent wilting point (θ<sub>WP</sub>) are needed for each considered field to estimate the daily soil available water content (AWC). This information can be derived either from *in situ* sampling and laboratory
analysis or can be inferred from soil texture products by applying pedotransfer functions (PTFs). In this study, we used both methodologies: data from *in situ* soil sampling at FR-Aur and data from the global soil databases GlobalSoilMap (GSM, Román Dobarco et al., 2019) and SoilGrids (SG, Hengl et al., 2017) at both FR-Aur and at the other 78 fields. These soil products were used to estimate the θ<sub>FC</sub>, θ<sub>WP</sub> and soil depth for all the considered fields in this study, including those at FR-Aur. The global soil databases provide the fraction of sand, clay and silt at a 90 m (GSM) or 250 m (SG) resolution as well as the soil depth. SG is based on machine learning and remote sensing data, while GSM uses a data-mining technique and a bottom-up approach (Minasny and McBratney, 2010).

<table><caption>**Table II-3 :** Soil characteristics derived from *in situ* measurements and GlobalSoilMap and SoilGrids data at FR-Aur corresponding to 3 layers simulated by the model.</caption><thead><tr><th rowspan="2">Source of data</th><th colspan="3">Field Capacity [m<sup>3</sup>.m<sup>-3</sup>]</th><th colspan="3">Wilting Point [m<sup>3</sup>.m<sup>-3</sup>]</th><th colspan="3">Available Water Capacity [m<sup>3</sup>.m<sup>-3</sup>]</th><th rowspan="2">Soil Depth [cm]</th></tr><tr><th>Surface layer</th><th>Root layer</th><th>Deep layer</th><th>Surface layer</th><th>Root layer</th><th>Deep layer</th><th>Surface layer</th><th>Root layer</th><th>Deep layer</th></tr></thead><tbody><tr><td>In situ</td><td>0.31</td><td>0.39</td><td>0.39</td><td>0.13</td><td>0.19</td><td>0.19</td><td>0.18</td><td>0.2</td><td>0.2</td><td>[100-200]</td></tr><tr><td>GSM</td><td>0.29</td><td>0.3</td><td>0.29</td><td>0.18</td><td>0.17</td><td>0.16</td><td>0.11</td><td>0.13</td><td>0.13</td><td>159</td></tr><tr><td>SoilGrids</td><td>0.34</td><td>0.33</td><td>0.33</td><td>0.25</td><td>0.21</td><td>0.21</td><td>0.09</td><td>0.12</td><td>0.12</td><td>188</td></tr></tbody></table>

<figure><img src="image_52.png" /><figcaption>**Figure II-7 :** Soil type distribution derived from GlobalSoilMap (left) and SoilGrids (right) products based on the USDA classification system. The red, black and blue circles respectively correspond to fields in the 2013, 2014 and 2015 seasons.</figcaption></figure>

We used the PTFs proposed by (Mercedes Román Dobarco et al., 2019) (eqs. 9-10) to calculate the volumetric soil moisture at field capacity (θ<sub>FC</sub>) and at permanent wilting point (θ<sub>WP</sub>). The estimations of the θ<sub>FC</sub>, θ<sub>WP</sub> and soil depth for FR-Aur using *in situ* samples, SG and GSM products are summarized in Table II-3 for the 3 soil layers considered in the proposed modeling approach (see section II-4-4.1).

The soil types of all the fields monitored in 2013, 2014 and 2015 were characterized based on the textures derived from the GSM and SG products (for the top 30 cm of the soil) and considering the United States Department of Agriculture (USDA) classification system (see Figure II-7). Considering all monitored fields, the GSM products suggest that 64% of the soils are clay loams, 23% are loams, 11% are silty clay loams, and the remaining soils vary between clay, silty clay and silty loams. Compared with the GSM products, the SG products suggest less heterogeneity among the soil types and classify 87% of the soils as clay loams, 12% as loams and less than 1% as silty clay loams.

### 4.3.3. Satellite data and products

#### 4.3.3.1. Multisatellite optical images

Liu et al., (2010) and Lobell et al., (2003) highlighted the usefulness of assimilating HSTR satellite products into crop models, as it increased the accuracy in reproducing crop dynamics and/or reduced the number of model parameters. Thus, the accuracy of our approach is partly determined by
the availability of remote sensing images. Indeed, to successfully calibrate some of the model parameters (see section II-4.4.2), a certain number of images at key dates had to be available. An extensive set of HSTR products derived from several sensors was thus used in this study. The main characteristics in terms of spatial resolution, temporal sampling and spectral bandwidth of each sensor are detailed in Table II-4, and the number of images (all in all, 147) from each sensor is shown in Figure II-8. For the purpose of the study, a temporal series of remote sensing products, averaged at plot scale, are needed. Products of different spatial (from 8 to 30 m) and temporal (from 1 to 16 days) scales are thus used in order to get the most complete time series (Fig. II-8).

<figure><img src="image_53.png" /><figcaption>Figure II-8 : Timeline of images for each year of study and each sensor.</figcaption></figure>

<table><caption>Table II-4 : Main characteristics of the satellite sensors used in this study and the bands used to derive the GAI.</caption><thead><tr><th></th><th>SPOT</th><th>Formosat-2</th><th>DEIMOS-1</th><th>Landsat-8</th><th>Sentinel-2A</th></tr></thead><tbody><tr><th>Spatial resolution</th><td>20 m (SPOT2-4) 10 m (SPOT5)</td><td>8 m</td><td>22 m</td><td>30 m</td><td>10 m</td></tr><tr><th>Revisit frequency</th><td>Between 1 and 3 days</td><td>1 day</td><td>Between 3 and 4 days</td><td>16 days</td><td>5 days</td></tr><tr><th>Footprint</th><td>60 km*60 km</td><td>24 km*24 km</td><td>Up to 625 km</td><td>(170 km*) 185 km</td><td>290 km</td></tr><tr><th>Wavelength (µm)</th><td>Green: [0.5-0.59];<br>Red: [0.61-0.68];<br>NIR: [0.79-0.89];<br>SWIR: [1.53-1.75]</td><td>Green: [0.52-0.6];<br>Red: [0.63-0.69];<br>NIR: [0.79-0.9]</td><td>Green: [0.52-0.6];<br>Red: [0.63-0.69];<br>NIR: [0.77-0.9]</td><td>Green: [0.53-0.59];<br>Red: [0.64-0.67];<br>NIR: [0.85-0.88];<br>SWIR1: [1.57-1.65];<br>SWIR2: [2.11-2.29]</td><td>Green: [0.54-0.58];<br>Red: [0.65-0.68];<br>VRE1: [0.7-0.71];<br>VRE2: [0.73-0.75];<br>VRE3: [0.77-0.79];<br>NIR: [0.72-0.95];<br>Narrow NIR: [0.85-0.87];<br>SWIR1: [1.57-1.66];<br>SWIR2: [2.11-2.29]</td></tr><tr><th>Viewing angle</th><td>Constant</td><td>Constant</td><td>Variable</td><td>Constant</td><td>Constant</td></tr></tbody></table>

(Battude et al., 2016a) showed that the advantages in terms of time sampling of using GAI estimates from multiple sensors outweigh the drawbacks in terms of sensor measurement mismatch. In their study, the authors compared the normalized difference vegetation index (NDVI) of different sensors and showed that, during the growing season (NDVI > 0.2) and over various crops (sunflower, maize, wheat, etc.), linear relationships [R values between 0.97 and 0.99 and relative root mean square error (rRMSE) values between 5.3 and 10.9%] exist between pairs of sensors [including DEIMOS-1, Formosat-2, Landsat-8, and SPOT4 (Take5) (Hagolle et al., 2015)]. The KALIDEOS processing chain [http://kalideos.cnes.fr] for atmospheric, radiometric and geometric corrections was used to derive
surface reflectance values from the satellite data. The mean geometric correction accuracy was close to
0.2 pixels (LaFrance, 2012), which is satisfactory considering the area of the studied fields.

4.3.3.2. From image reflectance to GAI estimates

The remotely sensed GAI estimates were derived from surface reflectance values using the Biophysical Variables Neural NETwork (BV-NNET) tool (Baret et al., 2007a). BV-NNET is a trained artificial neural network based on the inversion of the radiative transfer model PROSAIL (Jacquemoud et al., 2009b). The GAI estimates take into account the spectral and directional information of each sensor, which are of paramount importance when working with multiple satellites. The GAI value estimated by BV-NNET is an “effective GAI”, which does not consider the clumping effect and could lead to underestimation compared to that of the “true GAI”, mainly for high GAI values (Demarez et al., 2008a). For each simulated field, the mean GAI was calculated taking into account all the pixels of the considered plot after the application of an offset of 10 m to avoid edge effects.

**4.4. Methods**

**4.4.1. The SAFYE-CO₂ model**

The SAFYE-CO₂ model was constructed by coupling the SAFY-CO₂ model (Pique et al., 2020a) with the water module of the SAFY-WB model (Duchemin et al., 2015a) that has been adapted from the FAO-56 method (Allen et al., 1998). SAFY-CO₂ is a parsimonious agro-meteorological model running on a daily basis. It is driven by T_AIR and R_G as meteorological inputs. It simulates crop dynamics (GAI), crop production (DAM, YLD), CO₂ fluxes (GPP, R_ECO and NEE), and the annual net ecosystem C budget (NECB). Coupling with the water module, which can easily be added or removed, is needed to estimate crop TR and soil Ev and thus soil water content (SWC). From those data, soil water stress and plant water needs can be calculated. Adding this module requires two additional meteorological inputs, rain and ET₀.

Since all equations involved in the SAFYE-CO₂ model have already been described in detailed
in other studies (SAFY, (Duchemin et al., 2008a); SAFY-WB, (Duchemin et al., 2015b); SAFY-CO₂,
(Pique et al., 2020a)), we will summarize here the main methods and the interactions between the
water module and the SAFY-CO₂ model. All related parameters discussed in this section are
summarized in Table II-5.

In the SAFYE-CO₂ model, GPP is first estimated as a function of the R<sub>G</sub>, a climatic efficiency (ε<sub>C</sub>), the absorbed photosynthetically active radiation (APAR) calculating from the GAI and the light-interception coefficient (K<sub>EXT</sub>), a function of the light-use efficiency taking into account the fraction of the diffuse radiation (ELUE<sub>A</sub> and ELUE<sub>B</sub>; see (Pique et al., 2020a)), and a function accounting for the decrease in canopy photosynthesis capacity during senescence (sR<sub>10</sub>). The daily estimates of GPP are modulated by i) a thermal stress function defined by three cardinal temperatures limiting crop growth (T<sub>MIN</sub>, T<sub>MAX</sub>, T<sub>OPT</sub>) and ii) a water stress function (WtS<sub>v</sub>, see below). Net primary productivity (NPP) is estimated as the difference between the GPP and autotrophic respiration (R<sub>A</sub>). R<sub>A</sub> is estimated by the methods proposed by (McCree, 1974), which separates R<sub>A</sub> into two components: maintenance respiration (R<sub>M</sub>) and growth respiration (R<sub>GR</sub>). R<sub>M</sub> is computed as a function of the NPP of the previous day, a “Q<sub>10</sub>-type” R<sub>M</sub> equation (Van’t Hoff, 1898) (with two crop-specific parameters: Q<sub>10</sub> and R<sub>10</sub>), and the sR<sub>10</sub> coefficient. R<sub>GR</sub> is estimated considering the daily amount of biomass produced and using a R<sub>GR</sub> coefficient parameter (Y<sub>G</sub>) as proposed by (Amthor, 1989) and later improved by
(Choudhury, 2000). Heterotrophic respiration (R<sub>H</sub>) is estimated from a simple Q<sub>10</sub> first-order exponential equation (Q<sub>10_H</sub>) and a reference respiration parameter (Rh<sub>REF</sub>).

CO₂ fixed by the vegetation is converted into dry biomass via a plant C parameter (C<sub>VEG</sub>) and partitioned between aboveground (DAM) and belowground parts according to the methods proposed by (Baret et al., 1992). This process depends on three parameters (f<sub>0</sub>, f<sub>∞</sub> and c) that were adapted for sunflower crops based on data from (Ma et al., 2017a).

Like for SAFY, the daily increase in the GAI is estimated from the daily increase in DAM, the specific leaf area (SLA), and partition-to-leaf function (Pl<sub>A</sub> and Pl<sub>B</sub>). The decrease in the GAI starting with senescence is estimated based on a threshold from the growing degree days (starting from crop emergence as detected by the remote sensing data), from the senescence parameter (SEN<sub>A</sub>) and from the rate of senescence (SEN<sub>B</sub>). The final yield is simply deduced from the maximum value of DAM multiplied by a harvest index (HI). The set of equations used to describe the previous methods are detailed in the study by (Pique et al., 2020a).

In SAFY-CO₂, the soil is a single layer whose characteristics (parameters) are used to calculate R<sub>H</sub>. In the water module of SAFYE-CO₂, the soil is divided into 3 layers: a surface layer of constant depth (D<sub>EVAP</sub>), subject to Ev and TR: a root layer, subject to TR only, that increases in depth during the growing season at a rate (V<sub>R</sub>) that tracks crop growth and limited by a maximum root depth (D<sub>ROOTS</sub>); and a constant depth (D) deep layer. The last one is used to initialize the water content of the root layer during root growth and to control for the effects of percolation and drainage. The soil module is detailed in the paper of (Duchemin et al., 2015a). The soil AWC is essentially updated daily, taking into account rain (inputs), Ev and TR simulated based on the FAO-56 methods (Allen et al., 1998). The AWC is defined by the water fluxes (input and output) and the total water capacity (TWC) (i.e. TWC = θ<sub>FC</sub> - θ<sub>WP</sub>). TR depends on a maximal TR coefficient (Kcb<sub>MAX</sub>), a reduction coefficient (K<sub>TRP</sub>), a critical relative humidity parameter (Dfe), the GAI, and the maximum relative filling rate between the surface and root layer.

The water module also computes two water stress functions, one that drives Ev (WtS<sub>E</sub>, eq. II-19, considering the AWC of the surface layer) and another that drives the growth of vegetation (WtS<sub>V</sub>, eq. II-20, considering the surface and root layers) through GPP estimation. These functions represent the main link between the vegetation and the soil module and are defined as follows:

$$WtS_E = \left( \frac{SM - \theta_{WP}}{\theta_{FC} - \theta_{WP}} \right) \qquad \text{II-19}$$

$$WtS_V = \min \left\{ 1 - \frac{RH_{TOT} - Dfe}{Dfe}, 1 \right\} \qquad \text{II-20}$$

$$RH_{TOT} = \max \{RH_{evap}, RH_{root}\} \qquad \text{II-21}$$

where SM [m³.m⁻³] is the SWC of the surface layer and where θ<sub>FC</sub> [m³.m⁻³] and θ<sub>WP</sub> [m³.m⁻³] are the water content at field capacity and at the permanent wilting point, respectively. Dfe is the critical relative filling rate for water stress, and RH<sub>evap</sub> and RH<sub>root</sub> are the relative humidity of the surface and root layers, respectively.

In the present study, the calculation of Ev was modified from that of (Duchemin et al., 2015a) according to the methods of (Soares et al., 1988a) because differences were observed between simulated and observed ETR during bare soil periods and because of the dependence of Ev on ET<sub>0</sub>. The reason is that FAO-56 method for estimating Ev was developed to reproduce homogeneous grass
stands, whereas the soil is often completely bare during the fallow period. Additionally, (Soares et al., 1988a) asserted that, in addition to being ET₀ dependent, Ev depends on the hydraulic state of the soil surface layer. When ET₀ is high, the first several centimeters of the soil rapidly become dry, slowing the capillary rise of the remaining water. A new equation for estimating Ev was thus adopted :

$$Ev = K_{BS} * (Ev_{A} * ET_{0} + Ev_{B}) * ET_{0} * (1 - F_{cover}) * WtS_{E} \qquad II-22$$

where ET₀ [mm.d⁻¹] is the reference evapotranspiration, KBS [-] is a correction coefficient for bare soil, EvA [1/(mm.d⁻¹)] and EvB [-] are the evapotranspiration coefficients, and FCOVER is the vegetation cover fraction. FCOVER is estimated from the maximal vegetation fraction cover (KCOVER), an exponential coefficient (ECOV) and the GAI. FCOVER can be derived from remote sensing and thus drives the model, but previous quantitative experiments did not show significant enhancement of the outputs; thus, we calculated it here, from the GAI and two parameters (see Table II-5).

The water module is impacted by the vegetation module through the calculation of FCOVER (which considers the GAI), and in turn, the vegetation module is impacted by the water module through the water stress function (WtSv).

### 4.4.2. Model parameters

In addition to requiring meteorological data, the SAFYE-CO₂ model requires crop-specific parameters. These parameters came from several sources. Most of them could be obtained from the literature, some could be estimated via *in situ* measurements, and the rest were obtained from calibration (see Table II-5). The model could run using literature and *in situ* parameters only, but when the last category of parameters is not constrained by remote sensing, the model cannot cope with changing conditions in the fields, which renders the outputs inaccurate. The calibration procedure based on the satellite-derived GAI allows optimization of the phenological parameters and the A parameter of the light-use efficiency function. Once these parameters are calibrated, the outputs are representative of the GAI extraction area (e.g. pixels, fields, farms). To calculate the plot-scale C budget, data supplied by farmers, concerning the amount of C imported into fields as organic amendments are eventually needed.

### 4.4.3. Model parametrization and calibration

The calibration procedure incorporating remotely sensed GAI was detailed in the report by (Pique et al., 2020a). Six phenological parameters (PIA, PIB, SENA, SENB, SLA, D₀) as well as the A parameter of the light-use efficiency (see Table II-5) were calibrated during this process for each field and each year independently to that the simulated GAI matched the remotely sensed GAI. This procedure allows the model to accurately reproduce crop growth. This inversion process is adapted from the Nelder-Mead simplex method (Lagarias et al., 1998), in which a priori values and minimum/maximum boundaries are set for each retrieved parameter. A quadratic difference cost function—root mean square error (RMSE)—is minimized using the adapted simplex between the remotely sensed GAI and the simulated GAI. For minimization, a global approach based on 30 different a priori samples is used to avoid local minima while preventing unnecessary runs. Once completed, the inversion process results in a set of parameters allowing accurate reproduction of the remotely sensed GAI while avoiding biomass estimates as well as water and CO₂ flux estimates (see section II-4.4.1).
<table><caption>Table II-5: List of the SAFYE-CO2 model parameters for sunflower, values or a priori intervals of the parameters, method for setting them, and relevant sources.</caption><thead><tr><th>Model</th><th>Name</th><th>Notation</th><th>Unit</th><th>Value</th><th>Method</th><th>Source</th></tr></thead><tbody><tr><td></td><td>Climatic efficiency</td><td>ε<sub>C</sub></td><td>-</td><td>0.48</td><td>Literature</td><td>Varlet grancher 1982</td></tr><tr><td></td><td>Light interception coefficient</td><td>K<sub>EXT</sub></td><td>-</td><td>0.8</td><td>Literature</td><td>Flénet 96/Steer 96</td></tr><tr><td></td><td>Minimal temperature for growth</td><td>T<sub>MIN</sub></td><td>°C</td><td>8</td><td>Literature</td><td>STICS</td></tr><tr><td></td><td>Maximal temperature for growth</td><td>T<sub>MAX</sub></td><td>°C</td><td>42</td><td>Literature</td><td>STICS</td></tr><tr><td></td><td>Optimal temperature for growth</td><td>T<sub>OPT</sub></td><td>°C</td><td>25</td><td>Literature</td><td>STICS</td></tr><tr><td></td><td>Polynomial degree of thermal stress function</td><td>β</td><td>-</td><td>2</td><td>Literature</td><td>Duchemin 2008</td></tr><tr><td rowspan="7">SAFY model</td><td>Harvest index</td><td>HI</td><td>-</td><td>0.29</td><td>Literature</td><td>Andrade 1995</td></tr><tr><td>Specific leaf area</td><td>SLA</td><td>m<sup>2</sup>.g<sup>-1</sup></td><td>[10<sup>-3</sup> - 4.10<sup>-2</sup>]</td><td>Calibration</td><td>-</td></tr><tr><td>Partition-to-leaf function parameter A</td><td>Pl<sub>A</sub></td><td>-</td><td>[0.01 - 0.5]</td><td>Calibration</td><td>-</td></tr><tr><td>Partition-to-leaf function parameter B</td><td>Pl<sub>B</sub></td><td>-</td><td>[10<sup>-4</sup> - 10<sup>-2</sup>]</td><td>Calibration</td><td>-</td></tr><tr><td>Sum of temperature for senescence</td><td>SEN<sub>A</sub></td><td>°C</td><td>[200 - 1800]</td><td>Calibration</td><td>-</td></tr><tr><td>Rate of senescence</td><td>SEN<sub>B</sub></td><td>°C.d<sup>-1</sup></td><td>[10<sup>2</sup> - 3.10<sup>4</sup>]</td><td>Calibration</td><td>-</td></tr><tr><td>Day of plant emergence</td><td>Emerg</td><td>day</td><td>[100 - 200]</td><td>Calibration</td><td>-</td></tr><tr><td rowspan="10">Water module</td><td>Depth of evaporative layer</td><td>D<sub>EVAP</sub></td><td>m</td><td>0.1</td><td>Literature</td><td>FAO-56 Allen 1998</td></tr><tr><td>Maximum depth of the root system</td><td>D<sub>ROOTS</sub></td><td>m</td><td>1.4</td><td>Literature</td><td>Stone 2001</td></tr><tr><td>Soil depth</td><td>D</td><td>m</td><td>-</td><td>Soil Data</td><td>-</td></tr><tr><td>Humidity at field capacity</td><td>θ<sub>FC</sub></td><td>m<sup>3</sup>.m<sup>-3</sup></td><td>-</td><td>Soil Data</td><td>-</td></tr><tr><td>Humidity at the permanent wilting point</td><td>θ<sub>WP</sub></td><td>m<sup>3</sup>.m<sup>-3</sup></td><td>-</td><td>Soil Data</td><td>-</td></tr><tr><td>Root mean growth rate</td><td>V<sub>R</sub></td><td>m.<sup>-1</sup>°C<sup>-1</sup></td><td>0.0025</td><td>Literature</td><td>STICS</td></tr><tr><td>Evaporation coefficient A</td><td>Ev<sub>A</sub></td><td>mm<sup>-1</sup></td><td>0.09</td><td>Measurements</td><td>In situ ETR measurements</td></tr><tr><td>Evaporation coefficient B</td><td>Ev<sub>B</sub></td><td>-</td><td>1.09</td><td>Measurements</td><td>In situ ETR measurements</td></tr><tr><td>Evaporation coefficient C</td><td>K<sub>BS</sub></td><td>-</td><td>1.2</td><td>Measurements</td><td>In situ ETR measurements</td></tr><tr><td>Maximal transpiration coefficient</td><td>Kcb<sub>MAX</sub></td><td>-</td><td>1</td><td>Literature</td><td>FAO-56 Allen 98</td></tr><tr><td>Critical humidity parameter</td><td>Dfe</td><td>-</td><td>0.41</td><td>Literature</td><td>Casadebaig 2008</td></tr><tr><td>Maximal vegetation fraction cover</td><td>K<sub>COV</sub></td><td>-</td><td>0.92</td><td>Measurements</td><td>ICOS site data</td></tr><tr><td>Fraction cover exponential coefficient</td><td>E<sub>COV</sub></td><td>-</td><td>0.62</td><td>Measurements</td><td>ICOS site data</td></tr><tr><td>Transpiration reduction coefficient</td><td>K<sub>TRP</sub></td><td>1.01</td><td>Measure</td><td>Measurements</td><td>In situ ETR measurements</td></tr><tr><td rowspan="10">Carbon module</td><td>Effective light-use efficiency parameter A</td><td>ELUE<sub>A</sub></td><td>gC.MJ<sup>-1</sup></td><td>[0.75 - 0.95]</td><td>Calibration</td><td>-</td></tr><tr><td>Effective light-use efficiency parameter B</td><td>ELUE<sub>B</sub></td><td>-</td><td>1.14</td><td>Measurements</td><td>In-situ measurements</td></tr><tr><td>Q<sub>10</sub> maintenance respiration parameter</td><td>Q<sub>10</sub></td><td>-</td><td>2</td><td>Literature</td><td>Amthor 2000</td></tr><tr><td>Reference maintenance respiration at 10°C</td><td>R<sub>10</sub></td><td>gC.gDM<sup>-1</sup></td><td>0.0036</td><td>Literature</td><td>Szaniawski (1982)</td></tr><tr><td>Growth respiration conversion efficiency parameter</td><td>Y<sub>G</sub></td><td>-</td><td>0.76</td><td>Literature</td><td>Szaniawski (1982)</td></tr><tr><td>Root fraction parameters</td><td>f<sub>0</sub> - f<sub>∞</sub> - c</td><td>-</td><td>0.3 / 0.1 / 1.5</td><td>Literature</td><td>Baret 1992/Ma 2017</td></tr><tr><td>Carbon content coefficient</td><td>C<sub>VEG</sub></td><td>gC.gVEG<sup>-1</sup></td><td>0.46</td><td>Literature</td><td>Béziat 2009</td></tr><tr><td>Reference heterotrophic respiration at 0°C</td><td>Rh<sub>REF</sub></td><td>gC.m<sup>-2</sup>.d<sup>-1</sup></td><td>0.34</td><td>Literature</td><td>Suleau 2011</td></tr><tr><td>Q<sub>10</sub> heterotrophic respiration parameter</td><td>Q<sub>10_H</sub></td><td>-</td><td>2.3</td><td>Literature</td><td>Suleau 2011</td></tr><tr><td>Conversion factor of Tair in Tsoil</td><td>t</td><td>-</td><td>1.07</td><td>Measurements</td><td>ICOS site data</td></tr></tbody></table>

### 4.4.4. Model validation and sensitivity analysis

The validation strategy adopted in this study is summarized in Figure II-9. It relies on several types of datasets. First, we evaluated the ability of SAFYE-CO₂ to simulate crop growth and development [leaf area index (LAI) and DAM], final production (YLD) and CO₂ (GPP, R<sub>ECO</sub> and NEE) and water (Ev, TR) fluxes. To accomplish this, the model was run for the two crop seasons of sunflower at the FR-Aur site (CY-07 and CY-16) using *in situ* (IS) meteorological data. Additionally, *in situ* measurements were used to estimate θ<sub>FC</sub> and θ<sub>WP</sub> and the soil depth. These results are presented in section II-4.5.1.

For upscaling of the model, the impact of the soil products (GSM and SG were used to estimate θ<sub>FC</sub>, θ<sub>WP</sub> and soil depth) on the model's performance must be evaluated. For this purpose, the SAFYE-CO₂ model was run on the same two crop years at FR-Aur using the GSM or the SG products to estimate soil water retention parameters. The outputs were compared with those obtained with *in situ* measurements instead of soil products. Additionally, the results obtained with SAFYE-CO₂ were compared to those from the original SAFYE-CO₂ version to assess whether the coupling of SAFYE-CO₂ with the water module improves the model's outputs or not; section II-4.5.2.1 describes those results. Finally, the SAFYE-CO₂ and SAFYE-CO₂ models were used to simulate DAM and YLD over a large
area. In addition to the two crop years at FR-Aur, the fields from the 2013, 2014 and 2015 campaigns were simulated, and *in situ* measurements of DAM and YLD were used for validation. *In situ* meteorological variables were used to simulate the crop years at FR-Aur, while data from SAFRAN were used to simulate the other fields/crop years. The findings related to this exercise are described in section II-4.5.2.2. All simulations were performed from 15 November to 14 November of the following year, usually representing one cropping season in the area of study.

<figure><img src="image_54.png" /><figcaption>Figure II-9 : Summary of the validation strategy in terms of upscaling the SAFY-CO2 and SAFYE-CO2 models. IS is in situ, GSM is GlobalSoilMap and SG is SoilGrids.</figcaption></figure>

## 4.5. Results

### 4.5.1. Local assessment of the model's performances

The SAFYE-CO₂ model was first tested against the *in situ* data of the FR-Aur site. Figure II-10 presents the correlation between the simulated and observed GAI, DAM and CO₂ (GPP, R<sub>ECO</sub> and NEE) and water (ETR) fluxes for the two crop years of sunflower. Only data corresponding to the crop cycle (i.e. from emergence to harvest) were included to compute the GPP statistics, while the entire simulation period was considered to compute the model performances when estimating the other flux variables.

Since the model is calibrated by minimizing the quadratic difference between simulated and remotely sensed GAI, the estimation of GAI is used to verify that the inversion procedure is correctly performed rather than used to validate the modeled GAI. Consequently, the statistical results related to the GAI estimation are very good (R² = 0.99, RMSE = 0.099 m².m⁻² and rRMSE = 6.9 %), indicating that the sunflower's GAI dynamics are accurately reproduced by the model. The DAM is also robustly estimated in terms of correlation (R² = 0.96) and error (RMSE = 66 g.m⁻² and rRMSE = 16.2 %) data, and the yield (only 2 values at FR-Aur) is also estimated with a low error (RMSE = 53.8 g.m⁻²). Accurate estimation of biomass, and thus yield, is a prerequisite for assessing NECB with high precision. CO₂ fluxes are correctly simulated, although compared with CY-07 CO₂ fluxes, CY-16 CO₂ fluxes are estimated with lower accuracy due to the lack of satellite images throughout the agricultural
season. Indeed, for CY-16, only 4 images were available during sunflower development (it was the beginning of Sentinel's mission). In particular, no images were available from 30 July to 28 September (see Figure II-8), leading to a senescence period that was poorly documented by the remotely sensed GAI and a model-simulated senescence stage that was too early. The poor representation of the senescence phase in 2016 led to an underestimation of GPP and R<sub>ECO</sub>. This period is easily identifiable on the scatter plots (Figure II-10 C-D) and corresponds to the green points underestimated by the model. Indeed, the mean bias of the simulated GPP and R<sub>ECO</sub> are higher for CY-16 (respectively -2.21 and -0.69 gC.m<sup>-2</sup>.d<sup>-1</sup>) than for CY-07 (respectively -0.11 and -0.16 gC.m<sup>-2</sup>.d<sup>-1</sup>) as well as the slopes associated to these data (respectively 0.75 and 0.65 for CY-07 and 1.22 and 1.00 for CY-16). Nevertheless, when enough images are available to describe the crop dynamics (notably key dates such as those corresponding to emergence, maximum vegetation, start of senescence, harvest), the crop production as well as the CO<sub>2</sub> and water fluxes are estimated with acceptable accuracy (see (Pique et al., 2020a)).

<figure><img src="image_55.png" /></figure>

**Figure II-10 :** Scatter plots of the observed and simulated GAI (A), DAM (B), GPP (C), R<sub>ECO</sub> (D), NEE (E) and ETR (F) for the two crop years of sunflower at FR-Aur and mean statistic over the 2 crop years.

GPP and R<sub>ECO</sub> are estimated with relatively low errors (RMSE of 2.36 gC.m<sup>-2</sup>.d<sup>-1</sup> and 0.96 gC.m<sup>-2</sup>.d<sup>-1</sup>, respectively and rRMSE of 39.5 and 40.9 % respectively) and good correlations (R<sup>2</sup> of 0.81 and 0.83, respectively).
In the proposed approach, the errors associated with GPP and R<sub>ECO</sub> compensate for each other when estimating NEE (which is calculated from the difference between the two variables). This resulted in an accurate estimation of the NEE (RMSE of 0.97 gC.m<sup>-2</sup>.d<sup>-1</sup> and R<sup>2</sup> of 0.83) despite the lack of images during the 2016 senescence phase.

The model is also able to reproduce the ETR (Ev+TR) dynamics of the two crop years with good accuracy (RMSE of 0.6 mm.d⁻¹, rRMSE = 40 % and R² of 0.79). For the CO₂ fluxes, the ETR is underestimated by the model for CY-16 during the senescence period, which led to a decrease in TR before the end of the observed crop cycle (Figure II-10 F).

In the proposed approach, R<sub>H</sub> depends only on the soil temperature, while it is known to be impacted by many other factors (e.g. AWC and soil organic C content). In their study, (Pique et al., 2020a) discussed taking into account AWC in the estimation of R<sub>H</sub>. We thus estimated R<sub>H</sub> according to the (Delogu, 2013-eq. III.5) that accounts for AWC. During the two simulated crop years at FR-Aur (data not shown), the temperature-dependent-only formulation performed better (RMSE = 0.5 gC.m<sup>-2</sup>.d<sup>-1</sup> and R<sup>2</sup> = 0.36) than the combined temperature- and AWC-dependent equation (RMSE = 0.73 gC.m<sup>-2</sup>.d<sup>-1</sup> and R<sup>2</sup> = 0.06), even when considering *in situ* AWC to estimate R<sub>H</sub> (RMSE = 0.61 gC.m<sup>-2</sup>.d<sup>-1</sup> and R<sup>2</sup> = 0.15). We therefore kept the first-order exponential equation that depends only on soil temperature for the rest of the study.

### 4.5.2. Large scale application

#### 4.5.2.1. Impact of soil data and of coupling with the water module on model performance

The two sunflower crop years at FR-Aur were simulated with the SAFYE-CO₂ model using field capacity, the permanent wilting point and soil depth estimated from i) the *in situ* measurements, ii) the GSM products, and iii) the SG products. These two crop years were also simulated with iv) the SAFY-CO₂ model to evaluate the benefits of coupling with the water module compared to the previous version. The different evaluated configurations of the model are referred to as the following: i) SAFYE-CO₂<sup>IS</sup>, ii) SAFYE-CO₂<sup>GSM</sup>, iii) SAFYE-CO₂<sup>SG</sup> and iv) SAFY-CO₂.

Figure II-11 summarizes the statistics (RMSE and R²) obtained when testing the different configurations of the model in simulating the target variables (GAI, DAM, GPP, R<sub>ECO</sub>, NEE and ETR) for CY-07 and CY-16. There are no statistical data on ETR for SAFY-CO₂, as this version does not simulate water flux. The differences observed between the different tested configurations are mainly due to water stress. In addition to simulating water fluxes (Ev, TR), the water module impacts vegetation photosynthesis and growth through a water stress function (WtS<sub>V</sub>). Therefore, correctly estimating AWC is essential for accurate simulation of crop growth.

Regarding CY-07, the GAI and R<sub>ECO</sub> are estimated with nearly the same accuracy for each model configuration. Concerning the GAI, the good performance of each model configuration indicates that all the configurations were able to accurately fit the remotely sensed GAI regardless of the soil data used. The errors in R<sub>ECO</sub> are also similar between the different configurations for two reasons: i) the heterotrophic component of R<sub>ECO</sub> is estimated in the same way for all versions, and ii) the autotrophic component, depending on the GPP, differs between versions only at the end of the crop cycle when water stress occurs. However, greater differences occur for the DAM, GPP, NEE and ETR estimates, where SAFYE-CO₂<sup>GSM</sup> and SAFYE-CO₂<sup>SG</sup> show low accuracy. This is due to the overestimation of the water stress at the end of the season, which in turn underestimates the water fluxes and GPP and thus the DAM. Indeed, TWC is underestimated for both SAFYE-CO₂<sup>GSM</sup> and
SAFYE-CO₂SG during this period (see section II-4.3.2.4), leading to an underestimation of AWC, which induces stress at the end of the cropping season.

<figure><img src="image_56.png" /></figure>

**Figure II-11 :** RMSE and R² values for the main model's output for CY-07 (a) and CY-16 (b) at FR-Aur via SAFYE-CO₂ (orange bars) or SAFYE-CO₂ run with *in situ* soil data (SAFYE-CO₂IS, green bars), GlobalSoilMap data (SAFYE-CO₂GSM, red bars) and SoilGrids data (SAFYE-CO₂SG, blue bars).

For CY-16, the same observation can be made concerning the estimations of the GAI. Each version is able to accurately reproduce GAI dynamics. However, SAFYE-CO₂IS largely outperformed SAFYE-CO₂GSM and SAFYE-CO₂SG for all the other variables. SAFYE-CO₂SG and SAFYE-CO₂GSM show high errors when estimating GPP and ETR because, like CY-07, those versions largely overestimated water stress, inducing early senescence and thus an underestimation of the CO₂ and water fluxes.

Despite the observed water stress during both simulated years (at the end of the crop cycle), the best performances in estimating CO₂ fluxes were achieved by the SAFYE-CO₂ model that did not simulate water stress, indicating that the water stress is already taken into account through the
remotely sensed GAI. Additionally, the water module adds an additional level of complexity, while
the performance of SAFY-CO₂ is already very high.

Figure II-12 shows the components of the NECB estimated with the different configurations of the model compared with observed values: the net annually cumulated CO₂ fluxes [net annual ecosystem productivity (NEP)], the C input as seeds (C<sub>INP</sub>) and the C exported at harvest (C<sub>EXP</sub>). Because of the poor estimates of TWC deduced from GSM and SG products, SAFYE-CO₂<sup>GSM</sup> and SAFYE-CO₂<sup>SG</sup> overestimate NEP and underestimate C<sub>EXP</sub> for both simulated years. Both models' configurations overestimate water stress, especially at the end of the season when AWC is low, which leads to premature senescence and thus an early harvest. In the absence of vegetation, the only CO₂ flux is R<sub>H</sub> (flux towards the atmosphere), so shortening of the crop cycle by the two model configurations leads to an overestimation of NEP. The same reason explains the underestimation of C<sub>EXP</sub>. The more stressed the vegetation is, the lower the total biomass production and the yield are. On the other end, both SAFYE-CO₂<sup>IS</sup> and SAFY-CO₂ were able to accurately reproduce the NECB for CY-07 (-4.3 and -10.8% for SAFYE-CO₂<sup>IS</sup> and SAFY-CO₂, respectively) and for CY-16 (+15 and -6.2% for SAFYE-CO₂<sup>IS</sup> and SAFYE-CO₂, respectively), as both configurations were able to accurately simulate the components of the C budget.

<figure><img src="image_57.png" /><figcaption>Figure II-12 : Components of the net ecosystem carbon budget for CY-07 (a) and CY-16 (b) at FR-Aur for the SAFYE-CO2IS (green bars), SAFYE-CO2GSM (red bars), SAFYE-CO2SG (blue bars) and SAFY-CO2 (orange bars) models compared to observed value (black bars).</figcaption></figure>

4.5.2.2. Biomass and yield estimates

Apart from the two sunflower crop years available at FR-Aur, the model was evaluated for DAM and YLD over numerous fields in 2013, 2014 and 2015 (see section II-4.3.2.2) to assess its ability to reproduce those variables under a variety of climatic, soil and management conditions. Since no soil *in situ* measurements were available at the plot scale in 2013, 2014 and 2015, the SAFYE-CO₂<sup>IS</sup> configuration could not be tested.

Figure II-13 presents DAM and YLD simulated by SAFYE-CO₂GSM (A-B), SAFYE-CO₂SG (C-D) and SAFY-CO₂ (E-F). Although SAFYE-CO₂GSM and SAFYE-CO₂SG showed good overall performance at estimating DAM over the whole dataset (R² of 0.78 and 0.74, respectively; RMSE of 110 and 119 g.m⁻², respectively), these two models’ configurations had difficulty accurately estimating the final DAM for the 2013 field campaign. This could be explained by water stress simulated at the end of the season, which slows down crop growth and thus underestimates final biomass. This inaccurate estimation of water stress was the result of poor estimation of the TWC. In the previous sections, TWC derived from GSM and SG were compared to the value derived from *in situ* measurements performed at FR-Aur (Table II-3). The values associated with the global soil databases
were lower than those observed in *situ*, and only a few variations were observed in terms of texture (particularly for SG, Figure 2). The TWC distributions derived from GSM and SG for all the studied fields are presented in Figure II-14 and were compared to values derived from *in situ* measurements collected nearby (i.e. 45 fields located in the northwestern part of the study area—experiment described by (Pique et al., 2020a)). The values of TWC derived from the global soil databases are low compared to those obtained from the previous studied fields, and very little spatial variation is present. The TWC derived from GSM ranges between 0.11 and 0.14 m³.m⁻³, while the TWC derived from SG is even less variable—between 0.12 and 0.13 m³.m⁻³. Considering the *in situ* measurements, soils in the study area present more variation, with the TWC between 0.10 and 0.16 m³.m⁻³. These low and quite spatially constant TWC values derived from the global soil databases undoubtedly lead to an overestimation of water stress and thus to an underestimation of the final DAM. Nevertheless, SAFY-CO₂, which does not simulate water flux or water stress, better estimated DAM (R² of 0.91 and RMSE of 73 g.m⁻²) and especially better reproduced the final highest DAM values (red dots in Figure 8 E), confirming that SAFYE-CO₂GSM and SAFYE-CO₂SG simulated too much water stress in 2013.

<figure><img src="image_58.png" /><figcaption>Figure II-13 : Scatter plots of simulated and observed DAM (left) and YLD (right) obtained with SAFYE-CO₂GSM (A-B), SAFYE-CO₂SG (C-D) and SAFY-CO₂ (E-F) model.</figcaption></figure>

Regarding YLD estimates (Figure II-13 B-D-F), all versions have difficulty to reproduce the observed values accurately. This is possibly due to several factors. First, YLD is calculated from the maximum value of DAM using a constant HI. The HI of sunflower is known to vary with the length of the grain filling period and with temperature (Bange et al., 1998) but also with water stress (Ghelinezhad et al., 2009); moreover, there is some additional variations due to crop variety (Debaeke et al., 2010). These factors are not yet accounted for in the proposed approach. Additionally, the
observed yields were provided by the farmers, and may have by high uncertainty. Indeed, compared to
the “true” yield, those data are subject to significant grain losses at harvest.

<figure><img src="image_59.png" /></figure>

**Figure II-14 :** TWC of the fields simulated in this study derived from GSM and SG and of 45 fields within the study area estimated with *in situ* (IS) measurements.

**4.5.3. C budget assessment at the landscape scale**

Since the SAFYE-CO₂ model has difficulty reproducing sunflower growth if no accurate
information on soil TWC is available (i.e. SAFYE-CO₂GSM and SAFYE-CO₂SG, which is the case at
the landscape scale), the SAFY-CO₂ model was used to simulate the components of the C budget of
CY-15 at the landscape scale (Figure II-15) for more than 200 fields (black dashed square in Figure II-
6). Figure II-15 presents the corresponding NEP (A), CEXP (B) and resulting NECB (C).

<figure><img src="image_60.png" /></figure>

**Figure II-15 :** The annual carbon budget (NECB) and its component, NEP and Cexp, simulated for 200 fields with SAFY-CO₂ inside the study area for CY-15. Digital elevation model from *Shuttle Radar Topography Mission* is displayed in the background.

Considering all 200 fields, the simulated NEP is between -91 (sink) and 232 (source) gC.m⁻².yr⁻¹, while the amount of C exported at harvest is between 68 and 155 gC.m⁻².yr⁻¹. This results in a NECB between 64 and 300 gC.m⁻².yr⁻¹, which means that each simulated field behaves as a C source whose flux is towards the atmosphere. The topography seems to slightly affect sunflower development. This tendency can be seen in Figure II-15 A, plots where the C assimilation is higher (lowest NEP) are located mainly in plains and valleys. The same observation can be made for the CEXP term, meaning that the highest yields are observed in plains or in the lowest parts of hillsides. However, high yields are also observed in hilly areas, and poor yields are also observed in plains. Concerning the NECB, the tendency is less clear as the C budget is conditioned by the NEP and the CEXP term, whose signs may be opposite, with each favorably affected by good soil conditions (not short of water). Additionally, NEP can be affected by spontaneous regrowth events or weeds, increasing the CO₂ net assimilation, albeit with no effect on yield. Differences in the availability of water resources for crop management can also explain the spatial variability in the NECB and its components. All plots simulated here are sources of C because the length of the sunflower cycle is short, which leaves a long fallow period that
releases C from R<sub>H</sub> to the atmosphere. This observation is consistent with the results of a meta-analysis conducted for more than 15 flux crop sites in Europe by (Ceschia et al., 2010a).

**4.6. Discussion**

**4.6.1. SAFYE-CO₂ overall performances**

The SAFYE-CO₂ model for sunflower was validated in this study. Section 4.1 presents the results achieved with this approach using *in situ* soil measurements to calculate SWC parameters and soil depth. The performances for simulating the DAM and yield are within the range of those of other studies. As an example, (Todorovic et al., 2009a) compared the performances of three crop models (AquaCrop, CropSyst and WOFOST) for simulating aboveground biomass and yield of sunflower. Those authors showed that under rainfed conditions (similar to the conditions in the present study), the accuracy of the models ranged from 63 (WOFOST) to 204 (AquaCrop) g.m⁻², while SAFYE-CO₂ estimated DAM at FR-Aur with an error of 66 g.m⁻². Concerning yield, SAFYE-CO₂ (using *in situ* measurements to estimate the TWC) provided estimates with an accuracy (error of 54 g.m⁻²) similar to that shown by (Todorovic et al., 2009a), who reported errors ranging from 70 (AquaCrop) to 94 g.m⁻² (CropSyst). However, the performances for estimating yield with SAFYE-CO₂ should be taken with caution because of the small number of comparisons with field observations (only 2 corresponding to the Auradé site in 2007 and 2016), meaning that additional validation exercises are needed.

CO₂ fluxes are also robustly estimated by the SAFYE-CO₂ model and are within the range of other models. (Lehuger et al., 2010) evaluated the CERES-EGC model against CO₂ fluxes over several European flux sites, including FR-Aur. The model reproduced the net CO₂ fluxes over the growing period of CY-07 with an error of 30 gC.m⁻², while SAFYE-CO₂<sup>IS</sup> (using *in situ* measurements for soil description) estimated NEP (throughout the whole cropping year) with an error of -12.1 gC.m⁻². Additionally, in 2016, despite the lack of satellite images covering the period of senescence at FR-Aur, the model provided accurate estimates of NEE (mainly because the errors in GPP and R<sub>ECO</sub> compensate for each other). This accuracy in estimating CO₂ fluxes makes possible the realistic estimation of the NECB at the plot scale with the proposed approach. SAFYE-CO₂<sup>IS</sup> was able to reproduce NECB with very good accuracy (relative errors of -4.3% and +15% for CY-07 and CY-16, respectively). SAFYE-CO₂<sup>IS</sup> also performed well at estimating ETR. During the two simulated crop years at FR-Aur, the model reproduced ETR with a low error (RSME = 0.68 mm.d⁻¹), meaning that our approach can also be used to calculate annual water budgets for sunflower plots.

At the landscape scale, coupling with a soil water module tended to lower biomass and yield estimates. The poor estimates of TWC derived from soil products (GSM and SG) led to an underestimation of the final biomass and, thus, yield. In contrast, the biomass was accurately estimated by the SAFY-CO₂ model (R² = 0.91 and RMSE = 73 g.m⁻²), and the estimations were within the range of those of models specific to sunflower: Using the OILCROP-SUN model, (Villalobos et al., 1996) estimated DAM with an error of 170 g.m⁻², while using the SUNFLO model, (Casadebaig et al., 2011a) simulated DAM with an error of 114 g.m⁻². In contrast, the yield was poorly estimated by SAFY-CO₂ in terms of its correlation (R² = 0.31) but the error was still within the range or better than those other studies. The error achieved by SAFY-CO₂ in estimating YLD was 51 g.m⁻², while errors with OILCROP-SUN and SUNFLO were 80 and 50 g.m⁻², respectively (in (Villalobos et al., 1996) and (Casadebaig et al., 2011a)). The poor results achieved for yield with SAFY-CO₂ can be explained in part by differences in the HI between sunflower plots. Indeed, sunflower HI is known to change with genotype, management and environment ((Bange et al., 1998; Gholinezhad et al., 2009)) and depends on many factors that are not yet accounted for in the proposed approach. Accounting for
varietal potentials (Debaeke et al., 2010) and for abiotic stress after flowering as in CropSyst (Stöckle et al., 2003b) is an avenue for improvement with our approach.

### 4.6.2. Prerequisite to compute water fluxes over wide area

Estimating water fluxes and thus water stress over croplands requires an accurate description of soil properties ($\theta_{FC}$, $\theta_{WP}$ and soil depth). This need has been highlighted in the present study. Soil parameters were derived from *in situ* measurements when available or from soil texture products for upscaling. The two existing soil products tested in this study (GSM and SG) did not meet the needs for the expected accuracy. Indeed, we demonstrated (see section II-4.5.2) that the use of soil products led to an underestimation of the soil TWC, which in turn led to an overestimation of the water stress and an underestimation of the biomass and yield compared to the *in situ* data. However, this conclusion concerns only our study area, which is not representative of the spatial scale of GSM and SG products. Indeed, those products provide, among other things, soil texture and depth at a fine resolution (SG: 250 m; GSM: 90 m) all across the globe. As such, it would be difficult to expect those products to accurately describe differences in soil properties between plots over a small area of study. Additionally, taking water fluxes into account in our approach does not improve the estimation of biomass or CO₂ fluxes, even under water stressed conditions (see section II-4.5.2.2). Both sunflower crop years at FR-Aur experienced water stress (especially at the end of the season), and SAFY-CO₂ was still able to accurately reproduce the crop dynamics and production as well as the CO₂ flux, indicating that the plant stresses (i.e. water and nitrogen) were already correctly taken into account through the remotely sensed GAI. Vegetation stress slows down (or even stops) growth, which is observable through GAI. Therefore, we believe that coupling SAFY-CO₂ with a soil water module is useful for estimating ETR and water budgets or for assessing crop water-use efficiency (Tallec et al., 2013a) when safeguarding soil and water resources, but this coupling does not seem to improve the simulations of the annual C budget components. During the two sunflower crop years at FR-Aur, simulations were performed using SAFYE-CO₂, with no feedback of the SWC on GPP computation (WtS<sub>V</sub> = 1, eq. II-20). The accuracy in estimating biomass and CO₂ flux was similar to that of SAFY-CO₂; thus, the accuracy was slightly better than that of simulations considering the water stress function. Concerning the estimation of water fluxes, the same accuracy was obtained with or without the consideration of effects of water stress on vegetation development. We thus argue that SAFYE-CO₂ should be used in a decoupled way between vegetation growth and SWC since information about water stress is already provided by the remotely sensed GAI. However, further analysis that includes additional water flux data and the simulation of other crop species should reveal whether the water module should be used systematically in a decoupled way, allowing us to estimate water fluxes in the absence of vegetation stress.

### 4.6.3. Limitations and potential improvements

In addition to simulating crop dynamics and production, our approach fills the gap in assessing the C and water budget simultaneously at the plot scale over large areas. Notably, to achieve this, information on both organic fertilization and eventually irrigation at the plot scale are needed. This information is not yet available at a landscape spatial scale, which can lead, for example, to great uncertainties in the C budget calculations in areas with livestock farming. However, this drawback could be compensated by the combination of our approach with the use of data in farm management information systems (FMIS). Additionally, the representation of R<sub>H</sub> is quite simplistic in our approach and would require improvement. As shown in the present study and as discussed by (Pique et al., 2020a), this simple method allows for the accurate assessment of the annual C budget, but when