# Towards a global model for soil inorganic phosphorus dynamics: dependence of exchange kinetics and soil bioavailability on soil physicochemical properties

Ying-ping Wang, Yuanyuan Huang, Laurent Augusto, Daniel S Goll, Julian

Helfenstein, Enqing Hou

► To cite this version:

Ying-ping Wang, Yuanyuan Huang, Laurent Augusto, Daniel S Goll, Julian Helfenstein, et al. To-wards a global model for soil inorganic phosphorus dynamics: dependence of exchange kinetics and soil bioavailability on soil physicochemical properties. Global Biogeochemical Cycles, 2022, 36 (3), pp.e2021GB007061. 10.1029/2021GB007061. hal-03607493

HAL Id: hal-03607493

https://hal.science/hal-03607493

Submitted on 14 Mar 2022

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>
<figure><img src="image_2.png" /></figure>

# Towards a global model for soil inorganic phosphorus dynamics: dependence of exchange kinetics and soil bioavailability on soil physicochemical properties

Ying-Ping Wang¹*, Yuanyuan Huang¹, Laurent Augusto², Daniel S Goll³, Julian Helfenstein⁴, Enqing Hou⁵

1 CSIRO Oceans and Atmosphere, Aspendale, Victoria 3195 Australia

2 INRAE, Bordeaux Sciences Agro, UMR 1391 ISPA, 33882 Villenave d'Ornon, France

3 Université Paris Saclay, CEA-CNRS-UVSQ, LSCE/IPSL, Gif sur Yvette, France

4 Agroscope, Zurich 8046, Switzerland

5 Key Laboratory of Vegetation Restoration and Management of Degraded Ecosystems, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, 510650, China

To whom all correspondence should be sent. Yingping.wang@csiro.au

Submission to: Global Biogeochemical Cycles

Accepted Article
Key points

We developed and calibrated a model of soil inorganic P dynamics using the measured soil P fractions and fraction of added remaining over time of 147 soils;

We derived empirical relationships between model parameters and some soil chemical
properties;

Soil P bioavailability depends on soil P fractions, solution P concentration, desorption rate constants and the time scale.

Accepted Article
**Abstract** The representation of phosphorus cycling in global land models remains quite simplistic, particularly on soil inorganic phosphorus. For example, sorption and desorption remain unresolved and their dependence on soil physical and chemical properties is ignored. Empirical parameter values are usually based on expert knowledge or data from few sites with debatable global representativeness in most global land models. To overcome these issues, we compiled from data of inorganic soil phosphorus (P) fractions and calculated the fraction of added P remaining in soil solution over time of 147 soil samples to optimize three parameters in a model of soil inorganic P dynamics. The calibrated model performed well (r²>0.7 for 122 soil samples). Model parameters vary by several orders of magnitude, and correlate with soil P fractions of different inorganic pools, soil organic carbon and oxalate extractable metal oxide concentrations among the soil samples. The modelled bioavailability of soil P depends on, not only, the desorption rates of labile and sorbed pool, inorganic phosphorus fractions, the slope of P sorbed against solution P concentration, but also on the ability of biological uptake to deplete solution P concentration and the time scale. The model together with the empirical relationships of model parameters on soil properties can be used to quantify bioavailability of soil inorganic P on various timescale especially when coupled within global land models.

Key words phosphorus fractions, global modelling, Hedley fractionation, isotopic exchange kinetics, sorption, desorption, available phosphorus

**Plain Language Summary**

Phosphorus is a major nutrient limiting the productivity of many terrestrial ecosystems. About 20% to 60% of soil phosphorus is in inorganic form, and most inorganic soil P is sorbed or fixed on soil particles, leaving a small fraction (<1%) in soil solution available for direct uptake by plants. Sorption and desorption control inorganic P in solution and vary significantly with soil properties. However, sorption and desorption are not explicitly represented in most global land models. This study developed and calibrated a model of inorganic P dynamics using the observations from 147 soils worldwide. We found that the parameters in the model can vary by several orders of magnitude, and that a significant proportion of those variations can be explained by soil chemical properties, particularly soil P fractions, oxalate extractable metal oxide and soil organic carbon concentrations. The model

Accepted Article
and empirical relationships between model parameters and soil properties as developed in this study can be used to improve the representation of P cycle in land models.

Accepted Article
# 1. Introduction

Nitrogen (N) and phosphorus (P) are soil nutrients limiting primary production and other processes of terrestrial ecosystems (e.g. Lebauer and Treseder 2008; Fernandez-Martinez et al. 2014, Hou et al. 2020; and Jiang et al. 2021), and the limitation will likely intensify under future conditions (Penuelas et al. 2013; Zhang et al. 2014). This has far reaching implications for land carbon uptake (Goll et al 2012; Zhang et al., 2014), allowable CO₂ emissions (Wang et al. 2015), and carbon-climate feedback (Zhang et al. 2014; Arora et al. 2020). As a result, nutrient cycles have been implemented into an increasing number of global land models (see Achat et al. 2016a). While several land models (>10) have incorporated the N cycle (see Manzoni and Porporato 2009), the number of global land models with a P cycle remains quite small (<6) and the representation of P cycle is generally quite simplistic in both the parameter values and representation of some key processes (see Achat et al. 2016a; Fleischer et al. 2019).

The dynamics of soil P, in particular its availability for biological uptake (bioavailability) on time scales relevant for carbon cycling and land management, remains a major uncertainty in assessing phosphorus effects on the land carbon balance (Sun et al 2017). The representation of soil P in global models at present is largely based on the conceptual model developed by McGill and Cole (1981) who emphasized the importance of soil enzymes for soil P bioavailability, and the partitioning of different soil P pools based on Hedley fractions (Cross and Schlesinger 1995; Yang et al. 2013). Because of scarcity of global observational data, model equations for soil P dynamics were largely derived from a conceptual understanding with the model parameters chosen arbitrarily, or assumed to be globally constant (Wang and Goll 2021).

Different from inorganic soil N that only accounts for a few percent of total soil N for most mineral soils, inorganic soil P constitutes a significant fraction (20% to 60%) of total soil P with only a very small fraction (<1%) in soil solution (Morel et al. 2000; Mengel and Kirkby 2001). Although validity of chemical fractionation for soil phosphorus was debated (see Condron and Newman 2011), the Hedley fractionation technique (Hedley et al. 1982; Tiessen and Moir 1993), has become the widely used method to chemically separate soil inorganic P into a few distinct pools, named *soluble, labile, sorbed, occluded and primary mineral* (see Cross and Schlesinger 1995; Hou et al. 2018). Most of inorganic P is unavailable for direct uptake by organisms. A dynamic equilibrium is formed between the solution P, which is readily available for biotic uptake, and less available forms like labile P that is desorbed from

Accepted Article
the surface of soil particles and non-labile P that is strongly fixed by hydrous oxides and silicate minerals (Barrow 1999; Frossard et al. 2000). When solution P is depleted as a result of biological uptake, it is replenished by the desorption from other inorganic P pools and by mineralization of organic P. This dynamic equilibrium between inorganic P in soil solution and fixed P strongly depends on soil physical and chemical properties (Achat et al. 2016b).

Because of the large spatial variations of soil chemical and physical properties, dynamics of inorganic P exchanges in soil varies greatly among different soil types, terrain topography, vegetation, climate and so on (Achat et al. 2016b; Helfenstein et al. 2018a). These variations are largely ignored in most global land models. For example, in the CASA-CNP global biogeochemical model (Wang et al. 2010), the values of two parameters in the Langmuir equation for sorption were obtained by calibrating the soil P pool fractions at steady state against the mean fractions of different P pools for each of the twelve soil orders as estimated by Cross and Schlesinger (1995) based on a small number of observations (n=88), particularly in tropical soils (n=7). As a consequence, uncertainties of the estimated model parameters are very large and the influence of soil physical and chemical properties remains largely unresolved (see Achat et al. 2016b). Little advances have been made in the last decade (e.g. Kvakic et al. 2018), and several global models used similar parameter set (e.g. Goll et al 2012, Yang et al 2014, Zhu et al 2019).

Most global land models with P cycle use the Langmuir equation to simulate sorption of inorganic P in soil (see Wang et al. 2010; Goll et al. 2012; Yang et al. 2014). It has been argued that the Freundlich rather than the Langmuir equation should be used for simulating P sorption/desorption in soils (Barrow 2008), but only few global models applied it (e.g. Goll et al 2017). The Langmuir based type of approach for P sorption with fixed sorption capacity has been shown to lead to simulation of unrealistically strong P limitation under increasing CO2 (see Goll et al. 2012, their Figure 5). Theoretically, a Langmuir equation is valid when the surface of the sorbing particles is uniform and their properties are not affected by the sorption reaction. Both requirements are not satisfied for soil P sorption. Better suited is the Freundlich equation as it can be derived for heterogeneous soils for which the log of the binding constant decreases with the amount of sorption (Barrow 2008). Furthermore, most global land models assume that an equilibrium is reached within a model integration timestep (usually between 30 minutes to 1 day) between solution P and labile P (see Wang et al. 2010 and Yang et al. 2014 for example).

Accepted Article
Exchange of solution and labile or sorbed phosphorus can be measured in the laboratory using carrier-free isotopically labelled P (H₃³PO₄), or isotopic exchange kinetics (IEK) experiments (see Frossard et al. 1993, Fardeau et al. 1996). Recently, Helfenstein et al. (2018b) combined isotopic method and X-ray absorption spectroscopy to quantify the exchange kinetics of different inorganic soil P pools along a climate gradient in Hawaii. They found that the amount of P extracted using Hedley fractionation matched well the isotopically exchangeable amount of P at different timescale. This was further confirmed by Helfenstein et al. (2020) using 57 soil samples with both measured fractions of inorganic P in different pools and IEK data. Helfenstein et al. (2020) found that the mean residence time of inorganic P in each pool was quite variable. However, estimates of residence times for different pools alone are insufficient for constraining soil inorganic P dynamics. Furthermore, because of the hysteresis in the phosphate sorption and desorption isotherms (see Guedes et al. 2016), and soil's ability to replenish the depleted inorganic P in soil solution would be overestimated based on sorption isotherms (Okajima et al. 1983). Therefore, we need to represent the dynamics of sorption and desorption explicitly, and to take into account the effects of soil physical and chemical properties on sorption and desorption.

Several proxies for soil P bioavailability based on soil P fractionation were used for quantifying P limitation on ecosystem functioning on timescale relevant for carbon cycling: varying from the amount of labile P (see Sun et al. 2017) to the sum of labile and sorbed P (Ringeval et al. 2017). However, these proxies omit the dynamical nature of soil P where large fluxes connect a small stock of bioavailable P with less available forms. We need to quantify better the bioavailability of soil inorganic P based on the dynamics of inorganic P transformation, and by taking into account the dependency of P transformations on soil chemical and physical properties. To do so, we expanded the data compiled by Helfenstein et al. (2020) by including the 102 soil samples from the French forests from Achat et al. (2016b). The objectives of this study are to (1) develop and calibrate a model of inorganic P dynamics that is suitable for use in global land models; (2) quantify how model parameters vary with soil physical and chemical properties, and climate; (3) define soil bioavailable P (a<sub>P</sub>) mathematically, and quantify how a<sub>P</sub> varies with soil physical and chemical properties, and time scale.

## 2 Material and Methods

Accepted Article
## 2.1 A model of inorganic soil phosphorus dynamics

