Improved representation of phosphorus exchange on soil mineral surfaces reduces estimates of phosphorus limitation in temperate forest ecosystems

Phosphorus (P) availability affects the response of terrestrial ecosystems to environmental and climate change (e.g., elevated CO₂), yet the magnitude of this effect remains uncertain. This uncertainty arises mainly from a lack of quantitative understanding of the soil biological and geochemical P cycling processes, particularly the P exchange with soil mineral surfaces, which is often described by a Langmuir sorption isotherm.

We first conducted a literature review on P sorption experiments and terrestrial biosphere models (TBMs) using a Langmuir isotherm. We then developed a new algorithm to describe the inorganic P exchange between soil solution and soil matrix based on the double-surface Langmuir isotherm and extracted empirical equations to calculate the sorption capacity and Langmuir coefficient. We finally tested the conventional and new models of P sorption at five beech forest sites in Germany along a soil P stock gradient using the QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system) TBM.

We found that the conventional (single-surface) Langmuir isotherm approach in most TBMs largely differed from P sorption experiments regarding the sorption capacities and Langmuir coefficients, and it simulated an overly low soil P-buffering capacity. Conversely, the double-surface Langmuir isotherm approach adequately reproduced the observed patterns of soil inorganic P pools. The better representation of inorganic P cycling using the double-surface Langmuir approach also improved simulated foliar N and P concentrations as well as the patterns of gross primary production and vegetation carbon across the soil P gradient. The novel model generally reduces the estimates of P limitation compared with the conventional model, particularly at the low-P site, as the model constraint of slow inorganic P exchange on plant productivity is reduced.

1 Introduction

Nutrient availability is one of the key factors affecting the productivity of terrestrial ecosystems and their future carbon (C) balance (Fernández-Martínez et al., 2014; Wieder et al., 2015). Although nitrogen (N) is the main constraint of plant biomass responses to an elevated CO₂ (eCO₂) concentration in many terrestrial ecosystems, phosphorus (P) availability likely constrains the biomass responses to eCO₂ in major global biomes (Du et al., 2020; Elser et al., 2007; Lebauer and Treseder, 2008; Terrer et al., 2019). For instance, the tropical forests and forests grown on old soils are known to be strongly limited by P availability, and they might not be able to sequester additional C in the future as the CO₂ concentration continues to increase (Hubau et al., 2020; Jiang et al., 2020). While temperate and boreal forests are generally considered to be N limited (Lebauer and Treseder, 2008), recent studies have shown that a decreasing P nutritional status is concomitant with increasing atmospheric CO₂ (Penuelas et al., 2020; Jonard et al., 2015). These findings highlight the importance of representing C–P interactions in terrestrial biosphere models (TBMs) and including P cycle processes in the future global C balance predictions.

There has been a continuous effort to include P cycling processes into TBMs in the past decade (Goll et al., 2012, 2017; Sun et al., 2021; Thum et al., 2019; Yang et al., 2014; Zhu et al., 2019; Wang et al., 2010; Yu et al., 2018). Many current TBMs employ the scheme (Fig. 1a) developed by Wang et al. (2007) to describe soil geochemical processes. This model considers the soil inorganic P (Pi) as soluble Pi (P_sol), labile Pi (P_lab), sorbed Pi (P_sorb), occluded Pi (P_ocl), and primary Pi (P_primary), with few exceptions where labile and sorbed Pi pool are grouped into one pool (Zhu et al., 2016, 2019). The exchange between soil solution and soil matrix is described with a Langmuir adsorption isotherm assuming that P_sol quickly exchanges with P_lab on the mineral surface. The P_lab pool also slowly exchanges with P_sorb at a linear rate regardless of P_sol concentrations, but this relationship is parameterized very differently across TBMs (Helfenstein et al., 2020).

Many modeling studies have emphasized the significance of biological P processes (Fleischer et al., 2019; Jonard et al., 2010; Yu et al., 2018; Wang et al., 2010) and underestimated the role of geochemical P processes. In these models, organic P recycling processes are the major fluxes, whereas the geochemical P fluxes are small (Sun et al., 2021). Particularly, the effect of (ad)sorption kinetics has seldom been discussed in previous modeling studies (Fleischer et al., 2019; Yang et al., 2014), although they are known to directly and strongly regulate P_sol and P_lab, thereby affecting P bioavailability (Frossard et al., 2000; Shen et al., 2011).

In this paper, we first conducted an extensive literature review of important soil characteristics affecting soil phosphorus sorption kinetics and then developed a model concept to reconcile measured P stocks with model simulations. In this new model, we applied the “two-surface” modification (Holford and Mattingly, 1976; McGechan and Lewis, 2002) to the conventional single-surface Langmuir isotherm to formulate a novel algorithm for Pi exchange with mineral surfaces, namely a double-surface Langmuir isotherm. We hypothesized that both P_lab and P_sorb exchange with P_sol in the new model (Fig. 1b), following evidence from isotopic studies (Helfenstein et al., 2020; Frossard et al., 2011). We then compared and evaluated the performance of the novel and conventional models with measured soil Pi pools as well as foliar N and P concentrations for a gradient of soil P stocks (164–904 g P m⁻²) and availability in a similar climate (mean annual temperature, MAT, 4.5–8 ^∘C) and under similar vegetation conditions (mature German beech forests, 120–140 years). Lastly, we tested the sensitivity of the alternative Pi exchange schemes on ecosystem P and C cycling to changes in P cycling parameters, and we also tested their responses to changes in P availability (P fertilization), C availability (CO₂ fertilization), and changes in both P and C availability.

2 Methods

2.1 Literature review on the Langmuir isotherm

