Toward a Global Model for Soil Inorganic Phosphorus Dynamics: Dependence of Exchange Kinetics and Soil Bioavailability on Soil Physicochemical Properties.

The representation of phosphorus (P) 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 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 2 > 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 modeled 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.

Introduction

Nitrogen (N) and phosphorus (P) are soil nutrients limiting primary production and other processes of terrestrial ecosystems (e.g., Lebauer & 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 CO2 emissions (Wang et al., 2015), and carbon‐climate feedback (Arora et al., 2020; Zhang et al., 2014). As a result, nutrient cycles have been implemented into an increasing number of global land models (see Achat, Augusto, et al., 2016). While several land models (>10) have incorporated the N cycle (see Manzoni & 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, Augusto, et al., 2016; 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 & 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 & 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%–60%) of total soil P with only a very small fraction (<1%) in soil solution (Mengel & Kirkby, 2001; Morel et al., 2000). Although validity of chemical fractionation for soil phosphorus was debated (see Condron & Newman, 2011), the Hedley fractionation technique (Hedley et al., 1982; Tiessen & 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 & 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 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, Pousse, et al., 2016).

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, Pousse, et al., 2016; Helfenstein, Jegminat, et al., 2018). 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, Pousse, et al., 2016). 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 Goll et al., 2012; Wang et al., 2010; 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 min and 1 day) between solution P and labile P (see Wang et al., 2010 and Yang et al., 2014 for example).

Exchange of solution and labile or sorbed phosphorus can be measured in the laboratory using carrier‐free isotopically labelled P (H3 33PO4), or isotopic exchange kinetics (IEK) experiments (see Fardeau et al., 1996; Frossard et al., 1993). Recently, Helfenstein, Tamburini, et al. (2018) 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, Pousse, et al. (2016). The objectives of this study are to (a) develop and calibrate a model of inorganic P dynamics that is suitable for use in global land models; (b) quantify how model parameters vary with soil physical and chemical properties, and climate; (c) define soil bioavailable P (a P) mathematically, and quantify how a P varies with soil physical and chemical properties, and time scale.

Material and Methods

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

A schematic diagram showing structure of the inorganic soil P model and exchanges between solution P and labile, sorbed or occluded P. FWL 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.

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

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 & 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 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–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 & 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: and

Dynamics of the closed system of four inorganic soil P pools is given by with P  = vC .

The model of inorganic P dynamics consists of Equations 7–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 ∞ 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.

Table 1

A List of the Symbols Used in This Study

Symbol	Definition	Unit
Variabes
t	time	minute or day or year
v	Soil water content	L/kg
r	Radioactivity of soil solution	Bq
R	Radioactivity of added isotopically labeled P added at t = 0	Bq
Pin	Total inorganic P excluding primary inorganic P	mg P/kg soil
PW, PL, PS, PO	Amount of inorganic P in soil solution, labile, sorbed and occluded pool, respectively.	mg P/kg soil
PW∗,PL∗,PS∗,PO∗	Amount of added isotopically labeled P in soil solution, labile, sorbed and occluded pool, respectively.	mg P/kg soil
rW,rL,rS,rO	Fraction of isotopically labeled P in soil water, labile, sorbed and occluded pool, respectively.	___
C W and C ∞	Solution P concentration (C W ) and its steady state value (C ∞)	mg P (kg soil)−1 L−1
ρ W and ρ ∞	Fraction of the added isotopically labeled P remaining in solution per unit volume of soil water (ρ W ) and its steady state value (ρ ∞)	L−1
C W and C ∞	Solution P concentration (C W ) and its steady state value (C ∞)	mg P (kg soil)−1 L−1
ρ W and ρ ∞	Fraction of the added isotopically labeled P remaining in solution per unit volume of soil water (ρ W ) and its steady state value (ρ ∞)	L−1
Parameters
m, n	empirical parameters	___
kWL	Desorption rate constant from labile pool	day−1
kLW	Sorption coefficient for labile pool	mg P (kg soil)−1 day−1 (mg P/L)−b
kWS	Desorption rate constant from labile pool	day−1
kSW	Sorption coefficient for sorbed pool	mg P (kg soil)−1 day−1 (mg P/L)−b
kOS/kOS	Sorption/desorption rate constant for occluded pool	day−1
b	An empirical constant	___
J	Total cost	___
fW, fL, fS, fO	The modeled fractions of inorganic P in soil solution, labile, sorbed or occluded pool at steady state, respectively	___
fW∗,fL∗fS∗,fO∗	The measured fractions of inorganic P in soil solution, labile, sorbed or occluded pool at steady state, respectively.	___

