# Simulating phosphorus dynamics between the soil and the crop with the STICS model: Formalization and multi-site evaluation on maize in temperate area

Mounir Seghouani, Matthieu Bravin, Patrice Lecharpentier, Alain Mollier

► To cite this version:

Mounir Seghouani, Matthieu Bravin, Patrice Lecharpentier, Alain Mollier. Simulating phosphorus dynamics between the soil and the crop with the STICS model: Formalization and multi-site evaluation on maize in temperate area. European Journal of Agronomy, 2025, 164, pp.127475. 10.1016/j.eja.2024.127475. hal-04850769

HAL Id: hal-04850769

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

Submitted on 20 Dec 2024

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>

# Simulating phosphorus dynamics between the soil and the crop with the STICS model: Formalization and multi-site evaluation on maize in temperate area

Mounir Seghouani<sup>a,b,c</sup>, Matthieu N. Bravin<sup>b,c</sup>, Patrice Lecharpentier<sup>d</sup>, Alain Mollier<sup>a,*</sup>

## ARTICLE INFO

### Keywords:

Modeling

Plant

Soil-crop model

Field experiment

Fertilization

P response

## ABSTRACT

Soil-crop models are pertinent tools to study and manage phosphorus (P) in agroecosystems. However, P modeling is suffering a delay as compared to nitrogen and carbon. In this study, we extended the STICS model to simulate the P uptake and P feed-back by coupling it with a soil-plant P model. The paper aims at describing the P model and present the results showing the ability of the model to simulate contrasting P uptake and growth response pattern of maize submitted to different level of P inputs. in temperate area. The model simulates the soil P availability and the crop P demand, uptake, and partitioning. A major originality of this work is that it relies on soil solution P concentration and P sorption curves to simulate soil P availability, and critical P dilution curves to simulate crop P demand. We evaluated the model against a dataset coming from four field fertilization trials located at different site in mainland France. The trials consisted of fertilizing maize with a mineral fertilizer at three application rates (P0, P1, P2) which induced contrasted crop responses to P shortage. The model has shown great ability in predicting P uptake both dynamically and at the end of the cropping season for the entire dataset (EF >0.75). The model has satisfactory predictions of crop biomass accumulation (EF >0.5) and leaf area index. Considering each fertilization level separately, the evaluation has shown that the model had predicted the final P uptake of P1 and P2 treatments better than that of P0 treatment (EF of 0.74, 0.73 and 0.62 for P2, P1, and P0, respectively). The predictions made for the P0 treatment remained nonetheless satisfactory for both P uptake and plant growth. The good performance of the model is promising as it shows that the model is sufficiently robust to simulate maize P uptake across a range of soil P availability and P fertilization under contrasting temperate climatic conditions. Further validation on other crop species and soil and climatic conditions is discussed.

## 1. Introduction

Agroecosystems are complex systems that involve several components that may interact. While experimentation is a pertinent tool to understand components, it has difficulties to catch the complexity of such systems (Wallach et al., 2014). Crop models are useful tools which complement experimentation and allow to understand the overall performance of the system (Jones et al., 2017), while incorporating our knowledge on its components (Das et al., 2019). They allow notably to test in a short time, different scenarios that could not be tested experimentally (Hinsinger et al., 2011). They have notably been used extensively to optimize nitrogen (N) fertilization and management (Salo et al.,

2016). In contrast, there are little quantitative tools to monitor phosphorus (P) crop nutrition (Reid et al., 2011).

Modelers have been interested in modeling the P cycle since the 1970s. The models developed were mainly biophysical models focusing on specific mechanisms and their applications (Daroub et al., 2003). We can notably cite the mathematical models that have been developed to simulate P fluxes to roots (Claassen and Barber, 1974). However, such models are not adapted to handle the P cycle at the scale of a crop (Probert, 2004), and crop models that account for P remain scarce. Many reasons may explain the delay in the development of P crop modeling as compared to other elements such as N. According to Das et al., (2019), this is mainly due to the fact that most recent models (e.g. Daroub et al.,

<footnote><sup>a</sup> UMR1391 ISPA, INRAE, Bordeaux Sciences Agro, Villenave d'Ornon F-33882, France</footnote>

<footnote><sup>b</sup> CIRAD, UPR Recyclage et risque, Montpellier F-34398, France</footnote>

<footnote><sup>c</sup> Recyclage et risque, Univ Montpellier, CIRAD, Montpellier F-34398, France</footnote>

<footnote><sup>d</sup> INRAE, US1116 AgroClim, Avignon F-84914, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: alain.mollier@inrae.fr (A. Mollier).</footnote>
2003; Dzotsi et al., 2010; Probert, 2004) are based on the pioneering work of Jones et al. (1984) to integrate P in the EPIC model through a conceptual representation of soil P pools. While this approach may have good results, the difficulty in measuring, initializing, and evaluating this pool is a major limitation to the development of P modeling (Das et al., 2019). Yin et al. (2021) raise the concern that models grows in number but are not necessary enhanced from this point of view. Another break is the importance of rhizospheric processes in P crop nutrition and the difficulty in quantifying and formalizing them in soil-crop models (Hinsinger et al., 2011). In addition, most of the efforts invested in P modeling have focused on the soil compartment (Lewis and McGechan, 2002; Pferdmenges et al., 2020), while the plant compartment and its response to P shortage remains scarce.

The need for soil-crop models that address P nutrition and its effect on crop growth is particulary important in the context of agroecological transition to be able to simulate low-input systems (Probert, 2004). The current trend in soil-crop modeling is towards more mechanistic approaches (Yin et al., 2021). Semi-mechanistic models that handles N and water nutrition are numerous (Salo et al., 2016). In particular, STICS is a dynamic, generic, and semi-mechanistic soil-crop model. It handles both N and water nutrition of the crop and it is used to monitor both yield responses and environmental variables (Brisson et al., 2003). It has proven its accuracy in predicting both crop and soil variables, especially under European conditions (Coucheney et al., 2015). However such models may have limited usefulness in low-input systems as well as in tropical agroecosystems where other nutrients may limit crop growth (Delve et al., 2009). The need for models that manage both N and P stress has been reported in a recent review in order to address N-P co-limitation, which is common in agroecosystems (Seghouani et al., 2024). Integrating a P model into an existing soil-crop model is a common approach (e.g. Daroub et al., 2003; Dzotsi et al., 2010; Jones et al., 1984; Probert, 2004).

This work aims to introduce new soil and plant P modules in the STICS soil-crop model that allow the simulation of soil P availability, P demand, P uptake and the plant growth response to P. This study describes the formalisms of both soil and plant P modules and their integration in STICS, and evaluates of the performance of the model in predicting P uptake and crop growth using multi-site field measurements of maize in France.

## 2. Model description

The STICS P model is developed on the basis of the FUSSIM-P model (Mollier et al., 2008). It is structured as a set of modules (Fig. 1.). The P modules are coupled to the STICS soil-crop model, which was previously described by Beaudoin et al. (2023). The P modules run in the following order: (1) soil P availability (2) crop P demand (3) root P uptake (4) P

stress and (5) P partitioning. The P modules exchange variables among themselves and with the STICS core at a daily time step (Supplementary Information 1). Hereafter are described the P modules and their formalisms.

### 2.1. Soil P supply

Soil P supply is computed for each 1 cm layer (z) of the soil. It is represented both by the content of P in the soil solution ($Qw$, mg P kg$^{-1}$ DM) and by the amount of P ($Pr$, mg P kg$^{-1}$ DM) that can replenish the soil solution in 24 hours from the solid phase under a depletion gradient. The P content of the soil solution is calculated from the concentration ($C_p$, mg P L$^{-1}$) of orthophosphate ions ($H_2PO_4$ and $HPO_4^{2-}$ denoted as oP ions thereafter), the soil water content ($\theta$, L L$^{-1}$), and the bulk soil density ($D_a$, kg DM L$^{-1}$) as follows:

$$Q_w = \frac{\theta \times C_p}{D_a} \qquad (1)$$

$Pr$ is computed from $C_p$ by a Freundlich equation describing the solid-solution P equilibrium in soil as presented by (Schneider and Morel, 2000):

$$P_r = \nu \times C_p^w \times t^p \qquad (2)$$

Where, $\nu$, $w$, and $p$ are the Freundlich kinetic parameters that characterize the isotherm and, $t$ is the time (min). The available P in each layer ($E$, mg P kg$^{-1}$ DM) is obtained from the sum of both $Pr$ and $Qw$.

$$E = (P_r + Q_w) \qquad (3)$$

The model only accounts for the soil P mineral fraction and consequently neglects the contribution of organic P to plant nutrition. In fact, Raguet et al. (2023) showed that oP ions released by the mineralization of soil organic P are negligible compared to the oP ions released by desorption-diffusion processes at the solid/solution interface.

### 2.2. Plant P demand

Daily plant P demand is driven by biomass production. Aerial part demand is based on the dilution curve concept as described for major nutrients such as N (Lemaire et al., 2008) and P (Lemaire et al., 2019; Morel et al., 2021). The model computes the optimal P concentration ($P_{crit}$), g P kg$^{-1}$ DM) as a function of biomass (denoted as masec, t DM ha$^{-1}$ in STICS) as following:

$$[P_{crit}] = a \times masec^{-b} \qquad (4)$$

Where, a (g P kg$^{-1}$ DM) and b (dimensionless) are parameters characterizing the P dilution curve. The daily plant P demand (BasePdemand,

<figure><img src="image_3.png" /><figcaption>Fig. 1. Structure of the P module and its interactions with the STICS soil-crop model.</figcaption></figure>
kg P ha⁻¹) is derived from Eq. 4 according to the daily biomass growth increment (dltams, t DM ha⁻¹) as:

$$BasePdemand = dltams \times a(1-b) \times masec^{-b} \quad (5)$$

As the P dilution curves apply only from a minimum of biomass, the demand was considered as constant and maximal below this biomass (masecPmax, t DM ha⁻¹) and was computed as:

$$\text{if}(\text{masec} \le \text{masecPmax}) \text{ then } \text{BasePdemand} = dltams \times a \quad (6)$$

In order to simulate situations of P accumulation where the P concentration is higher than $[P_c]$ (i.e. luxury P uptake), we computed AerialPdemand (kg P ha⁻¹) from BasePdemand, taking into account the plant nutrition status and the surplus of P in the soil as compared to the BasePdemand as:

$$\text{AerialPdemand} = \text{BasePdemand} + (S_{sm} - \text{BasePdemand}) \times \text{Accumulationcoef} \quad (7)$$

Where, $S_{sm}$ is the plant's potential P uptake (kg P ha⁻¹), which is detailed below in Eq. 13. Eq. 7 is only computed if $S_{sm}$ is greater than Base P demand. Otherwise the AerialPdemand is equal to either BasePdemand according to Eqs. 5 and 6 or the amount of P required to reach $P_c$. In order to limit excessive P accumulation, a maximum P demand, Max P demand (kg P ha⁻¹), is computed according to Eq. 5, but taking into account for a maximum P concentration (Pmax) with $a_{max}$ and $b_{max}$ (Eq. 4). The daily AerialPdemand cannot be higher than MaxPdemand.

