Soil Cycles of Elements simulator for Predicting TERrestrial regulation of greenhouse gases: SCEPTER v0.9

The regulation of anthropogenic carbon dioxide (CO₂) is an urgent issue – continuously increasing atmospheric CO₂ from burning fossil fuels is leading to significant warming and acidification of the surface ocean. Timely and effective measures to curb CO₂ increases are thus needed in order to mitigate the potential degradation of natural ecosystems, food security, and livelihood caused by anthropogenic release of CO₂. Enhanced rock weathering (ERW) on croplands and hinterlands may be one of the most economically and ecologically effective ways to sequester CO₂ from the atmosphere, given that these soil environments generally favor mineral dissolution and because amending soils with crushed rock can result in a number of co-benefits to plant growth and crop yield. However, robust quantitative evaluation of CO₂ capture by ERW in terrestrial soil systems that can lead to coherent policy implementation will require an ensemble of traceable mechanistic models that are optimized for simulating ERW in managed systems. Here, we present a new 1D reactive transport model – SCEPTER. The model is designed to (1) mechanistically simulate natural weathering, including dissolution/precipitation of minerals along with uplift/erosion of solid phases, advection plus diffusion of aqueous phases and diffusion of gas phases, (2) allow targeted addition of solid phases at the soil–atmosphere interface, including multiple forms of organic matter (OM) and crushed mineral/rock feedstocks, (3) implement a range of soil mixing regimes as catalyzed by soil surface fauna (e.g., bioturbation) or humans (e.g., various forms of tilling), and (4) enable calculation of solid mineral surface area based on controlled initial particle size distributions coupled to a shrinking core framework. Here we describe the model structure and intrinsic thermodynamic/kinetic data, provide a series of idealized simulations to demonstrate the basic behavior of the code, and evaluate the computational and mechanistic performance of the model against observational data. We also provide selected example applications to highlight model features particularly useful for future prediction of CO₂ sequestration by ERW in soil systems.

1 Introduction

Continuously increasing emissions of CO₂ and other greenhouse gases (GHGs) from fossil fuel consumption and land use change have resulted in large changes to atmospheric chemistry and global temperature since the beginning of the industrial era (e.g., Le Quéré et al., 2018) and are expected to lead to significant climate perturbation and environmental degradation in the coming century (IPCC, 2006). Although reducing anthropogenic GHG emissions must be the linchpin for mitigating degradation of surface environments (e.g., Rogelj et al., 2018), all current pathways delineated by the Intergovernmental Panel on Climate Change (IPCC, 2006, 2018) as potentially limiting global warming to below 1.5 ^∘C by 2100 require carbon dioxide removal (CDR) on the order of 10²–10³ GtCO₂ (GtCO₂=10¹⁵ gCO₂) over the course of the next century (IPCC, 2018), and less severe CO₂ regulation trajectories will also likely require active CDR. As a result, various modes of CO₂ capture will likely be critical for achieving climate targets set by the international community (e.g., Fuss et al., 2014; Gasser et al., 2015; IPCC, 2018).

Enhanced rock weathering (ERW) at Earth's surface is one potential means of executing CDR on a gigaton scale (e.g., Köhler et al., 2010; Taylor et al., 2016; Beerling et al., 2020). Broadly, this class of CDR strategies involves the sequestration of atmospheric CO₂ as dissolved inorganic carbon (DIC) through reaction with silicate or carbonate minerals. In principle, this could be achieved across a range of marine (Rau et al., 2007; Köhler et al., 2013; Renforth and Henderson, 2017) and terrestrial (Köhler et al., 2010; Hartmann et al., 2013; Beerling et al., 2020) environments. However, terrestrial ERW in agricultural settings received particular recent attention as a potentially cost-effective strategy for CDR with a range of possible co-benefits, including increasing crop yields and neutralization of ongoing surface ocean acidification (Minx et al., 2018; Strefler et al., 2018; Beerling et al., 2020). Numerical tools will be essential to chart a path forward with CO₂ capture by ERW given environmental and economic constraints. Following more developed modeling fields (e.g., CMIP6; e.g., Liddicoat et al., 2021), the most robust estimates of CDR potential and operational cost will require an ensemble of traceable models that are optimized for addressing the feedbacks between soil systems and ERW intervention. Such model developments have become more readily feasible given the successful developments and applications of natural weathering models (e.g., Sverdrup and Warfvinge, 1993; Sverdrup et al., 1995; Goddéris et al., 2006, 2013; Maher et al., 2009; Brantley and Lebedeva, 2011; Beaulieu et al., 2012; Moore et al., 2012; Roland et al., 2013; Steefel et al., 2015; Zhi et al., 2022), whose theoretical frameworks and/or numerical schemes can be utilized as essential building blocks of a comprehensive modelling framework for ERW.

To help facilitate robust prediction of the CO₂ capture efficiency, environmental impacts, and operational costs of ERW in terrestrial soil systems, we have developed a new 1D reactive-transport model – SCEPTER. The model is designed for quantification of interactions between accelerated dissolution of applied rock/mineral powders and background natural weathering, including soil respiration and particle mixing by surface soil fauna. Soil mixing is implemented by adapting a transition matrix method, which allows the user to define their own transfer function for desired mixing regime. We track surface area differences between background rocks/soils and applied mineral/rock powders by tracking porosity evolution caused by addition of mineral/rock powders and/or by tracking size distributions of particles with an adapted shrinking particle model. Increased surface areas from milled grains can be retained or annealed with reaction progress. First, we describe the theoretical background and numerical implementation of SCEPTER (Sect. 2). Then, a series of idealized simulations are presented in order to illustrate the basic capability of the code as a natural weathering simulator and the utility of model components specifically designed to interrogate the impact of ERW (Sect. 3.1). Finally, we compare model results to soil depth profiles of pH and OM from site-specific US soil data (Sect. 3.2).

Table 1 Thermodynamic data of solid species^a.

