# Calibration of the STICS crop model combining remote sensing data and farmer information

Elena Pareja-Serrano, André Chanzy, Jose González-Piqueras

► To cite this version:

Elena Pareja-Serrano, André Chanzy, Jose González-Piqueras. Calibration of the STICS crop model combining remote sensing data and farmer information. European Journal of Agronomy, 2025, 170, pp.127756. 10.1016/j.eja.2025.127756. hal-05148593

HAL Id: hal-05148593

https://hal.inrae.fr/hal-05148593v1

Submitted on 7 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Calibration of the STICS crop model combining remote sensing data and farmer information

Elena Pareja-Serrano<sup>a, c</sup>, André Chanzy<sup>b, *</sup> <sup>D</sup>, Jose González-Piqueras<sup>c</sup>

## ARTICLE INFO

**Keywords:**

STICS crop model

Plant calibration

Spring wheat. Vegetation phenology

Yield

Optical remote sensing

## ABSTRACT

In climate change context, crop production assessment is a worldwide issue for agriculture that can be addressed using crop models like STICS, which simulate plant-environment interactions. Calibrating cultivar parameters is a key preliminary step that requires a substantial amount of field data, which might be a constraint for operational application. The objective of this work is to propose a methodology for calibrating STICS plant parameters using remote sensing and information from farm (soil, management, yield). The methodology was applied to spring wheat varieties (*Triticum aestivum* L.) grown in the *Central Plateau* of Spain. To avoid confusions between genetic traits and impact of edaphic stresses, the calibration was done on well irrigated and fertilized fields, while the validation encompassed a wider range of conditions. Green area index (GAI<sub>RS</sub>) time series computed on multispectral Sentinel-2 and Landsat 8 satellite images were used to characterize leaf cover development and phenology. The results led to *rmse* ranging between 6 and 14 days for the different phenological stages simulated by STICS, 1.3 m<sup>2</sup> m<sup>-2</sup> for GAI and 1.7 t ha<sup>-1</sup> for yield. These results were found better than using an existing parameterization established using large experimental data sets on a similar spring wheat cultivar. The evaluation of the calibrated crop model on water-stressed conditions was satisfactory showing the ability of the STICS crop model to represent environmental stresses even though the model was calibrated under optimal conditions. The adding value of using field phenological observations was marginal while calibration on biomass observations was required to accurately simulate biomass. However, it has been shown that a biased biomass simulation had little impact on yield and GAI. This work shows that by using readily available data such as data recorded by farmers or derived from remote sensing, it is possible to calibrate the plant parameters of a crop model leading to yield and GAI simulations better than those computed with existing set of plant parameter calibrated on comparable varieties.

## 1. Introduction

Crop production assessment is a global concern in the context of climate change. The Intergovernmental Panel on climate change Fifth Synthesis Assessment Report (IPCC, 2014) notes climate change impacts in agriculture include shifts in production areas around the world, which, combined with increasing food demand, pose significant risks to food security and intensify competition for water. In the past decades, several process-based models have been developed in agriculture to represent plant-environment-management interactions such as CERES (Ritchie and Otter, 1985), EPIC (J. R. Williams et al., 1984), DSSAT

(Jones et al., 2003), APSIM (Keating et al., 2003), CropSyst (Stöckle et al., 2003), STICS (Brisson and Launay, 2008), MONICA (Nendel et al., 2011), WOFOST (Boogaard et al., 2014), AquaCrop (Steduto et al., 2009), and DAISY (Abrahamsen and Hansen, 2000). Literature shows that crop models have the potential to understand the vulnerability of regions to climate change by predicting crop growth and yield under different growing conditions, including water and nutritional stress (Katerji et al., 2013; Jin et al., 2015; Kassie, 2016; Rahimi-Moghaddam et al., 2021) as well as variations in temperature, CO<sub>2</sub>, and irradiance (Asseng et al., 2013; Rosenzweig et al., 2013; Ewert et al., 2015; Tao et al., 2020). Other authors applied crop models for mitigation and

<footnote><sup>a</sup> AgriSat Iberia S.L., Paseo de la Innovación, 1, Albacete 02006, Spain</footnote>

<footnote><sup>b</sup> 1114 UMR INRAE-Avignon University EMMAH, Domaine St. Paul, Avignon 84914, France</footnote>

<footnote><sup>c</sup> Remote Sensing and GIS Group, Regional Research Institute, Campus of Albacete, University of Castilla-LaMancha, Albacete 02071, Spain</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: elena.pareja@agrisat.es (E. Pareja-Serrano), andre.chanzy@inrae.fr, andre.chanzy@inrae.fe (A. Chanzy), Jose.Gonzalez@uclm.es (J. González-Piqueras).</footnote>
adaptation purposes, exploring adapted management practices (Hadria
et al., 2007; Nendel et al., 2014; Ngwira et al., 2014; Bécel et al., 2015;
Kadiyala et al., 2015; Bai and Gao, 2021; Justino et al., 2022) and
supporting crop breeding (Rötter et al., 2015; He et al., 2017; Yang et al.,
2022).

Crop models represent plant-environment-management interactions
through equations driven by plant parameters, which may vary for each
variety of a given species (Manivasagam and Rozenstein, 2020).
Therefore, in addition to model structure errors, plant parameters as
well as the inputs describing crop growing conditions, i.e. weather, soil,
management, and initialisation, are the main sources of uncertainty
when running crop models (Jin et al., 2018; Hao et al., 2021). They are
usually adapted based on literature, in situ measurements, and cali-
bration procedures (Attia et al., 2016; Razzaghi et al., 2017; Seyoum
et al., 2018; Falconnier et al., 2019). In calibration procedures, simu-
lations of crop development and growth are matched to the corre-
sponding observations by fitting plant parameters. Jégo et al. (2011)
show some limitations in the calibration in dry conditions, highlighting
the importance of mastering the environmental conditions both in the
calibration phase and in the model implementation phase. Therefore, it
is necessary to mobilise large data sets describing crop development and
growth as well as growing conditions, which are not always easy to
collect, being one of the limiting factors for the calibration and, conse-
quently, the implementation of crop models. Furthermore, the fact that
the selection of parameters and the order of the steps in the calibration
procedures can have an impact on the results cannot be neglected.
Guillaume et al. (2011) showed that focusing on a variable such as yield,
gives better results for yield than attempting to represent both yield and
leaf area index development. Jégo et al. (2011), (2012) suggested that
calibration can be simplified to one variety per region, thus limiting to a
reduced number of parameters while maintaining good predictive
quality over the area considered. Therefore, when facing a calibration
process, we are confronted with the dual challenge of data availability
and standardisation of calibration protocols. The use of remote sensing
can be a way to overcome these two limitations, as will be shown in this
work.

Remote sensing gives access to key biophysical variables for crop
monitoring, such as the area of green leaves in a canopy (known as green
area index or GAI) (Weiss et al., 2020), the fraction of absorbed radiation
(fAPAR), biomass or phenological stages (Zhang et al., 2003; Sakamoto
et al., 2005; Lobell et al., 2013; González-Gómez et al., 2018; Zhou,
2019; Meroni et al., 2021; Liu et al., 2022). Recent advances in remote
sensing technology have significantly improved the potential for accu-
rate and reliable estimates of canopy biophysical variables (Jin et al.,
2018), such as remote sensors on-board the satellites Landsat 8 and
Sentinel-2, which offer a wide range of open-access spectral observa-
tions of the Earth's surface. The launch of Sentinel-2 satellites in 2015
and 2017 has boosted remote sensing capabilities by providing un-
precedented high spatial and temporal resolution (decametric and 5-day
revisit time) (Weiss et al., 2020), resulting in a wealth of open data sets
suitable for working at the field scale with a revisit time enabling
effective monitoring of crop dynamics throughout the growing season
simultaneously over wide areas, and offering a powerful alternative to
sampling data.

The coupling between remote sensing data and crop models has been
explored for a long time. There exists an abundant literature on the use
of remotely sensed data to improve crop models simulations by miti-
gating the impact of uncertainties associated to model plant parameters.
Reinitialisation of some of the parameters or the initial conditions of the
models is the most common approach. Maas (1988) proposed to improve
the yield predictions of a simple crop growth model, reinitializing the
model's state variable GAI through tuning the initial values to match
simulation with remote sensing observations. Guérif and Duke, (2000)
and Jégo et al. (2012) worked on the initial stages of the vegetation.
They use the time series of GAI to estimate the emergence date, which is
a critical stage for the simulation of the subsequent phenological stages.

Prévot et al. (2003), meanwhile, focused on the sowing date. An alter-
native way of using remotely sensed data is by forcing crop model
simulations (Huang et al., 2001). For instance, Casa et al. (2012) used
STICS model as a transformation function to infer performance from GAI
and showed the interest in having frequent and accurate measurements
overcoming some uncertainties in the management data. Data assimi-
lation schemes are an improvement of the forcing method. For instance,
the ensemble Kalman Filter (EnKF) is an optimal way to correct a model
simulation by taking into account model and observation errors as well
as the correlation between variables (de Wit and van Diepen, 2007;
Nearing et al., 2012; Zhao et al., 2013; Linker and Ioslovich, 2017;
Fattori Junior et al., 2022; Zhuo et al., 2022). Most of these studies used
a plant-calibrated crop model and have shown a gain in simulation ac-
curacy. Fattori Junior et al. (2022) underline the fact that the gain is
maximum with non-calibrated plant genotypes. However, some limita-
tions were highlighted when the link between the target variable and
observable quantities is not well simulated by the model due to either
the model parameters or process parameterization (Nearing et al., 2012;
Fattori Junior et al., 2022), highlighting the need of well-calibrated crop
models.

While the literature is abundant on the use of remote sensing data to
correct a simulation by attenuating the impact of uncertainties in the
model's inputs, there are relatively few studies on the use of remote
sensing for calibration purposes. For instance Prévot et al. (2003) used
airborne data to demonstrate the benefits of having large number of
measurements to calibrate two key STICS parameters regulating
phenology and leaf development, obtaining a good simulation of GAI.
Rodriguez et al. (2004) focused on the calibration of the parameter
describing GAI development in STICS and showed a positive impact of
this calibration on processes such as water transfer and surface energy
balance. Liao et al. (2019) calibrated parameters related to phenology
and leaf area development in the SAFY simple crop model. They used
phenological stages inferred from the application of the TSF shape
model developed by Sakamoto et al. (2005) to time series of vegetation
indices as observations of phenology. After calibration, they found a
positive impact on biomass and yield simulations even if the biomass to
yield conversion parameters were not calibrated. In the literature,
several works can be found that show the possibility of deriving
phenological stages from canopy reflectance dynamics such as Zhang
et al. (2003), Sakamoto et al. (2005) or Zhou (2019). These results
showed the potential of using remote sensing data for the calibration of
not only parameters related to canopy development, but also parameters
related to crop phenology. In these studies, a good estimate of leaf
development enables a good estimate of yield even if the parameters
enabling biomass and yield conversion are not calibrated. For this, in
situ yield data are required.

The implementation of crop models in an operational context, such
as their use to estimate agricultural production on a territorial scale or
their integration into decision support tools, requires the calibration of
plant parameters using biophysical and phenological data that are not
always available. The objectives of this work are 1) to propose a
simplified methodology for the calibration of plant parameters of the
STICS crop model based on remote sensing to infer the green leaf area
and plant phenology while assuming that crop practices and yield are
available at the farm level and 2) identify the added value of charac-
terizing *in situ* the biomass and plant phenology. The parsimonious use
of data imposed in this work does not allow us to describe the envi-
ronmental conditions that can have an impact on edaphic stresses, such
as water or nutritional stresses. Thus, there may be confusions between
the expression of plant parameters and those of these stresses. To avoid
such confusions, we have carried out the calibration on well irrigated
and fertilised agricultural plots, where we assume that the stresses can
be neglected. We also hypothesize that the effects of water stress are well
captured by the STICS model and associated generic parameters for
wheat. As a third objective 3), this assumption was evaluated using
validation data sets to assess the ability of the calibrated model to
simulate crop development and growth under different stress conditions. The study was implemented on a spring wheat variety (Triticum aestivum L.) grown in the semi-arid conditions of the Central Plateau of Spain.

# 2. Materials and methods

## 2.1. Study sites

