Gaétan Pique
email: piqueg@cesbio.cnes.fr

Rémy Fieuzal

Al Bitar Ahmad

Amanda Veloso

Tiphaine Tallec

Aurore Brut

Morgan Ferlicoq

Bartosz Zawilski

Jean-François Dejoux

Hervé Gibrin

Pique Gaétan

Fieuzal Rémy

Ahmad Al Bitar

Zawilski Bartosz

Gibrin Hervé

Ceschia Eric

Estimation of daily CO2 fluxes and of the components of the carbon budget for winter wheat by the assimilation of Sentinel 2-like remote sensing data into a crop model

Keywords: crop modelling, carbon budget, remote sensing monitoring, winter wheat

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Agricultural lands occupy nearly 12% of Earth's terrestrial surface. They not only contribute to but also affect climate change because climatic conditions and water resources affect crop production [START_REF] Smith | Carbon sequestration potential in European croplands has been overestimated[END_REF].

Additionally, the global food demand is increasing and may continue to increase for decades, driven by the increasing global population and per capita income that are anticipated through the middle of the next century [START_REF] Tilman | Global food demand and the sustainable intensification of agriculture[END_REF].

It is in this context that the '4 per mille Soils for Food Security and Climate' initiative was launched at COP21, with the aspiration to increase global soil organic matter stocks by 4 per 1000 (or 0.4%) per year as a compensation for part of the global emissions of greenhouse gases by anthropogenic sources and to increase food security [START_REF] Chabbi | Aligning agriculture and climate policy[END_REF][START_REF] Minasny | Soil carbon 4 per mille[END_REF]. Since then, this initiative has induced a wide debate in the scientific community concerning its feasibility [START_REF] Baveye | The "4 per 1000" initiative: A credibility issue for the soil science community?[END_REF][START_REF] Poulton | Major limitations to achieving "4 per 1000" increases in soil organic carbon stock in temperate regions: Evidence from long-term experiments at Rothamsted Research, United Kingdom[END_REF], and it has been recognized that such an increase in soil organic carbon (SOC) is likely achievable in soils that are being actively managed for agriculture at a rate of increase that may not be achievable everywhere [START_REF] Chabbi | Aligning agriculture and climate policy[END_REF][START_REF] Lal | Beyond COP 21: Potential and challenges of the "4 per Thousand" initiative[END_REF][START_REF] Minasny | Soil carbon 4 per mille[END_REF][START_REF] Pellerin | Stocker du carbone dans les sols français, quel potentiel au regard de l'objectif 4 pour 1000 et à quel coût[END_REF]. This debate illustrates the need for tools that can estimate changes in cropland SOC and identify potential levers to increase it. Currently, quantifying the net ecosystem carbon budgets (NECB) of croplands at regional or global scales remains difficult because of the heterogeneous character of agricultural landscapes, which have numerous plots with varied management practices and environmental conditions. This character results in uncertainties when assessing the impacts of specific management practices on the cropland NECB [START_REF] Osborne | Key questions and uncertainties associated with the assessment of the cropland greenhouse gas balance[END_REF] and when determining whether croplands are carbon sinks or sources [START_REF] Ciais | The European carbon balance. Part 2: croplands[END_REF][START_REF] West | Cropland carbon fluxes in the United States: increasing geospatial resolution of inventory-based carbon accounting[END_REF]. Indeed, the general biogeochemical models (such as SPA [START_REF] Williams | Modelling the soil-plant-atmosphere continuum in a Quercus-Acer stand at Harvard Forest: the regulation of stomatal conductance by light, nitrogen and soil/plant hydraulic properties[END_REF], Ecosys [START_REF] Grant | Net Biome Productivity of Irrigated and Rainfed Maize-Soybean Rotations: Modeling vs[END_REF], Isba-Ags [START_REF] Calvet | An interactive vegetation SVAT model tested against data from six contrasting sites[END_REF], ORCHIDEE [START_REF] Krinner | A dynamic global vegetation model for studies of the coupled atmospherebiosphere system[END_REF], and ORCHIDEE-STICS [START_REF] Gervois | Carbon and water balance of European croplands throughout the 20th century: CARBON BALANCE OF EUROPEAN CROPLANDS[END_REF]) that are commonly used to simulate the carbon cycle in terrestrial ecosystems are not suited to account for the specificities and complexities of agro-ecosystems, particularly the effect of management practices. In contrast, the agronomic models or so-called crop models (e.g., CERES [START_REF] Gabrielle | A model of leaf area development and senescence for winter oilseed rape[END_REF] or STICS [START_REF] Brisson | STICS: a generic model for the simulation of crops and their water and nitrogen balances. I. Theory and parameterization applied to wheat and corn[END_REF]) that account for management and pedoclimatic effects are primarily designed for simulating crop development and production (net primary production (NPP), yield) at the plot scale. However, unlike our approach, these models require information regarding management practices, which makes them less suitable for large spatial scale applications.

Several studies have demonstrated the benefit of assimilating remote sensing data into regional-scale crop models [START_REF] Sus | A data assimilation framework for constraining upscaled cropland carbon flux seasonality and biometry with MODIS[END_REF][START_REF] Wu | Assimilating MODIS-LAI into Crop Growth Model with EnKF to Predict Regional Crop Yield[END_REF]. In particular, the combination of high spatial and temporal resolution (HSTR) remote sensing data with crop models can provide, at the field scale over large areas, a timely and accurate picture of crop development [START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF][START_REF] Hadria | Potentiality of optical and radar satellite data at high spatio-temporal resolutions for the monitoring of irrigated wheat crops in Morocco[END_REF], cropland photosynthesis [START_REF] Wang | Vegetation primary production estimation at maize and alpine meadow over the Heihe River Basin, China[END_REF][START_REF] Wolanin | Estimating crop primary productivity with Sentinel-2 and Landsat 8 using machine learning methods trained with radiative transfer simulations[END_REF] and net CO2 fluxes [START_REF] Revill | Carbon cycling of European croplands: A framework for the assimilation of optical and microwave Earth observation data[END_REF][START_REF] Sus | A data assimilation framework for constraining upscaled cropland carbon flux seasonality and biometry with MODIS[END_REF]. Among others, the SAFY [START_REF] Duchemin | A simple algorithm for yield estimates: Evaluation for semi-arid irrigated winter wheat monitored with green leaf area index[END_REF][START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF] and SAFY-WB [START_REF] Battude | Modeling water needs and total irrigation depths of maize crop in the south west of France using high spatial and temporal resolution satellite imagery[END_REF][START_REF] Duchemin | Impact of Sowing Date on Yield and Water Use Efficiency of Wheat Analyzed through Spatial Modeling and FORMOSAT-2 Images[END_REF] crop models constitute coherent frameworks for estimating biomass, yield production and water requirements. These models describe the main biophysical processes underlying crop production by using climatic data and assimilate green area index (GAI) dynamic maps derived from remote sensing to avoid the need for management data, which makes them well suited for large-scale studies. In this work, we modified the SAFY model to simulate the components of the net ecosystem exchange (NEE) and to evaluate the potential of this approach for calculating cropland annual carbon budgets. The resulting model, called SAFY-CO2, is described and evaluated against in situ data. The objectives of this study are as follows:

1-To assess the potential of an approach combining HSTR remote sensing data and a simple crop model to quantify the components of the NEE and of the annual NECB for winter wheat plots in contrasting climatic and management conditions.

2-To address the potentialities and limitations of such an approach in the perspective of future regional-or global-scale applications.

To fulfil our objectives, GAI maps derived from HSTR optical data (Formosat-2 and SPOT satellites) from 2006 to 2014 in southwestern France were used to constrain the photosynthetic light-use efficiency and phenological parameters of the model. Consequently, the simulated crop phenology agreed well with the satellite observations, which is essential for correctly estimating CO2 fluxes and carbon budgets [START_REF] Grant | Net Biome Productivity of Irrigated and Rainfed Maize-Soybean Rotations: Modeling vs[END_REF][START_REF] Huang | Agro-C: A biogeophysical model for simulating the carbon budget of agroecosystems[END_REF][START_REF] Wattenbach | The carbon balance of European croplands: A cross-site comparison of simulation models[END_REF].