The daily P demand (Pdemand, kg P ha⁻¹) is equal to the sum of aerial and root P demand. The root P demand (RootPdemand, kg P ha⁻¹) is computed from a constant optimal root P content ($P_{root}$, g P kg⁻¹ DM), the daily root biomass production (dltamsrac, kg DM ha⁻¹) and the P nutrition index (PNI, dimensionless) as:

$$\text{RootPdemand} = P_{\text{root}} \times dltamsrac \times \text{PNI} \quad (8)$$

The PNI is computed according to Eq. 16 in order to account for variations in P allocation to roots during P shortage or over-accumulation (see below for rationale).

## 2.3. Root P uptake

The root P uptake module features the formalisms of the uptake module of FUSSIM-P (Mollier et al., 2008) which were adapted to the 1D description of soil processes in STICS. Accordingly, each soil layer is characterized by a height ($\Delta z$) equal to 1 cm and soil properties including a $C_p$ and a $P_r$ computed from Eq. 2. Each soil layer is also characterized by a root length density ($lracz$, cm cm⁻³) with a fixed mean root radius ($r_0$, cm). Phosphorus uptake is then computed for each soil layer using the microscopic nutrient uptake model developed by De Willigen and Van Noordwijk (1987) as:

$$\left(\frac{dQ}{dC_p} + \theta\right) \left(\frac{dC_p}{dt}\right) = \frac{D}{R} \frac{d}{dR} R \frac{dC_p}{dR} - ql \frac{dC_p}{dR} - S_s \quad (9)$$

Where, $dQ/dC_p$ is the soil buffering capacity (mL cm⁻³), $t$ is the time (day), $ql$ is the water flux towards the root (cm day⁻¹) computed by STICS, $R$ (cm) is the radial distance to the root center which varies between the mean root radius $r_0$ (cm) and the radius of the soil attributed to each soil $r_l$ (cm) calculated in Eq. 12 (below), $D$ is the diffusion coefficient of oP ions in the soil, and, $S_s$ is the P uptake rate (kg P ha⁻¹). The model given in Eq. 9 therefore accounts for the transport of oP ions to roots by both diffusion and mass flow.

The diffusion coefficient $D$ (m² s⁻¹) is calculated according to Barraclough and Tinker (1981) as:

$$D = D_0 f_l \quad (10)$$

Where, $D_0$ is the oP ions diffusion coefficient in free water (m² s⁻¹) and,

$f_l$ is the tortuosity factor of the soil (dimensionless). $f_l$ is computed from the broken line function of $\theta$ proposed by Barraclough and Tinker (1981) as:

$$f_l = \begin{cases} \frac{(f_l \theta + f_2)}{\theta(f_1 \theta + f_2)} & \text{for } \theta \ge \theta_l \\ \frac{\theta(f_1 \theta + f_2)}{\theta_l} & \text{for } \theta < \theta_l \end{cases} \quad (11)$$

Where, $f_1$ and $f_2$ are dimensionless parameters determined from soil properties and, $\theta_l$ is the water content at the intersection of the two lines.

The model assigns a soil to each root. This volume consists of a of soil characterized by the of the layer (1 cm) and on a cylinder of soil around each root characterized by the radius $r_l$ as:

$$r_l = \frac{1}{\sqrt{\pi L_{racz}}} \quad (12)$$

The maximum P uptake rate ($S_{sm}$, kg P ha⁻¹) is calculated as the sum of the maximum uptake rate of each layer. We assume that the root behaves as a zero sink and that the maximum uptake rate is equal to the maximum rate of transport to the root. $S_{sm}$ is thus derived from the steady-rate approximate solution for the concentration profile around the root (De Willigen and Van Noordwijk, 1994) as:

$$S_{sm} = \sum_{\Delta z=1}^{n} 2\pi \Delta z L_{racz} D \left( \frac{\rho^2 - 1}{2G(\rho, v)} \right) \bar{C}_p \quad (13)$$

Where, $\rho$ is the normalized radius computed as $r_l/r_0$ and $G(\rho, v)$ is the following geometry function:

$$G(\rho, v) = \frac{1}{2(v+1)} \left( \frac{1-\rho^2}{2} + \frac{\rho^2(\rho^{2v}-1)}{2v} + \frac{\rho^2(\rho^{2v}-1)(v+1)}{2v(\rho^{2v+2}-1)} + \frac{(1-\rho^{2v+4})(v+1)}{(2v+4)(\rho^{2v+2}-1)} \right) \quad (14)$$

Where, $v$ is the dimensionless uptake of water such as:

$$v = \frac{q_l}{4\Delta z\pi L_{racz} D} \quad (15)$$

If the maximum uptake rate $S_{sm}$ is greater than the crop P demand (PDemand), then the actual uptake $S_{sr}$ is equal to PDemand. If $S_{sm}$ cannot satisfy the crop P demand, then $S_{sr}$ is equal to $S_{sm}$. Effective uptake is then distributed across the layers, by reflecting the proportion of the root length of the layer relative to the total root length. If the roots cannot take up P at the required rate in one layer, the root system can compensate in other layers according to an iterative scheme. After accounting for root P uptake, $C_p$ is adjusted according to Eq. 2.

## 2.4. Phosphorus stress feed-back on crop growth and development

Plant P nutritional status is estimated through PNI as:

$$\text{PNI} = \frac{[P]}{[P_{\text{crit}}]} \quad (16)$$

If the PNI is less than 1, the crop is under P shortage. Conversely, a PNI greater than 1 indicates a situation of P luxury uptake. Accordingly, a stress index (Stress index, dimensionless) is computed with a scaling factor (stresscoef, dimensionless) that accounts for the crop sensitivity to P shortage as previously proposed by Beaudoin et al. (2023) in the STICS model:

$$\text{Stressindex} = \text{stresscoef} \times \min(\text{PNI}, 1) + 1 - \text{stresscoef} \quad (17)$$

The Stress index was used to impact crop development (Stress index dev), root elongation (Stress index root), and leaves production and senescence (Stress index leaf) by running Eq. 17 with stresscoef dev, stresscoef leaf, and stresscoef root, respectively.

The model simulates crop development through an accumulation of development units (Upvt, °Cd) which are computed using the sum of
effective temperature (Udevcult, °Cd), from one phenological stage to another. Under P shortage, the daily accumulation of development units is reduced, resulting in a delay in crop development, as:

$$Upvt = Udevcult \times Stressindexdev \qquad (18)$$

The reduction of daily new leaf production (Deltai, dimensionless), is simulated as:

$$Deltai = Deltai \times Stressindexleaf \qquad (19)$$

The effect of a P shortage on leaves senescence is simulated by a reduction of leaves lifespan (durvieF, °Cd) as:

$$durvieF = durevieF \times Stressindexleaf \qquad (20)$$

Plant adaptation to P shortage is further simulated by the production of longer roots as:

$$rlj = msrac \times Stressindexroots \qquad (21)$$

Where, *rlj* is the daily root length production (m d⁻¹) and, *msrac* is the root biomass (t ha⁻¹).

The roots produced are also thinner under P shortage. This effect is simulated by increasing the specific root length (longsperac, cm g⁻¹) when computing *msrac* as:

$$msrac(t) = msrac(t - 1) + \frac{rlj}{longsperac \times 100} \times Stressindexroots \qquad (22)$$

To simulate thinner roots, the model reduces the radius of newly produced roots (newrootradius, cm) as:

$$newrootradius = radius \times Stressindexroots \qquad (23)$$

Where, *radius* is a crop parameter describing the mean root radius (cm) under optimal conditions.

Finally, *r₀* which is the mean root radius of the crop, is updated by taking into account the radius of the newly produced roots (Eq. 24):

$$r_0(t) = \frac{r_0(t-1) \times rltot(t-1) + xrlj \times newrootradius}{rltot(t)} \qquad (24)$$

Where, *rltot* is the total root length of the crop per soil surface area (cm cm⁻²). and *xrlj* is the daily production of root length.

## 2.5. Plant P partitioning

Plant P partitioning is subdivided from total plant P (QPtot, kg ha⁻¹) in P into grain (QPgrain) and P in the vegetative parts (QPveg). The latter is further divided into P in aerials (QPaerial) and root parts (QProot). The initial P content of the seed is considered as an input parameter (QPinit).

QPgrain is computed from the beginning of the phenological stage “grain filling” (idrp) as:

$$QPgrain = PHI \times QPaerial \qquad (25)$$

Where, *PHI* is the P harvest index.

By analogy with the computation of the N harvest index in STICS (Beaudoin et al., 2023), the *PHI* is computed as:

$$PHI = vitirpho \times (day - idrp) + PHI_{ini} \qquad (26)$$

Where, *vitirpho* and *PHI ini* are, respectively, the rate of increase of *PHI* vs time and the *initial PHI* value. The P allocation to the grain is simulated from the day of the start of grain filling (idrp). The *PHI* cannot be greater than a maximum threshold (*PHImax*).

We account for a better of P remobilization to the grain when P nutrition is not optimal during the grain filling (eq.27):

$$PHI = PHI/min(PNI, 1) \qquad (27)$$

All P uptake and P demand parameters described in the model equations are listed in Table 1. Parameters were mostly derived from

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