A total of 12 commercial irrigated spring wheat fields grown during 2015 and 2016 are included in this study (Fig. 1). The fields are located on the Central Plateau of Spain and define a study area of 3600 km² (39°21' N 2°6' W; 38°40' N 1°32' W) with around 1000 km² of irrigated crops. The climate zone is characterised as Mediterranean with continental effects. The climatology of the area is characterized by the following average yearly values: 14° C for the air temperature, 1242 mm year⁻¹ for the cumulative evapotranspiration and 363 mm year⁻¹ for cumulative precipitation. These numbers were computed over a 20-year period (2002–2021) using the 4 weather stations located in the study area and operated by the Servicio Integral de Asesoramiento al Regante (SIAR) of Spanish Ministerio de Agricultura, Pesca y Alimentación (MAPA) (open data: http://eportal.mapama.gob.es/websiar/SelectionParametrosMap.aspx?dst=1). Soils are calcareous classified as Calcixerepts (Soil Survey Staff, 2014) typical of the Mediterranean region.

## 2.2. STICS crop model

STICS model (Brisson et al., 1998) is a generic crop model which simulates crop development, carbon, water and nitrogen dynamics in the soil–crop–atmosphere system for multiple crops (Guillaume et al., 2011). The version V8.50 (http://www6.paca.inrae.fr/stics_eng) is used in this work. The model is divided into modules that simulate crop development and biochemical and hydrological processes (Brisson and Launay, 2008) for a specie-variety growing under potential conditions that can be reduced according to thermal or nutritional constraints (water or nitrogen stress). It runs daily over a site-year-management unit, termed as the simulation unit, characterised by climate conditions, soil properties and management practices. The modules dealing with crop development, crop growth and yield formation on which this study is focused are briefly described. In order to make the names of the parameters and variables more readable, the names used in the STICS model have been changed. A correspondence table is provided in supplementary material (S5).

Crop development is variety dependent and is driven between phenological stages by the accumulation of daily photothermal units (gddPT, °Cd). gddPT corresponds to daily cardinal degree-days (°Cd) from emergence date ($d_{emerg}$) corrected by the required days for vernalisation $d_{vern}$ (days) and the sensitivity to the photoperiod *sensiphot* (dimensionless). The crop development is defined by the time intervals in gddPT units between wheat phenological stages, hereafter referred to as *pheno-phases*: $st_{embmg}$ (from emergence $d_{emerg}$ date to the beginning of maximal growth date $d_{bmg}$), $st_{bmgfl}$ (from $d_{bmg}$ to flag leaf date $d_{fl}$), $st_{emsgf}$ (from

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Study area: distribution of fields and weather stations considered in this study and Sentinel2 satellite tiles.</figcaption></figure>
$d_{emerg}$ to the onset of grain filling $d_{sgf}$) and $st_{sgfmat}$ (from $d_{sgf}$ to physiological maturity $d_{mat}$).

Crop growth formulation is divided into light interception driven by green leaf area expansion and aerial drymass production. The simulation of green leaf area expansion is divided into the growth and the senescence phases. The leaf expansion rate is simulated by a logistic function driven by plant development and scaled to plant population by $plai_{dens}$ (dimensionless). $plai_{dens}$ determines the ability of a variety to withstand increasing densities. Green leaf area senescence is based on the concept of leaf lifespan (Maas, 1993) that simulates the loss of greenness of leaves once the lifetime has elapsed driven by the variety parameter $plai_{scen}$ (cumulative Q10 unit).

In order to make GAI simulations comparable with the GAI as observed by the satellite (GAI<sub>RS</sub>), the contribution of the green area index of ears (EAI) need to be considered. EAI in STICS is taken into account to simulate the radiation interception while EAI is not considered in GAI computation. Moreover, the progressive onset of ear area and senescence, are not well represented, since the ear area is set to a single value applied between two phenological stages. Therefore, an EAI contribution was added to the leaf GAI computed by STICS in a post-processing step. In the EAI simulations three phases are considered. It is assumed that during the growing phase, green ear surface increases linearly from 0 to a maximum EAI value (eaimax). The growing phase starts at $d_{fl}$ and its duration is parameterized by $stmax$ given in growing degree day (gdd). Then, in the plateau phase period EAI = eaimax, the duration of which is parameterized by $stplat$. Finally, in the senescence phase, which duration is parameterized by $stsen$, EAI decreases linearly from EAI = eaimax to a value EAI = 0. An illustration and details are given in Section S3 of Supplementary Material.

The simulation of growth of aerial drymass, hereafter referred to as biomass, is computed daily based on specie maximum use efficiency of intercepted radiation. Three parameters are determined according to the development phases: $pbiom_{juv}$ (g MJ<sup>-1</sup>) during the juvenile phase (from $d_{emerg}$ to $d_{bmg}$), $pbiom_{veg}$ (g MJ<sup>-1</sup>) for the vegetative phase (from $d_{bmg}$ to $d_{sgf}$) and $pbiom_{rep}$ (g MJ<sup>-1</sup>) during the reproductive phase (from $d_{sgf}$ to $d_{mat}$).

Yield formation is simulated by assuming that a variety can reach a maximum grain weight ($pyield_{wgrmax}$) and develops a maximum number of grains ($pyield_{ngmax}$). Yield formation is represented assuming a progressive carbon allocation from biomass to grains by applying a harvest index ($pyield_{carb}$, weight of grain per weight of biomass per day) during the grain filling period, not exceeding the maximum yield being $pyield_{wgrmax} \cdot pyield_{ngmax}$ product. The number of grains is determined at the time of the onset of grain filling, considering the mean growth rate of the crop during the period preceding the start of grain filling driven by the $pyield_{cgrain}$ (g<sup>-1</sup> m<sup>2</sup>) parameter.

## 2.3. Field data

Field sampling was conducted on 75 one-hectare sampling areas (plots) located in the 12 fields considered in this study (Fig. 1). Three spring wheat varieties (Califa, Atomo, and Galera) were cropped during 2015 and 2016 in the frame of an experiment designed for the evaluation of the impact of the variable rate nitrogen fertilisation management on yield (Campos et al., 2019). In this study, the three varieties were considered as a unique variety since they have similar morphogenetic characteristics classified as short-cycle, semi-early spring wheat of medium-low height, and then a single variety was calibrated. The wheat season covered a period from mid-December (sowing period from 8th December to 22th January) to early July. Plant phenology, aerial biomass and soil properties were characterized on 3 sampling points (a 1.5 m x 2 row areas for vegetation observations) randomly located in each plot. Vegetation characteristics were observed from the 2 leaves stage to physiological maturity. Plot values were obtained from the average of the sampling points. The management of each plot was reported by farmers.

Phenological observations were made during 2015 and 2016 using the BBCH scale (Lancashire et al., 1991) at intervals of 7 – 20 days. Noted that, observations were not taken at the exact phenology stages as required by STICS and very few emergence dates were observed. Moreover, since observations were not frequently observed, some delay is expected between the day of the field visit and the date of the onset of the reported stage. The adaptation approach of the phenological field observations to the phenology stages as required by STICS is described in Section 2.4.

Biomass, was determined from the above-ground biomass field sampled and dried at 80 °C until a constant weight was reached. At maturity, yields were determined by weighing grains. The results were close to that given by the combine harvester ($R^2$ = 0.74 with a bias of 0.87 t ha<sup>-1</sup>, field sampling leading to larger yield) when yield map were collected, which was not the case in all fields. For this reason, we use our ground measurements to mimic the potential use of farmer data as yield monitored by the combine harvester. The difference in yield between the two sources of information can be explained by the yield sensor calibration and/or the sampling strategy that might underestimate the impact of missing rows where plant emergence failed.

Soils was sampled in the following layers: 0 – 20 and 20 – 40 cm depth. Granulometry (Bouyoucos method), calcium carbonate content (Bernard method), soil organic matter content (Walkley-Black method), organic and mineral nitrogen content (Kjeldahl method), pH (pH meter) and depth of surface calcareous horizons when it appears above 40 cm were determined. Soil characteristics are presented in Table 1 and details are given in supplementary material S6. Due to the low slope, we have assumed that runoff is negligible.

Climatic data as precipitation amount, air temperature, relative air humidity, wind velocity and incoming solar radiation were uploaded for the 2015 and 2016 crop seasons from the SIAR's website. The four meteorological sites in the area were considered while the technical characteristics of the sensors can be found in the website using the weather station label given in Fig. 1. We assumed that for each plot the climate was that of the closest meteorological station.

## 2.4. Adaptation of phenological field observations to STICS phenology stages

STICS phenology stages were derived from the phenological field observations (Section 2.3) based on the correspondences presented in Table 2. Instead of using a single BBCH index, a range of BBCH indices

<table><caption>**Table 1**

Pedoclimatic spring wheat growing conditions.</caption><thead><tr><th>Environmental indicator</th><th>Parameter/definition</th><th>Range</th></tr></thead><tbody><tr><td rowspan="10">Soil</td><td>Mean field elevation (m)</td><td>660 – 832</td></tr><tr><td>Mean field slope (%)</td><td>1.0 – 4.2</td></tr><tr><td>Soil clay content (%) (g kg<sup>-1</sup>)</td><td>15 – 36</td></tr><tr><td>Soil texture class (USDA)<sup>2</sup></td><td>SL, L, CL, SCL</td></tr><tr><td>Soil organic content (g kg<sup>-1</sup>)</td><td>19 – 38</td></tr><tr><td>Soil organic nitrogen content (g kg<sup>-1</sup>)</td><td>0.1 – 0.3</td></tr><tr><td>Water amount of field capacity (%w.)</td><td>13.3 – 22.7</td></tr><tr><td>Available water (%w.)</td><td>6.3 – 9.8</td></tr><tr><td rowspan="7">Climate<sup>1</sup></td><td></td><td><strong>2015</strong></td></tr><tr><td>Average daily mean temperature (°C)</td><td>11.9 ± 5.7</td></tr><tr><td>Cumulative rainfall (mm)</td><td>82 ± 13</td></tr><tr><td>Cumulative evaporative demand (mm)</td><td>305 ± 158</td></tr><tr><td>Mean daily global radiation (MJ m<sup>-2</sup> d<sup>-1</sup>)</td><td>18 ± 7</td></tr><tr><td></td><td><strong>2016</strong></td></tr><tr><td></td><td>11.5 ± 4.2</td></tr><tr><td></td><td>76 ± 14</td></tr><tr><td></td><td>281 ± 120</td></tr><tr><td></td><td>17 ± 7</td></tr></tbody></table>

1: arithmetic mean and cumulative variables of daily observations were computed over the wheat growing season; 2: SL: Sandy Loam; L: Loam; CL: Clay Loam; SCL: Sandy Clay Loam;
<caption>**Table 2**</caption>

Description of STICS wheat phenology stages and their correspondence with i) phenological field observations in BBCH scale, and ii) phenological vegetation indices (PVI) derived from remote sensing data.

<table><thead><tr><th rowspan="2">Phase</th><th colspan="2">STICS</th><th colspan="2">BBCH (from field observation)</th><th colspan="2">PVI (from satellite observation)</th></tr><tr><th>Acron.</th><th>Description¹</th><th>Index [range]⁴</th><th>Description²</th><th>Acron.</th><th>Description³</th></tr></thead><tbody><tr><td rowspan="3">Vegetative</td><td>d<sub>emerg</sub></td><td>Emergence</td><td>09<br>[Phylloc.]</td><td>Emergence: coleoptile penetrates soil surface.</td><td rowspan="3">BGR</td><td rowspan="3">Beginning growth rate (BGR) define the time when the quick increase of the GAI growth occur that is associated with the commonly known green-up.</td></tr><tr><td>d<sub>bmg</sub></td><td>Maximum acceleration of leaf growth, end of juvenile phase. It equates to the beginning of stem elongation and is generally not far from the end of ear initiation: it is the "1 cm ear" stage, just slightly later than the double-ridge stage for wheat.</td><td>21<br>[13 – 15,<br>21 – 23]</td><td>Beginning of tillering: first tiller detectable.</td></tr><tr><td>d<sub>fil</sub></td><td>Maximum leaf area index, end of leaf growth.</td><td>31<br>[26 – 33]</td><td>Stem elongation: First node at least 1 cm above tillering node.</td></tr><tr><td rowspan="3">Reproductive</td><td>d<sub>sig</sub></td><td>Onset of filling of harvested organs.</td><td>39<br>[37 – 42]</td><td>Flag leaf stage: flag leaf fully unrolled, ligule just visible.</td><td>MGR</td><td>Maximum growth rate (MGR) define the time of the maximum rate of crop growth.</td></tr><tr><td>d<sub>mat</sub></td><td>Physiological maturity, the harvested organs stop growing in dry matter terms.</td><td>65<br>[63 – 67]</td><td>(6) Flowering. (5) Full flowering, anthesis: 50 % of anthers mature.</td><td>BOP</td><td>Beginning of plateau (BOP) define the time when the slowdown of the GAI growth occur.</td></tr><tr><td>d<sub>gdd</sub></td><td></td><td>87<br>[85 – 89]</td><td>Hard dough: grain content solid. Fingernail impression held.</td><td>MOP</td><td>Mid of plateau (MOP) define the inflection point between the growth and the senescence that is associated with the onset of the leaves senescence.</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>EOS</td><td>End of season (EOS) define the time when GAI reaches the value of zero.</td></tr></tbody></table>