A validation of the simulated CO2 fluxes (photosynthesis, ecosystem respiration and NEE) was performed against eddy-covariance flux measurements that were carried out over two flux sites [START_REF] Béziat | Effet des conditions environnementales et des pratiques culturales sur les flux de carbone et d'eau dans les agrosystèmes[END_REF] in the study area.

In the next section of the manuscript, the study area, experimental datasets and satellite database are presented.

The following section describes the SAFY-CO2 mathematical formulations, the parameterization and calibration procedure and the method for computing the annual NECB. The inputs required to run the model as well as the validation procedure are also detailed. Section 4 is dedicated to the results. The biomass and yield results are presented first, followed by the flux estimate results and finally the annual C budgets. Section five discusses the potentialities, the limitations and potential improvements of such an approach in the perspective of future regional-or global-scale applications. The paper ends with a conclusion concerning the main results, limitations and insights into future developments.

MATERIALS

Study area

The study area is part of the Regional Spatial Observatory [http://www.cesbio.ups-tlse.fr/fr/osr.html] located next to Toulouse in southwest France which includes 2 instrumented agricultural sites, Auradé (FR-Aur) and Lamasquère (FR-Lam) (Figure 1). Those two sites belong to the Integrated Carbon Observation System (ICOS) network [https://www.icos-ri.eu /] for observations of surface fluxes (CO2, latent and sensible heat fluxes). The region has a temperate climate, with an annual mean precipitation of approximately 655 mm and an annual mean temperature of 12.9°C (measured by Meteo France at the Toulouse-Blagnac station between 1961 and 1990; see http://www.infoclimat.fr/climatologie/index.php). Agricultural activity occupies almost 90% of the landscape, and winter wheat is the main cultivated crop (covering approximately 20% of the total surface area). Sown from mid-October until the beginning of December, winter wheat is harvested from mid-June until the end of July, and straw is usually incorporated into the soil.

Soil characteristics

The nature of the soils of the study area is shaped by measurements are presented within the USDA triangle (United States Department of Agriculture, Figure 2).

With fractions between 9 and 50% for the clay, between 25 and 72% for the silt and between 12 and 55% for the sand, the observed contents cover wide ranges of each component. On average, the texture is composed of 48% of silt and 26% of clay and sand, illustrating the dominance of silt fraction within the study area. FR-Lam soil is more clayey than the FR-Aur one and thus less permeable. As FR-Lam is located on the terraces near a river flood can occur after heavy rainfall. As the FR-Aur site is located on a hillside, its soil is heterogeneous and its depth vary from 1 to more than 2 m.

In situ data

The FR- , 2005-2006[START_REF] Béziat | Effet des conditions environnementales et des pratiques culturales sur les flux de carbone et d'eau dans les agrosystèmes[END_REF]-2010, 2011-2012and 2013[START_REF] Veloso | Modélisation spatialisée de la production, des flux et des bilans de carbone et d'eau des cultures de blé à l'aide de données de télédétection : application au sud-ouest de la France[END_REF]at FR-Aur and 2006-2007, 2008[START_REF] Béziat | Effet des conditions environnementales et des pratiques culturales sur les flux de carbone et d'eau dans les agrosystèmes[END_REF], 2010-2011and 2012[START_REF] Delogu | Modélis ation de la respiration du sol dans les agrosystèmes[END_REF] at FR-Lam.

To facilitate the reading of this paper, we will identify each site-year by the three first letters of the site followed by the harvest year (e.g., AUR2006 for the site-year 2005-2006 at Auradé).

The FR-Aur field (23.5 ha) is located on a hillside area near the Garonne River terraces and is characterized by a rapeseed/winter-wheat/sunflower/winter-wheat rainfed rotation that only receives mineral fertilizers. Only the grain is exported. The FR-Lam field (23.8 ha) is part of an experimental farm for milk and chicken production owned by the Purpan engineering school EIP (Ecole d'Ingénieurs de Purpan). It is characterized by a maize/winter-wheat rotation that is used to feed livestock and provide litter. Therefore, nearly all aboveground biomass is exported as grain and straw for winter wheat, and irrigated maize is harvested when it is still green for silage. Both organic and mineral fertilizers are applied.

A field campaign was conducted in June-July 2011 in 16 winter wheat fields to obtain spatially distributed in situ biomass and yield data (referred to as the '2011 field campaign'). The selection of fields was based on the analysis of the intra-and inter-field variability of the Normalized Difference Vegetation Index (NDVI derived from the Formosat-2 and SPOT optical images of April) to monitor a wide range of vegetation development.

Crop biomass and yield measurements were performed just before the harvest. During the vegetative cycle, the crop development at the experimental sites was regularly monitored using destructive measurements of GAI and dry aboveground mass (DAM) (see [START_REF] Béziat | Effet des conditions environnementales et des pratiques culturales sur les flux de carbone et d'eau dans les agrosystèmes[END_REF] for protocol), while yield data were provided by the farmers that cultivate the two flux sites. Farmer's data are often a mean of the yield at several fields surrounding the instrumented sites.

GAI, biomass and yield data

During the 2011 field campaign, biomass and yield data were collected from 16 fields according to the VALERI sampling protocol [http://w3.avignon.inra.fr/valeri/]. The samples were collected from five homogeneous square subplots inside the 16 Elementary Sampling Units (ESUs) of 20×20 m 2 . The subplots sampled in each ESU were located in the ESU centre and corners. For each subplot, 4 rows with lengths of 50 cm and an inter-row distance of 13 cm were collected, which resulted in a sampling surface area of 0.25 m². The samples were dried and weighed, and the destructive grain yield was measured after threshing. The mean and associated standard deviations of the five subplots were calculated for each ESU.

Table 1 summarizes, for each year, the number and the kind of data that were acquired and it specifies if spontaneous regrowth, cover crop or weed development occurred or not after harvest for each year.

Flux and meteorological measurements

Turbulent fluxes of CO2, water vapor (evapotranspiration and latent heat), sensible heat and momentum were measured continuously using the eddy-covariance (EC) method [START_REF] Aubinet | Estimates of the Annual Net Carbon and Water Exchange of Forests: The EUROFLUX Methodology[END_REF][START_REF] Baldocchi | Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future[END_REF][START_REF] Moncrieff | A system to measure surface fluxes of momentum, sensible heat, water vapour and carbon dioxide[END_REF]. The EdiRe software (Robert Clement, © 1999, University of Edinburgh, UK) was used to calculate the turbulent fluxes. The NEE was calculated as the sum of turbulent CO2 fluxes and changes in CO2 storage under the EC devices. Flux filtering, quality controls and gap filling were performed following the CarboEurope-IP recommendations ([www.carboeurope.org], see [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF] for more details). The NEE was partitioned into gross primary production (GPP) and ecosystem respiration (Reco) components according to the method proposed by [START_REF] Reichstein | On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm[END_REF] and adapted by [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF] for croplands (a process that could lead to over-or underestimations of the two components of the NEE). During the periods of bare soil, the GPP was set to 0, and the measured NEE fluxes only represented the Reco component. Finally, the net ecosystem production (NEP) was derived from the annual integration of the NEE values. Synchronously, the standard meteorological variables were recorded at each experimental site and included different radiation components (i.e., direct and diffuse components of incoming global radiation). After pre-processing, the semihourly fluxes and meteorological data were integrated or averaged at a daily time scale to be consistent with the model time step. Note that because of instrument failure, there was a 3-month gap in the flux measurements at the beginning of 2011 at FR-Lam.

For the '2011 field campaign' simulations, the SAFRAN meteorological data produced by Météo-France [START_REF] Durand | A meteorological estimation of relevant parameters for snow models[END_REF] are used. The SAFRAN data provide the air temperature, incoming global radiation, precipitation, and relative air humidity 2 m above the ground and the wind speed 10 m above the ground based on weather station measurements and modelling. The data are available every 6 h over an 8 km spatial resolution grid. The daily means of these climatic variables are calculated for each 2011 campaign field using bilinear interpolation.

The other simulations are performed using the climatic data recorded at the instrumented sites.

During the studied years, the climatic conditions were very contrasted. Finally, the cropping year 2013-2014 was characterized by its warm winter (8.1°C against 6.6°C on average) allowing the crop to start early in the season.

Satellite data and products

Multi-satellite optical images

This study uses an extensive dataset of HSTR from several satellites. Because of the spatial and temporal resolution of this dataset, and also as the bands necessary for this study are available from Sentinel-2 we consider that their combined used in this modelling exercise is representative of what could be achieved with Sentinel 2. the optical domain (multispectral mode). The F2 images were characterized by a spatial resolution of 8 m (footprint of 24×24 km 2 ) and were acquired at the same viewing angle (±45°) in four narrow wavelengths (blue, green, red and near-infrared) [START_REF] Chern | Taiwan's second remote sensing satellite[END_REF]. The SPOT images were characterized by spatial resolutions of 20 m (Spot-2/4) and 10 m (Spot-5) covering an area of 60 × 60 km 2 . The images were acquired at two incidence angles (75° and 102°) in at least three wavelengths (green, red, and near-infrared), with the medium-infrared wavelength for SPOT-4/5 [START_REF] Arnaud | Performance of the APSIM-wheat model in Western Australia[END_REF].

Surface reflectances were derived from the satellite data using the KALIDEOS processing chain [http://kalideos.cnes.fr] for atmospheric, radiometric and geometric corrections. The mean geometric correction accuracy was close to 0.2 pixels (LaFrance, Lenot, Ruffel, Cao, & Rabaute, 2012), which is satisfactory for the surface area of the studied fields.

The combined use of images acquired by different satellites was important to increase the number observations per cropping year (see Figure 3). The comparison of the reflectances or of the NDVIs derived from different sensors, including Formosat-2 and SPOT, acquired at close dates over various crops, had highlighted the good performances of the processing chain and the limited effect of the sensor type [START_REF] Battude | Estimating maize biomass and yield over large areas using high spatial and temporal resolution Sentinel-2 like remote sensing data[END_REF][START_REF] Fieuzal | Estimation of corn yield using multi-temporal optical and radar satellite data and artificial neural networks[END_REF].

From image reflectance to GAI estimates

The seasonal dynamic maps of GAI were derived from the surface reflectances using the BV-NNET tool (Biophysical Variables Neural NETwork, [START_REF] Baret | LAI, fAPAR and fCover CYCLOPES global products derived from VEGETATION[END_REF]), which consists of a trained artificial neural network (ANN) using the outputs of a radiative transfer model (PROSAIL Jacquemoud et al., (2009)). ANNs were first trained with the wide range of conditions estimated by the radiative transfer model. Then, the trained network was used to predict the GAI from satellite reflectances. The GAI estimates derived from F2 and SPOT reflectances were compared to non-destructive measurements based on digital hemispherical photographs collected over a range of crops [START_REF] Demarez | Estimation of leaf area and clumping indexes of crops with hemispherical photographs[END_REF] and showed a determination coefficient of R² = 0.86 and an absolute root mean square error (RMSE) of approximately 0.5 m 2 .m -2 [START_REF] Veloso | Modélisation spatialisée de la production, des flux et des bilans de carbone et d'eau des cultures de blé à l'aide de données de télédétection : application au sud-ouest de la France[END_REF]. [START_REF] Battude | Estimating maize biomass and yield over large areas using high spatial and temporal resolution Sentinel-2 like remote sensing data[END_REF] also compared BV-NNET derived GAI from several satellites, including SPOT and Formosat-2. They showed very good correlation (R=0.92 and RRMSE=23%) and performances that were similar to the ones found in the literature [START_REF] Berjón | Retrieval of biophysical vegetation parameters using simultaneous inversion of high resolution remote sensing imagery constrained by a vegetation index[END_REF][START_REF] Bsaibes | Albedo and LAI estimates from FORMOSAT-2 data for crop monitoring[END_REF][START_REF] Duan | Inversion of the PROSAIL model to estimate leaf area index of maize, potato, and sunflower fields from unmanned aerial vehicle hyperspectral data[END_REF]. The BV-NNET procedure did not include the aggregation of the leaves, which can lead to the underestimation of GAI during periods of strong vegetation development [START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF]. GAI estimates were finally averaged considering all the pixels of the studied plots after the application of an offset of 10 m to avoid edge effects and to consider only the GAI of the considered crop.

METHODOLOGY

The SAFY-CO2 model

The SAFY-CO2 model (Figure 4) was adapted from the SAFY model (Simple Algorithm for Yield Estimates; [START_REF] Duchemin | A simple algorithm for yield estimates: Evaluation for semi-arid irrigated winter wheat monitored with green leaf area index[END_REF] to simulate the components of the net CO2 fluxes and the cropland annual carbon budget.

SAFY is a daily time step crop model that simulates the temporal evolution of GAI, DAM and final grain yield (YLD) by considering two climatic input variables: incoming global radiation and mean temperature. This approach is based on [START_REF] Monteith | Climate and the Efficiency of Crop Production in Britain [and Discussion[END_REF] light-use efficiency theory, which links the production of the total DAM with the photosynthetically active portion of the solar radiation (PAR) absorbed by the plant. In SAFY, the ratio of photosynthesis to autotrophic respiration is assumed to be constant when estimating the DAM from the absorbed PAR (APAR). The SAFY model has been extensively used for the estimation of biomass and yield in contrasting climatic conditions and crop types [START_REF] Battude | Estimating maize biomass and yield over large areas using high spatial and temporal resolution Sentinel-2 like remote sensing data[END_REF][START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF][START_REF] Duchemin | Impact of Sowing Date on Yield and Water Use Efficiency of Wheat Analyzed through Spatial Modeling and FORMOSAT-2 Images[END_REF][START_REF] Duchemin | A simple algorithm for yield estimates: Evaluation for semi-arid irrigated winter wheat monitored with green leaf area index[END_REF][START_REF] Fieuzal | Combined use of optical and radar satellite data for the monitoring of irrigation and soil moisture of wheat crops[END_REF].

Conversely, in SAFY-CO2, the GPP is first estimated as a function of the APAR. Then, the components of the biomass (above and below ground) and the corresponding components of the net CO2 fluxes and annual carbon budget are calculated. In this section, the main formalisms and equations of the model are presented and the parameters are detailed in the tables. Each table summarized the notations, the values or the ranges and the methods for estimating the parameters for the winter wheat crop and the post-harvest vegetative events. photosynthesis (fELUE) [eq.1.3], and a multiplicative coefficient (sR10) [eq. 1.4.1 and 1.4.2], which takes into account the decline in canopy photosynthetic capacity during the senescence phase (see [START_REF] Béziat | Effet des conditions environnementales et des pratiques culturales sur les flux de carbone et d'eau dans les agrosystèmes[END_REF]. sR10 is set to 1 until senescence begins [eq. 1.4.1] and then defined as the ratio between the GAI of the previous day and the maximum seasonal GAI value multiplied by the corrective factor CS. The senescence phase first acts on the lower portion of the plant (closer to the soil) and then acts on the higher canopy elements. Thus, the actual phenological senescence may be greater than the phenological senescence detected by satellite observations, which thus requires a corrective factor. Therefore, the CS coefficient is included in the computation of sR10 to correct for the effects of senescence over simulated fluxes. The effects of diffuse global radiation over canopy photosynthesis are not always considered in crop models when estimating crop productivity. However, measurements, including the measurements at our flux sites, have indicated that the efficiency is very sensitive to the diffuse components of incoming global radiation [START_REF] Béziat | Effet des conditions environnementales et des pratiques culturales sur les flux de carbone et d'eau dans les agrosystèmes[END_REF][START_REF] Hollinger | Forest-atmosphere carbon dioxide exchange in eastern Siberia[END_REF][START_REF] Roderick | On the direct effect of clouds and atmospheric particles on the productivity and structure of vegetation[END_REF]). An effective light-use efficiency function is thus defined to account for the fraction of diffuse global radiation (see 3.2). Because diffuse incoming radiation Rdf is not often measured in the field, [START_REF] De | Coupling the Soil-Vegetation-Atmosphere-Transfer Scheme ORCHIDEE to the agronomy model STICS to study the influence of croplands on the European carbon and water budgets[END_REF] approach was used to estimate the Rdf/Rg ratio over the study area from the top of the canopy and the top of the atmosphere radiation data. The NPP is then derived from the difference between the GPP and the autotrophic respiration (Ra) [eq.2], which was separated into two components: maintenance respiration (Rm) and growth respiration (Rgr) (McCree, 1974) [eq.3]. Rm is calculated from the NPP of the previous day and a maintenance coefficient mR [eq.3.1], which corresponds to the fraction of maintenance respiration per NPP unit. Because Rm responds strongly to the temperature [START_REF] Amthor | The McCree-de Wit-Penning de Vries-Thornley Respiration Paradigms: 30 Years Later[END_REF], it was estimated using a "Q10 type" equation ( Van't Hoff, 1898) [eq.3.1.1]. In this equation, R10 is the reference respiration at 10°C. Rgr is calculated using the method described by [START_REF] Amthor | Respiration and Crop Productivity[END_REF] and improved by [START_REF] Choudhury | A sensitivity analysis of the radiation use efficiency for gross photosynthesis and net carbon accumulation by wheat[END_REF], as shown in eq.3.2. The constant Yg is the growth conversion efficiency.

= * * * * * 10 1 = 1 - * ! 1.1 = 1 -" #$ % #$ &'( ) * for +,-< < /01 1.2.1 = 1 -" #$ % #$ &% ) * for /01 < < + 2 1.2.2 = 0 for < +,-or + 2 < 1.2.
B = - 2 = + + D 3 + = B * E 7 * 10 3.1 E 7 = FG * H FG I % FG FG J 3.1.1 D = K1 -L M * -+ 3.2
Finally, the total NPP is divided into root (NPPr,[eq.4.1]) and aerial (NPPa, [eq.4.2]) components, estimated by considering a root-to-shoot ratio (RtS) in accordance with the method proposed by Baret et al., (1992) [eq.5]. In this equation, SMT is the sum of temperature, D0 is the emergence date and DS is the first day of the spiking stage. N G is the extrapolated value of the root fraction N at emergence, N O is the asymptotic value of N, and c

is the relative rate of decrease. The DAM is estimated by dividing the NPPa by the coefficient Cveg, which represents the plant carbon content [eq.6].

B D = B * TU 4.1 B = B * 1 -TU 4.2 TU = N = N O + N G -N O * V WX Y Z Y < WX Y Z WX Y < WX Y Z [ 5 \ ] = B @ ^> 6
Once the biomass computed, the grain yield and GAI can be estimated, as in the SAFY version. The GAI is the sum of the GAI of the previous day and the positive and negative change in GAI of the current day [eq.7]. Leaf production and leaf senescence are controlled by a growing degree-day approach. The positive increment [eq.7.1] is the product of a function of leaf partitioning [eq.7.1.1], the specific leaf area parameter and the daily

∆ a = a + ∆ a b -∆ a 7 ∆ a b = ∆\ ] * c * U 7.1 c = 1 -c * de 6 * WX 7.1.1 ∆ a = a * U] -U f U f A 7.2 Lg ch = \ ] + 2 * ia 8
The NEE is calculated as the difference between the NPP and the carbon losses due to heterotrophic respiration (Rh) [eq.8]. Rh is calculated using a Q10 first-order exponential equation [START_REF] Delogu | Modélis ation de la respiration du sol dans les agrosystèmes[END_REF]. Rhref is the reference respiration at 0°C, and exp b*10 is equal to Q10.

The ecosystem respiration is defined as the sum of Ra and Rh [eq.10].

B = B -j 8 j = k * A * < 9 k = j l 9
9.1

H FG = ℎ n + 10 ℎ n = k * A * < bFG k * A * < = A * FG 9.2 > / = + j 10

Model parameterization and calibration

The parameters of the SAFY-CO2 model are set using one of the following three options: i) literature sources, ii) multi-site in situ measurements and iii) optimization using the time series of satellite-derived GAI.

The parameters in the first two categories are set as equal for all of the investigated fields and years of study (see equations 1 to 10 and Tables 2 to 6), while the parameter of the third category are optimized. Parameter b of the exponential function in the relationships between ELUE and Rdf/Rg [eq. 1.3] is fixed to 1.34 based on field data (i.e., CO2 fluxes and meteorological data) acquired over several years for contrasting climatic, soils, and management conditions at 5 European instrumented sites, including ours (Lonzee, BE; Auradé, FR; Grignon, The parameters in the third category are set using an iterative minimizing method with a RMSE objective function between satellite-based GAI and SAFY-CO2 GAI estimates. This procedure aims at determining the values of parameter "a" of the fELUE function [eq.1.3] and 6 phenological parameters, i.e., the day of plant emergence (D0), the specific leaf area (SLA), the two parameters of the partition-to-leave function (Pla, Plb), the sum temperature for senescence (Sena) and the rate of senescence (Senb).

The minimization procedure, applied to each simulations (i.e. each cropping year and each field), is based on an adapted version of the Nelder-Mead simplex method [START_REF] Lagarias | Convergence Properties of the Nelder--Mead Simplex Method in Low Dimensions[END_REF], which considers a priori boundaries for each parameter to constrain the solutions within realistic parameter intervals. The minimization process runs the model, computes a cost function (in this case the RMSE derived from the comparison between the estimated and the remote sensed GAI) and iteratively updates the values of parameters to converge to the best parameter combination, coinciding with the lowest cost function value. To reduce the probability of local minima, a global approach is applied that runs the optimization process 30 times, with different a priori conditions for each parameter. A set of the parameters with the best solution is considered (i.e., lower RMSE for the GAI estimates). The number of optimization runs is set to 30, based on a sensibility analysis, so that the best combination of parameters is always retrieved, while avoiding unnecessary runs.

The optimization process requires boundaries for the parameters to calibrate (Pla, Plb, Sena, Senb, SLA, ELUEa, D0). These boundaries are first estimated based on a literature review. Then, a sensitivity analysis of the model is conducted to adjust these boundaries using a grid search. The ranges of the parameters are discretized, and all possible combinations are simulated (more than 3 million simulations). Then, the outputs are compared to the outputs obtained by optimizations performed using the adapted simplex method described above and the same parameter boundaries. This comparison allows us to do the following.

1-Verify that the adjusted boundaries reproduce all plant development conditions, with the constraint of a limited dispersion in the outputs.

2-Validate the efficiency of the adapted version of the simplex in retrieving the best set of parameters.

Compared to SAFY, the SAFY-CO2 version considers the biomass production of re-growth, weeds and cover crops (hereafter called post-harvest vegetative events) and their effect on the CO2 fluxes and annual C budgets. Indeed, these post-harvest events can have an important impact on the NEE [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF][START_REF] Ceschia | Management effects on net ecosystem carbon and GHG budgets at European crop sites[END_REF][START_REF] Poeplau | Carbon sequestration in agricultural soils via cultivation of cover crops -A metaanalysis[END_REF] and thus on the NECB. Therefore, the model is adapted to simulate these events without the distinction of their nature (i.e., spontaneous re-growth, weeds or cover crop).

In a first attempt to estimate the effects of post-harvest vegetative events on the net CO2 flux components and ultimately on the annual NECB, the same parameterization considered for winter wheat is also considered for all post-harvest vegetative events. For the studied years, when the satellite-derived GAI indicates the presence of vegetation on the field after harvest (i.e., AUR2006, LAM2009, LAM2011, and LAM2013), the model is first run to simulate the main crop. Then, a second optimization is performed on the vegetation following the main crop, optimizing ELUEa and the 6 phenological parameters. The boundaries of parameters Sena, D0 and ELUEa are changed compared to those fixed for winter wheat, while the ranges of the 4 other parameters remain the same.

From daily net ecosystem CO2 fluxes, NEE to the annual net ecosystem carbon budget, NECB

To compute the annual NECB [eq.11], carbon input (Cinp) and export (Cexp) terms are added to the annual cumulated NEE (i.e., the NEP) [eq.11.1]. The NEP is the carbon absorbed or released by the field (through photosynthesis and respiration processes) over a cultural year, a positive NEP indicates that cumulated soil and autotrophic respiration are higher than cumulated photosynthesis, meaning that the field loses carbon towards the atmosphere and vice versa. The value of NEP is computed from October 1 st to September 30 th because this period usually corresponds to the agricultural cropping year in Europe [START_REF] Ceschia | Management effects on net ecosystem carbon and GHG budgets at European crop sites[END_REF]. The term Cinp represents the amount of C that is brought to the field as organic fertilizer (only at FR-Lam in this study) and as seeds. Since Cinp could not be simulated, the Cinp values provided by the farmer as well as analysis of the organic fertilizer C content [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF], were used to calculate the NECB at LAM.

Cexp generally corresponds to the yield [eq.11.2] in the study area, as typically, only grain is exported from the field (e.g., FR-Aur). However, in some cases, straws are also exported from the field (e.g., FR-Lam). From the perspective of regional-scale applications, this term (strawexp) is estimated as a function of the total straw biomass (strawtot), which corresponds to the final aboveground biomass (DAMmax) minus the final grain yield B @o = B + @ >20 -@ ,-0 11

B = p B qG r W>01 F < s 1
11.1

@ >20 = Lg ch * @ ^> 11.2 @ >20 = @ ^> * Lg ch + \ ] + 2 -Lg ch * _ 11.3 _ = TNkt >20 TNkt 1/1 = TNkt >20 \ ] + 2 -Lg ch 11.3.1

Model implementation and validation strategy

The proposed approach requires several types of input data, which are taken into account at different stages of the modelling process or during the calculation of the carbon budgets. To run the model, meteorological data (Rg and Tair), soil parameters (heterotrophic respiration parameters) and plant parameters (some being fixed and some being calibrated (see Table 2 to 6 for more details) are needed. The model can simulate the output variables only with those information, however, in the absence of calibration process, the estimates are very likely to be erroneous. To correctly simulate the crop development and thus get meaningful outputs, the satellite derived GAI is used to calibrate Pla, Plb, Sena, Senb, ELUEa, SLA and Do. Once the calibration process achieved, the parameters are set for each field and climatic year allowing to reproduce the vegetation dynamic (GAI, DAM, YLD) as well as the CO2 fluxes (GPP, Reco, NEE). Because at this stage, only remote sensing data and crop maps are needed, the approach can be up-scaled easily. Nevertheless assessing carbon budgets demand additional data that are not easily available at plot scale over large areas yet. The two additional input data needed are: the amount of C imported (Cinp) to the plot (as seeds and eventually as organic amendments), and an information on whether straw are exported from the plot or not. These two inputs allow to estimate the NECB [eq. 11].

The validation strategy relies on different types of datasets. Concerning the GAI, the comparison of the model output with the satellite derived GAI is not really a validation since the later is used to calibrate the model but rather a verification of the model's ability to reproduce vegetation dynamic. To evaluate the performances of the model to simulate biomass, data from the instrumented ICOS flux sites and from the 2011 field campaign are considered. For yield validation, only the results from the 2011 field campaigns are presented (covering a range of values from 1.87 to 8.93 t.ha -1 ) as the annual yields provided by the farmer at the two instrumented sites may be averages of several fields from the farm and therefore present high uncertainties (see [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF].

These validations are presented in section 4.1.1.

To evaluate the performance of the model in terms of CO2 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 2.2.2) makes our validation approach more robust. Those results are presented in section 4.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 4.2.1 presents the ability of the model to reproduce the cumulated NEE, which leads to the NEP. Then sections 4.2.2 and 5 present and discuss, respectively, the capability of the model to assess the NECB.

RESULTS

Evaluation of the overall model performances

GAI, DAM and yield estimates

Figure 5 shows an overview of the statistical performances of the model in estimating GAI and DAM by ). 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 [START_REF] Duchemin | Impact of Sowing Date on Yield and Water Use Efficiency of Wheat Analyzed through Spatial Modeling and FORMOSAT-2 Images[END_REF][START_REF] Fieuzal | Combined use of optical and radar satellite data for the monitoring of irrigation and soil moisture of wheat crops[END_REF][START_REF] Hadria | Potentiality of optical and radar satellite data at high spatio-temporal resolutions for the monitoring of irrigated wheat crops in Morocco[END_REF].