The two parameters in the Langmuir isotherm, the maximum sorption capacity (S_max) and the Langmuir coefficient (K_m), are statistically fitted from measured adsorption curves of batch sorption experiments (Barrow, 1978). Following the scheme of Lloyd et al. (2001), most TBMs adapted the constant-equilibrium assumption to calculate a partition coefficient, k_p (Eq. 2), which is the ratio of change in the P_sol concentration to the change in available P (P_sol and P_lab). The partition coefficient affects the rate of P_sol released from the P_lab pool as P_sol is consumed and, thus, eventually affects the soil P supply to plants and microbes.

$$\begin{array}{}\text{(1)}& {\displaystyle }{P}_{\mathrm{lab}}=S_{\max }{\displaystyle \frac{P_{\mathrm{sol}}}{K_{\mathrm{m}}+P_{\mathrm{sol}}}}\text{(2)}& {\displaystyle }{k}_{\mathrm{p}}={\displaystyle \frac{\frac{\mathrm{d}{P}_{\mathrm{sol}}}{\mathrm{d}t}}{\frac{\mathrm{d}{P}_{\mathrm{sol}}}{\mathrm{d}t}+\frac{\mathrm{d}{P}_{\mathrm{lab}}}{\mathrm{d}t}}}={\displaystyle \frac{{\left(P_{\mathrm{sol}}+K_{\mathrm{m}}\right)}^{\mathrm{2}}}{{\left(P_{\mathrm{sol}}+K_{\mathrm{m}}\right)}^{\mathrm{2}}+S_{\max }{K}_{\mathrm{m}}}}\end{array}$$

To better understand the characteristics of phosphate exchange between solution and the soil mineral surface, we conducted a literature review of batch phosphate sorption experiments that describe the P sorption curves using a Langmuir isotherm, and we converted the lab-fitted Langmuir parameters (Q_max, mg P per kg soil, and K_L, L per mg P; see Sect. S2 in the Supplement) to the parameters used in TBMs, sorption capacity and the Langmuir coefficient (S_max and K_m, g P m⁻²; Eq. 1). From the batch phosphate sorption experiments found in the literature, we calculated the partition coefficient (k_p) following Eq. (2). Similarly, we reviewed the values of S_max and K_m in modeling studies and calculated k_p accordingly.

To showcase the difference in the Langmuir isotherm parameter values in TBMs and experiments as well as the differences between single- and double-surface Langmuir isotherms, we compared the exchangeable soil Pi curves of different TBMs' and batch experiments' data to the double-surface Langmuir isotherm. We also simulated desorption experiments with the single- and double-surface Langmuir isotherms to demonstrate their different responses to P removal.

2.2 The QUINCY model

QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system) is a terrestrial biosphere model of coupled C, N, and P cycles as well as energy and water processes (Thum et al., 2019). The model represents the growth of vegetation and the turnover of litter and soil organic matter at 30 min timescales, coupled with the calculation of the terrestrial energy and water budgets. Vegetation is represented by average individuals of plant functional types (here, a temperate broadleaved tree; see Sect. 2.3). Gross carbon uptake as well as leaf area and vegetation biomass and structure are directly influenced by nutrient availability through their effect on photosynthesis, tissue growth, allocation, and mortality. The model explicitly considers depth profiles (discretized into 15 layers with a default total soil column depth of 9.5 m) of soil temperature, moisture, and biogeochemical pools, representing litter and soil organic matter compartments, as well as inorganic forms of N and P. Vertical transport processes include diffusion and bioturbation. The phosphate exchange between soil solution and soil matrix is described using a conventional single-surface Langmuir isotherm, which is implemented the same way as in other TBMs. Different from other TBMs, QUINCY estimates S_max and K_m using soil texture and organic matter content (Thum et al., 2019).

To test our hypothesis that both P_lab and P_sorb exchange with P_sol, we modified the QUINCY Pi pool structure to describe the Pi exchange between solution and the soil matrix (Fig. 1b), using a double-surface Langmuir isotherm (Eq. 3) (Holford and Mattingly, 1976; McGechan and Lewis, 2002).

$$\begin{array}{}\text{(3a)}& {\displaystyle }{P}_{\mathrm{lab}}=S_{\max ,\mathrm{1}}{\displaystyle \frac{P_{\mathrm{sol}}}{K_{\mathrm{m},\mathrm{1}}+P_{\mathrm{sol}}}};\text{(3b)}& {\displaystyle }{P}_{\mathrm{sorb}}=S_{\max ,\mathrm{2}}{\displaystyle \frac{P_{\mathrm{sol}}}{K_{\mathrm{m},\mathrm{2}}+P_{\mathrm{sol}}}};\text{(3c)}& {\displaystyle }{S}_{\max }=S_{\max ,\mathrm{1}}+S_{\max ,\mathrm{2}};\text{(3d)}& {\displaystyle }{k}_{\mathrm{m}}={\displaystyle \frac{\frac{S_{\max }}{\mathrm{hlp}}-\mathrm{2}{P}_{\mathrm{sol}}\pm \sqrt{\frac{S_{\max }^{\mathrm{2}}}{{\mathrm{hlp}}^{\mathrm{2}}}-\mathrm{4}\frac{S_{\max }{P}_{\mathrm{sol}}}{\mathrm{hlp}}}}{\mathrm{2}}},\end{array}$$

where

$$\begin{array}{}\text{(3e)}& \mathrm{hlp}={\displaystyle \frac{S_{\max ,\mathrm{1}}{K}_{\mathrm{m},\mathrm{1}}}{{\left(K_{\mathrm{m},\mathrm{1}}+P_{\mathrm{lab}}\right)}^{\mathrm{2}}}}+{\displaystyle \frac{S_{\max ,\mathrm{2}}{K}_{\mathrm{m},\mathrm{2}}}{{\left(K_{\mathrm{m},\mathrm{2}}+P_{\mathrm{lab}}\right)}^{\mathrm{2}}}}.\end{array}$$