Here we focus on the dynamics of inorganic P in soil and omit biological processes, such as plant uptake, microbial mineralization and immobilization as well as cleavage of organic phosphorus through phosphatase or organic acids. We assume a closed system without any external inputs or loss. Soil inorganic P is partitioned into different pools based on Hedley fractions, which are: solution inorganic P, labile inorganic P, sorbed inorganic P and occluded inorganic P. Solution inorganic P refers to water-extractable inorganic P, labile inorganic P as the resin- and bicarbonate-extractable inorganic P minus water-extractable inorganic P, sorbed inorganic P as NaOH-extractable inorganic P, and occluded inorganic P as residual inorganic P. HCl-extracted P, which includes P in apatite and some secondary compounds, is not considered (see Figure 1). The HCl-extracted P pool that can significantly contribute to biological P uptake in some P-poor soils (see Lambers et al. 2008) is excluded here, because data availability is yet too scarce to reliably extend the model to include HCl-P. In addition, data-driven formulations for P release from primary minerals for global modelling exist (e.g. Hartmann et al 2014) which can be coupled to our model within a land model. This study focuses on soil inorganic P, thereon P is used for soil inorganic phosphorus unless stated otherwise.

Figure 1 shows the exchanges between the solution and labile, or sorbed or between sorbed and occluded P pool. Exchange between solution P and labile or sorbed P pools are mediated by sorption and desorption. Sorption of solution P into labile P or sorbed P is described using a Freundlich equation (see Barrow 2008). That is

$$F_{LW} = k_{LW} C_{W}^{b} \qquad (1)$$

$$F_{SW} = k_{SW} C_{W}^{b} \qquad (2)$$

Here we use $F$ for flux in mg P (kg soil)$^{-1}$ day$^{-1}$, with the subscript LW in $F_{LW}$ representing a flow from solution P pool (W) to labile P pool (L), $F_{SW}$ is defined similarly. $C_{W}$ is the solution P concentration (mg P/L) and is time-dependent, $b$ is a dimensionless parameter (Barrow 2008), and is related to the chemical potential of phosphate in solution (Shayan and Davey 1978), varies between 0.1 and 2 for most soils. Parameter $k_{LW}$ and $k_{SW}$ are sorption coefficients for labile and sorbed P, respectively (see Chen et al. 1996); and both have a unit of mg P (kg soil)$^{-1}$ day$^{-1}$ (mg P/L)$^{-b}$.
Lookman et al. (1995) found that a double-exponential function adequately described the desorption of soil phosphorous, which is equivalent to a model of two pools with the first-order kinetics for desorption from each pool. Here we use the first-order kinetics to describe the desorption labile or sorbed P into solution. Desorption rate constants are different between labile and sorbed pools. That is

$$F_{WL} = k_{WL}P_{L} \qquad (3)$$

$$F_{WS} = k_{WS}P_{S} \qquad (4)$$

where $P_L$ and $P_S$ are the size of labile and sorbed P in mg P (kg soil)$^{-1}$, $k_{WL}$ and $k_{WS}$ are two desorption rate constants in day$^{-1}$, and $F_{WL}$ and $F_{WS}$ are the rates of desorption from labile and sorbed P pools in mg P (kg soil)$^{-1}$ day$^{-1}$.

Equations (1) to (4) were previously used by Chen et al. (1996) for describing sorption and desorption of soil phosphorus, and were also discussed by McGechan and Lewis (2002) (their equations 66 and 67).

Sorption of inorganic P into occluded P is quite slow and quite complicated (McGechan and Lewis 2002). For the sake of practicability, here we assume that sorption into the occluded P ($F_{OS}$) and desorption from occluded P ($F_{SO}$) are proportional to their respective pool sizes.

They are:

$$F_{OS} = k_{OS}P_{S} \qquad (5)$$

and

$$F_{SO} = k_{SO}P_{O} \qquad (6)$$

Dynamics of the closed system of four inorganic soil P pools is given by

$$\frac{dP_{W}}{dt} = \underbrace{k_{WL}P_{L} + k_{WS}P_{S}}_{\text{desorption}} - \underbrace{(k_{LW} + k_{SW})C_{W}}_{\text{sorption}} \qquad (7)$$

$$\frac{dP_{L}}{dt} = k_{LW}C_{W}^{b} - k_{WL}P_{L} \qquad (8)$$

$$\frac{dP_{S}}{dt} = k_{SW}C_{W}^{b} + k_{SO}P_{O} - k_{WS}P_{S} - k_{OS}P_{S} \qquad (9)$$

$$\frac{dP_{O}}{dt} = k_{OS}P_{S} - k_{SO}P_{O} \qquad (10)$$

with $P_W = vC_W$.
The model of inorganic P dynamics consists of eqn 7 to 10 with four state variables ($P_w$, $P_L$, $P_S$ and $P_O$) and seven parameters, $k_{WL}$, $k_{WS}$, $k_{LW}$, $k_{SW}$, $k_{OS}$, $k_{SO}$ and $b$, and two variables $v$ and $t$. $C_\infty$ is the value of $C_w$ at steady state, $v$ is soil water content in L/kg soil and at 10 L/kg, and $P_w+P_L+P_S+P_O=P_{in}$, where $P_{in}$ is the total inorganic P in mg/kg soil. Table 1 lists the definitions and units of all symbols. The model is linear when $b=1$, and nonlinear otherwise.

Exchange is relatively fast between labile and solution P, and slow between sorbed and solution P, or $k_{LW} \gg k_{SW}$, and $k_{WL} \gg k_{WS}$, the model as represented by eqn 7-10 can simulate both fast and slow processes governing soil inorganic P dynamics, as discussed by MeGechan and Lewis (2002), and Barrow (2008).

If exchanges between sorbed and occluded P pools are ignored, the steady state solution to eqn (7), (8) and (9) will give the following relationship:

$$P_L^* + P_S^* = \left( \frac{k_{LW}}{k_{WL}} + \frac{k_{SW}}{k_{WS}} \right) (C_W^*)^b \quad (11)$$

where $C_W^*$, $P_L^*$ and $P_S^*$ are solution P concentration, labile P and sorbed P at steady state. Eqn (11) is the Freundlich equation.

## 2.2 Isotopic exchange kinetics (IEK) experiments

We used the parameter values derived from isotopic exchange kinetics (IEK) experiments to calculate the time course of the depletion of added inorganic P in soil solution and measured fractions of inorganic P in different pools to estimate those seven model parameters, $k_{WL}$, $k_{WS}$, $k_{LW}$, $k_{SW}$, $k_{OS}$, $k_{SO}$ and $b$.

In an IEK experiment, plant roots and mycorrhizae in the soil samples were excluded, soil samples were then sterilized to inhibit microbial activities (see Achat et al. 2011). As a consequence biological P transformations were absent, in line with our modelling assumption. All measurements were conducted in laboratories in a soil-solution system with soil/solution ratio of 1:10 (see Fardeau et al. 1991). A very small amount of radioisotope ($^{32}$P or $^{33}$P) is added to soil solution at time $t=0$. The measured radioactivity remaining in soil solution decreases with time as a result of sorption processes. The data from the IEK experiments are used to fit the following equation (see Fardeau et al. 1991):

$$\frac{r(t)}{R} = m[t + m^{1/n}]^{-n} + \frac{r(\infty)}{R} \quad (11)$$
where $r(t)$ is the radioactivity in Bq measured in the solution at time $t$ after the addition (in minute), and $R$ is the total amount of radioactivity added at time $t=0$, $m$ and $n$ are two dimensionless model parameters that are estimated by fitting eqn (11) to the measured data. The term $r(\infty)/R$, represents the value of $r(t)/R$ at steady state, and is usually approximated by the ratio of solution inorganic P ($P_W$) to total soil inorganic P ($P_T$).

However, at $t=0$, eqn(11) gives: $\frac{r(t=0)}{R}=1+\frac{r(\infty)}{R}$. To correct for this small discrepancy, we modified eqn (11) as follows to make the equation correct also for $t = 0$:

$$ \frac{r(t)}{R} = m \left(1 - \frac{r(\infty)}{R}\right) \left[t + m^{1/n}\right]^{-n} + \frac{r(\infty)}{R} \quad (12) $$

It should be noted that adding a small dose of radioactive P to solution has no effect on the total inorganic P in solution, or $P_W$, since the added P mass is negligible (Frossard et al. 2011).

In the IEK data compiled (Helfenstein et al. 2020), parameters $m$ and $n$ were estimated by fitting eqn (11) rather than eqn (12) to the measurements. As $\frac{r(\infty)}{R}$ is usually less than 0.001, estimates of $m$ and $n$ should vary very little whether eqn (11) or eqn (12) was fitted to the measurements.

In theory, parameters $m$ and $n$ broadly account for fast and slow physical-chemical reactions, respectively (Fardeau et al., 1991; Frossard et al., 2011). In some studies, a value of $r/R$ at $t=1$ minute after the addition is used for $m$. The value of parameter $n$ usually varies between 0.1 and 0.7 (Achat et al. 2016b). Desorption of inorganic P into solution is generally slow for soil with a small $n$ value. This study did not fit the equation (11) to the isotopic data, but used the parameter values derived from the previous isotopic studies to calculate how $r/R$ varies with time. The calculated values of $r/R$ were then used to constrain the estimates of parameters in our model of inorganic P dynamics (eqns 13 to 16, see the below).2.3 Modelling the dynamics of added inorganic phosphorus in soil solution

To use $r(t)/R$ to constrain the parameters in the model of inorganic P dynamics, we further develop a model for simulating the dynamics of added inorganic P to soil solution, or. This requires making the following assumptions: (1) prior P addition, all inorganic pools are in a steady state; (2) when a new steady state is reached after P addition, the fractions of added P in different inorganic P pools are equal to their respective fractions of inorganic P in different pools before the addition.
As shown in Appendix 1, the equations governing the dynamics of the added P at any time t after addition are:

$$ \frac{dr_W}{dt} = k_{WL}r_L + k_{WS}r_S - (k_{LW} + k_{SW})(bC_{\infty}^{b-1}\rho_W) \quad (13) $$

$$ \frac{dr_L}{dt} = k_{LW}(bC_{\infty}^{b-1}\rho_W) - k_{WL}r_L \quad (14) $$

$$ \frac{dr_S}{dt} = k_{SW}(bC_{\infty}^{b-1}\rho_W) + k_{SO}r_o - k_{WS}r_S - k_{OS}r_S \quad (15) $$

$$ \frac{dr_o}{dt} = k_{OS}r_S - k_{SO}r_o \quad (16) $$

with initial condition of $r_W = 1$, and $r_L = r_S = r_o = 0$ at $t=0$.

Where $r_W$, $r_L$, $r_S$ and $r_o$ are modelled fractions of isotopically labelled inorganic P in solution, labile, sorbed and occluded pools at time t after addition, respectively. Therefore $r_W=r(t)/R$.

At steady state, $r_W(t \to \infty) = f_W$, $r_L(t \to \infty) = f_L$, $r_S(t \to \infty) = f_S$, and $r_o(t \to \infty) = f_o$, where $f_i$ is the steady state size fraction of pool i ($i$ = W, L, S or O) and $f_W+f_L+f_S+f_o=1$ (see the Appendix 1). The values of $k_{LW}/k_{WL}$, $k_{SW}/k_{WS}$ and $k_{OW}$ are determined by the inorganic P fractions given by

$$ \frac{k_{LW}}{k_{WL}} = \frac{f_L}{bC_{\infty}^{b-1}f_W/v}, \quad (17) $$

$$ \frac{k_{SW}}{k_{WS}} = \frac{f_S}{bC_{\infty}^{b-1}f_W/v}, \quad (18) $$

$$ \frac{k_{SO}}{k_{OS}} = \frac{f_S}{f_o}, \quad (19) $$

