Traditional static soil tests (STPs) for phosphorus (P) often fail to predict agronomic outcomes because they do not account for the kinetic nature of P supply to plant roots. This study investigated whether P desorption kinetic parameters could serve as more effective predictors. Using soils from the long-term Swiss agricultural experiment (STYCS), a sequential extraction method was refined and modeled with a non-linear approach to derive the maximum equilibrium P intensity (\(P_{desorb}\)) and a rate constant (\(k\)). The predictive power of these kinetic parameters was compared against standard Swiss STPs (\(P_{CO_2}\) and \(P_{AAE10}\)) using linear mixed-effects models. For predicting site-normalized yield, STPs were superior, while for national-normalized yield and P-uptake, both methods performed poorly due to dominant pedoclimatic factors. However, the kinetic model demonstrated exceptional success in predicting the long-term P-Balance, with \(P_{desorb}\) alone explaining 57% of the variance, whereas STP models showed no predictive power. We conclude that while traditional STPs remain adequate for short-term within-field fertility, kinetic parameters serve as a vastly superior functional proxy for assessing the long-term P status and sustainability of agricultural soils.
Phosphorus (P) is an essential macronutrient for all life, yet its management in agricultural systems presents a profound dilemma. In the soil solution, P exists almost exclusively as highly reactive orthophosphate anions (\(HPO_4^{2-}\) or \(H_2PO_4^{-}\)). These dissolved species are rapidly immobilized—either adsorbed onto the surfaces of clays and metal oxides or precipitated with cations like calcium, iron, and aluminum to form minerals of low solubility (Brady & Weil, 2016; Sposito, 2008). Consequently, the total soil P stock is vast, but the infinitesimally small fraction dissolved in the soil solution—the only form directly available to plant roots—is highly transient.
To conceptualize plant P availability, soil science relies on three fundamental thermodynamic and mass-balance parameters (Frossard et al., 2000): 1. P-Intensity (\(I\)): The thermodynamic energy state, or chemical activity, of P in the soil solution at any given moment. 2. P-Quantity or Capacity (\(Q\)): The physical mass of solid-phase P that is readily desorbable and can replenish the soil solution. 3. Phosphorus Buffering Capacity (PBC): The resistance of the soil to changes in P-Intensity, mathematically represented as the derivative of Quantity with respect to Intensity (\(dQ / dI\)).
While these concepts are theoretically elegant, measuring them in natural soils is notoriously difficult. In practice, standard agronomic soil tests do not measure absolute thermodynamic states; rather, they rely on operational definitions. A soil test value is strictly defined by the specific chemical extractant used, the extraction time, and the soil-to-solution ratio.
Different standard methods target different physicochemical pools of P by employing distinct extraction mechanisms: - Olsen P (\(NaHCO_3\), pH 8.5): Measures desorbable P. The high concentration of bicarbonate ions outcompetes orthophosphate for adsorption sites on soil surfaces (Olsen1954?). - P-AAE10 (Ammonium Acetate EDTA, pH 4.65): A strong extractant that attacks both sorbed and precipitated P forms. EDTA chelates background cations, actively dissolving apatites, while the acetate buffer facilitates desorption (Forschungsanstalt für Agrarökologie und Landbau (FAL), 1996). - Water Extractions (\(P_{H_2O}\) and \(P_{CO_2}\)): Weak extractions using distilled or \(CO_2\)-saturated water exert minimal chemical force, attempting to extract readily soluble forms without fundamentally altering the native soil matrix.
The reliance on operational definitions introduces a critical analytical flaw when interpreting soil tests, fundamentally tied to the principles of physical chemistry and plant physiology. Plant roots acquire P primarily via high-affinity transporter proteins (e.g., the PHT1 family). The binding kinetics of these transporters are governed by the true thermodynamic driving force of the substrate: its chemical activity (\(a\)), not its total stoichiometric concentration (\(c\)).
In the native soil solution, phosphate anions are electrostatically shielded by a cloud of background ions (\(Ca^{2+}, Mg^{2+}, K^+, NO_3^-, SO_4^{2-}\)). This native ionic strength (\(I\)) physically reduces the effective activity of the phosphate. The magnitude of this shielding is defined by the activity coefficient (\(\gamma\)), where true chemical activity is \(a = c \cdot \gamma\).
For strong extractants like Olsen or P-AAE10, the extractant solution has an overwhelmingly high ionic strength (e.g., \(I \approx 0.5 \text{ M}\) for Olsen). This artificial salt load obliterates the native ionic strength of the soil. Consequently, while these methods successfully quantify the physical mass of phosphorus (\(Q\)), they destroy the natural chemical environment, providing no information regarding the true thermodynamic Intensity (\(I\)) perceived by the root.
Conversely, while weak water extractions successfully preserve the native soil matrix, reading their variable activity against a standard \(I \approx 0\) calibration curve renders the values mathematically incomparable across different pedoclimatic sites. For example, in heavily fertilized Swiss agricultural soils (such as the long-term STYCS 167% P-treatments), the dissolution of accumulated amendments can elevate pore water ionic strength to \(0.02 \text{ M}\). At a neutral pH, a large portion of P exists as the divalent \(HPO_4^{2-}\) ion. Applying the Davies equation to this ionic strength yields an activity coefficient (\(\gamma\)) of approximately \(0.58\). Therefore, over 40% of the chemically extracted mass is thermodynamically “masked” from the root. To resolve this, raw stoichiometric concentrations derived from water extractions must be mathematically corrected via the Davies equation to yield true thermodynamic activity.
Once thermodynamic Intensity is correctly calculated, it can be paired with a standard mass extraction (\(Q\), such as \(P_{AAE10}\)) to construct a purely thermodynamic \(Q(I)\) framework. By parameterizing this relationship using a Freundlich isotherm (\(Q = K_f \cdot I^{1/n}\)), the true Phosphorus Buffering Capacity (PBC) can be mathematically deduced as the continuous first derivative (\(\frac{dQ}{dI}\)).
However, even with a correct thermodynamic understanding of \(Q\) and \(I\), plant uptake remains fundamentally limited by time. Standard conceptual models assume instantaneous solid-solution equilibrium. In reality, desorption is a rate-limiting kinetic bottleneck (Flossmann & Richter, 1982; Nye & Tinker, 2000). When a root depletes the rhizosphere, the soil cannot always resupply phosphate fast enough.
To address this, we define a novel dynamic predictor: the Initial Activity Flux (\(J_0\)). This flux shifts the focus from static mass to supply rate by multiplying the first-order desorption rate constant (\(k\)) by the maximum thermodynamic activity (\(a_{max}\)): \[J_0 = k \cdot a_{max}\]
?@tbl-terminology maps this transition from classic, operationally defined static tests to the proposed thermodynamic and kinetic parameters.
| Approach / Framework | Intensity Factor | Quantity (Capacity) Factor | Kinetic Factor |
|---|---|---|---|
| Theoretical Concept (Frossard et al. (2000)) | Chemical activity of P in the soil solution | Pool of solid-phase P that can replenish the solution | Rate of P transfer from the solid to the liquid phase |
| Traditional / Static Proxy (GRUD) | \(P_{CO_2}\) (Raw Water-soluble Concentration) | \(P_{AAE10}\) (Chelate-extractable Mass) | Not Measured |
| Dynamic / Thermodynamic Proxy (This Study) | Maximum Thermodynamic Activity (\(a_{max}\)) | \(P_{AAE10}\) | Initial Activity Flux (\(J_0 = k \cdot a_{max}\)) |
#|echo: false
library(magick)Linking to ImageMagick 6.9.12.98
Enabled features: fontconfig, freetype, fftw, heic, lcms, pango, raw, webp, x11
Disabled features: cairo, ghostscript, rsvgUsing 4 threadslibrary(pdftools)Using poppler version 24.02.0
This thesis hypothesizes that integrating thermodynamic activity and desorption kinetics provides a superior mechanistic framework for P-management compared to standard empirical testing. Using soils from the long-term Swiss agricultural experiment STYCS, this study aims to validate this framework through the following objectives:
The soil samples for this thesis originate from a set of six long-term field trials in Switzerland, established by Agroscope between 1989 and 1992. The primary objective of these experiments was to validate and re-evaluate Swiss phosphorus (P) fertilization guidelines by assessing long-term crop yield responses to varying P inputs across different pedoclimatic conditions. A detailed description of the experimental design and site characteristics can be found in Hirte et al. (2021).
The experiment was set up as a completely randomized block design with four field replications at each site. The core of the experiment consists of six fixed-plot treatments representing different P fertilization levels, which were applied annually as superphosphate before tillage and sowing. These levels were based on percentages of the officially recommended P inputs: 0% (Zero), 33% (Deficit), 67% (Reduced), 100% (Norm), 133% (Elevated), and 167% (Surplus).
The six experimental sites are located in the main crop-growing regions of Switzerland: Rümlang-Altwi (ALT), Cadenazzo (CAD), Ellighausen (ELL), Grabs (GRA), Oensingen (OEN), and Zurich-Reckenholz (REC). The key soil properties are summarized below.
| Site | Soil Type (WRB) | Clay (%) | Sand (%) | Organic C (g/kg) | pH (H2O) |
|---|---|---|---|---|---|
| ALT | Calcaric Cambisol | 22 | 48 | 21 | 7.9 |
| CAD | Eutric Fluvisol | 8 | 40 | 14 | 6.3 |
| ELL | Eutric Cambisol | 33 | 31 | 23 | 6.6 |
| GRA | Calcaric Fluvisol | 17 | 34 | 16 | 8.3 |
| OEN | Gleyic-calc. Cambisol | 37 | 32 | 24 | 7.1 |
| REC | Eutric Gleysol | 39 | 25 | 27 | 7.4 |
Soil samples for this thesis were collected in the year 2022 from the topsoil layer (0-20 cm).
The laboratory extraction of phosphorus was based on the sequential kinetic principles established by Flossmann & Richter (1982), subsequently modified to capture high-resolution temporal kinetics and corrected to reflect true thermodynamic activity.
To capture the desorption process with high temporal resolution while preserving the native soil matrix, a modified water-extraction protocol was employed:
Although the kinetic extraction utilized deionized water, the immediate dissolution of native soil salts (\(Ca^{2+}, Mg^{2+}, K^+\), etc.) generates a non-zero ionic strength (\(I\)) in the suspension. Because plant roots respond to thermodynamic energy rather than physical mass, the stoichiometric concentration of the extract (\(c_i\)) must be corrected to its true chemical activity (\(a_i\)).
The activity coefficient (\(\gamma\)) for the divalent phosphate ion (\(HPO_4^{2-}\), \(z = -2\)) was calculated for each site-specific matrix using the Davies Equation, which is highly robust for agricultural soil solutions up to \(I \approx 0.1 \text{ M}\):
\[\log_{10}(\gamma) = -0.5 \cdot z^2 \left( \frac{\sqrt{I}}{1+\sqrt{I}} - 0.3 I \right)\]
The true thermodynamic activity (\(a_i\)) of the phosphorus in the extract was then calculated as the product of the raw stoichiometric concentration and the activity coefficient:
\[a_i = c_i \cdot \gamma\]
This correction neutralizes pedochemical variances caused by differing native background salinities and pH levels across the STYCS sites, converting the raw data into a universal thermodynamic Intensity (\(I\)) metric.
All data processing, mathematical modeling, and visualization were conducted using the R programming language (v. 4.3.1) (R Core Team, 2022). Primary packages included nlme (Pinheiro et al., 2022) for non-linear kinetic modeling, lme4 (Bates et al., 2015) for linear mixed-effects modeling, and deSolve and ReacTran for partial differential equation solving.
To derive the kinetic parameters, a non-linear mixed-effects model (nlme) was fitted to the exact solution of the first-order rate equation, utilizing the thermodynamically corrected activities:
\[a(t) = a_{max} \times (1 - e^{-k \times t'})\]
Where \(a(t)\) is the P activity at time \(t\), \(a_{max}\) is the maximum thermodynamic activity of the desorbable P pool (the asymptote), \(k\) is the first-order desorption rate constant, and \(t'\) is an adjusted time (\(t_{min} + 3 \text{ min}\)) accounting for initial dissolution. Sample-specific deviations were modeled as random effects to capture the unique desorption characteristics of each soil.
From these parameters, a novel predictor for agronomic availability, the Initial Activity Flux (\(J_0\)), was defined mathematically as the product of the rate constant and the maximum thermodynamic activity:
\[J_0 = k \cdot a_{max}\]
To deduce the Phosphorus Buffering Capacity (PBC), a thermodynamic Quantity-Intensity (\(Q/I\)) framework was constructed. The \(P_{AAE10}\) extraction was utilized as the standard mass Quantity (\(Q\), \(\text{mg P kg}^{-1}\)), while the Davies-corrected \(a_{max}\) was utilized as the native Intensity (\(I\), \(\mu\text{M}\)).
The relationship between these variables across the varying fertilizer treatments at each site was parameterized using a log-linearized Freundlich Isotherm:
\[Q = K_f \cdot I^{\frac{1}{n}}\]
The true thermodynamic PBC for each specific site and state was then calculated analytically as the first derivative of the fitted Freundlich isotherm with respect to Intensity:
\[\text{PBC} = \frac{dQ}{dI} = K_f \cdot \frac{1}{n} \cdot I^{\left(\frac{1}{n} - 1\right)}\]
To statistically evaluate the predictive power of the chemical parameters across six highly diverse pedoclimatic sites, raw crop yield and P-uptake data required robust normalization to isolate chemical P-limitation from overarching climatic noise.
Standard normalization (dividing by the maximum yield of a given site within the specific experimental year) was rejected, as it artificially masks broad climatic limitations (e.g., drought) and forces the best plot of a poor year to equal 100%, mathematically erasing the biological penalty. Instead, a historical biological baseline strategy was employed.
Using the 40-year STYCS historical dataset, the true biological potential (\(Y_{max}\) and \(Uptake_{max}\)) was calculated for each specific site and crop species by taking the 95th percentile of all fully fertilized (100% and 167% P) treatments. Current experimental yields and uptake values were divided by this historical baseline to calculate Relative Historical Yield (\(Y_{rel\_hist}\)), providing a highly stable metric representing the percentage of the crop’s long-term biological potential achieved.
To test the core hypothesis, two competing sets of linear mixed-effects models (lmer) were constructed to predict \(Y_{rel\_hist}\) and Relative P-Uptake. Both sets shared an identical random effects structure—(1|year) + (1|Site) + (1|Site:block)—to absorb unmeasured pedoclimatic and meteorological variance.
| Abbreviation | Variable | Unit | Description |
|---|---|---|---|
| \(Y_{rel\_hist}\) | Relative Yield (Historic) | unitless | Plot yield normalized by the site/crop-specific 95th percentile of historical well-fertilized yields. |
| \(Uptake_{rel}\) | Relative P Uptake | unitless | Plot P-uptake normalized by historical maximum uptake for the specific site and crop. |
| \(P_{bal}\) | Cumulative P Balance | kg P ha⁻¹ | Long-term net P budget (inputs minus outputs). |
| \(k\) | Rate Constant | min⁻¹ | First-order rate constant of P desorption. |
| \(a_{max}\) | Maximum Activity (Intensity) | µM | Davies-corrected thermodynamic activity asymptote of the water-desorbable P pool. |
| \(J_0\) | Activity Flux | µM min⁻¹ | Initial thermodynamic supply rate, calculated as \(k \times a_{max}\). |
| \(P_{AAE10}\) | Standard Mass (Quantity) | mg P kg⁻¹ | Physical mass of extractable P (AAE10 method), acting as the Q-factor. |
| \(\text{PBC}\) | Phosphorus Buffering Capacity | mg kg⁻¹ µM⁻¹ | The mathematical derivative (\(dQ/dI\)) of the site-specific Freundlich isotherm. |
To mechanistically validate the field-scale statistical findings and visualize the limitations of static buffering capacity, a modified radial transport model based on Barber-Cushman principles was constructed. Standard iterations of the Barber-Cushman model calculate radial diffusion using stoichiometric concentration (\(C_l\)) and a static buffer power (\(b = \Delta Q / \Delta I\)), effectively assuming instantaneous soil-solution equilibrium.
To reflect the kinetic bottlenecks identified in this study, the partial differential equation (PDE) for radial transport in cylindrical coordinates was modified. The driving gradient was shifted to thermodynamic activity (\(a_l\)), and the static buffer power was replaced with a dynamic, first-order kinetic source term (\(R_{desorption}\)):
\[\frac{\partial a_l}{\partial t} = D_e \left( \frac{\partial^2 a_l}{\partial r^2} + \frac{1}{r} \frac{\partial a_l}{\partial r} \right) + \frac{r_0 v_0}{r} \frac{\partial a_l}{\partial r} + R_{desorption}\]
Where \(D_e\) is the effective diffusion coefficient, \(r\) is the radial distance, \(r_0\) is the root radius, and \(v_0\) is the water flux at the root surface. The kinetic source term resupplying the solution was defined using the experimentally derived parameters:
\[R_{desorption} = k \cdot (a_{max} - a_l)\]
Furthermore, the inner boundary condition at the root surface (\(r = r_0\)) was upgraded to reflect a dual-affinity Michaelis-Menten uptake system, accounting for both high-affinity (PHT1) and low-affinity (LATS) transporter kinetics responding directly to chemical activity:
\[J_r = \frac{V_{max1} \cdot (a_l - a_{min})}{K_{m1} + (a_l - a_{min})} + \frac{V_{max2} \cdot (a_l - a_{min})}{K_{m2} + (a_l - a_{min})}\]
The modified PDEs were solved numerically in R utilizing the ReacTran and deSolve packages to generate dynamic spatial profiles of the P-depletion zone over time.
Standard static soil test methods rely on the operational extraction of phosphorus. Strong extractants (e.g., Olsen P, P-AAE10) possess inherently high ionic strengths that dominate the soil matrix, rendering absolute thermodynamic comparisons invalid without complex modeling. In contrast, the \(P_{CO_2}\) method utilizes a very weak extractant: water saturated with carbon dioxide at atmospheric pressure (1 atm) and 25 °C.
To demonstrate that the \(P_{CO_2}\) extractant adopts the native ionic strength of the soil, we calculate the intrinsic ionic strength (\(I\)) of the pure extractant.
According to Henry’s Law, the concentration of dissolved \(CO_2\) (denoted as \(H_2CO_3^*\)) is: \[[H_2CO_3^*] = K_H \times P_{CO_2}\] Given Henry’s constant (\(K_H \approx 0.034 \text{ M atm}^{-1}\)) and \(P_{CO_2} = 1 \text{ atm}\), the molarity of dissolved \(CO_2\) is approximately \(0.034 \text{ M}\).
The primary ionic species generated arise from the first dissociation of carbonic acid: \[H_2CO_3^* \rightleftharpoons H^+ + HCO_3^-\] Using the acid dissociation constant (\(K_{a1} = 4.46 \times 10^{-7}\) at 25 °C): \[K_{a1} = \frac{[H^+][HCO_3^-]}{[H_2CO_3^*]}\] Because \([H^+] \approx [HCO_3^-]\), we can solve for the bicarbonate concentration: \[[HCO_3^-] = \sqrt{4.46 \times 10^{-7} \times 0.034} \approx 1.23 \times 10^{-4} \text{ M}\]
The intrinsic ionic strength (\(I\)) of the extractant is defined as: \[I = \frac{1}{2} \sum (C_i z_i^2)\] \[I = \frac{1}{2} \left( (1.23 \times 10^{-4} \times (+1)^2) + (1.23 \times 10^{-4} \times (-1)^2) \right) \approx 0.00012 \text{ M}\]
This intrinsic ionic strength (\(\approx 0.12 \text{ mM}\)) is infinitesimally small compared to native agricultural soil solutions (which typically range from \(10\) to \(50 \text{ mM}\)). Therefore, the final ionic strength of the \(P_{CO_2}\) extract is overwhelmingly dictated by the dissolved soluble salts of the specific soil being analyzed.
Because the \(P_{CO_2}\) extract adopts the native ionic composition of the soil, reading the orthophosphate concentrations against an \(I=0\) (distilled water) calibration curve introduces systematic error. The true thermodynamic driving force for P-availability is the chemical activity, not the stoichiometric concentration.
To correct for this, an activity coefficient (\(\gamma\)) must be applied. While the Pitzer equations are the gold standard for high-ionic-strength solutions (\(I > 0.5 \text{ M}\), such as the P-AAE10 extractant), the ionic strength of native agricultural soil solutions is typically much lower (\(I < 0.1 \text{ M}\)).
In this low-to-moderate ionic strength regime, the Davies equation (an empirical extension of the Debye-Hückel theory) is highly accurate and does not require complex, ion-specific interaction parameters.
The Davies equation calculates the activity coefficient (\(\gamma\)) as follows: \[\log_{10}(\gamma) = -A z^2 \left( \frac{\sqrt{I}}{1 + \sqrt{I}} - 0.3 I \right)\] Where: * \(A\) is a temperature-dependent constant (\(\approx 0.509 \text{ L}^{0.5} \text{mol}^{-0.5}\) for water at 25 °C). * \(z\) is the charge of the ion in question (for orthophosphate, \(z\) depends on pH, predominantly \(-1\) for \(H_2PO_4^-\) and \(-2\) for \(HPO_4^{2-}\)). * \(I\) is the total ionic strength of the extracted solution, calculated from the measured concentrations of major cations (\(Ca^{2+}, Mg^{2+}, K^+\)).
To assess the theoretical P-Intensity (\(I\)) and formulate the Phosphorus Buffering Capacity (PBC), the operationally defined \(P_{CO_2}\) concentrations are converted to thermodynamic activities using the following sequential steps:
By standardizing the measurements to their true thermodynamic activity, we remove the analytical bias introduced by varying soil salinity and pedoclimatic history, allowing for an absolute comparison of P-Intensity across all STYCS experimental sites. ## A.4 Supplementary Materials and Reproducibility
This Master’s thesis was produced using a fully reproducible workflow with the Quarto publishing system. All data, R scripts, and analytical notebooks used to generate the figures, tables, and results presented in this work are openly available.
The complete project can be cloned from the author’s GitHub repository, which contains the raw data, the R code for the kinetic and statistical models, and the Quarto source files.
GitHub Repository URL: [https://github.com/ptochosgeorgos/Article-P-desorption-kinetics]
A rendered version of the full project, including the analytical notebooks that document the development process, is also available as a GitHub Pages website at the following URL:
GitHub Pages Site URL: [https://ptochosgeorgos.github.io/Article-P-desorption-kinetics/]
This approach ensures full transparency and allows for the complete replication of the findings presented in this thesis.