Exchange is relatively fast between labile and solution P, and slow between sorbed and solution P, or k LW >> k SW, and k WL >> k WS, the model as represented by Equations (7), (8), (9), (10) can simulate both fast and slow processes governing soil inorganic P dynamics, as discussed by McGechan and Lewis (2002), and Barrow (2008).

If exchanges between sorbed and occluded P pools are ignored, the steady state solution to Equations (7), (8), (9) will give the following relationship: where are solution P concentration, labile P and sorbed P at steady state. Equation 11 is the Freundlich equation.

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 modeling 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 (32P or 33P) 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): 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 Equation 11 to the measured data. The term r(∞)/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, Equation 11 gives: r(t = 0)/R = 1 + r(∞)/R. To correct for this small discrepancy, we modified Equation 11 as follows to make the equation correct also for t = 0:

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

In the IEK data compiled (Helfenstein et al., 2020), parameters m and n were estimated by fitting Equation 11 rather than Equation 12 to the measurements. As r(∞)/R is usually less than 0.001, estimates of m and n should vary very little whether Equation 11 or Equation 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., 1993). In some studies, a value of r/R at t = 1 min after the addition is used for m. The value of parameter n usually varies between 0.1 and 0.7 (Achat, Pousse, et al., 2016). Desorption of inorganic P into solution is generally slow for soil with a small n value. This study did not fit 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 (Equations 13–16, see the below).

Modeling 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: (a) prior P addition, all inorganic pools are in a steady state; (b) 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: with initial condition of r  = 1, and r  = r  = r  = 0 at t = 0.

Where r , r ,r and r are modeled fractions of isotopically labeled 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 (t→∞) = f , r (t→∞) = f , r (t→∞) = f , and r (t→∞) = f , 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 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

As the uncertainty of the calculated the fraction of added P remaining in solution using Equation 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 years−1, a value used in CASA‐CNP (Wang et al., 2010) and calculated k OS using Equation 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 Equations (17), (18), (19).

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 and

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 et al., 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 and ΔP are the difference in the steady state pool sizes of labile and sorbed P pools, respectively, and are given by when b = 1, we have ΔP = (1−x)f P and ΔP = (1−x)f P .

The soil P bioavailability, a p, is calculated as

Because x ≤ 1, b > 0, then (1−x ) 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 Equation 11). That is equivalent to the soil P bioavailability calculated using Equation 25 with t→∞, x = 0.

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 modeled fraction of added P remaining in solution (r W) and the calculated value of r(t)/R using Equation 12 with the published values of m, n and r(∞)/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 and

Here J 1 is calculated as the sum of squared differences for the hourly values of the first 24 hr after labeling, and daily values of first 365 days and yearly values of the first 20 years f , f and f 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 Equations 13–16, respectively.  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 Equation 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.

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, Pousse, et al. (2016). For each soil sample, the data included estimates of r(∞)/R, m and n from the IEK experiments, together with measurements of C ∞ (in mg/L), total inorganic soil P (mg P/kg). For the 45 soil samples compiled by Helfenstein et al. (2020), measurements of the amounts of inorganic P extracted by resin (resin‐P; solution P), NaHCO3 (NaHCO3‐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, Tamburini, et al. (2018) 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 5,275 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, Tamburini, et al. (2018) 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, Pousse, et al., 2016), we have 147 soil samples representing a wide range of climate, soil physical and chemical properties (see Figure S1 in Supporting Information S1). Table S1 in Supporting Information S1 lists the ranges of variations for those variables or soil properties of the 147 soil samples.

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 a in oC or mean annual precipitation P a 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 a and P a), 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.

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 Equation 25 shows that, a P depends on the rate constants of desorption, parameter b, C ∞, 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 ∞ and k LW and k WS using Equations 17 and 18, then a P using Equation 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, β j, was calculated as where μ i and σ i are the mean and standard deviation of the elemental effect of T a or soil property on a P, and , 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).

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 (kWL in day𢈒1), Sorption Coefficient of Sorbed P (kSW in mg P (kg soil)𢈒1 day𢈒1 (mg P/L)−b), and Sorption Rate Constant for the Occluded P (kOS in day−1)