^a The thermodynamic constant for species θ (K_θ) is calculated as $K_{\mathit{\theta }}=K_{\mathit{\theta }}^{\mathrm{ref}}$exp[−ΔH_θ{ 1/($T_{\mathrm{C}}+\mathrm{273.15})-\mathrm{1}/(T_{\mathrm{C}}^{\mathrm{ref}}+\mathrm{273.15}\left)\mathit{\right\}}/\mathrm{\Re }$], where T_C is temperature in ^∘C and ℜ is the gas constant in units of kJ mol⁻¹ K⁻¹ ($\mathrm{\Re }=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹). Thermodynamic constants for SOMs and pyrite are not defined here as the saturation states for these species are kinetically defined (see Sect. 2.2.2 for more details): ${\mathrm{\Omega }}_{\mathrm{py}}=\mathrm{1}-p_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{0.5}}$ (pyrite) and ${\mathrm{\Omega }}_{\mathrm{g}}=\mathrm{1}-p_{{\mathrm{O}}_{\mathrm{2}}}/(p_{{\mathrm{O}}_{\mathrm{2}}}+K_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}$) (SOM Classes 1 to 3), where $p_{{\mathrm{O}}_{\mathrm{2}}}$ is the partial pressure of soil O₂ (atm) and $K_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}$ is the Michaelis constant for aerobic respiration (0.121 atm; see Kanzaki and Kump, 2017, and references therein). ^b Units change with y depending on minerals. ^c From Robie et al. (1978) except for illite and beidellites (Wolery and Jove-Colon, 2004), hedenbergite (calculated from the molar weight and density in https://www.mindat.org, last access: 15 June 2022) and SOMs (Mayer et al., 2004). Solid solutions are calculated from molar volumes of endmember minerals. ^d From Robie et al. (1978). When not available in Robie et al. (1978), calculated from chemical composition. ^e (1) From phreeqc.dat available in PHREEQC v.3 (Parkhurst and Appelo, 2013). (2) From minteq.v4.dat available in PHREEQC v.3 (Parkhurst and Appelo, 2013). (3) Data by Sugimori et al. (2012) from the Lawrence Livermore National Laboratory thermodynamic dataset, thermo.com.v8r6+.dat (Delany and Lundeen, 1990). (4) Summarized data by Kanzaki and Murakami (2018) from thermo.com.v8r6+.dat (Delany and Lundeen, 1990). (5) Data for disordered dolomite from minteq.v4.dat available in PHREEQC v.3 (Parkhurst and Appelo, 2013). (6) Wilkin and Barnes (1998). (7) Calculated assuming an ideal solid solution after Gislason and Arnorsson (1993). (8) Not defined (see caption a and Sect. 2.2.2).

Table 2 Dissolution/precipitation kinetic data of solid species^a.

^a The rate constant for species θ (mol m⁻² s⁻¹) is calculated as k_θ= [H${}^{+}{]}^{n_{\mathrm{acid}}}{k}_{\mathrm{acid}}^{\mathrm{ref}}$exp[$-E_{\mathrm{acid}}^{\mathrm{app}}\mathit{\{}$ 1/($T_{\mathrm{C}}+\mathrm{273.15})-\mathrm{1}/(T_{\mathrm{C}}^{\mathrm{ref}}+\mathrm{273.15}\left)\mathit{\right\}}/\mathrm{\Re }]+k_{\mathrm{ntrl}}^{\mathrm{ref}}$exp[$-E_{\mathrm{ntrl}}^{\mathrm{app}}\mathit{\{}$ 1/($T_{\mathrm{C}}+\mathrm{273.15})-\mathrm{1}/(T_{\mathrm{C}}^{\mathrm{ref}}+\mathrm{273.15}\left)\mathit{\right\}}/\mathrm{\Re }]+$ [H${}^{+}{]}^{n_{\mathrm{alk}}}{k}_{\mathrm{alk}}^{\mathrm{ref}}$exp[$-E_{\mathrm{alk}}^{\mathrm{app}}\mathit{\{}$ 1/(T_C+273.15) -1/($T_{\mathrm{C}}^{\mathrm{ref}}+\mathrm{273.15}\left)\mathit{\right\}}$ /ℜ], where T_C is temperature in ^∘C and ℜ is the gas constant in units of kJ mol⁻¹ K⁻¹ ($\mathrm{\Re }=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹), except for carbonates and SOMs. For carbonates, [H⁺] in the last term of the above equation (“base” mechanism) needs to be replaced by $p_{{\mathrm{CO}}_{\mathrm{2}}}$ (atm), the partial pressure of soil CO₂ (“carbonate” mechanism) according to Palandri and Kharaka (2004), and for SOMs, rate constants are assumed to be independent of pH and temperature. ^b Solid species that are unlikely to precipitate are not allowed to precipitate even when porewater is supersaturated with respect to them in the model, denoted with “n” in this column. Those that can precipitate are denoted with “y”. ^c (1) Palandri and Kharaka (2004). (2) “Acid” and “base” mechanisms from Brantley et al. (2008) (K_w is the water dissociation constant) and the “neutral” mechanism from Palandri and Kharaka (2004). Activation energy for the “neutral” mechanism is assumed to be applicable to other mechanisms. (3) Data for calcite are assumed. (4) The “base” mechanism is replaced by the “carbonate” mechanism (see caption a for details). (5) Data for disordered dolomite are assumed. (6) Data for enstatite are assumed. (7) Data for augite are assumed. (8) Data for nepheline are assumed (cf. Ragnarsdóttir, 1993). (9) Data for smectite are assumed (cf. Bibi et al., 2011). (10) $H_{{\mathrm{O}}_{\mathrm{2}}}$ is the solubility of O₂ (mol L⁻¹ atm⁻¹). Dependences on pH and O₂ are from Williamson and Rimstidt (1994) and apparent activation energy from McKibben and Barnes (1986). (11) See caption a. Corresponding to turnover times of 1, 8 and 1000 years for decomposition of SOM Classes 1, 2 and 3, respectively (cf. Chen et al., 2010).

2 Model description

2.1 Overview

The basic framework of SCEPTER is derived from previous models designed to simulate reactive transport and weathering in natural soil systems, with a focus on pyrite and organic matter oxidation (Kanzaki and Kump, 2017) and silicate mineral transformation (Kanzaki et al., 2020). SCEPTER can currently simulate up to 39 minerals and different classes of soil organic matter (SOM; currently configured with three classes), 10 independent aqueous species along with 48 dependent aqueous species, and 4 independent gas species (Tables 1–5). Reaction kinetics, especially dissolution/precipitation, are explicitly implemented for individual solid species and are fully coupled with the temporal and spatial evolution of aqueous and gas species (e.g., Table 2). One can further add/remove a series of additional reactions to the system associated with solid, aqueous, and/or gaseous species (referred to as “extra” reactions, e.g., Table 6). Particular features of SCEPTER designed to interrogate ERW in terrestrial soil systems include (1) implementation of bio-mixing in the upper soil layers using a transition matrix approach, (2) time-dependent application of crushed rock feedstock to the soil surface, (3) options for time-varying boundary conditions including porosity and advection rates of solids and porewater, and (4) a range of user options for calculation and specification of the surface area of solid species. These features of the model are specifically designed for aiding in robust prediction of enhanced weathering on terrestrial ecosystems, including croplands and hinterlands (e.g., Hartmann et al., 2013; Taylor et al., 2016; Beerling et al., 2020; Goll et al., 2021).

In the following we describe the theoretical framework of the model (Sect. 2.2), its numerical implementation (Sect. 2.3), and the user input in relation to model initialization and boundary conditions (Sect. 2.4). The model code is written in Fortran90 (see the code availability section).

Table 3 Thermodynamic data for aqueous species^a.

^a Thermodynamic constant (K_aq) is calculated as $K_{\mathrm{aq}}=K_{\mathrm{aq}}^{\mathrm{ref}}\exp [-\mathrm{\Delta }{H}_{\mathrm{aq}}\mathit{\{}\mathrm{1}/(T_{\mathrm{C}}+\mathrm{273.15})$ $-\mathrm{1}/(T_{\mathrm{C}}^{\mathrm{ref}}+\mathrm{273.15})\mathit{\}}/\mathrm{\Re }]$, where T_C is temperature in ^∘C and ℜ is the gas constant in units of kJ mol⁻¹ K⁻¹ ($\mathrm{\Re }=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹). ^b (1) From phreeqc.dat available in PHREEQC v.3 (Parkhurst and Appelo, 2013). (2) Maggi et al. (2008). Temperature dependence is assumed to be 0. (3) Kanzaki and Murakami (2016). (4) Thermodynamic constant for Fe${}^{\mathrm{2}+}+{\mathrm{HCO}}_{\mathrm{3}}^{-}=\mathrm{Fe}{\mathrm{HCO}}_{\mathrm{3}}^{+}$ from Kanzaki and Murakami (2016) divided by that for ${\mathrm{HCO}}_{\mathrm{3}}^{-}={\mathrm{H}}^{+}+{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (Table 4). (5) Thermodynamic constant for Na${}^{+}+{\mathrm{HCO}}_{\mathrm{3}}^{-}=$ NaHCO₃ from phreeqc.dat available in PHREEQC v.3 (Parkhurst and Appelo, 2013) divided by that for ${\mathrm{HCO}}_{\mathrm{3}}^{-}=$ H${}^{+}+$ ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (Table 4). (6) Kanzaki and Murakami (2015).

Table 4 Thermodynamic data for gaseous species^a.

^a The thermodynamic constant (K_gas) is calculated as $K_{\mathrm{gas}}=K_{\mathrm{gas}}^{\mathrm{ref}}\exp [-\mathrm{\Delta }{H}_{\mathrm{gas}}\mathit{\{}\mathrm{1}/(T_{\mathrm{C}}+\mathrm{273.15})$ $-\mathrm{1}/(T_{\mathrm{C}}^{\mathrm{ref}}+\mathrm{273.15})\mathit{\}}/\mathrm{\Re }]$, where T_C is temperature in ^∘C and ℜ is the gas constant in kJ mol⁻¹ K⁻¹ ($\mathrm{\Re }=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹). ^b (1) Kanzaki and Murakami (2016). (2) Kanzaki and Murakami (2015). (3) From wateq4f.dat available in PHREEQC v.3 (Parkhurst and Appelo, 2013). (4) Weiss and Price (1980).

2.2 Theoretical framework

2.2.1 Tracing solid, aqueous and gas species in soils

SCEPTER is based on previous models designed to represent fundamental aspects of natural weathering (e.g., Bolton et al., 2006; Brantley and Lebedeva, 2011; Li et al., 2014; Steefel et al., 2015; Kanzaki et al., 2020). Solid minerals are transported upward by continental uplift and eroded at the surface while reacting with solutes in porewater whose transport is dominated by molecular diffusion plus advection caused by downward infiltration. Some solutes are derived from the gas phase present in soil pore space whose composition is controlled by diffusive transport plus consumption/production through reactions within soil, including dissolution into and degassing from porewater. These reactions and transport of solid, aqueous and gaseous phases are simulated within a 1D model soil domain where materials are exchangeable at the top and bottom boundaries.

In addition to the above basic description of the model's framework as a natural soil/rock weathering simulator, SCEPTER implements bio-mixing of solid particles in soils and a rain of solid particles onto the soil surface (Sect. 2.1). Including these additional supply and mixing effects, a solid species θ is simulated according to the following generalized equation (cf. Boudreau, 1997; Kanzaki et al., 2020, 2021):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial m_{\mathit{\theta }}}{\partial t}}& ={\displaystyle \frac{\partial w{m}_{\mathit{\theta }}}{\partial z}}-R_{\mathit{\theta }}-\sum _{\mathit{\kappa }}^{n_{\mathrm{xrxn}}}{\mathit{\gamma }}_{\mathit{\kappa },\mathit{\theta }}{R}_{\mathit{\kappa }}+J_{\mathit{\theta }}\\ & -m_{\mathit{\theta }}\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 }}{m}_{\mathit{\theta }}\left(z^{\prime }\right)E_{\mathit{\theta }}(z^{\prime },z)\mathrm{d}{z}^{\prime },\end{array}\end{array}$$

where m_θ is the mole amount of solid species θ per unit bulk soil/rock volume (mol m⁻³), t is time (yr), z is the depth of the weathering profile (m), w is the advection rate of solid phases (m yr⁻¹), R_θ and J_θ are the net dissolution and rain rates of solid species θ, respectively (mol m⁻³ yr⁻¹), R_κ denotes the rate of the κth extra reaction (mol m⁻³ yr⁻¹) whose stoichiometry with respect to consumption of θ is defined by γ_κ,θ, n_xrxn is the total number of extra reactions, E_θ (z, z^′) is the rate of particle transfer between locations at z and z^′ by bio-mixing (m⁻¹ yr⁻¹) and z_ml is the mixed layer depth (m) within which bio-mixing occurs. The exchange function E_θ (z, z^′) expresses the bio-mixing rate in a generalized continuous form, whose discretized counterpart can be formulated with a transition matrix (e.g., Boudreau, 1997; Shull, 2001; Kanzaki et al., 2021). Various bio-mixing styles and corresponding transition matrices are elaborated in Sect. 2.2.5. As described above, SCEPTER also allows flexible addition of many additional reactions (with example “extra” reactions given in Table 6 and discussed in Sect. 2.2.2) beyond those parameterized with the thermodynamic and kinetic constants in Tables 1 and 2. See the following subsections for more details on individual transport and reaction terms.

Table 5 Diffusion coefficient for aqueous and gaseous species^a.

^a The diffusion coefficient is calculated as $D=D^{\mathrm{ref}}\exp [-E_{\mathrm{diff}}^{\mathrm{app}}\mathit{\{}\mathrm{1}/(T_{\mathrm{C}}+\mathrm{273.15})$ $-\mathrm{1}/(T_{\mathrm{C}}^{\mathrm{ref}}+\mathrm{273.15})\mathit{\}}/\mathrm{\Re }]$, where T_C is temperature in ^∘C and ℜ is the gas constant in units of kJ mol⁻¹ K⁻¹ ($\mathrm{\Re }=\mathrm{8.314}\times {\mathrm{10}}^{-\mathrm{3}}$ kJ mol⁻¹ K⁻¹), except for NH₃ gas. A power function of temperature is assumed for NH₃ gas from Massman (1998): D=624{ (T_C+273.15)/273.15}^1.81 (m² yr⁻¹). ^b (1) From Li and Gregory (1974). (2) From Rebreanu et al. (2008). (3) Kanzaki and Murakami (2016). (4) A value of 0.14 cm² s⁻¹ is assumed (e.g., Pritchard and Currie, 1982), and apparent activation energy is assumed to be the same as that for O₂ gas. (5) Massman (1998). See caption a. (6) Assumed to be the same as CO₂ diffusion (e.g., Pritchard and Currie, 1982). (7) The diffusion coefficient of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ is assumed to be the same as that of aqueous CO₂. (8) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ diffusion from Schulz and Zabel (2006). (9) From Schulz and Zabel (2006).

Table 6 Extra reactions and their kinetic laws.

^a Given in units of mol m⁻³ yr⁻¹. ^b (1) Singer and Stumm (1970). (2) After Maggi et al. (2008). A linear pH function is fitted to an exponential function.

When an element with a given redox state dissolves in porewater and is not sourced from the soil atmosphere, the element is regarded as an aqueous species in the model. The total amount of all dissolved forms of the element per unit porewater volume is traced as an independent variable for modeling solutes in porewater. Specific chemical forms (e.g., hydrolyzed forms and complexes with other ions) are assumed to be in equilibrium, and their concentrations are calculated based on the tracked total concentrations of individual dissolved elements and thermodynamics of association/dissociation reactions (Table 3, Sect. 2.2.2). For each dissolved element, all associated chemical forms are assumed to be transported uniformly (i.e., the same diffusion coefficients are applied to various aqueous forms of a dissolved element; Table 5, Sect. 2.2.4), and thus the governing equation for a dissolved element is given as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\varphi }\mathit{\sigma }\mathrm{\ell }{c}_{\mathit{\varsigma }}}{\partial t}}& =-{\displaystyle \frac{\partial \mathit{\varphi }\mathit{\sigma }\mathrm{\ell }v{c}_{\mathit{\varsigma }}}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\varphi }\mathit{\sigma }\mathrm{\ell }{\mathit{\tau }}_{\mathrm{aq}}{D}_{\mathit{\varsigma }}{\displaystyle \frac{\partial c_{\mathit{\varsigma }}}{\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 }},\end{array}\end{array}$$

where c_ς is the porewater concentration of dissolved element ς (mol L⁻¹), ϕ is the porosity, σ is the water saturation ratio, ℓ is a unit conversion factor (10³ L m⁻³), v is the porewater advection rate (m yr⁻¹), τ_aq is the tortuosity factor for solute diffusion in porewater, D_ς is the diffusion coefficient of dissolved element ς (m² yr⁻¹), n_sld is the total number of simulated solid species, γ_θ,ς is the mole amount of ς released upon dissolution of 1 mole of mineral θ and γ_κ,ς is the stoichiometry of ς production in the κth extra reaction.