On average, the model also reproduces the DAM with good precision. When considering all investigated fields across the studied years (Figure 6b), the RMSE, RRMSE and R 2 were 201 g.m -2 , 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 -2 , respectively, and RRMSEs of 24.7 and 39.1%, respectively, in [START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF], as well as the performances of other models such as the APSIM-wheat model (Asseng et al., 1998) and STICS [START_REF] Brisson | STICS: a generic model for simulating crops and their water and nitrogen balances. II. Model validation for wheat and maize[END_REF], which estimated winter wheat biomass with R² of 0.90 and 0.78, respectively, and RMSEs of 0.101 and 0.266 kg.m -2 , respectively. Nevertheless, the performances are more scattered when analysing annual statistics (Figure 5). 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 -2 and 2298 g.m -2 , 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 [START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF]. This issue will be further discussed in Section 5.2.

On average, yields (only data from the 2011 field campaign are analysed here, see 2.2.1) are estimated with good precision (RMSE=1.02 t.ha -1 , 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-CO2 are similar to those of SAFY [START_REF] Duchemin | A simple algorithm for yield estimates: Evaluation for semi-arid irrigated winter wheat monitored with green leaf area index[END_REF] and STICS [START_REF] Brisson | STICS: a generic model for simulating crops and their water and nitrogen balances. II. Model validation for wheat and maize[END_REF] for wheat, with R² values of 0.64 and 0.65, respectively, and RMSEs of 0.5 and 1.6 t.ha -1 , respectively. and correlation (R²) (see Figure 7). These statistical parameters are also calculated for the following periods of each cropping year (see Figure 8):

1-From the beginning of the cropping year (October 1 st ) 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 30 th )

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 -2 .d -1 .

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 1 st and March 31 th .