Equation	r2
C_∞=10^[−1.7425+1.7782(1−exp(−1.6653C__X))]	0.54
b=exp(0.532−0.13O_x−37.49fW*+0.81fL*+0.003s_s)	0.43
k_WL=−4.82+209fW*+14.64fO*+9.26fS*−0.008C−0.0003P_in−0.018s_i	0.28
k_SW=exp(0.002+4.0fS*+0.0008P_o+0.012C+0.108T_a−0.0002P_in)	0.64
K_OS=1.6×10^−5+0.0001fO*−0.0001fS*+1.3×10^−5p_H−9.0×10^−6O_X−4.48×10^−7s_c+1.1×10^−7C	0.90

Note. Ta is mean annual surface air temperature in °C, Ox 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, Cx is the ratio of total soil C and oxalate metal oxide concentrations in g/mmol (Cx = 0.001 C/Ox), , are the fractions of solution P, labile P, sorbed P and occluded P in soil, respectively; and . Ps is the amount of sorbed P in soil (mg P/kg soil), Pin is the total amount of inorganic P excluding primary mineral P in mg P/kg soil, Po is the organic P in mg P/kg soil, p H is soil pH measured in water, ss, sc and si are sand, clay and silt percentages, respectively. All correlations are significant (p < 0.05).

The elemental effect of x i (x i = T a, or soil property) is calculated as where Δ 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.

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(∞)/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.

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 Equations 1 and 2). Values of the optimized parameters (b, k WL, k WS) 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.

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 (rW) and calculated values of r/R using the soil‐specific parameters (m, n, r(∞)/R and C∞) 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.

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 r2 is good, with r2 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 years−1, 0.02 years−1 and 0.005 years−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 modeled time course of r/R with that calculated using Equation 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 <4 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 modeled 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).

Figure 3

Comparison of the variation of isotopically labeled 33P radioactivity in soil solution (r/R) with time as calculated using Equation 11 with the soil‐specific parameters (m, n, r(∞)/R and C∞) (dark gray curve) with that predicted by the inorganic P model (light gray curve) for soil sample 1 (a) or 19 (b); and the modeled 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.

An indication of a poor model fit is a large difference between the modeled and measured fractions of different P pools at steady state. The differences in the modeled 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 modeled 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.

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 WL, k LW) and sorbed (k WS, k SW) through Equations 17 and 18. The optimized sorption coefficient being >100 mg P (kg soil)−1 day−1 (mg P/L)−b for labile P, and between 1 and 20 mg P (kg soil)−1 day−1 (mg P/L)−b for most soils (see Figure S2a in Supporting Information S1), the desorption rates of labile P (k WL) vary mainly between 0.01 and 10 days−1 among soil samples, whereas the desorption rates for sorbed P range between 0.001 and 1 day−1 (Figure S2b in Supporting Information S1). The sorption rates of occluded P vary across a larger range than sorption rates of other pools, that is, from 10−6 to 10−3 day−1 (Figure S2b in Supporting Information S1).

In theory, parameters m and n in Equation 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 WL <1 day) and weekly to yearly (1/k WS > 7 days) timescale for most soil samples (Figure S2 in Supporting Information S1), 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 Figures 4a and 4b).

Figure 4

Variations of the optimized values of parameters in the inorganic P model with IEK parameter n for 147 soils. kWL is the desorption rate constant in day‐1 for labile, and kSW is the sorption coefficient in mg P (kg soil)−1 day−1 (mg P/L)−b for sorbed P, respectively, and kOS is the sorption rate constant in day−1 for the occluded P, b is the exponent in the Freundlich equation for sorption, f S and f O are the fractions of sorbed and occluded P, respectively. The best fitted regression equations are: k LW = 12.24exp(−3.48n), r 2 = 0.17*; kWS/fS = 8.37exp(6.57n), r 2 = 0.40*; kOW/f O = 2.29 × 10−4−2.0 × 10−4n, r 2 = 0.1*; and b = 1.09–1.22n, r 2 = 0.14*.

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 OS/f 0) decreases with an increase in n n (r 2 = 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 in Supporting Information S1). Parameters m and n are negatively correlated (see Figure S4a in Supporting Information S1). 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 in Supporting Information S1).

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 toward estimating parameter values at regional or global scales in global land modeling.