Soil and crop parameters description soil and values used during the simulations.</caption><thead><tr><th>Symbols</th><th>Explanation</th><th>Value</th><th>Units</th><th>Source</th></tr></thead><tbody><tr><td>Soil<br>D<sub>0</sub></td><td>Diffusion coefficient of oP ions in free water at 25 °C</td><td>0.76</td><td>cm<sup>2</sup><br>d<sup>-1</sup></td><td>Lide and Frederikse (1996)</td></tr><tr><td>f<sub>1</sub></td><td>Tortuosity parameters (Eq. 11)</td><td>1.58 (S1 and S2)<br>0.99 (S3 and S4)</td><td>-</td><td>Barraclough and Tinker (1981)</td></tr><tr><td>f<sub>2</sub></td><td>Tortuosity parameters (Eq. 11)</td><td>-0.17</td><td>-</td><td>Barraclough and Tinker (1981)</td></tr><tr><td>θ<sub>t</sub></td><td>Tortuosity parameters (Eq. 11)</td><td>0.12</td><td>-</td><td>Barraclough and Tinker (1981)</td></tr><tr><td>Crop<br>r<sub>0</sub></td><td>Mean root radius</td><td>0.02</td><td>cm</td><td>STICS corn file base value</td></tr><tr><td>QPinit</td><td>Initial P content (seeds P content)</td><td>0.065</td><td>kg P<br>ha<sup>-1</sup></td><td>Estimated from data (Sowing density and mean grain P content)</td></tr><tr><td>a</td><td>P critical dilution curve parameter</td><td>3.72</td><td>-</td><td>Estimated from data (evaluation and independent dataset)</td></tr><tr><td>b</td><td>P critical dilution curve parameter</td><td>0.254</td><td>-</td><td>Estimated from data (evaluation and independent dataset)</td></tr><tr><td>a<sub>max</sub></td><td>P maximum accumulation curve parameter</td><td>6</td><td>-</td><td>Estimated from data (evaluation and independent dataset)</td></tr><tr><td>b<sub>max</sub></td><td>P maximum accumulation curve parameter</td><td>0.27</td><td>-</td><td>Estimated from data (evaluation and independent dataset)</td></tr><tr><td>masecPmax</td><td>Shoot biomass threshold of P dilution curve</td><td>1</td><td>-</td><td>Gagnon et al. (2020)</td></tr><tr><td>Accumulationcoef</td><td>Accumulation Coef to account for P demand over the critical curve</td><td>0.06</td><td>-</td><td>trial-error optimisation (evaluation dataset)</td></tr><tr><td>[P<sub>root</sub>]</td><td>Optimal phosphorus concentration in the roots</td><td>1.2</td><td>g P<br>kg<sup>-1</sup></td><td>Li et al. (2017)</td></tr><tr><td>vitirPho</td><td>The rate of increase PHI vs time</td><td>0.0077</td><td>-</td><td>linear regression from data + trial-error optimisation (evaluation dataset)</td></tr><tr><td>PHI<sub>ini</sub></td><td>PHI value at beginning of grain filling</td><td>0.135</td><td>-</td><td>linear regression from data + trial-error optimisation (evaluation dataset)</td></tr><tr><td>PHI<sub>max</sub></td><td>Maximum value of the PHI</td><td>0.9</td><td>-</td><td>Estimated from data (evaluation dataset)</td></tr></tbody><tfoot><tr><td colspan="5" style="text-align: right;">(continued on next page)</td></tr></tfoot></table>
<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Symbols</th><th>Explanation</th><th>Value</th><th>Units</th><th>Source</th></tr></thead><tbody><tr><td>stresscoefleaf</td><td>Leaves growth and senescence sensitivity to P stress</td><td>0.4</td><td>-</td><td>trial-error optimisation (evaluation dataset)</td></tr><tr><td>stresscoefroots</td><td>Roots sensitivity to P stress (adaptation)</td><td>0.375</td><td>-</td><td>trial-error optimisation (evaluation dataset)</td></tr><tr><td>stresscoefdev</td><td>Development sensitivity to phenology</td><td>0 (insensitive)</td><td>-</td><td>Abiotic stresses (e.g Water and N) on phenology are disabled for maize crop plant file by default.</td></tr></tbody></table>

measured observations or from the literature. Few parameters were optimized by trial and error. In this work we have considered that the plant parameters are species-specific but they could be adapted to cultivars if this is justified and possible.

### 3. Material and methods

#### 3.1. Field trials description

The model was evaluated against observed crop data selected from four long-term P field experiments (hereafter referred to as sites; Table 2) of the French network on P fertilization experiments (Boniface and Trocme, 1988). The sites reflect a diversity of pedological and climatic conditions. Decades of different P fertilization rates resulted in contrasting soil and plant P contents and plant biomass. All experiments had three levels of P annual input as superphosphate ((Ca(H₂PO₄)₂ 2 H₂O), 45 % P₂O₅): a control with no P input (P0) and two levels of P input corresponding to either once (P1) or twice (P2) the grain P export. The treatments were arranged in a randomized complete block design with four replications at all sites. Nitrogen and potassium were homogeneously applied on the three P fertilization levels, with annual inputs of 180–200 kg N ha⁻¹ and 80–100 kg K ha⁻¹. Sites one (S1), two (S2), and four (S4) were irrigated and the site three (S3) was rainfed. Soils were sampled after the P fertilization application to measure the initial Cₚ at the beginning of the cropping cycle (Table 2). The Freundlich kinetic parameters are presented at the end of Table 2. They were determined from sorption/desorption or ³²P exchange experiments (Fardeau et al., 1991; Schneider and Morel, 2000). Morel et al. (2021) showed that E₂₄h was closely related to crop P status regardless of soils. The

<caption>**Table 2**</caption>

Sites characteristics, soil properties and cropping practices of the four field experiments. Mean values and standard errors are given between parentheses (n = 4).

<table><thead><tr><th>Site</th><th>Unit</th><th>S1</th><th>S2</th><th>S3</th><th>S4</th></tr></thead><tbody><tr><td>Site location</td><td>-</td><td>Caracarés Ste-Croix (40)</td><td>Mant (40)</td><td>Thierval-Grignon (78)</td><td>Saint Félix (16)</td></tr><tr><td>Site coordinates</td><td>-</td><td>43°52' N, 0°44''W</td><td>43°35' N, 0°30''W</td><td>48°50' N, 1°56''E</td><td>45°22' N, 0°00''E</td></tr><tr><td>Starting year of experimentation</td><td>year</td><td>1972</td><td>1975</td><td>1986</td><td>1978</td></tr><tr><td>Year of simulation</td><td>year</td><td>1995–1998</td><td>1991</td><td>1998</td><td>1998</td></tr><tr><td>P fertilization</td><td>kg P ha⁻¹ y⁻¹</td><td>P0: 0<br>P1: 44<br>P2: 96</td><td>P0: 0<br>P1: 21<br>P2: 79</td><td>P0: 0<br>P1: 26<br>P2: 52</td><td>P0: 0<br>P1: 30<br>P2: 45</td></tr><tr><td>Cultivar</td><td>-</td><td>Volga from 1995 to 1997<br>Cecilia in 1998</td><td>Cecilia</td><td>Nobilis</td><td>Cecilia</td></tr><tr><td>Mean annual T°</td><td>°C</td><td>13.2</td><td>13.2</td><td>10.9</td><td>12.6</td></tr><tr><td>Mean annual rainfall</td><td>mm</td><td>590</td><td>590</td><td>643</td><td>786</td></tr><tr><td>FAO soil Classification</td><td>-</td><td>Luvic Arenosol</td><td>Luvic Arenosol</td><td>Luvisol</td><td>Calcosol</td></tr><tr><td>Thickness of the ploughed layer</td><td>cm</td><td>25</td><td>28</td><td>30</td><td>20</td></tr><tr><td>Max depth of soil</td><td>cm</td><td>> 100</td><td>> 100</td><td>> 120</td><td>~30</td></tr><tr><td>Bulk density of the ploughed layer</td><td>g/cm3</td><td>1.5</td><td>1.5</td><td>1.3</td><td>1.3</td></tr><tr><td>Pebbles > 2 mm</td><td>% of total soil</td><td>negligible</td><td>negligible</td><td>24</td><td>5</td></tr><tr><td>Clay</td><td>g kg⁻¹</td><td>60 (10)</td><td>120 (7)</td><td>285 (15)</td><td>389 (53)</td></tr><tr><td>Silt</td><td>g kg⁻¹</td><td>135 (34)</td><td>520 (15)</td><td>609 (21)</td><td>290 (45)</td></tr><tr><td>Sand</td><td>g kg⁻¹</td><td>806 (71)</td><td>361 (17)</td><td>107 (7)</td><td>43 (6)</td></tr><tr><td>pHwater</td><td>-</td><td>5.9 (0.2)</td><td>7.3 (0.4)</td><td>8.1 (0.1)</td><td>8.1 (0.1)</td></tr><tr><td>CaCO₃</td><td>g kg⁻¹</td><td>< 2</td><td>< 2</td><td>18 (7)</td><td>268 (21)</td></tr><tr><td>C:N ratio</td><td>-</td><td>13.5</td><td>9.8</td><td>11.5</td><td>8.8</td></tr><tr><td>Olsen P</td><td>mg P kg⁻¹</td><td>P0: 22<br>P1: 63<br>P2: 82</td><td>P0: 5.3<br>P1: 22<br>P2: 68</td><td>P0: 5.8<br>P1: 34<br>P2: 55</td><td>P0: 7.4<br>P1: 20<br>P2: 30</td></tr><tr><td>CEC</td><td>cmol + kg⁻¹</td><td>3.6 (0.7)</td><td>5.9 (0.8)</td><td>18.2 (1.1)</td><td>23.3 (0.3)</td></tr><tr><td>Exch K</td><td>cmol + kg⁻¹</td><td>0.19 (0.01)</td><td>0.3 (0.1)</td><td>0.7 (0.1)</td><td>0.7 (0.2)</td></tr><tr><td>Alox</td><td>mmol kg⁻¹</td><td>22 (4)</td><td>44.1 (0.4)</td><td>56 (13)</td><td>73 (3)</td></tr><tr><td>Feox</td><td>mmol kg⁻¹</td><td>20 (5)</td><td>21.1 (0.2)</td><td>34 (5)</td><td>40 (2)</td></tr><tr><td>Freundlich kinetics parameters</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ν</td><td>-</td><td>1.76</td><td>5.70</td><td>20.7</td><td>27.9</td></tr><tr><td>ω</td><td>-</td><td>0.68</td><td>0.70</td><td>0.61</td><td>0.53</td></tr><tr><td>p</td><td>-</td><td>0.38</td><td>0.25</td><td>0.18</td><td>0.23</td></tr><tr><td>Cₚ of the ploughed layer</td><td>mg P L⁻¹</td><td>P0: 0.47 (0.13)<br>P1: 1.97 (0.37)<br>P2: 3.36 (0.44)</td><td>P0: 0.22 (0.06)<br>P1: 0.61 (0.19)<br>P2: 1.48 (0.66)</td><td>P0: 0.07 (0.01)<br>P1: 0.48 (0.08)<br>P2: 0.90 (0.20)</td><td>P0: 0.03 (0.02)<br>P1: 0.11 (0.05)<br>P2: 0.21 (0.08)</td></tr><tr><td>Cₚ of the subsoil</td><td>mg P L⁻¹</td><td>26–40 cm<br>P0: 0.17 (0.03)<br>P1: 0.47 (0.07)<br>P2: 1.44 (0.15)<br>41–60 cm<br>P0: 0.22 (0.03)<br>P1: 0.09 (0.04)<br>P2: 0.04 (0.01)</td><td>29–60 cm<br>P0: 0.08 (0.01)<br>P1: 0.11 (0.04)<br>P2: 0.47 (0.14)</td><td>31–50 cm<br>P0: 0.03 (0.002)<br>P1: 0.07 (0.005)<br>P2: 0.11 (0.02)<br>51–125 cm<br>P0: 0.03 (0.003)<br>P1: 0.04 (0.003)<br>P2: 0.06 (0.02)</td><td>20–30 cm<br>P0: 0.03 (0.02)<br>P1: 0.11 (0.05)<br>P2: 0.21 (0.08)</td></tr></tbody></table>
Freundlich kinetic parameters were considered identical for all the treatments and replicates of the same site. We used the *Cp* measured in each replicate (block) to simulate the variability of the initial *Cp* between the replicates of each treatment. Accordingly, each simulation unit (USM) is specific to the site, treatment, year, and replicate (block). We parameterized Freundlich kinetic isotherms and initial *Cp* of the subsoil for each layer. up to 60 cm for S1, S2, up to 125 cm for S3 and up to the maximum soil depth of 30 cm for S4. For the S1 and S2, which have high *Cp* in the ploughed soil layer, we considered only the 60 cm because most of the roots are found in these layers and the deeper layers do not have much impact on crop P nutrition. We considered the full depth of the S3 soil because it is a rainfed site with low *Cp* values in the soil. Thus, the contribution of deeper layers to water uptake and P nutrition is expected to be more important. The S4 was only characterized to 30 cm depth, as this is the maximum depth of the soil. The Freundlich kinetic isotherms parameters were considered homogeneous in a same site, and when the information on deeper layer sorption was not available, we considered it homogeneous over all the soil layers. Only the initial *Cp* values (both of the ploughing soil layer and of the subsoils) differed between blocks. *Cp* values were measured from soil samples collected after fertilizer application. Consequently, the effect of fertilization is not simulated explicitly but is accounted for via the initial differences of *Cp* values.