Consequently, the CO2 fluxes were gap filled, and thus, the simulated GPP is underestimated over this period, leading to low R² 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 -2 .d -1 . 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 in section 4.1.2).

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 2 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 3.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 postharvest vegetative event are considered.

Reco estimates

Overall, the ecosystem respiration follows the same dynamics as the GPP but in an attenuated way (i.e., the The mean annual performances of the model in estimating Reco are lower than those in estimating GPP but are consistent with other studies that aim to estimate ecosystem respiration in agricultural fields. [START_REF] Zhan | Modeled and Measured Ecosystem Respiration in Maize-Soybean Systems Over 10 Years[END_REF] found R² values between 0.84 and 0.87 and RMSEs between 1.52 and 1.65 gC.m -2 .d -1 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² between 0.71 and 0.79). A more detailed analysis of our results shows that the Reco is always well estimated during the crop growing periods, showing R² values between 0.88 and 0.95, while the statistical performances are lower in other periods than in the crop growing period (Figure 8). 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 5.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 4.1.2.4), 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 Reco during the crop season [START_REF] Béziat | Effet des conditions environnementales et des pratiques culturales sur les flux de carbone et d'eau dans les agrosystèmes[END_REF]. The modelled Reco estimates are thus satisfactory considering the simplicity of its representation as well as the limited number of inputs in this crop modelling approach.