A gas species is introduced into the model when it is produced/consumed by aqueous and/or solid species. In the current version of SCEPTER, soil CO₂, O₂, NH₃ and N₂O can be included as gas species. The independent variable is taken to be the soil partial pressure for a given gas species, and concentrations of all dissolved forms derived from the gas species are taken to be dependent variables. Transport occurs via diffusion for a gas species and via diffusion plus porewater advection for its dependent dissolved forms (Sect. 2.1). The governing equation for a gas species is accordingly given by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial {\mathit{\alpha }}_{\mathit{\epsilon }}{p}_{\mathit{\epsilon }}}{\partial t}}& =-{\displaystyle \frac{\partial \mathit{\varphi }\mathit{\sigma }\mathrm{\ell }v{H}_{\mathit{\epsilon }}{p}_{\mathit{\epsilon }}}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(D_{\mathit{\epsilon }}^{\mathrm{eff}}{\displaystyle \frac{\partial p_{\mathit{\epsilon }}}{\partial z}}\right)\\ & +\sum _{\mathit{\theta }}^{n_{\mathrm{sld}}}{\mathit{\gamma }}_{\mathit{\theta },\mathit{\epsilon }}{R}_{\mathit{\theta }}+\sum _{\mathit{\kappa }}^{n_{\mathrm{xrxn}}}{\mathit{\gamma }}_{\mathit{\kappa },\mathit{\epsilon }}{R}_{\mathit{\kappa }},\end{array}\end{array}$$

where p_ε is the partial pressure of soil gas ε (atm), γ_θ,ε is the mole amount of ε released upon dissolution of 1 mole of solid species θ, H_ε is the total solubility of ε (mol L⁻¹ atm⁻¹) and γ_κ,ε is the stoichiometry of ε production in the κth extra reaction. In Eq. (3), α_ε and $D_{\mathit{\epsilon }}^{\mathrm{eff}}$ are the unit conversion factor (mol m⁻³ atm⁻¹) and effective diffusion coefficient (mol atm⁻¹ m⁻¹ yr⁻¹) for soil gas ε, respectively, to include the reactive transport of dependent dissolved forms of ε (e.g., Elberling and Nicholson, 1996):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\alpha }}_{\mathit{\epsilon }}=\mathit{\eta }\mathit{\varphi }(\mathrm{1}-\mathit{\sigma })\mathrm{\ell }+\mathit{\varphi }\mathit{\sigma }\mathrm{\ell }{H}_{\mathit{\epsilon }},\text{(5)}& {\displaystyle }{D}_{\mathit{\epsilon }}^{\mathrm{eff}}=\mathit{\eta }\mathit{\varphi }(\mathrm{1}-\mathit{\sigma })\mathrm{\ell }{\mathit{\tau }}_{\mathrm{gas}}{D}_{\mathit{\epsilon }}^{\mathrm{gas}}+\mathit{\varphi }\mathit{\sigma }\mathrm{\ell }{H}_{\mathit{\epsilon }}{\mathit{\tau }}_{\mathrm{aq}}{D}_{\mathit{\epsilon }}^{\mathrm{aq}},\end{array}$$

where η is the unit conversion from atm to molarity ($={\mathrm{\Re }}^{-\mathrm{1}}{T}^{-\mathrm{1}}$, where ℜ is the gas constant, 0.08205 L atm mol⁻¹ K⁻¹, and T is the temperature in K), τ_gas is the tortuosity factor for gas diffusion in pore air space and $D_{\mathit{\epsilon }}^{\mathrm{gas}}$ and $D_{\mathit{\epsilon }}^{\mathrm{aq}}$ are the diffusion coefficients (m² yr⁻¹) for the gas and aqueous phases of ε, respectively. Individual reaction and transport terms are discussed below.

2.2.2 Reactions

Dissolution/precipitation of a solid species θ (R_θ) is formulated as a function of the saturation state of the species in porewater (Ω_θ), rate coefficient (k_θ, mol m⁻² yr⁻¹) and surface area of the species per unit bulk soil/rock volume available to porewater (S_θ, m² m⁻³):

$$\begin{array}{}\text{(6)}& R_{\mathit{\theta }}=S_{\mathit{\theta }}{k}_{\mathit{\theta }}(\mathrm{1}-{\mathrm{\Omega }}_{\mathit{\theta }}).\end{array}$$

Here, k_θ is a function of porewater pH and soil CO₂ (for carbonates) according to Palandri and Kharaka (2004) (Table 2), and S_θ is related to m_θ, with a number of potential scaling options (see Sect. 2.2.6). The saturation state Ω_θ is calculated based on thermodynamic data for solid species (Table 1). When decomposition of solid species occurs as a redox reaction, e.g., pyrite oxidation and SOM decomposition, Ω_θ is defined kinetically. For instance, pyrite oxidation proceeds according to, e.g., Williamson and Rimstidt (1994):

$$\begin{array}{}\text{(7)}& R_{\mathrm{py}}=S_{\mathrm{py}}{k}_{\mathrm{py}}{p}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{0.5}},\end{array}$$

where $p_{{\mathrm{O}}_{\mathrm{2}}}$ is the partial pressure of soil O₂ (atm). We conventionally define the saturation state for pyrite (Ω_py) as

$$\begin{array}{}\text{(8)}& {\mathrm{\Omega }}_{\mathrm{py}}=\mathrm{1}-p_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{0.5}},\end{array}$$

so that Eq. (6) can be applied to all solid species (Tables 1 and 2). In addition to the dissolution/precipitation reactions specific to individual solid species, one can add extra reactions to the system whose kinetics are explicitly considered. Currently implemented extra reactions and their kinetic laws are listed in Table 6.

All dependent aqueous species are calculated assuming that they are in equilibrium and satisfy the mass balance and charge balance in porewater at any model time step and at all model depths (cf. Steefel et al., 2015). Generally, the mass balance for dependent aqueous species derived from a given dissolved element is given as

$$\begin{array}{}\text{(9)}& c_{\mathit{\varsigma }}=\sum _i^{n_{\mathit{\varsigma }}}{c}_{\mathit{\varsigma }}^i,\end{array}$$

where $c_{\mathit{\varsigma }}^i$ is the porewater concentration of the ith dependent aqueous species derived from dissolved element ς (mol L⁻¹) and n_ς is the total number of dependent aqueous species of ς. Note that the current version of the model does not include any polymers (e.g., Si₂O(OH)₆), as aqueous species and thus 1 mole of dependent aqueous species of dissolved element ς contains 1 mole of ς as assumed in Eq. (9). We define the first dependent species of a given dissolved element as the free dissolved form, except for Si, whose first dependent species is defined as H₄SiO₄, and the other species as products formed via reactions of the first dependent species with other aqueous species in porewater, given, respectively, by

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{c}_{\mathit{\varsigma }}^{\mathrm{1}}& =c_{\mathit{\varsigma }}\left\{\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 }}}\right.\\ & {\left.\times \prod _{\mathit{\epsilon }}^{n_{\mathrm{gas}}}(p_{\mathit{\epsilon }}{)}^{{\mathit{\gamma }}_{\mathit{\varsigma },i,\mathit{\epsilon }}}\right\}}^{-\mathrm{1}},\end{array}\text{(11)}& {\displaystyle }\begin{array}{rl}{c}_{\mathit{\varsigma }}^{i\ne \mathrm{1}}& =c_{\mathit{\varsigma }}^{\mathrm{1}}{K}_{\mathit{\varsigma },i\ne \mathrm{1}}\left[{\mathrm{H}}^{+}{]}^{{\mathit{\gamma }}_{\mathit{\varsigma },i\ne \mathrm{1},\mathrm{p}}}\prod _{{\mathit{\varsigma }}^{\prime }\ne \mathit{\varsigma }}^{n_{\mathrm{aq}}}\right(c_{{\mathit{\varsigma }}^{\prime }}^{\mathrm{1}}{)}^{{\mathit{\gamma }}_{\mathit{\varsigma },i\ne \mathrm{1},{\mathit{\varsigma }}^{\prime }}}\\ & \times \prod _{\mathit{\epsilon }}^{n_{\mathrm{gas}}}(p_{\mathit{\epsilon }}{)}^{{\mathit{\gamma }}_{\mathit{\varsigma },i\ne \mathrm{1},\mathit{\epsilon }}},\end{array}\end{array}$$

where K_ς,i is the thermodynamic constant for production of the ith dependent aqueous species of dissolved element ς, [H⁺] is the porewater H⁺ concentration (mol L⁻¹), ${\mathit{\gamma }}_{\mathit{\varsigma },i,\mathrm{p}}$, ${\mathit{\gamma }}_{\mathit{\varsigma },i,{\mathit{\varsigma }}^{\prime }}$ and ${\mathit{\gamma }}_{\mathit{\varsigma },i,\mathit{\epsilon }}$ are the stoichiometry of H⁺, dissolved element ς^′ and gas species ε, respectively, in the reaction that produces the ith dependent aqueous species of ς, and n_aq and n_gas are the total numbers of independent aqueous and gas species, respectively.

Any aqueous species derived from a gas species is also assumed to be in equilibrium with the gas in soil pore space:

$$\begin{array}{}\text{(12)}& c_{\mathit{\epsilon }}^j=p_{\mathit{\epsilon }}{K}_{\mathrm{H},\mathit{\epsilon }}{K}_{\mathit{\epsilon },j}\left[{\mathrm{H}}^{+}{]}^{{\mathit{\gamma }}_{\mathit{\epsilon },j,\mathrm{p}}}\prod _{\mathit{\varsigma }}^{n_{\mathrm{aq}}}\right(c_{\mathit{\varsigma }}^{\mathrm{1}}{)}^{{\mathit{\gamma }}_{\mathit{\epsilon },j,\mathit{\varsigma }}},\end{array}$$

where $c_{\mathit{\epsilon }}^j$ is the porewater concentration of the jth dependent aqueous species derived from soil gas ε (mol L⁻¹), K_H,ε is Henry's constant for soil gas ε (mol L⁻¹ atm⁻¹), K_ε,j is the thermodynamic constant for production of jth dependent aqueous species of ε, and ${\mathit{\gamma }}_{\mathit{\epsilon },j,\mathrm{p}}$ and ${\mathit{\gamma }}_{\mathit{\epsilon },j,\mathit{\varsigma }}$ are the stoichiometry of H⁺ and the dissolved element ς, respectively, in the reaction that produces the jth dependent aqueous species of ε. The total solubility of ε (H_ε, Eq. 3) is then given by

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{H}_{\mathit{\epsilon }}\equiv \sum _j^{n_{\mathit{\epsilon }}}{c}_{\mathit{\epsilon }}^j/p_{\mathit{\epsilon }}& =K_{\mathrm{H},\mathit{\epsilon }}\sum _j^{n_{\mathit{\epsilon }}}{K}_{\mathit{\epsilon },j}[{\mathrm{H}}^{+}{]}^{{\mathit{\gamma }}_{\mathit{\epsilon },j,\mathrm{p}}}\\ & \times \prod _{\mathit{\varsigma }}^{n_{\mathrm{aq}}}(c_{\mathit{\varsigma }}^{\mathrm{1}}{)}^{{\mathit{\gamma }}_{\mathit{\epsilon },j,\mathit{\varsigma }}}.\end{array}\end{array}$$

Finally, porewater pH is obtained based on charge balance:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\left[{\mathrm{H}}^{+}\right]-\left[{\mathrm{OH}}^{-}\right]& +\sum _{\mathit{\varsigma }}^{n_{\mathrm{aq}}}\sum _i^{n_{\mathit{\varsigma }}}{Z}_{\mathit{\varsigma }}^i{c}_{\mathit{\varsigma }}^i\\ & +\sum _{\mathit{\epsilon }}^{n_{\mathrm{gas}}}\sum _j^{n_{\mathit{\epsilon }}}{Z}_{\mathit{\epsilon }}^j{c}_{\mathit{\epsilon }}^j(\ne c_{{\mathit{\varsigma }}^{\prime }}^{i^{\prime }})=\mathrm{0},\end{array}\end{array}$$

where $Z_{\mathit{\varsigma }}^i$ and $Z_{\mathit{\epsilon }}^j$ are the charges of the ith and jth dependent aqueous species derived from aqueous species ς and gas species ε, respectively, and [OH⁻] is the porewater OH⁻ concentration (mol L⁻¹). Note that double counting of dependent aqueous species is avoided in Eq. (14). Porewater concentrations of H⁺ and OH⁻ are related to one another via the thermodynamic constant of water dissociation K_w (mol² L⁻²):

$$\begin{array}{}\text{(15)}& K_{\mathrm{w}}=\left[{\mathrm{H}}^{+}\right]\left[{\mathrm{OH}}^{-}\right].\end{array}$$

Once p_ε and c_ς are known, Eqs. (9)–(15) can be solved to obtain porewater concentrations of all dependent aqueous species including [H⁺] or pH.

2.2.3 Porosity and advection of solids and porewater

The temporal and spatial evolution of porosity is calculated based on Eq. (1) and the constraint of volume conservation:

$$\begin{array}{}\text{(16)}& \sum _{\mathit{\theta }}^{n_{\mathrm{sld}}}{m}_{\mathit{\theta }}{V}_{\mathit{\theta }}=(\mathrm{1}-\mathit{\varphi }),\end{array}$$