1. Description from Brisson et al. (1998);

2. Description from Lancashire et al. (1991);

3. Description based on the interpretation of the algorithms used to smooth the shape of the GAI<sub>RS</sub> observations by Eq. 1;

4. BBCH index corresponding to each STICS stage (in bold) and the range of BBCH indices around the main index considered (in brackets) for the determination of STICS stages.

(given in brackets in Table 2) was considered around the BBCH index corresponding to each STICS phenological stage (highlighted in bold in Table 2) to enlarge the number of observations that match STICS phenological stages. By considering phenological stages before and after the target stage, bias is minimized at the cost of increased uncertainty, which can be reduced by a larger number of observations. It is important to note that *BBCH31* was also included, although it does not correspond to any STICS phenological stage. It was a well characterized phenological stage that enlarges the population to estimate the photoperiod sensitivity of the variety (refer to Section 2.6 for details).

Emergence dates, not observed during field visits, were inferred from the annotations of the appearance of the first leaves by applying the phyllochron concept. Phyllochron is defined as the thermal time interval between the appearance of two successive leaves (Gate, 1995; Slafer and Rawson, 1997), which is assumed constant up to 8 leaves (Jamieson et al., 1995). The phyllochron can be affected by both genetic and environmental factors (Slafer, 2015). Consequently, a phyllochron value was computed for each field. Values in the range of 70 – 100 growing degree-days (gdd, °Cd) were obtained, with lower values for more intensively managed fields that were unstressed.

## 2.5. Green area index and phenology from remote sensing data

In this work, GAI<sub>RS</sub> derived from satellite observations was used to monitor the growth dynamics of the canopy and the crop phenology. GAI<sub>RS</sub> was chosen among other vegetation indices to match the output of STICS model and to take profit of an existing analytical functions representing GAI dynamic which was used to determine phenological stages. GAI<sub>RS</sub> was computed using the BV-NET algorithm (Weiss and Baret, 1999) applied on multispectral images from Landsat 8 OLI (L8) in 2015 and Sentinel 2 (S2) in 2016. The BV-NET algorithm was implemented with the green, red and near infrared bands (B3, B4, and B5 bands for L8 with 30 m spatial resolution and B3, B4, and B8 bands for S2 with a 10 m spatial resolution). The time interval of acquisition dates varied between 3 days (path overlapping) to more than 1 month (cloudy periods), leading to 6 L8 dates over the wheat season in 2015 and from 8

to 14 S2 dates in 2016 according to the field locations. The satellite observations were regularly acquired between day 25 and 191 covering the crop cycle. The growth dynamics of the wheat canopy were monitored by considering the average GAI<sub>RS</sub> in each of 75 the observation plots leading to average between 9 (L8) and 100 (S2) pixels according to the image spatial resolution.

Crop development was monitored by detecting phenology through the timing of specific temporal features in the seasonal dynamics of GAI<sub>RS</sub>. Therefore, Phenological Vegetation Indices (PVI) were defined corresponding to the times when the features of interest occurred. A seasonal relationship was established by fitting the analytical model proposed by Baret (1986) on GAI<sub>RS</sub> time series. The model combines a logistic function for the growth phase, and an exponential function for the senescence phase, based on daily growing degree-days (gdd, °Cd) (Eq. 1).

$$GAI(gdd) \sim GAI_{\max} \left[ \frac{1}{1 + e^{r_s(gdd_m - gdd)}} - e^{r_s(gdd_s - gdd)} \right] \quad (1)$$

gdd were computed as the daily average air temperature from a base (0 °C) to a maximum (33 °C), which were also values used in STICS for wheat. In Eq. 1, GAI<sub>max</sub> represents the maximum GAI, r<sub>g</sub> is the parameter of the logistic function that describes the relative growth rate at the inflection point obtained when gdd=gdd<sub>m</sub>, and r<sub>s</sub> is the parameter of the exponential function that describes the relative senescence rate when all leaves have senesced (gdd=gdd<sub>s</sub>), i.e. when GAI = 0. A rmse = 0.30 m was obtained between the simulated GAI and GAI<sub>RS</sub> observations, which is comparable to the value reported by Casa et al. (2012). One can note that the impact of remote sensing acquisition dates, which could vary from one site to another or between years, is reduced by fitting an analytic function on GAI<sub>RS</sub> observations.

Subsequently, the following five PVIs were characterized using threshold and analytical methods: Beginning of Growth Rate (BGR), Mid Growth Rate (MGR), Beginning Of Plateau (BOP), Mid Of Plateau (MOP) and End Of Season (EOS) (Fig. 2).

BGR and EOS were derived using threshold methods to detect the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Example of the methodology to derive the phenological vegetation indices (PVI) (dotted) from GAI<sub>RS</sub> dynamics and their 1st (long dashed line) and 2nd (short dashed lines) derivatives.</figcaption></figure>

green-up and the end of the season, near the maturity, respectively. They were defined as the points when daily GAI values reach 10 % of the seasonal maximum, as used in previous studies by Lobell et al. (2013) and González-Gómez et al. (2018). MGR, BOP, and MOP were derived using the properties of analytic curve (Eq. 1) and its first and second derivatives to determine unambiguously critical days. MGR was determined as the time of maximum GAI growth rate, as proposed by Zhang et al. (2003) and Zhou (2019), indicating an intermediate stage during the vegetation growth, around the onset of the stem elongation. BOP was identified as the second local minimum of the second derivative of the daily GAI values, indicating a slowdown in vegetative growth rate after the period of rapid growth and just before vegetation stops growing around the time of the flag leaf stage. MOP was derived as the second local maximum of the second derivative of the daily GAI values, as proposed by Sakamoto et al. (2005). This point indicates the time when the crop reaches its maximum GAI values and undergoes a transition from the vegetative to the reproductive phase around flowering, accompanied by the onset of leaf senescence.

The correspondence between PVIs and STICS phenological stages (Table 2) were established based on the comparison of their time courses as a function of STICS green leaf area (GAI<sub>STICS</sub>). The process is detailed in of supplementary material (S1.1). The correspondence between PVIs and BBCH was done on interpolated BBCH times series as described in (S1.2).

## 2.6. Calibration scenarios

CAL1 is the main scenario of this study, in which we used:

* information from farm (crop management, crop yield and soil properties such as soil texture, soil organic matter and estimates of soil depth and soil mineral nitrogen content)

* existing information on yield potential found in commercial notes and/or trial reports delivered by technical institutes,

* GAI<sub>RS</sub> time series derived from satellite observations

* all crop parameters that were not calibrated were set to that of the existing STICS calibrated variety (Yecora Rojo). It corresponds to a spring wheat variety adapted to a semi-arid environment similar to that found in central Spain. It was calibrated as part of the crop model comparison within the AgMIP exercise, using an extensive database collected from an agronomic trial (Asseng et al., 2015).

In this study, to increase the number of calibration and validation points, we used plot-level averaged yield observations and the corresponding GAI<sub>RS</sub> time series.

In CAL2 we used *in situ* observations of wheat phenology to assess the impact of phenology when derived from remotely sensed PVIs. In CAL3, in addition to *in situ* phenology observations, we include biomass measurements to assess the impact of having an accurate description of the biomass dynamic in the calibration process (Table 3). Finally, the Yecora Rojo plant files was used in the YR scenario to evaluate the performances of an existing variety parameterization determined for a similar cultivar grown in similar climatic conditions but in another geographical context.

## 2.7. STICS calibration

The calibration approach considered in this work, is based on 5 or 6 successive steps. The order of the calibration steps following STICS guide line, which was often followed as for instance in Guillaume et al. (2011). The calibrated parameters and the target output variables in each step are presented in Table 3A total of 14 plant parameters were calibrated, selected for their sensitivity in simulating each respective target variable.

**Step1.** The parameters *d*<sub>vern</sub> and *sensiphot*, which play a role in the
<table><caption>**Table 3**

Calibration scenarios outline.</caption><thead><tr><th rowspan="2">Step</th><th rowspan="2">Process</th><th rowspan="2">Calibrated Parameter</th><th rowspan="2">Target variable</th><th colspan="3">Calibration scenarios</th></tr><tr><th>CAL1</th><th>CAL2</th><th>CAL3</th></tr></thead><tbody><tr><td>1</td><td>Photothermal time units (gdd<sub>PT</sub>)</td><td>d<sub>vern</sub></td><td>Phenology from RS data</td><td>X</td><td></td><td></td></tr><tr><td></td><td></td><td>sensiphot</td><td>Phenology from in-situ data</td><td></td><td>X</td><td>X</td></tr><tr><td rowspan="4">2</td><td rowspan="4">Crop development</td><td>s<sub>tembmg</sub></td><td>Phenology dates from RS data</td><td rowspan="4">X</td><td></td><td></td></tr><tr><td>s<sub>tpmgflb</sub></td><td>Phenology dates from in-situ data</td><td></td><td>X</td><td>X</td></tr><tr><td>s<sub>temsgfs</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>s<sub>zsgfmat</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>Biomass formation</td><td>pbiom<sub>veg</sub></td><td>pbiom<sub>rep</sub></td><td>Biomass</td><td></td><td>X</td></tr><tr><td>4</td><td>Green leaf area formation</td><td>plai<sub>dens</sub></td><td>plai<sub>scen</sub></td><td>GAI<sub>RS</sub></td><td>X</td><td>X</td></tr><tr><td>5</td><td>Potential yield</td><td>pyield<sub>wgrmax</sub></td><td>pyield<sub> agrmax</sub></td><td>Yield</td><td>X</td><td>X</td></tr><tr><td>6</td><td>Yield formation</td><td>pyield<sub>cgrain</sub></td><td>pyield<sub>carb</sub></td><td>Yield</td><td>X</td><td>X</td></tr></tbody></table>

simulation of photothermal units (gdd<sub>PT</sub>), were initially calibrated due to their significant impact on development parameters (step2). Optimal values were selected by minimizing the variability in gdd<sub>PT</sub> at various phenological stages across different climatic conditions of the plots within the calibration data set. The variability was quantified using the coefficient of variation (CV in %) defined as the ratio of the standard deviation divided by the average of gdd<sub>PT</sub> computed for each phenological stage and a given set of photothermal parameters. The d<sub>vern</sub> parameter was calibrated first considering a range of values between 0 and 38 days. The CV was evaluated exclusively at the d<sub>bmg</sub> phenological stage, since vernalization requirements are not expected beyond that stage. The sensiphot parameter was calibrated considering the BBCH31 stage. If BBCH31 stage was not a STICS phenological stage, it was taken as it was known to be sensitive to photoperiod (Slafer 1996 and Ochagavía et al. 2018). The CV was computed with sensiphot ranging from 0.1 to 1 and the values leading to the minimal value was selected.

**Step2.** The *pheno-phases* were derived by summing the daily gdd<sub>PT</sub> values between phenological stages. The computation was done on each plot in the calibration data set and then averaged.

**Step3.** pbiom<sub>veg</sub> and pbiom<sub>rep</sub> parameters were fitted to match simulated biomass with biomass measurements.

**Step4.** plai<sub>dens</sub> and plai<sub>scen</sub> parameters were fitted to match simulated GAI with GAI<sub>RS</sub> observations.

**Step5.** The parameters that define the maximum potential yield of the variety, pyield<sub> agrmax</sub> and pyield<sub> wgrmax</sub> were determined using data collected from official trials conducted by different organisations across Spain (refer to supplementary material S2).

**Step6.** pyield<sub>cgrain</sub> and pyield<sub>carb</sub> parameters were fitted to match simulated yields with measurements.

The fitting of parameters to match STICS simulations with the target observations (step3, step4, and step6) was performed by minimizing the goodness-of-fit criterion using the Nelder-Mead simplex algorithm. This algorithm is used for multidimensional minimisation of any parameters within the STICS crop model (Wallach et al., 2011). It is adapted to calibrate parameters of models with thresholds where gradient methods are not efficient. Moreover, the multistart mode minimizes the risk of reaching local minimum (Wallach et al., 2024).

The soil properties and the initial conditions for STICS simulations were determined from farmer knowledge. Soil profiles were characterized by three layers of 20 cm depth each, in line with the typical soils of the study area, which are mainly shallow, with a calcareous bedrock near the surface. Hydrological properties for each layer were derived using pedotransfer functions (Wösten et al., 1999). Simulations were initialised just after substantial rainfall events, allowing the initialisation of the simulations at field soil capacity, reducing the impact of initial soil conditions.