NEE estimates

The daily dynamics of the NEE are well reproduced by the SAFY-CO2 model even when GPP and Reco 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 Reco). 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 -2 .d -1 and R 2 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 Reco due to the reasons mentioned above).

In order to compare the performances of SAFY-CO2 to simulate NEE with other agronomical or land surface models, our results were confronted to those presented by [START_REF] Wattenbach | The carbon balance of European croplands: A cross-site comparison of simulation models[END_REF]. In their study, they compared the performances of DNDC [START_REF] Li | Modeling impacts of farming management alternatives on CO 2 , CH 4 , and N 2 O emissions: A case study for water management of rice agriculture of China: WATER MANAGEMENT AND CHINA PADDY GREENHOUSE GAS FLUXES[END_REF][START_REF] Li | Modeling carbon biogeochemistry in agricultural soils[END_REF][START_REF] Li | A model of nitrous oxide evolution from soil driven by rainfall events: 1. Model structure and sensitivity[END_REF], ORCHIDEE-STICS (de [START_REF] De | Coupling the Soil-Vegetation-Atmosphere-Transfer Scheme ORCHIDEE to the agronomy model STICS to study the influence of croplands on the European carbon and water budgets[END_REF][START_REF] Gervois | Carbon and water balance of European croplands throughout the 20th century: CARBON BALANCE OF EUROPEAN CROPLANDS[END_REF], SPA [START_REF] Williams | Modelling the soil-plant-atmosphere continuum in a Quercus-Acer stand at Harvard Forest: the regulation of stomatal conductance by light, nitrogen and soil/plant hydraulic properties[END_REF] and CERES-EGC [START_REF] Gabrielle | Process-based modeling of nitrous oxide emissions from wheat-cropped soils at the subregional scale: REGIONAL N 2 O EMISSIONS FROM ARABLE SOILS[END_REF][START_REF] Lehuger | Bayesian calibration of the nitrous oxide emission module of an agro-ecosystem model[END_REF] in reproducing the GPP, Reco and NEE for several site-years, including the AUR2006 crop season.

SPA and SAFY-CO2 outperformed the other approaches by reproducing the dynamics of the cumulated NEE, with R² values of 0.993 and 0.995, respectively. In terms of errors, SAFY-CO2 and CERES-EGC showed the lowest RMSEs, 33.6 gC.m -2 and 44.16 gC.m -2 , respectively.

Temporal evolution

We decided to show the performances of SAFY-CO2 in simulating the time courses of the GPP, Reco 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 Reco estimates that can affect the NEE and the net annual CO2 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.

• 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 9a), the model correctly reproduces the GPP and the Reco in terms of errors (RMSEs of 1.38 and 0.87 gC.m -2 .d -1 , respectively) and dynamics (R² of 0.92 and 0.78, respectively). Consequently, the NEE for this year is accurately estimated (RMSE and R² of 1.09 and 0.88, respectively). Nevertheless, the modelled GPP is slightly overestimated after maximum development and at the beginning of senescence (Figure 9a). This overestimation could be related to the water stress conditions observed after the maximum GAI was reached [START_REF] Béziat | Evaluation of a simple approach for crop evapotranspiration partitioning and analysis of the water budget distribution for several crop species[END_REF] showing RMSE values of 1.34, 0.82 and 1.14 gC.m -2 .d -1 , 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 postharvest GPP and Reco 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 -2 .d -1 ).

