Articles | Volume 17, issue 10
https://doi.org/10.5194/gmd-17-4515-2024
https://doi.org/10.5194/gmd-17-4515-2024
Model description paper
 | 
31 May 2024
Model description paper |  | 31 May 2024

In silico calculation of soil pH by SCEPTER v1.0

Yoshiki Kanzaki, Isabella Chiaravalloti, Shuang Zhang, Noah J. Planavsky, and Christopher T. Reinhard
Abstract

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 (pHs) and porewater pH (pHpw). 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 pHs. 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 pHs and pHpw 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 CO2 (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 CO2 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 CO2 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 pHs – which measures H+ residing in bulk soil samples that is measured in the laboratory as the pH of liquid extractants (deionized water or KCl/CaCl2 solution) of bulk soil samples taken from the field. The other is porewater pH – defined here as pHpw – 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 pHs (e.g., Thomas, 1996), while models that simulate biogeochemical reactions and transport within soils, including the dissolution of basalt during EW, typically calculate pHpw (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 (pHs) along with porewater pH (pHpw). 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 pHs or pHpw 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 1Default cation exchange capacity of solid speciesa.

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

Download Print Version | Download XLSX

2.1 Cation exchange in SCEPTER

The current version of SCEPTER allows for cation exchange involving H+, Na+, K+, Ca2+, Mg2+, and Al3+ 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:

(1) Z ς X ( θ ) - + ς Z ς + ς - X ( θ ) Z ς ,

where Zς is the valence number of cation ς; X(θ) denotes exchangeable surface sites of solid-phase exchanger θ; and ς-X(θ)Zς 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

(2) ( 1 / Z ς ) ς Z ς + + H - X ( θ ) ( 1 / Z ς ) ς - X ( θ ) Z ς + H + .

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

(3) K ς \ H , θ = f 1 / Z ς ( ς - X ( θ ) Z ς ) [ H + ] f ( H - X ( θ ) ) [ ς Z ς + ] 1 / Z ς ,

where f(i) denotes the charge-equivalent fraction of surface species i, and [j] represents the concentration of aqueous species j (mol L−1). The apparent equilibrium constant Kς\H,θ can vary as a result of surface fraction X(θ), and we adopt the formulation by Appelo (1994):

(4) K ς \ H , θ = η H , θ K ς \ H , θ .

Here, Kς\H,θ 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):

(5) log η H , θ = - α θ { 1 - f ( H - X ( θ ) ) } ,

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:

(6) CEC θ = ς Z ς ς - X ( θ ) Z ς ,

where CECθ is the cation exchange capacity of exchanger θ (eq g−1) and i is the concentration of surface species i (mol g−1). By definition, we have

(7) f ( ς - X ( θ ) Z ς ) Z ς ς - X ( θ ) Z ς CEC θ .

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

(8) 1 = ς f ( ς - X ( θ ) Z ς ) .

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

(9) 1 = f ( H - X ( θ ) ) + ς H η H , θ K ς \ H , θ f ( H - X ( θ ) ) [ H + ] Z ς [ ς Z ς + ] .

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 2Default thermodynamic data of cation exchangea.

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ς\Na=1.1, 0.507, and 0.665 for ς=K, Mg, and Ca, respectively, from Appelo (1994). (3) Calculated from log KAl\Na=0.41 from phreeqc.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013).

Download Print Version | Download XLSX

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 H4SiO4 concentration is tracked), denoted as cς1. These cς and cς1 are related to one another by the following equation (Kanzaki et al., 2022):

(10) c ς = c ς 1 + c ς 1 i = 2 n ς K ς , i [ H + ] γ ς , i , p ς ς n aq ( c ς 1 ) γ ς , i , ς ε n gas ( p ε ) γ ς , i , ε ,

where the second term on the right-hand side is the sum of the concentrations of dissolved element ς other than cς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 ς; γς,i,p, γς,i,ς, and γς,i,ε 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 naq and ngas 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ς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