We evaluated the model using data from four cropping seasons in S1 and one cropping season for S2, S3, and S4. Plants were periodically sampled during the selected cropping seasons to monitor crop growth and P accumulation. Leaf area index (LAI) was also measured in three out of the four cropping seasons in S1 as well as in the S3 and S4.

## 3.2. Model parameterization and simulations

The P modules were coupled with version 10 of STICS (http://www6.paca.inra.fr/stics/). Soil inputs were obtained from soil analyses (Table 1). If not provided in the datasets, hydrological soil parameters (field capacity and wilting point) were obtained using the pedo transfer function proposed by Dobarco et al. (2019). Meteorological data were provided daily from the nearest meteorological station at each site. Mean temperature and cumulated rainfall are presented in Table 1. Crop management information including sowing and N fertilization, were recorded. Automatic irrigation was run in S1, S2, and S4 to simulate optimal water supply and soil was initialized at field capacity. The emergence date was forced according to the observed date. The effects of nitrogen stress on crop growth were disabled because no nitrogen stress was observed during the experiment. Simulated maize cultivars were already calibrated in STICS v10. However, some parameters governing grain filling required to be parameterized as described in Supplementary Materials 2.

## 3.3. Model evaluation

We evaluated the model both graphically (i.e. visually) and quantitatively with statistics by using the CropPlotR package (Vezy et al., 2023). We performed the quantitative evaluation on both the dynamics and the maximum values of the variables (final biomass, grain yield, maximum P accumulation, and grain P content) using several statistical indicators (Moriasi et al., 2007). The root mean square error (RMSE) was used to quantify the mean absolute prediction error of the model. The relative root mean square error (rRMSE, i.e. an adimensional RMSE) was calculated by scaling residuals to actual values and was used to compare a variable across different scales (Guillaume et al., 2011). According to Jamieson et al. (1991), an rRMSE greater than 30 % or less than 30, 20, and 10 % was considered poor, fair, good, and excellent, respectively. The coefficient of determination (R²) was calculated to determine the proportion of the variance in the observed data explained by the model. The Nash-Sutcliffe model efficiency (EF), which ranges from -∞ to 1, was calculated to quantify the relative size of the residual variance compared

to the measured data variance. The EF indicates the degree of agreement between predicted and observed values. An EF of 1 would be achieved by a perfect of model and an EF of 0 is equivalent to the prediction of the mean of observations. According to Moriasi et al. (2007), > 0.5, > 0.65, and > 0.75 allowed to evidence a model performance rated as to be satisfactory, good, or very good, respectively.

# 4. Results

## 4.1. Model performance in predicting P accumulation and dry biomass in aerial parts and grain

We assessed the model performance in predicting the accumulation of P in the aerial parts and grain, as well as the aerial biomass and yield (Fig. 2). Overall, the model had an excellent performance in predicting the final P accumulation in the aerial parts (Fig. 2.a) with a model efficiency of 0.82 and a good prediction error (rRMSE = 16 %), which is equivalent to 6.3 kg P ha⁻¹. Simulated vs. observed data were homogeneously distributed around the 1:1 line, showing that there is no systematic over- or under-prediction of P accumulation in aerial parts. We also found a good correlation between observed and simulated data as the model was able to explain 82 % of the observed variability. Furthermore, the model was able to simulate the sixfold range of P accumulation, as the simulated P accumulation varied from 9.0 to 61 kg ha⁻¹. Regarding the prediction of the aerial biomass (Fig. 2.b), the model performance was overall satisfactory with an EF of 0.52, a good prediction error of 2.2 t ha⁻¹ (10 %) and a satisfactory correlation between observations and simulations (R² = 0.62). The model predicted grain P accumulation well, with a satisfactory EF of 0.6 and a good prediction error of 5.7 kg P ha⁻¹ (rRMSE = 20 %) (Fig. 2.c). The model also explained 70 % of the variation in the observed grain P accumulation. The evaluation shows that the model simulated the grain yield well with a low RMSE of 1.2 t ha⁻¹ (rRMSE = 11 %) (Fig. 2.d). However, the efficiency of the model in predicting the final grain yield is borderline (EF = 0.49) and the model explained only 56 % of the grain yield variation. This is mainly explained by the small grain yield, which only varied by a factor of two (from 7.5 to 16 t ha⁻¹).

The performance of the model was strongly influenced by the P treatments (Table 3). We found that the model performed better in simulating P accumulation for the P2 and P1 with good EF (0.78 and 0.73 respectively) compared to the P0, whose efficiency was satisfactory (EF = 0.62). Furthermore, the prediction error (< 20 %) and the correlation were good for all P treatments. The model performance in simulating crop biomass was also contrasted between P treatments. The model predicted the final biomass of P2 and P1 well with an EF of 0.65 and 0.57, respectively, and an excellent prediction error (< 10 %) and a good correlation (R² = 0.79 and 0.76). In contrast, the model prediction was less adequate for the P0 treatment, with an EF of -0.1. However, the model simulation remained good with a prediction error of 2.9 t ha⁻¹ (14 %). The model had a satisfactory to good prediction error in simulating grain P accumulation which was equivalent between P treatments (18–20 %). However, the model EF and R² were lower for the P1 treatment (EF = 0.12, R² = 0.57) because the observations had less variation compared to P2 (EF = 0.48, R² = 0.76) and P0 (EF = 0.65, R² = 0.67). We found that the model also performed better in predicting the grain yield for the P2 treatment with a satisfactory efficiency of 0.62 and a very good prediction error of 8.8 %. Compared to P2, the model performance was lower for P1 and P0 with a low EF (<0.5). However, the prediction error remains nonetheless good for P1 and P0 with a good rRMSE (< 20 %).

When comparing between sites (Table 3), we found that the prediction error for simulating aerial P accumulation was better in S1 and S2 with respectively a good rRMSE of 14 % and 10 %, respectively, as compared to S3 and S4, for which the rRMSE was only satisfactory at 25 % and 26 %, respectively. The model appeared to predict dry biomass accumulation equally well across the 4 sites with good rRMSE
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Simulated versus observed values at harvest for (a) phosphorus (P) accumulation and (b) dry biomass in the whole plant aerial parts and (c) P accumulation and (d) yield in grains for the three P fertilization treatments (P0, P1, and P2) across the four sites (S1, S2, S3, and S4). Four indicators of the prediction quality (RMSE, rRMSE, EF, and R², see the Section 3.3 for rationale) are given as insert.</figcaption></figure>

<caption>**Table 3**</caption>

Evaluation metrics assessing the quality of simulations (see the Section 3.3 for rationale) as compared to the observations in each P treatment and site for accumulated P in the aerial parts, dry biomass accumulation, accumulated P in the grain and grain yield at harvest.

<table><thead><tr><th rowspan="2"></th><th colspan="4">P in the aerial parts kg P ha<sup>-1</sup></th><th colspan="4">Dry aerial biomass t DM ha<sup>-1</sup></th><th colspan="4">P in the grain kg P ha<sup>-1</sup></th><th colspan="4">Grain yield t DM ha<sup>-1</sup></th></tr><tr><th>RMSE</th><th>rRMSE</th><th>R<sup>2</sup></th><th>EF</th><th>RMSE</th><th>rRMSE</th><th>R<sup>2</sup></th><th>EF</th><th>RMSE</th><th>rRMSE</th><th>R<sup>2</sup></th><th>EF</th><th>RMSE</th><th>rRMSE</th><th>R<sup>2</sup></th><th>EF</th></tr></thead><tbody><tr><td>P treatments</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>P0</td><td>5.2</td><td>18.9</td><td>0.81</td><td>0.62</td><td>2.9</td><td>14.4</td><td>0.49</td><td>-0.1</td><td>4.4</td><td>19.5</td><td>0.67</td><td>0.65</td><td>1.4</td><td>14.3</td><td>0.74</td><td>-0.02</td></tr><tr><td>P1</td><td>5.1</td><td>12.6</td><td>0.73</td><td>0.73</td><td>1.9</td><td>8.3</td><td>0.76</td><td>0.57</td><td>6.5</td><td>20.4</td><td>0.57</td><td>0.12</td><td>1.2</td><td>10.1</td><td>0.59</td><td>0.29</td></tr><tr><td>P2</td><td>8</td><td>16.5</td><td>0.78</td><td>0.74</td><td>1.8</td><td>7.7</td><td>0.79</td><td>0.65</td><td>6.2</td><td>18.3</td><td>0.76</td><td>0.48</td><td>1.1</td><td>8.8</td><td>0.6</td><td>0.62</td></tr><tr><td>Sites</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>S1</td><td>6.5</td><td>14.4</td><td>0.74</td><td>0.75</td><td>2.4</td><td>10.4</td><td>0.34</td><td>-0.02</td><td>5.5</td><td>17.2</td><td>0.64</td><td>-0.05</td><td>1.1</td><td>9.5</td><td>0.58</td><td>0.02</td></tr><tr><td>S2</td><td>4.8</td><td>10.1</td><td>0.75</td><td>0.69</td><td>2.4</td><td>10.3</td><td>0.76</td><td>-0.28</td><td>7.3</td><td>18.0</td><td>0.64</td><td>-0.1</td><td>1.3</td><td>9.4</td><td>0.14</td><td>0.01</td></tr><tr><td>S3</td><td>5.5</td><td>25.4</td><td>0.86</td><td>-0.89</td><td>2.3</td><td>13.6</td><td>0.12</td><td>-2.16</td><td>4.8</td><td>26.4</td><td>0.46</td><td>-0.73</td><td>1.0</td><td>11.5</td><td>0</td><td>-0.1</td></tr><tr><td>S4</td><td>6.5</td><td>26.1</td><td>0.68</td><td>0.3</td><td>3.0</td><td>14.7</td><td>0.42</td><td>-0.21</td><td>6.0</td><td>31.0</td><td>0.28</td><td>0.12</td><td>1.9</td><td>16.3</td><td>0.37</td><td>-0.1</td></tr></tbody></table>

($<$ 20 %). Similar to aerial P accumulation, the model performed better in predicting grain P accumulation in S1 and S2 with good rRMSE of 17 % and 18 %, respectively, compared to S3 and S4 which had only a satisfactory rRMSE of 26 and 31 %, respectively. For grain yield, the rRMSE was lower in S1 and S2 with an excellent prediction error ($<$ 10 %). The prediction error was still good in S3 and S4 ($<$ 20 %). Overall, the model seems to be less effective in explaining the variation of all the variables in S3 as compared to the other sites. The R² values were lower in S3. This is particularly true for the dry aerial biomass, for which the model captured only 12 % and none of the variation, respectively. The model was not efficient in predicting all the variables at each site.

When comparing each individual combination of site and P treatment, we found that in most situations the model had a good prediction error for all variables of interest with only few exceptions (see Supplementary Materials 3).

## 4.2. Temporal dynamics of the growth and P accumulation in aerial parts

The model underestimated the early accumulation of P in the aerial parts as compared to the observed data (Fig. 3 and the Supplementary Materials 4 and 5 for the dynamics of P accumulation and aerial biomass, respectively). This is particularly evident in S1 and S2 (Fig. 3a and b). The model fitted P accumulation from the middle to the end of the growing season depending on the site. Despite the tendency of the model to partially underestimate the observed data, it was able to adequately reproduce the response of P accumulation in aerial parts with increasing P fertilization. The model also simulated a higher variation in P accumulation in S2 and S4 as compared to S1 and S3 for which almost no variation was simulated.

The model predicted the temporal increase in the dry biomass of maize aerial parts throughout the cropping cycle, except for the P0 treatment in S3 and S4, for which the dry biomass was overestimated
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Temporal dynamics all along the cropping cycle of the simulated (lines) and observed (filled circles) phosphorus (P) accumulation in the whole aerial parts of maize for the sites a) S1, b) S2, c) S3, and d) S4 and the three P fertilization treatments (P0, P1, and P2). While the site S1 was monitored from 1995 to 1998, only the year 1997 is depicted in a) (see SM 4 for the observed and simulated data for the three other years). Shaded areas and error bars stand for standard deviation for simulations and observations, respectively. Green and blue dots overlap on some measurement dates.</figcaption></figure>

