In silico calculation of soil pH by SCEPTER v1.0

One of the soil properties most commonly measured to describe agronomic and biogeochemical conditions of soils is soil pH. Soil pH measures the concentration of exchangeable H⁺ that resides in bulk soil samples using extractants in the laboratory and thus differs from porewater pH, which we define here as an in situ measure of porewater H⁺ concentration in soil/weathering profiles. The difference between the two pH measurements is often not fully known for a given system but could lead to a misunderstanding of soil conditions if the two measurements are directly compared. Agricultural soils are one of the targeted loci for the application of enhanced weathering (EW), a technique aimed at counteracting increasing anthropogenic carbon dioxide from burning fossil fuels. An increase in pH is thought to be one of the key advantages of EW, given that the process can mitigate soil acidification and increase crop yields. As a result, fully evaluating the biogeochemical and agronomic consequences of EW approaches requires accurate simulation of both soil pH (pH_s) and porewater pH (pH_pw). This paper presents an updated version of the reactive transport code SCEPTER (Soil Cycles of Elements simulator for Predicting TERrestrial regulation of greenhouse gases), which enables simulation of bulk soil pH measurements in the laboratory, in addition to porewater pH, as measured in the field along with a more comprehensive representation of cation exchange with solid-phase constituents of bulk soil. We first describe the implementation of cation exchange in the SCEPTER model, then introduce conceptual modeling frameworks enabling the calculation of bulk pH_s. The validity of the model is examined through comparison of model results with soil pH measurements from mesocosm experiments on maize production with crushed basalt amendments. Finally, illustrative example simulations are shown, demonstrating that a difference between pH_s and pH_pw can lead to significantly different estimates of soil alkalinization and carbon capture by EW for a given targeted pH in cropland systems.

1 Introduction

Continuous harvesting and excess use of nitrogen fertilizers commonly lead to the acidification of agricultural soils, which may lead to soil degradation and food insecurity over the coming century (Kopittke et al., 2019). An addition of alkalinity to soils – traditionally through liming, the application of ground, relatively soluble (mostly carbonate) rock/mineral powder to soils (e.g., McLean, 1983; Thomas Sims, 1996; Rengel, 2003; Goulding, 2016) – is a widely utilized remedy to manage soil pH and stabilize crop yields. The addition of alkalinity to soil (including agricultural liming) has recently attracted attention because it can also sequester atmospheric CO₂ (e.g., Hamilton et al., 2007; Swoboda et al., 2022), an action that is urgently needed to help meet the climate targets delineated by the Intergovernmental Panel on Climate Change (IPCC, 2006, 2018). Indeed, enhanced weathering (EW) – the application of finely ground carbonate/silicate rock powder to soils – is one of a number of suggested schemes for actively removing anthropogenic CO₂ from the atmosphere at scale (e.g., Rau et al., 2007; Köhler et al., 2010; Taylor et al., 2016; Beerling et al., 2020; Vakilifard et al., 2021; Swoboda et al., 2022; Zhang et al., 2022; Kanzaki et al., 2023). In particular, applying basalt rock powder onto croplands/hinterlands has been suggested to be one of the most scalable, safe, and economically promising CO₂ removal schemes given the relatively low toxicity in basalt leachates, the sustainable availability of basalt rocks, and a range of potential co-benefits (e.g., Strefler et al., 2018; Beerling et al., 2020; Goll et al., 2021).

The pH change induced by the addition of basalt powder is central to the EW scheme because the resultant pH (reflecting, e.g., the soil buffer capacity, local climate, and particle size distribution of the milled rock) must be optimal for crop growth (e.g., Fernández and Hoeft, 2009), and the application rate of basalt feedstock and resultant carbon capture thus scale with the magnitude of the desired pH increase (e.g., Kelland et al., 2020; Kantzas et al., 2022; Zhang et al., 2022; Dietzen and Rosing, 2023). However, the interpretation of pH in soil is not always straightforward because two different types of pH measurements may potentially be regarded as a pH reference for evaluating soil acidity. One is referred to as soil pH – defined here as pH_s – which measures H⁺ residing in bulk soil samples that is measured in the laboratory as the pH of liquid extractants (deionized water or $\mathrm{KCl}/{\mathrm{CaCl}}_{\mathrm{2}}$ solution) of bulk soil samples taken from the field. The other is porewater pH – defined here as pH_pw – which measures in situ H⁺ concentrations in porewater flowing through or remaining in the soil/weathering profiles (e.g., Geibe et al., 2006; Steiner et al., 2018). In agricultural/agronomic situations, it is most common to measure pH_s (e.g., Thomas, 1996), while models that simulate biogeochemical reactions and transport within soils, including the dissolution of basalt during EW, typically calculate pH_pw (e.g., Kelland et al., 2020; Kanzaki et al., 2022). Potential differences between these distinct tracers of soil acidity are poorly explored, and in many cases the heterogeneous continuum that exists between dissolved H⁺ in pore fluids and exchangeable H⁺ on soil cation exchange sites is not discussed (cf. Nielsen et al., 2017).

Here, we present a newly developed numerical scheme in an attempt to fill this technical and knowledge gap and to develop a more mechanistic understanding of the difference between porewater pH and bulk soil pH. A numerical reactive transport model – SCEPTER (Soil Cycles of Elements simulator for Predicting TERrestrial regulation of greenhouse gases; Kanzaki et al., 2022) – has been updated to enable simulations of soil pH (pH_s) along with porewater pH (pH_pw). We first present the essential updates to the SCEPTER code (Sect. 2.1) and then describe potential modeling frameworks for simulating soil pH with the updated version of the model (Sect. 2.2). Then, the validity of the model is examined through comparison between simulated and observed soil pH in a set of mesocosm experiments amending a natural soil/maize system with crushed basalt (Sect. 3). We then discuss the implications of the difference between porewater and soil pH for EW and the associated impacts on soil acidity by showing example simulations in which basalt feedstock is added to cropland soil using either pH_s or pH_pw as a target pH for EW deployment (Sect. 4). Finally, we provide a summary of conclusions and touch briefly on future directions for model development (Sect. 5).

2 Model description

The SCEPTER model simulates the reactions and transport of solid, aqueous, and gas species within soil, including the dissolution/precipitation of minerals, three-phase biogeochemical reaction, bio-mixing and uplift/erosion of solid phases, advective and diffusive transport of aqueous species, and gaseous diffusion (Kanzaki et al., 2022). The model is developed for simulating not only natural weathering processes, but also EW with its specific features that allow for the explicit bio-mixing of soil, including tilling by farmers; addition of solid materials on the topsoil; and tracking of particle size distributions, which facilitates surface area calculations for individual solid species. This updated version of the code (v1.0) adds several new functions/options to the previously published version (v0.9; Kanzaki et al., 2022). Among them, the implementation of cation exchange is essential to the simulation of soil pH as the uptake of cations by solid-phase exchangers is a determining factor of the exchangeable acidity and nutrient cycling in soils. We first describe the implementation of cation exchange in SCEPTER (Sect. 2.1) and then the frameworks for the simulation of soil pH using the current version of the code (Sect. 2.2). All symbols used in this study and their definitions are summarized in Appendix A.

Table 1 Default cation exchange capacity of solid species^a.

^a Those solid species that are not listed here are assumed to have zero cation exchange capacity. ^b SOM: soil organic matter. ^c (1) Beerling et al. (2020). (2) Parfitt et al. (1995).

2.1 Cation exchange in SCEPTER

The current version of SCEPTER allows for cation exchange involving H⁺, Na⁺, K⁺, Ca²⁺, Mg²⁺, and Al³⁺ on any solid species specified by the user (on clay minerals and organic matter by default; Table 1). Cation exchange reactions are assumed to be in equilibrium, and their fundamental reactions can be written as reactions among the surfaces of solid-phase exchangers and the cations:

$$\begin{array}{}\text{(1)}& Z_{\mathit{\varsigma }}\mathrm{X}(\mathit{\theta }{)}^{-}+{\mathit{\varsigma }}^{Z_{\mathit{\varsigma }}+}\iff \mathit{\varsigma }\text{-}\mathrm{X}(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}},\end{array}$$

where Z_ς is the valence number of cation ς; X(θ)⁻ denotes exchangeable surface sites of solid-phase exchanger θ; and $\mathit{\varsigma }\text{-}\mathrm{X}(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}$ represents the cation ς adsorbed onto exchangeable sites of θ. Equation (1) should be regarded as a half reaction because the surface fraction of X(θ)⁻ must be very small compared to the surface sites where net local charge is zero because of adsorption under natural conditions (Appelo, 1994). Physically relevant net cation exchange can then be written as a combination of Eq. (1) for a given cation and Eq. (1) for a reference/competing cation so that the combined reaction equation does not have X(θ)⁻. As a reference cation, Na⁺ and Cs⁺ have been considered (e.g., Appelo, 1994; Steefel et al., 2003; Steefel, 2009). Here, we use H⁺ as a reference competing agent, with the net exchange reaction given as

$$\begin{array}{}\text{(2)}& (\mathrm{1}/Z_{\mathit{\varsigma }}){\mathit{\varsigma }}^{Z_{\mathit{\varsigma }}+}+\mathrm{H}\text{-}\mathrm{X}\left(\mathit{\theta }\right)\iff (\mathrm{1}/Z_{\mathit{\varsigma }})\mathit{\varsigma }\text{-}\mathrm{X}(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}+{\mathrm{H}}^{+}.\end{array}$$

The equilibrium constant for Eq. (2) can be defined as follows:

$$\begin{array}{}\text{(3)}& K_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}^{\prime }={\displaystyle \frac{f^{\mathrm{1}/Z_{\mathit{\varsigma }}}\left(\mathit{\varsigma }\text{-}\mathrm{X}\right(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}\left)\right[{\mathrm{H}}^{+}]}{f\left(\mathrm{H}\text{-}\mathrm{X}\right(\mathit{\theta }\left)\right)[{\mathit{\varsigma }}^{Z_{\mathit{\varsigma }}+}{]}^{\mathrm{1}/Z_{\mathit{\varsigma }}}}},\end{array}$$

where f(i) denotes the charge-equivalent fraction of surface species i, and [j] represents the concentration of aqueous species j (mol L⁻¹). The apparent equilibrium constant $K_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}^{\prime }$ can vary as a result of surface fraction X(θ)⁻, and we adopt the formulation by Appelo (1994):

$$\begin{array}{}\text{(4)}& K_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}^{\prime }={\mathit{\eta }}_{\mathrm{H},\mathit{\theta }}{K}_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}.\end{array}$$

Here, $K_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}$ is the intrinsic equilibrium constant and η_H,θ is formulated as a function of 1−f(H-X(θ)), assuming that f(X(θ)⁻) is proportional to 1−f(H-X(θ)) (Appelo, 1994):

$$\begin{array}{}\text{(5)}& \log {\mathit{\eta }}_{\mathrm{H},\mathit{\theta }}=-{\mathit{\alpha }}_{\mathit{\theta }}\mathit{\{}\mathrm{1}-f(\mathrm{H}\text{-}\mathrm{X}\left(\mathit{\theta }\right)\left)\mathit{\right\}},\end{array}$$

where α_θ is assumed to be 3.4 by default.

The solution for the fraction of surface species can be obtained by considering the mass balance at the exchangeable sites for each exchanger:

$$\begin{array}{}\text{(6)}& {\mathit{CEC}}_{\mathit{\theta }}=\sum _{\mathit{\varsigma }}{Z}_{\mathit{\varsigma }}\langle \mathit{\varsigma }\text{-}\mathrm{X}(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}\rangle ,\end{array}$$

where CEC_θ is the cation exchange capacity of exchanger θ (eq g⁻¹) and 〈i〉 is the concentration of surface species i (mol g⁻¹). By definition, we have

$$\begin{array}{}\text{(7)}& f\left(\mathit{\varsigma }\text{-}\mathrm{X}\right(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}})\equiv {\displaystyle \frac{Z_{\mathit{\varsigma }}\langle \mathit{\varsigma }\text{-}\mathrm{X}(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}\rangle }{{\mathit{CEC}}_{\mathit{\theta }}}}.\end{array}$$

Therefore, Eq. (6) can be alternatively written as

$$\begin{array}{}\text{(8)}& \mathrm{1}=\sum _{\mathit{\varsigma }}f\left(\mathit{\varsigma }\text{-}\mathrm{X}\right(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}).\end{array}$$

Further, with Eqs. (3) and (4), Eq. (8) can be transformed into

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{1}& =f\left(\mathrm{H}\text{-}\mathrm{X}\right(\mathit{\theta }\left)\right)\\ & +\sum _{\mathit{\varsigma }\ne \mathrm{H}}{\left({\displaystyle \frac{{\mathit{\eta }}_{\mathrm{H},\mathit{\theta }}{K}_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}f\left(\mathrm{H}\text{-}\mathrm{X}\right(\mathit{\theta }\left)\right)}{\left[{\mathrm{H}}^{+}\right]}}\right)}^{Z_{\mathit{\varsigma }}}\left[{\mathit{\varsigma }}^{Z_{\mathit{\varsigma }}+}\right].\end{array}\end{array}$$

Equation (9) is solved for f(H-X(θ)) once given a porewater chemistry and thermodynamic constants for exchange reactions (Table 2). Once f(H-X(θ)) is obtained, fractions of all surface species can be calculated using Eqs. (3)–(5).

Table 2 Default thermodynamic data of cation exchange^a.

^a The same set of thermodynamic data is assumed for any solid-phase exchanger. Therefore, the notation of solid-phase θ used in Sect. 2 is dropped in this table. ^b (1) From the modeled value at zero f(H-X) in Appelo (1994). (2) Calculated from log $K_{\mathit{\varsigma }\backslash \mathrm{Na}}=\mathrm{1.1}$, 0.507, and 0.665 for ς=K, Mg, and Ca, respectively, from Appelo (1994). (3) Calculated from log $K_{\mathrm{Al}\backslash \mathrm{Na}}=\mathrm{0.41}$ from phreeqc.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013).

In the previous version of SCEPTER, the key variables for tracking aqueous species are the total concentrations for individual dissolved elements (denoted as c_ς for dissolved element ς). In the updated model, the tracked independent variables have been changed to the concentrations of free dissolved species (except for Si, for which the H₄SiO₄ concentration is tracked), denoted as $c_{\mathit{\varsigma }}^{\mathrm{1}}$. These c_ς and $c_{\mathit{\varsigma }}^{\mathrm{1}}$ are related to one another by the following equation (Kanzaki et al., 2022):

$$\begin{array}{}\text{(10)}& c_{\mathit{\varsigma }}=c_{\mathit{\varsigma }}^{\mathrm{1}}+c_{\mathit{\varsigma }}^{\mathrm{1}}\sum _{i=\mathrm{2}}^{n_{\mathit{\varsigma }}}{K}_{\mathit{\varsigma },i}\left[{\mathrm{H}}^{+}{]}^{{\mathit{\gamma }}_{\mathit{\varsigma },i,\mathrm{p}}}\prod _{{\mathit{\varsigma }}^{\prime }\ne \mathit{\varsigma }}^{n_{\mathrm{aq}}}\right(c_{{\mathit{\varsigma }}^{\prime }}^{\mathrm{1}}{)}^{{\mathit{\gamma }}_{\mathit{\varsigma },i,{\mathit{\varsigma }}^{\prime }}}\prod _{\mathit{\epsilon }}^{n_{\mathrm{gas}}}(p_{\mathit{\epsilon }}{)}^{{\mathit{\gamma }}_{\mathit{\varsigma },i,\mathit{\epsilon }}},\end{array}$$