The first sorption sites are responsible for the fast exchange (Eq. 3a) between P_sol and P_lab and have much lower bonding strength than the second sorption sites that are responsible for the slower exchange (Eq. 3b) between P_sol and P_sorb. We estimated two soil sorption maxima (S_max,1 and S_max,2; Eq. 3c) and calculated two Langmuir coefficients (K_m,1 and K_m,2; Eqs. 3d and 3e) based on batch experiment data from the literature review, assuming both single- and double-surface isotherms can be fitted with batch experiment data. The detailed derivation is described in the Supplement (Sect. S3).

2.3 Sites and data

We performed analysis at five mature beech forest stands in Germany: Bad Brückenau (BBR), Mitterfels (MTF), Vessertal (VES), Conventwald (COM), and Lüß (LUE) (Table 1, based on Table 1 from Lang et al., 2017). Total soil P stocks (g P m⁻², up to 1 m depth) decrease strongly along the following gradient: BBR (904) > MTF (678) > VES (464) > COM (231) > LUE (164).

Figure 1 Model concept of (a) a conventional Langmuir model (siLang), where only the labile P pool (P_lab) is exchanging P with the soil solution (P_sol); (b) a double-surface Langmuir model (dbLang), where both P_lab and P_sorb are exchanging P with the soil solution but with different Langmuir parameters; and (c) a Langmuir model with only four inorganic P pools (4pool), where P_lab and P_sorb are combined to one pool (P_exchangeable) to which P_sol can get adsorbed via Langmuir sorption.

Soil was sampled up to 1 m depth at each site, with layer depths of 5–10 cm, for the measurement of total C, N, organic and inorganic P, and other physiochemical properties, such as soil texture, pH, and oxalate-extractable Fe and Al. Modified Hedley fractionations (Tiessen and Moir, 2008) were conducted on soils from all depths for the measurements of labile Pi (P extractable with anionic resin and Pi extractable with 0.5 M NaHCO₃), sorbed Pi (0.1 M NaOH extractable Pi), occluded Pi (P residual, acid digestion), and primary P (1 M HCl extractable P). Beech leaves were sampled in July and August from five trees for the measurement of leaf N and P concentrations.

Table 1 Site characteristics of the Bad Brückenau (BBR), Mitterfels (MTF), Vessertal (VES), Conventwald (COM), and Lüß (LUE) study sites, based on Table 1 from Lang et al. (2017). In the soil measurement section of the table, values before the slash refer to the topsoil, and values after the slash refer to the subsoil.

2.4 Model setup, experiments, and evaluation

2.4.1 Model setup

To compare the performance of different phosphate exchange schemes, we applied the conventional single-surface Langmuir isotherm here (siLang; Fig. 1a) as well as the novel double-surface Langmuir isotherm (dbLang; Fig. 1b) as model variants. To further test the causes of differences between these model variants, we introduced a Control simulation, which serves as a reference run without P limitation, and a four-pool simulation (4pool; Fig. 1c), which serves as a special case of siLang that does not include slow Pi exchange. All four models employed Eq. (3) to calculate the S_max and K_m in the Langmuir isotherms. In siLang and 4pool, total sorption capacity (S_max; Eq. 3c) only refers to P_lab (or P_exchangeable), whereas the two sorption maxima (S_max,1 and S_max,2) refer to P_lab and P_sorb, respectively, in dbLang. In the Control model, P_sol was kept at concentrations not limiting plant uptake or soil organic matter (SOM) decomposition. All models were initialized with a 15-layer soil column of 9.5 m which has a decreasing SOM content as soil depth increases (Thum et al., 2019). As for the inorganic P pools, the initial values of the top 1 m of soil were prescribed from the Hedley fractionation measurements and extrapolated to the deeper soil assuming an increased fraction of primary P and a decreasing fraction of exchangeable soil Pi pools (P_lab and P_sorb). The initialization of soil inorganic P pools deeper than 1 m is described in the Supplement (Sect. S4).

2.4.2 Model experiment protocol

The models were spun up for 500 years with meteorology and other atmospheric forcing (atmospheric CO₂ as well as N and P deposition), which were randomly drawn from the years from 1901 to 1930. During the model spin-up, the P cycle was simulated dynamically, but the more stable Pi pools, i.e., P_ocl, and P_primary, were kept constant to ensure that all of the models initialized at the same field P status. After spin-up, all models were run from 1901 to 2015 using the annual values for atmospheric CO₂, N (Lamarque et al., 2010, 2011) and P (Brahney et al., 2015; Chien et al., 2016) deposition, and the meteorology (Viovy, 2018) of the respective year. The C : N and C : P ratios of the slow SOM pool were calibrated per site per soil depth to match the measured soil C : N and C : P ratios so that the simulated SOM cycling adequately represented the measured site condition. All of the other parameter values were either taken from Thum et al. (2019) or modified based on the development in this study, except for the maximum biological N fixation rate, which is set to be 2.5 kg N ha⁻¹ yr⁻¹ for all of the study sites.

2.4.3 Model evaluation

Trend analyses were carried out with the Mann–Kendall test (M–K test) of the Kendall R package (McLeod, 2011). In the Mann–Kendall test, the tau (τ) value varies between −1 and 1, where −1 represents a decreasing trend and 1 represents an increasing trend. The modeled soil profile against the measured soil profile was evaluated with a normalized root-mean-square ratio term, K_nrmsr, which is modified to represent the average proportions between modeled and measured values (Yu et al., 2020b).

$$\begin{array}{}\text{(4)}& K_{\mathrm{nrmsr}}=\sqrt{{\displaystyle \frac{{\sum }_{\mathrm{1}}^n{K}_i^{\mathrm{2}}}{n}}}\text{where}{K}_i=min\left({\displaystyle \frac{{\mathrm{Mod}}_i}{{\mathrm{Meas}}_i}},{\displaystyle \frac{{\mathrm{Meas}}_i}{{\mathrm{Mod}}_i}}\right)\end{array}$$

K_i is the variable representing the ratio between simulated and measured values (in volumetric units) at the measured ith layer. A paired t test was conducted on the K_nrmsr of all study sites between siLang and dbLang to verify if one model was statistically better than the other.