One such parameter is the solution P concentration at steady state (C ∞). Figure 5a shows C ∞ increases exponentially with C/O x but reaches a saturation around C/Ox ∼ 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 in Supporting Information S1 and Table 2). All those four parameters depend on the fractions of inorganic P pools. In addition, k SW varies with T a, 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 in Supporting Information S1). The best fitted regression predicts the optimized value of k OS very well (see Figure S5 in Supporting Information S1).

Figure 5

(a) Variation of solution P concentration at steady state (C∞ in mg/L) with Cx, 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∞ by the best‐fitted regression and the observed values. The solid line represents 1:1 line.

Soil P Bioavailability

Different from previous studies, our definition of soil available inorganic P (a P) 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 p as defined in this study, we compared the estimated a p 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, Columbia. 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 Equation 25 and the optimized parameters for the four soils (see Table 3), we estimated the available soil P (A) as a p P T 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).

Table 3

Inorganic P Fractions, Total Inorganic P, Parameter b in Freundlich Equation, Desorption Rate Constant for Labile P (kWL), or Sorbed P (kWS) for the Top 0–20 cm Soils at the Four Sites in Carigmague, Columbia

	SAV	GL	CR	RGM
Fractions (labile/sorbed/occluded)	0.03/0.33/0.63	0.06/0.32/0.62	0.08/0.60/0.32	0.157/0.47/0.38
Total inorganic P (mg P/kg)	66	85	171	213
b	1.252	0.29	0.20	0.77
kWL (1/day)	0.0195	6.23	4.93	3.97
kWS (1/day)	0.0132	0.47	0.12	0.10
Plant available P (observed) (mg P/kg soil)	23	29	91	134
Our estimate (mg P/kg)	25	30	97	131

Note. 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). kWL and kWS were estimated from model calibration. The measured plant available was taken from Table 6 of Buhler et al., 2002) (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).

Furthermore, we explore why a p differs among the four different land uses at Carmagua, Columbia. As Equation 25 shows, a P depends on desorption rate, b and soil P fractions. The terms in Equation 25, (1‐exp(‐k WL t)) and (1‐exp(‐k WS t)) are very close to 1 at the four sites except (1‐exp(‐k WS 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 b) = 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.

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 (gray bar) to oxalate extractable metal oxide concentration (Ox), mean annual temperature (Ta), sand (ss), clay (sc) or silt (si) percentage, total soil C (C), total inorganic P (Pin), fractions of solution P (), labile P () or sorbed P (); (b) importance of different soil properties. Properties as listed in Table S1 that have little important (<10−3) or low sensitivity (<10−3) are not included here.

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 , and quite small (<0.1) for soil pH (not shown), T and soil texture variables. The absolute relative sensitivity increased for the fraction of solution P, T a 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.

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.

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 a (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 2 = 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(αT a)) with α 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 , and that the desorption rate for sorbed pool (k WS) that is proportional to k SW increases with fractions of sorbed P .

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 Achat, Pousse, et al., 2016; Borggaard et al., 2004; Hinsinger, 2001). 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 Regelink et al., 2015), which explains the negative sensitivity of sorption of labile P, k LW that is proportional to k WL by Equation 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 Equation 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 modeling P cycle in global land models.

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., Augusto et al., 2017; Ringeval 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, Pousse, et al. (2016), who identified the importance of soil C concentration and metal oxide concentration controlling inorganic P bioavailability. This is not unexpected, as two thirds of data used in this study were from Achat, Pousse, et al. (2016). But, by including the additional data (n = 45), we confirmed the findings of Achat, Pousse, et al. (2016).

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

Applications to Global Land Modeling

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 modeling. 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 data set and implementation of soil inorganic P model into global land models, a significant advance can be made in modeling of terrestrial biogeochemical cycles globally (see Wang & Goll, 2021).

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.

Supporting Information S1

Click here for additional data file.

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Figure S3

Click here for additional data file.

Figure S4

Click here for additional data file.

Figure S5

Click here for additional data file.