where the second term on the right-hand side is the sum of the concentrations of dissolved element ς other than $c_{\mathit{\varsigma }}^{\mathrm{1}}$, denoted as the ith species of dissolved element ς where i≠1, with K_ς,i being the thermodynamic constant for production of ith aqueous species of dissolved element ς; ${\mathit{\gamma }}_{\mathit{\varsigma },i,\mathrm{p}}$, ${\mathit{\gamma }}_{\mathit{\varsigma },i,{\mathit{\varsigma }}^{\prime }}$, and ${\mathit{\gamma }}_{\mathit{\varsigma },i,\mathit{\epsilon }}$ the stoichiometry of H⁺, dissolved element ς^′, and gas species ε, respectively, in the reaction that produces ith aqueous species of ς; p_ε the partial pressure (atm) of gas species ε; and n_aq and n_gas the total numbers of independent aqueous and gas species, respectively (see Kanzaki et al., 2022, for more details). This modification of tracked independent variables (from c_ς to $c_{\mathit{\varsigma }}^{\mathrm{1}}$) facilitates our implementation of cation exchange.

In accordance with the implementation of cation exchange and modification of independent variables to track for aqueous species described above, the governing equation for aqueous species has been updated to

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\varphi }\mathit{\sigma }\mathrm{\ell }{\mathit{\beta }}_{\mathit{\varsigma }}^{\mathrm{aq}}{c}_{\mathit{\varsigma }}^{\mathrm{1}}}{\partial t}}& +{\displaystyle \frac{\partial B_{\mathit{\varsigma }}^{\mathrm{ads}}{c}_{\mathit{\varsigma }}^{\mathrm{1}}}{\partial t}}=-{\displaystyle \frac{\partial \mathit{\varphi }\mathit{\sigma }\mathrm{\ell }v{\mathit{\beta }}_{\mathit{\varsigma }}^{\mathrm{aq}}{c}_{\mathit{\varsigma }}^{\mathrm{1}}}{\partial z}}\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\varphi }\mathit{\sigma }\mathrm{\ell }{\mathit{\tau }}_{\mathrm{aq}}{D}_{\mathit{\varsigma }}{\displaystyle \frac{\partial {\mathit{\beta }}_{\mathit{\varsigma }}^{\mathrm{aq}}{c}_{\mathit{\varsigma }}^{\mathrm{1}}}{\partial z}}\right)\\ & +\sum _{\mathit{\theta }}^{n_{\mathrm{sld}}}{\mathit{\gamma }}_{\mathit{\theta },\mathit{\varsigma }}{R}_{\mathit{\theta }}+\sum _{\mathit{\kappa }}^{n_{\mathrm{xrxn}}}{\mathit{\gamma }}_{\mathit{\kappa },\mathit{\varsigma }}{R}_{\mathit{\kappa }}\\ & +{\displaystyle \frac{\partial w{B}_{\mathit{\varsigma }}^{\mathrm{ads}}{c}_{\mathit{\varsigma }}^{\mathrm{1}}}{\partial z}}-B_{\mathit{\varsigma }}^{\mathrm{ads}}{c}_{\mathit{\varsigma }}^{\mathrm{1}}\underset{\mathrm{0}}{\overset{z_{\mathrm{ml}}}{\int }}{E}_{\mathit{\theta }}(z,z^{\prime })\mathrm{d}{z}^{\prime }\\ & +\underset{\mathrm{0}}{\overset{z_{\mathrm{ml}}}{\int }}{B}_{\mathit{\varsigma }}^{\mathrm{ads}}{c}_{\mathit{\varsigma }}^{\mathrm{1}}\left(z^{\prime }\right)E_{\mathit{\theta }}(z^{\prime },z)\mathrm{d}{z}^{\prime }.\end{array}\end{array}$$

The first and second terms on the left-hand side of Eq. (11) denote the time rate of change of dissolved and adsorbed forms of ς, respectively, with ${\mathit{\beta }}_{\mathit{\varsigma }}^{\mathrm{aq}}$ and $B_{\mathit{\varsigma }}^{\mathrm{ads}}$ (m⁻³ L) defined as the factors to convert $c_{\mathit{\varsigma }}^{\mathrm{1}}$ to c_ς and to the total concentration of element ς adsorbed onto solid phases, respectively. The first and second terms on the right-hand side show the advective and diffusive transport rates of dissolved forms of ς, respectively; the third and fourth terms the net supply of ς through dissolution/precipitation of solid phases and extra biogeochemical reactions, respectively; and the rest of the terms the advective transport (fifth term) and bio-mixing (sixth and seventh terms) rates of adsorbed forms of ς along with solid-phase exchangers. The parameters to formulate the individual terms, reactions, and transport in Eq. (11) described above are tabulated in Appendix A. See Kanzaki et al. (2022) for further details on the reactions and transport schemes implemented in SCEPTER. The values of ${\mathit{\beta }}_{\mathit{\varsigma }}^{\mathrm{aq}}$ and $B_{\mathit{\varsigma }}^{\mathrm{ads}}$ can be calculated from Eqs. (10) and (3)–(9), respectively:

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}{\mathit{\beta }}_{\mathit{\varsigma }}^{\mathrm{aq}}& \equiv {\displaystyle \frac{c_{\mathit{\varsigma }}}{c_{\mathit{\varsigma }}^{\mathrm{1}}}}=\mathrm{1}\\ & +\sum _{i=\mathrm{2}}^{n_{\mathit{\varsigma }}}{K}_{\mathit{\varsigma },i}\left[{\mathrm{H}}^{+}{]}^{{\mathit{\gamma }}_{\mathit{\varsigma },i,\mathrm{p}}}\prod _{{\mathit{\varsigma }}^{\prime }\ne \mathit{\varsigma }}^{n_{\mathrm{aq}}}\right(c_{{\mathit{\varsigma }}^{\prime }}^{\mathrm{1}}{)}^{{\mathit{\gamma }}_{\mathit{\varsigma },i,{\mathit{\varsigma }}^{\prime }}}\prod _{\mathit{\epsilon }}^{n_{\mathrm{gas}}}(p_{\mathit{\epsilon }}{)}^{{\mathit{\gamma }}_{\mathit{\varsigma },i,\mathit{\epsilon }}},\end{array}\text{(13)}& {\displaystyle }{B}_{\mathit{\varsigma }}^{\mathrm{ads}}\equiv \left\{\begin{array}{cc}\frac{\mathrm{1}}{c_{\mathit{\varsigma }}^{\mathrm{1}}}\sum _{\mathit{\theta }}{m}_{\mathit{\theta }}{M}_{\mathit{\theta }}\langle \mathit{\varsigma }\text{-}\mathrm{X}(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}\rangle & \\ =\sum _{\mathit{\theta }}\frac{m_{\mathit{\theta }}{M}_{\mathit{\theta }}{\mathit{CEC}}_{\mathit{\theta }}}{Z_{\mathit{\varsigma }}}{\left(\frac{{\mathit{\eta }}_{\mathrm{H},\mathit{\theta }}{K}_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}f\left(\mathrm{H}\text{-}\mathrm{X}\right(\mathit{\theta }\left)\right)}{\left[{\mathrm{H}}^{+}\right]}\right)}^{Z_{\mathit{\varsigma }}}& (\mathit{\varsigma }\in \mathit{\{}\mathrm{Na},\mathrm{K},\mathrm{Ca},\mathrm{Mg},\mathrm{Al}\mathit{\left\}}\right)\\ \mathrm{0}& \text{(else)}\end{array}\right.,\end{array}$$