2.5 Sensitivity analysis

To evaluate the response of alternative Langmuir isotherms to changes in P cycling processes, we further tested the sensitivity of the siLang and dbLang models to the P cycling parameterization at one low-P site, Conventwald (COM), and one high-P site, Mitterfels (MTF), using a standard Latin hypercube design (LHS; Saltelli et al., 2004). We selected 16 parameters that directly control flux rates of inorganic or organic P cycling processes and varied each parameter by ± 20 % of its default value (Table S2) using LHS sampling of a uniform distribution in order to form a set of 1000 LHS samples. The model outputs were evaluated in terms of gross primary productivity (GPP), foliar N and P content, vegetation C stock, soil organic carbon (SOC) content, total soil organic P (Po) and inorganic P, and the ratio between P_lab and exchangeable Pi. We measured parameter importance as the rank-transformed partial correlation coefficients (RPCCs) to account for potential nonlinearities in the association between model parameters and output (Zaehle et al., 2005; Saltelli et al., 2004).

2.6 Simulated CO₂ and P fertilization experiments

To test the effects of phosphate exchange schemes on environmental changes, we conducted a CO₂ fertilization model experiment, a P fertilization model experiment, and a CO₂–P (CP) fertilization model experiment using the siLang, dbLang, and 4pool models at each study site. In the CO₂ fertilization experiment, the atmospheric CO₂ concentration was increased by 200 ppm from 2006 to 2015. In the P fertilization experiment, 50 kg ha⁻¹ KH₂PO₄, i.e., 1139.7 mg P m⁻², was added once to the soil as soluble phosphate on 16 September 2006.

3 Results

3.1 Langmuir isotherm parameters in batch sorption experiments and TBMs

In batch sorption experiments, the maximum soil P sorption capacity S_max ranges from 187 to 829 g P m⁻² (25th quantile and 75th quantile values, median 390 g P m⁻²) and the Langmuir coefficient K_m ranges from 0.21 to 4.5 g P m⁻² (median 0.93 g P m⁻²); therefore, the calculated partition coefficient k_p varies between 0.005 and 0.022 (median 0.01; Table 2, Fig. S1). Regarding the parameterization of S_max and K_m in TBMs, a few studies have directly retrieved values from sorption experiments (Wang et al., 2007; Yang et al., 2014), although most modeling studies estimate S_max and K_m based on soil types (Wang et al., 2010; Goll et al., 2012; Zhu et al., 2019) or soil texture (Thum et al., 2019).

The two studies using lab-derived Langmuir parameters fitted well in the range of experimental values; however, for those modeling studies that estimate Langmuir parameters, only the QUINCY model (Thum et al., 2019) produced reasonable Langmuir parameters and k_p values, while most other TBMs greatly overestimated K_m, thereby generating much higher k_p values than batch experiments (Table 2).

Table 2 The maximum sorption capacity (S_max), the Langmuir coefficient (K_m), and the calculated partition coefficient (k_p) of the phosphorus Langmuir (ad)sorption isotherm in TBMs and batch sorption experiment data (BED).

n/a: not applicable.

In Fig. 2a, it is also shown that most TBMs are at or below the lower range of curves from batch sorption experiment data (BED) except the two isotherms from QUINCY that are close to the BED-median curve. All of the other TBMs either have very low S_max or overly high K_m, which can both lead to flattened soil Pi curves, indicating that very low amounts of P are stored as exchangeable Pi in the soil, even at high P_sol concentrations.

However, the responses to P removal are noticeably different between the single- and double-surface Langmuir models (siLang and dbLang) given the similar soil Pi curves (Fig. 2b). The P_sol in siLang almost reached zero after a 30 d desorption experiment, whereas dbLang still maintained a P_sol level of 0.2 g m⁻², suggesting a much stronger buffering capacity in dbLang compared with siLang. This is mainly because the replenishment of P_sol in siLang is strongly limited by the rate of desorption from P_sorb to P_lab, but P_sol directly exchanges with P_sorb when the P is removed in dbLang.

Figure 2 (a) Exchangeable soil inorganic P (P_lab and P_sorb) curves based on different terrestrial biosphere model (TBM) parameters and batch experiment data (BED). In TBMs, exchangeable soil Pi is defined as the sum of P_lab and P_sorb. In BED, soil exchangeable Pi only refers to P_lab in Eq. (1) due to a destructive experimental procedure that mobilizes the soil Pi well (Barrow and Shaw, 1979). For the conventional Langmuir model (siLang), P_sorb is calculated as 9 : 8 (global average ratio between P_sorb and P_lab, reported in Yang et al., 2013) of P_lab at all soluble P concentrations; for the double-surface Langmuir model (dbLang), P_sorb is calculated following Eq. (3b). (b) Simulated desorption curves for single- and double-surface Langmuir isotherms. The two desorption curves start with the same exchangeable soil Pi (152.2 g P m⁻², P_lab+ P_sorb) and soluble Pi (0.769 g P m⁻², P_sol) content, and the same amount of P (2.5 g P m⁻²) is removed from soil solution every day. Both Langmuir isotherms are assumed to be in equilibrium at a daily time step.

3.2 Simulated and measured ecosystem properties

The siLang, dbLang, and 4pool models can adequately reproduce the measured SOC content as well as the SOM C : N and C : P ratios along soil profiles at all of the study sites after site- and depth-specific calibration (Fig. 3a, b, c; Table S1). The performance of the three models did not differ much regarding the simulation of SOM profiles due to our calibration strategy (Tables 2, S1). However, the novel dbLang model better reproduced the ratio between P_lab and exchangeable Pi (Lab-to-Exchangeable P) than siLang (K_nrmsr improvement of 0.130 ± 0.077, p<0.05) (Table 3). The improvements of dbLang in modeling the labile (K_nrmsr improvement of 0.030 ± 0.147, p>0.05) and sorbed Pi (K_nrmsr improvement of 0.096 ± 0.137, p>0.05) pools were not as significant as that in modeling the Lab-to-Exchangeable P ratio (Fig. 3d, e, f; Table 3).