• 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 21 st of August and incorporated in the soil on the 4 th 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 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 -2 ) compared to that of other years (mean and standard deviation over the 8 studied years: 1566 +/-453 g.m -2 ). 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 CO2 fluxes (particularly on the GPP and the Reco) observed just before and after harvest. The difference between the observed and simulated NEE dynamics is small, either because the "observed" GPP and Reco partly balanced each other or because the increases in the "observed" Reco and GPP are caused by errors in the NEE partitioning process.

From the cumulated NEE to yearly carbon budget

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 1 st until September 30 th 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-CO2 allowed for the identification of the CO2 net assimilation and release phases (Figure 10). A negative slope in the cumulated NEE curve (i.e., corresponding to net assimilation, with GPP>Reco) is observed during the growing season and during crop re-growth or weed or cover crop development. A positive slope (i.e., corresponding to CO2 release, with GPP<Reco) 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 Reco 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 Reco 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 -2 to 33 and 16 gC.m -2 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.

Carbon budget over 8 agricultural seasons of winter wheat

The NECB model estimates and its components (NEP, Cexp) are compared with those of the eight site-year measurements (Table 7). In this table, the effects of post-harvest vegetative events on NEE and NECB are considered. For FR-Lam, the Cinp 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, Cinp 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-CO2, which vary from -191 gC.m -2 .yr -1 (LAM2009) to -486 gC.m -2 .yr -1 (LAM2011), are generally close to the measured values, which vary from -208 gC.m -2 .yr -1 (LAM2009) to -410 gC.m -2 .yr -1 (LAM2011). Table 7 shows that for all site-years, winter wheat is a CO2 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 Rh 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 Rh (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%).

The analysis of the amount of carbon exported from the ecosystem at harvest (Cexp) shows that this amount varies considerably from one site-year to another. The simulated Cexp varies from 253 gC.m -2 .yr -1 (AUR2010) to 436 gC.m -2 .yr -1 (LAM2013), while the observations range from 204 gC.m -2 .yr -1 (AUR2010) to 488 gC.m -2 .yr -1 (LAM2013). Cexp 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 Cexp. The differences between the observations and model estimates of Cexp vary between 5% and 46%. However, the comparison of the modelled Cexp and the observed Cexp 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 2.2.1.). For this reason, [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF] concluded that the uncertainty of in situ Cexp is often higher than the uncertainty of the NEP at our sites. For AUR2012, for instance, the Cexp estimated from the farmer's data is 223 gC.m -2 .yr -1 , while our destructive measurements encompassing over 30 subplots of 3.75*10 -2 m 2 in the field suggest a Cexp of 406+/-53 gC.m -2 .yr -1 . When analysing the performance of the model against our own destructive field samples, the model generally performs much better (see Figure 6b and6c). This indicates that our modelling approach may perform better for estimating Cexp (and therefore NECB) than what is presented in Table 7.

The NECB estimated from in situ data or from SAFY-CO2 characterize all of the site-years as carbon sinks. The modelled NECB estimates vary from -304 gC.m -2 .yr -1 (LAM2007) to -6 gC.m -2 .yr -1 (AUR2010), while the in situ NECB estimates vary from -439 gC.m -2 .yr -1 (LAM2007) to -29 gC.m -2 .yr -1 (AUR2014). The FR-Lam site-years present the greatest variations between years, and the carbon inputs (NEP + Cinp) are stronger than those of FR-Aur, partly because of organic fertilization.

As shown in Table 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 -2 .yr -1 (AUR2014) to 135 gC.m -2 .yr -1 (LAM2007) but the sign of NECB is always similar between both approaches. Regarding all simulated years, the model shows RMSE of 77 gC.m -2 .yr -1 . 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, Cinp, Cexp) range from 40 to 105 gC.m -2 .d -1 . 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).

DISCUSSION

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 CO2 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-CO2 achieves equivalent or better performances regarding estimates of key components of the C budget: CO2 fluxes, biomass and yield (see sections 4. 1.1. and 4.1.2.).