where m_θ and M_θ are the concentration (mol m⁻³) and molar weight (g mol⁻¹) of solid species θ, respectively.

The updated version of the governing equation for aqueous species (Eq. 11) is solved together with those for solid and gaseous species as described by Kanzaki et al. (2022), except that the calculation of surface speciation via cation exchange is additionally included during each update of porewater pH and aqueous speciation. The default capacities and thermodynamic constants of cation exchange are tabulated in Tables 1 and 2, respectively. Cation exchange can be switched on and off by specifying so in the switches.in input file. One can also modify the cation exchange parameters for any solid species using another input file, cec.in; for example, it is possible to assign different cation exchange parameters to different classes of organic matter that differ from the default values in Tables 1 and 2. Instructions for running example simulations from this paper are given in the “Code and data availability” section.

2.2 Soil pH simulation by SCEPTER

The in silico calculation of bulk soil pH (pH_s) imitates the procedure to measure soil pH in the laboratory: sampling bulk soils, mixing them with an extractant solution (e.g., deionized water or $\mathrm{KCl}/{\mathrm{CaCl}}_{\mathrm{2}}$ solution) at a given soil-to-solution ratio (e.g., 1:1 or 1:5 g mL⁻¹), bringing the mixtures to a short-term equilibrium, and measuring extractant solution pH (e.g., McLean, 1983; Thomas, 1996; Jones, 1999; Kissel and Sonon, 2008). Soil buffer pH – a measure of the resistance of bulk soil to a pH change – can be calculated in silico using the same procedure but with a specified buffer solution (e.g., Thomas Sims, 1996; Sikora, 2006) instead of the extractant solutions implemented for measuring bulk agronomic pH. Our procedure for calculating soil (buffer) pH can be summarized as follows:

• 1.

  A field simulation is run, which can be fed by field observations.

• 2.

  Data from the field run are retrieved at a given model field depth and/or averaged over a given depth interval, including output for

  • a.

    concentrations and cation exchange properties (e.g., Tables 1 and 2) of unextractable solid phases (e.g., silicates)

  • b.

    concentrations of exchangeable (i.e., dissolved plus adsorbed) cations and anions

  • c.

    concentrations of cations and anions in extractable solid phases (e.g., salts).

• 3.

  Boundary conditions for a laboratory simulation are determined based on Step 2 in order to realize a hypothetical laboratory beaker/flask, where a bulk soil sample and an extractant solution (deionized water or electrolyte solution) are mixed homogeneously at a given soil-to-solution ratio.

  • a.

    Concentrations of unextractable solid species obtained in Step 2 are given as the initial/boundary concentrations in an input file (parentrock.in) for the laboratory run. Those solid species are not allowed to dissolve/precipitate in the laboratory run because of the short duration for soil pH measurements (e.g., Thomas, 1996), realized by setting their rate constants at zero in an input file (kinspc.in). Meanwhile, cation exchange properties of the unextractable solid species are assumed to be the same as those in the field run (these can be specified in the corresponding input file, cec.in).

  • b.

    Exchangeable/extractable cations and anions are added to the calculation domain of the laboratory beaker/flask as an appropriate combination of oxides and salts whose complete dissolution is allowed (Table 3). Note that dissolved inorganic carbon (DIC) is an exception and is instead added as the most labile class of organic matter (Table 3) so that carbon can be added without additional cations (compare, e.g., carbonates). The amount of solid species added is calculated as $z_{\mathrm{lab}}(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{lab}})C_{\mathit{\varsigma }}{M}_{\mathit{\theta }}/{\mathit{\gamma }}_{\mathit{\theta },\mathit{\varsigma }}$ (g m⁻²), where z_lab (m) is the depth of the beaker/flask filled with the mixture of a soil sample and a solution; ϕ_lab is the volume ratio of fluid against solid plus fluid phases calculated as ${\mathit{\varphi }}_{\mathrm{lab}}={\mathit{\psi }}^{-\mathrm{1}}({\mathit{\rho }}^{-\mathrm{1}}+{\mathit{\psi }}^{-\mathrm{1}}{)}^{-\mathrm{1}}$ with the soil-to-solution ratio used in the laboratory (ψ; g cm⁻³) and bulk soil particle density (ρ; g cm⁻³) observed in the in silico field; C_ς is the concentration of exchangeable/extractable cation/anion ς (mol m⁻³); M_θ is the molar weight of the added solid species, θ (g mol⁻¹); and γ_θ,ς is the mole of ς contained in 1 mol of θ. When soil pH is measured in the mixture of the bulk soil sample and an electrolyte solution, the corresponding salt is additionally added in the amount of $z_{\mathrm{lab}}{\mathit{\varphi }}_{\mathrm{lab}}\mathrm{\ell }{c}_{\mathrm{\Theta }}{M}_{\mathit{\theta }}/{\mathit{\gamma }}_{\mathit{\theta },\mathrm{\Theta }}$ (g m⁻²), where c_Θ and γ_θ,Θ are the solution concentration (mol L⁻¹) of electrolyte Θ and the mole of electrolyte Θ in 1 mol of salt θ, respectively (e.g., c_Θ=0.01 mol L⁻¹ and ${\mathit{\gamma }}_{\mathit{\theta },\mathrm{\Theta }}=\mathrm{1}$ if θ=CaCl₂; else ${\mathit{\gamma }}_{\mathit{\theta },\mathrm{\Theta }}=\mathrm{0}$, when Θ=CaCl₂). When simulating soil buffer pH, the added salt corresponding to the electrolyte described above must be replaced by a series of solid phases corresponding to solute ingredients, according to the recipe of the buffer solution (e.g., Table 4 for a buffer solution by Sikora, 2006), enabling at the same time the tracking of corresponding aqueous species with relevant aqueous diffusion coefficients and association/dissociation thermodynamics (e.g., Tables 5 and 6, respectively, for the Sikora buffer solution). These constituents are added to the beaker/flask only once at the beginning of a laboratory simulation.

  • c.

    The beaker/flask domain of the laboratory simulation is assumed to be a closed system for solid, aqueous, and gaseous species, except for the addition of solid/salt phases at the beginning of the run described in Step 3b above; i.e., there is no advective transport for solid, aqueous, and gaseous phases and no diffusive in- and outfluxes of aqueous and gaseous species at the boundaries (specified in input files frame.in and switches.in).

• 4.

  The laboratory simulation is run long enough to achieve equilibrium.

Figure 1 schematically illustrates the procedure described above, from running a field simulation and sampling data from the in silico field to the soil pH measurement in the laboratory. As implied by the schematic (e.g., compare aqueous compositions illustrated for porewater in Step 2 and extractant solution in Step 4 of Fig. 1), porewater and soil pH can differ. Indeed, under the conditions considered in our analysis, a significant offset between pH_s and pH_pw is confirmed to be a general phenomenon, as discussed below. In the next section, we discuss the validity of our approach toward simulating pH_s with SCEPTER using observed soil and porewater pH data from a mesocosm experiment along with other observed soil chemical characteristics. After examining the validity of the model (Sect. 3), we present examples of the model application to EW and discuss how the difference between porewater and soil pH can potentially lead to significant differences in the prediction of the amount of basalt feedstock required to achieve a given agronomic target pH in agricultural soils (Sect. 4).

Table 3 Solid species to be dissolved in laboratory simulations^a.