As the uncertainty of the calculated the fraction of added P remaining in solution using eqn (12) is likely to be large, and will not provide a good constraint on the estimates of two rate constants for the occluded P pool ($k_{OS}$, $k_{SO}$) that will be of the order of $10^{-1}$ to $10^{-2}$ year$^{-1}$ (see Helfenstein et al. 2020), we fixed $k_{SO}$ at 0.02 year$^{-1}$, a value used in CASA-CNP (Wang et al. 2010) and calculated $k_{OS}$ using eqn (19).

Only three of the seven parameters in the inorganic P model, $k_{WL}$, $k_{WS}$ and $b$ are to be estimated, $k_{LW}$, $k_{SW}$ and $k_{OS}$ can be calculated using eqns (17)-(19).

## 2.4 Soil P bioavailability at different timescale
Following Buehler et al. (2003), we define soil P bioavailability as the fraction of inorganic P that can leave the solid phase of the soil and arrive in the soil solution for uptake by plants or soil microbes within a given time. As biological uptake drives down the solution P, desorption of P from labile and sorbed pools will replenish solution P and provide additional P available for biological uptake. That additional available P, A, can be calculated as

$$A(x, t) = P_L(0,0) - P_L(x, t) + P_S(0,0) - P_S(x, t) \quad (20)$$

and

$$P_L(0,0) - P_L(x, t) = (1 - \exp(-k_{WL}t))\Delta P_L \quad (21)$$

$$P_S(0,0) - P_S(x, t) = (1 - \exp(-k_{WS}t))\Delta P_S \quad (22)$$

Here we ignored the contribution of desorption of occluded P to soil solution via the sorbed P pool. In practice, occluded P can be desorbed directly into soil solution (see Schubert, Steffens and Ashraf 2020), and that rate usually is quite small, and is therefore excluded here, given the duration of IEK experiments being a few hours to days and the limited data that hamper the parameterization of additional model parameters. Where x is the fraction of solution P concentration reduced by biological P uptake, and t is the time over which the additional P becomes available by desorption. ΔP_L and ΔP_S are the difference in the steady state pool sizes of labile and sorbed P pools, respectively, and are given by

$$\Delta P_L = \frac{k_{LW}}{k_{WL}} (1 - x^b) C_\infty^b = (1 - x^b) f_L P_T \quad (23)$$

$$\Delta P_S = \frac{k_{SW}}{k_{WS}} (1 - x^b) C_\infty^b = (1 - x^b) f_S P_T \quad (24)$$

when b=1, we have ΔP_L = (1 - x)f_L P_T and ΔP_S = (1 - x)f_S P_T.

The soil P bioavailability, a_p, is calculated as

$$a_p = \frac{A(x,t)}{P_T} = (1 - \exp(-k_{WL}t))(1 - x^b)f_L + (1 - \exp(-k_{WS}t))(1 - x^b)f_S \quad (25).$$

Because $x \le 1$, $b>0$, then $(1 - x^b)$ increases with an increase in b or a decrease in x.

Augusto et al. (2017) defined the available inorganic P in soil for plant uptake as the sum of labile and sorbed P (see their eqn 11). That is equivalent to the soil P bioavailability calculated using eqn 25 with $t \to \infty$, $x=0$.
## 2.5 Model parameter optimization

To optimize the three unknown model parameters $k_{WL}$, $k_{WS}$ and $b$ for each soil sample, we minimized the squared differences between the modelled fraction of added P remaining in solution ($r_W$) and the calculated value of $r(t)/R$ using eqn (12) with the published values of $m$, $n$ and $r(\infty)/R$ at different times, and between the simulated fractions of added P in different inorganic P pools at steady state and the measured fractions of different inorganic P pools. The cost function ($J$) for parameter optimization was constructed as

$$J = J_1 + J_2 \qquad (26)$$

And

$$J_1 = \sum_{hour=1,24} (r_W - \frac{r}{R})^2 + \sum_{day=2,365} (r_W - \frac{r}{R})^2 + \sum_{year=2,20} (r_W - \frac{r}{R})^2 \quad (27)$$ $$J_2 = (f_L - f_L^*)^2 + (f_S - f_S^*)^2 + (f_O - f_O^*)^2 \qquad (27)$$

Here $J_1$ is calculated as the sum of squared differences for the hourly values of the first 24 hours after labelling, and daily values of first 365 days and yearly values of the first 20 years. $f_L$, $f_S$ and $f_O$ are also taken as the steady state values of $r_L$, $r_S$ and $r_O$ obtained by setting $dr_W/dt = dr_L/dt = dr_S/dt = dr_O/dt =0$ in eqns (13) to (16), respectively. $f_L^*$, $f_S^*$ and $f_O^*$ are the measured fractions of labile, sorbed and occluded P, respectively.

To estimate $b$, $k_{WL}$ and $k_{WS}$, we used a highly efficient and effective optimization method based on shuffled complex evolution, or SCE-UA (see Duan et al. 1993) to minimize $J$ in eqn (26). This method is highly robust and efficient (Duan 2003), and is well suited for estimating parameters in highly nonlinear models (see Wang et al. 2021 for example). As this study is to quantify how model parameter values vary with soil physical and chemical properties, and the uncertainties of the observed fractions of different inorganic P pools are poorly quantified, we therefore only estimated the values of model parameters, rather than the posterior parameter distributions.

## 2.6 Data for model calibration

We used the data compiled by Helfenstein et al. (2020) for 45 soil samples and 102 French soil samples compiled by Achat et al. (2016b). For each soil sample, the data included estimates of $r(\infty)/R$, $m$ and $n$ from the IEK experiments, together with measurements of $C_\infty$ (in mg/L), total inorganic soil P (mg P/kg). For the 45 soil samples compiled by Helfenstein et al. (2002), measurements of the amounts of inorganic P extracted by resin (resin-P) (solution
P), NaHCO₃ (NaHCO₃-P) (labile P), NaOH (NAOH-P) (sorbed P) and HCl (HCl-P) (primary P), and the amount of residual P in mg P/kg were made by Helfenstein et al. (2018b) or extracted from the published papers by Helfenstein et al. (2020). As inorganic P fractions were not measured for the 102 French forest soils, we used estimates based on a machine-learning approach by He et al. (in preparation) that explains 66% of the variance in the observed P fractions among 5275 soil samples representing all major vegetation and soil types in the world.

We partitioned the residual P into an organic and inorganic fraction for all 147 soil samples using an empirical model based on the soil phosphorus data compiled in Augusto et al. (2017). The inorganic P in the residual P is the occluded P. The fractions of solution, labile, sorbed and occluded inorganic P were then calculated for each of 147 soil samples.

Data of 45 soil samples compiled by Helfenstein et al. (2020) were from the following four studies. The first study is the one by Helfenstein et al. (2018b) who collected soil samples at two different soil depths of six sites along a precipitation gradient in Hawaii. The second study by Chen et al. (2003) measured 15 different grassland soils in New Zealand. The third study by Lang et al. (2017) measured the soils at two different soil depths from five mature beech forests in Germany along a geo-sequence (different parental material). The fourth study by Buehler et al. (2002, 2003) measured highly weathered soils with very low P under different land uses in Columbia. Together with the data from 102 French forest soils (Achat et al. 2016b), we have 147 soil samples representing a wide range of climate, soil physical and chemical properties (see Figure S1). Table S1 lists the ranges of variations for those variables or soil properties of the 147 soil samples.

## 2.7 Multivariate regression between the model parameters and soil properties and climate variables

To relate the soil specific optimized values of model parameters to soil properties (see below) and climate variables (mean annual air temperature *T*<sub>a</sub> in °C or mean annual precipitation *P*<sub>a</sub> in mm/year), we fitted a stepwise multivariate linear regression using python v3.8.6 statsmodels package. The independent variables we used include: climate variables (*T*<sub>a</sub> and *P*<sub>a</sub>), soil texture properties (fractions of clay, silt or sand), USDA soil order (see www.nrcs.usda.gov), soil chemical properties (organic carbon, total inorganic P and its fractions in soil water, labile, sorbed and occluded pools, soil pH as measured in water, oxalate extractable metal (Al and Fe) oxide concentration). In the final regression, only
independent variables that make significant contributions (p <0.05) to the explained variables were included.

## 2.8 Importance of climate, soil physical and chemical properties for soil P bioavailability

Following Buehler et al. (2003), here we define soil P bioavailability ($a_P$) as the fraction of inorganic P that can leave the solid phase of the soil and arrive in the soil solution for taking up by plants or soil microbes within a given time. As Eqn (25) shows that, $a_P$ depends on the rate constants of desorption, parameter $b$, $C_\infty$, soil P fractions and time.

Using the multivariate regression equations listed in Table 2, we calculated the rate constants ($k_{WL}$, $k_{SW}$), $b$ or $C_\infty$ and $k_{LW}$ and $k_{WS}$ using eqns 17-18, then $a_P$ using eqn (25) for given values of soil physical and chemical properties. To study the importance of the different soil properties for $a_P$, we used the Morris method to compute the elemental effects of different soil properties on $a_P$. Following Lu et al. (2013), the importance of $T_a$, or soil property, $\theta_j$, was calculated as

$$ \beta_j = \frac{\sqrt{\mu_i^2 + \sigma_i^2}}{\sigma} \qquad (28) $$

where $\mu_i$ and $\sigma_i$ are the mean and standard deviation of the elemental effect of $T_a$ or soil property on $a_P$, and $\sigma^2 = \sum_{1}^{N} \sigma_i^2$, where N is the number of independent variables (N=13 in this study, e.g. all variables are listed in Table S1 except mean annual precipitation and soil order).

The elemental effect of $x_i$ ($x_i=T_a$, or soil property) is calculated as

$$ E_i = \frac{a_P(x_1, x_2, \dots, x_i + \Delta, \dots, x_k) - a_P(x_1, x_2, \dots, x_i, \dots, x_k)}{\Delta} \qquad (29) $$

where $\Delta$ is the step size of relative change in $x_i$ within its range. See Lu et al. (2013) and Saltelli et al. (2004) for further details.

## 3 Results

In the following we will present the results of model calibration, and how the estimated model parameters vary with IEK parameters ($m$, $n$ and $r(\infty)/R$), or soil properties. Then we will estimate how the soil P bioavailability ($a_P$) varies with soil physical and chemical properties, and the importance of different soil properties on the fraction of $a_P$ at different time scales.
## 3.1 Performance of the calibrated model

We first compared the performance of linear and nonlinear inorganic P models with optimized model parameters. Model equations for the linear model are the same as the nonlinear model except that *b* is set to 1 (see eqns 1 and 2). Values of the optimized parameters (*b*, *k*<sub>WL</sub>, *k*<sub>WS</sub>) differ between the linear and nonlinear models for most soils (data not shown). Figure 2a shows that the linear model has larger RMSEs than the nonlinear model for most of the 147 soil samples. Our results support the use of the nonlinear model to describe P sorption in soil in line with previous findings (Barrow 2008). In the following, only the results of the nonlinear P model are presented.

As shown in Figure 2b, RMSE of the calibrated model is less than 0.003 for 35 soil samples, exceeding 0.1 for only 17 out of the 147 soil samples. The model performance as measured by $r^2$ is good, with $r^2$ being greater for 0.7 for 122 soils and <0.5 for only 7 soils (see Figure 2c). Further analysis shows that large RMSE are concentrated at soils with small (<0.2) values of *n* or large (>0.8) values of *m* (see Figure 2d).

We also assessed the impact of three different fixed values of $k_{SO}$ (0.01 year$^{-1}$, 0.02 year$^{-1}$ and 0.005 year$^{-1}$) on the estimated parameters, $k_{WL}$, $k_{WS}$ and *b*, and found that impact were small (<5%).