The plant parameters obtained from each calibration scenario were assessed using prediction errors metrics, namely *rmse* (Root Mean Square Error) and *bias*.

## 2.8. Calibration and validation data sets

To address the uncertainties arising from the inaccurate characterization of soil and initial conditions, specific calibration and validation data sets were established. As a result, plots with the highest yield, well-irrigated and well-fertilised throughout the growing season, were selected for calibration purposes. The plots included in the calibration data set (DBOpti) were identified based on both the maximum GAI<sub>RS</sub> during the whole season and highest yield. A total of 6 plots were selected for DBOpti from fields 9 and 10 (see Table S6). The remaining 63 plots, independent from DBOpti, were employed for evaluating the calibrated plant parameters (DBValid). Furthermore, DBValid was divided into three validation subsets (DBValid1 – 3) also based on GAI<sub>RS</sub> dynamics and yield to represent three different levels of water and nitrogen stress attributed to crop intensification. A comprehensive description of the data sets, encompassing the range of GAI<sub>RS</sub> values, yields, as well as the total nitrogen and water doses applied, is provided in Table 4.

# 3. Results

## 3.1. Phenological vegetation indices (PVI)

The PVIs, used in CAL1 as phenology observations, were evaluated against the BBCH observations to determine their accuracy in predicting STICS phenology stages (Table 2). The prediction errors of each PVI in DBOpti showed a moderate bias of < 5 days, except for BOP, which showed the largest deviation of –9 days (Table 5). BOP, intended to detect the d<sub>PL</sub> stage (flag leaf stage BBCH39), occurred earlier than expected, during stem elongation (BBCH34). This advance can be attributed to the fact that remote sensors introduce uncertainties when GAI<sub>RS</sub> values exceed 6 m<sup>2</sup> m<sup>-2</sup>, owing to the saturation of the remote sensing signals in densely vegetated areas (Weiss and Baret, 2016). This phenomenon accelerates the rapid development phase and leads to the slowdown of the vegetative growth rate before the flag leaf, as observed in DBOpti plots which had dense canopy cover. The accuracy identified

<table><caption>**Table 4**

Calibration (DBOpti) and validation (DBValid) data sets.</caption><thead><tr><th></th><th>DBOpti</th><th>DBValid1</th><th>DBValid2</th><th>DBValid3</th></tr></thead><tbody><tr><td>n</td><td>6</td><td>24</td><td>25</td><td>14</td></tr><tr><td>Sowing date (Julian days<sup>1</sup>)</td><td>358±8</td><td>369±17</td><td>372±16</td><td>359±15</td></tr><tr><td>N fertilisation (kgN ha<sup>-1</sup>)</td><td>304±32</td><td>226±69</td><td>211±56</td><td>137±35</td></tr><tr><td>Irrigation (mm)</td><td>499±4</td><td>440±42</td><td>421±54</td><td>298±71</td></tr><tr><td>Yield (t ha<sup>-1</sup>)</td><td>9.6±2.0</td><td>9.7±1.4</td><td>7.4±1.3</td><td>5.1±0.9</td></tr><tr><td>GAI<sub>RS</sub> (m<sup>2</sup> m<sup>-2</sup>)</td><td>7.6±0.5</td><td>6.8±0.8</td><td>5.6±0.6</td><td>4.7±0.5</td></tr></tbody></table>

1: Julian days of two consecutive years.

n: number of simulation units included in each group; arithmetic mean±standard deviation of each group of the total amount of nitrogen fertilized (N fertilization) and irrigation water (Irrigation) applied during the season, yield field measured (Yield) and maximum GAI<sub>RS</sub> (GAI<sub>RS</sub>).
<caption>**Table 5**</caption>

Evaluation of the PVIs: prediction errors (rmse and bias) of PVIs for detecting STICS wheat phenology stages against the dates derived from phenological field observations and phenological stages detected by PVI in terms of BBCH (arithmetic mean ± standard deviation).

<table><thead><tr><th>PVI (STICS stages)</th><th>BGR (d<sub>bmg</sub>)</th><th>MGR (BBCH31)<sup>1</sup></th><th>BOP (d<sub>fl</sub>)</th><th>MOP (d<sub>sgf</sub>)</th><th>EOS (d<sub>mat</sub>)</th></tr></thead><tbody><tr><td><strong>DBOpti</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><em>rmse</em>, days</td><td>18</td><td>4</td><td>10</td><td>6</td><td>1</td></tr><tr><td><em>bias</em>, days</td><td>4</td><td>1</td><td>-9</td><td>-3</td><td>&lt; 1</td></tr><tr><td><strong>BBCH</strong></td><td>20 ± 4</td><td>30 ± 1</td><td>34 ± 2</td><td>63 ± 4</td><td>89 ± 1</td></tr><tr><td><strong>DBValid1</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><em>rmse</em>, days</td><td>10</td><td>7</td><td>15</td><td>9</td><td>4</td></tr><tr><td><em>bias</em>, days</td><td>3</td><td>-3</td><td>-10</td><td>-1</td><td>-1</td></tr><tr><td><strong>BBCH</strong></td><td>20 ± 4</td><td>29 ± 2</td><td>34 ± 5</td><td>66 ± 8</td><td>88 ± 1</td></tr><tr><td><strong>DBValid2</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><em>rmse</em>, days</td><td>8</td><td>9</td><td>15</td><td>14</td><td>11</td></tr><tr><td><em>bias</em>, days</td><td>4</td><td>3</td><td>-4</td><td>3</td><td>-9</td></tr><tr><td><strong>BBCH</strong></td><td>20 ± 3</td><td>31 ± 4</td><td>40 ± 11</td><td>66 ± 12</td><td>88 ± 1</td></tr><tr><td><strong>DBValid3</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><em>rmse</em>, days</td><td>33</td><td>16</td><td>12</td><td>15</td><td>8</td></tr><tr><td><em>bias</em>, days</td><td>30</td><td>11</td><td>10</td><td>11</td><td>-5</td></tr><tr><td><strong>BBCH</strong></td><td>25 ± 4</td><td>35 ± 6</td><td>45 ± 10</td><td>69 ± 7</td><td>86 ± 2</td></tr></tbody></table>

¹: additional stage to the STICS stages included in this study for the determination of the sensitivity to photoperiod.

in this study is consistent with results presented in Section 2.5, where the PVIs derived from simulated STICS GAI<sub>STICS</sub> time series were compared with simulated phenological stages. Hence, the discrepancy between PVIs and phenology, particularly with BOP, could potentially influence the calibration of plant parameters. This impact will be further evaluated.

The results presented in Table 5 also highlight the importance of paying attention to BGR, as a significant bias was observed in DBValid3, which contains the less irrigated plots. These particular plots exhibited the shortest canopy greenness dynamics throughout the season, which was further influenced by lower crop density. With low leaf development, the fitting of the beginning of the curve to determine BGR might be delayed. This indicates that the relationship between the PVI as determined in this study and phenology requires well developed wheat canopy to be significant.

## 3.2. Calibration of plant parameters

The parameters obtained in each calibration scenario, CAL1, CAL2 and CAL3, are presented in this section.