where V_θ is the molar volume of solid species θ (m³ mol⁻¹) (Table 1). Summing Eq. (1) (multiplied by V_θ for all the considered solid species in the model) and using Eq. (16) yields

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& {\displaystyle \frac{\partial (\mathrm{1}-\mathit{\varphi })}{\partial t}}={\displaystyle \frac{\partial (\mathrm{1}-\mathit{\varphi })w}{\partial z}}+\sum _{\mathit{\theta }}^{n_{\mathrm{sld}}}{V}_{\mathit{\theta }}\left\{-R_{\mathit{\theta }}-\sum _{\mathit{\kappa }}^{n_{\mathrm{xrxn}}}{\mathit{\gamma }}_{\mathit{\kappa },\mathit{\theta }}{R}_{\mathit{\kappa }}\right.\\ & \left.+J_{\mathit{\theta }}+m_{\mathit{\theta }}\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 }}{m}_{\mathit{\theta }}\left(z^{\prime }\right)E_{\mathit{\theta }}(z^{\prime },z)\mathrm{d}{z}^{\prime }\right\}.\end{array}\end{array}$$

To satisfy the volume conservation (Eq. 16), porosity (ϕ) and advection rate (w) are calculated simultaneously by Eq. (17) assuming a scaling relationship. In the default setting, w is assumed to be constant and independent of ϕ, but two other user options are available in which wϕ or w(1−ϕ) can be assumed to be constant (cf. Wang et al., 2014; Chen et al., 2020).

Advection of porewater is assumed to be caused by steady-state porewater flow (e.g., Stonestrom et al., 1998). Thus, advection terms for aqueous species (the first terms on the right-hand side of Eqs. 2 and 3) are calculated based on a net water flux to the soil system, q (m yr⁻¹), and a water saturation profile (σ):

$$\begin{array}{}\text{(18)}& q=\mathit{\varphi }\mathit{\sigma }v.\end{array}$$

By default, SCEPTER fixes the profiles of q and σ, but these can be changed along model time as a user option. The depth profile of σ is assumed to increase linearly from the surface level to the water table, where σ=1, and is thus controlled by two variables – the value of σ at the surface and the depth of the water table (cf. Sect. 2.4; Kanzaki et al., 2020).

2.2.4 Diffusion

Molecular diffusion becomes slower in porous media because of tortuosity (Clennell, 1997). The tortuosities of pore air- and water-filled spaces are represented by factors τ_gas and τ_aq, both of which are parameterized with porosity and water content in soils according to Aachib et al. (2004):

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{gas}}={\mathit{\varphi }}^{\mathrm{1.4}}(\mathrm{1}-\mathit{\sigma }{)}^{\mathrm{2.4}},\text{(20)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{aq}}={\mathit{\varphi }}^{\mathrm{1.4}}{\mathit{\sigma }}^{\mathrm{2.4}}.\end{array}$$

See Table 5 for the molecular diffusion coefficients for gas and aqueous species where tortuosity factors are not accounted for.

2.2.5 Bio-mixing

SCEPTER parameterizes various styles of soil mixing (bio-mixing) within a transition matrix framework. Here, we provide a short description of the transition matrix framework for parameterization of bio-mixing. For more details, the reader is referred to, e.g., Boudreau (1997), Shull (2001) and Kanzaki et al. (2021).

The particle transport rate for solid species θ from layer i to layer j is defined here as P_θ,ij (yr⁻¹), given by

$$\begin{array}{}\text{(21)}& P_{\mathit{\theta },ij}={\displaystyle \frac{N_{\mathit{\theta },ij}}{\sum _{j=\mathrm{1}}^{n_{\mathrm{ml}}}{N}_{\mathit{\theta },ij}}}{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}},\end{array}$$

where N_θ,ij is the number of particles of species θ moved from layer i to layer j, n_ml is the total number of layers within the bio-mixed zone and τ is the time (yr) required for the particle displacements. Note that the particle transport probability, given as $P_{\mathit{\theta },ij}\times \mathit{\tau }$, corresponds to components at (i, j) of the transition matrix (Trauth, 1998; Shull, 2001). The time rate of change in soil concentration of solid species θ at layer i caused by bio-mixing ([$\partial m_{\mathit{\theta },i}/\partial t$]_mix) can then be described with P_θ,ij:

$$\begin{array}{}\text{(22)}& {\left[{\displaystyle \frac{\partial m_{\mathit{\theta },i}}{\partial t}}\right]}_{\mathrm{mix}}=-m_{\mathit{\theta },i}\sum _{j=\mathrm{1}}^{n_{\mathrm{ml}}}{P}_{\mathit{\theta },ij}+\sum _{j=\mathrm{1}}^{n_{\mathrm{ml}}}{\displaystyle \frac{\mathit{\delta }{z}_j}{\mathit{\delta }{z}_i}}{m}_{\mathit{\theta },j}{P}_{\mathit{\theta },ji},\end{array}$$

where m_θ,i is the concentration of θ (mol m⁻³) at soil layer i and δz_i is the thickness of layer i (m) (see Kanzaki et al., 2021, for a more detailed derivation). To simplify Eq. (22), a modified transition matrix for species θ is introduced, whose components at (i, j) are denoted as K_θ,ij and calculated based on the particle transport rate P_θ,ij:

$$\begin{array}{}\text{(23)}& K_{\mathit{\theta },ij}=\left\{\begin{array}{cc}\mathit{\delta }{z}_i{P}_{\mathit{\theta },ij}/\mathit{\delta }{z}_j& (i\ne j)\\ -\sum _{j\ne i}^{n_{\mathrm{ml}}}{P}_{\mathit{\theta },ij}& (i=j).\end{array}\right.\end{array}$$

From Eqs. (22) and (23), we obtain

$$\begin{array}{}\text{(24)}& {\left[{\displaystyle \frac{\partial m_{\mathit{\theta },i}}{\partial t}}\right]}_{\mathrm{mix}}=\sum _j^{n_{\mathrm{ml}}}{m}_{\mathit{\theta },j}{K}_{\mathit{\theta },ji}.\end{array}$$

Equation (24) can be regarded as the discretized form of the bio-mixing term shown in Eq. (1).

SCEPTER can implement a range of different transition matrices for parameterizing different soil mixing regimes, including Fickian mixing (e.g., bulk bioturbation), homogeneous mixing (e.g., deep rotary tilling prior to sowing row crops), inversion tilling, or particle-tracking automata-based simulators (see Boudreau et al., 2001; Choi et al., 2002; Kanzaki et al., 2019). Simulations shown here implement either Fickian or homogeneous mixing (see below). The transition matrix for Fickian bioturbation (parameterized with a biodiffusion coefficient D_b, Goldberg and Koide, 1962) can be expressed by

$$\begin{array}{}\text{(25)}& K_{\mathit{\theta },ij}=\left\{\begin{array}{ll}-K_{\mathit{\theta },ij}(j=i+\mathrm{1})& (i=j=\mathrm{1})\\ -K_{\mathit{\theta },ij}(j=i+\mathrm{1})-K_{\mathit{\theta },ij}(j=i-\mathrm{1})& (\mathrm{1}<i=j<n_{\mathrm{ml}})\\ -K_{\mathit{\theta },ij}(j=i-\mathrm{1})& (i=j=n_{\mathrm{ml}})\\ \mathit{\left\{}\right(\mathrm{1}-{\mathit{\varphi }}_i)D_{\mathrm{b},i}+(\mathrm{1}-{\mathit{\varphi }}_j\left)D_{\mathrm{b},j}\mathit{\right\}}/& (\mathrm{2}\le j=i+\mathrm{1}\le n_{\mathrm{ml}}\\ \mathit{\left\{}\mathit{\delta }{z}_i\right(\mathrm{1}-{\mathit{\varphi }}_i\left)\right(\mathit{\delta }{z}_i+\mathit{\delta }{z}_j\left)\mathit{\right\}}& \mathrm{or}\mathrm{1}\le j=i-\mathrm{1}\le n_{\mathrm{ml}}-\mathrm{1})\\ \mathrm{0}& \left(\mathrm{else}\right),\end{array}\right.\end{array}$$

where D_b,i represents the biodiffusion coefficient at soil layer i. SCEPTER adopts a depth-independent coefficient by default: $D_{\mathrm{b}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m² yr⁻¹ (cf. Jarvis et al., 2010; Astete et al., 2016).

The transition matrix for homogeneous mixing is given by

$$\begin{array}{}\text{(26)}& K_{\mathit{\theta },ij}=\left\{\begin{array}{ll}\mathit{\delta }{z}_i{P}_{\mathrm{h}}/\mathit{\delta }{z}_j& (i\ne j\mathrm{and}\mathrm{1}\le i,j\le n_{\mathrm{ml}})\\ -(n_{\mathrm{ml}}-\mathrm{1})P_{\mathrm{h}}& (\mathrm{1}\le i=j\le n_{\mathrm{ml}})\\ \mathrm{0}& \left(\mathrm{else}\right),\end{array}\right.\end{array}$$

where P_h (yr⁻¹) is the homogeneous transport rate of solid particles between soil layers. The model assumes P_h=0.01 yr⁻¹ as the default parameterization (cf. Kanzaki et al., 2021).

2.2.6 Surface area

Surface areas of solid species θ available for reaction with porewater (S_θ) significantly impact dissolution/precipitation rates of solid phases (Eq. 6), and thus the parameterization of surface area is of critical importance for predicting the behavior of elements in soils. However, there are significant differences in the parameterization of S_θ between models (e.g., Bolton et al., 2006; Steefel, 2009; Li et al., 2014) and/or between solid species (e.g., minerals vs. SOM; e.g., Jia et al., 2021). SCEPTER provides multiple options for parameterization of S_θ.

(1) Surface area parameterization based on hydraulic radius

The hydraulic radius of a given porous medium (r_H, m) can be defined as the ratio of pore volume against pore surface area (Fanchi, 2018). Conceptually, the value of r_H can be thought of as the average effective radius of particles within a weathering/soil system. In this formulation, the overall surface area of pores in a unit bulk soil/rock volume can be described as $\mathit{\varphi }/r_{\mathrm{H}}$ (m² m⁻³). The area available for mineral-phase θ can then be calculated as the product of the soil volume fraction of mineral θ and $\mathit{\varphi }/r_{\mathrm{H}}$, i.e.,

$$\begin{array}{}\text{(27)}& S_{\mathit{\theta }}=\mathit{\varphi }{m}_{\mathit{\theta }}{V}_{\mathit{\theta }}{r}_{\mathrm{H}}^{-\mathrm{1}}.\end{array}$$

Equation (27) is specified as the default option in SCEPTER, following Kanzaki and Kump (2017) and Kanzaki et al. (2020). A linear relationship between S_θ and m_θ as in Eq. (27) is also widely assumed in other models that calculate the surface area based on measured specific surface area (m² g⁻¹), soil concentration (mol m⁻³) and molar weight (g mol⁻¹) of minerals (see (4) below; Brantley and Lebedeva, 2011; Li et al., 2014).

However, several other reactive-transport models assume that the pore surface area available to mineral θ is not directly proportional to its volume fraction m_θV_θ but is instead proportional to ($m_{\mathit{\theta }}{V}_{\mathit{\theta }}{)}^{\mathrm{2}/\mathrm{3}}$ with an intension to convert volume fraction to surface area fraction (e.g., Steefel, 2009; see also Bolton et al., 2006, who implemented a similar dependence with a shrinking particle model). To realize this S_θ–m_θ relation, SCEPTER can also adopt an alternative function for S_θ:

$$\begin{array}{}\text{(28)}& S_{\mathit{\theta }}=\mathit{\varphi }(m_{\mathit{\theta }}{V}_{\mathit{\theta }}{)}^{\mathrm{2}/\mathrm{3}}{r}_{\mathrm{H}}^{-\mathrm{1}}.\end{array}$$

In Eq. (28), the hydraulic radius r_H is still used to specify the total surface area available to minerals.

With the options presented above (Eqs. 27 and 28), surface area normalized to porosity and mineral fraction (i.e., $\mathrm{1}/r_{\mathrm{H}}$) does not change with time or depth. To enable evolving surface area while using only average effective radius of particles, SCEPTER adopts two optional scaling relationships between surface area and porosity (e.g., Emmanuel and Berkowitz, 2005):

$$\begin{array}{}\text{(29)}& {\displaystyle }{r}_{\mathrm{H}}^{-\mathrm{1}}\propto (\mathrm{1}-\mathit{\varphi }{)}^{\mathrm{2}/\mathrm{3}},\text{(30)}& {\displaystyle }{r}_{\mathrm{H}}^{-\mathrm{1}}\propto {\mathit{\varphi }}^{\mathrm{2}/\mathrm{3}}.\end{array}$$

Equation (29) is adopted as the default option in the current version of SCEPTER.

(2) Parameterization based on particle size distribution

The surface area of a porous medium can be explicitly calculated if the shapes of component particles and the particle size distribution (PSD) are known, assuming that aggregation of particles does not reduce the surface area available to porewater (which can occur in natural systems depending on the assumed particle shapes). Tracking the PSD can be facilitated by adopting a shrinking particle model in which all particles maintain their shapes as their volumes change with ongoing dissolution/precipitation (e.g., Nicholson et al., 1990; Safari et al., 2009). For a batch solution system in which minerals only dissolve without being transported, surface area can be calculated relatively easily in a shrinking particle framework (e.g., LeBlanc and Fogler, 1987). For example, Beerling et al. (2020) adopted a shrinking particle model to calculate the particle size distribution and particle surface area in their “performance” simulations of basalt powder application onto croplands assuming a uniform spherical shape. The problem becomes more complex with solid-phase transport (e.g., in 1-D), where solid phases not only precipitate/dissolve but are also added at the top of the model domain and transported vertically (e.g., Eq. 1). To facilitate inter-model comparison and to add flexibility for representing reactive surface area, SCEPTER contains a numerical scheme and corresponding subroutine for explicit PSD tracking within a shrinking particle framework. The essential framework is provided below, while numerical implementation is detailed in Sect. 2.3.

Particles are all assumed to be spherical (e.g., Beerling et al., 2020) and chemically and mineralogically homogeneous regardless of the size. A shrinking particle scheme is then applied where net dissolution/precipitation occurs within soils. First, SCEPTER defines the PSD as f(r, t, z) – the number of particles per unit bulk soil volume for a given radius bin as a function of radius r, time t and soil depth z. Applying the general equation for solid species to f(r, t, z) yields

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& {\displaystyle \frac{\partial f(r,t,z)}{\partial t}}={\displaystyle \frac{\partial wf(r,t,z)}{\partial z}}-{\left[{\displaystyle \frac{\partial f(r,t,z)}{\partial t}}\right]}_{\mathrm{diss}}\\ & +{\left[{\displaystyle \frac{\partial f(r,t,z)}{\partial t}}\right]}_{\mathrm{rain}}-f(r,t,z)\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 }}f(r,t,z^{\prime })E_{\mathit{\theta }}(z^{\prime },z)\mathrm{d}{z}^{\prime },\end{array}\end{array}$$

where [∂f(r, t, z)/∂t]_diss is the time rate of change in f(r, t, z) caused by net dissolution and [∂f(r, t, z)/∂t]_rain is the supply rate of particles at the soil surface with a specified PSD, which satisfies the following:

$$\begin{array}{}\text{(32)}& {\displaystyle }\sum _{\mathit{\theta }}^{n_{\mathrm{sld}}}{V}_{\mathit{\theta }}(R_{\mathit{\theta }}+\sum _{\mathit{\kappa }}^{n_{\mathrm{xrxn}}}{\mathit{\gamma }}_{\mathit{\kappa },\mathit{\theta }}{R}_{\mathit{\kappa }})=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{4}\mathit{\pi }{r}^{\mathrm{3}}}{\mathrm{3}}}{\left[{\displaystyle \frac{\partial f(r,t,z)}{\partial t}}\right]}_{\mathrm{diss}}\mathrm{d}r,\text{(33)}& {\displaystyle }\sum _{\mathit{\theta }}^{n_{\mathrm{sld}}}{V}_{\mathit{\theta }}{J}_{\mathit{\theta }}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{4}\mathit{\pi }{r}^{\mathrm{3}}}{\mathrm{3}}}{\left[{\displaystyle \frac{\partial f(r,t,z)}{\partial t}}\right]}_{\mathrm{rain}}\mathrm{d}r.\end{array}$$