(11) ϕ σ β ς aq c ς 1 t + B ς ads c ς 1 t = - ϕ σ v β ς aq c ς 1 z + z ϕ σ τ aq D ς β ς aq c ς 1 z + θ n sld γ θ , ς R θ + κ n xrxn γ κ , ς R κ + w B ς ads c ς 1 z - B ς ads c ς 1 0 z ml E θ ( z , z ) d z + 0 z ml B ς ads c ς 1 ( z ) E θ ( z , z ) d z .

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 βςaq and Bςads (m−3 L) defined as the factors to convert cς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 βςaq and Bςads can be calculated from Eqs. (10) and (3)–(9), respectively:

(12)βςaqcςcς1=1+i=2nςKς,i[H+]γς,i,pςςnaq(cς1)γς,i,ςεngas(pε)γς,i,ε,(13)Bςads1cς1θmθMθς-X(θ)Zς=θmθMθCECθZςηH,θKς\H,θf(H-X(θ))[H+]Zς(ς{Na,K,Ca,Mg,Al})0(else),

where mθ and Mθ are the concentration (mol m−3) and molar weight (g mol−1) 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 (pHs) imitates the procedure to measure soil pH in the laboratory: sampling bulk soils, mixing them with an extractant solution (e.g., deionized water or KCl/CaCl2 solution) at a given soil-to-solution ratio (e.g., 1:1 or 1:5 g mL−1), 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 zlab(1-ϕlab)CςMθ/γθ,ς (g m−2), where zlab (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 ϕlab=ψ-1(ρ-1+ψ-1)-1 with the soil-to-solution ratio used in the laboratory (ψ; g cm−3) and bulk soil particle density (ρ; g cm−3) observed in the in silico field; Cς is the concentration of exchangeable/extractable cation/anion ς (mol m−3); Mθ is the molar weight of the added solid species, θ (g mol−1); 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 zlabϕlabcΘMθ/γθ,Θ (g m−2), where cΘ and γθ are the solution concentration (mol L−1) of electrolyte Θ and the mole of electrolyte Θ in 1 mol of salt θ, respectively (e.g., cΘ=0.01 mol L−1 and γθ,Θ=1 if θ=CaCl2; else γθ,Θ=0, when Θ=CaCl2). 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 pHs and pHpw is confirmed to be a general phenomenon, as discussed below. In the next section, we discuss the validity of our approach toward simulating pHs 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 3Solid species to be dissolved in laboratory simulationsa.

a Thermodynamic constants for solid species θ (Kθ) are calculated as Kθ=Kθrefexp(-ΔHθ(T-1-298-1)R-1), where T is the temperature in kelvins and is the gas constant in units of kJ mol−1 K−1 (R=8.314×10-3 kJ mol−1 K−1). 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θref and ΔHθ from llnl.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013). (3) Kθref and ΔHθ are assumed to be the same as those for amorphous Si. (4) Kθ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−3 for NaOH, NH4NO3, SOM Class 1, triethanolamine, imidazole, MES (2-(N-morpholino)ethanesulfonic acid) monohydrate, and acetic acid, respectively.

Download Print Version | Download XLSX

Table 4Sikora 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.

Download Print Version | Download XLSX

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

 (1) The diffusion coefficient (m2 yr−1) is calculated as D=0.4415 (μw-1.1a0.6)−1, where μw is the water viscosity (mPa s) and a is the molar volume of solute (cm3 mol−1) (Othmer and Thakar, 1953; La-Scalea et al., 2005). The water viscosity μw is calculated as μw=0.024152exp(4.7428(T-139.86)-1R-1), where R=8.314×10-3 kJ mol−1 K−1 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 (m2 yr−1) is calculated as D=a×exp(-b(T-1-288-1)R-1) where a is the pre-exponential factor (m2 yr−1) and b is the apparent activation energy (kJ mol−1). (5) a and b from Schulz and Zabel (2006). (6) a and b from Li and Gregory (1974).

Download Print Version | Download XLSX

Table 6Thermodynamic data for aqueous species in the Sikora buffera.

a The thermodynamic constant (Kaq) is calculated as Kaq=Kaqrefexp(-ΔHaq(T-1-298-1)R-1), where R=8.314×10-3 kJ mol−1 K−1 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) Kaqref from Sikora (2006) and ΔHaq from Goldberg et al. (2002). (2) From llnl.dat, available in PHREEQC v3.0 (Parkhurst and Appelo, 2013).