**Step 1 (Photothermal time unit).** In the case of CAL2, the minimum CV% of accumulated gdd<sub>PT</sub> at the d<sub>bmg</sub> stage, calculated from the plots of DBOpti was obtained with d<sub>vern</sub> values of 0 and 10 days (Table 6). Similarly, with CAL1, the minimum CV% was also achieved with d<sub>vern</sub> = 0 (Table 6. These results indicate that the variety does not require days of vernalization, as expected for spring wheat varieties and are consistent with the values characterized for the Yecora Rojo (YR) variety and with the results reported by Jégo et al. (2010) for spring wheat in Canada. The consistency of results between CAL1 and CAL3 was obtained despite the larger values of CV% attributed to the broader range of dates corresponding to BGR stage.

In all calibration cases, sensiphot = 1 give the best results (Table 6),

<caption>**Table 6**</caption>

gdd<sub>PT</sub> parameterisation: coefficients of variance (CV%) for the tested d<sub>vern</sub> and sensiphot values (in bold) at the considered phenology stages.

<table><thead><tr><th></th><th colspan="4">CAL1</th><th colspan="4">CAL3</th></tr></thead><tbody><tr><td>d<sub>vern</sub></td><td>0</td><td>10</td><td>20</td><td>38</td><td>0</td><td>10</td><td>20</td><td>38</td></tr><tr><td>d<sub>bmg</sub></td><td>48</td><td>53</td><td>61</td><td>69</td><td>17</td><td>15</td><td>16</td><td>30</td></tr><tr><td><strong>Sensiphot</strong></td><td><strong>0.1</strong></td><td><strong>0.4</strong></td><td><strong>0.7</strong></td><td><strong>1.0</strong></td><td><strong>0.1</strong></td><td><strong>0.4</strong></td><td><strong>0.7</strong></td><td><strong>1.0</strong></td></tr><tr><td><strong>BBCH31</strong></td><td>8</td><td>5</td><td>4</td><td>4</td><td>10</td><td>7</td><td>5</td><td>5</td></tr></tbody></table>

suggesting no substantial impact of photoperiod. This value is different to that determined for Yecora Rojo variety.

**Step 2 (Pheno-phases).** The durations of the pheno-phases obtained with CAL2 are significantly different from those with YR (Table 7), due to the different sensiphot values (0.1 for YR while 1 in our study) which have a strong impact on gdd<sub>PT</sub> computation. Similar results were reported by Jégo et al., (2010) for st<sub>embmg</sub> and st<sub>sgfmat</sub>, where values of 235 °Cd and 700 °Cd were reported respectively. With CAL1, the pheno-phases obtained were generally similar to CAL2, except for st<sub>bmgflb</sub> with a reduction in thermal time of 116 °Cd. This shift is the consequence of the bias observed when estimating d<sub>fl</sub> with the BOP PVI (Section 3.1).

**Step 3 (Biomass).** The calibrated values of pbiom<sub>rep</sub> and pbiom<sub>veg</sub> are 3.5 g MJ<sup>-1</sup> and 3.1 g MJ<sup>-1</sup> respectively (Table 7). These values are significantly smaller than the default values from the YR scenario, which used the generic default value proposed by the STICS model for wheat. However, the calibrated values look realistic as they are within the range of radiation use efficiency measurements reported in the literature under optimal growing conditions (2.3 g MJ<sup>-1</sup> in Xu et al., 2024, between 2 and 3.5 g MJ<sup>-1</sup> in Robles-Zazueta et al., 2021 and 4.5 g MJ-1 in Ahmed et al., 2024).

**Step 4 (Green leaf area).** The plai<sub>scen</sub> are consistent between CAL1, CAL2 and YR while a lower value is obtained with CAL3 (Table 7). The common point between CAL1, CAL2 and YR was the absence of biomass parameters calibration, leading to overestimated biomass simulations (Section 3.3). Some compensation might exist between biomass production and leave senescence. Regarding plai<sub>dens</sub>, the impact of biomass parameters (pbiom<sub>rep</sub> and pbiom<sub>veg</sub>) calibration was limited, as the results are similar with CAL2, CAL3 and YR scenarios. In CAL1, a lower plai<sub>dens</sub> value was obtained, which can be attributed to the shorter duration of the vegetative period indicated by the st<sub>emsgf</sub> parameter.

**Step 5 (Potential yield).** The maximum yield potential threshold, parameterised by pyield<sub>wgrmax</sub> and pyield<sub>ngmax</sub>, was increased to 12.5 t ha<sup>-1</sup> to match the maximum yield measurements (Table 7). These revised values were uniformly applied across all calibration scenarios, as they were established based on accessible reports from variety trials, assumed accessible at all levels.

**Step 6 (Yield formation).** The calibrated pyield<sub>cgrain</sub> values are variable between scenarios (Table 7). In STICS, pyield<sub>cgrain</sub> is determined as dependent on biomass during 30 days preceding the onset of grain filling, d<sub>sgf</sub>. Consequently, the lower value obtained in CAL2 in comparison to CAL3 might balance the higher biomass simulations obtained in CAL2. The lower value obtained with CAL2 compared to CAL1 might be the consequence of shifts in pheno phases during the vegetative period (see Step 2) that led to compensations. Regarding pyield<sub>carb</sub> similar values were obtained for all scenarios.

<caption>**Table 7**</caption>

Plant parameters obtained with each calibration scenario (CAL1, CAL2 and CAL3) and the initial parameters of an existing spring wheat in STICS, Yecora Rojo (YR).

<table><thead><tr><th>Step</th><th>Module</th><th>Parameter</th><th>YR</th><th>CAL1</th><th>CAL2</th><th>CAL3</th></tr></thead><tbody><tr><td rowspan="2">1</td><td rowspan="2">Photothermal time units (gdd<sub>PT</sub>)</td><td>d<sub>vern</sub></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td><strong>sensiphot</strong></td><td>0.1</td><td>1.0</td><td>1.0</td><td>1.0</td></tr><tr><td rowspan="5">2</td><td rowspan="5">Crop development</td><td>st<sub>embmg</sub></td><td>100</td><td>229</td><td>237</td><td>237</td></tr><tr><td>st<sub>bmgfl</sub></td><td>294</td><td>419</td><td>535</td><td>535</td></tr><tr><td>st<sub>emsgf</sub></td><td>563</td><td>1113</td><td>1199</td><td>1199</td></tr><tr><td>st<sub>sgfmat</sub></td><td>824</td><td>796</td><td>736</td><td>736</td></tr><tr><td>pbiom<sub>veg</sub></td><td>4.25</td><td>4.25</td><td>4.25</td><td>3.08</td></tr><tr><td rowspan="2">3</td><td rowspan="2">Potential aerial dry biomass</td><td>pbiom<sub>rep</sub></td><td>4.25</td><td>4.25</td><td>4.25</td><td>3.53</td></tr><tr><td>plai<sub>dens</sub></td><td>-0.60</td><td>-0.53</td><td>-0.58</td><td>-0.60</td></tr><tr><td rowspan="2">4</td><td rowspan="2">Potential green leaf area</td><td>plai<sub>scen</sub></td><td>250.0</td><td>249.9</td><td>249.7</td><td>237.7</td></tr><tr><td>pyield<sub>wgrmax</sub></td><td>0.0480</td><td>0.0489</td><td>0.0489</td><td>0.0489</td></tr><tr><td rowspan="2">5</td><td rowspan="2">Potential yield</td><td>pyield<sub>ngmax</sub></td><td>22500</td><td>25562</td><td>25562</td><td>25562</td></tr><tr><td>pyield<sub>cgrain</sub></td><td>0.036</td><td>0.053</td><td>0.040</td><td>0.056</td></tr><tr><td rowspan="2">6</td><td rowspan="2">Yield formation</td><td>pyield<sub>carb</sub></td><td>0.011</td><td>0.010</td><td>0.011</td><td>0.014</td></tr></tbody></table>
### 3.3. Scenarios evaluation

The plant parameters calibrated for the different scenarios were evaluated using prediction error (*rmse*, *bias*) of key variables linked to the phenology, GAI, biomass, and yield. The evaluation was made on both DBOpti and DBValid data sets. The results are presented in Table 8.

**Crop development.** The evaluation on the calibration data set shows errors in the determination of the STICS phenological stages that do not exceed 7 days (Table 8). The results are comparable from one scenario to another with a few small differences if we refer to a given stage. It should be noted that the parameters for the YR variety gave results comparable to those given in the calibrations scenarios, which shows similarities in crop development for varieties adapted to the semi-arid context. It is significant to point out that the errors structure is mainly due to the bias which can be explained by the low variability of the conditions used to calibrate the parameters. For the validation data set, the errors are much larger with a very strong error on the calculation of the day of emergence, irrespective of the scenario. STICS calculates the date of emergence as a function of soil moisture and temperature conditions and involves plant parameters that have not been modified in the scenarios. The error observed reflects the difficulty in accurately representing the interactions between soil physical state and emergence date when crop conditions are not optimal. This is a major difficulty encountered with STICS, which bases the whole crop development on the emergence date. In STICS we have the option of forcing the emergence date, which was done by running an additional scenario based on CAL3. In this case, the *rmse* was divided by two (results not shown) in the validation dataset, particularly at the beginning and end of the crop cycle, with errors comparable to those obtained in the calibration dataset. However, at the flag leaf stage and at the start of grain filling, the errors remain significant (8 and 11 days for $d_{fl}$ and $d_{sgf}$, respectively). This shows that it is difficult to represent various cropping conditions with a single set of parameters.

**Biomass.** The simulations of biomass dynamics throughout the growing season, after calibration of the radiation use efficiency parameters (CAL3), led on DBOpti plots to halve the *rmse*. This improvement is mainly attributed to *bias* reduction from values around 2.4–2.5 t ha⁻¹ with the other scenarios using the default STICS parameters to −0.7 t ha⁻¹ (Table 8, Fig. 3a,d). It is noteworthy from Table 8 that a better characterization of biomass has little impact on other simulated characteristics, such as GAI and yields. In the DBValid plots, the results are comparable to those obtained with DBOpti plots. Although the *rmse* increased slightly, the *bias* was nearly zero, showing that the calibrated

parameters were effective over a wide range of environmental conditions.

**Green leaf area.** The simulations of green leaf area dynamics throughout the growing season using CAL1 and CAL3 on DBOpti plots showed a reduction in both *rmse* and *bias* compared to those obtained with YR parameters (Table 8). When considering the DBValid plots the advantage is less clear on *rmse* while *bias* is reduced when using calibrated parameters. These results are comparable to previous calibration exercises (Prévot et al., 2003; Guillaume et al., 2011).

The comparison between CAL3 simulations and *in situ* observations showed good estimations along the growing season, as presented in Fig. 3e with a scatter plot around the 1:1 line. The temporal dynamics is well reproduced (Fig. 4a).

In the case of DBOpti, the prediction errors with CAL1 were higher than those obtained with CAL2 and CAL3 and comparable to those with YR (Table 8) while the reverse was found with DBValid. The comparison between CAL1 simulations and *in situ* observations indicated a tendency to underestimate GAI, particularly when intermediate values were reached (Fig. 3b). This feature is illustrated in Fig. 4a, which shows an underestimation during the senescence phase. CAL1 simulations showed an overall advancement of the growing season compared to the simulations of the other Calibration scenarios. This might result from the approximation used to determine the phenology and its subsequent impact on STICS *pheno-phase*.

When considering DBValid plots, GAI results showed a reduction of the *bias* across the three calibration scenarios compared to YR, while all scenarios resulted in comparable *rmse* values (Table 8, Fig. 4b). In all cases, STICS performance on GAI was only slightly degraded for validation plots, even when the stress conditions encountered were much more significant than those in the calibrated plots. This is very encouraging for further STICS implementations.

**Yield.** The evaluation of yield formation simulations using CAL1, CAL2, CAL3 on DBOpti yielded similar results to those obtained with YR (Table 8). However, significant differences were found in terms of *bias*, which can be attributed to the increase in the maximum yield potential threshold of the variety under consideration. This increase is highlighted by the strong negative *bias* observed with YR. Fig. 3f shows that the observed yield variability was not well simulated by STICS. The factors explaining such variability may not be fully accounted for in the observed STICS inputs (e.g. local variations in irrigation doses) or modelled by STICS (crop health).

The evaluation of yield simulations on DBValid significantly increased the *rmse* to around 1.7 t ha⁻¹, regardless of the calibration

<table><caption>**Table 8**

Evaluation of the plant parameters: prediction errors (*rmse* and *bias*) of the plant parameters and simulated variables obtained on each calibration scenario: CAL1, CAL2 and CAL3, and with Yecora Rojo (YR), on the calibration (DBOpti) and validation (DBValid) data sets.</caption><thead><tr><th rowspan="2">calibration scenarios</th><th colspan="4">rmse</th><th colspan="4">Bias</th></tr><tr><th>YR</th><th>CAL1</th><th>CAL2</th><th>CAL3</th><th>YR</th><th>CAL1</th><th>CAL2</th><th>CAL3</th></tr></thead><tbody><tr><td><strong>DBOpti</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>d<sub>emerg</sub>, day</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><td>d<sub>bmg</sub>, day</td><td>5</td><td>7</td><td>6</td><td>6</td><td>5</td><td>7</td><td>5</td><td>5</td></tr><tr><td>d<sub>fls</sub>, day</td><td>6</td><td>5</td><td>7</td><td>6</td><td>6</td><td>-3</td><td>6</td><td>5</td></tr><tr><td>d<sub>sgf</sub>, day</td><td>6</td><td>3</td><td>7</td><td>6</td><td>-6</td><td>3</td><td>7</td><td>5</td></tr><tr><td>d<sub>mat</sub>, day</td><td>2</td><td>6</td><td>6</td><td>5</td><td>2</td><td>6</td><td>6</td><td>5</td></tr><tr><td><strong>Biomass</strong>, t ha<sup>-1</sup></td><td>3.5</td><td>3.4</td><td>3.6</td><td>1.8</td><td>2.4</td><td>2.4</td><td>2.5</td><td>-0.7</td></tr><tr><td><strong>GAI<sub>RS</sub></strong>, m<sup>2</sup> m<sup>-2</sup></td><td>1.4</td><td>1.3</td><td>1.1</td><td>1.0</td><td>-0.7</td><td>-0.5</td><td>-0.2</td><td>-0.2</td></tr><tr><td><strong>Yield</strong>, t ha<sup>-1</sup></td><td>1.3</td><td>1.4</td><td>1.4</td><td>1.2</td><td>-1.0</td><td>1.0</td><td>1.0</td><td>0.8</td></tr><tr><td><strong>DBValid</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>d<sub>emerg</sub>, day</td><td>14</td><td>14</td><td>14</td><td>14</td><td>13</td><td>13</td><td>13</td><td>13</td></tr><tr><td>d<sub>bmg</sub>, day</td><td>10</td><td>13</td><td>12</td><td>12</td><td>5</td><td>10</td><td>8</td><td>8</td></tr><tr><td>d<sub>fls</sub>, day</td><td>10</td><td>6</td><td>12</td><td>11</td><td>9</td><td>3</td><td>11</td><td>10</td></tr><tr><td>d<sub>sgf</sub>, day</td><td>6</td><td>10</td><td>13</td><td>12</td><td>-3</td><td>9</td><td>12</td><td>11</td></tr><tr><td>d<sub>mat</sub>, day</td><td>7</td><td>7</td><td>7</td><td>7</td><td>-2</td><td>4</td><td>4</td><td>4</td></tr><tr><td><strong>Biomass</strong>, t ha<sup>-1</sup></td><td>4.1</td><td>4.1</td><td>4.3</td><td>2.3</td><td>2.5</td><td>2.7</td><td>2.7</td><td>-0.1</td></tr><tr><td><strong>GAI<sub>RS</sub></strong>, m<sup>2</sup> m<sup>-2</sup></td><td>1.3</td><td>1.2</td><td>1.3</td><td>1.5</td><td>-0.6</td><td>-0.1</td><td>0.1</td><td>0.3</td></tr><tr><td><strong>Yield</strong>, t ha<sup>-1</sup></td><td>1.9</td><td>1.7</td><td>1.7</td><td>1.7</td><td>1.3</td><td>0.4</td><td>0.7</td><td>0.7</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Evaluation of calibrated model: simulation of biomass (a, d), green leaf area (b, e), and yield (c, f) with CAL3 (a – c) and CAL1 (d – f). The grey solid line is the 1:1 line. Color corresponds to the data set: DBOpti (black), DBValid1 (green), DBValid2 (orange), DBValid3 (violet).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** STICS green leaf area daily simulations and GAI<sub>RS</sub> observations (x) from a) a calibration plot and b) a validation plot.</figcaption></figure>

scenario considered, and these results were better than those of YR. In Fig. 3c and Fig. 3f, the three circled points of DBValid1 data set (comprising intensively cultivated plots) show a strong bias. These points correspond to a fertilisation trial where the final nitrogen dressing was omitted. STICS was very sensitive to such a reduction in N

fertilization, while the actual impacts were moderate. In this case STICS seems to overestimate the nitrogen stress and its impact on yield which may results from the difficulties to determine environmental properties (soil depth, water storage capacity, organic matter, initial conditions) that may mitigate the nitrogen deficit in the actual conditions.
# 4. Discussion

## 4.1. Impact of the assumptions underlying the calibration process

The calibration scenarios were based on several assumptions that need to be discussed. The strongest hypothesis given in the introduction part was to focus the calibration on unstressed crops to mitigate the impact of soil conditions that are generally not monitored in farm context. Soils conditions might generate edaphic stresses leading to errors in plant parameter determination induced by confusion between genetic traits and the impact of stresses on crop development. In this study, the proposed approach was made possible by the availability of fields with intensive irrigation and fertilisation, which enabled optimum conditions to be achieved or approached. Admittedly, this is a rather favourable situation, which cannot be generalised. Overcoming this limitation without having to carry out detailed monitoring of the crops and the water and nitrogen status of the soil remains a major challenge that will have to be addressed in the future. However, we believe that this can be achieved by carefully selecting the calibration years or by working on large datasets covering a wide range of crop conditions. A positive point shown in this study is that the default parameters involved in the impact of edaphic stresses on crop development and yield production seem to give good results with non-optimal growing conditions.