Shifts in the particle size distribution caused by net dissolution/precipitation ([∂f(r, t, z)/∂t]_diss) can be formulated with a population balance equation (e.g., LeBlanc and Fogler, 1987; Iggland and Mazzotti, 2012):

$$\begin{array}{}\text{(34)}& {\left[{\displaystyle \frac{\partial f(r,t,z)}{\partial t}}\right]}_{\mathrm{diss}}=-{\displaystyle \frac{\partial g(r,t,z)f(r,t,z)}{\partial r}},\end{array}$$

where g(r, t, z) is the particle growth rate (m yr⁻¹). Within a shrinking particle framework, the particle growth rate is specified as (e.g., LeBlanc and Fogler, 1988; Safari et al., 2009)

$$\begin{array}{}\text{(35)}& g(r,t,z)={\displaystyle \frac{\partial r}{\partial t}}=-k^{\prime }.\end{array}$$

Here, k^′ is the particle dissolution rate in units of m yr⁻¹. Loss of particles is allowed at the lower boundary, i.e., particles are assumed to be completely dissolved when r becomes smaller than the minimum radius considered in the model. Also, when dissolution (imposed by Eq. 32) is intense enough to consume all existing particles, particles are allowed to be lost. Meanwhile, particles are not allowed to grow over the maximum radius set by the model. At each time step, Eqs. (32), (34) and (35) are solved at once to obtain the values for [∂f(r, t, z)/∂t]_diss that satisfy the volume balance (Sect.  2.3).

Given a PSD for particles applied at the soil surface, defined here as f_rain(r) that is only a function of radius but not of time or depth, one can obtain a function k_rain(r, t, z) in units of yr⁻¹ satisfying Eq. (33) and the following:

$$\begin{array}{}\text{(36)}& {\left[{\displaystyle \frac{\partial f(r,t,z)}{\partial t}}\right]}_{\mathrm{rain}}=k_{\mathrm{rain}}(r,t,z)f_{\mathrm{rain}}\left(r\right).\end{array}$$

By solving Eqs. (31)–(36), one can obtain the full PSD reflecting both shrinking particles via net precipitation/dissolution and addition/transport of particles. The above scheme generally satisfies the volume balance:

$$\begin{array}{}\text{(37)}& \mathrm{4}\mathit{\pi }/\mathrm{3}{\int }_{r=r_{min}}^{r=r_{max}}{r}^{\mathrm{3}}f\left(r,t,z\right)\mathrm{d}r=\sum _{\mathit{\theta }}^{n_{\mathrm{sld}}}{m}_{\mathit{\theta }}{V}_{\mathit{\theta }}.\end{array}$$

The total surface area is then calculated based on the PSD as follows:

$$\begin{array}{}\text{(38)}& \mathrm{4}\mathit{\pi }{\int }_{r=r_{min}}^{r=r_{max}}{r}^{\mathrm{2}}f\left(r,t,z\right)\mathrm{d}r=\mathit{\varphi }{r}_{\mathrm{H}}^{-\mathrm{1}}.\end{array}$$

Here, the total pore surface area per unit soil/rock volume is again represented by the hydraulic radius r_H to facilitate comparison with the default surface area calculation where r_H is directly specified as the average effective radius of particles (see (1) above). The surface areas for individual solid species (S_θ) can then be calculated by either Eq. (27) (the default in SCEPTER) or (28). When adopting the PSD-based formulation, the surface area evolves temporally and spatially in a dynamic fashion without any further parameterization.

(3) Surface roughness

Mineral surfaces are not necessarily smooth and can have complicated geometry and resultant surface roughness that increases the surface area per unit volume/mass (e.g., White and Peterson, 1990). One can introduce a roughness factor λ to correct the surface area calculated for a smooth geometry of solid particles with an assumed shape. This additional parameterization is especially useful when one considers a particle size distribution based on ideal shapes ((2) above). The roughness factor is calculated assuming a dependency on particle size following Beerling et al. (2020) (cf. Navarre-Sitchler and Brantley, 2007):

$$\begin{array}{}\text{(39)}& \mathit{\lambda }=({\mathrm{10}}^{\mathrm{10}}r{)}^{\mathrm{0.33}}.\end{array}$$

As a user option, one can include λ in the PSD calculation method in (2); in this case, the k^′ and r² terms in Eqs. (35) and (38) are replaced with λk^′ and λr², respectively.

Implementation of a roughness factor is of limited use for formulations that calculate surface area directly from hydraulic radius (Eqs. 27–30), because the hydraulic radius should in principle already account for roughness in the pore surface. As a result, in the default surface area calculation (Eqs. 27 and 29), a roughness factor is not included, although SCEPTER contains a user option to add a roughness term to Eqs. (27)–(30).

(4) Specific surface area of solid species

The options for calculating surface area in SCEPTER described above are based on the total surface area of pores and the fractions of individual solid species in soils. These approaches all assume that every patch of pore surface is mineralogically and chemically homogeneous. However, individual solid species can have unique particle size distributions and thus specific surface areas that are different from the surface areas of the bulk soil multiplied by solid species fractions. Based on the concept of specific surface areas for individual solid species, SCEPTER includes an additional user option that enables the surface area calculation for individual solid species.

Defining the specific surface area of solid species θ as A_θ (m² g⁻¹), the surface area of θ available to porewater (S_θ, m² m⁻³) can alternatively be represented by (e.g., Brantley and Lebedeva, 2011; Li et al., 2014)

$$\begin{array}{}\text{(40)}& S_{\mathit{\theta }}=A_{\mathit{\theta }}{M}_{\mathit{\theta }}{m}_{\mathit{\theta }},\end{array}$$

where M_θ (g mol⁻¹) is the molar weight of solid species θ (Table 1). To facilitate comparison with the surface area parameterization using a hydraulic radius (Eq. 27), we introduce an apparent hydraulic radius for solid species θ as r_H,θ (m):

$$\begin{array}{}\text{(41)}& r_{\mathrm{H},\mathit{\theta }}^{-\mathrm{1}}={\displaystyle \frac{A_{\mathit{\theta }}{\mathit{\rho }}_{\mathit{\theta }}}{\mathit{\varphi }}},\end{array}$$

where ρ_θ (g m⁻³) is the particle density of solid species θ and ${\mathit{\rho }}_{\mathit{\theta }}\equiv M_{\mathit{\theta }}/V_{\mathit{\theta }}$. The apparent hydraulic radius r_H,θ (m) can be thought of as the average effective radius of particles which are composed solely of solid species θ. SCEPTER has the option of calculating individual surface areas according to Eqs. (40) and (41), with a specified r_H,θ the evolution of which can be further specified by replacing r_H with r_H,θ in Eq. (29) or (30).

SCEPTER can evaluate specific surface areas for individual solid phases by tracking PSDs for individual species. The scheme introduced in (2) still holds, but in this case the PSD is defined and calculated for individual solid species. Defining the PSD and particle growth and dissolution rates for solid species θ as f_θ(r, t, z), g_θ(r, t, z) and $k_{\mathit{\theta }}^{\prime }$, respectively, the equations to solve PSD for bulk soil, i.e., Eqs. (31)–(38), can be used to solve f_θ(r, t, z) by replacing f(r, t, z), g(r, t, z) and k^′ with f_θ(r, t, z), g_θ(r, t, z) and $k_{\mathit{\theta }}^{\prime }$, respectively, and dropping summation symbols and notations (i.e., replacing ${\sum }_{\mathit{\theta }}^{n_{\mathrm{sld}}}{V}_{\mathit{\theta }}(R_{\mathit{\theta }}+{\sum }_{\mathit{\kappa }}^{n_{\mathrm{xrxn}}}{\mathit{\gamma }}_{\mathit{\kappa },\mathit{\theta }}{R}_{\mathit{\kappa }}$), ${\sum }_{\mathit{\theta }}^{n_{\mathrm{sld}}}{V}_{\mathit{\theta }}{J}_{\mathit{\theta }}$ and ${\sum }_{\mathit{\theta }}^{n_{\mathrm{sld}}}{m}_{\mathit{\theta }}{V}_{\mathit{\theta }}$ with $V_{\mathit{\theta }}(R_{\mathit{\theta }}+{\sum }_{\mathit{\kappa }}^{n_{\mathrm{xrxn}}}{\mathit{\gamma }}_{\mathit{\kappa },\mathit{\theta }}{R}_{\mathit{\kappa }}$), V_θJ_θ and m_θV_θ, respectively). Then, the specific surface area for solid species θ can be calculated as