Figure 3 Simulated and measured (a) soil organic carbon (SOC), soil organic matter (SOM) (b) C : N and (c) C : P ratios, (d) ratio of labile to exchangeable Pi, (e) labile Pi, and (f) sorbed Pi of the study sites along the following soil P availability gradient: BBR > MTF > VES > COM > LUE. The black dotted line is the field measurements, the lime line represents the Langmuir model with only four inorganic P pools (4pool), the blue line represents the single-surface Langmuir model (siLang), the red line represents the double-surface Langmuir model (dbLang), and the light sky-blue line represents the Control model. The 4pool and Control simulations are not applicable in panels (d), (e), and (f).

Table 3 Results for the paired t test between the normalized root-mean-square ratios (K_nrmsr; Table S1) of the siLang and dbLang models of measured soil profile properties at all study sites. The null hypothesis is that dbLang performs no differently than siLang (i.e., dbLang has the same K_nrmsr values as siLang), and the hypothesis is rejected when the p value is smaller than 0.05 (values in bold with an asterisk).

The simulated average foliar N content values of the four models (23.35 ± 1.59, 23.05 ± 1.84, 22.45 ± 2.35, 22.74 ± 2.27 mg N g⁻¹ DW (where DW denotes dry weight) for Control, dbLang, siLang, and 4pool, respectively) were within the average range of measured values (24.32± 1.43 mg N g⁻¹ DW). However, the simulated decreasing trend in the foliar N content (M–K test, $\mathit{\tau }=-\mathrm{0.6}$, −0.8, −1, −0.8 for Control, dbLang, siLang, and 4pool, respectively) along the soil P gradient was not found in the measured data (M–K test, τ = −0.2) (Fig. 4a). Decreasing trends in simulated foliar P content in all models (M–K test, $\mathit{\tau }=-\mathrm{0.8}$, −0.95, −1, −1 for Control, dbLang, siLang, and 4pool, respectively) were instead also seen in measurements, but the decreasing trend was much weaker (M–K test, $\mathit{\tau }=-\mathrm{0.6}$) (Fig. 4b). The simulated foliar P content was highest in the Control model (1.40 ± 0.09 mg P g⁻¹ DW) and lowest in siLang (0.95 ± 0.21 mg P g⁻¹ DW) at each study site, as Control is not P limited and the strongest P limitation occurs in siLang. The simulated difference in foliar P content across models reflects both plant P uptake and productivity (Figs. 4, S2–5). For example, at the P-rich BBR site, although the simulated GPP, leaf area index (LAI), aboveground C, and fine-root C for the four models were almost identical (Fig. S2), the foliar P content values of four models were different due to differing plant P uptake (0.99, 0.95, 0.89, 0.95 g P m⁻² yr⁻¹ for Control, dbLang, siLang, and 4pool, respectively; Fig. S2). In contrast, at the P-poor LUE site, the differences among siLang, dbLang and 4pool for GPP, LAI, and plant C were more pronounced than that for foliar P content, as the limited P uptake (0.67, 0.53, 0.55 g P m⁻² yr⁻¹ for Control, dbLang, and siLang, respectively; Fig. S2) drastically changed the plant C allocation (e.g., leaf-to-root ratio) and led to subsequent changes in plant properties (Figs. 4, S2).

Figure 4 Simulated and measured foliar N (a) and P (b) content along the soil P availability gradient. The simulated values are the yearly average of the 2006–2015 period, and the measured values were sampled in 2014. The black dots are the field measurements, the lime line represents the Langmuir model with only four inorganic P pools (4pool), the blue line represents the single-surface Langmuir model (siLang), the red line represents the double-surface Langmuir model (dbLang), and the light sky-blue line represents the Control model.

3.3 Sensitivity analysis

Our sensitivity tests at the low-P COM site (231 g P m⁻²) and the high-P MTF site (678 g P m⁻²) showed a diverging effect of model choice and parameterization on selected ecosystem properties, such as GPP, plant C, and foliar P content (Figs. 5, S6). Both models were very consistent in maintaining the vertical pattern of the simulated ratio between labile and exchangeable Pi (Lab-to-Exchangeable P) across different parameterizations (Fig. S7). In other words, the better fit of dbLang in capturing the measured decreasing trend in Lab-to-Exchangeable P was robust against model parameters choice (Figs. S6, S7).

The P cycling parameters with the strongest effect on the selected ecosystem properties were very different between the two models and moderately different between the two sites (Fig. 5, Table S3). In siLang, the parameters with the highest impact for most selected outputs (e.g., foliar P, GPP, and plant C, as in Fig. 5) were the absorption and desorption rates between P_lab and P_sorb (k_abs and k_des), i.e., the slow exchange process, whereas they were the turnover rates and N : P ratios of slow and fast SOM pools (τ_slow, τ_fast, SOM_np, and microbial_np) in dbLang. One common feature of both models is that the low-P site (COM) was more affected by SOM_np, the maximum weathering rate (k_weath_mineral) and the maximum plant P uptake rate (vmax_uptake_p_p4) than the high-P site (MTF), inferring a tighter P cycle in a low-P ecosystem than in a high-P ecosystem.

Figure 5 Partial correlation coefficient (pcc) values of eight main parameters of foliar P content, gross primary production (GPP), and vegetation C stock at the COM and MTF sites of the single-surface Langmuir model (siLang) and the double-surface Langmuir model (dbLang) in the Latin hypercube design (LHS) sensitivity runs. The eight parameters are as follows: the maximum plant P uptake rate (vmax_uptake_p_p4), the turnover rates of slow and fast soil organic matter (SOM) pools (tau_slow and tau_fast, respectively), the N : P ratio of the slow SOM pool (SOM_NP), the phosphate sorption capacity of fine soil (qmax_po4_silt), the correction coefficient of pH on Langmuir K_m (km_po4_ph), and the absorption and desorption rate between P_lab and P_sorb (k_ adsorpt_po4 and k_desorpt_po4, respectively). The significance levels shown in figure (using asterisks) are as follows: p<0.001 (${}^{***}$), p<0.01 (${}^{**}$), and p<0.05 (^*).