One might have thought that the official agronomic trials on the varieties studied could have been used for the calibration thanks to better documented crops. However, the results of these trials show that the growing conditions were far from the agronomic optimum, with yields always lower than those recorded on our calibration plots. Moreover, such trials focussed on yield elaboration and did not provide soil monitoring as required to characterize the edaphic stresses, which would have been necessary to determine plant parameters.

Technically, in the proposed calibration we made some choices as the correction function used to account for the ears in STICS green leaf area simulation or the use of within-field yield maps to enlarge the range of variability and take profit of more consistent data set with leaf development and yield observed at the same spatial scale.

Green leaf area as derived from remote sensing incorporates all the green contributors of the canopy (Weiss and Baret, 2016), which in the case of wheat, includes different organs like leaves, stems and ears. Ears, which develop above the canopy, contribute significantly to the green area, even though they are not considered in STICS simulations (Weiss et al., 2001). The proposed EAI module is driven by parameters such as *eaimax*, *stmax*, *stplat* and *stsen*, whose values were obtained from measurements by Casals (1996) on a different variety from the one considered in this study. EAI simulations modify the GAI simulations of STICS and may impact the calibration process. The sensitivity analysis performed on the well-developed canopies of the DBOpti database shows that *eaimax* presented the strongest impact on GAI simulations (Figure S4 in S3). Surprisingly, the best results were obtained when the EAI correction was minimised or absent. The STICS simulations were almost unbiased when the addition of the ears was not considered, *eaimax* = 0. It can be inferred from this result that the integration of the EAI module into STICS is not necessary to reproduce the *GAI*<sub>RS</sub> as observed by remote sensors. However, this might be a consequence of dense vegetation cover with a saturation of the remote sensing signal, as reported by Weiss and Baret (2016) who showed that GAI product derived from Sentinel 2 has uncertainties with GAI values > 6 m<sup>2</sup> m<sup>-2</sup>. Additionally, Li et al. (2015) reported that the spectral reflectance of the canopy in the near-infrared (760–1150 nm) showed a limited response to spike removal under conditions of excessive nitrogen content. This is consistent with the selection criteria applied to the DBOpti plots, which had highest nitrogen fertilisation and irrigation leading to large *GAI*<sub>RS</sub>. *GAI*<sub>RS</sub> dynamics of the DBOpti plots reached their maximum at the flag leaf with values of around 7.5 m<sup>2</sup> m<sup>-2</sup>, and remained constant and insensitive to full ear development. (Figure S5 in S3). In contrast, with the DBValid plots when the maximum *GAI*<sub>RS</sub> ≤ 6.5 m, an improvement

in GAI simulation was obtained. Consequently, the EAI module is not necessary when *GAI*<sub>RS</sub> values are over 6.5 m<sup>2</sup> m<sup>-2</sup>, while it seems to have a positive impact when *GAI*<sub>RS</sub> values < 6.5 m. For the calibration process, which is done on well-developed crops, hypothesis made on ear contribution to GAI had little impact on the results. The calibrated parameters were comparable with variations of about 1 %, which are considerably smaller than the differences observed between the CAL1 and CAL3 scenarios (Table S5 in S3).

In this work we considered within field spatial variation in yield at the spatial resolution of *GAI*<sub>RS</sub> observation using local measurement. In operational context, yield map recorded by the combine harvester could be used. However, such a map might not be available leading to consider average field yield as an alternative, which might have an impact on the calibration process. In our case this led to reduce the considered yield from 9.0 – 11.8 t ha<sup>-1</sup> (plots) to 10.0 – 10.6 t ha<sup>-1</sup> (field average). Considering such average yields in Step 6 of the calibration process following the CAL3 scenario procedure, we got slightly lower values for *pyield*<sub>cgrain</sub>, 0.053 (instead of 0.056 with CAL3), while no change was observed for *pyield*<sub>carb</sub>. The new set of parameters had little impact on prediction errors for yield simulations in DBValid with a rmse of 1.6 t ha<sup>-1</sup> and a bias of 0.6 t ha<sup>-1</sup>, comparable to those obtained with CAL3, 1.7 t ha<sup>-1</sup> and 0.7 t ha<sup>-1</sup> respectively (Table 8). Therefore, it can be concluded that the average field yield commonly known by farmers can be used in the calibration process, making the adaptation of STICS to a variety more accessible.

## 4.2. Equifinality of CAL1 and CAL3 plant parameters

The calibration strategies yielded notably different plant parameter sets while simulation performances were comparable. This brings up the concept of equifinality, where multiple parameter combinations can produce similar model outputs. One may wonder to what extent STICS simulations might diverge under climatic conditions other than those considered in this study. Addressing this question is of great importance for climate change studies. The aim is to determine whether the differences in plant parameters induced by the calibration strategy generate comparable results whatever the climatic conditions, or whether, on the contrary, the simulations can diverge significantly under certain climatic conditions. For this purpose, the discrepancies between simulations using CAL1 and CAL3 plant parameters were examined for the climate period 2002–2020 (see S4). On average, the differences in GAI simulations showed little bias (-0.11 m<sup>2</sup> m<sup>-2</sup>) and a difference that remain lower than the prediction errors computed against observed *GAI*<sub>RS</sub> (Table 8). The interannual variability is small compared to the average differences (Figure S6 in S4), thus both plant files can be used to study the behaviour of GAI under different climatic conditions.

Regarding the simulations of yield, the differences using CAL1 and CAL3 parameters showed comparable simulation results in average than those obtained during the calibration processes. However, the interannual variability of bias between the two sets of simulations (Figure S6, S4) is large, with an average value of -0.45 t ha<sup>-1</sup> and yearly variation ranging from -1 t ha<sup>-1</sup> in 2006–0 in 2013. The processes that led to such differences are complex and difficult to analyse. The outlier as in 2013 can be related to anomaly in temperature during the grain filling period. Such a result leads to recommend a careful analysis of the range of validity of the plant parameters that may be impacted by the calibration strategy. This is essential for climate change studies.

# 5. Conclusions

The methodology proposed for the calibration of a wheat variety in STICS, based on farmers’ knowledge and remote sensing information, gave satisfactory results with accuracy comparable to earlier studies and better than using existing set of plant parameters. The advantage of using remote sensing data is twofold: on the one hand, it overcomes the difficulties of access to data and, on the other hand, it ensures
consistency between the calibration data and the data that can then be
used in assimilation schemes, thus minimising the effects of differences
in determination methods, such as in the determination of leaf area data
(remote sensing versus in situ). To mitigate the impact of inaccurately
described environmental crop growing conditions on the calibrated
parameters, such as soil properties like soil water content, which could
lead to the simulation of uncontrolled stress conditions, we performed
the calibration on well-irrigated and well-fertilised plots. This was
possible in the area of *La Macha* (Spain), where we found intensively
irrigated wheat fields. This is the main limitation of this work and, as
mentioned in the discussion, new calibration strategies should be pro-
posed to take into account limiting environmental factors while
remaining with a limited number of field observations.

The implementation of STICS using the calibrated plant parameters was successfully performed under different environmental conditions spanning a wide range of water and nitrogen stresses. This means that the impact of stress on plant development and yield is well addressed by the STICS formalism once the plant parameters are calibrated. The obtained accuracy (rmse) in GAI ranged between 1.2 and 1.5 m² m⁻² and in yield between 1.7 and 1.8 t ha⁻¹, across the calibration scenarios considered. These error values include STICS model errors, the uncertainties of uncalibrated parameters set to default values and errors of environmental crop growing conditions (soil and climate). These results represent an improvement compared to the simulations obtained with the Yecora Rojo plant parameters, which is the closest variety already calibrated using an extensive field data set. In this study, we therefore show that it is possible to improve crop model simulation by calibrating the varietal parameters for wheat using data that are easily accessible in an operational context. This strengthens the possibility of using crop models in decision-support tools for crop management (irrigation, fertilization) or in territorial monitoring of agricultural production.

Some technical issues concerning phenology, biomass and GAI were
also considered. Firstly, we compared two calibration scenarios, with
and without field phenological observations. When phenological ob-
servations were not available, we proposed to use PVIs that were derived
from analytical properties of the GAI<sub>RS</sub> time series and were therefore
unambiguously defined. Some shifts were observed compared to the
expected STICS key phenology stages, leading to differences in the
calibrated parameters. However, the accuracy of the simulations of yield
and GAI were comparable, being also consistent between them over a
wide range of climatic conditions encountered over the last 20 years.
Regarding biomass, while the calibration of the parameter governing
biomass dynamics proved necessary to accurately simulate biomass
(strong biases were found when it was not included in the calibration
process), not including it had low impact on the other variables such as
GAI, yield and phenology stages. Finally, we considered an empirical
module for estimating the EAI. We showed that with well-developed
crop canopies, such as those used for calibration, the inclusion of the
EAI module had little impact on the GAI simulations and the calibration
process. However, when canopies were less developed (maximum GAI
<6.5 m² m⁻²), the inclusion of the ears in GAI simulations improved the
agreement between STICS and remote sensing observations. This is
essential for assimilating remotely sensed data into crop models, high-
lighting the importance of obtaining data sets on the contribution of
ears, which remain scarce at present.

CRediT authorship contribution statement

**André Chanzy:** Writing – original draft, Validation, Supervision, Methodology, Formal analysis, Conceptualization. **Elena Pareja-Serrano:** Writing – original draft, Validation, Methodology, Formal analysis, Data curation, Conceptualization. **Jose González-Piqueras:** Writing – review & editing, Supervision, Methodology, Funding acquisition, Data curation, Conceptualization.