$$\begin{array}{}\text{(42)}& A_{\mathit{\theta }}={\displaystyle \frac{\mathrm{3}{\int }_{r=r_{min}}^{r=r_{max}}{r}^{\mathrm{2}}{f}_{\mathit{\theta }}\left(r,t,z\right)\mathrm{d}r}{{\mathit{\rho }}_{\mathit{\theta }}{\int }_{r=r_{min}}^{r=r_{max}}{r}^{\mathrm{3}}{f}_{\mathit{\theta }}\left(r,t,z\right)\mathrm{d}r}},\end{array}$$

and the corresponding apparent hydraulic radius can be evaluated as

$$\begin{array}{}\text{(43)}& r_{\mathrm{H},\mathit{\theta }}^{-\mathrm{1}}={\displaystyle \frac{\mathrm{3}{\int }_{r=r_{min}}^{r=r_{max}}{r}^{\mathrm{2}}{f}_{\mathit{\theta }}\left(r,t,z\right)\mathrm{d}r}{\mathit{\varphi }{\int }_{r=r_{min}}^{r=r_{max}}{r}^{\mathrm{3}}{f}_{\mathit{\theta }}\left(r,t,z\right)\mathrm{d}r}}.\end{array}$$

The roughness factor λ can further be included by replacing $k_{\mathit{\theta }}^{\prime }$ and r² with $\mathit{\lambda }{k}_{\mathit{\theta }}^{\prime }$ and λr², respectively, as described for the PSD calculation for bulk soil ((3) above). For instance, ground minerals have been characterized by significant surface roughness based on measured surface areas that are larger than expected from their grain sizes (e.g., Brantley and Mellot, 2000; Renforth, 2012; Renforth et al., 2015).

(5) Reactions not limited by surface area

Decomposition of some solid species can proceed without being affected by the surface area available to porewater. In the current version of SCEPTER, three classes of SOM are assumed to decay depending on their concentrations but independently of their surface areas (e.g., Jia et al., 2021):

$$\begin{array}{}\text{(44)}& S_{\mathit{\theta }}=m_{\mathit{\theta }}.\end{array}$$

For all the other solid species, their dissolution/precipitation kinetics are assumed to be dependent on the surface areas which are calculated based on r_H or r_H,θ, either directly specified by the user or calculated from tracked PSD(s), as described above.

Table 7 Boundary conditions for example simulations.

^a Only IDs of solid species and extra reactions are denoted (Tables 1 and 6). ^b Composition of dust is assumed to be the basalt composition by Beerling et al. (2020). ^cN is the number of grid points, z_tot is the bottom depth of simulated soil, z_sat is the depth of the water table, σ₀ is the surface water saturation ratio and ϕ₀ is the initial porosity. ^d When surface area is calculated based on tracked particle size distributions, “PSD” is denoted. ^e (1) Rainwater composition is assumed to be that of pure water saturated with respect to atmospheric CO₂ and O₂. (2) Atmospheric CO₂ and O₂ are assumed to be $\mathrm{3.16}\times {\mathrm{10}}^{-\mathrm{4}}$ and 0.21 atm, respectively. (3) Parent rock concentrations of albite and kaolinite are 10 and 10⁻³ wt %, respectively. (4) Parent rock concentrations of pyrite and goethite are 0.56 and 10⁻³ wt %, respectively. (5) Parent rock concentrations of solid species are from those for Site 1 in Table 8. (6) Abiotic weathering, biotic weathering and basalt application experiments in Figs. 3, 4 and 6–9 are repeated with surface area parameterization based on tracked particle size distributions for bulk soil (Fig. 10). The same basalt application experiment was conducted except with tracking PSDs for individual solid species and calculating corresponding specific surface areas (Figs. 11 and 12).

Figure 1 Simulation of abiotic weathering of albite. Soil concentrations, saturation states and dissolution/precipitation rates of solid species are plotted in panels (a)–(c), respectively; porewater concentrations of aqueous species and pH in (d) and (e), respectively; soil porosity and particle density in (f) and (g), respectively; ratio of total volume of solid species against solid space specified with porosity in (h). See Table 7 for the details on the boundary conditions of the simulation.

2.3 Numerical implementation

Initialization of SCEPTER involves loading input data, including chosen species to be simulated, and initial and boundary conditions (see Sect. 2.4). In the case of a restart experiment from a previous simulation, initial conditions are overwritten with the results from the previous (restart) experiment. After initialization, SCEPTER begins time integration of the governing equations (see Sect. 2.2). For a given time step, the boundary conditions and time step duration can be modified if the user selects time-varying changes to water flux, temperature, particle rain rate or water saturation ratio. Time step duration also evolves adaptively depending on the time to convergence in the previous integration step (from 10⁻¹⁸ up to 10³ years). Kinetic and thermodynamic constants are then updated, and the concentrations of all chosen species are solved via Newton–Raphson iteration in a fully coupled way (see below). Finally, porosity, surface area and advection rate(s) of solid phases are updated. By default, SCEPTER updates porosity, surface area and advection rates iteratively at each time step as verified by porosity convergence, but there is also a user option to bypass the porosity iteration and simply update porosity, surface area and advection rate(s) once per time step. The latter option is computationally less expensive and will yield the same solutions when the time step duration is relatively small. The default criterion for numerical convergence is a solution difference of $<{\mathrm{10}}^{-\mathrm{9}}$ (mol m⁻³ or m³ m⁻³) from the previous iteration in each time-integration step (cf. Steefel, 2009).

Figure 2 Simulation of abiotic weathering of pyrite. Soil concentrations and saturation states of solid species are plotted in (a) and (b), respectively; rate profile of reactions within soil in (c); porewater concentrations of aqueous species and pH in (d) and (e), respectively; concentrations of soil gases in (f); soil porosity and particle density in (g) and (h), respectively; ratio of the total volume of solid species against solid space specified with porosity in (i). See Table 7 for details on the boundary conditions of the simulation.

The governing equations of SCEPTER are differentiated via a finite-difference method. First-order upwind and second-order central differencing schemes are applied to advection and diffusion terms, respectively (Eqs. 1–3). Equation (24) is used as a difference form of the bio-mixing term in Eq. (1). Time derivatives (left-hand side of Eqs. 1–3) are discretized in accordance with a backward Euler scheme. The finite-difference expressions are then solved via a fully coupled Newton–Raphson method (Steefel and Lasaga, 1994). Parent rock values are enforced below the bottom of the model domain as a boundary condition for solid species. For aqueous and gas species, the bottom of the model domain is assumed to be impermeable with respect to molecular diffusion (i.e., zero concentration gradient), while the compositions of rainwater and the overlying atmosphere are enforced as boundary conditions at the top of the model domain on dissolved and gaseous species, respectively. When there is no user input, a value of 10⁻²⁰ mol m⁻³ or mol L⁻¹ is assumed for the boundary parent-rock or rainwater concentration, respectively, and 0.21, 10^−3.5, 10⁻⁹ and $\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{7}}$ atm for atmospheric O₂, CO₂, NH₃ and N₂O, respectively. Concentrations of all dependent aqueous species as well as associated rate constants are always updated (Sect. 2.2.2, Tables 1–4), including during Newton–Raphson iteration.

The calculations of porosity and advection rate of solid species are conducted by differentiating Eq. (17), again using an implicit finite-difference method. Because the reaction term (the summation on the right-hand side of Eq. 17) is fixed by the solution of Eqs. (1)–(3), the discretized equations become linear with respect to porosity and advection rate, and thus the Newton–Raphson method is not used to solve the finite-difference form of Eq. (17). Boundary conditions are imposed consistently with those for solid species – e.g., the porosity and uplift rate of the parent rock at the bottom of the model domain.

Table 8 Boundary conditions for natural weathering sites.

^a Solid species denoted as amorphous, hornblende, pyroxene, 10Å clay, 14Å clay, plagioclase and zeolite phases in the original data are assumed to be amorphous Si, tremolite, diopside, illite, Ca-beidellite, albite and analcime, respectively. Solid species are denoted with ID in Table 1. ^b w is the uplift rate in units of 10⁻⁵ m yr⁻¹ (Larsen et al., 2014), q is the net water flux to sites in units of m yr⁻¹ (Fekete et al., 2002), σ₀ is the surface water saturation ratio, T_C is temperature in ^∘C (Fick and Hijimans, 2017) and OM rain value is provided with units of g C m⁻² yr⁻¹. ^c Obtained as 1.5 times the net primary production following Beerling et al. (2020).

Table 9 Observed soil depth profiles of OM and pH^*.

^* From Hengl et al. (2017).

Figure 3 Soil concentrations of solid species in simulation of abiotic weathering at Site 1 in Table 8. See Table 7 for the details on the boundary conditions of the simulation.

In the default version of SCEPTER in which PSDs are not tracked, surface area is calculated by combining Eqs. (27)–(30), (40) and (41). When PSD tracking is enabled in the surface area calculation, the governing equation for PSD (Eq. 31) is discretized by a finite-difference method and solved time-implicitly as for the solution of the governing equation for solid species (Eq. 1). The PSD shifts caused by net volume change by dissolution/precipitation (the second term on the right-hand side of Eq. 31) are enforced from solutions of Eqs. (32), (34) and (35) by the first-order upwind finite-difference scheme and the Newton–Raphson method. A PSD for parent rock is imposed as the lower boundary condition. This procedure is repeated for individual solid species when calculating species-specific PSDs and surface areas ((4) in Sect. 2.2.6). By default, SCEPTER considers the PSD calculation converged when the difference in particle number from the previous iteration is less than 10⁻¹² times the maximum number of particles. The calculated PSD is truncated below one particle per cubic meter of bulk soil for a given radius bin.

Figure 4 Simulation of abiotic weathering at Site 1 in Table 8. Porewater concentrations of aqueous species and pH are plotted in (a)–(d) and (e), respectively; concentrations of soil gases in (f); soil porosity and particle density in (f) and (g), respectively; ratio of the total volume of solid species against solid space specified with porosity in (h). See Table 7 for details on the boundary conditions of the simulation.

2.4 User input

2.4.1 Independent variables

Solid, aqueous and gas species to be tracked in a simulation are listed in individual input files (slds.in, solutes.in and gases.in). If a user wants to include reactions other than dissolution/precipitation reactions specified for individual solid species, the ID string of the extra reaction (e.g., fe2o2 in Table 6 for aqueous Fe(II) oxidation by O₂) can be specified in another input file (exrxns.in).

Figure 5 Locations of natural weathering sites whose boundary conditions are tabulated in Table 8 and utilized for example simulations in Sect. 3.

Figure 6 Simulation of biotic weathering at Site 1 in Table 8. Shown are concentrations of simulated solid species. Note that the model configuration is the same as that for Fig. 3 except that soil respiration and bio-mixing are included.

Figure 7 Simulation of biotic weathering at Site 1 in Table 8. Panels show the depth profiles of aqueous and gas species and soil physical properties as in Fig. 4 except that simulation includes soil respiration and bio-mixing.

2.4.2 Boundary conditions

Fundamental variables such as grid size, total depth of simulated soil, water flux, water table depth, surface water saturation, mixing layer thickness, application rate of powdered rock feedstock, OM rain flux, temperature, initial/bottom uplift rate and initial soil porosity are specified in the input file frame.in. Modifying options for, e.g., mixing regime and method of surface area calculation can be made by modifying another input file, switches.in. One can also specify whether to do a re-start experiment in switches.in, and, when doing a re-start experiment, the previous experiment from which the current simulation should be restarted can be specified in frame.in.