3.4 Modeled ecosystem responses to C, P, and CP fertilizations

The simulated ecosystem responses to CO₂ fertilization, P fertilization, and CP fertilization differed greatly among study sites, models, and fertilization types due to the differences in soil P availability and model schemes of Pi exchange (Figs. 6, 7, S3–S5). At the P-rich BBR site, the simulated GPP, LAI, aboveground and fine-root C, and plant uptake of N and phosphate increased after CO₂ and CP fertilizations but did not change after P fertilization in all three models (siLang, dbLang, and 4pool; Fig. 6a). This indicates that the simulations at BBR were not limited by P, which is also supported by the small difference in the responses between the CO₂ and CP fertilization experiments. In contrast, at the P-poor LUE site, the simulated responses to P and CP fertilization were much stronger than those after CO₂ fertilization in all three models (Fig. 6c), indicating a strong P limitation at LUE.

This simulated P stress among models can be quantified when compared to the Control model (Fig. 6d). The siLang model simulated the lowest plant P uptake at all three sites with high, moderate and low P concentrations in soil. This difference in plant P uptake was not only reflected in the foliar P concentration at the high-P site but also in other vegetation properties at moderate- to low-P sites. The 4pool model, which is a special case of the single-surface Langmuir isotherm, simulated similar P stress to dbLang at high- and moderate-P sites, whereas it behaved similarly to siLang at the low-P site and simulated a high P stress (Fig. 6d).

Figure 6 Simulated changes in gross primary production (GPP), aboveground C, fine-root C, leaf area index (LAI), uptake of N and P, and foliar N and P content due to CO₂, P, and CP fertilizations at the (a) BBR, (b) VES, and (c) LUE sites, and (d) the ratio of unfertilized scenarios of each model to the Control model. Bars are calculated as the percentage change in subplots (a) to (c), and the bars in subplot (d) are fractions representing how much the models deviate from a non-P-limited condition.

Figure 7 Simulated temporal responses of model P pools (labile Pi; sorbed Pi; litter Po; Po in fast and slow soil organic matter (SOM) pools – fast-SOM_Po and slow-SOM_Po, respectively; and P in the plant – Veg_P) to CO₂, P, and CP fertilizations at the BBR, VES, and LUE sites. In the C fertilization experiment, the atmospheric CO₂ concentration was increased by 200 ppm from 2006 to 2015. In the P fertilization experiment, 1139.7 mg P m⁻² was added once to the soil as soluble phosphate on 16 September 2006. The CP fertilization is the combination of C and P fertilizations. Each point in the x axis represents a P budget change at the specific time, with the size change in all P pools equal to zero (C fertilization) or the added P (1139.7 mg P m⁻², P and CP fertilization). The change in a specific P pool at any given time is equal to the total change in related P fluxes since the fertilization (e.g., the plant P pool change relates to plant P uptake and P litterfall, the labile Pi pool change relates to P adsorption (and P absorption) (Fig. 1), the Po in the fast SOM pool change relates to litter and fast SOM decomposition, and so on).

The changes in the model P pools after CO₂, P, or CP fertilizations further illustrated the effects of Pi exchange schemes on ecosystem responses under varying soil P availability (Fig. 7). In the CO₂ fertilization experiment, the chronic increase of CO₂ led to increases in plant biomass P at all study sites in all three models (Fig. 7, top row). The concurrent increases in plant P pools were compensated for by decreases in labile and sorbed Pi as well as Po in the microbial pool (fast SOM). At the high-P BBR site, the increases in plant P (631.4 ± 1.4 mg P m⁻²), litter P (107.8 ± 0.8 mg P m⁻²), and slow Po in SOM (62.8 ± 0.3 mg P m⁻²) were similar in all models, and the mobilization from soil exchangeable Pi (−417.4 ± 16.1 mg P m⁻²) and fast Po SOM (−354.2 ± 2.5 mg P m⁻²) contributed similarly to P sources. Conversely, at the P-poor LUE site, dbLang simulated a significantly higher plant P increase (188.3 mg P m⁻²) than siLang (69.5 mg P m⁻²) and 4pool (18.5 mg P m⁻²), due to its stronger capacity to mineralize P from the microbial P pool. dbLang can maintain a much higher SOM pool in topsoil (Fig. 3a), as P_sol is better buffered under low soil Pi (Fig. 2b) compared with siLang; thus, more P can be mineralized under eCO₂.

In the P fertilization experiments, the fate of added P largely differed between siLang and dbLang, as the added P was preferably and quickly transferred to P_sorb in dbLang compared with siLang (Fig. 7, bottom row).

In the CP fertilization experiments (Fig. 7, middle row), plant P did not gain more P at high- and moderate-P sites (BBR and VES) compared to the C fertilization, indicating both sites might not be P limited under eCO₂. Surprisingly, microbes (fast SOM) also did not benefit from P addition on top of eCO₂ at BBR and VES, as most added P was transferred to the soil Pi pools. However, at the low-P LUE site (Fig. 7, last three columns), the combination of C and P fertilization produced higher plant P increases (928.5, 1202.3, and 971.1 mg P m⁻² for 4pool, dbLang and siLang, respectively) than P fertilization alone in all three models, inferring a very strong P limitation on C sequestration at the low-P LUE site.

4 Discussion

4.1 Model representation of soil Pi cycling and Langmuir sorption