We explored the underlying cause for the relatively larger RMSE for soils with small values of *n*. Parameter *n* in the IEK equation is closely related to slow processes controlling sorption and desorption of P (Fardeau 1993). When *n* is close to 0, the rate of sorption or desorption is very slow, therefore such a system takes a very long time (i.e. > decade) to reach steady state. To explore why the model did not fit the IEK data well for soils with low *n*, we compared the modelled time course of *r/R* with that calculated using eqn (12) for two soils with different values of *n*. As shown in Figure 3, sorption and desorption take place rapidly for soil 1 with *n*=0.437, as indicated by the time reaching their respective equilibrium of all the four pools being < four days (see Figure 3c), and the model fits very well the calculated values of *r/R* using IEK equation (see Figure 3a). The differences in the modelled and measured fractions of different inorganic pools are negligible for soil 1, with high *n* value. For soil sample 19 with a low value of *n* (0.12) and very high fraction of sand (84% in texture), the calculated value of *r/R* continues to decline slowly after 10 years (see Figure 3b), which was poorly fit by our model (see Figure 3d).
An indication of a poor model fit is a large difference between the modelled and measured
fractions of different P pools at steady state. The differences in the modelled and measured
pool size fractions are negligible (<0.001) for the high-n soil, but quite large for the low-n
soil (Figure 3). The fraction of inorganic P pools for the low-n soil as modelled at steady state
are 0.018, 0.12, 0.22 and 0.65, as compared the observed fractions of 0.046, 0.30, 0.56 and
0.1 for the fractions of inorganic P in solution, labile, sorbed and occluded pools,
respectively. The model underestimates the fractions of labile and sorbed pools, but
overestimates the fraction of the occluded P for low-n soils.

## 3.2 Dependence of the kinetic parameters in the inorganic P model on isotopic exchange kinetics

Among the six model parameters, the parameters for sorption are related to the ones for
desorption in case of labile (k<sub>WL</sub>, k<sub>LW</sub>) and sorbed (k<sub>WS</sub>, k<sub>SW</sub>) through eqns (17) and (18). The
optimized sorption coefficient being > 100 mg P (kg soil)<sup>-1</sup> day<sup>-1</sup> (mg P/L)<sup>-b</sup> for labile P, and
between 1 and 20 mg P (kg soil)<sup>-1</sup> day<sup>-1</sup> (mg P/L)<sup>-b</sup> for most soils (see Figure S2a), the
desorption rates of labile P (k<sub>WL</sub>) vary mainly between 0.01 to 10 day<sup>-1</sup> among soil samples,
whereas the desorption rates for sorbed P range between 0.001 to 1 day<sup>-1</sup> (Figure S2b). The
sorption rates of occluded P vary across a larger range than sorption rates of other pools, i.e.
from 10<sup>-6</sup> to 10<sup>-3</sup> day<sup>-1</sup> (Figure S2b).

In theory, parameters *m* and *n* in the eqn 12 are broadly related to the fast and slow rates of
processes controlling sorption and desorption of P, respectively (Fardeau 1993). In our model
of P dynamics, labile and sorbed P are used to represent P that is exchanged with the solution
P at hourly to daily (1/k<sub>WL</sub> <1 day) and weekly to yearly (1/k<sub>WS</sub> > 7 days) timescale for most
soil samples (Figure S2), therefore their respective sorption/desorption rates should be related
to parameters *m* and *n*, respectively. This is also supported by our results. Although the
correlation is not statistically significant between sorption/desorption rate constant of labile P
and *m* (data not shown), the desorption/sorption rate constant decreases for labile P, but
increases for sorbed P with *n* (see Figure 4a and 4b).

Different from the desorption of sorbed P into soil solution, Figure 4c shows that the rate
constant of transfer from sorbed to occluded P pool per unit fraction of occluded P (k<sub>OS</sub>/f<sub>0</sub>)
decreases with an increase in *n* (*r*<sup>2</sup>=0.1). Parameter *b* in the Freundlich equation for sorption
decreases with an increase in *n* (see Figure 4d), and is <1 for most soil samples, therefore the
sorption rate into labile or sorbed P pool increases with solution P concentration nonlinearly.
The sorption coefficient in the Freundlich equation is poorly correlated with *m* or *n* for the labile P, but significantly decreases with *m*, or increases with *n* for sorbed P (see Figure S3). Parameters *m* and *n* are negatively correlated (see Figure S4a). Therefore, a soil with a lower value of *n* or higher value of *m* usually has faster sorption and desorption rates of labile P, slower sorption and desorption rates of sorbed P, as shown in Figure 2 for soil sample 19. As a result, both *C*∞ and *r*(∞)/R significantly decrease with an increase in *n* (see Figure S4).

## 3.2 Variations of model parameters with soil properties

One objective of this study was to develop relationships between model parameters and soil properties or climate as the first step towards estimating parameter values at regional or global scales in global land modelling.

One such parameter is the solution P concentration at steady state (*C*∞). Figure 5a shows *C*∞ increases exponentially with C/Oₓ but reaches a saturation around C/Oₓ ~ 2 g C/mmol. The best fitted regression explains 54% of the variance in the observed *C*∞ (see Figure 5b). We also fitted regression equations for the other four model parameters, *b*, *k*WL, *k*SW, and *k*OS (see Figure S5 and Table 2). All those four parameters depend on the fractions of inorganic P pools. In addition, *k*SW varies with *T*ₐ, and *k*OS varies with soil texture and soil pH (see Table 2). For *k*WL, *k*SW and *b*, the best fitted regressions overestimate the parameter values when parameter values are low, and underestimate the parameter values otherwise (see Figure S5). The best fitted regression predicts the optimized value of *k*OS very well (see Figure S5).

## 3.3 Soil P bioavailability

Different from previous studies, our definition of soil available inorganic P (*a*<sub>P</sub>) is based on a theoretical model which goes beyond existing approaches which rely solely on data from soil P fractionation, whereas our approach also takes into consideration the turnover of soil P fractions. To assess the plausibility of *a*<sub>P</sub> as defined in this study, we compared the estimated *a*<sub>P</sub> based on our definition with those by Buehler et al. (2003) for the soils from four different land uses in Carmagua, Columbia.

Buehler et al. (2002) measured the soil P fraction and obtained IEK kinetics (*m*, *n*, *r*(∞)/R) from the soils from four different land uses in Carimagua, Colombia. The data were used as part of the 147 soils for our model calibration. Buehler et al. (2003) used the same soils to
measure the uptake of isotopically labelled P by plants from the harvest biomass at 12 weeks
after adding the isotopically labelled P. They also compared the measured isotopically
labelled P taken up by plants from biomass harvest (L) with the amount of exchangeable P as
calculated using the IEK kinetics (E), found that E is about 20% greater than L. Given the
likely experimental errors (not provided by Buehler et al. (2003)), E can be considered to be
not significantly different from L.

Using the eqn (25) and the optimized parameters for the four soils (see Table 3), we
estimated the available soil P (A) as a<sub>p</sub>P<sub>T</sub> for plants grown on each of the four soils over the
12-week period. Our estimates of available soil P agree with the measured values for the four
soils within ±10% (see Table 3).

Furthermore, we explore why *a*<sub>p</sub> differs among the four different land uses at Carmagua,
Columbia. As Eqn 25 shows, *a*<sub>p</sub> depends on desorption rate, *b* and soil P fractions. The terms
in eqn 25, (1-exp(-*k*<sub>WL</sub>*t*)) and (1-exp(-*k*<sub>WST</sub>*t*)) are very close to 1 at the four sites except (1-
exp(-*k*<sub>WST</sub>*t*)) for RGM site (=0.89), therefore available soil P is largely dependent on soil P
fractions and parameter *b*. The two sites, SAV and GL have lower amount of inorganic P and
smaller fractions of labile and sorbed P, therefore lower amount of soil available P and soil P
bioavailability than the other two sites. Between the two fertile sites, CR and RGM, the total
fractions of labile and sorbed P are quite similar (0.68 at CR, and 0.62 at RGM), but the soil
available P at CR was lower than RGM (see Table 3). This difference is largely contributed
by the difference in *b* between the two sites: (1-*x*<sup>b</sup>)=0.84 at CR and =1 at RGM.

To assess how important different soil properties (physical and chemical) and climate (T_a) on
a_p are, we calculated the sorption coefficients, desorption rate constant for labile and sorbed
P, and sorption rate constant for occluded P, b and c_w at steady state for a short timescale
(t=1 day) or a growing season (t=180 days). We quantified the importance of a given soil
property or T_a using the Morris method, and the results are shown in Figure 6.

Among the 14 variables (see Table S1), relative sensitivities are >0.01 or <-0.01 for 10
variables (see Figure 6a). A negative sensitivity indicates that a_p decreases with an increase in
that variable. Among those 10 variables, a_p decreases with an increase in oxalate extractable
metal oxide concentration, or clay percentage the amount of inorganic P, and fraction of
inorganic P in soil solution, and increases with an increase in all other eight variables (see
Figure 6a). The absolute value of the relative sensitivity is greatest for f<sub>W</sub><sup>*</sup>, and quite small
(<0.1) for soil pH (not shown), T_a and soil texture variables. The absolute relative sensitivity
increased for the fraction of solution P, Tₐ and decreased for sand fraction, when t is varied from 1 day to 180 days, suggesting an increase in importance of soil chemical properties (oxalate extractable metal oxide concentration, inorganic P fractions).

As shown in Figure 6b, fractions of inorganic P pools are more important than all other variables at both timescales, and soil pH and soil texture percentage are the least important variables. The importance of oxalate extractable metal oxide concentration, soil P fractions increases when t is varied from 1 day to 180 days, which is consistent with the responses of their relative sensitivities to t.

## 4. Discussion

A process-based model of inorganic P dynamics was developed in this study. This model uses the Freundlich equation to represent inorganic P sorption, which was considered to be more suitable than the widely used Langmuir equation in global land models (Barrow 2008). We also showed that the simulations by the model with the Freundlich equation (b≠1) agreed much better than the model applying a linear equation. Therefore, nonlinear sorption of soil inorganic P should be used in global land models with P cycle.

Using the measured fractions of different inorganic P pools and results from IEK experiments from 147 soils worldwide, we found that the model simulated the measurements quite well for most soil samples, except those sandy soils (sand fraction >0.6) with low b (<0.2) and very slow exchanges between the different inorganic P pools. A low value of b corresponds to small sensitivity of sorption rate to solution P concentration based on Freundlich equation. For those sandy soils, our model tends to simulate a faster decline in the first few hours, then slower decline in the added P in soil solution than the IEK model suggests. Those biases likely result from errors in model structure, further studies are needed to identify the source of errors and develop a better model for those sandy soils.

### 4.1 Soil properties controlling soil inorganic P transformation

The optimized model parameters vary widely among the different soil samples, and those variations can be well explained empirically by soil physical and chemical properties and Tₐ (see Table 2). Those empirical relationships can be used to estimate model parameters for studying soil P cycle at regional or global scales. One of those model parameters is b in the Freundlich equation. Theoretically, as soil pH increases, and the electrostatic potential
decreases, parameter *b* should increase (see Figure 1 of Barrow 2008). The value of *b* is <1 for neutral or acidic soils (Barrow 1983), can be >1 for alkaline soils (see Bertrand et al. 2003). Among the 147 soil samples, only 7 soils were with pH >7, and the value of *b* for those 7 soil samples varied from 0.2 to 1.6 with a mean value of 0.7. We found that the correlation between *b* and soil pH was not statistically significantly (*r*²=0.001). More data from alkaline soils are needed to assess the dependence of *b* on soil pH.