Download Print Version | Download XLSX

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f01

Figure 1Schematic 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 (pHs) of the in silico field soil in Steps 1 and 2.

Download

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f02

Figure 2Comparison 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 NO3 in equilibrium with a 1.1 meq L−1 cation exchanger. (b) A 0.6 mM CaCl2 solution replacing a solution of 1 mM Na, 0.2 mM K and 1.2 mM NO3, initially equilibrated with a 1.1 meq L−1 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.

Download

Table 7Boundary 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, ztot: total depth of the calculation domain, w: uplift/erosion rate, rH: hydraulic radius of particles for solid phases, q: annual runoff, σ0: water saturation ratio at the surface, zsat: water table depth, CECinrt: cation exchange capacity assumed for bulk species, cς0: concentration of ς at the surface (ς= Na, K, Ca, NO3, Cl), Kς\H,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 cNO30=1.2 mmol L−1 and N=100. To satisfy the Péclet number of 40, diffusion coefficients for all aqueous species are set at 6.5975×10-2 (m2 yr−1) 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−1, 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).

Download Print Version | Download XLSX

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−1 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 NO3 is flushed by a 0.6 mM CaCl2 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 8Compositional data measured for mesocosm soil sample.

Download Print Version | Download XLSX

Table 9Boundary 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, ztot: total depth of the calculation domain, w: uplift/erosion rate, zml: mixed layer depth, rH: hydraulic radius of particles for solid phases, q: annual runoff, σ0: water saturation ratio at the surface, zsat: water table depth, CECθ: cation exchange capacity for solid species θ, cCl0: 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: NH4NO3, g1: SOM Class 1 (most labile class), g2: SOM Class 2 (second most labile class), na2o: Na2O, k2o: K2O, mgo: MgO, cao: CaO, kcl: KCl, cacl2: CaCl2, 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 CaCl2 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.
i Those values are applied only to bulk and SOM Class 2 species.

Download Print Version | Download XLSX

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f03

Figure 3Comparison 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).

Download

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f04

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.

Download

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), NO3, and Cl, as well as CO2 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−1, respectively), with their average constrained by the observed bulk soil CEC (8.9 cmol kg−1). 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 NO3 is added as the NH4NO3 fertilizer at the upper boundary in the field simulation at the same rate of total N supply as the urea–NH4–NO3 fertilizer is applied to the mesocosm (24.210 g N m−2 yr−1). 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−3) and in 0.0025, 0.005, and 0.01 M CaCl2 solutions at a 1:1 soil-to-solution ratio (g cm−3). 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−3) and in 0.0025, 0.005, and 0.01 M CaCl2 solutions at a 1:1 soil-to-solution ratio (g cm−3) 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 1:1:1 ratio (g : cm3: cm3), following the recipe by Sikora (2006). The observation shows significant amounts of extractable NO3 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.

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f05

Figure 5(a) Soil pH in deionized water at different soil-to-solution ratios and in a CaCl2 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−3 ratio between in deionized water and a 0.01 M CaCl2 solution (ΔpH1:1) plotted against electrical conductivity for both simulated and observed mesocosm, along with the ΔpH1: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.

Download

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 CaCl2 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−1 mol−1 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 CaCl2 solution at the same soil-to-solution ratio of 1:1 (g cm−3), defined as ΔpH1: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 10Porewater and soil (buffer) pH of the mesocosm.

Download Print Version | Download XLSX

Table 11Boundary 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, ztot: total depth of the calculation domain, w: uplift/erosion rate, zml: mixed layer depth, rH: hydraulic radius of particles for solid phases, q: annual runoff, σ0: water saturation ratio at the surface, zsat: water table depth, CEC: cation exchange capacity assumed for bulk species and SOM, cCa0: concentration of Ca at the surface, KCa\H: thermodynamic constant for Ca–H exchange.
b Only IDs of solid species are denoted; inrt: bulk species, amnt: NH4NO3, g1: SOM Class 1 (most labile class), g2: SOM Class 2 (second most labile class), gbas: glassy basalt, na2o: Na2O, k2o: K2O, mgo: MgO, cao: CaO. The chemical composition of glassy basalt is given by the stoichiometry of γgbas,ς/γgbas,Si=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θ,ij=δziPinv/δzj if i=j-1 or i=j+1 or i=nml+1-j, else 0, where Pθ,ij is assumed to have a phase- and location-independent value of Pinv=0.1 yr−1; δzi is the thickness (m) of soil layer i; and nml 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 KCa\H; for example, log Kς\H=-4.389, −3.289, and −7.764 for ς= Na, K, and Mg, respectively.