The majority of terrestrial biosphere models (TBMs) describe the soil inorganic P (Pi) cycling processes with a similar structure of pools and fluxes (Fig. 1a or c), and they mostly describe the key soil Pi exchange process, i.e., the (ad)sorption and desorption between soil solution and matrix, using the Langmuir isotherm. In this study, we have shown that the values of Langmuir parameters in most TBMs, namely maximum P sorption capacity (S_max; Eq. 1) and the half-saturation Langmuir coefficient (K_m; Eq. 1), largely differ from those reported in the P sorption batch experiments (Table 2). This disagreement between TBMs and experiments leads to much lower exchangeable soil Pi (P_lab and P_sorb) content in TBMs than experiments under the same soluble Pi (P_sol) conditions (Fig. 2a). This is in line with the much higher partition coefficient (k_p; Eq. 2) values in most TBMs compared with experiments (Table 2). Additionally, the reported S_max values have a comparable size to the total soil Pi in experiments (Fig. S1), indicating that most of the soil-associated Pi could be exchanged with the soil solution at a certain point. This implies that the Langmuir isotherm can describe not only the fast exchange between P_sol and P_lab but also the slower exchange with P_sorb. This hypothesis is supported by the classic sorption model of Barrow (Barrow, 1978, 1983) and the isotopic exchange kinetics conceptual model (Fardeau, 1995; Frossard and Sinaj, 1997; Morel et al., 2000) that Pi ions located in the soil matrix are distributed along a continuum of solubility, with some being in rapid equilibrium with Pi ions in soil solution and some being in slow equilibrium.

Following the concept of the soil P solubility continuum, we have developed the double-surface Langmuir isotherm (dbLang; Fig. 1b, Eq. 3) in the QUINCY TBM and parameterized it with batch sorption experiment data. The double-surface Langmuir isotherm simulated a higher P-buffering capacity than the single-surface Langmuir isotherm (siLang) when P is removed from the soil (Fig. 2b), which is in agreement with the experimental results of Roberts and Johnston (2015). By comparing dbLang with the conventional siLang model and a special four-Pi-pool structure of siLang (4pool model, used in Zhu et al., 2016, 2019), we further investigated the effects of different sorption–desorption schemes on soil P under the QUINCY TBM framework (Fig. 3). The model performance in simulating the relationship between soil labile and sorbed Pi pools (Lab-to-Exchangeable P) has been improved by dbLang at our study sites of varying soil P stocks, soil texture, and stoichiometry (Lang et al., 2017; Yu et al., 2020b). Such an improvement of dbLang with respect to reproducing the Lab-to-Exchangeable P ratio much better captures the vertical patterns of soil inorganic pools (Figs. 3d, e, f), which is also supported by statistical analysis that the measured soil P profiles were slightly better reproduced (Tables 3, S1). This was not simply caused by improvements in modeling the individual labile or sorbed Pi pool but rather by a better representation of the Pi exchange among P_sol, P_lab, and P_sorb, as such an improvement is independent of site conditions and model parameterization (Figs. 3, S7). Additionally, the siLang and dbLang models in this study shared the same set of QUINCY parameters; therefore, they had the same model complexity but differed in model performance. However, the new scheme did require more input than other TBMs (listed in Table 2), as the new parameterizations were dependent on soil texture, pH, and extractable Al or Fe, whereas most other TBMs use prescribed values.

In this study, we have included the Hedley soil P data as part of the model validation, which, to our knowledge, has seldom been done before, although the Hedley soil P data (Yang and Post, 2011; Yang et al., 2013) have been widely used as initialization data for TBMs. The depth-explicit and quantitative evaluation of model Pi pools with field measurements provides us a unique tool to diagnose different sorption–desorption schemes, and it warrants future applications in model–data comparison or model intercomparison studies.

4.2 Model performance and confidence

We have shown that dbLang performed better than siLang (and 4pool) in reproducing the measured soil Pi pools and foliar P content (Figs. 3, 4; Table 3). Several observations suggest that the novel dbLang model produced more realistic GPP and vegetation C patterns, although there is no direct evidence. First of all, Lang et al. (2017) found that there are no trends in the measured tree heights nor stand biomass along our soil P gradient (Table 1). As all sites are unmanaged beech stands with a similar age structure, we can assume that GPP and vegetation C should also be similar between the sites in the absence of nutrient limitation. This was better reproduced using dbLang than siLang, as less P limitation is simulated using dbLang. Secondly, studies have shown that biological processes dominate the C–P interactions in P-poor beech forests (Bünemann et al., 2016; Pistocchi et al., 2018); this was supported by the sensitivity results of only dbLang (Fig. 5, Table S3), whereas slow sorption/desorption process dominated GPP and vegetation C regardless of P availability in siLang. Most importantly, siLang overestimated P_lab and P_sorb in topsoil but still greatly underestimated the plant productivity and biomass at low-P sites (COM and LUE). Collectively, these findings suggest that a double-surface Langmuir model of Pi sorption–desorption better described forest C and nutrient dynamics at these sites than the conventional single-surface Langmuir model.

Yang et al. (2014) used the siLang scheme to study P cycling in P-limited tropical forests. They showed that a doubled rate of P transfer from P_sorb to P_lab leads to similar increases in plant productivity and biomass as a doubled P mineralization rate. This is in line with our sensitivity analysis of siLang, showing that slow sorption/desorption of P has similar impacts to SOM turnover on GPP and plant biomass (Fig. 5, Table S3). However, this seems unlikely in tropical forest sites depleted of soil Pi, and we argue that a dbLang model could also better capture such processes in tropical forest ecosystems.

From the perspective of Pi cycling, the parameterization of linear P exchange between P_lab and P_sorb in siLang lacked clear foundations and experimental evidence; hence, it largely differed among TBMs (Helfenstein et al., 2020). Similarly, we have shown that the exchange between P_sol and P_lab in most TBMs also largely varies and deviates from the reported ranges of experiments (Table 2). Most importantly, the assumption that P_sol only exchanges with P_lab has been heavily questioned by experimentalists (Helfenstein et al., 2020; Morel et al., 2000), and it also differs from earlier models (Barrow, 1983; Devau et al., 2009; van der Zee and Gjaltema, 1992). The double-surface Langmuir isotherm much better represents the conceptual sorption–desorption scheme that is supported by experimental data (Frossard et al., 2000; Frossard and Sinaj, 1997; Fardeau, 1995), and it also implicitly considers some slower Pi exchange processes that are currently ignored in TBMs, such as chelation and dissolution/precipitation.