Next, we demonstrated the ability of the model to reproduce winter wheat dynamics, production and CO2 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 CO2 fluxes. As showed by [START_REF] Ceschia | Management effects on net ecosystem carbon and GHG budgets at European crop sites[END_REF] 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-CO2 model for those post-harvest vegetative events, the overall performances of the model for simulating CO2 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).

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 CO2 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 [START_REF] Ameline | Estimation of Corn Yield by Assimilating SAR and Optical Time Series Into a Simplified Agro-Meteorological Model: From Diagnostic to Forecast[END_REF][START_REF] Baup | Temporal Evolution of Corn Mass Production Based on Agro-Meteorological Modelling Controlled by Satellite Optical and SAR Images[END_REF][START_REF] Betbeder | Assimilation of LAI and Dry Biomass Data From Optical and SAR Images Into an Agro-Meteorological Model to Estimate Soybean Yield[END_REF][START_REF] Fieuzal | Estimation of Sunflower Yield Using a Simplified Agrometeorological Model Controlled by Optical and SAR Satellite Data[END_REF][START_REF] Revill | Carbon cycling of European croplands: A framework for the assimilation of optical and microwave Earth observation data[END_REF]. 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 [START_REF] Brown | High-resolution measurements of scattering in wheat canopies-implications for crop parameter retrieval[END_REF].

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 scenarii on the C budgets (e.g. accounting or not for the effects of the postharvest 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 [START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF][START_REF] Duchemin | Impact of Sowing Date on Yield and Water Use Efficiency of Wheat Analyzed through Spatial Modeling and FORMOSAT-2 Images[END_REF][START_REF] Duchemin | A simple algorithm for yield estimates: Evaluation for semi-arid irrigated winter wheat monitored with green leaf area index[END_REF] respectively) and for different crop species (corn, soybean, sunflower and wheat). Also, the parametrisation of the Q10 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, [START_REF] Delogu | Improved methodology to quantify the temperature sensitivity of the soil heterotrophic respiration in croplands[END_REF] obtained good Rh estimates (RMSE comprised between 0.15 and 0.73 gC.m -2 .d -1 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 Rh which should be considered as a first step for estimating Rh. In the future, this method could be improved, with little changes in the formalisms, by considering a R10 parameter that depends on top soil slow carbon content as it is proposed in [START_REF] Delogu | Improved methodology to quantify the temperature sensitivity of the soil heterotrophic respiration in croplands[END_REF]. Another step of improvement could be to account for the priming effect following the incorporation of fresh organic matter into the soil [START_REF] Kuzyakov | Review of mechanisms and quanti®cation of priming e ects[END_REF]. 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 [START_REF] Eugster | Management effects on European cropland respiration[END_REF] and no clear effect on SOC mineralisation (e.g. in [START_REF] Dimassi | Effect of nutrients availability and long-term tillage on priming effect and soil C mineralization[END_REF][START_REF] Powlson | Does conservation agriculture deliver climate change mitigation through soil carbon sequestration in tropical agro-ecosystems?[END_REF][START_REF] Virto | Carbon input differences as the main factor explaining the variability in soil organic C storage in no-tilled compared to inversion tilled agrosystems[END_REF]. Also, a potential drawback of our approach for estimating Rh is that the effect of the soil water content is not accounted for. As for the SAFY model [START_REF] Battude | Modeling water needs and total irrigation depths of maize crop in the south west of France using high spatial and temporal resolution satellite imagery[END_REF][START_REF] Duchemin | Impact of Sowing Date on Yield and Water Use Efficiency of Wheat Analyzed through Spatial Modeling and FORMOSAT-2 Images[END_REF], we have already tested the coupling of a soil water module (FAO56; [START_REF] Allen | Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56[END_REF] to the SAFY-CO2 model (see [START_REF] Veloso | Modélisation spatialisée de la production, des flux et des bilans de carbone et d'eau des cultures de blé à l'aide de données de télédétection : application au sud-ouest de la France[END_REF]. 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-CO2 -soil water module coupled approach at plot scale (e.g. the performances of GlobalSoilMap in France are R² = 0.27 and RMSE = 128 g.kg -1 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 Q10 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-CO2 model with a soil organic matter model (e.g. RothC, [START_REF] Coleman | RothC-26.3 -A Model for the turnover of carbon in soil[END_REF]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 [START_REF] Castaldi | Evaluating the capability of the Sentinel 2 data for soil organic carbon prediction in croplands[END_REF][START_REF] Vaudour | Sentinel-2 image capacities to predict common topsoil properties of temperate and Mediterranean agroecosystems[END_REF] 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 [START_REF] Claverie | Maize and sunflower biomass estimation in southwest France using high spatial and temporal resolution remote sensing data[END_REF]. 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 [START_REF] Betbeder | Assimilation of LAI and Dry Biomass Data From Optical and SAR Images Into an Agro-Meteorological Model to Estimate Soybean Yield[END_REF] 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 CO2 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 -2 .yr -1 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 insitu NECB calculated with more accurate yield data. Still, even if SAFY-CO2 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 [START_REF] Coleman | RothC-26.3 -A Model for the turnover of carbon in soil[END_REF] and CENTURY [START_REF] Parton | Analysis of Factors Controlling Soil Organic Matter Levels in Great Plains Grasslands[END_REF]. 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% [START_REF] Falloon | Simulating SOC changes in long-term experiments with RothC and CENTURY: model evaluation for a regional scale application[END_REF][START_REF] Guo | Application of the RothC model to the results of long-term experiments on typical upland soils in northern China[END_REF][START_REF] Smith | A comparison of the performance of nine soil organic matter models using datasets from seven long-term experiments[END_REF] for Roth-C model and between 1.8-16.4% for CENTURY [START_REF] Cong | Evaluation of the CENTURY Model Using Long-Term Fertilization Trials under Corn-Wheat Cropping Systems in the Typical Croplands of China[END_REF][START_REF] Falloon | Accounting for changes in soil carbon under the Kyoto Protocol: need for improved long-term data sets to reduce uncertainty in model projections[END_REF]. These results, which are more accurate than those achieved with SAFY-CO2 model should be tempered by the fact that they represent two different approaches, requiring different input and designed for different purposes and timescales. Also, as in other studies [START_REF] Ceschia | Management effects on net ecosystem carbon and GHG budgets at European crop sites[END_REF][START_REF] Schmidt | The carbon budget of a winter wheat field: An eddy covariance analysis of seasonal and inter-annual variability[END_REF] our results showed that for all cropping years the plots behave as a net CO2 sinks and our results concerning the potential C storage of winter wheat crop are consistent with other studies (e.g. [START_REF] Aubinet | Carbon sequestration by a crop over a 4-year sugar beet/winter wheat/seed potato/winter wheat rotation cycle[END_REF][START_REF] Ceschia | Management effects on net ecosystem carbon and GHG budgets at European crop sites[END_REF]. 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 [START_REF] Kaye | Using cover crops to mitigate and adapt to climate change[END_REF][START_REF] Pellerin | Stocker du carbone dans les sols français, quel potentiel au regard de l'objectif 4 pour 1000 et à quel coût[END_REF][START_REF] Poeplau | Carbon sequestration in agricultural soils via cultivation of cover crops -A metaanalysis[END_REF][START_REF] Tribouillois | Cover crops mitigate direct greenhouse gases balance but reduce drainage under climate change scenarios in temperate climate with dry summers[END_REF].

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.

CONCLUSION

In this work, we demonstrate the potential of high-resolution remote sensing data assimilation in a semi-physical crop model (SAFY-CO2) to successfully provide estimates of some of the main components of cropland annual carbon budgets (i.e., net CO2 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 CO2 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 Q10 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 [START_REF] Veloso | Understanding the temporal behavior of crops using Sentinel-1 and Sentinel-2-like data for agricultural applications[END_REF]. In addition, our results show that the performance of the model in estimating net CO2 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.

  the Garonne River. The Garonne River flows from the South to the North on the east side of the study area. It has spread sediments over a 15 km wide area along its western side resulting in vast terraces of heterogeneous soils called "boulbènes" and "terrefort" characterized by low-permeability and composed of a silt layer of variable thickness over stony clay soils. The geology is old quaternary and the main lithology is old alluviums. The area west of the terraces is characterized by a hilly landscape, consisting of hills and slopes resulting from the erosion of the oldest terraces. Further west the landscape is hilly over hundreds of kilometres and the soils become more calcareous with deposits formed of marl and clayey molasses with limestone. The heterogeneous character of the soils of the study area is illustrated (Figure 2) by the texture measurements (fractions of clay, silt and sand) collected on the flux sites and on a network of fields within the footprint of the satellite images during the year 2018 (see also https://soilgrids.org/ for predicted soil classification of the study area). The texture

Figure 1

 1 Figure 1: The upper right corner shows the location of the study site in southwestern France, as well as the footprint of Formosat-2 (green square) and Spot (orange square) images (in 2014 and 2012, respectively). The SPOT-4 false color image used as the background shows the flux sites of FR-Aur and FR-Lam (zoomed areas), the network of fields sampled for biomass and yield during the 2011 field campaign (yellow points), and the SAFRAN meteorological grid (black crosses).

Figure

  Figure 2: Surface texture measurements (FR-Lam in red, FR-Aur in blue and field campaign in black) displayed on USDA classification, with the following classes: clay (Cl), silty clay (SiCl), sandy clay (SaCl), clay loam (ClLo), silty clay loam (SiClLo), sandy clay loam (SaClLo), loam (Lo), silty loam (SiLo), sandy loam (SaLo), silt (Si), loamy sand (LoSa), sand (Sa).

Figure 3 Figure 3 :

 33 Figure 3 presents a chronogram of the satellite images used in this study between 2006 and 2014. The images from those different satellites were combined to better monitor crop development and to reduce the gaps between successive observations. Nevertheless, the presence of clouds and/or shadows reduced the number of useful images. For instance, only one cloud-free SPOT image (April 26 th ) was available from mid-February until mid-June 2008; consequently, this site-year was not processed in the present study.The high spatial resolution images provided by Formosat-2 (F2, 155 images) and Spot-2/4/5 (80 images) were in

FirstFigure 4 :

 4 Figure 4: Schematic representation of the assimilation procedure of GAI derived from high resolution satellite optical images for the calibration of the agro-meteorological model SAFY-CO2 by minimizing difference between satellite derived (SAT) and simulated (SIM) GAI. Also the figure shows the procedure for estimating the crop biomass, the components of the net CO2 fluxes (GPP, RECO, NEE) and the annual carbon budgets (NECB) over a cropping season (see Tables 2 to 6 and equations 1 to 11.3.1 for more details concerning the processes simulated and the parameters).

[

  eq.11.3].The sc parameter [eq.11.3.1] is estimated from in situ data during the 2011 field campaign and set to 0.3.

Figure 5 :

 5 Figure 5: 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 6: 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.

4. 1 . 2 .

 12 Components of the net CO2 fluxes: daily GPP, Reco and NEE In this section, the components of the net CO2 fluxes simulated by SAFY-CO2 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)

Figure 7 :

 7 Figure 7: 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 8 :

 8 Figure 8: Boxplot of the pluriannual determination coefficients (R²) and the pluriannual root mean square errors (RMSEs) of the simulated GPP, Reco and NEE in the before-crop, crop-growth, cropsenescence, 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.

Figure 9a :

 9a Figure 9a: 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 9b :

 9b Figure 9b: Temporal evolution of the measured (in blue) and estimated (in red or red dashed lines) cumulated GPP, Reco 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.

Figure 10 :

 10 Figure 10: Temporal evolution of the measured (in blue) and estimated (in red and red dashed lines) daily values of the GPP, Reco 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).

Table 7 :

 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).

Figure

  Figure S.3.2: Temporal evolution of the measured (in blue) and estimated (in red) GPP, Reco 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

  Figure S.3.3: Temporal evolution of the measured (in blue) and estimated (in red or red dashed lines) GPP, Reco 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

  Figure S.3.4: Temporal evolution of the measured (in blue) and estimated (in red lines) GPP, Reco 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

  Figure S.3.5: 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

  Figure S.3.6: Temporal evolution of the measured (in blue) and estimated (in red lines) GPP, Reco 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.

  Aur and FR-Lam sites have been intensively monitored since 2005. Micrometeorological, meteorological, soil and vegetation measurements are performed since then (see Béziat et al., 2009 for more details). Both sites have similar climatic conditions but different soil properties (see 2.1.1), topography and agricultural management practices. Winter wheat was cultivated throughout 8 cropping years

Table 1 : Overview of the in situ data collected from 2005 until 2014, agricultural practices (dates of sowing

 1

	Year	Site/ESU	Date of sowing Date of harvest	Vegetation after harvest (nature)	GAI	Biomass	Grain yield	Flux and meteorological measurements
	2006	AUR	27/10/2005	29/06/2006 Yes (weeds/re-growth)	9	9	F	Yes
	2007	LAM	18/10/2006	15/07/2007	No	11	11	F	Yes
	2009	LAM	19/11/2008	13/07/2009	Yes (weeds)	8	8	F	Yes
	2010	AUR	19/11/2009	12/07/2010	No	5	5	F	Yes
	2011	LAM Field Campaign	03/11/2010 -	02/07/2011 -	Yes (re-growth) -	5 -	5 16	F D	Yes No
	2012	AUR	21/10/2011	14/07/2012	No	5	5	F	Yes
	2013	LAM	29/10/2012	22/07/2013	Yes (cover-crop)	5	5	F	Yes
	2014	AUR	26/10/2013	10/07/2014	No	5	5	F	Yes

harvest, etc.), and presence and type of vegetation during the fallow period. F: provided by farmer, D: destructive measurements.

Table 2 . List of SAFY-CO2 model parameters for calculating the GPP

 2

	3

Table 3 . List of SAFY-CO2 model parameters for calculating autotrophic respiration

 3

	Description	Notation	Unit	Value/Range	Method	Source
	Maintenance respiration parameter: Q10 Maintenance respiration parameter: R10 Growth respiration	H FG FG	-gC/gDM	2 0.0025	Literature literature	Amthor 2000 Béziat 2009
	conversion efficiency parameter	L	-	0.74	Literature	Amthor 1989

Table 4 . List of SAFY-CO2 model parameters for calculating the aboveground and the belowground mass

 4 The negative increment, which is only evaluated from the beginning of senescence, depends on the two senescence parameters Sena and Senb [eq.7.2]. The grain yield estimation [eq.8] depends on the total biomass production at the end of the vegetative period (DAMmax) and a constant harvest index HI (see 3.2).

	Description	Notation	Unit	Value/Range	Method	Source
	Root fraction parameters Carbon content coefficient	G -O -_ @ ^>	-gC/gveg	0.63 / 0.11 / 1.48 0.46 [20 th Oct-15 th Jan]	Literature Literature	Baret et al., 1992 Béziat 2009
				[end of main
	Day of plant emergence	\ G	day	(post-harvest vegetative crop-31 th Dec]	Calibration	-
				events)

Table 5 . List of SAFY-CO2 model parameters for calculating the leaf biomass production, the yield and the senescence

 5

	Description	Notation	Unit	Value/Range	Method	Source
	Harvest index Specific leaf area Partition-to-leaf function parameter a Partition-to-leaf function parameter b	ia U c c A	-m 2 .g -1 --	0.45 [0.005 -0.04] [0.01-0.5] [0.0001-0.02] [1045-2000]	Literature Calibration Calibration Calibration	(Dai et al., 2016) ---
	Sum of temperature for senescence	U f	°C	[100-900] (post-vegetative events) harvest	Calibration	-
	Rate of senescence	U f A	°C.day -1	[10 3 -2.10 4 ]	Calibration	-

Table 6 . List of SAFY-CO2 model parameters for calculating the heterotrophic respiration Description Notation Unit Value/Range Method Source

 6 Lamasquère, FR; Oensingen, CH; see Supplementary Material). This parameter is set to a generic value to facilitate the large-scale application of the approach.

	Heterotrophic respiration parameter: Rhref Heterotrophic respiration parameter: Q10 Conversion factor of Ta into Ts	j l 9 H FG T	gC.m -2 .d -1 --	0.34 2.3 1.07	Literature Literature ICOS sites data	(Suleau et al., 2011) (Suleau et al., 2011) -

  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 Reco 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

9a

); as a consequence, GPP and Reco are overestimated. NEE estimates during this period are improved (RMSE from 1.22 to 0.93 gC.m -2 .d -1 ) 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, Reco and NEE present R² values of 0.78 and 0.81, respectively, over the cultivation year. The errors of the GPP, Reco and NEE are also low,

Acknowledgements

This work was made possible through the support of the Agence De l'Environnement et de la Maîtrise de l'Energie (ADEME), which financed the project "Couverts Intermédiaires pour l'atténuation du Changement Climatique" (CICC) and half of Gaétan's and Morgan's PhDs, the Centre National d'Etudes Spatiales (CNES), which financed half of Gaétan's PhD, and the Agence de l'Eau Adour Garonne (AEAG), which financed the Bag'ages project and the support of the European Union that financed the Sensagri project (Horizon 2020

Research and Innovation Programme, Grant Agreement n°730074).

Data acquisition at FR-Lam and FR-Aur were mainly funded by the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS-INSU) through the ICOS and OSR SW observatories. Facilities and staff were also funded and supported by the University Toulouse III -Paul Sabatier, the CNES and IRD (Institut de Recherche pour le Développement). We are grateful to Franck Granouillac, Nicole Claverie, Bernard Marciel and Pascal Keravec for their technical support, advice, and valuable assistance in the field and with site management. We extend special thanks to Mr. Andréoni (farmer) and to the Ecole d'Ingénieur de Purpan for accommodating our measurement devices in their fields at FR-Aur and FR-Lam respectively.

We are also very grateful to Pauline Buysse and Benjamin Loubet (respectively Grignon's site manager and PI), Bernard Heinesch (Lonzee's PI), Nina Buchmann (Oensingen's PIs) and Christian Bernhofer (Klingenberg's PI) for granting and facilitating our access to their flux site data (destructive GAI, CO2 flux, radiation and temperature measurement and management data), which allowed us to account for diffuse radiation effects in GPP estimates.

Supplementary material S.1 Effect of diffuse radiation on the ELUE

The fraction of diffuse radiation affects photosynthesis [START_REF] Béziat | Carbon balance of a three crop succession over two cropland sites in South West France[END_REF] 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 [START_REF] Roderick | On the direct effect of clouds and atmospheric particles on the productivity and structure of vegetation[END_REF]. Thus, the photosynthetic efficiency is expected that 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 nonlinearly with the ratio of diffuse over total global radiation. 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.