Download Print Version | Download XLSX

4 Example EW application

To illustrate the potential importance of distinguishing between pHs and pHpw 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 pHs with an equivalent ensemble in which it is assumed to be pHpw. Here, we consider another simple soil system which enables us to focus on any potential difference between pHs and pHpw: tracked solid species include the bulk and SOM species, aqueous species are Ca2+ and NO3, and CO2 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 Ca2+ 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 Ca2+ 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 pHs as the operative target pH and one with pHpw as the target, allowing us to compare the estimated basalt feedstock application required to reach the identical target pH when using pHs or pHpw 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−3 ratio, while the average over the top 15 cm is considered for porewater pH, calculated as -log(00.15[H+]dz/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.

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f06

Figure 6Basalt 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.

Download

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f07

Figure 7Evolution 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.

Download

https://gmd.copernicus.org/articles/17/4515/2024/gmd-17-4515-2024-f08

Figure 8Evolution 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.

Download

Depending on the pH reference (i.e., either soil pH, pHs, or porewater pH, pHpw), 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 pCO2 while soil pH has lower re-equilibrated pCO2 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 pHpw as the index target requires higher alkalinity input via basalt dissolution for a given target pH value because pHpw is lower than pHs in the background, and it requires pHpw to reach higher values at shallower depths to compensate for the lower background pHpw 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., PO4) 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 CO2 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ςads factor to convert cς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ς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θ(z,z) 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ς\H,θ, Kς\H,θ intrinsic and apparent equilibrium constants for cation exchange, respectively mol1-1/ZςL1/Zς-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(θ))
βςaq factor to convert cς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
γς,i,p, γς,i,ς, γς,i,ε 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ς\H,θ
ρ bulk soil particle density g cm−3
σ water saturation ratio
ς-X(θ)Zς 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.

Code and data availability