(Fig. 4). The reduction in dry biomass observed in all sites for the P0 treatment was adequately reproduced by the model for S1 and S2. However, the model tended to neglect or at least underestimate the intensity of this reduction in S3 and S4. The model simulated very little variations around dry biomass and notably less than the observed dry biomass at the four sites.

### 4.3. Temporal dynamics of the variables used to simulate the stress induced by P

The model was able to simulate the feedbacks of P deficiency on crop growth (see the Supplementary Materials 6 and 7 for the dynamics of LAI and root surface, respectively). In particular, it reproduced the overall reduction in LAI growth during the cropping season. The model was also able to simulate a clear feedback on root growth, with the production of a higher root surface at lower P levels. The simulated feedbacks were observed at all sites and were particularly pronounced for the P0 treatment.

As a consequence of the simulated feedbacks for both the LAI and the relative root surface area, the model simulated for the four sites and from the second part of the growing cycle, an increasing leaf-to-root ratio with the increasing P treatment, i.e. from P0 to P2 (Fig. 5 and the Supplementary Materials 8). Thus, the model simulated a higher root development at the expense of leaf development under P stress.

## 5. Discussion

### 5.1. Model formalisms are consistent with agronomic knowledge

We have proposed a novel approach to predict P uptake in soil-crop models. We followed general guidelines of parsimony as many parts of the proposed model were close or analogous to those developed for N nutrition in STICS, in particular the use of dilution curves and harvest index (Beaudoin et al., 2023). This was done according to the of Probert (2004), who argues that models with the same level of complexity between their components perform better. We also considered the specificity of P in both soil and plant.

Since the pioneering work of Jones et al. (1984), which proposed a module for predicting soil P dynamics in the EPIC model, modeling of the P cycle has been less successful compared to N and carbon modeling (Das et al., 2019). This is because models have difficulty in linking soil P measurements to the P pool that is readily available to the crop (Das et al., 2019). Here, we propose an approach based on the notion of P sorption-desorption curve described by a Freundlich isotherm, which allows a more mechanistic and robust representation of the soil-P physicochemical interactions (Morel et al., 2014, 2021; Ziadi et al., 2013). Another important difference of our soil P availability formalism, compared to other models is its relative simplicity as we consider only inorganic P. This choice is supported by the findings of Raguet et al. (2023), who showed that over a single cropping cycle the mineralization of organic P contributes only marginally to crop P nutrition in temperate agroecosystems. Finally, the model requires only a few parameters to describe the soil P availability (oP ions), the concentration of ions in the
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Temporal dynamics all along the cropping cycle of the simulated (lines) and observed (filled circles) dry biomass in the whole aerial parts of maize for the sites a) S1, b) S2, c) S3, and d) S4 and the three P fertilization treatments (P0, P1, and P2). While the site S1 was monitored from 1995 to 1998, only the year 1997 is depicted in a) (see SM 5 for the observed and simulated data for the three other years). Shaded areas and error bars stand for standard deviation for simulations and observations, respectively.</figcaption></figure>

soil solution and the Freundlich isotherm parameters. This is an
advantage as overparameterization has been identified as a major lim-
itation in modeling P nutrition (Das et al., 2019).

The proposed model simulates both the diffusive and convective P
fluxes to the root surface. Thus, crop P uptake is highly dependent on soil
moisture. This effect well was captured in the model and was clearly
highlighted in S3, where soil moisture was highly variable because it
was rainfed. Therefore, such a model would be particularly relevant to
study the interaction between crop P nutrition and soil water avail-
ability. No P uptake limitation by root length unit has been included
because crop P uptake is thought to be limited by soil P availability
rather than root uptake capacity (Rengel, 1993).

The major difference from the older nutrient uptake model that we
have incorporated into the formalism is the consideration of crop P
demand in computing the effective P uptake. Although determining the
optimal and maximum P concentration for each developmental stage is a
common approach in models (Daroub et al., 2003; Dzotsi et al., 2010),
we preferred to use a more dynamic estimation of the optimal and
maximum P concentration in the plant by linking it to biomass growth
according to the dilution curves (Lemaire et al., 2019, 2008). This
approach has already demonstrated its utility and reliability (Cadot
et al., 2018; Fontana et al., 2021; Morel et al., 2021). We obtained good
results using a single dilution curve for three different maize varieties.
This is in line with the results of Gómez et al. (2019), who found that the
of the dilution curves are specific parameters and do not vary much
between cultivars of the same group for potato crop. The effect of spe-
cies, cultivar, root age, and plant nutritional status in other major nu-
trients (especially N) on P uptake was not considered in our model.
These effects are rarely considered in nutrient uptake models because

their consideration would complicate the models too much (Mollier
et al., 2008), while nutrient uptake models perform reasonably well
without considering them (Barber, 1995).

The use of phosphorus nutrition index (PNI) on a daily basis allows the diagnosis of P shortage periods. This would support the evaluation of management practices for P nutrition. We used a simple linear relationship between time and the phosphorus harvest index (PHI) to represent the P allocation to the grain. This PHI has been reported to be effective in quantifying P allocation to the grain (El Mazlouzi et al., 2020). The PNI effect on the PHI was added to account for the greater remobilization of P when crops are suffering from P deficiency. This is due to the fact that the remobilization of P is independent from that of biomass (Rodriguez and Goudriaan, 1995).

It is well known that P deficiency mainly leads to a reduced above-
ground biomass (Assuero et al., 2004; Plenet et al., 2000; Shi et al.,
2020). We have reproduced this effect through feedbacks on interme-
diate variables in accordance with the current knowledge of plant
physiology. The main simulated effect of P on crop growth is a reduction
in the leaf area expansion in line with the findings from field experi-
ments (Assuero et al., 2004; Plenet et al., 2000; Rodriguez et al., 2000)
and also an increase in leaf senescence through a reduction of leaf life-
span, which has also been reported in crops suffering from P shortage as
an adaptation mechanism to ensure P remobilization during the grain
filling stage (El Mazlouzi et al., 2020). Because the STICS model simu-
lates biomass production through a radiation use efficiency (RUE)
approach, we did not simulate a reduction in photosynthesis due to P
deficiency, unlike other mechanistic approaches. While a severe P
deficiency reduces photosynthetic activity (Naeem et al., 2010; Shi
et al., 2020), field experiments show that radiation use efficiency (RUE)
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Temporal dynamics all along the cropping cycle of the simulated leaf to root surface ratio (leaf/root) of maize for the sites a) S1, b) S2, c) S3, and d) S4 and the three P fertilization treatments (P0, P1, and P2). While the site S1 was monitored from 1995 to 1998, only the year 1997 is depicted in a) (see SM 8 for the observed and simulated data for the three other years). Shaded areas stand for standard deviation.</figcaption></figure>

is not affected by soil P levels (Fletcher et al., 2008; Plenet et al., 2000). P deficiency is also responsible for multiple adaptations of the root system, including changes in root depth, density, architecture and radius (Lopez et al., 2023; Niu et al., 2013; Vance et al., 2003). As a result of the P shortage feedback, a reduction in the shoot:root ratio is commonly reported under field conditions (Malhotra et al., 2018). The model was able to simulate the main effect of a P shortage, including the reduction in leaf expansion and the increase in the shoot:root ratio. The reduction in leaf expansion ultimately leads to a reduction in the root surface area due to reduced photosynthesis. This effect has been reported in the literature (Mollier and Pellerin, 1999) and it was captured by the model without a direct formalization for this effect. Another important emergent effect of the model is the reduction of the total P demand of the deficient crop as a result of the decrease in biomass due to feedbacks. This response is due to the fact that plant P demand is driven by crop biomass production. It is also noteworthy that the model was able to simulate clear differences in the main (i.e. biomass and yield) and intermediate (i.e. LAI and root surface area) variables between but also within P treatments, while only the initial soil oP ions availability differed between simulation units.

## 5.2. The *P* module predicts well the final biomass and *P* content of maize under contrasting soil *P* availability

The evaluation of the developed model shows a good performance in simulating the final P accumulation in the aerial parts and its allocation to the grain as well as the final biomass and grain yield. The simulations were generally of higher prediction quality at higher levels of P fertilization. The model showed contrasting performance between sites. Two

situations deserve to be highlighted: (i) in the rainfed site (S3), it appears that the model overestimated the intensity of water stress on crop growth. This would explain the lack of response to P levels and the fact that the model underestimated the final biomass even at the highest P levels. (ii) In S4, which has a high degree of variability in soil conditions. This variability included key parameters such as the proportion of pebbles in the soils. We could not account for this variability because we did not have access to data on the distribution by plot, and we changed values of initial oP ions concentration in the soil solution between simulations. As a result, the model was not able to reproduce all of the observed variation compared to other sites.

The STICS model has been extensively evaluated under temperate conditions. An evaluation of the model on 15 crops across France shows that the model has a good prediction of final N uptake (rRMSE = 33 %, EF = 0.67) (Coucheney et al., 2015). Older work on maize crops reported an even better performance (Brisson et al., 2002). Our work suggests that the prediction of plant uptake could be as good for P as for N in STICS. However, further validation is needed to assess its robustness and genericity. When compared to other dynamic semi-mechanistic soil-crop models, we found that our model had good performance in predicting final aerial accumulation of P. The DSSAT P module proposed by Daroub et al. (2003) shows low performance in predicting P accumulation of different crops grown on both calcareous and acidic soils (18–107 % rRMSE). The reported rRMSE for maize P uptake was 41 %. A more recent version proposed by Dzotsi et al. (2010) was evaluated by Amouzou et al. (2018) for maize and sorghum. They reported good simulations of P uptake. Similar to our study, they reported a lower performance of the model at lower P levels under rainfed conditions.