4.3 Ecosystem responses to perturbations

The different Pi sorption–desorption schemes have not only led to different simulated Pi pool sizes but have also caused varying ecosystem responses to elevated CO₂ and P addition (Figs. 6, 7). The main difference between siLang and dbLang was the speed and extent of exchange between the soluble/labile and sorbed Pi pools. The added P was quickly transferred to the P_lab in both models, but the speed of Pi transfer to P_sorb was much faster in dbLang. This was because Pi directly transferred from P_sol to P_sorb in dbLang, whereas the transfer to P_sorb only occurred from P_lab and at a much slower rate in siLang. As a consequence, P_sorb in dbLang stored more P after P addition, and it also released more P under higher CO₂ concentrations and elevated plant P demand. The fast transfer of P from P_sol to P_lab in both models is supported by evidence from ³³P tracer studies showing that the radioactivity of P_sol and P_lab converged quickly (<3 months) after ³³P tracer addition (Pistocchi et al., 2018). However, the speed and extent of transfer to P_sorb cannot be easily confirmed with isotope signals because both ³³P and ³²P have short half-lives, but it is estimated to vary from days to weeks (Buehler et al., 2002; Helfenstein et al., 2021). Nevertheless, the mathematical and conceptual descriptions of isotopic exchange kinetics tend to support the dbLang model in that the P transfer to P_sorb took place at a timescale of months rather than years in acidic soils (Bünemann et al., 2016; Frossard and Sinaj 1997; Helfenstein et al., 2020). Hou et al. (2019), using a data assimilation approach, also concluded that P transfer to P_sorb might happen at a much faster rate than the conventional TBM parameter values suggest.

Our advance in capturing soil Pi exchange has altered the responses of QUINCY to changing P availability, particularly at low-P sites. The main difference in the soil Pi responses to P additions between siLang (also 4pool) and dbLang is that the novel dbLang model can store or release more P from the sorbed Pi pool at a faster rate than siLang. Particularly, at sites with high P stress (e.g., LUE), a further increase in P stress (e.g., CO₂ fertilization) led to an increased P mineralization in dbLang rather than increased P desorption in siLang (Fig. 7), as dbLang can establish a much larger (implicit) microbial pool (fast SOM) in the topsoil to mineralize P (Fig. 3a). It was confirmed in the EucFACE (Eucalyptus Free Air CO₂ Enrichment) experiment that both N and P mineralization increased by as much as 200 % under eCO₂ (Hasegawa et al., 2016). Additionally, evidence from isotopic studies has shown that biological processes, rather than geochemical processes, are dominant in both topsoil and subsoil of the P-poor LUE site (Pistocchi et al., 2018; Bünemann et al., 2016).

4.4 Limitations and future directions

The novel dbLang model much better represented soil Pi exchange and, thus, improved the C–P interactions in QUINCY, but there was a caveat when the P stress was high. At the P-poor site (LUE), the GPP and vegetation C in dbLang were almost twice as high as those in siLang (Fig. S2a, c). However, such a huge release of P stress by dbLang was still not enough for LUE to reach similar plant productivity and biomass values to BBR or VES, as indicated by field evidence (Lang et al., 2017). In ecosystems with low soil P or high P stress, CO₂ fertilization enhances the organic P cycling processes, such as root exudation, phosphatase production, and microbial P mining (Ellsworth et al., 2017; Jiang et al., 2020; Lang et al., 2017; Pistocchi et al., 2020). Such responses in SOM cycling have not yet been implemented in QUINCY, leading to much lower GPP and plant biomass values at the low-P site than at other sites. These mechanisms have also not been described in other TBMs, leading to poor model performance in low-P ecosystems, e.g., Amazon forests (Fleischer et al., 2019) and Eucalyptus forests (Medlyn et al., 2016). Recent developments in soil models/modules have endeavored to improve the model description of SOM cycling by including parameters such as organomineral association process, explicit microbial dynamics, enzyme dynamics or allocation, and mycorrhizal association (Huang et al., 2021; Yu et al., 2020a; Ahrens et al., 2020; Wutzler et al., 2017; Tang and Riley, 2014; Sulman et al., 2014; Wang et al., 2013). We believe that coupling the novel dbLang scheme with a more mechanistic representation of SOM cycling will open the door to model the different forest ecosystem P cycling strategies (Lang et al., 2017) as well as to simulate the responses of low-P ecosystems to elevated CO₂ (Jiang et al., 2020).

5 Conclusions

In this study, we first reviewed the implementation of the soil Pi exchange process in terrestrial biosphere models and compared the model implementations of P sorption with batch sorption experiment data. We found that the parameterizations used by most TBMs strongly underestimate the soil's P sorption capacity and overestimate the half-saturation Langmuir coefficient compared with the experimental data. In the QUINCY model, such a formulation leads to a much lower soil P binding capacity in TBMs than in reality, causing the unrealistically large constraint of a slow Pi exchange process on GPP and plant biomass, regardless of the actual soil P availability. We presented a novel model scheme, based on a double-surface Langmuir isotherm, to describe the soil Pi exchange in better accordance with data from sorption experiments. This model parameterization better simulated the measured soil Pi pools and the GPP and plant C patterns of our study sites. It also better reproduced the topsoil SOM pools at the low-P site, thereby better simulating the responses of P pools/fluxes to elevated CO₂ than conventional TBMs. Thus, the novel double-surface Langmuir approach can serve as a better modeling tool to understand the ecosystem response to global change.