**Declaration of Competing Interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was performed in the framework of the projects HER-
MANA (HERramientas para el MAneko sostenible de fertilización Nitro-
genada y Agua, funded by the Spanish Ministry Science and Innovation,
Grant No. AGL2015-68700-R) and FATIMA (FArming Tools for external
nutrient Inputs and water MAnagement, European Union's Horizon 2020
Research and Innovation Programme, Grant No. 633945). We would like
to thank the researchers and technicians of the EMMAH unit of INRAE
(Institut National de Recherche pour l'Agriculture l'alimentation et
l'Environnement, Avignon, France) for their support and advice on
STICS and for providing us the satellite products specially to Samuel Buis
and Guillaume Pouget. We also thank all members of the department of
"Teledetección y SIG" of IDR-UCLM (Instituto de Desarrollo Regional-
Universidad de Castilla-La Mancha), ITAP (Instituto Técnico Agro-
nómico Provincial) and AgriSat Iberia S.L. who participated intensively
in the field measurements and the farmers who let us perform the ex-
periments on their fields. We thanks reviewers for their kind and wise
suggestions that led to improve the text.

**Appendix A. Supporting information**

Supplementary data associated with this article can be found in the
online version at doi:10.1016/j.eja.2025.127756.

Data Availability

Data will be made available on request.

References

Abrahamsen, P., Hansen, S., 2000. Daisy: an open soil-crop-atmosphere system model. Environ. Model. Softw. 15, 313–330. https://doi.org/10.1016/S1364-8152(00)00003-7.

Ahmed, M., Bilal, M., Ahmad, S., 2024. Simulation of source sink partitioning in wheat under varying nitrogen regimes using DSSAT-CERES-wheat model. Agric. Water Manag. 303, 109028. https://doi.org/10.1016/j.agwat.2024.109028.

Asseng, S., Ewert, F., Rosenzweig, C., Jones, J.W., Hatfield, J.L., Ruane, A.C., Boote, K.J., Thorburn, P.J., Rötter, R.P., Cammarano, D., Brisson, N., Basso, B., Martre, P., Aggarwal, P.K., Angulo, C., Bertuzzi, P., Biernath, C., Challinor, A.J., Doltra, J., Gayler, S., Goldberg, R., Grant, R., Heng, L., Hooker, J., Hunt, L.A., Ingwersen, J., Izaurralde, R.C., Kersebaum, K.C., Müller, C., Naresh Kumar, S., Nendel, C., O'Leary, G., Olesen, J.E., Osborne, T.M., Palosuo, T., Priesack, E., Ripoche, D., Semenov, M.A., Shcherbak, I., Steduto, P., Stöckle, C., Stratonovitch, P., Streck, T., Supit, I., Tao, F., Travasso, M., Waha, K., Wallach, D., White, J.W., Williams, J.R., Wolf, J., 2013. Uncertainty in simulating wheat yields under climate change. Nat. Clim. Change 3, 827–832. https://doi.org/10.1038/nclimate1916.

Asseng, S., Ewert, F., Martre, P., Rötter, R.P., Lobell, D.B., Cammarano, D., Kimball, B.A., Ottman, M.J., Wall, G.W., White, J.W., Reynolds, M.P., Alderman, P.D., Prasad, P.V. V., Aggarwal, P.K., Anothai, J., Basso, B., Biernath, C., Challinor, A.J., De Sanctis, G., Doltra, J., Fereres, E., Garcia-Vila, M., Gayler, S., Hoogenboom, G., Hunt, L.A., Izaurralde, R.C., Jabloun, M., Jones, C.D., Kersebaum, K.C., Koehler, A.-K., Müller, C., Naresh Kumar, S., Nendel, C., O'Leary, G., Olesen, J.E., Palosuo, T., Priesack, E., Eyshi Rezaei, E., Ruane, A.C., Semenov, M.A., Shcherbak, I., Stöckle, C., Stratonovitch, P., Streck, T., Supit, I., Tao, F., Thorburn, P.J., Waha, K., Wang, E., Wallach, D., Wolf, J., Zhao, Z., Zhu, Y., 2015. Rising temperatures reduce global wheat production. Nat. Clim. Change 5, 143–147. https://doi.org/10.1038/
nclimate2470.

Attia, A., Rajan, N., Xue, Q., Nair, S., Ibrahim, A., Hays, D., 2016. Application of DSSAT-CERES-Wheat model to simulate winter wheat response to irrigation management in the Texas High Plains. Agric. Water Manag. 165, 50–60. https://doi.org/10.1016/j.agwat.2015.11.002.

Bai, Y., Gao, J., 2021. Optimization of the nitrogen fertilizer schedule of maize under drip irrigation in Jilin, China, based on DSSAT and GA. Agric. Water Manag. 244, 106555. https://doi.org/10.1016/j.agwat.2020.106555.

Baret, F., 1986. Contribution au suivi radiometrique de cultures de cereales. Université de Paris-Sud, Paris.
Bécel, C., Munier-Jolain, N.M., Nicolardot, B., 2015. Assessing nitrate leaching in cropping systems based on integrated weed management using the STICS soil-crop model. Eur. J. Agron. 62, 46–54. https://doi.org/10.1016/j.eja.2014.09.005.

Boogaard, H.L., de Wit, A.J.W., te Roller, J.A., Van Diepen, C.A., 2014. WOFOST CONTROL CENTRE 2.1; User's guide for the WOFOST CONTROL CENTRE 2.1 and the crop growth simulation model WOFOST 7.1.7. Wageningen University & Research Centre, Wageningen (Netherlands), Alterna.

Brisson, N., Launay, M., 2008. Conceptual Basis, Formalisations and Parameterization of the Stics Crop Model. Editions Quae, Paris.

Brisson, N., Mary, B., Ripoche, D., Jeuffroy, M.H., Ruget, F., Nicoulaull, B., Gate, P., Devienne-Barret, F., Antonioletti, R., Durr, C., Richard, G., Beaudoin, N., Recous, S., Tayot, X., Plenet, D., Cellier, P., Machet, J.-M., Meynard, J.M., Delécolle, R., 1998. STICS: a generic model for the simulation of crops and their water and nitrogen balances. 1. Theory aramingerization Appl. wheat Corn. Agron. 18, 311–346.

Campos, I., González-Gómez, L., Villodre, J., Calera, M., Campoy, J., Jiménez, N., Plaza, C., Sánchez-Prieto, S., Calera, A., 2019. Mapping within-field variability in wheat yield and biomass using remote sensing vegetation indices. Precis. Agric. 20, 214–236. https://doi.org/10.1007/s11119-018-9596-z.

Casa, R., Varella, H., Buis, S., Guérif, M., De Solan, B., Baret, F., 2012. Forcing a wheat crop model with LAI data to access agronomic variables: evaluation of the impact of model and LAI uncertainties and comparison with an empirical approach. Eur. J. Agron. 37, 1–10. https://doi.org/10.1016/j.eja.2011.09.004.

Casals, M.L., 1996. Introduction des mécanismes de résistance à la sécheresse du Blé dur au fonctionnement phénologique et trophique de la plante dans un modèle dynamique de croissance. (Thèse de doctorat INA-PG).

Ewert, F., Rötter, R.P., Bindi, M., Webber, H., Trnka, M., Kersebaum, K.C., Olesen, J.E., van Ittersum, M.K., Janssen, S., Rivington, M., Semenov, M.A., Wallach, D., Porter, J. R., Stewart, D., Verhagen, J., Gaiser, T., Palosuo, T., Tao, F., Nendel, C., Roggero, P. P., Bartosová, L., Asseng, S., 2015. Crop modelling for integrated assessment of risk to food production from climate change. Environ. Model. Softw. 72, 287–303. https://doi.org/10.1016/j.envsoft.2014.12.003.

Falconier, G.N., Journet, E.-P., Bedoussac, L., Vermue, A., Chlédowski, F., Beaudoin, N., Justes, E., 2019. Calibration and evaluation of the STICS soil-crop model for faba bean to explain variability in yield and N2 fixation. Eur. J. Agron. 104, 63–77. https://doi.org/10.1016/j.eja.2019.01.001.

Fattori Junior, I.M., dos Santos Vianna, M., Marin, F.R., 2022. Assimilating leaf area index data into a sugarcane process-based crop model for improving yield estimation. Eur. J. Agron. 136, 126501. https://doi.org/10.1016/j.eja.2022.126501.

Gate, P., 1995. Ecophysiologie du blé: de la plante à la culture. Paris: Technique et Doc.-Lavoisier, DL 1995, France.

González-Gómez, L., Campos, I., Calera, A., 2018. Use of different temporal scales to monitor phenology and its relationship with temporal evolution of normalized difference vegetation index in wheat. J. Appl. Remote Sens. 12, 1–18. https://doi.org/10.11171/1.JRS.12.026010.

Guérif, M., Duke, C.L., 2000. Adjustment procedures of a crop model to the site specific characteristics of soil and crop using remote sensing data assimilation. Agric. Ecosyst. Environ. 81, 57–69. https://doi.org/10.1016/S0167-8809(00)00168-7.

Guillaume, S., Bergez, J.-E., Wallach, D., Justes, E., 2011. Methodological comparison of calibration procedures for durum wheat parameters in the STICS model. Eur. J. Agron. 35, 115–126. https://doi.org/10.1016/j.eja.2011.05.003.

Hadria, R., Khabba, S., Lahrouni, A., Duchemin, B., Chehbouni, A., Carriou, J., Ouzine, L., 2007. Calibration valid. stetics crop model managing wheat Irrig. Semiarid Marrakech/al Haouz plain 16.

Hao, S., Ryu, D., Western, A., Perry, E., Bogena, H., Franssen, H.J.H., 2021. Performance of a wheat yield prediction model and factors influencing the performance: a review and meta-analysis. Agric. Syst. 194, 103278. https://doi.org/10.1016/j.agsy.2021.103278.

He, D., Wang, E., Wang, J., Lilley, J.M., 2017. Genotype × environment × management interactions of canola across China: a simulation study. Agric. For. Meteorol. 247, 424–433. https://doi.org/10.1016/j.agrformet.2017.08.027.

Huang, Jingfeng, Tang, Shuchuan, Ousama, Abou-Ismail, Wang, Renchao, 2001. Integration of remote sensing data and simulation model to estimate rice yield. In: 2001 International Conferences on Info-Tech and Info-Net. Proceedings (Cat. No.01EX479). Presented at the 2001 International Conferences on Info-tech and Info-net. IEEE, Beijing, China, pp. 101–107. https://doi.org/10.1109/ICII.2001.982729.

IPCC (Ed.), 2014. Climate change 2014: synthesis report. Contribution of Winking Groups I, II, and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Core Writing Team, R.K. Pachauri and L.A. Meyer (eds.). ed. IPCC, Geneva, Switzerland.

Jamieson, P.D., Brookling, I.R., Porter, J.R., Wilson, D.R., 1995. Prediction of leaf appearance in wheat: a question of temperature. Field Crops Res. 41, 35–44. https://doi.org/10.1016/0378-4290(94)00102-1.

Jégo, G., Pattey, E., Bourgeois, G., Morrison, M.J., Drury, C.F., Tremblay, N., Tremblay, G., 2010. Calibration and performance evaluation of soybean and spring wheat cultivars using the STICS crop model in Eastern Canada. Field Crops Res. 117, 183–196. https://doi.org/10.1016/j.fcr.2010.03.008.

Jégo, G., Pattey, E., Bourgeois, G., Drury, C.F., Tremblay, N., 2011. Evaluation of the STICS crop growth model with maize yield parameters calibrated for Eastern Canada. Agron. Sust. Dev. 31, 557–570. https://doi.org/10.1007/s13593-011-0014-4.

Jégo, G., Pattey, E., Liu, J., 2012. Using Leaf Area Index, retrieved from optical imagery, in the STICS crop model for predicting yield and biomass of field crops. Field Crops Res. 131, 63–74. https://doi.org/10.1016/j.fcr.2012.02.012.

Jin, M., Liu, X., Wu, L., Liu, M., 2015. An improved assimilation method with stress factors incorporated in the WOFOST model for the efficient assessment of heavy

metal stress levels in rice. Int. J. Appl. Earth Obs. Geoinf. 41, 118–129. https://doi.org/10.1016/j.jag.2015.04.023.

Jin, X., Kumar, L., Li, Z., Feng, H., Xu, X., Yang, G., Wang, J., 2018. A review of data assimilation of remote sensing and crop models. Eur. J. Agron. 92, 141–152. https://doi.org/10.1016/j.eja.2017.11.002.

Jones, J.W., Hoogenboom, B., Porter, C.H., Boote, K.J., Batchelor, W.D., Hunt, L.A., Wilkens, P.W., Singh, U., Gijsman, A.J., Ritchie, J.T., 2003. The DSSAT cropping system model. Eur. J. Agron. 18, 235–265. https://doi.org/10.1016/S1161-0301 (02)00107-7.

Justino, L.F., Battisti, R., Stone, L.F., Heinemann, A.B., 2022. In silico assessment of sowing dates and nitrogen management in common bean crops. Eur. J. Agron. 133, 126434. https://doi.org/10.1016/j.eja.2021.126434.

Kadiyala, M.D.M., Jones, J.W., Mylavarapu, R.S., Li, Y.C., Reddy, M.D., 2015. Identifying irrigation and nitrogen best management practices for aerobic rice–maize cropping system for semi-arid tropics using CERES-rice and maize models. Agric. Water Manag. 149, 23–32. https://doi.org/10.1016/j.agwat.2014.10.019.

Kassie, B.T., 2016. Performance of DSSAT-Nwheat across a wide range of current and future growing conditions 10. https://doi.org/10.1016/j.eja.2016.08.012.

Katerji, N., Campi, P., Mastrorilli, M., 2013. Productivity, evapotranspiration, and water use efficiency of corn and tomato crops simulated by AquaCrop under contrasting water stress conditions in the Mediterranean region. Agric. Water Manag. 130, 14–26. https://doi.org/10.1016/j.agwat.2013.08.005.

Keating, B.A., Carberry, P.S., Hammer, G.L., Probert, M.E., Robertson, M.J., Holzworth, D., Huth, N.I., Hargreaves, J.N.G., Meinke, H., Hochman, Z., McLean, G., Verburg, K., Snow, V., Dimes, J.P., Silburn, M., Wang, E., Brown, S., Bristow, K.L., Asseng, S., Chapman, S., McCown, R.L., Freebairn, D.M., Smith, C.J., 2003. An overview of APSIM, a model designed for farming systems simulation. Eur. J. Agron. 18, 267–288. https://doi.org/10.1016/S1161-0301(02)00108-9.

Lancashire, P.D., Bleiholder, H., Boom, T.V.D., Langelüdeke, P., Stauss, R., Weber, E., Witzenberger, A., 1991. A uniform decimal code for growth stages of crops and weeds. Ann. Appl. Biol. 119, 561–601. https://doi.org/10.1111/j.1744-7348.1991.tb04895.x.

Li, H., Zhao, C., Yang, G., Feng, H., 2015. Variations in crop variables within wheat canopies and responses of canopy spectral characteristics and derived vegetation indices to different vertical leaf layers and spikes. Remote Sens. Environ. 169, 358–374. https://doi.org/10.1016/j.rse.2015.08.021.

Liao, C., Wang, J., Dong, T., Shang, J., Liu, J., Song, Y., 2019. Using spatio-temporal fusion of Landsat-8 and MODIS data to derive phenology, biomass and yield estimates for corn and soybean. Sci. Total Environ. 650, 1707–1721. https://doi.org/10.1016/j.scitotenv.2018.09.308.

Linker, R., Ioslovich, I., 2017. Assimilation of canopy cover and biomass measurements in the crop model AquaCrop. Biosyst. En. 162, 57–66. https://doi.org/10.1016/j.biosystemseng.2017.08.003.

Liu, L., Cao, R., Chen, J., Shen, M., Wang, S., Zhou, J., He, B., 2022. Detecting crop phenology from vegetation index time-series data by improved shape model fitting in each phenological stage. Remote Sens. Environ. 277, 113060. https://doi.org/10.1016/j.rse.2022.113060.

Lobell, D.B., Ortiz-Monasterio, J.I., Sibley, A.M., Sohu, V.S., 2013. Satellite detection of earlier wheat sowing in India and implications for yield trends. Agric. Syst. 115, 137–143. https://doi.org/10.1016/j.agsy.2012.09.003.

Maas, S.J., 1988. Using satellite data to improve model estimates of crop yield. Agron. J. 80, 655–662. https://doi.org/10.2134/agronj.1988.00021962008000040002.

Maas, S.J., 1993. Parameterized model of gramineous crop growth: I. Leaf Area and dry mass simulation. Agron. J. 85, 348–353. https://doi.org/10.2134/agronj1993.00021962008500020034x.

Manivasagam, V.S., Rozenstein, O., 2020. Practices for upscaling crop simulation models from field scale to large regions. Comput. Electron. Agric. 175, 105554. https://doi.org/10.1016/j.compag.2020.105554.

Meroni, M., d'Andriment, R., Vrieling, A., Fasbender, D., Lemoine, G., Rembold, F., Seguini, L., Verhegghen, A., 2021. Comparing land surface phenology of major European crops as derived from SAR and multispectral data of Sentinel-1 and -2. Remote Sens. Environ. 253, 112232. https://doi.org/10.1016/j.rse.2020.112232.

Nearing, G.S., Crow, W.T., Thorp, K.R., Moran, M.S., Reichle, R.H., Gupta, H.V., 2012. Assimilating remote sensing observations of leaf area index and soil moisture for wheat yield estimates: an observing system simulation experiment: ENKF AND SIRF OSSES FOR WHEAT YIELD. Water Resour. Res. 48. https://doi.org/10.1029/2011WR011420.

Nendel, C., Berg, M., Kersebaum, K.C., Mirschel, W., Specka, X., Wegehenkel, M., Wenkel, K.O., Wieland, R., 2011. The MONICA model: Testing predictability for crop growth, soil moisture and nitrogen dynamics. Ecol. Model. 222, 1614–1625. https://doi.org/10.1016/j.ecolmodel.2011.02.018.

Nendel, C., Kersebaum, K.C., Mirschel, W., Wenkel, K.O., 2014. Testing farm management options as climate change adaptation strategies using the MONICA model. Eur. J. Agron. 52, 47–56. https://doi.org/10.1016/j.eja.2012.09.005.

Ngwira, A.R., Aune, J.B., Thierfelder, C., 2014. DSSAT modelling of conservation agriculture maize response to climate change in Malawi. Soil Tillage Res. 143, 85–94. https://doi.org/10.1016/j.still.2014.05.003.

Ochagavía, H., Prieto, P., Savin, R., Griffiths, S., Safer, G., 2018. Dynamics of leaf and spikelet primordia initiation in wheat as affected by PdP-1a alleles under field conditions. J. Exp. Bot. 69, 2621–2631. https://doi.org/10.1093/jxb/ery104.

Prévot, L., Chauki, H., Trouffleau, D., Weiss, M., Baret, F., Brisson, N., 2003. Assimilating optical and radar data into the STICS crop model for wheat. Agronomie 23, 297–303. https://doi.org/10.1051/agro:200303.

Rahimi-Moghaddam, S., Deihimfard, R., Azizi, K., Roostaei, M., 2021. Characterizing spatial and temporal trends in drought patterns of rainfed wheat (Triticum aestivum
L.) across various climatic conditions: a modelling approach. Eur. J. Agron. 129, 126333. https://doi.org/10.1016/j.eja.2021.126333.

Razzaghi, F., Zhou, Z., Andersen, M.N., Plauborg, F., 2017. Simulation of potato yield in temperate condition by the AquaCrop model. Agric. Water Manag. 191, 113–123. https://doi.org/10.1016/j.agwat.2017.06.008.

Ritchie, J., Otter, S., 1985. Descr. Perform. CERESWheat A Use Oriented Wheat Yield Model 38, 12.

Robles-Zazueta, C.A., Molero, G., Pinto, F., Foulkes, M.J., Reynolds, M.P., Murchie, E.H., 2021. Field-based remote sensing models predict radiation use efficiency in wheat. J. Exp. Bot. 72 (10), 3756–3773. https://doi.org/10.1093/jxb/erab115.

Rodriguez, J.C., Duchemin, B., Hadria, R., Watts, C., Garatuza, J., Chehbouni, A., Khabba, S., Boulet, G., Palacios, E., Lahrouni, A., 2004. Wheat yield estimation using remote sensing and the STICS model in the semiarid Yaqui valley, Mexico. Agronomie 24, 295–304. https://doi.org/10.1051/agro:2004037.

Rosenzweig, C., Jones, J.W., Hatfield, J.L., Ruane, A.C., Boote, K.J., Thorburn, P., Antle, J.M., Nelson, G.C., Porter, C., Janssen, S., Asseng, S., Basso, B., Ewert, F., Wallach, D., Baigorria, G., Winter, J.M., 2013. The Agricultural Model Intercomparison and Improvement Project (AgMIP): protocols and pilot studies. Agric. For. Meteorol. 170, 166–182. https://doi.org/10.1016/j.agrformet.2012.09.011.

Rötter, R.P., Tao, F., Höhn, J.G., Palosuo, T., 2015. Use of crop simulation modelling to aid ideotype design of future cereal cultivars. EXBOTJ 66, 3463–3476. https://doi.org/10.1093/jxb/erv098.

Sakamoto, T., Yokozawa, M., Toritani, H., Shibayama, M., Ishitsuka, N., Ohno, H., 2005. A crop phenology detection method using time-series MODIS data. Remote Sens. Environ. 96, 366–374. https://doi.org/10.1016/j.rse.2005.03.008.

Seyoum, S., Rachaputi, R., Chauhan, Y., Prasanna, B., Fekybelu, S., 2018. Application of the APSIM model to exploit G × E × M interactions for maize improvement in Ethiopia. Field Crops Res. 217, 113–124. https://doi.org/10.1016/j.fcr.2017.12.012.

Slafer, G.A., 1996. Differences in phasic development rate amongst wheat cultivars independent of responses to photoperiod and vernalization. A viewpoint of the intrinsic earliness hypothesis. J. Agric. Sci. 126, 403–419. https://doi.org/10.1017/S0021859600075493.

Slafer, G.A., 2015. Genetic and environmental effects on crop development determining adaptation and yield. Crop Physiol.

Slafer, G.A., Rawson, H.M., 1997. Phyllochron in wheat as affected by photoperiod under two temperature regimes. Funct. Plant Biol. 24, 151. https://doi.org/10.1071/PP96021.

Soil Survey Staff, 2014. Keys to Soil Taxonomy. 12th edn.

Steduto, P., Hsiao, T.C., Raes, D., Fereres, E., 2009. AquaCrop—The FAO Crop Model to Simulate Yield Response to Water: I. Concepts and Underlying Principles. Agron. J. 101, 426–437. https://doi.org/10.2134/agronj2008.0139s.

Stöckle, C.O., Donatelli, M., Nelson, R., 2003. CropSyst, a cropping systems simulation model. Eur. J. Agron. 18, 289–307. https://doi.org/10.1016/S1161-0301(02)00109-0.

Tao, F., Palosuo, T., Rötter, R.P., Díaz-Ambrona, C.G.H., Inés Mínguez, M., Semenov, M. A., Kersebaum, K.C., Cammarano, D., Specka, X., Nendel, C., Srivastava, A.K., Ewert, F., Padovan, G., Ferrise, R., Martre, P., Rodríguez, L., Ruiz-Ramos, M., Gaiser, T., Höhn, J.G., Salo, T., Dibari, C., Schulman, A.H., 2020. Why do crop models diverge substantially in climate impact projections? A comprehensive

analysis based on eight barley crop models. Agric. For. Meteorol. 281, 107851. https://doi.org/10.1016/j.agrformet.2019.107851.

Wallach, D., Buis, S., Lecharpentier, P., Bourges, J., Clastre, P., Launay, M., Bergez, J.-E., Guerif, M., Soudais, J., Justes, E., 2011. A package of parameter estimation methods and implementation for the STICS crop-soil model. Environ. Model. Softw. 26, 386–394. https://doi.org/10.1016/j.envsoft.2010.09.004.

Wallach, D., Buis, S., Seserman, D.-M., Palosuo, T., Thornburn, P.J., Mielenz, H., Justes, E., Kersebaum, K.-C., Dumont, B., Launay, M., Seidel, S.J., 2024. A calibration protocol for soil-crop models. Environ. Model. Softaware 180, 106147. https://doi.org/10.1016/j.envsoft.2024.106147.

Weiss, M., Baret, F., 1999. Evaluation of canopy biophysical variable retrieval performances from the accumulation of large swath satellite data. Remote Sens. Environ. 70, 293–306. https://doi.org/10.1016/S0034-4257(99)00045-0.

Weiss, M., Troufleau, D., Baret, F., Chauki, H., Prévot, L., Olioso, A., Bruguier, N., Brisson, N., 2001. Coupling canopy functioning and radiative transfer models for remote sensing data assimilation. Agric. For. Meteorol. 108, 113–128. https://doi.org/10.1016/S0168-1923(01)00234-9.

Weiss, M., Jacob, F., Duveiller, G., 2020. Remote sensing for agricultural applications: a meta-review. Remote Sens. Environ. 236, 111402. https://doi.org/10.1016/j.rse.2019.111402.

Williams, J.R., Jones, C.A., Dyke, P.T., 1984. A modeling approach to determining the relationship between erosion and soil productivity. Trans. ASAE 27, 0129–0144. https://doi.org/10.13031/2013.32748.

de Wit, A.J.W., van Diepen, C.A., 2007. Crop model data assimilation with the Ensemble Kalman filter for improving regional crop yield forecasts. Agric. For. Meteorol. 146, 38–56. https://doi.org/10.1016/j.agrformet.2007.05.004.

Wösten, J.H.M., Lilly, A., Nemes, A., Le Bas, C., 1999. Development and use of a database of hydraulic properties of European soils. Geoderma 90, 169–185. https://doi.org/10.1016/S0016-7061(98)00132-3.

Xu, H., Liu, M., Li, C., Tang, Y., Xue, Q., Xiao, W., Gao, D., Peng, D., Dai, X., 2024. Optimizing agronomic management practices for enhanced radiation capture and improved radiation use efficiency in winter wheat. Plants 13, 2036. https://doi.org/10.3390/plants13152036.

Yang, Y., Ogle, S., Del Grosso, S., Mueller, N., Spencer, S., Ray, D., 2022. Regionalizing crop types to enhance global ecosystem modeling of maize production. Environ. Res. Lett. 17, 014013. https://doi.org/10.1088/1748-9326/ac3f06.

Zhang, X., Friedl, M.A., Schaaf, C.B., Strahler, A.H., Hodges, J.C.F., Gao, F., Reed, B.C., Huete, A., 2003. Monitoring vegetation phenology using MODIS. Remote Sens. Environ. 84, 471–475. https://doi.org/10.1016/S0034-4257(02)00135-9.

Zhao, Y., Chen, S., Shen, S., 2013. Assimilating remote sensing information with crop model using Ensemble Kalman Filter for improving LAI monitoring and yield estimation. Ecol. Model. 270, 30–42. https://doi.org/10.1016/j.ecolmodel.2013.08.016.

Zhou, Y., 2019. Asymmetric behavior of vegetation seasonal growth and the climatic cause: evidence from long-term NDVI dataset in Northeast China. Remote Sens. 11, 2107. https://doi.org/10.3390/rs11182107.

Zhuo, W., Huang, J., Xiao, X., Huang, H., Bajgain, R., Wu, X., Gao, X., Wang, J., Li, X., Wagle, P., 2022. Assimilating remote sensing-based VPM GPP into the WOFOST model for improving regional winter wheat yield estimation. Eur. J. Agron. 139, 126556. https://doi.org/10.1016/j.eja.2022.126556.