Our model is based on the modeling work of Mollier et al. (2008),
who used the FUSSIM-P to simulate P uptake and biomass accumulation of early stages of maize subjected to three levels of P fertilization. We used a similar approach to simulate P availability and uptake. The model evaluations of Mollier et al. (2008) on P uptake in early stages of maize show similar results to our study with a good overall prediction of all P levels (EF > 0.65). Similar to our study, they reported a better simulation of high and medium level of P fertilization (EF > 0.9) as compared to low level of P fertilization which was underestimated at 336° Cd (EF = 0.66). Our work thus confirms the fact that the simulation of soil P and plant uptake as proposed by Mollier et al. (2008) is robust and performs well under temperate and moderate P shortage conditions.

Few works have evaluated soil-crop models for their ability to predict the final P accumulation in the grain. Das et al. (2022) reported low variation in the simulation of P allocation to the grain and they found unsatisfactory results (EF = 0.36). They suggested that the model overestimated P uptake and P allocation from moderately available P pools at high levels of P fertilization. They highlighted the need for better measurements of soil P availability, and its balance with crop P demand as well as a possible omission of feedbacks due to P limitation. In contrast, the simulation of N allocation to the grain has been more intensively evaluated. In comparison, our model shows good performance in predicting P allocation to the grain, which that could be equivalent to the STICS simulation of the N allocation to the grain.

Our model performed satisfactorily in predicting final grain yield. However, the relative error was still small. We can also argue that the relatively low EF is explained by the fact that a mean of the observations remains a good indicator of grain yield across treatments rather than by the fact that the prediction error is high. This is mainly due to the fact that there was little variation in observed grain yield between different P fertilization levels. These results are in line with literature, as yield reduction due to P shortage is reported to be less than the reduction in the total aerial biomass (Plenet et al., 2000; Fontana et al., 2021). Although the model predicted yield satisfactorily, its sensitivity to the initial soil P availability should be improved. Biomass was well predicted with slight overestimation at lower levels of P nutrition.

When compared with other modeling works such as APSIM or DSSAT, we found that our performance was similar. For example, Das et al. (2022) reported similar Nash-Sutcliff efficiency (EF = 0.46) in predicting grain yield of several crops in a 35-years field trial of contrasting N and P treatments with the APSIM model. They also reported an overestimation of yield under greater P shortage. Previous work with the APSIM model shows good simulations of biomass. However, they found satisfactory results when considering only higher and lower P levels (EF = 0.52). The lower efficiency found in the Das et al. (2022) study as compared to other APSIM work may be due to N-P co-limitations that make the response more complex and unaddressed factors such as pests and diseases. Work done with the DSSAT model on maize, peanuts and sorghum shows that the model performed well in predicting biomass and yield under contrasting P inputs. Accounting for P stress improved the performance of the DSSAT model (Dzotsi et al., 2010; Naab et al., 2015; Halder et al., 2017; Amouzou et al., 2018). Similar to our results, the FUSSIM-P model performed better at high and medium levels of P fertilization (EF) compared to low level of P fertilization (EF). However, our model overestimated the biomass of P0 while the FUSSIM-P model underestimated it (Mollier et al., 2008). They explained this underestimation mainly by an overestimation of P stress due to the lack of mechanisms that may be more important for P nutrition under P limitation. While we also neglected these mechanisms, we overestimated the final biomass of P0. These may be due to an underestimation of the feedback on biomass growth as compared to the more direct C:P equilibrium represented in FUSSIM-P model.

Overall, the model had good performance for predicting all variables of interest. It had a performance that was at least equal to most of the reported evaluations of soil-crop model dealing with P accumulation and P feedback to crop growth. We also found that the module was well implemented in the STICS model, as it had a performance in predicting

final P uptake that was at least equal to the performance of STICS model in predicting N uptake. It is also noteworthy that our work complements the work of Mollier et al. (2008) with the FUSSIM-P model by extending it to the whole cropping season. Furthermore, we highlight the fact that we obtained good results with a minimal calibration of the model, as most of the parameters were determined by measurements, literature, or by using a trial-error optimization for empirical parameters as compared to recent modeling work that relies mostly on mathematical optimization (Kherif et al., 2022; Ravelojaona et al., 2023). We believe that the model performance could be improved by using mathematical optimization of the parameters, as suggested in the literature (Wallach et al., 2014).

### 5.3. Limitations of the model and future paths of improvement

Despite the overall good performance of the model we developed, the simulation results show two contradictory issues. The first issue is related to the underestimation of the LAI temporal dynamics. The second issue is related to the final prediction of aerial biomass and grain yield. While the model had good prediction error over all P treatments, the performance decreased with soil P availability. Consequently the model tends to overestimate the final aerial biomass and grain yield of P despite a clear underestimation of LAI. With the current formalization, it is noteworthy that solving the first issue would exacerbate the second. In fact, parameterizing the model to increase plant P uptake at the beginning of the crop cycle would lead to a decrease in the feedback due to P shortage and, ultimately, to a further overestimation of the final plant aerial biomass in low P availability.

One of the main way to improve the model is to improve the root representation in the crop model. Despite the fact that we have formalized a feedback on root growth, the implemented formalisms do not allow to reflect all root adaptation strategies to P shortage, such as the changes in spatial configurations or the contrasting effects on primary and secondary roots (Niu et al., 2013; Vance et al., 2003; Jia et al., 2018).

A simple representation of the root system is sufficient to simulate N and water uptake (Coucheney et al., 2015). However, for a better simulation of the P uptake during the cycle, a more detailed representation of the root system may be required. Further improvement of these aspects of crop adaptation to P stress is limited by the 1D representation of the root-soil system. This issue has been highlighted in the literature, as it interferes not only with good P uptake but also with the simulation of certain agronomic practices such as banded P fertilization (Naab et al., 2015). However, there is still room for improvement. For example, the FUSSIM-P model which simulates a carbon allocation equilibrium between shoot and root showed a better simulation of P accumulation at the beginning of the crop cycle as compared to STICS model, which computes roots based on shoot growth, but without a proper simulation of carbon allocation with a source-sink relationship (Beaudoin et al., 2023). This is due to the fact that when the crop is limited by P, the growth of the crop leaves simulated by FUSSIM-P model is reduced and so is their demand for assimilates. Consequently, a larger fraction of assimilates is allocated to the roots (Mollier et al., 2008). This behavior is consistent with commonly reported crop adaptation mechanisms that allow relatively more carbon allocation to roots when subjected to P shortage (Vance et al., 2003). Such improvements in STICS model would allow simulation of a faster root establishment, and thus correct, at least in part, the underestimation of P at the beginning of the crop cycle. Furthermore, it would help to solve the two aforementioned issues, as it would also reduce the final aerial biomass by allocating more carbon to roots. Therefore, a better simulation of the roots would allow for a better representation of the C:P stoichiometry of the crop. This would be considered as a more intense feedback on the final aerial biomass, which would allow to reduce the intensity of the P shortage feedbacks on LAI, thus improving its simulation.

Other mechanisms involved in P nutrition that are not considered in
the model are the rhizospheric mechanisms such as proton and phosphatase release by plant roots (George et al., 2011; Hinsinger, 2001; Vance et al., 2003). Some of these mechanisms may even be involved in C:P equilibrium because they involve carbon investment. This is the case of mycorrhizal association, which increases soil exploration and contributes to P solubilization (Balemi and Negisho, 2012), and organic anion exudation, which increases P mobility (Hinsinger, 2001; Malhotra et al., 2018; Vance et al., 2003). This later involves a high carbon especially for citrate and malate, which can account for up to 20 % of carbon assimilates (Pueyo et al., 2021; Vance et al., 2003). All these considerations led to the good predictions of final P uptake in maize that we obtained under all contrasting P treatments to moderate this statement. However, we validated the STICS-P model against experimental data obtained under a relatively high level of soil P availability and consequently a moderate level of plant P shortage. Model performance therefore needs to be evaluated under more extreme P shortage conditions, such as those found in the acid soils of Sub-Saharan Africa (van der Velde et al., 2014).

In the current study, no effect of N on P nutrition was considered. The model formalisms follow a strict law of the minimum when considering N-P co-limitation. Field evidence shows that N nutrition affects P nutrition (Dotaniya and Meena, 2015), such as changes in the P dilution curve (Bélanger et al., 2015; Bélanger and Ziadi, 2008; Ziadi et al., 2007). In a recent review, we showed that in most cases (i.e., 84 % of the fertilization field experiments), the crop responds to N-P co-limitation by following the multiple limitation hypothesis (MLH) rather than the law of the minimum (LM) which is however the law considered by most crop models. Future work should evaluate the ability of the model to simulate situations of N-P co-limitations using real data sets, with formalisms that consider either the usual LM or the MLH.

## 6. Conclusion

The aim of this work was to extend the STICS soil crop model to simulate P uptake and its feedback on plant growth. We make a first evaluation of the model for maize under temperate conditions. The evaluations showed good performance of the model in predicting final P accumulation and partitioning in the grain as well as final biomass and yield over a range of soil P availability. The fact that such results were obtained with a formalization that neglects many mechanisms, such as organic P mineralization and rhizospheric processes, supports the hypothesis that these mechanisms do not play a major role in plant P nutrition under the conditions. We have also shown that a good simulation of P uptake is highly dependent on a good simulation of soil water status. In order to ensure a good simulation of P accumulation dynamics, special attention must be paid to the simulation of root development and its response to soil P availability especially at the beginning of the crop cycle. The model needs to be evaluated for other crop species and under other experimental conditions to assess both its robustness and its genericity. These experimental conditions should include soils with higher sorption capacity and lower soil P availability, greater water shortage and N-P co-limitation.

## CRediT authorship contribution statement

**Mounir Seghouani:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Matthieu N. Bravin:** Writing – review & editing, Methodology, Conceptualization.
**Patrice Lecharpentier:** Validation, Software.
**Alain Mollier:** Writing – review & editing, Software, Methodology, 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. The research done by the Recyclage et Risque group, to which Mounir Seghouani and Matthieu N. Bravin belong, aims at suggesting methods and models to recycle agricultural, agro-industrial, and urban organic residues in agriculture by minimizing agro-environmental issues and managing recycling at the territory scale.

No conflict of interest exits in the submission of this manuscript, and manuscript is approved by all authors for publication. I would like to declare on behalf of my co-authors that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part. All the authors listed have approved the manuscript that is enclosed.

## Acknowledgements

The Ph.D. of M. Seghouani was funded by the French National Institute for Agriculture, Food and Environment (INRAE, Agro-EcoSystem Division) and the French Agricultural Research Center for International Development (CIRAD, Persyst Division). We also would like to thank Mathias Christina (UPR Aïda, CIRAD), Gatien Falconnier (UPR Aïda, CIRAD), David Houben (AGHYLE, UniLaSalle), Gaetan Louarn (URP3F, INRAE), Christophe Nguyen (UMR ISPA, INRAE), and Bastien Castagneyrol (UMR Biogeco) for their interest in this work and their valuable advice.