Surface area of parent rock is calculated from the average effective radius of particles specified in frame.in. When one chooses to calculate specific surface areas for individual solid species, different average radius values can be assigned to different solid species in the input file sa.in. The particle size distributions for parent rock are calculated assuming log-normal distributions with 1 SD centered at the radius value(s) input from frame.in or sa.in. Currently, the PSD for OM rain is assumed to be the same as that for parent rock, while the PSD for rock feedstock needs to be specified within the source codes by providing mean radius(es) and standard deviation(s).

Figure 8 Basalt powder application at Site 1 in Table 8. Panels show the depth profiles of simulated solid species as in Fig. 6, except that basalt powder is continuously added at a rate of 40 t ha⁻¹ yr⁻¹ and plots are focused on the soil mixed layer (0–25 cm).

Boundary conditions for solid, aqueous and gas species need to be provided as parent rock, rainwater and atmospheric concentrations, respectively, in the corresponding input files (parentrock.in, rain.in and atm.in). Compositions of solids being introduced at the soil surface must also be specified using individual input files (dust.in and OM_rain.in, respectively). In the case of time-varying changes to water flux, temperature and water saturation ratio, a series of additional input files is necessary (T_temp.in, q_temp.in and Wet_temp.in). See the full README information included in the code repository for further details (code availability).

Figure 9 Basalt powder application at Site 1 in Table 8. Panels show the depth profiles of aqueous and gas species and soil physical properties as in Fig. 7 except that basalt powder is continuously added at a rate of 40 t ha⁻¹ yr⁻¹.

2.4.3 Initial conditions

At the beginning of a simulation, concentrations of solid, aqueous and gas species are assumed to be the same as those of parent rock, rainwater and atmosphere, respectively, which can be specified by the user with the corresponding input files (see above). Initial surface areas as well as PSDs at all depths are assumed to be the same as those of parent rock specified by the user. Porosity similarly takes the initial value provided by the user at all depths. When the volume sum of all chosen solid species in the parent rock is less than the solid fraction implied by the initial/bottom porosity, then an imaginary “bulk” species is additionally tracked to occupy the void space whose physical properties are assumed to be the same as those of kaolinite but with no reactions allowed (R_θ≡0 with θ= “bulk”). When the volume sum of all chosen solid species exceeds the assigned value from the initial/bottom porosity, the parent rock concentrations of chosen solid species are rescaled to be consistent with the initial/bottom porosity. Volume conservation (Eq. 16) is thus always satisfied by SCEPTER (cf. Archer et al., 2002; Munhoven, 2021; Kanzaki et al., 2021).

Figure 10 Simulations of abiotic weathering (a–c), biotic weathering (d–f) and basalt powder application (g–i) at Site 1 in Table 8 with surface area calculated based on tracked particle size distribution (PSD) for bulk soil. PSDs at surface and bottom of soil mixed layer are shown in (a), (d) and (g). Pore surface area per unit pore volume is plotted against depth (b, e, h) or porosity (c, f, i) and compared with the calculation with the default surface area parameterization. See Sect. 3.1.4 and Table 7 for the details on the boundary conditions of the simulations.

3 Application examples

To illustrate the features and capabilities of SCEPTER, we present a series of idealized example experiments. First, the basic features of SCEPTER are illustrated by simulating abiotic weathering (e.g., without SOM and mixing; Sect. 3.1.1), biotic weathering (with SOM and mixing; Sect. 3.1.2), and basalt application (with SOM, mixing, and additional mineral supplied at the soil surface; Sect. 3.1.3) scenarios with the default surface area calculation method. The same set of simulations is then repeated using the PSD-based surface area calculation (Sect. 3.1.4). Finally, we show a series of SCEPTER runs driven by observational data from a subset of relatively pristine natural systems coupled with the USGS soil chemistry database (U.S. Geological Survey, 2016), and track time-dependent CO₂ capture across a range of timescales (Sect. 3.2).

Figure 11 Particle size distributions (PSDs) at surface and bottom of soil mixed layer for individual solid species in simulation of basalt powder application at Site 1 in Table 8 with the surface area calculated based on tracked PSDs for individual solid species. See Sect. 3.1.4 and Table 7 for the details on the boundary conditions of the simulation.

3.1 Illustration of basic features of SCEPTER

3.1.1 Abiotic weathering

The simplest configuration of SCEPTER explored here is the simulation of abiotic weathering of albite and pyrite, respectively (Figs. 1, 2; Table 7). Dominant controls on and locations of reaction fronts of albite and pyrite are consistent with simpler models (e.g., water flow for albite weathering, and water table depth on pyrite weathering; e.g., Kanzaki et al., 2020). SCEPTER allows tracking of time-dependent changes to gradients in solid and solute species, which in the simple abiotic cases evolve as expected – in the first case, with gradual conversion of albite to kaolinite and progressive release of Si and Na to soil pore fluids (Fig. 1), and in the second case a sharp reaction front along which pyrite is converted to goethite and O₂ is drown down to negligible levels at the water table depth (Fig. 2). Overall depth-dependent changes to porosity and particle density are relatively small in both of these simplified abiotic cases (Figs. 1 and 2). Note that, in these examples input parent-rock concentrations of minerals (10 wt % albite and 0.56 wt % pyrite) are smaller than inferred from assumed initial/bottom porosity (0.1). Thus, an imaginary bulk species is simulated along with the above minerals (not shown in Figs. 1 and 2), ensuring that volume conservation is satisfied (Eq. 16).

Figure 12 Surface area of individual solid species per unit pore volume plotted against depth (focused on mixed layer; 0–25 cm) in simulation of basalt powder application at Site 1 in Table 8 with the surface area calculated based on tracked PSDs for individual solid species (“Full PSD”). The pore surface areas per unit pore volume calculated with the default parameterization (“Def.”) and based on PSD for bulk soil (“Bulk PSD”) are also plotted for comparison. See Sect. 3.1.4 and Table 7 for details on the boundary conditions of the simulation.

Results of a more complicated abiotic weathering simulation are shown in Figs. 3 and 4, in which SCEPTER is fed by the bulk mineralogy and climatological boundary conditions for a natural weathering site (Site 1 in Fig. 5 and Table 8) (see also Sect. 3.2). Here, we assume zero OM rain to the system and no mixing to exclude biotic aspects of weathering (cf. Sect. 3.1.2) and run an abiotic weathering simulation to reach a steady state. Because production of soil CO₂ is minimized without a flux of OM to the soil surface, soil CO₂ drops at depth as a result of cation release (and alkalinity production) from mineral dissolution. Primary silicates such as albite and diopside are largely dissolved. Initially, carbonate phases dissolve at the surface but precipitate at depth. However, as the system approaches steady state (∼10⁴ years), carbonate phases redissolve and secondary clays accumulate (Fig. 3). Solute profiles evolve in accordance with mineral profiles, and once again overall changes to soil porosity and particle density are relatively small (Fig. 4). In this case, SCEPTER does not include the imaginary bulk species (Fig. 3; Tables 7 and 8), but volume conservation is nevertheless always satisfied.

3.1.2 Biotic weathering: addition of organic matter flux to soil surface and soil mixing

Following the expansion of vascular plant ecosystems across Earth's land surface, natural weathering generally occurs in the presence of SOM and soil respiration (cf. Volk, 1987; Berner, 1992; Kanzaki and Kump, 2017; Ibarra et al., 2019; Wen et al., 2021). Indeed, the recycling of SOM represents a critical component of soil acid–base chemistry and CO₂ cycling and is likely to change significantly across a range of ecosystems in the coming century (e.g., Brovkin et al., 2013; McGuire et al., 2018). The current version of SCEPTER can simulate up to three classes of SOM of varying lability. By default, the lability of each class of SOM progressively decreases, with turnover times of 1, 8 and 10³ years, respectively. However, the intrinsic SOM labilities can be scaled arbitrarily for a range of applications (e.g., manure application, biochar amendment, etc.), provided that reliable kinetic data are implemented. The OM flux to the soil surface can either be specified arbitrarily or can be scaled to above-ground net primary production; by default, SCEPTER scales the OM flux to the soil surface as 1.5 times above ground primary production, following Beerling et al. (2020). Note that additional biotic factors relating to the soil OM cycling, such as secretion of organic acids and uptake of nutrients and/or cations by plants, are not treated in the current version of the model but will be included in a future release.

Figure 13 Model fit to observed depth profiles of soil OM and pH for Sites 2 and 3 in Table 8. See Sect. 3.2.1 for the details.

Figures 6 and 7 show results of a simple biotic weathering experiment in which Class 2 SOM (here taken to represent a generalized “natural” SOM; e.g., Chen et al., 2010) and soil bioturbation (here termed Fickian mixing) are added to the experimental conditions for the abiotic weathering experiment shown in Figs. 3 and 4. As a result of SOM respiration, soil CO₂ builds up as high as >0.01 atm as weathering proceeds (Fig. 7). This buildup in CO₂ significantly lowers soil pH, which promotes dissolution of primary mineral phases and leads to an increase in porosity around the mid-depth of the soil mixed layer (∼0.1 m). A range of silicate mineral phases becomes unstable and dissolves with the presence of SOM respiration (e.g., compare Figs. 3 and 6). On a timescale of ∼10² years, carbonate phases dissolve at the surface of the model domain and precipitate at the bottom of the soil column, but all carbonate phases ultimately dissolve at steady state.

At steady state, soil CO₂ is around 0.006 atm, still more than 10 times higher than the ambient atmospheric level. SOM dominates amongst the solid phases near the soil surface, below which clays and Si-oxide phases become more dominant (Fig. 6). Dissolved cations are concentrated only at the deepest depths of the model domain, consistent with the observed distributions of mineral phases (Fig. 7). Dissolved Al shows higher concentrations in the biotic weathering scenario relative to the abiotic weathering regime due to increased solubility at lower pH (e.g., compare Figs. 4 and 7). Average grains (aggregates) of soil become less dense closer to the surface because they consist more of SOM whose particle density is relatively small (1.5 g cm⁻³) compared with those of other solid phases (Fig. 7). Note that volume conservation is again satisfied throughout the experiment even with a continuous flux of OM to the soil surface (Fig. 7i).

3.1.3 Basalt powder application: enhanced rock weathering (ERW)

A useful user option for simulating enhanced rock weathering in SCEPTER is the ability to conduct a re-start experiment from the previously conducted simulation. This feature allows the user to first run a spinup experiment – for instance, to fit the model to current observations at a given location – and then examine the impacts of adding crushed rock feedstock in a transient experiment branched from the spinup simulation. Here, we illustrate this procedure by conducting a re-start experiment including basalt powder application to the soil surface branched from the biotic weathering experiment for Site 1 (Table 8 and Fig. 5) presented in Sect. 3.1.2, which has been spun up to steady state by running the model for ∼10⁵ years.

Figure 14 CO₂ capture predicted in simulations of continuous basalt powder application at Sites 2 and 3 in Table 8 with the surface area calculation based on the assumed hydraulic radius (“Def.”), tracked PSD for bulk soil (“Bulk PSD”) and tracked PSDs for individual solid species (“Full PSD”). See Sect. 3.2.2 for details.

Figure 15 CO₂ capture predicted in simulations of pulsative basalt application at Sites 2 and 3 in Table 8. Milled basalt (dominated by 5, 20, 50 and 70 µm particles) is applied non-continuously and homogeneously mixed at the soil surface. See Sect. 3.2.2 for details.