^a Thermodynamic constants for solid species θ (K_θ) are calculated as $K_{\mathit{\theta }}=K_{\mathit{\theta }}^{\mathrm{ref}}$exp($-\mathrm{\Delta }{H}_{\mathit{\theta }}(T^{-\mathrm{1}}-{\mathrm{298}}^{-\mathrm{1}}){\mathcal{R}}^{-\mathrm{1}}$), where T is the temperature in kelvins and ℛ is the gas constant in units of kJ mol⁻¹ K⁻¹ ($\mathcal{R}=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹). Solid species listed here are assumed to have decomposition rates that are represented by a short turnover time (≤1 year) and do not depend on surface areas but their concentrations (see Kanzaki et al., 2022). The variation in kinetic constants does not affect the soil pH simulations as long as they are run long enough to attain equilibrium.
^b Units change with x depending on the solid species.
^c (1) M_θ and V_θ from Robie et al. (1978). (2) $K_{\mathit{\theta }}^{\mathrm{ref}}$ and ΔH_θ from llnl.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013). (3) $K_{\mathit{\theta }}^{\mathrm{ref}}$ and ΔH_θ are assumed to be the same as those for amorphous Si. (4) $K_{\mathit{\theta }}^{\mathrm{ref}}$ and ΔH_θ from minteq.v4.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013). (5) Assumed to be undersaturated unconditionally. M_θ is calculated from the chemical formula, and V_θ is based on M_θ assuming densities of 2.13, 1.725, 1.5, 1.124, 1.23, 0.56, and 1.27 g cm⁻³ for NaOH, NH₄NO₃, SOM Class 1, triethanolamine, imidazole, MES (2-(N-morpholino)ethanesulfonic acid) monohydrate, and acetic acid, respectively.

Table 4 Sikora buffer composition^∗.

^∗ From Sikora (2006), except that the NaOH concentration is modified so that the mixture of the Sikora buffer with deionized water at a 1:1 volume ratio has a pH of 7.5.

Table 5 Diffusion coefficients for aqueous species in the Sikora buffer.

^∗ (1) The diffusion coefficient (m² yr⁻¹) is calculated as D=0.4415 (${\mathit{\mu }}_{\mathrm{w}}^{-\mathrm{1.1}}{a}^{\mathrm{0.6}}$)⁻¹, where μ_w is the water viscosity (mPa s) and a is the molar volume of solute (cm³ mol⁻¹) (Othmer and Thakar, 1953; La-Scalea et al., 2005). The water viscosity μ_w is calculated as ${\mathit{\mu }}_{\mathrm{w}}=\mathrm{0.024152}\exp \left(\mathrm{4.7428}\right(T-\mathrm{139.86}{)}^{-\mathrm{1}}{\mathcal{R}}^{-\mathrm{1}})$, where $\mathcal{R}=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹ and T is temperature in kelvins, according to Likhachev (2003). (2) a from La-Scalea et al. (2005). (3) a is assumed to be equivalent to that of MES monohydrate. (4) The diffusion coefficient (m² yr⁻¹) is calculated as $D=a\times \exp (-b(T^{-\mathrm{1}}-{\mathrm{288}}^{-\mathrm{1}}\left){\mathcal{R}}^{-\mathrm{1}}\right)$ where a is the pre-exponential factor (m² yr⁻¹) and b is the apparent activation energy (kJ mol⁻¹). (5) a and b from Schulz and Zabel (2006). (6) a and b from Li and Gregory (1974).

Table 6 Thermodynamic data for aqueous species in the Sikora buffer^a.

^a The thermodynamic constant (K_aq) is calculated as $K_{\mathrm{aq}}=K_{\mathrm{aq}}^{\mathrm{ref}}$exp($-\mathrm{\Delta }{H}_{\mathrm{aq}}(T^{-\mathrm{1}}-$298${}^{-\mathrm{1}}\left){\mathcal{R}}^{-\mathrm{1}}\right)$, where $\mathcal{R}=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹ and T is temperature in kelvins.
^b TEA: triethanolamine, TEA(H)⁺: H⁺-associated triethanolamine, IM: imidazole, IM(H)⁺: H⁺-associated imidazole, MES: 2-(N-morpholino)ethanesulfonic acid, MES(−H)⁻: H⁺-dissociated MES, AcO⁻ : acetate anion, AcOH: acetic acid.
^c Units change with x depending on reaction.
^d (1) $K_{\mathrm{aq}}^{\mathrm{ref}}$ from Sikora (2006) and ΔH_aq from Goldberg et al. (2002). (2) From llnl.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013).

Figure 1 Schematic of the soil pH calculation procedure. After a field simulation is run to represent a specific field soil (Step 1), in silico field data are obtained (Step 2) for the concentrations of solid phases (left), adsorbed cations (middle), and dissolved cations and anions (right). In Step 3, sampled in silico field data are converted to input data for a laboratory simulation in which extractable/exchangeable cations/anions are converted to a combination of salt/oxide phases to be added to the laboratory beaker/flask, with additional phases depending on the extractant (or buffer) solution. In Step 4, these added phases are dissolved in the laboratory beaker/flask to reach equilibrium, after which the calculated solution pH corresponds to the soil pH (pH_s) of the in silico field soil in Steps 1 and 2.

Figure 2 Comparison of cation exchange simulations between SCEPTER (v1.0) and PHREEQC (v3.0) with the thermodynamic dataset Tipping_Hurley.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013). (a) Solution of 1 mM Na, 0.2 mM K, and 1 to 15 mM NO₃ in equilibrium with a 1.1 meq L⁻¹ cation exchanger. (b) A 0.6 mM CaCl₂ solution replacing a solution of 1 mM Na, 0.2 mM K and 1.2 mM NO₃, initially equilibrated with a 1.1 meq L⁻¹ cation exchanger homogeneously distributed along the soil column, through advection and dispersion with a Péclet number of 40 as in Appelo (1994). See Table 7 for the details on the experimental setups.

Table 7 Boundary conditions for cation exchange simulations.

^a J_θ: addition rate of solid species θ at the upper boundary of the calculation domain, N: number of grid cells in the calculation domain, z_tot: total depth of the calculation domain, w: uplift/erosion rate, r_H: hydraulic radius of particles for solid phases, q: annual runoff, σ₀: water saturation ratio at the surface, z_sat: water table depth, CEC_inrt: cation exchange capacity assumed for bulk species, $c_{\mathit{\varsigma }}^{\mathrm{0}}$: concentration of ς at the surface (ς= Na, K, Ca, NO₃, Cl), $K_{\mathit{\varsigma }\backslash \mathrm{H},\mathrm{inrt}}$: thermodynamic constant for ς–H exchange (ς= Na, K, Ca) for bulk solid species.
^b Run as a restart from the spin-up whose boundary condition is given in the second column with $c_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{0}}=\mathrm{1.2}$ mmol L⁻¹ and N=100. To satisfy the Péclet number of 40, diffusion coefficients for all aqueous species are set at $\mathrm{6.5975}\times {\mathrm{10}}^{-\mathrm{2}}$ (m² yr⁻¹) because the tortuosity factor is calculated to be 0.3789, with assumed porosity and water saturation (Kanzaki et al., 2022).
^c Only the IDs of solid species are denoted (inrt: bulk species).
^d Equivalent to 1.1 meq L⁻¹, with assumed porosity and water saturation.
^e Calculated to be consistent with the thermodynamic dataset Tipping_Hurley.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013).

3 Model validation

Before we examine the validity of our framework for soil pH calculation, the model's capacity to simulate cation exchange is compared with that of PHREEQC v3.0 (https://www.usgs.gov/software/phreeqc-version-3, last access: 7 June 2023; Fig. 2). First, a series of experiments (Fig. 2a) is conducted with the two models with common cation exchange thermodynamics (Table 7) in order to compare the relationship between solution pH and base saturation at equilibrium when a solution with fixed concentrations of cations (1 mM Na and 0.2 mM K) and different concentrations of nitrate (1 to 15 mM) is brought to equilibrium with a 1.1 meq L⁻¹ cation exchanger. Second, we perform a cation exchange simulation (Fig. 2b) in which a cation-exchanging soil column initially homogeneously equilibrated with porewater consisting of 1 mM Na, 0.2 mM K, and 1.2 mM NO₃ is flushed by a 0.6 mM CaCl₂ solution through advection and dispersion with a Péclet number of 40, as in Appelo (1994), again using the same cation exchange properties (Table 7). We find negligible differences in the equilibria and dynamics of solutes and exchangeable cations between the two models (Fig. 2), indicating that the capacity of the current version of SCEPTER to simulate cation exchange is comparable to that of PHREEQC v3.0 (see also the Supplement for the effect of cation exchange in some of the previous example simulations run by Kanzaki et al., 2022).