This study found that sorption coefficient ($k_{LW}$) and desorption rate constant ($k_{WL}$) of labile P did not vary with $T_a$, whereas sorption coefficient or desorption rate constant of sorbed pool did. Theoretically, sorption of inorganic P in soil varies with temperature with an activation energy of 90.7 kJ/mol (see Barrow 2008), which would give an exponential dependence on temperature (or $y=\exp(\alpha T_a)$) with $\alpha$ value of 0.091, which is close to the value of 0.108 from our regression (see Table 2). Here $T_a$ represents the mean annual temperature of the sampling sites rather than the ambient temperature at which IEK experiments were conducted, therefore sensitivity of sorption coefficient into sorbed P ($k_{SW}$) to $T_a$ may include the effects of different environmental conditions (e.g. climate, soil, lithology and ecosystem), which requires further studies.

Theoretically, desorption of adsorbed or fixed P occurs mostly via ligand exchange reaction (Hinsinger 2001), which was not explicitly represented in the empirical relationships for $k_{WL}$ or $k_{WS}$ (see Table 2). Instead, the empirical relationship indicates that the desorption rate constant for labile P ($k_{WL}$) increases with a decrease in the labile fraction ($f_L^*=1-f_W^* - f_O^*$), and that the desorption rate for sorbed pool ($k_{WS}$) that is proportional to $k_{SW}$ increases with fractions of sorbed P ($f_S^*$).

Phosphate can be adsorbed to the surface of metal oxides, clay minerals and organic matter (Hinsinger 2001). The empirical equations in Table 2 predict that the sorption rate into the sorbed and occluded P increases with metal oxide concentration and total soil organic C (see Table 2), which is in line with previous studies (see Hinsinger 2001, Achat et al. 2016b, and Borggaard et al. 2004). Increasing soil carbon concentration could increase P sorption indirectly by inhibiting iron oxide crystallisation (Borggaard 1986) or by increasing the concentration of carboxyl, a strong sorbent of inorganic P in soil (Hinsinger 2001). Such effects could explain the positive dependence of P sorption into the occluded pool on the soil C concentration we found. On the other hand, increasing soil C concentration can reduce P sorption by increasing competition of organic anion with inorganic P for reactive surface (see

Accepted Article
Regelink et al. 2015), which explains the negative sensitivity of sorption of labile P, $k_{LW}$ that is proportional to $k_{WL}$ by eqn (17) to soil C concentration (see Table 2).

Although an empirical relationship (see Table 2) can explain 90% of the variance of the estimated desorption rate of the occluded P (see Table 2), the uncertainties are likely to be large, as we used the calculated fraction of added P remaining in soil solution with the parameter values from the relatively short-term (hours to days) IEK measurements. The rate constant of sorption and desorption of occluded P are generally about $10^{-5}$ day$^{-1}$, will therefore not be significantly constrained by the calculated fractions using eqn (12).

Finally, contrary to the assumption made in global land models with P cycle (see Wang et al. 2010 for example), our results show that parameters controlling sorption and desorption dynamics vary by several orders of magnitude, and metal oxide and other soil chemical properties rather than soil order are significant contributors to those variations. Our results consequently do not support using parameter values based on soil order for modelling P cycle in global land models.

## 4.2 Soil P bioavailability: new insights from this study

The quantification of soil bioavailable P on time scales relevant for carbon and nutrient cycles which span from days to decades is challenging. Here we defined the bioavailable fraction as the inorganic P which is exchangeable in both solid and solution phases for a given time. Our definition is similar to that by Hamon et al. (2002) and Buehler et al. (2003) but differs significantly from other previous studies which are often based on one or several soil P fractions on various time scales (e.g. Ringeval et al. 2017, Augusto et al. 2017, Sun et al 2017). Our estimate of soil P bioavailability depends not only on P fractions, but also on the desorption rate, parameter *b* in the Freundlich equation, and on the extent to which active biological P uptake can drive down the solution P concentration. Our theoretical derivation of available P is supported by Barros et al. (2005) who found that available inorganic P in highly weathered soils in Brazil strongly depended on the sorption and desorption characteristics of the soils. Soil P fractions alone are poor indicators for bioavailability (see also Johnson et al. 2003).

In addition to the fractions of different soil P pools, soil P bioavailability increases with soil C concentration but decreases with an increase in oxalate extractable metal oxide concentration (see Figure 6). This is consistent with the results of Achat et al. (2016b), who identified the importance of soil C concentration and metal oxide concentration controlling

Accepted Article
inorganic P bioavailability. This is not unexpected, as two thirds of data used in this study were from Achat et al. (2016b). But, by including the additional data (n=45), we confirmed the findings of Achat et al. (2016b).

Different from previous studies, soil P bioavailability as defined in this study also depends on the capacity of biological uptake to drive down the solution P concentration (x) (see eqn 25). This is supported by Hinsinger (2001). Therefore, soil available P is not just a soil property, as it also depends on plant types. The transfer of inorganic dissolved P from the surrounding soil towards the rhizosphere is dominated by diffusion (see Barber 1995), and is a major constraint on plant P uptake due to the low mobility of phosphate ions in soils. The commonly observed depletion of inorganic P concentration within the rhizosphere is, on the one hand, the result of the low mobility of P in soils, and on the other hand, enhances the diffusion of P by increasing the concentration gradient between root surface and the surrounding soil. The latter is particularly important for plants growing in soils with low P concentration (Lambers et al. 2008). Although we did not explicitly simulate rhizosphere processes here, the modelling framework will capture these processes when it is coupled to a model that resolves the rhizosphere (e.g. Goll et al 2018).

## 4.3 Applications to global land modelling

One of the key objectives for this study was to develop a model of inorganic P dynamics for global land models. The model of inorganic P dynamics developed here is suitable for being implemented into global land models. By calibrating the model using observations, we identified significant variations of all model parameters with soil physical and chemical properties, which have so far been ignored in all global land models. The empirical relationships developed in this study (see Table 2) are useful for deriving spatially varying model parameter values for global modelling. It is highly desirable to assess the robustness of those empirical relationships and the inorganic P model using additional field observations.

However, one of the key variables in the empirical relationships is the oxalate extractable metal oxide concentration ($O_x$), for which large scale data sets are missing from the available global datasets (e.g. Shangguan et al. 2014). Given the importance of this variable on inorganic P dynamics in soil, a concerted effort is needed to develop a global database for oxalate extractable metal oxide concentration. An alternative would be to develop an empirical function to relate oxalate extractable metal oxide concentration with some other
soil variables that are available spatially at a global scale, such as the *GlobalSoilMap* initiative (see Arrouays et al. 2014).

Biological activities can play significant role in soil phosphorus dynamics through immobilization, mineralization (both biological and biochemical) and uptake. When the inorganic P model developed is coupled to a global land model, such as CABLE (Wang et al. 2011) or ORCHIDEE (Goll et al. 2012), effects of biological activities on soil P dynamics will be accounted for.

Recently, He et al. (2021) published a global data set of soil phosphorus, which significantly improved the previous estimates of total soil P by Yang et al. (2013). With the improved global dataset and implementation of soil inorganic P model into global land models, a significant advance can be made in modelling of terrestrial biogeochemical cycles globally (see Wang and Goll 2021).

## 5. Conclusions

Our model of inorganic P transformations provides the basis for a data-driven quantification of soil P availability over various timescales which is needed for the assessment of P limitation on food & fiber production, carbon cycling, and climate change. By accounting of inorganic P transformations, it goes well beyond the current approaches that do not explicitly simulate the dynamics of inorganic and organic P pools with arbitrarily chosen parameter values.

The model, which is calibrated to observed soil P fractions and isotopic exchange kinetics, produces predictions in line with our theoretical understanding of soil P transformations. Model parameters vary significantly among soils, which highlights the need for accounting for the effects of soil chemical properties on inorganic soil P dynamics in global P models. The use of empirical formulation relating model parameters to soil properties seems promising to do so, but some of the needed information (i.e. oxalate extractable metal oxide concentration) is currently missing.

## Acknowledgment

Contributions from YPW is partially funded by the National Environmental Science Program funded by the Australian government. DG is supported by the ANR CLAND Convergence
Institute and the European Research Council Synergy project SyG-2013-610028
IMBALANCE-P. We also acknowledge support from a National Science Foundation
Research Coordination Grant (INCyTE; DEB-1754126) to investigate nutrient cycling in
terrestrial ecosystems. Data of IEK kinetics and soil P fractions are available from
Helfenstein et al. (2020) and Achat et al. (2016b).

## References

Achat, D.L., Augusto, L., Morel, C. & Bakker, M.R. (2011). Predicting available phosphate
ions from physical-chemical soil properties in acidic sandy soils under pine forests. *Journal of
Soils Sediments*, 11, 452-466, https://doi.org/10.1007/s11368-010-02329-9.

Achat, D.L., Augusto, L., Gallet-Budynek, A. & Loustau, D. (2016a). Future challenges in coupled C-N-P cycle models for terrestrial ecosystems under global change: a review. *Biogeochemistry*, 131, 173-202. https://doi.org/10.1007/s10533-016-0274-9.

Achat, D. L., Pousse, N., Nicolas, M., Brédoire, F., & Augusto, L. (2016b). Soil properties controlling inorganic phosphorus availability: general results from a national forest network and a global compilation of the literature. *Biogeochemistry*, 127, 255-272, https://doi.org/10.1007/s10533-015-0178-0.

Arrouays, D., Grundy, M. G., Hartemink, A. E., Hempel, J. W., Heuvelink, G. B et al. (2014). GlobalSoilMap: Toward a fine-resolution global grid of soil properties. *Advances in Agronomy*, 125, 93-134.

Arora, V. K., Katavouta, A., Williams, R.G., et al. (2020). Carbon0concentration and carbon-climate feedbacks in CMIP6 models and their comparison to CMIP5 models. *Biogeosciences* 17, 4173-4222.

Augusto, L., Achat, D., Jonard, M., Vidal, D. and Ringeval, B. (2017). Soil parent material – a major driver of plant nutrient limitations in terrestrial ecosystems. *Global Change Biology*, doi:10.1111/gcb.13691.

Barber, S.A. (1995). *Soil nutrient bioavailability: a mechanistic approach*. 2nd Edition, John Wiley, New York., 414p.
Barros, N.F., Comerford, N.B. & Barros, N.F. (2005). Phosphorus sorption, desorption, and resorption by soils of the Brazilian Cerrado supporting eucalypt. *Biomass and Bioenergy*, 28, 229-236.

Barrow, N.J. (1979). The description of desorption of phosphate from soil. *European Journal of Soil Science*, 30,259-270, http://doi.org/10.1111/j.1365-2389.1979.tb00983.x

Barrow, N.J. (1983). A mechanistic model for describing the sorption and desorption of phosphate by soil. *Journal of Soil Science*, 34:733-750.

Barrow, N.J. (1984). Modelling effects of pH on phosphate sorption by soils. *Journal of Soil Sciences*, 35, 283-297.

Barrow, N.J. (1999). The four laws of soil chemistry: the Leeper lecture 1998. *Australian Journal of Soil Sciences*, 37, 787-829.

Barrow, N.J. (2008). The description of sorption curves. *European Journal of Soil Science*, 59, 900-910.

Bertrand, I., Hollway, R.E., Armstrong, R.D. & McLaughlin, M.J. (2003). Chemical characteristics of phosphorus in alkaline soils from southern Australia. *Australian Journal of Soil Research*, 41, 61-76.

Borda, T., Celi, L., Bunemann, E.K., Oberson, A., Frossard, E. & Barberis, E. (2014). Fertilization strategies affect phosphorus forms and release from soils and suspended solids. *Journal Environmental Quality*, 43, 1024-1031, https://doi.org/10.2134/jeq2013.11.0436.

Borggaard, O.K. (1986). Iron oxides in relation to phosphate adsorption by soils. *Acta Agriculturae Scandinavica*, 36, 107-118.