## Appendix A. Supporting information

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

## Data availability

Data will be made available on request.

## References

* Amouzou, K.A., Naab, J.B., Lamers, J.P.A., Becker, M., 2018. CERES-Maize and CERES-Sorghum for modeling growth, nitrogen and phosphorus uptake, and soil moisture dynamics in the dry savanna of West Africa. Field Crops Res 217, 134–149. https://doi.org/10.1016/j.fcr.2017.12.017.

* Assuero, S.G., Mollier, A., Pellerin, S., 2004. The decrease in growth of phosphorus-deficient maize leaves is related to a lower cell production. Plant Cell Environ. 27, 887–895. https://doi.org/10.1111/j.1365-3040.2004.01194.x.

* Balemi, T., Negisho, K., 2012. Management of soil phosphorus and plant adaptation mechanisms to phosphorus stress for sustainable crop production: a review. J. Soil Sci. Plant Nutr. 12, 547–561. https://doi.org/10.4067/S0718-95162012005000015.

* Barber, S.A., 1995. Soil nutrient bioavailability: a mechanistic approach. John Wiley & Sons.

* Barraclough, P.B., Tinker, P.B., 1981. The determination of ionic diffusion coefficients in field soils. I. Diffusion coefficients in sieved soils in relation to water content and bulk density. J. Soil Sci. 32, 225–236. https://doi.org/10.1111/j.1365-2389.1981.tb01702.x.

* Beaudoin, N., Lecharpentier, P., Ripoche-Wachter, D., Strullu, L., Mary, B., Léonard, J., Launay, M., Justes, E. (Eds.), 2023. STICS soil-crop model: Conceptual framework, equations and uses. Quae Editions. https://doi.org/10.35690/978-2-7592-3679-4.

* Bélanger, G., Ziadi, N., 2008. Phosphorus and nitrogen relationships during spring growth of an aging timothy sward. Agron. J. 100, 1757–1762. https://doi.org/10.2134/agronj2008.0132.

* Bélanger, G., Ziadi, N., Pageau, D., Grant, C., Lafond, J., Nyiraneza, J., 2015. Shoot Growth, Phosphorus–Nitrogen Relationships, and Yield of Canola in Response to Mineral Phosphorus Fertilization. Agron. J. 107, 1458–1464. https://doi.org/10.2134/agronj15.0050.

* Boniface, R., Trocme, S., 1988. Enseignements fournis par des essais de longue durée sur la fumure phosphatée et potassique. II-Essai sur la fumure phosphatée. In: Gachon, L. (Ed.), Phosphore et potassium dans les relations sol-plante: conséquences sur la fertilisation. INRA, pp. 279–402.

* Brisson, N., Gary, C., Justes, E., Roche, R., Mary, B., Ripoche, D., Zimmer, D., Sierra, J., Bertuzzi, P., Burger, P., 2003. An overview of the crop model STICS. Eur. J. Agron. 18, 309–332. https://doi.org/10.1016/S1161-0301(02)00110-7.

* Brisson, N., Ruget, F., Gate, P., Lorgeou, J., Nicoullaud, B., Tayot, X., Plenet, D., Jeuffroy, M.-H., Bouthier, A., Ripoche, D., 2002. STICS: a generic model for simulating crops and their water and nitrogen balances. II. Model validation for wheat and maize. Agronomie 22, 69–92. https://doi.org/10.1051/agro:2001005.

* Cadot, S., Bélanger, G., Ziadi, N., Morel, C., Sinaj, S., 2018. Critical plant and soil phosphorus for wheat, maize, and rapeseed after 44 years of P fertilization. Nutr. Cycl. Agroecosyst. 112, 417–433. https://doi.org/10.1007/s10705-018-9956-0.
Claassen, N., Barber, S.A., 1974. A Method for Characterizing the Relation between Nutrient Concentration and Flux into Roots of Intact Plants. Plant Physiol. 54, 564-568. https://doi.org/10.1104/pp.54.4.564.

Coucheney, E., Buis, S., Launay, M., Constantin, J., Mary, B., de Corta, I.G., Ruget, F., Andrianarisoa, K.S., Bas, C.L., Justes, E., 2015. Accuracy, robustness and behavior of the STICS soilecile model for plant, water and nitrogen outputs: evaluation over a wide range of agro-environmental conditions in France, 2015 Environ. Model. Softw. 64, 177-190. https://doi.org/10.1016/j.envsoft.2014.11.024.

Daroub, S.H., Gerakis, A., Ritchie, J.T., Friesen, D.K., Ryan, J., 2003. Development of a soil-plant phosphorus simulation model for calcareous and weathered tropical soils. Agric. Syst. 76, 1157-1181. https://doi.org/10.1016/S0308-521X(02)00082-3.

Das, B., Huth, N., Probert, M., Condon, L., Schmidt, S., 2019. Soil phosphorus modeling for modern agriculture requires balance of science and practicality: a perspective. J. Environ. Qual. 48, 1281-1294. https://doi.org/10.2134/jeq2019.05.0201.

Das, B.T., Schmidt, S., Biggs, J., Lester, D.W., Bourne, N., Huth, N.I., 2022. Simulating long-term phosphorus, nitrogen, and carbon dynamics to advance nutrient assessment in dryland cropping. Field Crops Res 285, 108590. https://doi.org/10.1016/j.fcr.2022.108590.

De Willigen, P., Van Noordwijk, M., 1987. Roots, plant production and nutrient use efficiency. Wageningen Agricultural University, Wageningen, The Netherlands, p. 281.

De Willigen, P., Van Noordwijk, M., 1994. Mass flow and diffusion of nutrients to a root with constant or zero-sink uptake. II. Zero-sink uptake. Soil Sci. 157, 171-175.

Delve, R.J., Probert, M.E., Cobo, J.G., Ricaurte, J., Rivera, M., Barrios, E., Rao, I.M., 2009. Simulating phosphorus responses in annual crops using APSIM: model evaluation on contrasting soil types. Nutr. Cycl. Agroecosyst. 84, 293-306. https://doi.org/10.1007/s10705-008-9243-6.

Dobarco, M.R., Cousin, I., Le Bas, C., Martin, M.P., 2019. Pedotransfer functions for predicting available water capacity in French soils, their applicability domain and associated uncertainty. Geoderma 336, 81-95. https://doi.org/10.1016/j.geoderma.2018.08.022.

Dotaniya, M.L., Meena, V.D., 2015. Rhizosphere effect on nutrient availability in soil and its uptake by plants: a review. Proc. Natl. Acad. Sci., India, Sect. B Biol. Sci. 85, 1-12. https://doi.org/10.1007/s40011-013-0297-0.

Dzotsi, K.A., Jones, J.W., Adiku, S.G.K., Naab, J.B., Singh, U., Porter, C.H., Gijsman, A.J., 2010. Modeling soil and plant phosphorus within DSSAT. Ecol. Modell. 221, 2839-2849. https://doi.org/10.1016/j.ecolmodel.2010.08.023.

El Mazlouzi, M., Morel, C., Robert, T., Yan, B.-F., Mollier, A., 2020. Phosphorus uptake and partitioning in two durum wheat cultivars with contrasting biomass allocation as affected by different P supply during grain filling. Plant Soil 449, 179-192. https://doi.org/10.1007/s11104-020-04444-0.

Fardeau, J.C., Morel, C., Boniface, R., 1991. Phosphate ion transfer from soil to soil solution: kinetic parameters. Agronomie 11, 787-797.

Fletcher, A.L., Moot, D.J., Stone, P.J., 2008. Radiation use efficiency and leaf photosynthesis of sweet corn in response to phosphorus in a cool temperate environment. Eur. J. Agron. 29, 88-93. https://doi.org/10.1016/j.eja.2008.04.002.

Fontana, M., Bélanger, G., Hirte, J., Ziadi, N., Elfouki, S., Bragazza, L., Liebisch, F., Sinaj, S., 2021. Critical plant phosphorus for winter wheat assessed from long-term field experiments. Eur. J. Agron. 126, 126263. https://doi.org/10.1016/j.eja.2021.126263.

Gagnon, B., Ziadi, N., Bélanger, G., Parent, G., 2020. Validation and use of critical phosphorus concentration in maize. Eur. J. Agron. 120, 126147. https://doi.org/10.1016/j.eja.2020.126147.

George, T.S., Fransson, A.-M., Hammond, J.P., White, P.J., 2011. Phosphorus nutrition: rhizosphere processes, plant response and adaptations. In: Bünemann, E., Oberson, A., Frossard, E. (Eds.), Phosphorus in Action, Soil Biology. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 245-271. https://doi.org/10.1007/978-3-642-15271-9_10.

Gómez, M.I., Magnitskiy, S., Rodríguez, L.E., 2019. Critical dilution curves for nitrogen, phosphorus, and potassium in potato group andigenum. Agron. J. 111, 419-427. https://doi.org/10.2134/agronj2018.05.0357.

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.

Halder, D., Panda, R.K., Srivastava, R.K., Kheroar, S., 2017. Evaluation of the CROPGRO-Peanut model in simulating appropriate sowing date and phosphorus fertilizer application rate for peanut in a subtropical region of eastern India. Crop J. 5, 317-325. https://doi.org/10.1016/j.cj.2017.02.005.

Hinsinger, P., 2001. Bioavailability of soil inorganic P in the rhizosphere as affected by root-induced chemical changes: a review. Plant Soil 237, 173-195. https://doi.org/10.1023/A:1013351617532.

Hinsinger, P., Brauman, A., Devau, N., Gérard, F., Jourdan, C., Laclau, J.-P., Le Cadre, E., Jaillard, B., Plassard, C., 2011. Acquisition of phosphorus and other poorly mobile nutrients by roots. Where do plant nutrition models fail? Plant Soil 348, 29-61. https://doi.org/10.1007/s11104-011-0903-y.

Jamieson, P.D., Porter, J.R., Wilson, D.R., 1991. A test of the computer simulation model ARCWHEAT1 on wheat crops grown in New Zealand. Field Crops Res., 27 337-350. https://doi.org/10.1016/0378-4290(91)90040-3.

Jia, X., Liu, P., Lynch, J.P., 2018. Greater lateral root branching density in maize improves phosphorus acquisition from low phosphorus soil. J. Exp. Bot. 69, 4961-4970.

Jones, C.A., Cole, C.V., Sharpley, A.N., Williams, J.R., 1984. A simplified soil and plant phosphorus model: I. Documentation. Soil Sci. Soc. Am. J. 48, 800-805. https://doi.org/10.2136/ssaj1984.03615995004800040020x.

Jones, J.W., Antle, J.M., Basso, B., Boote, K.J., Conant, R.T., Foster, I., Godfray, H.C.J., Herrero, M., Howitt, R.E., Janssen, S., Keating, B.A., Munoz-Carpena, R., Porter, C.

H., Rosenzweig, C., Wheeler, T.R., 2017. Brief history of agricultural systems modeling. Agric. Syst. 155, 240-254. https://doi.org/10.1016/j.agsy.2016.05.014.