Table 8 Compositional data measured for mesocosm soil sample.

Table 9 Boundary conditions for mesocosm simulations.

^a J_θ: addition rate of solid species θ at the upper boundary of the calculation domain, N: number of grid cells in the calculation domain, z_tot: total depth of the calculation domain, w: uplift/erosion rate, z_ml: mixed layer depth, r_H: hydraulic radius of particles for solid phases, q: annual runoff, σ₀: water saturation ratio at the surface, z_sat: water table depth, CEC_θ: cation exchange capacity for solid species θ, $c_{\mathrm{Cl}}^{\mathrm{0}}$: concentration of Cl at the surface, K_ς\H: intrinsic thermodynamic constant for ς–H exchange (ς= Na, K, Mg, and Ca), α_θ: coefficient to describe surface charge effect on cation exchange thermodynamics for solid species θ (Sect. 2.1; Appelo, 1994).
^b Only IDs of solid species are denoted; inrt: bulk species, amnt: NH₄NO₃, g1: SOM Class 1 (most labile class), g2: SOM Class 2 (second most labile class), na2o: Na₂O, k2o: K₂O, mgo: MgO, cao: CaO, kcl: KCl, cacl2: CaCl₂, teas: triethanolamine, ims: imidazole, mesmh: MES (2-(N-morpholino)ethanesulfonic acid) monohydrate, gac: acetic acid, naoh: NaOH.
^c Added only when simulating soil pH in CaCl₂ solution.
^d Added only when simulating soil buffer pH by Sikora (2006).
^e Some of aqueous species in Sikora buffer are abbreviated; TEA: triethanolamine, IM: imidazole, MES: 2-(N-morpholino)ethanesulfonic acid, AcO: Acetate anion.
^f Parameter values optimized for the reproduction of observations (Sect. 3).
^g CEC_θ=0 for solid species not listed here.
^h See Sect. 3 for base cation concentrations at the upper boundary.
ⁱ Those values are applied only to bulk and SOM Class 2 species.

Figure 3 Comparison of soil composition between our model simulation and observations from the mesocosm experiments. (a) Porewater chemistry at 15 cm soil depth. (b) Concentrations of exchangeable cations over top 15 cm. A uniform 10 % error is assumed for observational measurements (see Table 8).

Figure 4 (a) Comparison of porewater and soil buffer pH between mesocosm observations and model simulation. (b) Data-model comparison of Sikora buffer pH (2006) as a function of neutralized acidity.

In order to validate our approach toward calculating bulk agronomic soil pH in the reaction–transport model, we compare a series of soil pH simulations fed by a field simulation with observed boundary conditions to results from a mesocosm experiment. The mesocosm has been monitored since July 2022 at a greenhouse controlled under average growing season conditions. The field simulation is constrained using detailed measurements conducted in August 2022 (Table 8) as boundary conditions (Table 9). The field simulation is simplified as much as possible as the focus of this paper is the simulation of soil pH (see Kanzaki et al., 2022, for some additional examples of field simulations fitted to observations and the Supplement). A detailed description of the mesocosm setup can be found in Chiaravalloti et al. (2023). Its tracked solid species are limited to soil organic matter (SOM) and a bulk solid-phase species (i.e., a hypothetical species representing the solid phases other than soil organic matter dumped together as a whole) treated as two cation exchangers; tracked aqueous species include base cations (Na, K, Ca, and Mg), NO₃, and Cl, as well as CO₂ gas. The tracked solid species (i.e., SOM and bulk species) are assumed to have the same values for thermodynamic parameters for cation exchange except that they have different cation exchange capacity (CEC) values (120 and 3.176 cmol kg⁻¹, respectively), with their average constrained by the observed bulk soil CEC (8.9 cmol kg⁻¹). Measured porewater composition at 15 cm depth is used as the upper boundary condition for aqueous base cations in the field simulation so that simulated porewater composition at 15 cm depth is consistent with observations (Fig. 3a). Aqueous NO₃ is added as the NH₄NO₃ fertilizer at the upper boundary in the field simulation at the same rate of total N supply as the urea–NH₄–NO₃ fertilizer is applied to the mesocosm (24.210 g N m⁻² yr⁻¹). Upper aqueous Cl concentration takes a fitted value (Table 9) so that the simulated porewater pH at 15 cm depth is consistent with observations (Fig. 4a). SOM input is fixed at the value (Table 9) with which simulated average organic matter concentration over the top 15 cm is consistent with observations (4.9 wt %). See Table 9 for more details on the boundary conditions for the field simulation.

Soil samples from the mesocosm experiments were homogenized from the top 15 cm of soil (dried at 60 °C overnight and sieved at 2 mm), and measured soil pH values and electrical conductivity values were obtained from a series of solutions in deionized water at soil-to-solution ratios of 1:5, 1:2, 1:1, and 1:0.5 (g cm⁻³) and in 0.0025, 0.005, and 0.01 M CaCl₂ solutions at a 1:1 soil-to-solution ratio (g cm⁻³). The pH for each soil/solution slurry was measured with a Thermo Scientific Orion ROSS Ultra pH/ATC Triode paired with a Thermo Scientific Orion STARA2215 Orion Star A221 Portable pH Meter (ThermoFisher Scientific, Massachusetts). Electrical conductivity was measured by placing a few drops of the liquid from the soil/solution slurry on a HOBO U24 Conductivity Logger (U24-002-C) (Onset Computer Corporation, Massachusetts). We also measured buffer pH from a soil split using the method and recipe developed by Sikora (2006).

Soil pH simulations are conducted based on averaged data over the top 15 cm of bulk soil from the field simulation described above, supplemented with the mesocosm observations according to the procedure described in Sect. 2.2. A series of soil pH values is calculated in deionized water at soil-to-solution ratios of 1:5, 1:2, 1:1, and 1:0.5 (g cm⁻³) and in 0.0025, 0.005, and 0.01 M CaCl₂ solutions at a 1:1 soil-to-solution ratio (g cm⁻³) following Miller and Kissel (2010). We also calculate soil buffer pH, where the bulk soil over the upper 15 cm, deionized water, and Sikora buffer solution are mixed at a $\mathrm{1}:\mathrm{1}:\mathrm{1}$ ratio (g : cm³ : cm³), following the recipe by Sikora (2006). The observation shows significant amounts of extractable NO₃ and Cl (Table 8), which probably exist as some forms of salts, given measured electrical conductivities, and are not explicitly simulated in the field run. Therefore, those extractable anions are added to the laboratory runs so that all major extractable/exchangeable elements measured in the mesocosm samples are consistent between the laboratory simulation and observations.

Figure 5 (a) Soil pH in deionized water at different soil-to-solution ratios and in a CaCl₂ solution at different concentrations plotted against electrical conductivity for both simulations and mesocosm observations. (b) Difference in soil pH at a 1:1 soil-to-solution g cm⁻³ ratio between in deionized water and a 0.01 M CaCl₂ solution (ΔpH_1:1) plotted against electrical conductivity for both simulated and observed mesocosm, along with the ΔpH_1:1 relationship with electrical conductivity derived for US soils by Miller and Kissel (2010). In panels (a) and (b), measured pH is assumed to have a uniform error of 0.02.