The source codes of the model are available at GitHub (https://github.com/cdr-laboratory/SCEPTER, last access: 11 March 2024) under the GNU General Public License v3.0. The specific version of the model used in this paper is tagged as v1.0.1 and has been assigned a DOI (https://doi.org/10.5281/zenodo.10805268, Kanzaki, 2024). A readme file on the web provides the instructions for executing the simulations. All underlying research data can be assessed in the references cited in this paper.

Supplement

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-17-4515-2024-supplement.

Author contributions

YK and CTR designed and implemented the model, with significant contributions from the other authors. IC obtained the mesocosm data. SZ and YK compiled and filtered soil data from the literature. All authors contributed to the experimental design, and YK conducted the experiments and analyzed the results. All authors contributed to the writing of the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We are grateful to two anonymous reviewers for their useful comments and to Sandra Arndt for editorial handling.

Review statement

This paper was edited by Sandra Arndt and reviewed by two anonymous referees.

References

Alva, A. K., Sumner, M. E., and Miller, W. P.: Relationship between ionic strength and electrical conductivity for soil solutions, Soil Sci., 152, 239–242, https://doi.org/10.1097/00010694-199110000-00001, 1991. 

Appelo, C. A. J.: Cation and proton exchange, pH variations, and carbonate reactions in a freshening aquifer, Water Resour. Res., 30, 2793–2805, https://doi.org/10.1029/94WR01048, 1994. 

Aradóttir, E. S. P., Sonnenthal, E. L., and Jónsson, H.: Development and evaluation of a thermodynamic dataset for phases of interest in CO2 mineral sequestration in basaltic rocks, Chem. Geol., 304, 26–38, https://doi.org/10.1016/j.chemgeo.2012.01.031, 2012. 

Beerling, D. J., Kantzas, E. P., Lomas, M. R., Wade, P., Eufrasio, R. M., Renforth, P., Sarkar, B., Andrews, M. G., James, R. H., Pearce, C. R., Mercure, J.-F., Pollitt, H., Holden, P. B., Edwards, N. R., Khanna, M., Koh, L., Quegan, S., Pidgeon, N. F., Janssens, I. A., Hansen, J., and Banwart, S. A.: Potential for large-scale CO2 removal via enhanced rock weathering with croplands, Nature, 583, 242–248, https://doi.org/10.1038/s41586-020-2448-9, 2020. 

Brantley, S. L., Kubicki, J. D., and White, A. F.: Kinetics of Water-Rock Interaction, Springer, ISBN 978-0-387-73562-7, 2008. 

Chiaravalloti, I., Theunissen, N., Zhang, S., Wang, J., Sun, F., Ahmed, A. A., Pihlap, E., Reinhard, C. T., and Planavsky, N. J.: Mitigation of soil nitrous oxide emissions during maize production with basalt amendments, Front. Clim., 5, 1203043, https://doi.org/10.3389/fclim.2023.1203043, 2023. 

Dietzen, C. and Rosing, M. T.: Quantification of CO2 uptake by enhanced weathering of silicate minerals applied to acidic soils, Int. J. Greenhouse Gas Contr., 125, 103872, https://doi.org/10.1016/j.ijggc.2023.103872, 2023. 

Fernández, F. G. and Hoeft, R. G.: Managing soil pH and crop nutrients, in: Illinois Agronomy Handbook, 24, 91–112, 2009. 

Fick, S. E. and Hijmans, R. J.: WorldClim 2: new 1km spatial resolution climate surfaces for global land areas, Int. J. Climatol., 37, 4302–4315, 2017. 

Geibe, C. E., Danielsson, R., van Hees, P. A., and Lundström, U. S.: Comparison of soil solution chemistry sampled by centrifugation, two types of suction lysimeters and zero-tension lysimeters, Appl. Geochem., 21, 2096–2111, https://doi.org/10.1016/j.apgeochem.2006.07.010, 2006. 

Goldberg, R. N., Kishore, N., and Lennen, R. M.: Thermodynamic quantities for the ionization reactions of buffers, J. Phys. Chem. Ref. Data, 31, 231–370, https://doi.org/10.1063/1.1416902, 2002. 

Goll, D. S., Ciais, P., Amann, T., Buermann, W., Chang, J., Eker, S., Hartmann, J., Janssens, I., Li, W., Obersteiner, M., Penuelas, J., Tanaka, K., and Vicca, S: Potential CO2 removal from enhanced weathering by ecosystem responses to powdered rock, Nat. Geosci., 14, 545–549, https://doi.org/10.1038/s41561-021-00798-x, 2021. 

Goulding, K. W. T.: Soil acidification and the importance of liming agricultural soils with particular reference to the United Kingdom, Soil Use Manage., 32, 390–399, https://doi.org/10.1111/sum.12270, 2016. 

Grybos, M., Davranche, M., Gruau, G., Petitjean, P., and Pédrot, M.: Increasing pH drives organic matter solubilization from wetland soils under reducing conditions, Geoderma, 154, 13–19, https://doi.org/10.1016/j.geoderma.2009.09.001, 2009. 

Hamilton, S. K., Kurzman, A. L., Arango, C., Jin, L., and Robertson, G. P.: Evidence for carbon sequestration by agricultural liming, Global Biogeochem. Cycles, 21, GB2021, https://doi.org/10.1029/2006GB002738, 2007. 

IPCC: 2006 IPCC Guidelines for National Greenhouse Gas Inventories, IPCC, ISBN 4-88788-032-4, 2006. 

IPCC: Global Warming of 1.5 °C, IPCC, https://doi.org/10.1017/9781009157940, 2018. 

ISRIC: WISE Soil Property Databases, https://www.isric.org/explore/wise-databases, last access: 19 January 2023, 2022. 

Jones Jr, J. B. (Ed.): Soil Analysis Handbook of Reference Methods, CRC press, https://doi.org/10.1201/9780203739433, 1999. 

Kantzas, E. P., Val Martin, M., Lomas, M. R., Eufrasio, R. M., Renforth, P., Lewis, A. L., Taylor, L. L., Mecure, J. F., Pollitt, H., Vercoulen, P. V., and Vakilifard, N.: Substantial carbon drawdown potential from enhanced rock weathering in the United Kingdom. Nat. Geosci., 15, 382–389, https://doi.org/10.1038/s41561-022-00925-2, 2022. 

Kanzaki, Y.: cdr-laboratory/SCEPTER: v1.0.1 (v1.0.1), Zenodo [code], https://doi.org/10.5281/zenodo.10805268, 2024. 

Kanzaki, Y., Zhang, S., Planavsky, N. J., and Reinhard, C. T.: Soil Cycles of Elements simulator for Predicting TERrestrial regulation of greenhouse gases: SCEPTER v0.9, Geosci. Model Dev., 15, 4959–4990, https://doi.org/10.5194/gmd-15-4959-2022, 2022. 

Kanzaki, Y., Planavsky, N. J., and Reinhard, C. T.: New estimates of the storage permanence and ocean co-benefits of enhanced rock weathering, PNAS Nexus, 2, pgad059, https://doi.org/10.1093/pnasnexus/pgad059, 2023. 

Kelland, M. E., Wade, P. W., Lewis, A. L., Taylor, L. L., Sarkar, B., Andrews, M. G., Lomas, M. R., Cotton, T. A., Kemp, S. J., James, R. H., and Pearce, C. R.: Increased yield and CO2 sequestration potential with the C4 cereal Sorghum bicolor cultivated in basaltic rock dust-amended agricultural soil, Global Change Biol., 26, 3658–3676, https://doi.org/10.1111/gcb.15089, 2020. 

Kissel, D. E. and Sonon, L. S.: Soil Test Handbook for Georgia, University of Georgia College of Agricultural and Environmental Sciences, Athens, GA, 2008. 

Kopittke, P. M., Menzies, N. W., Wang, P., McKenna, B. A., and Lombi, E.: Soil and the intensification of agriculture for global food security, Environ. Int., 132, 105078, https://doi.org/10.1016/j.envint.2019.105078, 2019. 

Köhler, P., Hartmann, J., and Wolf-Gladrow, D. A.: Geoengineering potential of artificially enhanced silicate weathering of olivine, P. Natl. Acad. Sci. USA, 107, 20228–20233, https://doi.org/10.1073/pnas.1000545107, 2010. 

La-Scalea, M. A., Menezes, C. M. S., and Ferreira, E. I.: Molecular volume calculation using AM1 semi-empirical method toward diffusion coefficients and electrophoretic mobility estimates in aqueous solution, J. Mol. Struct.-THEOCHEM, 730, 111–120, https://doi.org/10.1016/j.theochem.2005.05.030, 2005. 

Li, Y.-H. and Gregory, S.: Diffusion of ions in sea water and in deep-sea sediments, Geochim. Cosmochim. Ac., 38, 703–714, https://doi.org/10.1016/0016-7037(74)90145-8, 1974. 

Likhachev, E. R.: Dependence of water viscosity on temperature and pressure, Tech. Phys., 48, 514–515, https://doi.org/10.1134/1.1568496, 2003. 

McGechan, M. B. and Lewis, D. R.: SW – soil and water: sorption of phosphorus by soil, part 1: principles, equations and models, Biosyst. Eng., 82, 1–24, https://doi.org/10.1006/bioe.2002.0054, 2002. 

McLean, E. O.: Soil pH and lime requirement, In: Methods of Soil Analysis: Part 2 Chemical and Microbiological Properties, 9, edited by: Page, A. L., 199–224, https://doi.org/10.2134/agronmonogr9.2.2ed.c12, 1983. 

Miller, R. O. and Kissel, D. E.: Comparison of soil pH methods on soils of North America, Soil Sci. Soc. Am. J., 74, 310–316, https://doi.org/10.2136/sssaj2008.0047, 2010. 

Nambu, K. and Yonebayashi, K.: Acidic properties of dissolved organic matter leached from organic layers in temperate forests, Soil Sci. Plant Nutr., 45, 65–77, https://doi.org/10.1080/00380768.1999.10409324, 1999. 

Nielsen, K. E., Irizar, A., Nielsen, L. P., Kristiansen, S. M., Damgaard, C., Holmstrup, M., Petersen, A. R., and Strandberg, M.: In situ measurements reveal extremely low pH in soil, Soil Biol. Biochem., 115, 63–65, https://doi.org/10.1016/j.soilbio.2017.08.003, 2017. 

Othmer, D. F. and Thakar, M. S.: Correlating diffusion coefficient in liquids, Ind. Eng. Chem., 45, 589–593, https://doi.org/10.1021/ie50519a036, 1953. 

Pan, B., Lam, S. K., Wang, E., Mosier, A., and Chen, D.: New approach for predicting nitrification and its fraction of N2O emissions in global terrestrial ecosystems, Environ. Res. Lett., 16, 034053, https://doi.org/10.1088/1748-9326/abe4f5, 2021. 

Parfitt, R. L., Giltrap, D. J., and Whitton, J. S.: Contribution of organic matter and clay minerals to the cation exchange capacity of soils, Commun. Soil Sci. Plant Anal., 26, 1343–1355, https://doi.org/10.1080/00103629509369376, 1995. 

Parkhurst, D. L. and Appelo, C. A. J.: Description of input and examples for PHREEQC version 3: a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, U. S. Geological Survey, https://doi.org/10.3133/tm6A43, 2013. 

Poggio, L., de Sousa, L. M., Batjes, N. H., Heuvelink, G. B. M., Kempen, B., Ribeiro, E., and Rossiter, D.: SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty, SOIL, 7, 217–240, https://doi.org/10.5194/soil-7-217-2021, 2021. 

Pollyea, R. M. and Rimstidt, J. D.: Rate equations for modeling carbon dioxide sequestration in basalt, Appl. Geochem., 81, 53–62, https://doi.org/10.1016/j.apgeochem.2017.03.020, 2017. 

Ponnamperuma, F. N., Tianco, E. M., and Loy, T. A.: Ionic strengths of the solutions of flooded soils and other natural aqueous solutions from specific conductance, Soil Sci., 102, 408–413, https://doi.org/10.1097/00010694-196612000-00009, 1966. 

Rau, G. H., Knauss, K. G., Langer, W. H., and Caldeira, K.: Reducing energy-related CO2 emissions using accelerated weathering of limestone, Energy, 32, 1471–1477, https://doi.org/10.1016/j.energy.2006.10.011, 2007. 

Reitz, M., Sanford, W. E., Senay, G. B., and Cazenas, J.: Annual estimates of recharge, quick-flow runoff, and evapotranspiration for the contiguous US using empirical regression equations. J. Am. Water Resour. Assoc., 53, 961–983, https://doi.org/10.1111/1752-1688.12546, 2017. 

Rengel, Z. (Ed.): Handbook of Soil Acidity, CRC Press, https://doi.org/10.1201/9780203912317, 2003. 

Robie, R. A. Hemingway, B. S., and Fisher, J. R.: Thermodynamic properties of minerals and related substances at 298.15 K and 1 bar (105 pascals) pressure and at higher temperatures, U. S. Geological Survey, https://doi.org/10.3133/b1452, 1978. 

Schulz, H. D. and Zabel, M.: Marine Geochemistry, Springer, https://doi.org/10.1007/3-540-32144-6, 2006. 

Sikora, F. J.: A buffer that mimics the SMP buffer for determining lime requirement of soil, Soil Sci. Soc. Am. J., 70, 474–486, https://doi.org/10.2136/sssaj2005.0164, 2006. 

Steefel, C. I.: CrunchFlow Software for Modeling Multicomponent Reactive Flow and Transport USER'S MANUAL, 2009. 

Steefel, C. I., Carroll, S., Zhao, P., and Roberts, S.: Cesium migration in Hanford sediment: a multisite cation exchange model based on laboratory transport experiments, J. Contam. Hydrol., 67, 219–246, https://doi.org/10.1016/S0169-7722(03)00033-0, 2003. 

Steiner, Z., Lazar, B., Erez, J., and Turchyn, A. V.: Comparing Rhizon samplers and centrifugation for pore-water separation in studies of the marine carbonate system in sediments, Limnol. Oceanogr.-Meth., 16, 828–839, https://doi.org/10.1002/lom3.10286, 2018. 

Strefler, J., Amann, T., Bauer, N., Kriegler, E., and Hartmann, J.: Potential and costs of carbon dioxide removal by enhanced weathering of rocks, Environ. Res. Lett., 13, 034010, https://doi.org/10.1088/1748-9326/aaa9c4, 2018. 

Swoboda, P., Döring, T. F., and Hamer, M.: Remineralizing soils? The agricultural usage of silicate rock powders: A review, Sci. Total Environ., 807, 150976, https://doi.org/10.1016/j.scitotenv.2021.150976, 2022. 

Taylor, L. L., Quirk, J., Thorley, R. M. S., Kharecha, P. A., Hansen, J., Ridgwell, A., Lomas, M. R., Banwart, S. A., and Beerling, D. J.: Enhanced weathering strategies for stabilizing climate and averting ocean acidification, Nat. Clim. Change, 6, 402–406, https://doi.org/10.1038/NCLIMATE2882, 2016. 

Thomas, G. W.: Soil pH and soil acidity, in: Methods of Soil Analysis: Part 3 Chemical Methods, 5, edited by: Sparks, D. L., Page, A. L., Helmke, P. A., Loeppert, R. H., Soltanpour, P. N., Tabatabai, M. A., Johnston, C. T., and Sumner, M. E., 475–490, https://doi.org/10.2136/sssabookser5.3.c16, 1996.  

Thomas Sims, J.: Lime requirement, in: Methods of Soil Analysis: Part 3 Chemical Methods, 5, edited by: Sparks, D. L., Page, A. L., Helmke, P. A., Loeppert, R. H., Soltanpour, P. N., Tabatabai, M. A., Johnston, C. T., and Sumner, M. E., 491–515, https://doi.org/10.2136/sssabookser5.3.c17, 1996. 

Vakilifard, N., Kantzas, E. P., Edwards, N. R., Holden, P. B., and Beerling, D. J.: The role of enhanced rock weathering deployment with agriculture in limiting future warming and protecting coral reefs, Environ. Res. Lett., 16, 094005, https://doi.org/10.1088/1748-9326/ac1818, 2021. 

Van der Zee, S. E. and van Riemsdijk, W. H.: Model for long-term phosphate reaction kinetics in soil, J. Environ. Qual., 17, 35–41, https://doi.org/10.2134/jeq1988.00472425001700010005x, 1988. 

Walkinshaw, M., O'Geen, A. T., and Beaudette, D. E.: Soil Properties, California Soil Resource Lab, http://casoilresource.lawr.ucdavis.edu/soil-properties/ (last access: 22 May 2023), 2022. 

Wang, Y., Mao, J., Jin, M., Hoffman, F. M., Shi, X., Wullschleger, S. D., and Dai, Y.: Development of observation-based global multilayer soil moisture products for 1970 to 2016, Earth Syst. Sci. Data, 13, 4385–4405, https://doi.org/10.5194/essd-13-4385-2021, 2021. 

Zhang, S., Planavsky, N. J., Katchinoff, J., Raymond, P. A., Kanzaki, Y., Reershemius, T., and Reinhard, C. T.: River chemistry constraints on the carbon capture potential of surficial enhanced rock weathering, Limnol. Oceanogr., 67. S148–S157, https://doi.org/10.1002/lno.12244, 2022. 

Download
Short summary
Soil pH is one of the most commonly measured agronomical and biogeochemical indices, mostly reflecting exchangeable acidity. Explicit simulation of both porewater and bulk soil pH is thus crucial to the accurate evaluation of alkalinity required to counteract soil acidification and the resulting capture of anthropogenic carbon dioxide through the enhanced weathering technique. This has been enabled by the updated reactive–transport SCEPTER code and newly developed framework to simulate soil pH.