Borggaard, O.K., Szilas, C., Gimsing, A.L. & Rasmussen, L.H. (2004). Estimation of soil phosphate adsorption capacity by means of a pedotransfer function. *Geoderma*, 118,55-61.

Buehler, S., Oberson, A., Rao, I.M., Frisen, D.K. & Frossard, E. (2002). Sequential phosphorus extraction of a 33P-labeled oxisol under contrasting agricultural systems. *Soil Science Society of America Journal*, 66, 868-877.

Buehler, S., Oberson, A., Sinaj, S., Friesen, D.K. & Frossard, E. (2003). Isotope methods for assessing plant available phosphorus in acid tropical soils. *European Journal of Soil Sciences*, 54, 605-616.

Accepted Article
Chen, J.S., Mansell, R.S., Nkedi-Kizza, P., & Burgoa, B.A. (1996) Phosphorus transport during transient, unsaturated water flow in an acid sandy soil. Soil. Sci. Am. J., 60:42-48.

Chen, C.R., Condron, L.M., Sinaj, S., Davis, M.R., Sherlock, R.R. & Frossard, E. (2003). Effects of plant species on phosphorus availability in a range of grassland soils. *Plant and Soil*, 256, 115-130.

Condron, L.M. and Newman, S. (2011) Revisiting the fundamentals of phosphorus
fractionation of sediments and soils. *Journal of Soils and Sediments*, 11, 830-840,
https://doi.org/10.01007/s11368-011-0363-2.

Cross, A.F. & Schlesinger, W.H. (1995). A literature review and evaluation of the Hedley
fractionation: Applications to the biogeochemical cycle of soil phosphorus in natural
ecosystems. *Geoderma*, 64, 197-214.

Davies, J.A.C., Tipping, E., Rowe, E.C., Boyle, J.F., Pannatier, E.G. & Martinsen, V. (2016). Long-term P weathering and recent N deposition control contemporary plant-soil C, N, and P. *Global Biogeochemical Cycles*, 30, 231-249, https://doi.org/10.1002/2015GB005167.

Duan, Q. Y., Gupta, V. K., & Sorooshian, S. (1993). Shuffled complex evolution approach for effective and efficient global minimization. *Journal of Optimization Theory and Applications*, 76, 501-521. https://doi.org/10.1007/bf00939380

Duan Q. (2003) Global optimization for watershed model calibration. In *Advances in Calibration of Watershed Models, Water Science and Application*, Series 6, American Geophysical Union, Washington p 89-104. Edited by Q Duan, H. Gupta, S. Sorooshian, A. Rousseau and R. Turcotte.

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

Fardeau, J.C., Guiraud, G. & Morel, C. (1996). The role of isotopic techniques on the evaluation of the agronomic effectiveness of P fertilizers. *Fertilizer Research*, 45, 101-109

Fernandez-Martinez, M., Vicca, S., Janssens, I.A. et al (2014). Nutrient availability as the key regulator of global forest carbon balance. *Nature Climate Change*, 4, 471-476.

Fleischer, K., Ramming, A., de Kauwe M.G., et al. (2019). Amazon forest response to CO2 fertilization dependent on plant phosphorus acquisition. *Nature Geoscience*, 12, 736-741.
https://doi.org/10.1038/s41561-019-0404-9.

Accepted Article
Frossard, E., Feller, C., Tiessen, H., Stewart, J.W.B., Fardeau, J. C. & Morel, J. L. (1993).

Can an isotopic method allow for the determination of the phosphate-fixing capacity of soils?

*Communication in Soil Science and Plant Analysis*, 24, 367–377.

Frossard, E., Condron, L.M., Oberson, A., et al (2000). Processes governing phosphorus availability in temperate soils. *Journal of Environmental Quality*, 29, 15–23.

Gimsing, A.L. & Borggaard, O.K. (2002). Competitive adsorption and desorption of glyphosate and phosphate on clay silicates and oxides. *Clay Mineralogy*, 37, 509–515.

Goll, D., Brovkin, V., Parida, B., Reick, C.H., Kattge, J., Reich, P., Van Bodegom, P. & Niinemets, Ü. (2012). Nutrient limitation reduces land carbon uptake in simulations with a model of combined carbon, nitrogen and phosphorus cycling, *Biogeosciences*, 9(9), 3547–3569.

Guedes, R.S., Melo, L.C.A., Vergutz, L., Rodrigues-Vila, A., Covelo, E.F. & Gernandes, A.R. (2016). Adsorption and desorption kinetics and phosphorus hysteresis in highly weathered soil by stirred flow chamber experiments. *Soil Tillage Research*, 162, 46-54.

Hamon, R.E., Bertrand, I. & McLaughlin, M.J. (2002). Use and abuse of isotopic exchange data in soil chemistry. *Australian Journal of Soil Sciences*, 40, 1371-1381.

Hartmann, J., Moosdorf, N., Lauerwald, R., Hinderer, M., and West A.J. (2014). Global chemical weathering and associated P-release- The role of lithology, temperature and soil properties. *Chemical Geology*, 363, 145-163.

http://dx.doi.org/10.1016/j.chemgeo.2013.10.025/

He, X.J., Augusto, L., Goll, D.S., Ringeval, B., Wang, Y.P., Helfenstein, J., Huang, Y.Y., Yu, K.L., Wang, Z.Q., Yang, Y.C. and Hou, E.Q. (2021) Global patterns and drivers of soil total phosphorus concentration. *Earth System Science Data Discussions*, https://doi.org/10.5194/essd-2021-166.

Hedley, M.J., Stewart, J.W.B. & Chauhan, B.S. (1982). Changes in inorganic and organic phosphorus fractions induced by cultivation practices and by laboratory incubation. *Soil Science Society of America Journal*, 46, 970–976

Helfenstein, J., Jegminat, J., McLaren, T.I. & Frossard, E. (2018a). Soil solution phosphorus turnover: derivation, interpretation, and insights from a global compilation of isotope exchange kinetics studies. *Biogeosciences*, 15, 105-114, https://doi.org/10.5194/bg-15-105-2018.
Helfenstein, J., Tamburini, F., von Sperber, C., Massey, M., Pistocchi, C., Chadwick, O., Vitousek, P., Kretzschmar, R., and Frossard, E. (2018b). Combining spectroscopic and isotopic techniques gives a dynamic view of phosphorus cycling in soil. *Nature Communications*, 9, 3226, https://doi.org/10.1038/s41467-018-05731-2.

Helfenstein, J., Pistocchi, C., Obrson, A., Tamburini, F., Goll, D.S. & Frossard, E. (2020). Estimates of mean residence times of phosphorus in commonly considered inorganic soil phosphorus pools. *Biogeosciences*, 17, 441-454, https://doi.org/10.5194/bg-17-441-2020.

Hinsinger, P. (2001). Bioavailability of soil inorganic P in the rhizosphere as affected by root-induced chemical changes: a review. *Plant and Soil*, 237, 173-195.

Hou, E., Luo, Y., Kuang, Y., Chen, C., Lu, X. et al. (2020). Global meta-analysis shows pervasive phosphorus limitation of aboveground plant production in natural terrestrial ecosystems. *Nature communications*, 11(1), 1-9.

Hou, E., Tan, X., Heenan, M. & Wen, D. (2018). A global dataset of plant available and unavailable phosphorus in natural soils derived by Hedley method. *Scientific Data*, 5, 180166.

Huang, W.J., Zhou, G.Y., Liu, J.X., Duan, H.L., Liu, X.Z., Fang, X. & Zhang, D.Q. (2014). Shifts in soil phosphorus fractions under elevated CO2 and N addition in model forest ecosystems in subtropical China. *Plant Ecology*, 215, 1373-1384.

Jiang, J., Wang, Y-P., Liu, F.C., Du, Y., Zhuang, W., Chang, Z.B., Yu, M.X. & Yan, J.H. (2021). Antagonistic and additive interactions dominate the response of belowground carbon-cycling processes to nitrogen and phosphorus additions. *Soil Biology and Biochemistry*, 156, 108216, https://doi.org/10.1016/j.soilbio.2021.108216.

Johnson, A., Frizano, J., and Vann, D. (2003). Biogeochemical implications of labile phosphorus in forest soils determined by the Hedley fractionation procedure, *Oecologia*, 135, 487-499, https://doi.org/10.1007/s00442-002-1164-5.

Kvakic, M., Pellerin, S., Ciais, P., Achat, D.L., Augusto, L., Denoroy, P., et al. (2018). Quantifying the limitation to world cereal production due to soil phosphorus status. *Global Biogeochemical Cycles*, 32, 143-157, https://doi.org/10.1002/2017GB005754.

Lambers, H., Raven, J.A., Shaver, G.R. & Smith, S.E. (2008). Plant nutrient-acquisition strategies change with soil age. *Trends in Ecology & Evolution*, 23, 95-103.

Accepted Article
Lang, F., Kruger, J., Amelung, W. et al. (2017). Soil phosphorus supply controls P nutrition strategies of beech forest ecosystems in central Europe. *Biogeochemistry*, 136, 5-29, https://doi.org/10.1007/s10533-0375-0.

Lebauer, D. S., & Treseder, K. K. (2008). Nitrogen limitation of net primary productivity in terrestrial ecosystems is globally distributed. *Ecology*, 89, 371–379.

Lookman, R., Freese, D., Merckx, R., Vlassek, K., Van Riemsdijk, W.H. (1995). Long term kinetics of phosphate release from soil. *Environmental Science Technology*, 29, 1569-1575.

Lu, X.J., Wang, Y.P., Ziehn, T. & Dai, Y.J. (2013). An efficient method for global parameter sensitivity analysis and its applications to the Australian community land surface models (CABLE). *Agricultural and Forest Meteorology*, 182-183, 292-303.

Manzoni, S. & Porporato, A. (2009). Soil carbon and nitrogen mineralization: theory and model across scales. *Soil Biology and Biochemistry*, 41, 1355-1379.

McGechan, M.B. and Lewis, D.R. (2002). Sorption of phosphorus by soil, Part 1: Principles, equations and models. *Biosystems Engineering*, 82, 1-24. https://doi.org/10.1006/bioe.2002.0054.

McGill, W.B. & Cole, C.V. (1981). Comparative aspects of cycling of organic C, N, S and P through soil organic matter. *Geoderma*, 26, 267-286.

Mengel, K. & Kirkby, E.A. (2001). *Principles of plant nutrition*, Dordrecht, Kluwer Academic Publishers. Fifth edition, 849 pp.

Morel, C., Tunney, P., Plenet, D. & Pellerin, S. (2000). Transfer of phosphate ions between soil and solution: Perspectives in soil testing. *Journal of Environment Quality*, 29, 50–59.

Okajima, H., Kubota, H. & Sakuma, T. (1983). Hysteresis in the phosphorus sorption and desorption processes of soils. *Soil Science and Plant Nutrition*, 29, 3, 271-283, https://doi.org/10.1080/00380768.1983.10434628.

Parton, W. J., D. S. Schimel, C. V. Cole, & D. S. Ojima (1987). Analysis of factors controlling soil organic matter levels in Great Plains grassland. *Soil Science Society of America Journal*, 51, 1173–1179.

Penuelas, J., Poulter, B., Sardans, J., Ciais, P. et al. (2013). Human-induced nitrogen phosphorus imbalances alter nature and managed ecosystem across the globe. *Nature Communications*, 4, 2934, https://doi.org/10.1038/ncomms3934.
Regelink, I.C., Weng, L., Lair, G.J. & Comans, R.N.J. (2015). Adsorption of phosphate and organic matter on metal (hydr)oxides in arable and forest soil: a mechanistic modelling study. *European Journal of Soil Science*, 66, 867-875.