The simulated field run shows an abundance of exchangeable cations over the top 15 cm that matches well with observations (Fig. 3b) with the optimized thermodynamic parameters for cation exchange (Table 9). Slight offsets might be attributable to chemical gradients developed especially for relatively strongly bound Ca and Mg, caused by CEC variation with depth. The simulated soil buffer pH is also consistent with the observed buffer pH for the topsoil of the mesocosm (Fig. 4a). Although soil buffer pH was measured using the Sikora buffer, we can confirm that the model can effectively reproduce the relationship between Sikora buffer pH and neutralized acid measured by the Sikora method (Fig. 4b). Therefore, the in silico measurement of soil buffer pH should be directly comparable with the observational data. Simulated soil pH varies as a function of dilution by deionized water and/or the concentration of CaCl₂ in solution, a trend especially obvious when soil pH is plotted against electrical conductivity as shown in Fig. 5a (in silico electrical conductivity is calculated from ionic strength, assuming a conversion factor of 0.016 dS m⁻¹ mol⁻¹ L from Ponnamperuma et al., 1966; cf. Alva et al., 1991). This trend is also consistent with observations (Table 10 and Fig. 5a). The difference in soil pH in deionized water from that in a 0.01 M CaCl₂ solution at the same soil-to-solution ratio of 1:1 (g cm⁻³), defined as ΔpH_1:1 (Miller and Kissel, 2010), is also consistent with the mesocosm observations and the trend observed for US soils by Miller and Kissel (2010) (Fig. 5b). Overall, with optimized thermodynamics of cation exchange, the model can very closely reproduce observed porewater and soil (buffer) pH results for both our mesocosm experiments and previously published data (Miller and Kissel, 2010).

Table 10 Porewater and soil (buffer) pH of the mesocosm.

Table 11 Boundary conditions for EW simulations.

^a J_θ: addition rate of solid species θ at the upper boundary of the calculation domain, N: number of grid cells in the calculation domain, z_tot: total depth of the calculation domain, w: uplift/erosion rate, z_ml: mixed layer depth, r_H: hydraulic radius of particles for solid phases, q: annual runoff, σ₀: water saturation ratio at the surface, z_sat: water table depth, CEC: cation exchange capacity assumed for bulk species and SOM, $c_{\mathrm{Ca}}^{\mathrm{0}}$: concentration of Ca at the surface, K_Ca\H: thermodynamic constant for Ca–H exchange.
^b Only IDs of solid species are denoted; inrt: bulk species, amnt: NH₄NO₃, g1: SOM Class 1 (most labile class), g2: SOM Class 2 (second most labile class), gbas: glassy basalt, na2o: Na₂O, k2o: K₂O, mgo: MgO, cao: CaO. The chemical composition of glassy basalt is given by the stoichiometry of ${\mathit{\gamma }}_{\mathrm{gbas},\mathit{\varsigma }}/{\mathit{\gamma }}_{\mathrm{gbas},\mathrm{Si}}=\mathrm{0.0809}$, 0.0084, 0.2439, 0.2722, 0.1251, 0.4683, and 1 for ς= Na, K, Ca, Mg, Fe, Al, and Si, respectively.
^c Only enabled when basalt is applied in a field run or soil pH is simulated for basalt-applied soils.
^d See Sect. 4 for the calculation of those parameter values.
^e See Fig. 6.
^f Bio-mixing is defined using a modified transition matrix (K_θ,ij), which is a discretized form of the continuous exchange function E_θ in Eq. (11) and can be formulated based on transport probability between soil layers i and j (P_θ,ij). Inversion mixing in this paper is implemented as $K_{\mathit{\theta },ij}=\mathit{\delta }{z}_i{P}_{\mathrm{inv}}/\mathit{\delta }{z}_j$ if $i=j-\mathrm{1}$ or $i=j+\mathrm{1}$ or $i=n_{\mathrm{ml}}+\mathrm{1}-j$, else 0, where P_θ,ij is assumed to have a phase- and location-independent value of P_inv=0.1 yr⁻¹; δz_i is the thickness (m) of soil layer i; and n_ml is the total number of mixed layers. See Kanzaki et al. (2022) for the formulation for Fickian mixing.
^g Other thermodynamic constants for cation exchange are modified from their default values in Table 2 consistently with the change in K_Ca\H; for example, log $K_{\mathit{\varsigma }\backslash \mathrm{H}}=-\mathrm{4.389}$, −3.289, and −7.764 for ς= Na, K, and Mg, respectively.

4 Example EW application

To illustrate the potential importance of distinguishing between pH_s and pH_pw and modeling both accurately, we present example simulations in which the alkalinity addition to soils through EW for 1 year is limited by an assumed target pH and compare cases in which the target value is assumed to be pH_s with an equivalent ensemble in which it is assumed to be pH_pw. Here, we consider another simple soil system which enables us to focus on any potential difference between pH_s and pH_pw: tracked solid species include the bulk and SOM species, aqueous species are Ca²⁺ and NO₃, and CO₂ is the only tracked gaseous species. Boundary conditions are those of an arbitrarily chosen field site from the midwestern USA (Table 11), and the cation exchange thermodynamics, soil respiration, and base saturation are correspondingly constrained from the observation at the site. More specifically, the system is tuned by varying Ca²⁺ concentration at the upper boundary, thermodynamic coefficients for Ca–H exchange, and organic matter input to soil (three unknown variables) until the system satisfies the observed soil pH, exchangeable acidity, and SOM wt % (6.058, 20.980 % CEC, and 2.052 wt %, respectively; three observed known variables) at steady state. Mechanistically, the non-zero value of Ca²⁺ concentration at the upper boundary can be taken to reflect the net result of historical liming at the site. We then add a glassy basalt solid species iteratively to meet a range of target pH values (6.2, 6.5, and 6.8) after 1 year of basalt application and use the model to estimate the rate of basalt application required to achieve a given target pH value. We run two ensembles, one with pH_s as the operative target pH and one with pH_pw as the target, allowing us to compare the estimated basalt feedstock application required to reach the identical target pH when using pH_s or pH_pw as an index. For comparison with the target values, soil pH is calculated in a mixture of top 15 cm bulk soil and deionized water at a 1:1 g cm⁻³ ratio, while the average over the top 15 cm is considered for porewater pH, calculated as $-\log \left({\int }_{\mathrm{0}}^{\mathrm{0.15}}\right[{\mathrm{H}}^{+}]\mathrm{d}z/\mathrm{0.15})$. The observed data used for the initial spin-up/tune-up is from Fick and Hijmans (2017) for temperature, Wang et al. (2021) for soil moisture, Reitz et al. (2017) for runoff, Poggio et al. (2021) for soil pH and organic matter, Walkinshaw et al. (2022) for cation exchange capacity, Pan et al. (2021) for nitrification rate, and ISRIC (2022) for base saturation. Basalt application simulations are all conducted as re-starts from the end of the same spin-up/tune-up described above, where glassy basalt is applied and mixed with bulk soil via tilling during the initial 0.005 years (∼2 d). For glassy basalt, we use the kinetic law formulated by Brantley et al. (2008) and the thermodynamic calculation method used by Aradóttir et al. (2012) and Pollyea and Rimstidt (2017) and assume a log-normal distribution centered at 10 µm with a 0.2 log unit standard deviation for the initial particle size distribution, and the chemical composition in the footnote of Table 11. See Table 11 for additional details on model boundary conditions.

Figure 6 Basalt requirements for different target pH values after the first year following feedstock application using either bulk soil or porewater pH averaged over 0–15 cm as a pH reference value.

Figure 7 Evolution of (a–c) soil and porewater pH and (d–f) exchangeable acidity during the first year following basalt feedstock application at target pH values (vertical dotted lines) of (a, d) 6.2, (b, e) 6.5, and (c, f) 6.8 using soil pH averaged over 0–15 cm as a pH reference. Note that the curves of soil pH and porewater pH depicted here show values at individual depth points, while target pH values are integrated across the top 0–15 cm.