Kherif, O., Seghouani, M., Justes, E., Plaza-Bonilla, D., Bouhenache, A., Zemmouri, B., Dokukin, P., Latati, M., 2022. The first calibration and evaluation of the STICS soil-crop model on chickpea-based intercropping system under Mediterranean conditions. Eur. J. Agron. 133, 126449. https://doi.org/10.1016/j.eja.2021.126449.

Lemaire, G., Jeuffroy, M.-H., Gastal, F., 2008. Diagnosis tool for plant and crop N status in vegetative stage. Eur. J. Agron. 28, 614-624. https://doi.org/10.1016/j.eja.2008.01.005.

Lemaire, G., Sinclair, T., Sadras, V., Bélanger, G., 2019. Allometric approach to crop nutrition and implications for crop diagnosis and phenotyping. A review. Agron. Sustain. Dev. 39, 27. https://doi.org/10.1007/s13593-019-0570-6.

Lewis, D.R., McGechan, M.B., 2002. A Review of Field Scale Phosphorus Dynamics Models. Biosyst. Eng. 82, 359-380. https://doi.org/10.1006/bioe.2002.0102.

Li, H., Mollier, A., Ziadi, N., Shi, Y., Parent, L.-É., Morel, C., 2017. The long-term effects of tillage practice and phosphorus fertilization on the distribution and morphology of corn root. Plant Soil 412, 97-114. https://doi.org/10.1007/s11104-016-2925-y.

Lide, D.R., Frederikse, H.P.R., 1996. CRC handbook of chemistry and physics: a ready-reference book of chemical and physical data, 77th ed. CRC Press, Boca Raton, p. 2491.

Lopez, G., Ahmadi, S.H., Amelung, W., Athmann, M., Ewert, F., Gaiser, T., Gocke, M.I., Kautz, T., Postma, J., Rachmilevitch, S., Schaaf, G., Schnepf, A., Stoschus, A., Watt, M., Yu, P., Seidel, S.J., 2023. Nutrient deficiency effects on root architecture and root-to-shoot ratio in arable crops. Front. Plant Sci. 13. https://doi.org/10.3389/fpls.2022.1067498.

Malhotra, H., Vandana, Sharma, S., Pandey, R., 2018. Phosphorus nutrition: plant growth in response to deficiency and excess. In: Hasanuzzaman, M., Fujita, M., Oku, H., Nahar, K., Hawrylak-Nowak, B. (Eds.), Plant Nutrients and Abiotic Stress Tolerance. Springer Singapore, Singapore, pp. 171-190. https://doi.org/10.1007/978-981-10-9044-8_7.

Mollier, A., De Willigen, P., Heinen, M., Morel, C., Schneider, A., Pellerin, S., 2008. A two-dimensional simulation model of phosphorus uptake including crop growth and P-response. Ecol. Modell. 210, 453-464. https://doi.org/10.1016/j.ecolmodel.2007.08.008.

Mollier, A., Pellerin, S., 1999. Maize root system growth and development as influenced by phosphorus deficiency. J. Exp. Bot. 50, 487-497. https://doi.org/10.1093/jxb/50.333.487.

Morel, C., Plenet, D., Mollier, A., 2021. Calibration of maize phosphorus status by plant-available soil P assessed by common and process-based approaches. Is it soil-specific or not? Eur. J. Agron. 22, 126174. https://doi.org/10.1016/j.eja.2020.126174.

Morel, C., Ziadi, N., Messiga, A., Bélanger, G., Denoroy, P., Jeangros, B., Jouany, S., Fardeau, J.-C., Mollier, A., Parent, L.-E., Proix, N., Rabeharisoa, L., Jouany, S., 2014. Modeling of phosphorus dynamics in contrasting agroecosystems using long-term field experiments. Can. J. Soil. Sci. 94, 377-387. https://doi.org/10.4141/cjss2013-024.

Moriasi, D.N., Arnold, J.G., Van Liew, M.W., Bingner, R.L., Harmel, R.D., Veith, T.L., 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50, 885-900. https://doi.org/10.13031/2013.231531.

Naab, J.B., Boote, K.J., Jones, J.W., Porter, C.H., 2015. Adapting and evaluating the CROPGRO-peanut model for response to phosphorus on a sandy-loam soil under semi-arid tropical conditions. Field Crops Res 176, 71-86. https://doi.org/10.1016/j.fcr.2015.02.016.

Naeem, M., Khan, M.M.A., Moinuddin, Idrees, M., Aftab, T., 2010. Phosphorus ameliorates crop productivity, photosynthetic efficiency, nitrogen-fixation, activities of the enzymes and content of nutraceuticals of Lablac purpureus L. Sci. Hortic. 126, 205-214. https://doi.org/10.1016/j.scienta.2010.07.009.

Niu, Y.F., Chai, R.S., Jin, G.L., Wang, H., Tang, C.X., Zhang, Y.S., 2013. Responses of root architecture development to low phosphorus availability: a review. Ann. Bot. 112, 391-408. https://doi.org/10.1093/aob/mcs285.

Pferdmenges, P., Breuer, L., Julich, S., Kraft, P., 2020. Review of soil phosphorus routines in ecosystem models. Environ. Model. Softw. 126, 104639. https://doi.org/10.1016/j.envsoft.2020.104639.

Plenet, D., Etchebest, S., Mollier, A., Pellerin, S., 2000. Growth analysis of maize field crops under phosphorus deficiency. Plant Soil 223, 119-132. https://doi.org/10.1023/A:1004877111238.

Probert, M.E., 2004. A capability in APSIM to model phosphorus responses in crops: Model. Nutr. Manag. Trop. Crop. Syst., ACIAR Proc. Ser. 92-100.

Pueyo, J.J., Quinones, M.A., Coba de la Pena, T., Fedorova, E.E., Lucas, M.M., 2021. Nitrogen and phosphorus interplay in lupin root nodules and cluster roots. Front. Plant Sci. 12, 644218. https://doi.org/10.3389/fpls.2021.644218.

Raguet, P., Cade-Menun, B., Mollier, A., Abdi, D., Ziadi, N., Karam, A., Morel, C., 2023. Mineralization and speciation of organic phosphorus in a sandy soil continuously cropped and phosphorus-fertilized for 28 years. Soil Biol. Biochem. 178, 108938. https://doi.org/10.1016/j.soilbio.2022.108938.

Ravelojaona, N., Jégo, G., Ziadi, N., Mollier, A., Lafond, J., Karam, A., Morel, C., 2023. STICS soil-crop model performance for predicting biomass and nitrogen status of spring barley cropped for 31 years in a gleysoitic soil from northeastern Quebec (Canada). Agronomy 13, 2540. https://doi.org/10.3390/agronomy1302540.

Reid, J.B., Searle, B.P., Sinton, S., 2011. Does my potato crop need fertiliser? A modelling approach. Agron. N. Z. 41. https://www.agronomysociety.org.nz/files/2011_8_Doe_s_potato_crop need_fertiliser.pdf.

Rengel, Z., 1993. Mechanistic simulation models of nutrient uptake: a review. Plant Soil 152, 161-173.
Rodriguez, D., Andrade, F.H., Goudriaan, J., 2000. Does assimilate supply limit leaf expansion in wheat grown in the field under low phosphorus availability? Field Crops Res. 12.

Rodriguez, D., Goudriaan, J., 1995. Effects of phosphorus and drought stresses on dry matter and phosphorus allocation in wheat. J. Plant Nutr. 18, 2501-2517. https://doi.org/10.1080/01904169509365080.

Salo, T.J., Palosuo, T., Kersebaum, K.C., Nendel, C., Angulo, C., Ewert, F., Bindi, M., Calanca, P., Klein, T., Moriondo, M., Ferrise, R., Olesen, J.E., Patil, R.H., Ruget, F., Takáč, J., Hlavinka, P., Trnka, M., Rötter, R.P., 2016. Comparing the performance of 11 crop simulation models in predicting yield response to nitrogen fertilization. J. Agric. Sci. 154, 1218-1240. https://doi.org/10.1017/S0021859615001124.

Schneider, A., Morel, C., 2000. Relationship between the isotopically exchangeable and resin-extractable phosphate of deficient to heavily fertilized soil: isotopically exchangeable and resin-extractable phosphate. Eur. J. Soil Sci. 51, 709-715. https://doi.org/10.1046/j.1365-2389.2000.00351.x.

Seghouani, M., Bravin, M.N., Mollier, A., 2024. Crop response to nitrogen-phosphorus colimitation: theory, experimental evidences, mechanisms, and models. A review. Agron. Sustain. Dev. 44, 3. https://doi.org/10.1007/s13593-023-00939-z.

Shi, Q., Pang, J., Yong, J.W.H., Bai, C., Pereira, C.G., Song, Q., Wu, D., Dong, Q., Cheng, X., Wang, F., Zheng, J., Liu, Y., Lambers, H., 2020. Phosphorus-fertilisation has differential effects on leaf growth and photosynthetic capacity of Arachis hypogaea L. Plant Soil 447, 99-116. https://doi.org/10.1007/s11104-019-04041-w.

Vance, C.P., Uhde-Stone, C., Allan, D.L., 2003. Phosphorus acquisition and use: critical adaptations by plants for securing a nonrenewable resource. N. Phytol. 157, 423-447. https://doi.org/10.1046/j.1469-8137.2003.00695.x.

Vezy, R., Munz, S., Gaudio, N., Launay, M., Lecharpentier, P., Ripoche, D., Justes, E., 2023. Modeling soil-plant functioning of intercrops using comprehensive and generic formalisms implemented in the STICS model. Agron. Sustain. Dev. 43, 61. https://doi.org/10.1007/s13593-023-00917-5.

van der Velde, M., Folberth, C., Balkovič, J., Ciais, P., Fritz, S., Janssens, I.A., Obersteiner, M., See, L., Skalský, R., Xiong, W., Peñuelas, J., 2014. African crop yield reductions due to increasingly unbalanced Nitrogen and Phosphorus consumption. Glob. Change Biol. 20, 1278-1288. https://doi.org/10.1111/gcb.12481.

Wallach, D., Makowski, D., Jones, J.W., Brun, F., 2014. Working with dynamic crop models: methods, tools and examples for agriculture and environment, Second edition. ed. Elsevier/Academic Press, Amsterdam.

Yin, X., Struik, P.C., Goudriaan, J., 2021. On the needs for combining physiological principles and mathematics to improve crop models. Field Crops Res 271, 108254. https://doi.org/10.1016/j.fcr.2021.108254.

Ziadi, N., Bélanger, G., Cambouris, A.N., Tremblay, N., Nolin, M.C., Claessens, A., 2007. Relationship between P and N concentrations in corn. Agron. J. 99, 833-841. https://doi.org/10.2134/agronj2006.0199.

Ziadi, N., Whalen, J.K., Messiga, A.J., Morel, C., 2013. Assessment and modeling of soil available phosphorus in sustainable cropping systems. Adv. Agron. 122, 85-126. https://doi.org/10.1016/B978-0-12-417187-9.00002-4.