Ringeval, B., Augusto, L., Monod, H., Apeldoorn, D.V. et al. (2017). Phosphorus in agricultural soils: drivers of its distribution at the global scale. *Global Change Biology*, 23, 3418-3432, https://doi.org/10.1111/gcb.13618

Schubert, S., Steffens, D. and Ashraf, I. (2020). Is occluded phosphate plat-available? *Journal of Plant Nutrition and Soil Science*, 183:338-344, https://doi.org/10.1002/jpin.201900402.

Shang, C., Steward, J.W.B. & Huang, P.M. (1992). pH effect on kinetics of adsorption of organic and inorganic phosphates by short-range ordered aluminum and iron precipitates. *Geoderma*, 53, 1-14

Shayan, A. & Davey, B.G. (1978). A universal dimensionless phosphate adsorption isotherm for soil. *Soil Science Society Journal of America Journal*, 42, 878-882.

Sun, Y., Peng, S.S., Goll, D.S., Ciais, P. et al. (2017). Diagnosing phosphorus limitations in natural terrestrial ecosystems in carbon cycle models. *Earth's Future*, 5, 730-749, https://doi.org/10.1002/2016EF000472.

Shangguan, W., Dai, Y., Duan, Q., Liu, B., & Yuan, H. (2014). A global soil data set for earth system modeling. *Journal of Advances in Modeling Earth Systems*, 6(1), 249-263.

Tiessen, H., Moir, J. (1993). Characterization of available P by sequential extraction. In *Soil Sampling and Methods of Analysis* (ed. Carter, M.R), 75-86 (Lewis Publishers, 1993).

Wang, Y.P., Kowalczyk, E., Leuning, R., Abramowitz, G., Raupach, M.R., Pak, B., Gorsel, E.V. and Luhar, A. (2011) Diagnosing errors in a land surface model (CABLE) in the time and frequency domains. *Journal of Geophysical Research*, 116, G01034, doi:10.1029/2010JG001385.

Wang, Y.P. & Goll, D.S. (2021). Modelling of land nutrient cycles: recent progresses and future development. *F1000Research, Faculty Review*, 10:(53) https://doi.org/10.12703/r/10-53.

Wang, Y.P., Zhang, Q., Pitman, A.J. & Dai, Y.J. (2015). Nitrogen and phosphorus limitation reduces the effects of land use change on net atmosphere-land CO2 exchange. *Environmental Research Letters*, 10, 014001, https://doi.org/10.1088/1748-9326/10/014001.
Wang, Y.P., Law, R.M. & Pak, B. (2010). A global model of carbon, nitrogen and phosphorus cycles for the terrestrial biosphere. *Biogeosciences*. 7, 2261-2282.

Wang, Y.P. and Goll, D.S. (2021) Modelling of land nutrient cycles: recent progresses and future development. *Faculty Review*, 10:(53) https://doil.org/10.12703/r/10-53.

Wang, Y.P., Zhang, H.C., Ciais, P., Goll, D., Huang Y.Y., Qood, J.D., Ollinger, S.V., Tang, X.L. and Prescher, A.K. (2021). Microbial activity and root carbon inputs are more important than soil carbon diffusion in simulating soil carbon profiles. *Journal of Geophysical Research: Biogeosciences*, 126, e2020JG006205, https://doi.org/10.1029/2020JG006205.

Weng, L.P., Vega, F.A. & Riemsdijk, W.H.V. (2011). Comparative and synergistic effects in pH dependent phosphate adsorption in soils: LCD modelling. *Environmental Science and Technology*, 45,8420-8428, dx.doi.org/10.1021/es201844d.

Yang, X., W. M. Post, P. E. Thornton, and A. Jain (2013). The distribution of soil phosphorus for global biogeochemical modeling, *Biogeosciences*, 10(4), 2525-2537.

Yang, X., Thornton, P.E., Ricciuto, D.M. and Post, W.M. (2014) The role of phosphorus dynamics in tropical forests – a modeling study using CLM-CNP. *Biogeosciences*, 11:1667-1681.

Zhang, Q., Wang, Y.P., Matear, R.J., Pitman, A.J. & Dai, Y.J. (2014). Nitrogen and phosphorous limitations significantly reduce future allowable CO2 emissions. *Geophysical Research Letters*, 41, 632-637, https://doi.org/1002/2013GL058352.

Zhu, Q., Riley, W.J., Tang, J., Collier, N., Hoffman, F.M. et al. (2019). Representing nitrogen, phosphorus, and carbon interactions in the E3SM land model: development and global benchmarking. *Journal of Advances in Modeling Earth Systems*, 11, 2238-2258, https://doi.org/10.1029/2018MS001571.
Appendix 1: A model of the isotopically labelled P dynamics

Assuming that soil inorganic P among the four pools are at steady state, and that a minute amount of P was added to the soil solution at time t=0, the dynamics of the added P can be described by the following equations:

$$ \frac{d(P_W + P_W^*)}{dt} = k_{WL}(P_L + P_L^*) + k_{WS}(P_S + P_S^*) - (k_{LW} + k_{SW})(C_\infty + C_W^*)^b \quad (A1) $$

$$ \frac{d(P_L + P_L^*)}{dt} = k_{LW}(C_\infty + C_W^*)^b - k_{WL}(P_L + P_L^*) \quad (A2) $$

$$ \frac{d(P_S + P_S^*)}{dt} = k_{SW}(C_\infty + C_W^*)^b - k_{WS}(P_S + P_S^*) + k_{SO}(P_O + P_O^*) - k_{OS}(P_S + P_S^*) \quad (A3) $$

$$ \frac{d(P_O + P_O^*)}{dt} = k_{OS}(P_S + P_S^*) - k_{SO}(P_O + P_O^*) \quad (A4) $$

With the conditions of

$$ P_W^* = R, P_L^* = P_S^* = P_O^* = 0 \text{ at } t=0. $$

where $P_W^*, P_L^*, P_S^*$ and $P_O^*$ are the isotopically labelled P in solution, labile, sorbed and occluded P pools in mg P/kg, and R is the amount of isotopically labelled P added to soil solution at t=0, $C_W^* = \frac{P_W^*}{v}$, and $C_\infty = \frac{P_W(t \to \infty)}{v}$. After addition, the labelled P in soil solution decrease over time, as a result of the exchanges between solution P and other inorganic P pools (labile, sorbed and occluded P) in soil in the absence of biological activities.

Substituting eqn 7-10 into eqn (A1) to (A4), and make use of the approximation of $(C_\infty + C_W^*)^b \approx C_\infty^b + bC_\infty^{b-1}C_W^*$, we have

$$ \frac{dP_W^*}{dt} = k_{WL}P_L^* + k_{WS}P_S^* - (k_{LW} + k_{SW})(bC_\infty^{b-1}C_W^*) \quad (A5) $$

$$ \frac{dP_L^*}{dt} = k_{LW}(bC_\infty^{b-1}C_W^*) - k_{WL}P_L^* \quad (A6) $$

$$ \frac{dP_S^*}{dt} = k_{SW}(bC_\infty^{b-1}C_W^*) - k_{WS}P_S^* - k_{OS}P_S^* + k_{SO}P_O^* \quad (A7) $$

$$ \frac{dP_O^*}{dt} = k_{OS}P_S^* \quad (A8) $$

By normalizing all the pool size by the amount of isotopically labelled P added to the solution at t=0, R, we have

$$ r_W = \frac{P_W^*}{R}, r_L = \frac{P_L^*}{R}, r_S = \frac{P_S^*}{R}, r_O = \frac{P_O^*}{R}, \rho_W = \frac{f_W}{v}, \rho_\infty = \frac{C_\infty^*}{R}. $$

This article is protected by copyright. All rights reserved.
We can rewrite eqn A5-A8 as

$$ \frac{dr_W}{dt} = k_{WL}r_L + k_{WS}r_S - (k_{LW} + k_{SW})(bC_\infty^{b-1}\rho_W) \quad (A9) $$

$$ \frac{dr_L}{dt} = k_{LW}(bC_\infty^{b-1}\rho_W) - k_{WL}r_L \quad (A10) $$

$$ \frac{dr_S}{dt} = k_{SW}(bC_\infty^{b-1}\rho_W) - k_{WS}r_S - k_{OS}r_S + k_{SO}r_O \quad (A11) $$

$$ \frac{dr_O}{dt} = k_{OS}r_S - k_{SO}r_O \quad (A12) $$

with initial conditions of $r_W = 1$, and $r_L = r_S = r_O = 0$ at $t=0$.

and $r_W = \frac{r(t)}{R}$ at any time $t$, where $r(t)$ is the isotopically labelled $^{33}$P in solution.

At steady state, we also have $r_W(t \to \infty) = f_W = \frac{r(\infty)}{R}$, $r_L(t \to \infty) = f_L$, $r_S(t \to \infty) = f_S$, $r_O(t \to \infty) = f_O$. Setting the left-hand sides of eqns A9- A12 and substituting all steady-state pool sizes into those four equations, we have:

$$ \frac{k_{LW}}{k_{WL}} = \frac{f_L}{bC_\infty^{b-1}f_W/v}, \quad (A13) $$

$$ \frac{k_{SW}}{k_{WS}} = \frac{f_S}{bC_\infty^{b-1}f_W/v}, \quad (A14) $$

$$ \frac{k_{SO}}{k_{OS}} = \frac{f_S}{f_O}, \quad (A15) $$