We set the application rate of basalt powder at 40 t ha⁻¹ yr⁻¹, whose mineral composition is assumed to be the same as that adopted for basalt by Beerling et al. (2020). Because basalt powder will generally contain significant amounts of reduced Fe, we implement aqueous Fe(II) oxidation by O₂ as an extra reaction in the re-start experiment (Table 7). Application of basalt powder adds a large amount of primary minerals to surface soils, where they have generally been depleted during spinup to steady state (Fig. 8). One outcome of this is that carbonate phases precipitate during basalt application (Fig. 8d) together with kaolinite and goethite and later with smectite (Fig. 8c, n and o). There is a drop in SOM near the soil surface (Fig. 8p), despite the OM rain flux remaining identical to that in the spinup simulation (Sect. 3.1.2). Rained minerals with high solubility and reactivity result in similar depth profiles to that of SOM (e.g., Fig. 8h, k, l and p), while those with lower solubility/reactivity tend to accumulate at depth (e.g., Fig. 8f and g). Soil CO₂ accumulates as high as >0.075 atm because surface porosity drops significantly following basalt application (Fig. 9g). Soil O₂ becomes depleted over time, in part due to the porosity change, but also because of the introduction of Fe(II) oxidation as an additional O₂ sink that does not exist during the spinup simulation. Once again, volume conservation is satisfied even with a significant external flux of mineral phases to the system (Fig. 9i).

3.1.4 The impact of particle size tracking

In order to illustrate the potential use of optional PSD tracking for calculating the mineral surface area in SCEPTER, we repeated the same experiments presented above for Site 1 (Fig. 5 and Table 8; Sects. 3.1.1–3.1.3) while implementing mineral surface area based on tracked PSDs for bulk soil minerals, including a dynamic roughness factor (see (2)–(4) in Sect. 2.2.6). For parent rock, we assume a log-normal size distribution centered at 10 µm with a standard deviation of 1 log unit (i.e., f_PR(r)∝ exp(−(log r+ 5)²/2), where f_PR is defined as the PSD of parent rock and r is the particle radius in meters). We assign the same distribution to the OM rain flux in the biotic weathering experiment, while basalt powder is assumed to have a PSD that mixes equally four log-normal distributions centered at 5, 20, 50 and 70 µm with standard deviations of 0.2 log units. Note that because our intention is to illustrate the basic behavior of the model, the PSD parameterizations adopted here are not necessarily realistic for a given application (cf. Eberl et al., 1998; Sklar et al., 2017; Beerling et al., 2020) and should be modified within the code by the user when necessary.

Concentration profiles for solid, aqueous, and gas species are largely similar to those obtained with the default surface area scheme, particularly for the abiotic and biotic natural weathering experiments, because the resultant surface areas are fairly similar (Fig. 10b, c, e and f). In the abiotic weathering experiment, particles with small radii dissolve and disappear rapidly, particularly at shallower depths (Fig. 10a). In the biotic weathering experiment, the soil surface receives an OM rain flux whose particle distribution is assumed to be the same as that of parent rock, and thus the relative depletion of small particles in the PSD is most prominent near the bottom of the mixed layer (dotted curves in Fig. 10d), where mineral dissolution is significant but the particle input from the OM rain flux is attenuated.

In contrast, mineral surface area in the basalt application case is significantly different between the default and PSD-tracking schemes (Fig. 10h and i). This is because the crushed rock particles, whose size distributions are centered at 5–70 µm with relatively small standard deviations, become the dominant constituent of the surface mixed layer of soil (e.g., solid and dotted curves in Fig. 10g). The dominance of particles with a relatively small radius results in a large surface area per unit mass, leading to a large surface area difference between the PSD-based and default methods (Fig. 10h and i).

To further explore the impact of PSD tracking on mineral surface area, we conducted another set of experiments identical to those above but with PSD tracking for individual solid species (rather than the bulk solid phase). We present the results of a basalt application experiment in which PSD tracking has the most prominent effects on the surface area calculation (Figs. 11 and 12) and compare these with the default surface area parameterization and the calculation based on the PSD for bulk soil. Individually tracked PSDs are different between solid species (Fig. 11) and also differ from the PSD for bulk soil (Fig. 10g). For example, solid species that are introduced to the system only through basalt powder application show patterns similar to the PSD specified for basalt powder, with smaller particles being produced via dissolution as time proceeds (e.g., Fig. 11k), and those precipitated as a result of basalt powder application show significant particle growth (Fig. 11c, d, n and o). Accordingly, surface area depth profiles vary between solid species and are significantly different from those simulated with the default method and the PSD tracked only for bulk soil (Fig. 12). In the basalt application scenario, surface areas for individual solid species increase more significantly when tracking species-specific PSDs relative to those calculated based on PSD for bulk soil (e.g., Fig. 12k). On the other hand, surface areas of solid species that are not supplied through basalt application are generally smaller when tracking species-specific PSDs than tracking only PSD for bulk soil (e.g., Fig. 12j and m). Note that surface areas of secondary minerals (e.g., Fig. 12d) are not significantly reduced compared with those of dissolving primary minerals that are not supplied through basalt application (Fig. 12j) because surface roughness is accounted for. Overall, when conducting experiments involving addition of particles whose size distributions are significantly different from the PSD of the parent rock, we recommend PSD tracking for individual solid species though computationally more demanding.

3.2 Example ERW application – time-dependent CO₂ capture following basalt addition at contrasting pristine hinterland sites

As a final illustration of model capability, we estimate time-dependent CO₂ capture during continuous and pulsed basalt application at two pristine hinterland watersheds in the US. The two watersheds are obtained by delineating all watersheds corresponding to the USGS river gauges (U.S. Geological Survey, 2016) using GRASS GIS (GRASS Development Team, 2017), followed by selecting watersheds that have experienced minimum human interferences (watersheds that have a population density of less than one person per square kilometer; land proportion of cultivated vegetation less than 1 %; land proportion of urban land less than 1 %). The mineral compositions for these two pristine watersheds are derived from the USGS soil chemistry database (U.S. Geological Survey, 2016). We emphasize that the sites used for illustration are not optimized for CO₂ capture, given that they are characterized by relatively low runoff and temperature (Table 8). The only aim here is to demonstrate the ability of SCEPTER to examine first-order controls on CO₂ capture during enhanced rock weathering, and these sites, given the lack of human intervention, have a more predictable history than sites that have been strongly directly influenced by human activity. Experimental conditions are the same as those adopted in Sect. 3.1.3 and 3.1.4 except for the used sites (Sites 2 and 3; Table 8 and Fig. 5).

3.2.1 Tuned spinup to observational data

When detailed observational data are available for a specific individual site of potential enhanced rock weathering, one can perform site-specific tuning at spinup prior to crushed rock application. As an illustration, we tune SCEPTER to observed OM and pH profiles at two pristine hinterland sites (Sites 2 and 3 in Table 8; see Fig. 5 for their locations) assuming Fickian mixing (e.g., bioturbation) (Fig. 13). For Site 2, optimized OM and pH profiles can be obtained by introducing an additional class of SOM (e.g., Class 3 SOM introduced at 0.1016 times the Class 2 SOM rain flux) and assuming slightly different labilities for the two classes of SOM relative to the default values (5 and 85 years for turnover of SOM Classes 2 and 3, respectively). In tuning to Site 3, no extra SOM class is added, but the rain flux and turnover time are modified for Class 2 SOM to 0.64 g C m⁻² yr⁻¹ and 1600 years, respectively. The simulated soil thickness and mixed layer depth are extended to 2 m to fit to the observed OM and pH profiles, and the total grid number is doubled (60). Surface area is calculated during spinup with the default method assuming 100 and 7 µm of average particle radii for Sites 2 and 3, respectively. The remaining boundary conditions are defined in the same way as detailed in the biotic weathering experiment for Site 1 (Table 7; Sect. 3.1.2) but using the parameter values specified for Sites 2 and 3 in Table 8.

Observed soil OM and pH depth profiles are reproduced well by SCEPTER with minimal tuning (Fig. 13). Fickian mixing is effective only where solid concentrations are relatively high and is thus particularly intense close to the surface. Therefore, Fickian mixing introduces some difficulty in reproducing the SOM profile for Site 2 at depth (∼1 wt %, Table 9). In contrast, because a relatively low SOM reactivity is assumed, Fickian mixing closely reproduces SOM concentrations throughout the soil column for Site 3, except for the low surface value. In any case, this exercise demonstrates that SCEPTER can be readily tuned using site-specific observations in order to optimize ERW capture for local boundary conditions and that in general minimal effort should be required to perform robust site-specific tuning across a very wide range of soil pore fluid pH and organic matter content (Fig. 13).

3.2.2 Time-dependent variation in CO₂ capture

We perform two sets of time-dependent CO₂ capture experiments, applying crushed basalt at a rate of 40 ton ha⁻¹ yr⁻¹ to the tuned spinups discussed in Sect. 3.2.1. In the first set, basalt is applied continuously throughout the simulation, and Fickian mixing is applied across all solid phases. In the second set, we perform pulsed basalt addition in which crushed basalt is applied at the same overall annual application rate and homogeneously mixed into the soil, but only during the initial part of each year (0.1 year), while mixing occurs in the Fickian regime throughout the remainder of the year. We perform all simulations across the range of PSD options in SCEPTER (see above), including the default setting, bulk PSD tracking and full PSD tracking.

We find that the surface area scheme significantly impacts predicted CO₂ capture (Fig. 14). In particular, significant differences between the default parameterization and PSD-based methods (Fig. 14) suggest that explicit differentiation of particle sizes between added powder and existent soil/rock is crucial for accurate predictive quantification of carbon capture during enhanced rock weathering. The results also suggest that species-specific PSD tracking is fundamental particularly for the prediction of short-term CO₂ capture efficiency (<10 years; Fig. 14a and d). It is also clear that the CO₂ capture rate depends strongly on the chosen site and associated background boundary conditions (compare Fig. 14a–c with d–f; see also Tables 8 and 9), particularly over longer time horizons.

Pulsed basalt addition leads to a significant short-term enhancement in CO₂ capture, but on century timescales overall capture rates for both sites are comparable regardless of whether feedstock amendment is continuous or pulsed (Figs. 14 and 15). The long-term CO₂ capture rates estimated here are broadly comparable with previous estimates (e.g., Beerling et al., 2020), but our model predicts significant induction time before achieving maximum CO₂ capture at a given feedstock application rate (e.g., 0.16 capture efficiency (≡ mass ratio of captured CO₂ against deployed basalt) in 100 years at Site 2 vs. 0.15–0.25 capture efficiency in 10 years by Beerling et al., 2020). However, we emphasize that these simulations are only meant to be illustrative of the capabilities of the model and that the sites and ERW deployment procedure have not been optimized for maximizing CO₂ capture rates. As a result, these results should not in any way be interpreted as conclusive or generalizable with respect to CO₂ capture rate, efficiency and induction time. In addition, the potential difference in induction time for CDR should be evaluated in a more comprehensive intercomparison study in which model parameterizations, boundary conditions and experimental setups are aligned as much as possible. The purpose here is simply to illustrate the differences between parameterizations of particle size and surface area.

4 Conclusions

SCEPTER is a traceable code with the capability to comprehensively realize phenomena occurring within soil weathering systems, including abiotic/biotic weathering of minerals, mixing of soil particles and addition of OM/minerals under natural or managed conditions. The model is equipped with options to calculate surface areas of soil particles by tracking particle size distributions through time and space. This specific feature may be of particular importance for calculating the cost performance of terrestrial enhanced rock weathering, as a significant component of both economic cost and secondary CO₂ emissions is the grinding and transport of rock feedstocks (e.g., Renforth, 2012). Application of the model to US soil data indicates that it is well suited to capturing the major mineral transformations and solute fluxes attendant to natural weathering.

The current version of the model can serve several purposes relevant for natural weathering, carbon cycle dynamics and managed soil chemistry. Ongoing model developments include a full mechanistic implementation and validation of nutrient (P, N, and K) cycles and nutrient uptake and secretion of organic acids by plants (e.g., Lawrence et al., 2014; Perez-Fodich and Derry, 2019), inclusion of size-dependent changes in the particle solubility and nucleation kinetic barrier when tracking species-specific PSDs (e.g., Hochella, 2003; Perez et al., 2008; Emmanuel and Ague, 2011), implementation of a wider range of potential management practices and coupling to Earth system model frameworks in order to more fully explore the dynamics of CDR (e.g., Ridgwell et al., 2007; Holden et al., 2016; Taylor et al., 2016). We suggest that SCEPTER can be a powerful tool for predicting and diagnosing the effects of enhanced rock weathering on existing soil systems, particularly following robust site-specific tuning against background observations.