Figure 8 Evolution of (a–c) soil and porewater pH and (d–f) exchangeable acidity during the first year following basalt feedstock application at target pH values (vertical dotted lines) of (a, d) 6.2, (b, e) 6.5, and (c, f) 6.8 using porewater pH averaged over 0–15 cm as a pH reference. Note that the curves of soil pH and porewater pH depicted here show values at individual depth points, while target pH values are integrated across the top 0–15 cm.

Depending on the pH reference (i.e., either soil pH, pH_s, or porewater pH, pH_pw), the required amount of basalt is significantly different at any of the target pH values examined here (Fig. 6), though all are within comparable ranges to previous EW deployments (e.g., Swoboda et al., 2022). Comparison of soil and porewater pH (Figs. 7 and 8) shows that variation in soil pH is more limited compared to that of porewater pH because soil pH largely reflects exchangeable acidity, which can more effectively buffer the input of alkalinity compared to acidity of porewater although the total exchangeable acidity is dependent on the cation exchange capacity and initial base saturation of soil. Porewater pH is lower than soil pH at a relatively low alkalinity input (e.g., at an earlier time after basalt deployment and/or at deep depths; Figs. 7 and 8), given that in situ porewater pH reflects higher soil pCO₂ while soil pH has lower re-equilibrated pCO₂ from conserved DIC because of dilution by deionized water. With a higher alkalinity input (e.g., at a later time after basalt deployment and/or at shallower depths; Figs. 7 and 8), porewater pH is higher than soil pH because soil pH has a maximum value set by the cation exchange capacity at 100 % base saturation. In general, using pH_pw as the index target requires higher alkalinity input via basalt dissolution for a given target pH value because pH_pw is lower than pH_s in the background, and it requires pH_pw to reach higher values at shallower depths to compensate for the lower background pH_pw at deeper depths. Though only meant to be illustrative, the example simulations shown here demonstrate the importance of distinguishing between soil and porewater pH in numerical frameworks for representing soil pH regulation. Our results indicate that care must be taken in the reporting and validation of simulated pH values in reaction–transport models, particularly when comparing to analytical data for a bulk (multiphase) parameter such as soil pH.

5 Conclusions

We update the SCEPTER model (v1.0) to simulate the mechanics of cation exchange, and an associated, newly developed framework that enables the calculation of soil pH in silico. By comparing it to observational measurements from mesocosm experiments, we demonstrate that the soil pH simulation in SCEPTER can accurately reproduce systematic variations in observed porewater pH, soil pH, and soil buffer pH so long as a field simulation can be validated by accessory soil chemistry. We also present example simulations which focus on the application of the model to the estimation of required basalt for agricultural soils to reach different target pH values through EW. We observe significant differences in response to an alkalinity input via basalt dissolution between porewater and soil pH, with important implications for diagnosing agricultural soils with respect to an optimal basalt deployment rate/style through EW and managing crop yields. Future model developments include an extension of cation exchange to a more generalized suite of sorption reactions, e.g., the implementation of anion (e.g., PO₄) adsorption onto oxides (e.g., van der Zee and van Riemsdijk, 1988; McGechan and Lewis, 2002) and nutrient uptake by plants to comprehensively predict nutrient cycling and productivity in cropland soils in parallel with anthropogenic alkalinity modification and CO₂ removal through EW. Also, tracking additional potential aqueous pH-affecting agents (e.g., dissolved organic matter; e.g., Nambu and Yonebayashi, 1999; Grybos et al., 2009) will widen the soil conditions to which the model can be applied to estimate shifts in porewater and soil pH as a result of EW. Further and wider use of the current/future code coupled with mesocosm/field observations of porewater pH and soil pH is expected to enhance our mechanistic understanding of the agronomic and climatic impacts of EW in croplands.

Appendix A: List of symbols

Symbol	Definition	Units
$B_{\mathit{\varsigma }}^{\mathrm{ads}}$	factor to convert $c_{\mathit{\varsigma }}^{\mathrm{1}}$ to total concentration of element ς adsorbed onto solid phases	m−3 L
cΘ	solution concentration of electrolyte Θ	mol L−1
cς	total concentration of dissolved element ς	mol L−1
$c_{\mathit{\varsigma }}^{\mathrm{1}}$	concentration of free dissolved species for ς or H4SiO4 if ς= Si	mol L−1
Cς	concentration of exchangeable/extractable cation/anion ς	mol m−3
CECθ	cation exchange capacity of exchanger θ	eq g−1
Dς	diffusion coefficient of dissolved element ς	m2 yr−1
$E_{\mathit{\theta }}(z,z^{\prime })$	rate of particle transfer between locations at z and z′ by bio-mixing	m−1 yr−1
f(i)	charge-equivalent fraction of surface species i	–
〈i〉	concentration of surface species i	mol g−1
[j]	concentration of aqueous species j	mol L−1
Kς,i	thermodynamic constant for production of ith aqueous species of dissolved element ς	variablea
$K_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}$, $K_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}^{\prime }$	intrinsic and apparent equilibrium constants for cation exchange, respectively	${\mathrm{mol}}^{\mathrm{1}-\mathrm{1}/Z_{\mathit{\varsigma }}}{\mathrm{L}}^{\mathrm{1}/Z_{\mathit{\varsigma }}-\mathrm{1}}$
ℓ	unit conversion factor	L m−3
mθ	concentration of solid species θ	mol m−3
Mθ	molar weight of solid species θ	g mol−1
naq, ngas, nsld	total number of simulated aqueous, gaseous, solid species, respectively	–
nxrxn	total number of extra reactions	–
pε	partial pressure of gas species ε	atm
Rθ	net dissolution rate of solid species θ	mol m−3 yr−1
Rκ	rate of κth extra reaction	mol m−3 yr−1
t	time	years
X(θ)−	exchangeable surface sites of solid-phase exchanger θ	–
v	porewater advection rate	m yr−1
w	advection rate of solid phases	m yr−1
z	depth of weathering profile	m
zlab	depth of laboratory beaker/flask filled with the mixture of soil sample and solution	m
zml	mixed layer depth	m
Zς	valence number of cation ς	–
αθ	factor to represent relation of ηH,θ to f(X(θ)−)	–
${\mathit{\beta }}_{\mathit{\varsigma }}^{\mathrm{aq}}$	factor to convert $c_{\mathit{\varsigma }}^{\mathrm{1}}$ to total concentration of dissolved element ς	–
γθ,Θ	mole of electrolyte Θ in 1 mole of solid species θ	–
γθ,ς	mole amount of ς released upon dissolution of 1 mole of solid species θ	–
γκ,ς	stoichiometry of ς production in κth extra reaction	–
${\mathit{\gamma }}_{\mathit{\varsigma },i,\mathrm{p}}$, ${\mathit{\gamma }}_{\mathit{\varsigma },i,{\mathit{\varsigma }}^{\prime }}$, ${\mathit{\gamma }}_{\mathit{\varsigma },i,\mathit{\epsilon }}$	stoichiometry of H+, dissolved element ς′ and gas species ε, respectively, in production ith aqueous species of ς	–
ηH,θ	factor to reflect the effect of surface potential on $K_{\mathit{\varsigma }\backslash \mathrm{H},\mathit{\theta }}^{\prime }$	–
ρ	bulk soil particle density	g cm−3
σ	water saturation ratio	–
$\mathit{\varsigma }\text{-}\mathrm{X}(\mathit{\theta }{)}_{Z_{\mathit{\varsigma }}}$	cation ς adsorbed onto exchangeable sites of θ	–
τaq	tortuosity factor for solute diffusion in porewater	–
ϕ	porosity	–
ϕlab	volume ratio of fluid against solid plus fluid phases	–
ψ	soil-to-solution ratio used in the laboratory	g cm−3

^a Units can vary depending on reactions considered.