Therefore parameter the ratios of $\frac{k_{LW}}{k_{WL}}$, $\frac{k_{SW}}{k_{WS}}$ and $\frac{k_{SO}}{k_{OS}}$ are known if the fractions of different pools at steady state, and initial concentration of solution P at steady state before ($C_\infty$) $t=0$ ($C_\infty$) are given. Because short-term (< a few days) IEK measurements are likely to constrain $k_{SO}$ or $k_{OS}$ well, we fixed $k_{SO}$ at 0.02 year$^{-1}$, and calculated $k_{OS}$ using eqn A15. As a result, we only need to estimate the values of $k_{LW}$, $k_{SW}$ and $b$ from the results of IEK experiments.
<figure><img src="image_3.png" /><figcaption>Figure 1. A schematic diagram showing structure of the inorganic soil P model and exchanges between solution P and labile, sorbed or occluded P. $F_{WL}$ represents a flux from labile (L) pool to solution pool (W), subscript S and O represent sorbed and occluded P pools, respectively. and There is no exchange between solution P and primary mineral P or organic P.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2. Results of model calibration. (a) comparison of the root mean square error (RMSE) of the optimized linear and nonlinear models for the 147 sites; (b) frequency distribution of the RMSE of the nonlinear model for the 147 sites; (c) frequency distribution of the squared correlation coefficient of the linear regression between the predicted ($r_W$) and calculated values of $r/R$ using the soil-specific parameters ($m, n, r(\infty)/R$ and $C_\infty$) at different times and between the predicted and measured fractions of different pools at steady state and (d) variations of mean values of IEK parameters $m$ and $n$ for each class of RMSE of the nonlinear model.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3. Comparison of the variation of isotopically labelled $^{33}$P radioactivity in soil solution ($r/R$) with time as calculated using the eqn 11 with the soil-specific parameters ($m, n, r(\infty)/R$ and $C_\infty$) (dark grey curve) with that predicted by the inorganic P model (light grey curve) for soil sample 1 (a) or 19 (b); and the modelled simulated fractions of added P isotopes in different inorganic P pools as a function of time (c) for soil 1 and d for soil 19. Parameter $n=0.437$ and sand:silt:clay =40:40:20 for soil 1; $n=0.12$ and sand:silt:clay=84:11:5 for soil 19.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4. Variations of the optimized values of parameters in the inorganic P model with IEK parameter *n* for 147 soils. *k*<sub>WL</sub> is the desorption rate constant in day<sup>-1</sup> for labile, and *k*<sub>SW</sub> is the sorption coefficient in mg P (kg soil)<sup>-1</sup> day<sup>-1</sup> (mg P/L)<sup>-b</sup> for sorbed P, respectively, and *k*<sub>OS</sub> is the sorption rate constant in day<sup>-1</sup> for the occluded P, *b* is the exponent in the Freundlich equation for sorption, *f*<sub>S</sub> and *f*<sub>O</sub> are the fractions of sorbed and occluded P, respectively. The best fitted regression equations are: *k*<sub>LW</sub>=12.24exp(-3.48*n*), *r*<sup>2</sup>=0.17*; *k*<sub>WS</sub>/*f*<sub>S</sub>=8.37exp(6.57*n*), *r*<sup>2</sup>=0.40*; *k*<sub>OW</sub>/*f*<sub>O</sub>=2.29x10<sup>-4</sup>-2.0x10<sup>-4</sup>*n*, *r*<sup>2</sup>=0.1*; and *b*=1.09-1.22*n*, *r*<sup>2</sup>=0.14*.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5. (a) Variation of solution P concentration at steady state ($C_\infty$ in mg/L) with $C_x$, or the ratio of total soil carbon concentration and oxalate extracted Al and Fe concentration in g C/mmol. The best fitted equation: $y=10^{-[-1.7425+1.7782(1-\exp(-1.6653C_x))]}$, $r^2=0.54$; (b) comparison of the predicted $C_\infty$ by the best-fitted regression and the observed values. The solid line represents 1:1 line.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. (a) mean (bar height) and one standard deviation (error bar) of the relative sensitivity of the soil P bioavailability over 7 days (black bar) or 180 days (grey bar) to oxalate extractable metal oxide concentration ($O_x$), mean annual temperature ($T_a$), sand ($s_s$), clay ($s_c$) or silt ($s_i$) percentage, total soil C (C), total inorganic P ($P_{in}$), fractions of solution P ($f_W^*$), labile P ($f_L^*$) or sorbed P ($f_S^*$); (b) importance of different soil properties. Properties as listed in Table S1 that have little important (<10⁻³) or low sensitivity (<10⁻³) are not included here.</figcaption></figure>
<table><caption>Table 1. A list of the symbols used in this study</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Unit</th></tr></thead><tbody><tr><td colspan="3"><strong>Variables</strong></td></tr><tr><td>t</td><td>time</td><td>minute or day or year</td></tr><tr><td>v</td><td>Soil water content</td><td>L/kg</td></tr><tr><td>r</td><td>Radioactivity of soil solution</td><td>Bq</td></tr><tr><td>R</td><td>Radioactivity of added isotopically labelled P added at t=0</td><td>Bq</td></tr><tr><td>P<sub>in</sub></td><td>Total inorganic P excluding primary inorganic P</td><td>mg P/kg soil</td></tr><tr><td>P<sub>W</sub>, P<sub>L</sub>, P<sub>S</sub>, P<sub>O</sub></td><td>Amount of inorganic P in soil solution, labile, sorbed and occluded pool, respectively.</td><td>mg P/kg soil</td></tr><tr><td>P<sub>W</sub><sup>*</sup>, P<sub>L</sub><sup>*</sup>, P<sub>S</sub><sup>*</sup>, P<sub>O</sub><sup>*</sup></td><td>Amount of added isotopically labelled P in soil solution, labile, sorbed and occluded pool, respectively.</td><td>mg P/kg soil</td></tr><tr><td>r<sub>W</sub>, r<sub>L</sub>, r<sub>S</sub>, r<sub>O</sub></td><td>Fraction of isotopically labelled P in soil water, labile, sorbed and occluded pool, respectively.</td><td>—</td></tr><tr><td>C<sub>W</sub> and C<sub>∞</sub></td><td>Solution P concentration (C<sub>W</sub>) and its steady state value (C<sub>∞</sub>)</td><td>mg P (kg soil)<sup>-1</sup> L<sup>-1</sup></td></tr><tr><td>ρ<sub>W</sub> and ρ<sub>∞</sub></td><td>Fraction of the added isotopically labelled P remaining in solution per unit volume of soil water (ρ<sub>W</sub>) and its steady state value (ρ<sub>∞</sub>)</td><td>L<sup>-1</sup></td></tr><tr><td>C<sub>W</sub> and C<sub>∞</sub></td><td>Solution P concentration (C<sub>W</sub>) and its steady state value (C<sub>∞</sub>)</td><td>mg P (kg soil)<sup>-1</sup> L<sup>-1</sup></td></tr><tr><td>ρ<sub>W</sub> and ρ<sub>∞</sub></td><td>Fraction of the added isotopically labelled P remaining in solution per unit volume of soil water (ρ<sub>W</sub>) and its steady state value (ρ<sub>∞</sub>)</td><td>L<sup>-1</sup></td></tr><tr><td colspan="3"><strong>Parameters</strong></td></tr><tr><td>m, n</td><td>empirical parameters</td><td>—</td></tr><tr><td>k<sub>WL</sub></td><td>Desorption rate constant from labile pool</td><td>day<sup>-1</sup></td></tr><tr><td>k<sub>LW</sub></td><td>Sorption coefficient for labile pool</td><td>mg P (kg soil)<sup>-1</sup> day<sup>-1</sup> (mg P/L)<sup>-b</sup></td></tr><tr><td>k<sub>WS</sub></td><td>Desorption rate constant from labile pool</td><td>day<sup>-1</sup></td></tr><tr><td>k<sub>SW</sub></td><td>Sorption coefficient for sorbed pool</td><td>mg P (kg soil)<sup>-1</sup> day<sup>-1</sup> (mg P/L)<sup>-b</sup></td></tr><tr><td>k<sub>OS</sub>/k<sub>OS</sub></td><td>Sorption/desorption rate constant for occluded pool</td><td>day<sup>-1</sup></td></tr><tr><td>b</td><td>An empirical constant</td><td>—</td></tr><tr><td>J</td><td>Total cost</td><td>—</td></tr><tr><td>f<sub>W</sub>, f<sub>L</sub>, f<sub>S</sub>, f<sub>O</sub></td><td>The modelled fractions of inorganic P in soil solution, labile, sorbed or occluded pool at steady state, respectively</td><td>—</td></tr><tr><td>f<sub>W</sub><sup>*</sup>, f<sub>L</sub><sup>*</sup>, f<sub>S</sub><sup>*</sup>, f<sub>O</sub><sup>*</sup></td><td>The measured fractions of inorganic P in soil solution, labile, sorbed or occluded pool at steady state, respectively.</td><td>—</td></tr></tbody></table>
Table 2. The best fitted nonlinear regressions for solution P concentration at steady state ($C_∞$, in mg P/L), parameter *b* in the Freundlich equation for P sorption, desorption rate constant of labile P ($k_{WL}$ in day⁻¹), sorption coefficient of sorbed P ($k_{SW}$ in mg P (kg soil)⁻¹ day⁻¹ (mg P/L)⁻⁶), and sorption rate constant for the occluded P ($k_{OS}$ in day⁻¹). $T_a$ is mean annual surface air temperature in °C, $O_x$ is oxalate metal oxide (Al and Fe) concentration in mmol/kg soil, and *C* is soil carbon concentration in g C per kg of soil, $C_x$ is the ratio of total soil C and oxalate metal oxide concentrations in g/mmol ($C_x=0.001C/O_x$), $f_W^*$, $f_L^*$, $f_S^*$ and $f_O^*$, are the fractions of solution P, labile P, sorbed P and occluded P in soil, respectively; and $f_W^*+f_L^*+f_S^*+f_O^*=1$. $P_s$ is the amount of sorbed P in soil (mg P/kg soil), $P_{in}$ is the total amount of inorganic P excluding primary mineral P in mg P/kg soil, $P_o$ is the organic P in mg P/kg soil, $p_H$ is soil pH measured in water, $s_s$, $s_c$ and $s_i$ are sand, clay and silt percentages, respectively. All correlations are significant (p<0.05).

<table><thead><tr><td>Equation</td><td>r<sup>2</sup></td></tr></thead><tbody><tr><td>C<sub>∞</sub>=10<sup>[-1.7425+1.7782(1-exp(-1.6653C<sub>x</sub>))]</sup></td><td>0.54</td></tr><tr><td>b = exp(0.532 - 0.130O<sub>x</sub> - 37.49f<sub>W</sub><sup>*</sup> + 0.81f<sub>L</sub><sup>*</sup> + 0.003s<sub>s</sub>)</td><td>0.43</td></tr><tr><td>k<sub>WL</sub> = -4.82 + 209f<sub>W</sub><sup>*</sup> + 14.64f<sub>O</sub><sup>*</sup> + 9.26f<sub>S</sub><sup>*</sup> - 0.008C - 0.0003P<sub>in</sub> - 0.018s<sub>i</sub></td><td>0.28</td></tr><tr><td>k<sub>SW</sub> = exp(0.002 + 4.0f<sub>S</sub><sup>*</sup> + 0.0008P<sub>o</sub> + 0.012C + 0.108T<sub>a</sub> - 0.0002P<sub>in</sub>)</td><td>0.64</td></tr><tr><td>k<sub>OS</sub>=1.6×10<sup>-5</sup>+0.0001f<sub>O</sub><sup>*</sup>-0.0001f<sub>S</sub><sup>*</sup>+1.3×10<sup>-5</sup> p<sub>H</sub>-9.0×10<sup>-6</sup> O<sub>x</sub>-4.48×10<sup>-7</sup>s<sub>c</sub>+1.1×10<sup>-7</sup>C</td><td>0.90</td></tr></tbody></table>
<table><caption>Table 3. Inorganic P fractions, total inorganic P, parameter b in Freundlich equation, desorption rate constant for labile P ($k_{WL}$), or sorbed P ($k_{WS}$) for the top 0-20cm soils at the four sites in Carigmague, Columbia. SAV, GL, CR and RGM represent native savannah, grass-legume, continuous rice or rice green manure rotation, respectively. Data of P fractions and total amount of inorganic P were taken from Buehler et al. (2002). $k_{WL}$ and $k_{WS}$ were estimated from model calibration. The measured plant available was taken from Table 6 of Buhler et al. 2003) (third harvest at 12 weeks). In calculating plant available P, $x$=0.0001 and $t$=84 days for all four sites. Error estimate of plant available P based on observations were not provided by Buehler et al. (2002).</caption><thead><tr><th></th><th>SAV</th><th>GL</th><th>CR</th><th>RGM</th></tr></thead><tbody><tr><td>Fractions (labile/sorbed/occluded)</td><td>0.03/0.33/0.63</td><td>0.06/0.32/0.62</td><td>0.08/0.60/0.32</td><td>0.157/0.47/0.38</td></tr><tr><td>Total inorganic P (mg P/kg)</td><td>66</td><td>85</td><td>171</td><td>213</td></tr><tr><td>b</td><td>1.252</td><td>0.29</td><td>0.20</td><td>0.77</td></tr><tr><td>k<sub>WL</sub> (1/day)</td><td>0.0195</td><td>6.23</td><td>4.93</td><td>3.97</td></tr><tr><td>k<sub>WS</sub> (1/day)</td><td>0.0132</td><td>0.47</td><td>0.12</td><td>0.10</td></tr><tr><td>Plant available P (observed) (mg P/kg soil)</td><td>23</td><td>29</td><td>91</td><td>134</td></tr><tr><td>Our estimate (mg P/kg)</td><td>25</td><td>30</td><td>97</td><td>131</td></tr></tbody></table>