the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
SedTrace 1.0: a Juliabased framework for generating and running reactivetransport models of marine sediment diagenesis specializing in trace elements and isotopes
Trace elements and isotopes (TEIs) are important tools in studying ocean biogeochemistry. Understanding their modern ocean budgets and using their sedimentary records to reconstruct paleoceanographic conditions require a mechanistic understanding of the diagenesis of TEIs, yet the lack of appropriate modeling tools has limited our ability to perform such research. Here I introduce SedTrace, a modeling framework that can be used to generate reactivetransport code for modeling marine sediment diagenesis and assist model simulation using advanced numerical tools in Julia. SedTrace enables mechanistic TEI modeling by providing flexible tools for pH and speciation modeling, which are essential in studying TEI diagenesis. SedTrace is designed to solve one particular challenge facing users of diagenetic models: existing models are usually casespecific and not easily adaptable for new problems such that the user has to choose between modifying published code and writing their own code, both of which demand strong coding skills. To lower this barrier, SedTrace can generate diagenetic models only requiring the user to supply Excel spreadsheets containing necessary model information. The resulting code is clearly structured and readable, and it is integrated with Julia's differential equation solving ecosystems, utilizing tools such as automatic differentiation, sparse numerical methods, Newton–Krylov solvers and preconditioners. This allows efficient solution of large systems of stiff diagenetic equations. I demonstrate the capacity of SedTrace using case studies of modeling the diagenesis of pH as well as radiogenic and stable isotopes of TEIs.
 Article
(2802 KB)  Fulltext XML
 BibTeX
 EndNote
Trace elements and isotopes (TEIs) play critical roles in the heathy functioning of the marine ecosystem (SCOR Working Group, 2007). For example, Fe, Zn, Ni and other transition metals serve as cofactors in metalloenzymes that are essential for phytoplankton physiology such that the low concentrations of these metals can limit phytoplankton productivity with implications for the global carbon cycle (Martin et al., 1987; Vance et al., 2017; Tagliabue et al., 2017; Morel et al., 2020; Lemaitre et al., 2022). TEIs are also powerful tracers of geological, physical and biogeochemical processes in the ocean (Lam and Anderson, 2018). For instance, the radiogenic isotope of Nd has been used to study ocean circulation and the role of continental weathering and tectonics in past climate changes (Piepgras et al., 1979; Palmer and Elderfield, 1985; Frank, 2002; Goldstein and Hemming, 2003; Du et al., 2020); the stable isotopes of Mo have been used to study past ocean oxygenation (Archer and Vance, 2008; Anbar and Rouxel, 2007; Lyons et al., 2009). Since the dawn of the international GEOTRACES program, there have been increasing efforts aimed at characterizing the distributions, global oceanic budgets, internal cycling and external sources of TEIs in the modern ocean and their utility in paleoceanography (Anderson, 2020; Schlitzer et al., 2018).
In recently years, however, two issues have emerged as major challenges to the understanding and application of marine TEIs. First, the modern oceanic budgets of many TEIs, such as Fe, Cu, Ni, Zn and Nd, cannot be balanced by known sources at the riverine, atmospheric and hydrothermal interfaces, and increasing evidence suggests that fluxes across the sediment–water interface (SWI) may play important, if not dominant, roles in setting the global oceanic budgets of TEIs (Homoky et al., 2016; Jeandel, 2016; Little et al., 2014; Du et al., 2020; Haley et al., 2017; Elrod et al., 2004; Cameron and Vance, 2014; Little et al., 2020). Second, applying TEIs as proxies to study the geological evolution of the ocean system is hampered by the lack of quantitative knowledge of the sedimentary cycling of TEIs, which may alter or completely erase the original proxy signals preserved in sedimentary archives (Horner et al., 2021; Crusius and Thomson, 2000). Models of the sedimentary cycling of TEIs are thus needed to resolve these issues. Quantitative modeling is particularly indispensable in this case because the complexity, heterogeneity, and highly coupled nature of biogeochemical and physical processes in marine sediments make straightforward interpretation of measured modern and paleoTEI data difficult. In addition, TEI data in sediment pore water are extremely scarce because of analytical and sampling challenges. I thus hope that models can assist our understanding of the sedimentary TEI cycling in the presence of large data gaps.
A diagenetic model is typically a 1D reactivetransport model that includes physical transport and biogeochemical reactions to simulate the distributions of dissolved substances in pore water and solid substances in sediments (Berner, 1980; Boudreau, 1997; Burdige, 2006). Such models have traditionally been used to study the sedimentary cycling of carbon, oxygen and nutrients (Wang and Van Cappellen, 1996; Boudreau, 1996; Meysman et al., 2003; Soetaert et al., 1996), as well as to a limited extent TEIs like Fe, Mn and U (Dale et al., 2015; Burdige, 1993; Lau et al., 2020; Maher et al., 2006). However, to date there has been little concerted effort to develop diagenetic models specifically for TEIs due to some particular challenges of TEI biogeochemistry. First, a realistic digenetic TEI model needs to include a pH module to enable speciation modeling. Pore water pH is difficult to model because of the large network of biogeochemical and physical processes involved (Boudreau and Canfield, 1988, 1993; Jourabchi et al., 2008; Reimers et al., 1996). Many studies thus ignore TEI speciation in diagenetic models. Moreover, TEIs are commonly influenced by more biogeochemical reactions than the major and minor constituents. Together with the necessity for pH and speciation modeling, a diagenetic model may need to include a large system of differential equations that are also highly stiff because the span of the reaction timescales can be large (Boudreau, 1997). Thus, unlike traditional diagenetic models of carbon, oxygen and nutrients, diagenetic TEI models also need to take into account numerical efficiency, especially if the goal is to couple diagenetic models to global ocean biogeochemical models (Hülse et al., 2018; Archer et al., 2002).
The greatest challenge from the user's point of view is that most diagenetic models are too specialized to be adaptable. The strong heterogeneity of marine sediments implies that no single diagenetic model can be created to include all sedimentary processes applicable to all environmental settings (Paraska et al., 2014). Often, modelers create specific diagenetic models with fixed selection of substances and processes for their own studies. The model code is often not opensource, and it could be timeconsuming and errorprone for other users to adapt the code for new studies. Writing and modifying code require strong programming skills, which is a significant hurdle to the general user. Thus, rather than creating one allencompassing diagenetic model, it is preferable to create a modeling framework that can generate custom models, allowing users with limited programming skills to create and run models satisfying their own needs (Soetaert and Meysman, 2012).
In this study, I describe SedTrace 1.0, a modeling framework that automates the generation of Julia code for diagenetic models and provides highperformance computing tools to assist model simulation, only requiring the user to supply information using a Microsoft Excel input file. SedTrace specializes in modeling TEIs and can accommodate a wide range of custom pH and speciation modeling choices. The design principle is to give as much control to the user as possible when generating the model such that SedTrace is only meant to help convert user ideas to code rather than make model decisions for the user. Julia is an opensource, dynamically typed programming language for highperformance scientific computing (Bezanson et al., 2017). It offers the high productivity of scripting languages like Python, but can also match the performance of statically typed languages such as C and FORTRAN. It is increasingly being adopted by the climate and ocean modeling community (Pasquier et al., 2022; Sridhar et al., 2022; Sulpis et al., 2022), with wellsupported ecosystems relevant to solving differential equations (Rackauckas and Nie, 2017).
This paper is structured as follows. First I describe the model equations and framework in Sect. 2. I then discuss the code generation procedure for physical processes and biogeochemical reactions in Sects. 3 to 5. The numerical solution of the model is discussed in Sects. 6 and 7. Finally, I present a few case studies in Sect. 8 and briefly touch on future development in Sect. 9.
SedTrace uses the 1D diagenetic equation (Boudreau, 1997; Meysman et al., 2003):
where ${C}_{i}^{\mathit{\xi}}$ is the concentration of model substance i in the phase ξ (f for pore water or s for solid sediment), ϕ^{ξ} is the phase volume fraction, T_{i} is the transport due to advection, diffusion and bioirrigation, R_{i} is the net source or sink due to biogeochemical reactions, t is time, and x is the sediment depth coordinate staring from the SWI pointing downward. In SedTrace, the default units are millimoles (mmol) for mass, years for time and centimeters (cm) for length. For example, reaction rates are in units of millimoles per cubic centimeter per year (mmol cm^{−3} yr^{−1}).
The transport terms have wellestablished forms in diagenetic models (Berner, 1980; Boudreau, 1997). It is largely the reaction terms which are casespecific that create the diversity of diagenetic models and thus cause their poor adaptability. SedTrace fixes the general forms of the transport terms and only requires users to supply casespecific parameters, while letting the user set the reaction terms freely.
Using userspecified spatial grids, which can be nonuniform, SedTrace discretizes Eq. (1) using the method of lines in the time dimension and finitevolume method in the space dimension (Boudreau, 1996), resulting in a system of ordinary different equations (ODEs) of time only:
where j is the index of the grid cell. ${F}_{i,j\frac{\mathrm{1}}{\mathrm{2}}}^{\mathit{\xi}}$ and ${F}_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathit{\xi}}$ are the fluxes across the left and right boundaries of the cell, respectively. Δx_{j} is the volume of the cell.
The diffusive and advective fluxes are discretized using the finitevolume–complete flux (FVCF) scheme (Llorente et al., 2020; ten Thije Boonkkamp and Anthonissen, 2010, 2011). The finitevolume method is preferred in diagenetic models because it ensures mass conservation, whereas the finitedifference method may not.
The finitevolume flux ${F}_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathit{\xi}}$ is derived based on the analytical solution of the local twopoint boundary value problems and is partitioned into a homogeneous flux and an inhomogeneous flux. Presently SedTrace only uses the homogeneous flux term. The symbolic derivation of these fluxes using Mathematica can be found in /math/fvcf.pdf
.
The resulting discretization may be viewed as a weighted average of the centered scheme and the upwind scheme, and the weight is controlled by the cell Péclet number $\frac{{v}^{\mathit{\xi}}\mathrm{\Delta}x}{\mathrm{2}{D}^{\mathit{\xi}}}$, where v^{ξ} is the advective velocity and D^{ξ} is the diffusion coefficient. At high Péclet numbers (advection dominant), the scheme approaches the upwind scheme, while at low Péclet numbers (diffusion dominant) it approximates the centered scheme. This ensures at least firstorder uniform convergence regardless of the Péclet number. Similar Pécletnumberdependent schemes have been applied to diagenetic and ocean modeling (Fiadeiro and Veronis, 1977; Soetaert and Meysman, 2012; Boudreau, 1996). The numerical diffusion introduced in advectiondominated cases has a greater impact on transient than steadystate simulations (ten Thije Boonkkamp and Anthonissen, 2011, 2010). Thus, SedTrace 1.0 is best suited for steadystate modeling of tracer profiles in sediments rather than transient simulations of paleoproxy evolution.
SedTrace collects the discretized model substances into an MN vector:
where M is the number of model substances and N is the number of grid cells. The system of ODEs including boundary conditions can be written in the matrix form:
A is an MN×MN matrix including the linear diffusion and advection operators. B is an MN×MN matrix containing the homogeneous part of the boundary conditions. b is an MN vector containing the nonhomogeneous part of the boundary conditions. S is an MN vector, generally a nonlinear function of C, incorporating the reaction and bioirrigation sources and sinks, as well as nonlinear coupling of transport.
To generate code for Eq. (5) the user supplies an Excel file model_config.xlsx
to SedTrace. An example of a simple Fe cycle model (SimpleFe
) can be found in /examples/SimpleFe/model_config_simpleFe
.xlsx
, and the spreadsheets are shown in Tables 1 to 7. The substances
sheet (Table 1) lists the modeled substances, their types (e.g., solid or dissolved), chemical formula and boundary conditions. The reactions
sheet (Table 2) lists the kinetic reactions, their chemical equations and rate expressions. The speciation
sheet (Table 3) lists aqueous speciation reactions. The adsorption
sheet (Table 4) lists the adsorbed species. The diffusion
sheet (Table 5) lists information to compute the diffusion coefficients of dissolved substances. The parameters
sheet (Table 6) lists the parameters required by the model. The output
sheet (Table 7) is used to formulate output and plotting. The data
sheet (not shown here) includes observational data that will be plotted together with model outputs.
The workflow of generating and running diagenetic models using SedTrace is shown in Fig. 1. In Sects. 3 and 4, I will describe the mathematical formulation of each of the terms in Eqs. (1), (2) and (5), as well as how SedTrace generates the corresponding code. I will use SimpleFe
to illustrate this process. In this paper, I use monospaced font to indicate computer code; variable names in the code, Excel spreadsheet and column names; and Julia file names.
^{a} Code name for the model substance. ^{b} Chemical formula that can be used to write chemically balanced equations in the reactions
sheet. ^{c} For model substances of type dissolved_speciation
, the boundary conditions listed here are for the total dissolved concentration.
^{a} Chemical balance is checked if check
is set to 1. ^{b} Equations can be written using the code name or the chemical formula of model substances or the code names and formulae of their dissolved and adsorbed species. ^{c} Rate and omega expressions can be written only using the code name of the model substances or the code names of their dissolved and adsorbed species.
^{a} Code name of the dissolved species; it must be different from the names of model substances in the substances
sheet. ^{b} Chemical formula that can be used to write the chemical equations; it must be different from the formulae of model substances in the substances
sheet. ^{c} Equations should be written using the formula, not code name. ^{d} Set to 1 if the concentration of the species is required.
^{a} Code name of dissolved species that appears in the speciation
sheet or the code name of the total dissolved concentration. ^{b} Code name of adsorbed species. ^{c} Surface to be adsorbed onto; it can be either the code name of a solid substance from the substances
sheet or left empty. ^{d} Expression can be written only using the code names of the dissolved and adsorbed species listed in the same row; adsorption parameters need to be supplied to the parameters
sheet. ^{e} The boundary conditions of the adsorbed species, which should be the same for all adsorbed species of the same substance.
^{a} Code name of the dissolved model substance or its total dissolved concentration. ^{b} The corresponding name listed in SedTrace's database of the diffusion coefficients.
^{a} Class must be specified using one of the key words shown here following the same order. ^{b} Type can be constant or a function of depth x. ^{c} Code name for parameters. ^{d} To enter the value for the parameter, either supply a numerical value, a function of depth, or a string expression that computes the value using other parameters, in which case the parameters being depended on must appear earlier in the table. The values shown here for the SimpleFe
model are only for illustration. ^{e} Optional.
^{a} Code name for output variables; this can be the same as the model substances or new variable names. ^{b} Expressions to compute new variables if needed; the expression can use model parameters listed in parm.jl
. ^{c} Multiplication factors that convert default SedTrace units to custom units in unit_profile
; these can use model parameters listed in parm.jl
– must not be empty (use 1 instead). ^{d} TA is computed internally and does not need an expression.
3.1 Grid
To generate the model grid, the user specifies the number of grids (Ngrid
), the length of the sediment domain (L
) and a grid transformation function (gridtran
) in the parameters
sheet. SedTrace internally creates a uniform grid between 0 and L
cm and uses gridtran
to transform it to the userdesired grid (using function x>x
will preserve the original uniform grid). A nonuniform grid is used to reduce numerical error at locations where finer spacing is needed (Meysman et al., 2003): for example, close to the SWI where biogeochemical gradients are often the sharpest. SedTrace allows the user to supply parameters as constants or functions of depth x, as labeled in the type
column of the parameters
sheet (Table 6): for example, gridtran
as a function in this case. SedTrace will convert the function string to Julia function and use broadcast()
to vectorize the function if necessary. The generated model grid code for SimpleFe
using a gridtran
that has small grid spacing closer to the SWI is as follows.
#
# grid parameters
#
L = 50.0 # cm # model sediment section thickness
Ngrid = 200 # integer # number of model grid
ξ = range(0, step = L / (Ngrid), length = Ngrid + 1)
# cm # uniform grid
x_{v} = broadcast(x > L * (exp(x / 5)  1) / (exp(L / 5)  1), ξ) # cm # nonuniform grid transformation
x = (x_{v}[2:(Ngrid + 1)] . +x_{v}[1:Ngrid]) / 2 # cm # cell center
dx = x_{v}[2:(Ngrid+1)] . x_{v}[1:Ngrid] # cm # cell volume
3.2 Advection
The advective flux in Eq. (2) is
where v^{ξ} is the phase velocity. To calculate v^{ξ}, SedTrace makes two assumptions (Boudreau, 1997; Meysman et al., 2003; Berner, 1980): sediment compaction is at steady state and therefore the volume fractions of fluid ϕ^{f} and solid ${\mathit{\varphi}}^{\mathrm{s}}\left(=\mathrm{1}{\mathit{\varphi}}^{\mathrm{f}}\right)$ are functions of depth but not time, and the final burial velocities of fluid and solid are the same, ${v}_{\mathrm{\infty}}^{\mathrm{f}}={v}_{\mathrm{\infty}}^{\mathrm{s}}$, without externally forced pore water advection. Using the usersupplied porosity function ϕ^{f} (phif
), porosity at burial depth ${\mathit{\varphi}}_{\mathrm{\infty}}^{\mathrm{f}}$ (phif_Inf
), density of dry sediments ρ^{s} (ds_rho
) and solid burial flux ${F}_{\mathrm{total}}^{\mathrm{s}}$ (Fsed
) in parameters, SedTrace computes the phase velocities:
In SimpleFe
the resulting code is as follows.
#
# porosity parameters
#
phi_Inf = 0.7884 # dimensionless # porosity at infinite depth
phif = broadcast(x > 0.8 + (0.9  0.8) * exp(x / 2), x) # dimensionless # fluid volume fraction
phis = 1.0 . phif # dimensionless # solid volume fraction
#
# phase velocity parameters
#
Fsed = 0.073 # g cm^2 yr^1 # total sediment flux
w_Inf = Fsed / ds_rho / (1  phi_Inf) # cm yr^1 # solid sediment burial velocity at infinite depth
uf = phi_Inf * w_Inf ./ phif # cm yr^1 # pore water burial velocity
us = Fsed / ds_rho ./ phis # cm yr^1 # solid sediment burial velocity
3.3 Diffusion
For solid substances SedTrace considers the diffusive flux due to bioturbation (Boudreau, 1997):
where D^{b} (Ds
) is the bioturbation coefficient as a function of depth specified in parameters
.
#
# bioturbation parameters
#
Ds = broadcast(x > 10 * exp(x / 3), x) # cm^2 yr^1 # Bioturbation coefficient
For dissolved substances SedTrace considers molecular diffusion ${D}_{i}^{\mathrm{md}}$ corrected for the tortuosity factor θ^{2} (Boudreau, 1997):
Boudreau (1997) parameterized ${D}_{i}^{\mathrm{md}}$ at infinite dilution and atmospheric pressure (P_{atm}) as linear functions (m_{0}+m_{1}T) of temperature (T in Celsius) for selected dissolved substances. In this case SedTrace computes ${D}_{i}^{\mathrm{md}}$ at userspecified in situ T, salinity (S) and pressure (P) using the Stokes–Einstein relationship (Li and Gregory, 1974):
where μ is the dynamic viscosity of pore water as a function of T, S and P. However, if only the diffusion coefficient ${D}_{i}^{{\mathrm{md}}^{\mathrm{25}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}}}$ at infinite dilution, 25 ^{∘}C and atmospheric pressure is known, SedTrace computes
The values of m_{0}, m_{1} and ${D}_{i}^{{\mathrm{md}}^{\mathrm{25}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}}}$ of selected dissolved substances compiled by Boudreau (1997) are stored in the database file /src/diffusion.xlsx
. If a model substance is in the database, then the user simply needs to list it in the diffusion
sheet, where their model_name
given by the user needs to match the SedTrace_name
in the database (Table 5). SedTrace will compute ${D}_{i}^{\mathrm{md}}$ automatically. If the model substance is not in the database, the user can either modify the database file to include it or directly supply ${D}_{i}^{\mathrm{md}}$ to parameters
.
The example code for diffusion coefficient is as follows.
#
# solute diffusivity
#
DSO4 = 1.8034511184582917E+02 ./ (1.0 . 2log.(phif)) # cm^2 yr^{1} # Sediment diffusion coefficient
3.4 Total physical transport
SedTrace collects the discretized advective and diffusive fluxes as an interior transport matrix A (Am
) by calling the fvcf()
function, which performs the FVCF discretization. The AC term in Eq. (5) is then computed. The code for the transport of SO_{4} (AmSO4
) is as follows.
#
# Interior transport matrix
#
AmSO4 = fvcf(phif, DSO4, uf, dx, Ngrid) # # Interior transport matrix of SO4
#
# Transport term A*C
#
mul!(dSO4, AmSO4, SO4) # dSO4 is the rate of change
3.5 Bioirrigation
Bioirrigation sources and sinks of a dissolved substance are modeled as a nonlocal exchange (Boudreau, 1997):
where α (alpha
) is the bioirrigation coefficient as a function of depth x and ${C}_{i}^{\mathrm{BW}}$ is the bottomwater concentration, both supplied to parameters
. Previous studies have suggested that the bioirrigation coefficient α is also substancespecific (Meile et al., 2005), but it is unknown how this effect should be parameterized. SedTrace thus does not consider it at the moment.
The code for the biological transport term of SO_{4} is as follows.
#
# bioirrigation parameters
#
alpha = broadcast(x > 10 * exp(x / 2), x) # yr^{1} # Bioirrigation coefficient
#
# biological transport
#
@.. dSO4 += alpha * (SO4BW  SO4)
3.6 Boundary conditions
At the SWI, for a solid substance (${C}_{i}^{\mathrm{s}})$ users can specify a Robin boundary condition using the incoming flux ($F{C}_{i}^{{\mathrm{s}}^{\mathrm{0}}}$) (Boudreau, 1997),
or a Dirichlet boundary condition using the concentration at the SWI (${C}_{i}^{{\mathrm{s}}^{\mathrm{0}}}$):
At the SWI, for a dissolved substance (${C}_{i}^{\mathrm{f}}$) users can specify a Robin boundary condition assuming the existence of a diffusive boundary layer (DBL) of thickness δ and the overlying bottomwater concentration of ${C}_{i}^{\mathrm{BW}}$ (Boudreau, 1997),
or a Dirichlet boundary condition using the concentration at the SWI (${C}_{i}^{{\mathrm{f}}^{\mathrm{0}}}$):
At the bottom of the model domain (x=L), users can specify either a Neumann or a Dirichlet boundary condition for model substances:
Here Eq. (20) assumes no concentration gradient at x=L.
The user specifies the type of boundary conditions in substances
and provides the necessary parameters in parameters
(Table 6). SedTrace will format the boundary conditions as
For example, setting the Robin boundary condition in Eq. (18) at the SWI for SO_{4} (Table 1), SedTrace will compute $\mathit{\beta}=\frac{{D}_{j}^{i,\mathrm{md}}}{\mathit{\delta}}$ (beta
), which is the mass transfer velocity, as well as ${\mathit{\alpha}}_{\mathrm{1}}^{\mathrm{0}}=\mathit{\beta}+{\mathit{\varphi}}^{\mathrm{f}}{v}^{\mathrm{f}}$, ${a}_{\mathrm{2}}^{\mathrm{0}}={\mathit{\varphi}}^{\mathrm{f}}\frac{{D}_{j}^{i,\mathrm{md}}}{{\mathit{\theta}}^{\mathrm{2}}}$ and ${a}_{\mathrm{3}}^{\mathrm{0}}=\mathit{\beta}{C}_{i}^{\mathrm{BW}}$ using δ (delta
) and ${C}_{i}^{\mathrm{BW}}$ (SO4BW
) from the parameters
. The Neumann bottom boundary condition is simply ${\mathit{\alpha}}_{\mathrm{1}}^{L}=\mathrm{0}$, ${a}_{\mathrm{2}}^{L}=\mathrm{1}$ and ${a}_{\mathrm{3}}^{L}=\mathrm{0}$. SedTrace collects these coefficients into a Tuple
((${a}_{\mathrm{1}}^{\mathrm{0}},{a}_{\mathrm{2}}^{\mathrm{0}},{a}_{\mathrm{3}}^{\mathrm{0}})$, (${a}_{\mathrm{1}}^{L},{a}_{\mathrm{2}}^{L},{a}_{\mathrm{3}}^{L}\left)\right)$, which is passed to fvcf_bc()
to generate the homogeneous (BcAM
, B) and nonhomogeneous (BcCm
, b) parts of the boundary condition based on the FVCF discretization. SedTrace then updates the righthand side of Eq. (5) by adding the BC+b term. The code for SO_{4} is as follows.
#
# assemble boundary conditions
#
BcSO4 =
(
(betaSO4 + phif[1]uf[1], phif[1]DSO4[1], betaSO4 * SO4BW),
(0.0, 1.0, 0.0),
) # # Boundary condition of SO4
#
# Boundary transport matrix
#
BcAmSO4, BcCmSO4 = fvcf_bc(phif, DSO4, uf, dx, BcSO4, Ngrid) # # Boundary transport matrix of SO4
#
# Boundary terms
#
dSO4[1] += BcAmSO4[1] * SO4[1] + BcCmSO4[1]
dSO4[Ngrid] += BcAmSO4[2] * SO4[Ngrid] + BcCmSO4[2]
Biogeochemical reactions in diagenetic models can be classified as kinetic or equilibrium reactions (Boudreau, 1997; Meysman et al., 2003). The difference mainly lies in the reaction timescale: reactions that happen on much shorter timescales than physical transport are often treated as equilibrium rather than kinetic reactions. The user is free to add the kinetic reactions, while SedTrace has a specific user interface for two types of equilibrium reactions: acid dissociation for pH modeling and aqueous complexation–sorption for speciation modeling. SedTrace uses the direct substitution approach (DSA) to handle the equilibrium reactions, which performs better than other approaches when dealing with a highly coupled and stiff biogeochemical reaction network (Meysman et al., 2003).
4.1 Kinetic reactions
The summed rate of kinetic reactions for substance ${C}_{i}^{\mathit{\xi}}$ is
where r^{η,k} is the rate of the kth reaction in units per volume of phase η, and ${\mathit{\nu}}_{i}^{\mathit{\eta},k}$ is the stoichiometry of the ith substance in this reaction. The unit of R_{i} is per volume of phase ξ. The convention of SedTrace is that homogeneous reactions between dissolved substances are in units per volume of fluid, while heterogeneous reactions between dissolved and solid substances and homogenous reactions between the solid substances are in units per volume of solid. Conversion between fluid and solid concentration units is done using the conversion factor ${\mathit{\varphi}}^{\mathrm{f}}/{\mathit{\varphi}}^{\mathrm{s}}$ (fluid to solid, pwtods
) or ${\mathit{\varphi}}^{\mathrm{s}}/{\mathit{\varphi}}^{\mathrm{f}}$ (solid to fluid, dstopw
). Kinetic reactions are added to the reactions
sheet, including their chemical equations and rate expressions. SedTrace parses the reactions and reassembles them into Julia code using Julia's Perlcompatible regular expression engine and the symbolic computing utility from the SymPy.jl package (Meurer et al., 2017).
4.1.1 Chemical equations
In SedTrace each model substance and its dissolved and adsorbed species need to have a code name and a chemical formula, supplied to the Excel sheets as noted in the tables. The code names are used in the code, and I suggest only using the Latin alphabet and “_” in the name. The formulae are necessary only when writing chemically balanced reaction equations. The code name and the chemical formula can be the same if (and only if) the chemical formula is an allowed Julia variable name, such as FeOOH
but not Fe{2+}
. The chemical formula of a substance should be written in the form of (E)[e](F)[f](G)[g]{h+}
, where E/F/G
are compounds, e/f/g
are their stoichiometric coefficients and h
is the number of charge. For example, the user can name organic matter POM
and write its chemical formula as (CH2O)(NH3)[rNC](H3PO4)[rPC]
, where rNC
and rPC
are the $\mathrm{N}/\mathrm{C}$ and $\mathrm{P}/\mathrm{C}$ ratios. SedTrace allows parameterized stoichiometry, and the parameters rNC
and rPC
should be supplied to parameters
. This feature makes it easy to perform model sensitivity tests since natural chemical compounds such as organic matter in marine sediments rarely have fixed stoichiometry.
SedTrace requires the user to write chemical equations in the form of a*A + b*B = c*C + d*D
, where a/b/c/d
are the stoichiometry coefficients of substances A/B/C/D
, and the reactants are on the lefthand side. The user has two options when writing these equations, decided by the check
column in the reactions
sheet (Table 2), which informs SedTrace if the chemical balance of the equation should be checked. The user can write a chemically balanced equation using the chemical formulae of the model substances A/B/C/D
and set check
to 1. Or the user can leave check
empty and write a heuristic equation without considering chemical balance using code names only or a mixture of code names and chemical formulae. SedTrace offers the second choice because it is common in the diagenetic literature to write heuristic reactions. This happens because of the complexity of biogeochemical reactions in marine sediment such that it is often not possible to know the exact chemical formula or reaction mechanism. An example heuristic chemical equations is POM = Carbon + rNC*NH3 + rPC*PO4
for organic matter remineralization.
SedTrace uses regular expressions to split the chemical equation into reactants and products, as well as identify the stoichiometry coefficients ${\mathit{\nu}}_{i}^{\mathit{\eta},k}$ and charges of the model substances. SedTrace further parses the equation down to the level of individual elements, identifying the stoichiometry of each element. The results are saved in model_parsing_diagnostics.xlsx
that the user should check to make sure the parsing is correct. Here I show the parsing result of a more complex version of the Fe reduction reaction in Table 2. (CH2O)(NH3)[rNC](H3PO4)[rPC] + 4*FeOOH + (8+rNCrPC)*H{+} = CO2 + rNC*NH4{+} + rPC*H2PO4{} + 4*Fe{2+} + 7*H2O
:
#
# Parsing the chemical equation into reactants/products
#
Row 
species 
stoichiometry 
charge 
role 

String 
String 
String 
String 

1 
(CH2O)(NH3) 
1 
0 
reactant 
[rNC](H3PO4) 

[rPC] 

2 
FeOOH 
4 
0 
reactant 
3 
H 
(8+rNCrPC) 
+1 
reactant 
4 
CO2 
1 
0 
product 
5 
NH4 
rNC 
+1 
product 
6 
H2PO4 
rPC 
1 
product 
7 
Fe 
4 
+2 
product 
8 
H2O 
7 
0 
product 
#
# Parsing the chemical equation into elements
# negative value indicates reactant
#
Row 
element 
coef 

String 
String 

1 
C 
1 
2 
H 
2 
3 
O 
1 
4 
N 
(rNC) 
5 
H 
(3*rNC) 
6 
H 
(3*rPC) 
7 
P 
(rPC) 
8 
O 
(4*rPC) 
9 
Fe 
4 
10 
O 
4 
11 
O 
4 
12 
H 
4 
13 
H 
(rNC + rPC  8) 
14 
C 
1 
15 
O 
2 
16 
N 
(rNC) 
17 
H 
(4*rNC) 
18 
H 
(2*rPC) 
19 
P 
(rPC) 
20 
O 
(4*rPC) 
21 
Fe 
4 
22 
H 
14 
23 
O 
7 
Since the stoichiometry may be parameterized, SedTrace uses symbolic computing to check the chemical balance if check = 1
. Parameters like rNC
and rPC
are converted to Julia symbols of type SymPy.Sym
. The sums of charge and mass are computed symbolically, and errors are thrown if the sums are not zero. If check
is empty SedTrace will parse the equation and identify the stoichiometric coefficients, but it will not check the chemical balance.
4.1.2 Reaction rates
The string expressions of reaction rates, r^{η,k} in Eq. (23), supplied to the reactions sheet are copied and pasted as Julia code directly, and therefore they should be written only using the code names. For dissolution and precipitation reactions, the user can add the definition of saturation state in a separate column: Omega
. SedTrace will allow dissolution or precipitation only when Omega < 1
or Omega > 1
, respectively. However, an if Omega > 1
, then the precipitate statement, i.e., a Heaviside step function, in the code induces numerical discontinuity, which can hurt numerical performance, especially when applying automatic differentiation to the code. SedTrace provides the option of approximating the Heaviside step function using the logistic function, which is continuously differentiable:
where r_{pre} is a precipitation rate and the parameter τ controls how close the approximation is. By default τ=10^{3}, which results in a sharp transition near saturation (Fig. 2). For example, the code for FeS
precipitation rate RFeS_pre
is as follows.
@.. Omega_RFeS_pre = Fe_eq * HS / (H * KspFeS)
# saturation state for precipitation
@.. RFeS_pre =
(tanh(1e3 * (Omega_RFeS_pre  1.0)) / 2 + 0.5) *
(kFeSpre * (Omega_RFeS_pre  1))
# precipitation rate
This option can be disabled during code generation (see Sect. 7).
To compute the summed reaction rate R_{i} in Eq. (23), SedTrace collects the stoichiometry coefficients ${\mathit{\nu}}_{i}^{\mathit{\eta},k}$ of the model substances in each kinetic reaction after parsing the chemical equations and applies the appropriate unit conversion factors (dstopw
or pwtods
). For example, code for the summed reaction rate S_FeOOH
of Fe oxide is as follows.
@.. S_FeOOH = 4 * RFeOOHPOC + 2
* RFeOOHH2S
4.2 pH modeling
Traditionally diagenetic models enable pH modeling by combining the differential equations describing the time evolution of model substances together with the pHdependent nonlinear algebraic equations of charge or alkalinity balance, forming a system of differential algebraic equations (Wang and Van Cappellen, 1996; Boudreau, 1996). However, this approach is not only numerically difficult, but it also does not give direct information on how model processes affect pH. SedTrace models pH using the DSA outlined by Hofmann et al. (2008, 2009, 2010). Under this approach, proton concentration is modeled dynamically, avoiding the challenge of solving DAEs while allowing detailed partitioning of pH changes to individual reaction and transport processes. The dynamic equation for proton concentration ([H^{+}], free scale) is
where TA is the total alkalinity (TA), and EI_{l} is the lth equilibrium invariant (EI). The EIs are composite variables, such as the total dissolved inorganic carbon (TCO_{2}). They are so named because they are invariant with respect to the equilibrium reaction rates.
The full definition of TA in seawater (Dickson et al., 2007) is
To include all these species, SedTrace provides the following EIs:
Except for protons, the concentration in these equations refers to the total concentration, i.e., the sum of the concentrations of free and complexed species. It is usually unnecessary to include the entire set in diagenetic models. The user is free to choose any subset of this collection. The user adds the selected EIs to the substances
sheet and specifies the type as dissolved_pH
. Apart from supplying the boundary conditions, all other aspects of pH modeling are handled internally by SedTrace, requiring no user input. Based on the user's choice, SedTrace defines TA as a subset of the full definition in Eq. (26). For example, in SimpleFe
TCO_{2} and TH_{2}S are selected (Table 1), and then $\text{TA}=\left[{\mathrm{HCO}}_{\mathrm{3}}^{}\right]+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}\right]+\left[{\mathrm{HS}}^{}\right]+\left[{\mathrm{OH}}^{}\right]\left[{\mathrm{H}}^{+}\right]$. H^{+} and OH^{−} are always included by default.
SedTrace uses EquilibriumInvariant
to store the information of the EIs, including the analytical expressions to compute the concentrations of the individual species listed in Eqs. (27) to (34), their coefficients in Eq. (26), and the expressions to compute $\frac{\partial \text{TA}}{\partial {\text{EI}}_{l}}$ and $\frac{\partial \text{TA}}{\partial \left[{\mathrm{H}}^{+}\right]}$, for example TCO_{2}.
struct EquilibriumInvariant
name::String # name of EI
species::Vector{String} # species
charge::Vector{String} # charges of the species
expr::Vector{String} # expression to compute the species concentration
coef::Vector{String} # coefficient of species in TA definition
dTAdEI::String # expression to compute dTA/dEI
dTAdH::String # expression to compute this EI's
contribution to dTA/dH
diss_const::VectorString # acid dissociation
constants
end
# TCO2 for example
EquilibriumInvariant(
"TCO2",
["HCO3", "CO3", "CO_{2}"],
["{}","{2}",""],
[ "H * KCO2 * TCO2 / (H^2 + H * KCO2
+ KCO2 * KHCO3)",
"KCO2 * KHCO3 * TCO2 / (H^2 + H * KCO2 + KCO2 * KHCO3)",
"H^2 * TCO2 / (H^2 + H * KCO2 + KCO2 * KHCO3)",
],
["1", "2", "0"],
"KCO2*(H + 2*KHCO3)/(H^2 + H*KCO2 + KCO2
* KHCO3)",
"KCO2*TCO2*(H^2 + 4*H*KHCO3 + KCO2*KHCO3)/
(H^2 + H*KCO2 + KCO2*KHCO3)^2",
["KCO2","KHCO3"]
)
SedTrace stores precomputed dissociation constants (on the free proton scale) on highresolution grids of salinity, temperature and pressure in /src/dissociation_constant.jld2
following the Guide to best practices for ocean CO_{2} measurements (Dickson et al., 2007) as implemented in the seacarb package (Gattuso et al., 2021). During code generation, SedTrace will compute the dissociation constants at the in situ salinity, temperature and pressure by interpolation. In the case of SimpleFe
, the dissociation constant of H_{2}O (KH2O
), the first (KCO2
) and second (KHCO3
) dissociation constants of H_{2}CO_{3}, and the first dissociate constant of H_{2}S (KH2S
) are computed as follows.
#————————————–
# Acid dissociation constants
#————————————–
KH2O = 7.9445878598109790E15 #H 1th dissociation constant
KCO2 = 8.3698755808176183E07 #TCO2 1th dissociation constant
KHCO3 = 4.6352156109843975E10 #TCO2 2th dissociation constant
KH2S = 1.2845618784784923E07 #H2S 1th dissociation constant
Given the concentrations of EIs and protons, SedTrace computes the concentrations of the individual species, as well as their transport and boundary condition terms following Sect. 3. SedTrace uses speciesspecific diffusion coefficients computed following Sect. 3.3. Many diagenetic models transport EIs and TA using the diffusion coefficients of the dominant species (for example, using the diffusion coefficient of HCO${}_{\mathrm{3}}^{}$ to represent all TCO_{2} species) or a fixed weighted average of the diffusion coefficients of the individual species (for example, computing a bulk diffusion coefficient for TCO_{2} using the weighted average of the diffusion coefficients of HCO${}_{\mathrm{3}}^{}$, CO${}_{\mathrm{3}}^{\mathrm{2}}$ and CO_{2} where the weight is fixed based on bottomwater speciation). However, studies have shown that these practices may lead to modeled pore water pH being different than using the speciesspecific diffusion coefficient (Luff et al., 2001; Jourabchi et al., 2005). In the SimpleFe
example, code for the transport and boundary conditions of HCO${}_{\mathrm{3}}^{}$ is as follows.
@.. HCO3 = H * KCO2 * TCO2 / (H^2 + H
* KCO2 + KCO2 * KHCO3)
mul!(HCO3_tran, AmHCO3, HCO3)
HCO3_tran[1] += BcAmHCO3[1] * HCO3[1]
+ BcCmHCO3[1]
HCO3_tran[Ngrid] += BcAmHCO3[2]
* HCO3[Ngrid] + BcCmHCO3[2]
@.. HCO3_tran += alpha * (HCO30  HCO3)
SedTrace then computes the transport and kinetic reaction terms of TA and EIs by summing over the species:
where ${\text{EI}}_{l}^{m}$ is the mth species of EI_{l}, TA^{n} is the nth species in the definition of TA and ζ_{n} is its coefficient in Eq. (26). ${\mathit{\nu}}_{{\mathrm{EI}}_{l}^{m}}^{\mathit{\eta},k}$ and ${\mathit{\nu}}_{{\mathrm{TA}}^{n}}^{\mathit{\eta},k}$ are the stoichiometric coefficients of these species in the kinetic reaction r^{η,k}, respectively.
The user can use the individual species when writing the chemical equations of the kinetic reaction rates. SedTrace will parse the equations and identify ${\mathit{\nu}}_{{\mathrm{EI}}_{l}^{m}}^{\mathit{\eta},k}$ and ${\mathit{\nu}}_{{\mathrm{TA}}^{n}}^{\mathit{\eta},k}$ automatically. For example, SedTrace recognizes that for the reaction RFeS_pre
in Table 2, the stoichiometric coefficient of TH_{2}S is ${\mathit{\nu}}_{{\mathrm{HS}}^{}}=\mathrm{1}$, and the stoichiometric coefficient of TA is ${\mathit{\nu}}_{{\mathrm{HS}}^{}}{\mathit{\nu}}_{{\mathrm{H}}^{+}}=\mathrm{2}$. For the SimpleFe
model the reactiontransport code for TA and EIs is as follows.
# Transport of EIs
@.. dTCO2 = HCO3_tran + CO3_tran + CO2_tran
@.. dTH2S = H2S_tran + HS_tran
# Transport of TA
@.. TA_tran = 1 * H_tran + 1 * OH_tran
@.. TA_tran += 1 * HCO3_tran + 2 * CO3_tran + 0 * CO2_tran
@.. TA_tran += 0 * H2S_tran + 1 * HS_tran
# Kinetic reaction rates of EIs
@.. S_TCO2 = 1 * RFeOOHPOC * dstopw + 1 * RSO4POC * dstopw
@.. S_TH2S =
1 / 2 * RSO4POC * dstopw + 1 * RFeOOHH2S * dstopw + 1 * RFeS_pre
# Kinetic reaction rates of TA
@.. S_TA =
8 * RFeOOHPOC * dstopw + 1 * RSO4POC * dstopw + 4 * RFeOOHH2S * dstopw + 1 * RFeS_pre
@.. S_TA += 1 * RFeS_pre
SedTrace does not explicitly model TA; rather, Eqs. (35) and (36) are substituted back into Eq. (25) to eliminate the TA terms to arrive at a diagenetic equation of [H^{+}]:
In the terminology of Hofmann et al. (2010), $\frac{\partial \text{TA}}{\partial \left[{\mathrm{H}}^{+}\right]}$ is the buffer factor and ${\mathit{\nu}}_{\mathrm{H}}^{\mathit{\eta},k}=(\sum _{n}{\mathit{\zeta}}_{n}{\mathit{\nu}}_{{\mathrm{TA}}^{n}}^{\mathit{\eta},k}\sum _{l}\frac{\partial \text{TA}}{\partial {\text{EI}}_{l}}\sum _{m}{\mathit{\nu}}_{{\mathrm{EI}}_{l}^{m}}^{\mathit{\eta},k})$ is the fractional stoichiometric coefficient of protons in the kinetic reaction r^{η,k}. Equations (37) to (39) show that the advantage of DSA is that the rate of pH change can be clearly partitioned at the level of individual species and reactions.
The user only needs to supply the boundary conditions for the EIs and pH. SedTrace will internally compute the boundary conditions for the individual species and form the transport terms in Eq. (37) according to Sect. 3. The user can use the EIs or the individual species when writing the chemical reactions. SedTrace will set the stoichiometry coefficients in Eq. (38) automatically according to Sect. 4.1. For the SimpleFe
model the code for Eqs. (37) to (39) is as follows.
# dTA/dEIs
@.. dTA_dTCO2 = KCO2 * (H + 2 * KHCO3)
/ (H^2 + H * KCO2 + KCO2 * KHCO3)
@.. dTA_dTH2S = KH2S / (H + KH2S)# dTA/dH
@.. dTA_dH = (H^2 + KH2O) / H^2
@.. dTA_dH +=
KCO2 * TCO2 * (H^2 + 4 * H * KHCO3
+ KCO2 * KHCO3) /
(H^2 + H * KCO2 + KCO2 * KHCO3)^2
@.. dTA_dH += KH2S * TH2S / (H + KH2S)^2# transport of individual species
mul!(HCO3_tran, AmHCO3, HCO3)
HCO3_tran[1] += BcAmHCO3[1] * HCO3[1]
+ BcCmHCO3[1]
HCO3_tran[Ngrid] += BcAmHCO3[2]
* HCO3[Ngrid] + BcCmHCO3[2]
@.. HCO3_tran += alpha * (HCO30  HCO3)# ... other species
# transport of EIs
@.. dTCO2 = HCO3_tran + CO3_tran
+ CO2_tran
@.. dTH2S = H2S_tran + HS_tran# transport of TA
@.. TA_tran = 1 * H_tran + 1 * OH_tran
@.. TA_tran += 1 * HCO3_tran + 2
* CO3_tran + 0 * CO2_tran
@.. TA_tran += 0 * H2S_tran + 1 * HS_tran# transport of proton
@.. dH = TA_tran
@.. dH = dTCO2 * dTA_dTCO2
@.. dH = dTH2S * dTA_dTH2S
@.. dH = dH / dTA_dH# kinetic reaction rates of proton
@.. S_H = S_TA
@.. S_H = S_TCO_2 * dTA_dTCO_2
@.. S_H = S_TH2S * dTA_dTH2S
@.. S_H = S_H / dTA_dH
The apparent dissociation constants depend not only on ionic strength but also composition (Zeebe and WolfGladrow, 2001). The assumption of the pH modeling laid out above is that the ionic medium of interest has a similar major ion composition as modern seawater; otherwise, the dissociation constants will not be applicable. Thus, SedTrace 1.0 is not suitable for modeling pore water with a major ion composition considerably different from modern seawater, which may happen in late diagenesis or in extreme or ancient environments.
4.3 Speciation modeling
Speciation of TEIs in diagenetic models is highly diverse and there is a lack of universally accepted practice. As such, SedTrace leaves the model decision to the user, while only requiring the user to format the input consistently. Thus, the responsibility of choosing a suitable set of aqueous and adsorbed species, equilibrium constants, and speciation equations lies with the user, whereas SedTrace is used to convert user input to code.
In the DSA of speciation modeling, SedTrace models the total concentration of a model substance (M^{T} in unit of per volume fluid), which is the sum of the total dissolved (M^{f} in units per volume of fluid) and total adsorbed (M^{s} in units per volume of solid) concentrations:
and M^{f}and M^{s} are themselves sums of the concentrations of individual dissolved and adsorbed species, respectively.
To indicate that speciation modeling is required for a model substance, the user needs to specify dissolved_speciation
in the type
column in the substances
sheet (Table 1) and provide the dissolved and adsorbed speciation information in the speciation
(Table 3) and adsorption
(Table 4) sheets. The name given in the substances
sheet is the code name of the total concentration M^{T}. Internally SedTrace will set the code names of the total dissolved M^{f} and total adsorbed M^{s} concentrations by appending the postfix _dis
and _ads
to the code name of the model substance. For example, the user names the total Fe TFe
in SimpleFe
, and SedTrace will name the total dissolved and adsorbed Fe TFe_dis
and TFe_ads
, respectively. The code names (column dissolved
) and chemical formulae (column formula
) of the individual dissolved species should be supplied to speciation
(Table 3). The code names (column adsorbed
) of the individual adsorbed species should be supplied to adsorption
(Table 4), but no chemical formula for the adsorbed species is needed.
The user can add the dissolved speciation reaction of a trace element M of the following format to the column equation
in speciation
(Pierrot and Millero, 2017):
which describes complexation with the pth dissolved ligand L_{p} to form aqueous species M(L_{p})_{q}, and q is the number of the ligand in the complex. Assuming local equilibrium, the concentration of the complexed species is
where ${K}_{M({L}_{p}{)}_{q}}$ is the apparent equilibrium constant supplied to column logK
, and [M] is the concentration of the base species (e.g., free ion) M. [M^{f}] is the sum of the concentrations of the base and complexed species:
In the speciation
sheet, the base species is indicated by writing a trivial speciation equation
such as Fe{2+} = Fe{2+}
with logK = 0
in the SimpleFe
model where free Fe^{2+} is the chosen base species of TFe
(Table 3). SedTrace will parse and check the chemical balance of the equations as discussed in Sect. 4.1.1. The dissolved ligands have to be modeled or specified by users. In SimpleFe
[HS^{−}] is computed as part of the pH model, and the concentration of Cl^{−} is not modeled but deemed a constant and supplied to parameters
. SedTrace 1.0 does not model the speciation of the ligands, and thus [L_{p}] in Eq. (42) refers to the total dissolved ligand concentrations. Therefore the apparent equilibrium constant supplied by the user needs to be precomputed using the total activity coefficients of model species based on ion pairing or Pitzer models (Pierrot and Millero, 2017; Millero and Schreiber, 1982; Zeebe and WolfGladrow, 2001). This is a limitation of SedTrace 1.0 which I plan to resolve in the future.
The formulation of adsorption in the literature of diagenetic modeling is also diverse (Boudreau, 1997; Berner, 1980; Wang and Van Cappellen, 1996; Katsev et al., 2006). The challenge lies in the fact that the surface properties of sedimentary particles are poorly known. SedTrace thus does not constrain the formulation of adsorption and lets the user specify the concentrations of the adsorbed species directly.
In the adsorption
sheet (Table 4), each row should list the adsorption of one dissolved species onto one surface. The user names the adsorbed species in the adsorbed
column, the dissolved species to be adsorbed in the dissolved
column and the surface to be adsorbed onto in the surface
column. All three columns should contain code names only. The user then supplies a mathematical expression f^{ads} to compute the adsorbed species concentration in the expression
column as a function of the concentrations of the dissolved species and surface:
where $[{M}^{\mathrm{s}}\equiv {S}_{\mathit{\kappa}}^{\mathit{\lambda}}]$ is the concentration of the λth dissolved species adsorbed onto the κth particle surface ≡S_{κ}. The dissolved species M^{dis} can be one of M, M(L_{p})_{q} or M^{f}.
The term “surface” is used heuristically here and can refer to any modeled solid substance. For example, in the SimpleFe
model (Table 4), Fe is adsorbed onto POC
, and the expression
for the adsorbed species is Fe_ads = KFe_ads*POC*TFe_dis
, where the adsorption is assumed to be indifferent to the dissolved speciation, and thus the total dissolved Fe concentration TFe_dis
(M^{f}) is used. The adsorption constant KFe_ads
here is an apparent constant specific to the user's model decision and needs to be provided to the parameters
sheet. The concentration of adsorbed species could sometimes be independent of any surface: for example, in the classic linear isothermal used by many diagenetic models (Berner, 1980) such that Fe_ads = KFe_ads*TFe_dis
. In this case the surface
column should be left empty.
SedTrace will group the adsorbed species by surface and compute the total concentration adsorbed onto the surface ≡S_{κ} ($\left[{M}_{\mathit{\kappa}}^{\mathrm{s}}\right]$, “empty surface” is a special surface) and the total concentration adsorbed onto all surfaces ([M^{s}]):
Internally SedTrace creates a code name for ${M}_{\mathit{\kappa}}^{\mathrm{s}}$ by appending the postfix ads_surface
to the substance name: for example, TFe_ads_POC
for the total Fe adsorbed onto POC. If the surface is empty, the code name for this adsorbed species is created by appending the postfix _ads_nsf
to the substance name: for example, TFe_ads_nsf
for the linear isothermal case.
To generate the speciation code, SedTrace solves Eqs. (40) to (46) using the solveset
function of solving systems of symbolic nonlinear equations from SymPy.jl (Meurer et al., 2017). To do so the usersupplied expression f^{ads} needs to be analytically invertible. The result is a set of symbolic expressions to compute [M], [M(L_{p})_{q}], [M^{f}], $[{M}^{\mathrm{s}}\equiv {S}_{\mathit{\kappa}}^{\mathit{\lambda}}$], $\left[{M}_{\mathit{\kappa}}^{\mathrm{s}}\right]$ and [M^{s}] using [M^{T}], [L_{p}], [≡S_{κ}] and the equilibrium constants, which are converted to Julia code.
The symbolic derivation starts by rewriting Eq. (44) to substitute [M] and [M(L_{p})_{q}] by [M^{f}] using Eqs. (42) and (43). Together with Eqs. (40), (45) and (46), SedTrace derives an expression to compute [M^{f}] using [M^{T}] and [≡S_{κ}]. In the case of TFe_dis
in the SimpleFe
model the code is as follows.
# Concentrations of total dissolved species
@.. TFe_dis = TFe / (KFe_ads * POC * dstopw + 1)
Substituting this expression back to Eqs. (42) and (43), SedTrace derives the expressions to compute [M] and [M(L_{p})_{n}] using [M^{f}] and [L_{p}]. The code for the individual dissolved Fe species is as follows.
# Concentrations of the individual dissolved species
@.. Fe_aq =
The user may not need to compute the concentrations of all species; for example, in this example only the concentration of the free species
3.98107170553497e6 * TFe_dis / (
0.01778279410038921 * CO3
+ 3.019951720402014e6 * Cl + 1.0
* HS +
3.630780547701011e5 * SO4
+ 3.98107170553497e6
)
@.. FeCl_aq =
3.019951720402014e6 * Cl * TFe_dis / (
0.01778279410038921 * CO3
+ 3.019951720402014e6 * Cl + 1.0
* HS +
3.630780547701011e5 * SO4
+ 3.98107170553497e6
)
@.. FeSO4_aq =
3.630780547701011e5 * SO4 * TFe_dis / (
0.01778279410038921 * CO3
+ 3.019951720402014e6 * Cl + 1.0
* HS
+
3.630780547701011e5 * SO4
+ 3.98107170553497e6
)
@.. FeCO3_aq =
0.01778279410038921 * CO3 * TFe_dis / (
0.01778279410038921 * CO3
+ 3.019951720402014e6 * Cl + 1.0
* HS
+
3.630780547701011e5 * SO4
+ 3.98107170553497e6
)
@.. FeHS_aq =
1.0 * HS * TFe_dis / (
0.01778279410038921 * CO3
+ 3.019951720402014e6 * Cl + 1.0
* HS
+
3.630780547701011e5 * SO4
+ 3.98107170553497e6
)Fe_aq
is needed. To choose the species whose concentrations are required by the model, set the value to 1 in the code
column in the speciation
sheet.
SedTrace then computes $\left[{M}^{\mathrm{s}}\equiv {S}_{\mathit{\kappa}}^{\mathit{\lambda}}\right]$, $\left[{M}_{\mathit{\kappa}}^{\mathrm{s}}\right]$ and [M^{s}] using [M], [M^{f}], [M(L_{p})_{n}] and [≡S_{κ}] based on Eqs. (44) to (46); for example, the code for the individual and total adsorbed Fe species is as follows.
# Concentrations of the individual adsorbed species
@.. Fe_ads = KFe_ads * POC * TFe_dis
# Concentrations of the total adsorbed species onto each surface
@.. TFe_ads_POC = Fe_ads
# Concentrations of the total adsorbed species
@.. TFe_ads = Fe_ads
In the next step, SedTrace computes the transport and reaction terms for M^{T} using M^{f} and M^{s}:
${T}_{{M}^{\mathrm{f}}}$ is the transport of total dissolved species summed over individual species. SedTrace uses the same diffusion coefficients for the dissolved species, as the diffusivities of the aqueous complexes are generally not known. ${T}_{{M}^{\mathrm{s}}}$ is the transport of the total adsorbed species summed over individual species, which are subject to the same transport mechanism like normal solid substances. The user needs to supply two sets of boundary conditions when modeling M^{T}. The boundary conditions of M^{f} are supplied to the substances
sheet, and the boundary conditions of M^{s} are supplied to the adsorption
sheet.
When writing the chemical equations and rate expressions of the kinetic reactions, the user can use the code name or chemical formula of any of the following: M, M(L_{p})_{q}, M^{f}, ${M}^{\mathrm{s}}\equiv {S}_{\mathit{\kappa}}^{\mathit{\lambda}}$, ${M}_{\mathit{\kappa}}^{\mathrm{s}}$, M^{s} or M^{T}. SedTrace will parse the reactions and set the stoichiometric coefficient for M^{T} (${\mathit{\nu}}_{{M}^{\mathrm{T}}}^{\mathit{\eta},k})$. The reactivetransport code for TFe
in the SimpleFe
model is as follows.
# Transport and boundary condition of total dissolved concentration
mul!(TFe_dis_tran, AmTFe_dis, TFe_dis)
TFe_dis_tran[1] += BcAmTFe_dis[1] * TFe_dis[1] + BcCmTFe_dis[1]
TFe_dis_tran[Ngrid] += BcAmTFe_dis[2] * TFe_dis[Ngrid] + BcCmTFe_dis[2]
# Transport and boundary condition of total adsorbed concentration
mul!(TFe_ads_tran, AmTFe_ads, TFe_ads)
TFe_ads_tran[1] += BcAmTFe_ads[1] * TFe_ads[1] + BcCmTFe_ads[1]
TFe_ads_tran[Ngrid] += BcAmTFe_ads[2] * TFe_ads[Ngrid] + BcCmTFe_ads[2]
# Transport of total concentration
@.. dTFe = TFe_dis_tran * 1 + TFe_ads_tran * dstopw
@.. dTFe += alpha * (TFe_dis0  TFe_dis)
# Reaction rate of total concentration
@.. S_TFe = 4 * RFeOOHPOC * dstopw + 2 * RFeOOHH2S * dstopw + 1 * RFeS_pre
4.4 Sediment age
The user may sometimes need to know the sediment age. Use cases may rise, for example, when using the continuum reactivity model for organic carbon remineralization (Boudreau et al., 2008):
where k_{POC} is the reaction rate constant of organic carbon, ζ is the parameter controlling the shape of the gamma distribution of reactivity, a_{0} is the initial age of organic carbon, and a_{t} is the duration of remineralization. Users may use sediment age in placement of a_{t} (Freitas et al., 2021).
The diagenetic equation for sediment age is (Meile and Van Cappellen, 2005)
The top boundary condition of sediment age can be specified as Robin or Dirichlet type. If specifying a zero incoming flux, the modeled age can be interpreted with respect to the incoming sediments, the age of which is zero. If specifying a zero age at the SWI, the modeled age is set to zero whenever sediments make contact with the SWI. The bottom boundary condition of sediment age is specified as ${\left(\frac{\partial \text{Age}}{\partial x}\right}_{x=L}=\frac{\mathrm{1}}{{\left({v}^{\mathrm{s}}\right}_{x=L}}$ to be consistent with the burial velocity. The use can add Age
to the substances
sheet with proper boundary conditions to enable sediment age modeling.
The code generated by SedTrace is assembled into five Julia files (Fig. 1): parm.jl
and parm.struct.jl
containing the model parameters, cache.jl
containing the model cache, reactran.jl
containing the reactivetransport code, and jactype.jl
containing the sparsity pattern of the Jacobian matrix.
5.1 Parameters
The usersupplied parameters in the parameters
sheet, and those computed internally by SedTrace, are included in parm.jl
. Parameters that are directly needed by the ODE function in Eq. (5) are further packed into a container, called ParamStruct
, within the module Param
inside parm.struct.jl
using the Julia package Parameters.jl (Werder, 2022a). For instance, the diffusion coefficients are not needed by Eq. (5) directly and are thus only included in parm.jl
. The transport matrix A, which incorporates the diffusion coefficients, is needed directly by Eq. (5) and is thus included in ParamStruct
. This container makes it easier to pass and modify parameters during model simulations. For the SimpleFe
model the code of the ParamStruct
is as follows.
@with_kw mutable struct ParamStruct{ T} # T is a Parametric Type
# only showing a few entries
TFeID::StepRange{ Int64,Int64} = TFeID # index of TFe
AmTFe_dis::Tridiagonal{ T,Vector{ T} } = AmTFe_dis # transport matrix of TFe
TFe_dis0::T = TFe_dis0 # top boundary condition of TFe
kFeSpre::T = kFeSpre # FeS precipitation rate constant
end
5.2 Cache
The code generation process creates many intermediate variables, which could cause repeated memory allocation during model simulation. Technically, they can be eliminated by substitution using Symbolic computing, but it would render the code unreadable. SedTrace preserves them for code clarity and preallocates memories for them when initiating the model so the memories can be reused.
During code generation, SedTrace keeps track of which variables are intermediate variables. SedTrace adds them to a container Reactran
within the module Cache
inside cache.jl
and preallocates their memories using the package PreallocationTools.jl (Rackauckas and Nie, 2017). The cache is compatible with the dual number type used by the package ForwardDiff.jl, so automatic differentiation can be applied to the code (Revels et al., 2016).
mutable struct Reactran{ T} # T is a Parametric Type
TFe_dis::PreallocationTools.DiffCache{ Array{ T,1} ,Array} T,1} }
FeCl_aq::PreallocationTools.DiffCache{ Array{ T,1} ,Array{ T,1} }
...... # other intermediate variables
end
5.3 ODE function
The reactran.jl
file contains a function function(f::Cache.Reactran)(dC, C, parms:: Param.ParamStruct, t)
, which includes the reactivetransport code to update $\frac{\mathrm{d}}{\mathrm{d}t}\mathit{C}$ (dC
) in Eq. (5) at time t in place, given the current model state vector C
and parameters parms
. This function is compatible with the ODE solvers from DifferentialEquations.jl (Rackauckas and Nie, 2017) and the automatic differentiation tools in ForwardDiff.jl (Revels et al., 2016).
This function is assembled in the following sequence: (1) unpack the parameters contained in parms using the package UnPack.jl (Werder, 2022b); (2) retrieve preallocated cache; (3) compute the transport and boundary condition terms of model substances that do not require speciation; (4) compute pH, EI speciation, and related transport and boundary condition terms; (5) compute the speciation of model substances that require speciation, as well as their transport and boundary condition terms; (6) compute the individual kinetic reaction rates and the summed rates for model substances. A code skeleton for SimpleFe
model is shown here.
function (f::Cache.Reactran)(dC, C, parms::Param.ParamStruct, t)
#
# Parameters
#
@unpack TFeID,
#...... other parameters
kFeSpre = parms
#
# Cache
#
TFe_dis = PreallocationTools.get_tmp (f.TFe_dis, C)
#...... other caches
#
# Model state
#
TFe = @view C[TFeID]
dTFe = @view dC[TFeID]
#...... other model substances
#
# Transport of solid and dissolved substances see Sect. 3.4
#
#
# pH code see Sect. 4.2
#
#
# Speciation code see Sect. 4.3
#
# Concentrations of adsorbed/dissolved species
# Transport of adsorbed/dissolved species
@.. dTFe = TFe_dis_tran * 1 + TFe_ads_tran * dstopw
#
# Reaction code see Sect. 4.1
#
# Individual reaction rates
# Summed rates for model substances
@.. S_TFe = 4 * RFeOOHPOC * dstopw + 2 * RFeOOHH2S * dstopw + 1 * RFeS_pre
# sum transport and reaction
@.. dTFe += S_TFe
return nothing
end
5.4 Jacobian pattern
The Jacobian of the righthand side of Eq. (5) is
which is needed to improve numerical performance, and knowing its sparsity pattern (i.e., which elements are nonzero) can accelerate numerical computation considerably (Rackauckas and Nie, 2017). The sparsity pattern of A+B is set by the discretization scheme, and in SedTrace it is a Julia Tridiagonal
matrix. The sparsity pattern of $\frac{\partial \mathit{S}}{\partial \mathit{C}}$ is modelspecific, and SedTrace detects it during code generation. Without pH and speciation modeling, $\frac{\partial \mathit{S}}{\partial \mathit{C}}$ can be treated as a Julia BandedMatrix
, the upper and lower bandwidths of which are equal to the number of model substances. However, pH and speciation modeling introduces additional coupling to the Jacobian and increases the bandwidths of $\frac{\partial \mathit{S}}{\partial \mathit{C}}$. Thus, it may be better to treat $\frac{\partial \mathit{S}}{\partial \mathit{C}}$ as a Julia SparseMatrixCSC
, especially when the size of the Jacobian is large, given that most of the elements on the diagonals of $\frac{\partial \mathit{S}}{\partial \mathit{C}}$ are zero.
Nonzero elements can be introduced by direct coupling through kinetic reactions. Such coupling happens within the same grid cell. When parsing the chemical equations in Sect. 4.1.1 SedTrace knows which reactions affect a model substance. And by further parsing the rate expressions of the reactions, SedTrace knows which model substances the reaction rates of these reactions depend on. The parsing result is saved in model_parsing_diagnostics.xlsx
. SedTrace then establishes a dependency relationship for all model substances and sets the corresponding elements of the Jacobian to nonzero. For example, in SimpleFe
the summed reaction rate of TFe
(S_TFe
, see code in Sect. 5.3) includes the reaction RFeOOHPOC
, the kinetic rate of which depends on FeOOH
, O2
and POC
. Thus, the [${i}_{\mathrm{TFe}}+(j\mathrm{1})M$, ${i}_{\mathrm{FeOOH}}+(j\mathrm{1})M$], [${i}_{\mathrm{TFe}}+(j\mathrm{1})M$, ${i}_{{\mathrm{O}}_{\mathrm{2}}}+(j\mathrm{1})M$] and [${i}_{\mathrm{TFe}}+(j\mathrm{1})M$, ${i}_{\mathrm{POC}}+(j\mathrm{1})M$] elements of the Jacobian are nonzero. The dependency relationship in the SimpleFe
model is as follows.
Row 
substance 
dependence 

String 
String 

1 
H 
FeOOH,POC,SO4,TH2S,TCO2,H,TFe 
2 
POC 
FeOOH,POC,SO4 
3 
FeOOH 
FeOOH,POC,TH2S 
4 
TCO2 
FeOOH,POC,SO4 
5 
TFe 
FeOOH,POC,TH2S,TCO2,H,TFe,SO4 
6 
SO4 
SO4,FeOOH,POC 
7 
TH2S 
SO4,FeOOH,POC,TH2S,TCO2,H,TFe 
8 
FeS 
TCO2,H,TH2S,TFe,POC,SO4 
Speciation further introduces transport coupling to the Jacobian that does not happen inside the same grid cell. In the SimpleFe
model, the adsorbed Fe concentration depends on the surface POC
. Thus, transport causes TFe
to depend not only on POC
inside the same grid cell, but also the two neighboring cells. Therefore, the [${i}_{\mathrm{TFe}}+(j\mathrm{1})M$, ${i}_{\mathrm{POC}}+(j\mathrm{1})M$], [${i}_{\mathrm{TFe}}+(j\mathrm{1})M$, ${i}_{\mathrm{POC}}+(j\mathrm{2})M$] and [${i}_{\mathrm{TFe}}+(j\mathrm{1})M$, i_{POC}+(j)M] elements of the Jacobian are nonzero. SedTrace keeps a record of such relationships during code generation.
Similar coupling outside the same grid cell is introduced by pH modeling. The proton concentration depends on the EIs not only inside the same grid cell, but also the two neighboring cells because of coupled transport. Therefore, in SimpleFe
the [${i}_{\mathrm{H}}+(j\mathrm{1})M$, ${i}_{{\mathrm{TCO}}_{\mathrm{2}}}+(j\mathrm{1})M$], [${i}_{\mathrm{H}}+(j\mathrm{1})M$, ${i}_{{\mathrm{TCO}}_{\mathrm{2}}}+(j\mathrm{2})M$] and [${i}_{\mathrm{H}}+(j\mathrm{1})M$, ${i}_{{\mathrm{TCO}}_{\mathrm{2}}}+\left(j\right)M$] elements of the Jacobian are nonzero. SedTrace knows such dependency relationships internally when generating the pH code.
SedTrace assembles the nonzero elements and produces the Jacobian sparsity pattern, which is saved in jactype.jl
.
Equation (5) is a system of ODEs that are coupled, nonlinear (in S) and generally speaking stiff. Since TEIs are sensitive to many biogeochemical processes, a model for TEI diagenesis likely needs to include a large reaction network, together with pH and speciation modeling, which introduces additional nonlinear coupling. And to capture the sharp chemical gradient close to the SWI, a fine grid is often needed. The number of equations can thus reach greater than 10^{4} as in the case studies presented below. Our experience shows that the backward differential formula (BDF), a family of implicit linear multistep methods of time stepping, is among the most efficient for solving large systems of stiff diagenetic equations. SedTrace uses the BDF method from the CVODE solver in the SUite of Nonlinear and DIfferential/ALgebraic equation Solvers (SUNDIALS) package (Hindmarsh et al., 2005; Gardner et al., 2022), which is written in C but made accessible to Julia by the Sundials.jl package as part of DifferentialEquations.jl (Rackauckas and Nie, 2017).
At the nth time step of integration a system of nonlinear equations resulting from time discretization of Eq. (5) needs to be solved by CVODE (Hindmarsh et al., 2005):
where f is the righthand side of Eq. (5), and γ_{n} and a_{n} are the coefficients set by CVODE. SedTrace solves Eq. (51) using the Newton–Krylov method, which is efficient for large sparse systems (Knoll and Keyes, 2004).
The mth Newton iteration step to update C^{n} involves solving a system of linear equations:
where $\mathrm{\Delta}{\mathit{C}}^{n}={\mathit{C}}^{n,m+\mathrm{1}}{\mathit{C}}^{n,m}$ is the increment, I is the identity matrix and $\mathbf{J}=\frac{\partial \mathit{f}}{\partial \mathit{C}}$ is the Jacobian of the righthand side of Eq. (5).
The Krylov space iterative method of solving Eq. (52) requires proper preconditioning to be numerically fast (Knoll and Keyes, 2004). SedTrace applies a right preconditioner by default:
SedTrace uses the incomplete LU factorization (ILU) of I−γ_{n}J as the preconditioner P. Currently two options are available: the ILU with zerolevel fill (ILU0) from the ILUZero.jl (Covalt, 2022) package and the Crout version of ILU with drop tolerance from the IncompleteLU.jl package (Stoppels, 2022). The advantage of ILU0 is memory efficiency. Since the sparsity pattern of J does not change during iteration, the sparsity pattern of P is also fixed. SedTrace can reuse the preallocated memory when updating P. In comparison, the ILU with drop tolerance uses more memory because during each iteration the sparsity pattern of P may change and new memory needs to be allocated, but it has the advantage that the resulting factorization is a better approximation of I−γ_{n}J than ILU0.
Creating the preconditioner requires computing J. SedTrace computes J using the forwardmode automatic differentiation tools from ForwardDiff.jl (Revels et al., 2016). This computation is accelerated using a matrix coloring algorithm (Gebremedhin et al., 2005) from SparseDiffTools.jl (Gowda et al., 2022), taking advantage of the knowledge of the sparsity pattern detected by SedTrace in Sect. 5.4. The preconditioned system in Eq. (53) is then solved using an iterative method: for example, the generalized minimal residual method (GMRES) from the CVODE package.
The user generates model code and performs a simulation in the main.jl
file. The example of SimpleFe
is shown here.
using SedTrace
# model configuration
modeldirectory = (@__DIR__) * "$\mathrm{\backslash}\mathrm{\backslash}$"
modelfile = "model_config.SimpleFe.xlsx"
modelname = "SimpleFe"
modelconfig = ModelConfig(modeldirectory, modelfile, modelname, AllowDiscontinuity = false)
# generate a parameter sheet template
@time generate_parameter_template (modelconfig)
# generate model code
@time generate_code(modelconfig)
# load model code files
IncludeFiles(modelconfig)
# initial values
C0 = Param.C0;
# initialize parameters
parm = Param.ParamStruct();
# initialize cache and ODE function
OdeFun = Cache.init(C0, parm.Ngrid);
# initialize Jacobian
JacPrototype = JacType(Param.IDdict);
# test if the Jacobian is correct
TestJacobian(JacPrototype, OdeFun, C0, parm)
# benchmark the ODE function performance
BenchmarkReactran(OdeFun, C0, parm)
# benchmark the Jacobian performance
BenchmarkJacobian(JacPrototype, OdeFun, C0, parm)
# benchmark the preconditioner performance
BenchmarkPreconditioner(JacPrototype, OdeFun, C0, parm,:ILU0)
# configure the solver
solverconfig = SolverConfig(:GMRES, :ILU0, 2)
# configure the solution
solutionconfig = SolutionConfig
(
C0, # initial values
(0.0, 3000.0), # time span
reltol = 1e6, # relative tolerance
abstol = 1e18, # absolute tolerance
saveat = 100.0, # save time steps
callback = TerminateSteadyState(1e16, 1e6),
# terminate when steady state is reached
);
# run the model
solution = modelrun(OdeFun, parm, JacPrototype, solverconfig, solutionconfig);
# generate output and plot
generate_output(modelconfig, solution, showplt = true, saveplt=true)
The workflow (Fig. 1) begins with configuring the model using modelconfig
and then supplying information of the model directory, the Excel file and model name. If the user prefers to use the logistic approximation of the Heaviside function discussed in Sect. 4.1.2, set AllowDiscontinuity = false
. If needed the user can call generate_parameter_template(modelconfig)
, which parses the substances
, reactions
, speciation
and adsorption
sheets to identify which parameters are needed by the model and outputs a template model_parameter_template.xlsx
to assist the creation of the parameters
sheet. Once the Excel model configuration file is created, code generation is done by calling generate_code(modelconfig)
, creating the five Julia files discussed in Sect. 5. The Julia files needs to be loaded into the workspace module Main
by calling IncludeFiles(modelconfig)
. The parameters are loaded into the module Param
, and the cache and ODE function are loaded into the module Cache
. During code generation, SedTrace collects the results of parsing the Excel sheets and creates a file called model_parsing_diagnostics.xlsx
, which can help the user to diagnose potential issues of code generation.
In the next step the user initializes the initial values, parameters, ODE function and Jacobian pattern. Internally SedTrace generates a set of initial values C0
that are constant with respect to depth based on usersupplied boundary conditions (Meysman et al., 2003). Users can also supply their own initial values, for example using the output from previous model runs. The parameters are initialized by parm = Param.ParamStruct()
. The ODE function is initialized using OdeFun = Cache.init(C0, Ngrid)
. The Jacobian sparsity pattern is initialized by JacPrototype = JacType(Param.IDdict)
, where IDdict
is a Julia dictionary that stores the indices of the model substances.
SedTrace also provides functions for code testing. The function TestJacobian()
computes the Jacobian assuming it is dense, which is timeconsuming but accurate. The result is then compared with the Jacobian computed using JacPrototype
. This serves as a check on code generation. BenchmarkJacobian()
, BenchmarkReactran()
and BenchmarkPreconditioner()
are used to benchmark the performance and memory allocations of the Jacobian, ODE and preconditioner functions, respectively.
The user configures the numerical solver using solverconfig = SolverConfig(:method, :preconditioner, prec_side)
, where :method
is the numerical method such as :GMRES
; :preconditioner
is the name of the preconditioner, by default :ILU0
; and prec_side
controls whether it is left (1) or right (2) preconditioning in Eq. (53). The numerical solution is configured using solutionconfig = SolutionConfig(C0, tspan, reltol, abstol; callback)
to set the initial values (C0
), the time span (tspan
), and the relative (reltol
) and absolute (abstol
) tolerances for numerical convergence. Any callback function compatible with DifferentialEquations.jl can be supplied too. For example, the user can use TerminateSteadyState(rtol, atol)
from the DiffEqCallbacks.jl (Rackauckas and Nie, 2017) library to terminate the simulation once steady state is reached given the relative and absolute tolerances of $\frac{\mathrm{d}}{\mathrm{d}t}\mathit{C}$ (rtol
and atol
).
Model simulation is performed by calling solution = modelrun(OdeFun, parm, JacPrototype, solverconfig ,solutionconfig)
. Internally SedTrace creates the Jacobian, solver and preconditioner functions and formats the ODE function to be compatible with DifferentialEquations.jl, which carries out the numerical solution.
Finally, outputs are created by calling generate_output(modelconfig, solution; site = nothing, showplt = true, saveplt = false)
. SedTrace will compute the output variables listed in the output
sheet (Table 7). New variables can be computed by supplying their mathematical expressions as functions of the model substances in the expression
column. For example, if the user wants to output pH on the free proton scale, then an expression log10(H)
should be supplied. SedTrace also converts the default units to units specified by the user in the unit_profile
column by multiplying the conversion factors in the conversion_profile
column. SedTrace can also compute the benthic fluxes of the output variables at the SWI. This is enabled by setting the flux_top
column to 1. Similar to the model profiles, unit conversion for the flux is done using the conversion_flux
and unit_flux
columns.
SedTrace then plots the profiles and the fluxes of the output variables. The user can supply the measured profiles of these variables in the data sheet and the measured fluxes in the flux_top_measured
column in the output sheet. SedTrace will plot the measurements with the model output. To do so the user needs to specify the site name in the site
column and supply this name to generate_output
. The name of the substance
and the unit
in the data
sheet must exactly match those in the output
sheet for SedTrace to match the model results with measurements. Measured data are supplied to the depth
and value
columns, with optional error
values that will be used to create error bars on the plots. SedTrace will save the output in model_output.xlsx
, containing internal output of modeled substances (sheet Substances
), reaction rates (ReacRates
), saturation state (Omega
), pHrelated species (pHspecies
), speciation results (Speciation
) and intermediate variables (IntermediateVar
) in the default SedTrace units. The userspecified output is reported in the OutputProfile
and OutputFlux
sheets with the custom units. If saveplt = true
SedTrace will save the plots in a plot directory inside modeldirectory
.
In the /examples
folder I provide case studies of the generation and running of diagenetic models at different levels of complexity. These examples include models with analytical solutions that are used to validate SedTrace's code generation and numerical methods. I also present advanced cases studies to demonstrate SedTrace's capacity for modeling pH, speciation, radiogenic and stable isotopes.
8.1 Models with analytical solutions
A few simple models of the diagenesis of carbon and nutrients that have analytical solutions are presented in examples/analytical
. These include the POC1G
model for POC remineralization using the 1G kinetics, the Ammonia
model for organic N remineralization and NH${}_{\mathrm{4}}^{+}$ adsorption, the SulfateRedcution
model for oxidation of POC by sulfate, and the Phosphate
model for organic P remineralization and authigenic phosphate precipitation. These examples come from Berner's classic textbook (Berner, 1980).
Here I discuss the model pHBB1991
in more detail. Boudreau (1991) created a diagenetic model with an analytical solution to explain the pH change across the mat of sulfideoxidizing bacteria Beggiatoa in sediments from Danish lagoons (Jørgensen and Revsbech, 1983). This model is now generated here using SedTrace. It includes one kinetic reaction, the oxidation of HS^{−} by O_{2}. The kinetic rate is ${k}_{\mathrm{OS}}{e}^{a{\left(x{x}_{\mathrm{0}}\right)}^{\mathrm{2}}}$, where k_{OS} is the rate constant, and x is depth. The reaction is assumed to happen close to the mat at x_{0}=0.005 cm where dissolved O_{2} disappears and H_{2}S starts to increase; a controls the sharpness of this interface. The model substances are dissolved O_{2} and H^{+} as well as the EIs TCO_{2}, TH_{2}S and TH_{3}BO_{3}. Their Dirichlet boundary conditions are specified at the top (−0.05 cm) and bottom (0.15 cm) of the model domain. Porosity is assumed to be constant and equal to 1, and thus no distinction is made between seawater above the SWI and the pore water below. Molecular diffusion is the only transport mechanism.
I tested the sensitivity of the model error (with respect to the analytical solution) to the model grid. I created the model with either a uniform or a nonuniform grid. The gridtran
for the uniform grid is simply x>x0.05
, where the 0.05 cm offset is to take into account the fact that SedTrace's internally uniform grid starts from 0 cm. The gridtran
for the nonuniform grid is constructed using hyperbolic functions (Hoffmann and Chiang, 2000):
where L=0.2 cm is the length of the model domain and 0.05 cm is the offset. The resulting grid points are concentrated near x_{0}=0.005 cm, the degree of which is controlled by b. I tested both types of grid with Ngrid
of 100 and 200.
The model results, the analytical solution and the observations made using profiling microelectrodes are shown in Fig. 3. The nonuniform grid captures the sharp biogeochemical gradient near x_{0}=0.005 cm better than the uniform grid when Ngrid
is the same. The L^{∞} norms of the model errors of O_{2} are 0.086, 0.021, 0.015 and 0.0039 µM for the uniform grid (Ngrid=100
), uniform grid 2X (Ngrid=200
), nonuniform grid (Ngrid=100
) and uniform grid 2X (Ngrid=200
), respectively. Thus, using a suitable grid can considerably improve the model accuracy.
8.2 Oregon margin diagenetic Nd cycle
Neodymium is one of the rareearth elements (REEs), which are important tracers in chemical oceanography (Elderfield and Greaves, 1982). Its radiogenic isotope composition, expressed as ${\mathit{\epsilon}}_{\mathrm{Nd}}=(\frac{{}^{\mathrm{143}}\mathrm{Nd}/{}^{\mathrm{144}}{\mathrm{Nd}}_{\mathrm{Sample}}}{{}^{\mathrm{143}}\mathrm{Nd}/{}^{\mathrm{144}}{\mathrm{Nd}}_{\mathrm{CHUR}}}\mathrm{1})\times {\mathrm{10}}^{\mathrm{4}}$ where CHUR is the chondritic uniform reservoir, has been used to study modern and past ocean circulation, as well as marine and continental weathering (Goldstein and Hemming, 2003; Haley et al., 2017; Lacan and Jeandel, 2005; Frank, 2002). It is designated as a “key parameter” by the GEOTRACES program so that it needs to be measured by all affiliated cruises. One of the greatest challenges facing the study of the modern ocean Nd cycle is that its sedimentary cycle is poorly constrained. Recently studies suggest that a benthic flux from marine sediments, particularly in the deep sea, is likely the dominant source of seawater Nd, far exceeding that of riverine and dust inputs at the ocean surface (Abbott et al., 2015b; Du et al., 2018, 2020; Haley et al., 2017). Consequently, pore water ε_{Nd}, subject to diagenetic processes such as marine silicate weathering, can affect the water column ε_{Nd} (Abbott et al., 2015a, 2016; Du et al., 2016). Because of the low Nd concentration, ∼1 L of pore water is typically required to make an isotope measurement. Thus, pore water ε_{Nd} analysis to date has only been done at three sites on the Oregon margin in the northeastern Pacific (Abbott et al., 2016, 2015a). Modeling is needed to make efficient use of the data and extrapolate to other regions once the diagenetic processes are well understood.
Our group recently published a reactivetransport model for the early diagenesis of Nd at a deepsea site (3000 m water depth) on the Oregon margin (Du et al., 2022), which is regenerated here using SedTrace. This model has 41 kinetic reactions, including the organic matter remineralization sequence, secondary redox reactions, and the diagenesis of carbonate, sulfide and opal. Pore water pH is modeled by including the EIs TCO_{2}, TH_{2}S and TH_{3}BO_{3}. Aqueous speciation of Fe^{2+} and Al^{3+} is included to model mineral dissolution and precipitation. Adsorption of Fe^{2+} and Mn^{2+} onto Fe and Mn oxides is also included. ^{144}Nd (Ndnr
) and the radiogenic ^{143}Nd (Ndr
) are modeled as two tracers. The following aqueous and solidphase speciation is included for both isotopes: complexation with CO${}_{\mathrm{3}}^{\mathrm{2}}$, HCO${}_{\mathrm{3}}^{}$, OH^{−}, Cl^{−}, SO${}_{\mathrm{4}}^{\mathrm{2}}$ and H_{3}SiO${}_{\mathrm{4}}^{}$; adsorption onto Fe and Mn oxides; incorporation into Fe and Mn oxides by coprecipitation; and formation of the authigenic phosphate mineral rhabdophane, released as trace constituents from primary silicates. By including cocycling with Fe, Mn and phosphate as well as release by marine silicate weathering coupled to reverse weathering, the model successfully simulated the pore water Nd concentration and ε_{Nd} data (Fig. 4). This model is relatively large, containing 10 200 equations of 34 model substances on 300 grids between 0 and 50 cm. Interested readers should consult Du et al. (2022) for details of the model description. This model can be found in /examples/OregonNd
.
8.3 Santa Barbara Basin sediment biogeochemistry, pH and Mo
Santa Barbara Basin (SBB) is one of the California borderland basins. Its seasonally anoxic condition leads to organicrich and laminated sediments (Reimers et al., 1996). SBB is among the most studied locations for sediment diagenesis and has perhaps one of the most complete pore water datasets to offer in the literature (Reimers et al., 1996). Highquality pH measurements by in situ profiling microelectrodes and the availability of various TEI data make it ideal for benchmarking diagenetic models (Meysman et al., 2003). Using SedTrace, I generated a diagenetic model for SBB that includes sediment biogeochemistry, pH and Mo cycling. This example is included in /examples/SBB
.
The biogeochemical reaction network includes the classic redox sequence of aerobic respiration, denitrification, Mn and Fe reduction, sulfate reduction, and methanogenesis. The model also includes secondary redox reactions and the authigenesis of carbonates, sulfide, opal and carbonate fluorapatite (CFA). The adsorption of NH${}_{\mathrm{4}}^{+}$, Fe^{2+} and Mn^{2+} is treated using the linear isothermal. To model pH, I include the following equilibrium invariants: TCO_{2}, TH_{2}S, THSO_{4}, TH_{3}BO_{3}, TH_{3}PO_{4} and THF. The biogeochemical model includes 26 substances and 36 kinetic reactions. I use a 500point nonuniform grid with finer spacing close to the SWI. The total number of equations in the biogeochemical model is thus 13 000. Figure 5 shows the modeled sediment biogeochemistry.
The top 50 cm of the sediment represents ∼150 years of sedimentation. The model successfully reproduced the measurements of pore water constituents, which respond much faster to perturbations of the seasonal cycle and other variabilities, while the solid sediment components show nonsteadystate behaviors (Fig. 5). During the sampling time bottom water was low (∼9 µM) but not anoxic, and the oxygen penetration depth was ∼1 cm, below which H_{2}S is immediately detectable. There was active formation of authigenic minerals (Reimers et al., 1996). Intense Fe cycling leads to high concentrations of FeS and FeS_{2}. Decreasing pore water Ca concentration is evidence of authigenic carbonate precipitation, and the model shows that Ca, Fe and Mn carbonates are likely formed. Decreasing pore water F concentration and a relatively low PO_{4} concentration are explained by CFA precipitation in the model (Jahnke, 1984; Reimers et al., 1996).
The model also captures the measured pH profile (Fig. 6). The pH measurements were reported on the seawater scale, and the model computes pH on this scale by summing the concentrations of free protons, HSO${}_{\mathrm{4}}^{}$ and HF. Using DSA (Hofmann et al., 2008), the model can easily partition the changes in pH to relevant transport and reaction terms (Fig. 6b–e). The overall increase in pH is driven by Fe reduction in the model (Fig. 6d), as suggested by previous studies (Reimers et al., 1996; Boudreau and Canfield, 1988). This increase in pH is responsible for the saturation and precipitation of authigenic carbonate and CFA. The slight decrease in pH at ∼3.5 cm (Fig. 6a) is modeled by FeS precipitation, which releases protons to pore water (Fig. 6d). The reaction rates of proton are balanced by the transport rate largely attributed to by total alkalinity, TCO_{2} and TH_{2}S (Fig. 6b). This example highlights the capacity of SedTrace for mechanistic studies of what controls pore water pH following the DSA approach in Sect. 4.2.
Molybdenum is sensitive to sedimentary redox conditions, and its stable isotope composition, expressed as ${\mathit{\delta}}^{\mathrm{98}}\text{Mo}=(\frac{{}^{\mathrm{98}}\mathrm{Mo}/{}^{\mathrm{95}}{\mathrm{Mo}}_{\mathrm{sample}}}{{}^{\mathrm{98}}\mathrm{Mo}/{}^{\mathrm{95}}{\mathrm{Mo}}_{\mathrm{standard}}}\mathrm{1})\times {\mathrm{10}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}+\mathrm{0.25},$ where NIST SRM3134 is the commonly used standard. Its δ^{98}Mo is 0.25 ‰ by convention, which makes it an important proxy to study past ocean deoxygenation (Kendall et al., 2017). SBB provides a useful analogy for an anoxic ocean, and modeling the sedimentary Mo cycle here may help understand how the δ^{98}Mo works as a redox proxy. Here I present a test model for Mo diagenesis in SBB to demonstrate the capability of SedTrace to model stable isotope fractionation, complementing the radiogenic isotope example above.
In this model, I consider five dissolved Mo species, MoO${}_{\mathrm{4}}^{\mathrm{2}}$ and four thiomolybdate species (Erickson and Helz, 2000):
K_{i} represents the apparent equilibrium constants.
I include ^{98}Mo (Moh
) and ^{95}Mo (Mol
) as two tracers. I assume that the equilibrium isotope fractionation is induced during thiolation:
where ${\mathit{\alpha}}_{i}^{\mathrm{98}/\mathrm{95}}$ represents the fractionation factors, which are 1.0014, 1.0028, 1.00455 and 1.0063 (or 1.40 ‰, 2.80 ‰, 4.54 ‰ and 6.28 ‰, respectively, in the δ notation) for i=1 to 4 respectively estimated by ab initio calculation (Tossell, 2005) and recalculated by Kendall et al. (2017). In SedTrace, I add Eq. (55) to the speciation
sheet, choosing MoO${}_{\mathrm{4}}^{\mathrm{2}}$ as the base species. Presently SedTrace does not provide special treatment of isotope fractionation, so the user needs to incorporate the fractionation factor in the parameters (e.g., K_{i}) before supplying them to SedTrace. I use the constants from Erickson and Helz (2000) as ${K}_{i}^{\mathrm{95}}$ and then multiply them by ${\mathit{\alpha}}_{i}^{\mathrm{98}/\mathrm{95}}$ to get ${K}_{i}^{\mathrm{98}}$.
I test the model sensitivity to the Mo removal mechanism. In Case 1, I assume all thiomolybdate species can be removed by scavenging:
In Case 2 I assume that only tetrathiomolybdate can be removed by scavenging:
In this test I do not consider kinetic isotope fractionation during removal and diffusion, and reactions in Eq. (55) are assumed to be fast enough to reach local equilibrium. In Eq. (57) I assume all species are removed at the same rate constant. I assume the bottomwater δ^{98}Mo to be the same as the seawater, which is globally uniform at 2.34 ‰ (Kendall et al., 2017). The only source of authigenic Mo accumulation is pore water Mo removal supported by diffusion of seawater into sediment. In the model I also supply a lithogenic Mo flux to account for the reported 2 ppm lithogenic Mo in sediments (Zheng et al., 2000). The lithogenic δ^{98}Mo is assumed to be the same as the upper continental crust (UCC) of ∼0.3 ‰ (Kendall et al., 2017).
In SedTrace I further add six Mo tracers to the biogeochemical model described above, including the two Mo isotopes in pore water, as well as lithogenic and authigenic fractions, increasing the total number of model equations to 16 000.
In the two model cases I vary k_{rm1} and k_{rm2} to fit the measured pore water Mo concentration profile (Fig. 7b) (Zheng et al., 2000). The model results of the two cases are roughly the same. As H_{2}S concentration increases, pore water Mo speciation is dominated by the thiomolybdate species (Fig. 7d). The aqueous speciation in the two cases is identical since it only depends on the concentration of H_{2}S. The modeled total solid (sum of authigenic and lithogenic) sediment Mo concentrations in the two cases are also similar, but both are much lower than measured (Fig. 7c) (Zheng et al., 2000). Modeled authigenic Mo enrichments, relative to the 2 ppm lithogenic background, do not start until below ∼5 cm. The model results suggest that the Mo diffusion and removal rate at the time of pore water sampling cannot explain the history of sediment Mo enrichment in the top 10 cm. Nonsteadystate accumulation and other Mo sources, such as Mo carried into sediments by settling particles, are needed to close the sedimentary budget (Zheng et al., 2000).
Modeled thiomolybdate species all have lighter δ^{98}Mo than MoO${}_{\mathrm{4}}^{\mathrm{2}}$ (Fig. 7e). The removal of thiomolybdate thus makes the residual pore water δ^{98}Mo heavier than seawater (Fig. 7f). In Case 1 where all thiomolybdate species are removed, the apparent fractionation between pore water and authigenic sediment (<1.4 ‰) is much smaller than in Case 2 where only tetrathiomolybdate is removed (as much as 6.5 ‰). In contrast, the δ^{98}Mo of total sediment has little sensitivity to the removal mechanism (Fig. 7f). This illustrates the underappreciated challenge of applying detrital corrections to TEI data: the most useful data to differentiate the removal mechanism in this test are the authigenic δ^{98}Mo at nearzero authigenic enrichment, which unfortunately will have the greatest uncertainty when applying detrital correction in reality (Ciscato et al., 2018).
As pore water Mo removal becomes quantitative below 25 cm, authigenic δ^{98}Mo approaches the seawater δ^{98}Mo. However, unlike in a closed system, modeled authigenic δ^{98}Mo never reaches the seawater value even at depths of quantitative removal. This is because in a reactivetransport system, the memory of the authigenic enrichment of light Mo isotope in the zone of partial removal will persist.
Future developments of SedTrace will improve the user interface, the speciation modeling capacity and parameter selection. Currently SedTrace has no specialized interface for isotopes, and the user needs to add the isotopes as individual tracers and supply parameters that already incorporate the stable isotope fractionation factor and radiogenic isotope source. A more userfriendly interface for modeling isotopes is planned. Also, the current implementation of adsorption and aqueous speciation can be further improved to allow more flexible and mechanistic formulations of speciation: for example, including the speciation of dissolved and surface ligands. At the moment, SedTrace only precomputes selected parameters such as diffusion coefficients and acid dissociation constants, while it requires the user to supply the rest. Integration with the Miami ion interaction model for seawater (Pierrot and Millero, 2017) and the Pitzerequationbased seawater speciation model being developed by the SCOR working group (Humphreys et al., 2022) would allow SedTrace to precompute more parameters and lessen the burden on the user.
Future developments will also aim to add more choices of numerical methods and improve numerical performance. Currently SedTrace only uses CVODE with the iterative method of solving linear systems (Hindmarsh et al., 2005; Gardner et al., 2022). Although I have found it to be more efficient when solving a large stiff system of equations, in other cases other ODE solvers or the direct method of solving linear systems may be preferable. Further release will enhance the integration of SedTrace with DifferentialEquations.jl (Rackauckas and Nie, 2017) and provide the user with more numerical choices. Similarly, more options of preconditioners can also be added. Also, SedTrace currently has no parallel computing capacity (aside from those done internally by Julia), and model simulation cannot take advantage of the multicore architecture of modern computers. Future development will focus on enabling SedTrace for parallel computing, for example, through integration with the PETSc.jl package, the Julia interface for the Portable, Extensible Toolkit for Scientific Computation system that uses the Messaging Passing Interface (MPI) for a scalable solution of differential equations (Balay et al., 2022). Such high performance is needed if SedTrace is to be used to solve globalscale problems.
Currently SedTrace's capacity is limited to the forward modeling of steadystate diagenesis. A special interface can be added to enable the user to add timedependent forcing in parameters or boundary conditions for fully transient modeling. Moreover, since SedTrace already enables automatic differentiation, capacities of adjoint sensitivity analysis, parameter estimation, inverse modeling and data assimilation can be added by integration with relevant existing Julia packages (Rackauckas and Nie, 2017; Dunning et al., 2017).
Code of SedTrace 1.0, along with the Excel sheets and Julia scripts used in the case studies, can be found online in the GitHub repository (https://github.com/JianghuiDu/SedTrace.jl, last access: 26 July 2023) and the Zenodo repository (https://doi.org/10.5281/zenodo.7225861) (Du, 2022).
No original data are generated by this study. Source data from the pore water and sediment observations used for model evaluation in the case studies can be found in their original references, as cited in the main text, and are also contained in the SedTrace/examples/ case study directory (https://doi.org/10.5281/zenodo.7225861) (Du, 2022).
The author has declared that there are no competing interests.
I thank Derek Vance for suggestions that improved the paper. I thank Clare Reimers for sharing the data from the Santa Barbara Basin. SedTrace was inspired by many published diagenetic models, including CANDI (Boudreau, 1996), STEADYSED (Wang and Van Cappellen, 1996), MEDIA (Meysman et al., 2003) and ReacTran (Soetaert and Meysman, 2012). I thank the topical editor Andrew Yool for handling the paper and two anonymous reviewers for providing helpful comments.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No 891489. This work was supported by an ETH Zurich Postdoctoral Fellowship 192 FEL32.
This paper was edited by Andrew Yool and reviewed by two anonymous referees.
Abbott, A. N., Haley, B. A., and McManus, J.: Bottoms up: Sedimentary control of the deep North Pacific Ocean's ϵ_{Nd} signature, Geology, 43, 1035–1035, https://doi.org/10.1130/G37114.1, 2015a.
Abbott, A. N., Haley, B. A., McManus, J., and Reimers, C. E.: The sedimentary flux of dissolved rare earth elements to the ocean, Geochim. Cosmochim. Ac., 154, 186–200, https://doi.org/10.1016/j.gca.2015.01.010, 2015b.
Abbott, A. N., Haley, B. A., and McManus, J.: The impact of sedimentary coatings on the diagenetic Nd flux, Earth Planet. Sc. Lett., 449, 217–227, https://doi.org/10.1016/j.epsl.2016.06.001, 2016.
Anbar, A. D. and Rouxel, O.: Metal Stable Isotopes in Paleoceanography, Annu. Rev. Earth Planet. Sci., 35, 717–746, https://doi.org/10.1146/annurev.earth.34.031405.125029, 2007.
Anderson, R. F.: GEOTRACES: Accelerating Research on the Marine Biogeochemical Cycles of Trace Elements and Their Isotopes, Annu. Rev. Mar. Sci., 12, 49–85, https://doi.org/10.1146/annurevmarine010318095123, 2020.
Archer, C. and Vance, D.: The isotopic signature of the global riverine molybdenum flux and anoxia in the ancient oceans, Nat. Geosci., 1, 597–600, https://doi.org/10.1038/ngeo282, 2008.
Archer, D. E., Morford, J. L., and Emerson, S. R.: A model of suboxic sedimentary diagenesis suitable for automatic tuning and gridded global domains, Global Biogeochem. Cycles, 16, 171–1721, https://doi.org/10.1029/2000GB001288, 2002.
Balay, S., Abhyankar, S., Adams, M. F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E., Dalcin, L., Dener, A., Eijkhout, V., Faibussowitsch, J., Gropp, W. D., Hapla, V., Isaac, T., Jolivet, P., Karpeev, D., Kaushik, D., Knepley, M. G., Kong, F., Kruger, S., May, D. A., McInnes, L. C., Mills, R. T., Mitchell, L., Munson, T., Roman, J. E., Rupp, K., Sanan, P., Sarich, J., Smith, B. F., Zampini, S., Zhang, H., Zhang, H., and Zhang, J.: PETSc/TAO Users Manual, technical report, Argonne National Laboratory, https://doi.org/10.2172/1968587, 2022.
Berner, R. A.: Early Diagenesis: A Theoretical Approach, Princeton University Press, 260 pp., ISBN 9780691082608, 1980.
Bezanson, J., Edelman, A., Karpinski, S., and Shah, V.: Julia: A Fresh Approach to Numerical Computing, SIAM Rev., 59, 65–98, https://doi.org/10.1137/141000671, 2017.
Boudreau, B. P.: Modelling the sulfideoxygen reaction and associated pH gradients in porewaters, Geochim. Cosmochim. Ac., 55, 145–159, https://doi.org/10.1016/00167037(91)90407V, 1991.
Boudreau, B. P.: A methodoflines code for carbon and nutrient diagenesis in aquatic sediments, Comput. Geosci., 22, 479–496, https://doi.org/10.1016/00983004(95)001158, 1996.
Boudreau, B. P.: Diagenetic models and their implementation: modelling transport and reactions in aquatic sediments, Springer, Berlin, New York, ISBN 3540611258, 1997.
Boudreau, B. P. and Canfield, D. E.: A provisional diagenetic model for pH in anoxic porewaters: Application to the FOAM Site, J. Marine Res., 46, 429–455, https://doi.org/10.1357/002224088785113603, 1988.
Boudreau, B. P. and Canfield, D. E.: A comparison of closed and opensystem models for porewater pH and calcitesaturation state, Geochim. Cosmochim. Ac., 57, 317–334, https://doi.org/10.1016/00167037(93)90434X, 1993.
Boudreau, B. P., Arnosti, C., Jørgensen, B. B., and Canfield, D. E.: Comment on “Physical Model for the Decay and Preservation of Marine Organic Carbon”, Science, 319, 1616, https://doi.org/10.1126/science.1148589, 2008.
Burdige, D. J.: The biogeochemistry of manganese and iron reduction in marine sediments, EarthSci. Rev., 35, 249–284, https://doi.org/10.1016/00128252(93)90040E, 1993.
Burdige, D. J.: Geochemistry of Marine Sediments, Princeton University Press, 646 pp., ISBN 9780691095066, 2006.
Cameron, V. and Vance, D.: Heavy nickel isotope compositions in rivers and the oceans, Geochim. Cosmochim. Ac., 128, 195–211, https://doi.org/10.1016/j.gca.2013.12.007, 2014.
Ciscato, E. R., Bontognali, T. R. R., and Vance, D.: Nickel and its isotopes in organicrich sediments: implications for oceanic budgets and a potential record of ancient seawater, Earth Planet. Sc. Lett., 494, 239–250, https://doi.org/10.1016/j.epsl.2018.04.061, 2018.
Covalt, M.: ILUZero.jl, GitHub [code], https://github.com/mcovalt/ILUZero.jl (last access: 26 July 2023), 2022.
Crusius, J. and Thomson, J.: Comparative behavior of authigenic Re, U, and Mo during reoxidation and subsequent longterm burial in marine sediments, Geochim. Cosmochim. Ac., 64, 2233–2242, https://doi.org/10.1016/S00167037(99)004330, 2000.
Dale, A. W., Nickelsen, L., Scholz, F., Hensen, C., Oschlies, A., and Wallmann, K.: A revised global estimate of dissolved iron fluxes from marine sediments, Global Biogeochem. Cycles, 29, 691–707, https://doi.org/10.1002/2014GB005017, 2015.
Dickson, A. G., Sabine, C. L., and Christian, J. R. (Eds.): Guide to best practices for ocean CO_{2} measurements, PICES Special Publication 3, 191 pp., https://doi.org/10.25607/OBP1342, 2007.
Du, J.: SedTrace.jl: a Julia package for generating and running reactivetransport models of marine sediment diagenesis, Zenodo [code], https://doi.org/10.5281/zenodo.7225861, 2022.
Du, J., Haley, B. A., and Mix, A. C.: Neodymium isotopes in authigenic phases, bottom waters and detrital sediments in the Gulf of Alaska and their implications for paleocirculation reconstruction, Geochim. Cosmochim. Ac., 193, 14–35, https://doi.org/10.1016/j.gca.2016.08.005, 2016.
Du, J., Haley, B. A., Mix, A. C., Walczak, M. H., and Praetorius, S. K.: Flushing of the deep Pacific Ocean and the deglacial rise of atmospheric CO_{2} concentrations, Nat. Geosci., 11, 749–755, https://doi.org/10.1038/s4156101802056, 2018.
Du, J., Haley, B. A., and Mix, A. C.: Evolution of the Global Overturning Circulation since the Last Glacial Maximum based on marine authigenic neodymium isotopes, Quaternary Sci. Rev., 241, 106396, https://doi.org/10.1016/j.quascirev.2020.106396, 2020.
Du, J., Haley, B. A., Mix, A. C., Abbott, A. N., McManus, J., and Vance, D.: Reactivetransport modeling of neodymium and its radiogenic isotope in deepsea sediments: The roles of authigenesis, marine silicate weathering and reverse weathering, Earth Planet. Sc. Lett., 596, 117792, https://doi.org/10.1016/j.epsl.2022.117792, 2022.
Dunning, I., Huchette, J., and Lubin, M.: JuMP: A Modeling Language for Mathematical Optimization, SIAM Rev., 59, 295–320, https://doi.org/10.1137/15M1020575, 2017.
Elderfield, H. and Greaves, M. J.: The rare earth elements in seawater, Nature, 296, 214–219, https://doi.org/10.1038/296214a0, 1982.
Elrod, V. A., Berelson, W. M., Coale, K. H., and Johnson, K. S.: The flux of iron from continental shelf sediments: A missing source for global budgets, Geophys. Res. Lett., 31, L12307, https://doi.org/10.1029/2004GL020216, 2004.
Erickson, B. E. and Helz, G. R.: Molybdenum(VI) speciation in sulfidic waters:: Stability and lability of thiomolybdates, Geochim. Cosmochim. Ac., 64, 1149–1158, https://doi.org/10.1016/S00167037(99)004238, 2000.
Fiadeiro, M. E. and Veronis, G.: On weightedmean schemes for the finitedifference approximation to the advectiondiffusion equation, Tellus, 29, 512–522, https://doi.org/10.3402/tellusa.v29i6.11385, 1977.
Frank, M.: Radiogenic isotopes: Tracers of past ocean circulation and erosional input, Rev. Geophys., 40, 1001, https://doi.org/10.1029/2000RG000094, 2002.
Freitas, F. S., Pika, P. A., Kasten, S., Jørgensen, B. B., Rassmann, J., Rabouille, C., Thomas, S., Sass, H., Pancost, R. D., and Arndt, S.: New insights into largescale trends of apparent organic matter reactivity in marine sediments and patterns of benthic carbon transformation, Biogeosciences, 18, 4651–4679, https://doi.org/10.5194/bg1846512021, 2021.
Gardner, D. J., Reynolds, D. R., Woodward, C. S., and Balos, C. J.: Enabling New Flexibility in the SUNDIALS Suite of Nonlinear and Differential/Algebraic Equation Solvers, ACM Trans. Math. Softw., 48, 31:1–31:24, https://doi.org/10.1145/3539801, 2022.
Gattuso, J.P., Orr, J., Epitalon, J.M., Baldry, K., Hoshijima, U., and Abulos: seacarb: seawater carbonate chemistry. R package version 3.3.0., Zenodo [code], https://doi.org/10.5281/zenodo.4600014, 2021.
Gebremedhin, A. H., Manne, F., and Pothen, A.: What Color Is Your Jacobian? Graph Coloring for Computing Derivatives, SIAM Rev., 47, 629–705, https://doi.org/10.1137/S0036144504444711, 2005.
Goldstein, S. L. and Hemming, S. R.: Longlived Isotopic Tracers in Oceanography, Paleoceanography, and Icesheet Dynamics, in: Treatise on Geochemistry, edited by: Holland, H. D. and Turekian, K. K., Pergamon, Oxford, 453–489, https://doi.org/10.1016/B9780080959757.006173, 2003.
Gowda, S., Ma, Y., Churavy, V., Edelman, A., and Rackauckas, C.: Sparsity Programming: Automated SparsityAware Optimizations in Differentiable Programming, Program Transformations for ML Workshop at NeurIPS 2019, 2022.
Haley, B. A., Du, J., Abbott, A. N., and McManus, J.: The Impact of Benthic Processes on Rare Earth Element and Neodymium Isotope Distributions in the Oceans, Front. Marine Sci., 4, 426, https://doi.org/10.3389/fmars.2017.00426, 2017.
Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., and Woodward, C. S.: SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers, ACM Trans. Math. Softw., 31, 363–396, https://doi.org/10.1145/1089014.1089020, 2005.
Hoffmann, K. A. and Chiang, S. T.: Computational Fluid Dynamics, Engineering Education System, 502 pp., ISBN 0962373109, 2000.
Hofmann, A. F., Meysman, F. J. R., Soetaert, K., and Middelburg, J. J.: A stepbystep procedure for pH model construction in aquatic systems, Biogeosciences, 5, 227–251, https://doi.org/10.5194/bg52272008, 2008.
Hofmann, A. F., Middelburg, J. J., Soetaert, K., and Meysman, F. J. R.: pH modelling in aquatic systems with timevariable acidbase dissociation constants applied to the turbid, tidal Scheldt estuary, Biogeosciences, 6, 1539–1561, https://doi.org/10.5194/bg615392009, 2009.
Hofmann, A. F., Middelburg, J. J., Soetaert, K., A.WolfGladrow, D., and Meysman, F. J. R.: Proton cycling, buffering, and reaction stoichiometry in natural waters, Marine Chem., 121, 246–255, https://doi.org/10.1016/j.marchem.2010.05.004, 2010.
Homoky, W. B., Weber, T., Berelson, W. M., Conway, T. M., Henderson, G. M., Hulten, M. van, Jeandel, C., Severmann, S., and Tagliabue, A.: Quantifying trace element and isotope fluxes at the ocean–sediment boundary: a review, Phil. T. Roy. Soc. A, 374, 20160246, https://doi.org/10.1098/rsta.2016.0246, 2016.
Horner, T. J., Little, S. H., Conway, T. M., Farmer, J. R., Hertzberg, J. E., Janssen, D. J., Lough, A. J. M., McKay, J. L., Tessin, A., Galer, S. J. G., Jaccard, S. L., Lacan, F., Paytan, A., Wuttig, K., and Members, G.P. B. P. W. G.: Bioactive Trace Metals and Their Isotopes as Paleoproductivity Proxies: An Assessment Using GEOTRACESEra Data, Global Biogeochem. Cycles, 35, e2020GB006814, https://doi.org/10.1029/2020GB006814, 2021.
Hülse, D., Arndt, S., Daines, S., Regnier, P., and Ridgwell, A.: OMENSED 1.0: a novel, numerically efficient organic matter sediment diagenesis module for coupling to Earth system models, Geosci. Model Dev., 11, 2649–2689, https://doi.org/10.5194/gmd1126492018, 2018.
Humphreys, M. P., Waters, J. F., Turner, D. R., Dickson, A. G., and Clegg, S. L.: Chemical speciation models based upon the Pitzer activity coefficient equations, including the propagation of uncertainties: Artificial seawater from 0 to 45 ^{∘}C, Marine Chem., 244, 104095, https://doi.org/10.1016/j.marchem.2022.104095, 2022.
Jahnke, R. A.: The synthesis and solubility of carbonate fluorapatite, Am. J. Sci., 284, 58–78, https://doi.org/10.2475/ajs.284.1.58, 1984.
Jeandel, C.: Overview of the mechanisms that could explain the “Boundary Exchange” at the land–ocean contact, Phil. Trans. R. Soc. A, 374, 20150287, https://doi.org/10.1098/rsta.2015.0287, 2016.
Jørgensen, B. B. and Revsbech, N. P.: Colorless Sulfur Bacteria, Beggiatoa spp. and Thiovulum spp., in O_{2} and H_{2}S Microgradients, Appl. Environ. Microbiol., 45, 1261–1270, https://doi.org/10.1128/aem.45.4.12611270.1983, 1983.
Jourabchi, P., Cappellen, P. V., and Regnier, P.: Quantitative interpretation of pH distributions in aquatic sediments: A reactiontransport modeling approach, Am. J. Sci., 305, 919–956, https://doi.org/10.2475/ajs.305.9.919, 2005.
Jourabchi, P., Meile, C., Pasion, L. R., and Van Cappellen, P.: Quantitative interpretation of pore water O_{2} and pH distributions in deepsea sediments, Geochim. Cosmochim. Ac., 72, 1350–1364, https://doi.org/10.1016/j.gca.2007.12.012, 2008.
Katsev, S., Sundby, B., and Mucci, A.: Modeling vertical excursions of the redox boundary in sediments: Application to deep basins of the Arctic Ocean, Limnol. Oceanogr., 51, 1581–1593, https://doi.org/10.4319/lo.2006.51.4.1581, 2006.
Kendall, B., Dahl, T. W., and Anbar, A. D.: Good Golly, Why Moly? THE STABLE ISOTOPE GEOCHEMISTRY OF MOLYBDENUM, in: NonTraditional Stable Isotopes, edited by: Teng, F.Z., Watkins, J., and Dauphas, N., De Gruyter, 683–732, https://doi.org/10.1515/9783110545630017, 2017.
Knoll, D. A. and Keyes, D. E.: Jacobianfree Newton–Krylov methods: a survey of approaches and applications, J. Comput. Phys., 193, 357–397, https://doi.org/10.1016/j.jcp.2003.08.010, 2004.
Lacan, F. and Jeandel, C.: Neodymium isotopes as a new tool for quantifying exchange fluxes at the continent–ocean interface, Earth Planet. Sc. Lett., 232, 245–257, https://doi.org/10.1016/j.epsl.2005.01.004, 2005.
Lam, P. J. and Anderson, R. F.: GEOTRACES: The Marine Biogeochemical Cycle of Trace Elements and Their Isotopes, Elements, 14, 377–378, https://doi.org/10.2138/gselements.14.6.377, 2018.
Lau, K. V., Lyons, T. W., and Maher, K.: Uranium reduction and isotopic fractionation in reducing sediments: Insights from reactive transport modeling, Geochim. Cosmochim. Ac., 287, 65–92, https://doi.org/10.1016/j.gca.2020.01.021, 2020.
Lemaitre, N., Du, J., de Souza, G. F., Archer, C., and Vance, D.: The essential bioactive role of nickel in the oceans: Evidence from nickel isotopes, Earth Planet. Sc. Lett., 584, 117513, https://doi.org/10.1016/j.epsl.2022.117513, 2022.
Li, Y.H. and Gregory, S.: Diffusion of ions in sea water and in deepsea sediments, Geochim. Cosmochim. Ac., 38, 703–714, https://doi.org/10.1016/00167037(74)901458, 1974.
Little, S. H., Vance, D., WalkerBrown, C., and Landing, W. M.: The oceanic mass balance of copper and zinc isotopes, investigated by analysis of their inputs, and outputs to ferromanganese oxide sediments, Geochim. Cosmochim. Ac., 125, 673–693, https://doi.org/10.1016/j.gca.2013.07.046, 2014.
Little, S. H., Archer, C., McManus, J., Najorka, J., Wegorzewski, A. V., and Vance, D.: Towards balancing the oceanic Ni budget, Earth Planet. Sc. Lett., 547, 116461, https://doi.org/10.1016/j.epsl.2020.116461, 2020.
Llorente, V. J., ten Thije Boonkkamp, J. H. M., Pascau, A., and Anthonissen, M. J. H.: Similarities and differences of two exponential schemes for convectiondiffusion problems: The FVCF and ENATE schemes, Appl. Math. Comput., 365, 124700, https://doi.org/10.1016/j.amc.2019.124700, 2020.
Luff, R., Haeckel, M., and Wallmann, K.: Robust and fast FORTRAN and MATLAB^{®} libraries to calculate pH distributions in marine systems, Comput. Geosci., 27, 157–169, https://doi.org/10.1016/S00983004(00)000972, 2001.
Lyons, T. W., Anbar, A. D., Severmann, S., Scott, C., and Gill, B. C.: Tracking Euxinia in the Ancient Ocean: A Multiproxy Perspective and Proterozoic Case Study, Annu. Rev. Earth Planet. Sci., 37, 507–534, https://doi.org/10.1146/annurev.earth.36.031207.124233, 2009.
Maher, K., Steefel, C. I., DePaolo, D. J., and Viani, B. E.: The mineral dissolution rate conundrum: Insights from reactive transport modeling of U isotopes and pore fluid chemistry in marine sediments, Geochim. Cosmochim. Ac., 70, 337–363, https://doi.org/10.1016/j.gca.2005.09.001, 2006.
Martin, J. H., Knauer, G. A., Karl, D. M., and Broenkow, W. W.: VERTEX: carbon cycling in the northeast Pacific, DeepSea Res. Pt. A, 34, 267–285, https://doi.org/10.1016/01980149(87)900860, 1987.
Meile, C. and Van Cappellen, P.: Particle age distributions and O_{2} exposure times: Timescales in bioturbated sediments, Global Biogeochem. Cycles, 19, GB3013, https://doi.org/10.1029/2004GB002371, 2005.
Meile, C., Berg, P., Van Cappellen, P., and Tuncay, K.: Solutespecific pore water irrigation: Implications for chemical cycling in early diagenesis, J. Mar. Res., 63, 601–621, https://doi.org/10.1357/0022240054307885, 2005.
Meurer, A., Smith, C. P., Paprocki, M., Èertík, O., Kirpichev, S. B., Rocklin, M., Kumar, Am., Ivanov, S., Moore, J. K., Singh, S., Rathnayake, T., Vig, S., Granger, B. E., Muller, R. P., Bonazzi, F., Gupta, H., Vats, S., Johansson, F., Pedregosa, F., Curry, M. J., Terrel, A. R., Rouèka, Š., Saboo, A., Fernando, I., Kulal, S., Cimrman, R., and Scopatz, A.: SymPy: symbolic computing in Python, PeerJ Comput. Sci., 3, e103, https://doi.org/10.7717/peerjcs.103, 2017.
Meysman, F. J. R., Middelburg, J. J., Herman, P. M. J., and Heip, C. H. R.: Reactive transport in surface sediments. II. Media: an objectoriented problemsolving environment for early diagenesis, Comput. Geosci., 29, 301–318, https://doi.org/10.1016/S00983004(03)000074, 2003.
Millero, F. J. and Schreiber, D. R.: Use of the ion pairing model to estimate activity coefficients of the ionic components of natural waters, Am. J. Sci., 282, 1508–1540, https://doi.org/10.2475/ajs.282.9.1508, 1982.
Morel, F. M. M., Lam, P. J., and Saito, M. A.: Trace Metal Substitution in Marine Phytoplankton, Annu. Rev. Earth Planet. Sci., 48, 491–517, https://doi.org/10.1146/annurevearth053018060108, 2020.
Palmer, M. R. and Elderfield, H.: Sr isotope composition of sea water over the past 75 Myr, Nature, 314, 526–528, https://doi.org/10.1038/314526a0, 1985.
Paraska, D. W., Hipsey, M. R., and Salmon, S. U.: Sediment diagenesis models: Review of approaches, challenges and opportunities, Environ. Model. Softw., 61, 297–325, https://doi.org/10.1016/j.envsoft.2014.05.011, 2014.
Pasquier, B., Primeau, F. W., and John, S. G.: AIBECS.jl: A tool for exploring global marine biogeochemical cycles, J. Open Source Softw., 7, 3814, https://doi.org/10.21105/joss.03814, 2022.
Piepgras, D. J., Wasserburg, G. J., and Dasch, E. J.: The isotopic composition of Nd in different ocean masses, Earth Planet. Sc. Lett., 45, 223–236, https://doi.org/10.1016/0012821X(79)901250, 1979.
Pierrot, D. and Millero, F. J.: The Speciation of Metals in Natural Waters, Aquat. Geochem., 23, 1–20, https://doi.org/10.1007/s1049801692924, 2017.
Rackauckas, C. and Nie, Q.: DifferentialEquations.jl – A Performant and FeatureRich Ecosystem for Solving Differential Equations in Julia, J. Open Res. Softw., 5, 15, https://doi.org/10.5334/jors.151, 2017.
Reimers, C. E., Ruttenberg, K. C., Canfield, D. E., Christiansen, M. B., and Martin, J. B.: Porewater pH and authigenic phases formed in the uppermost sediments of the Santa Barbara Basin, Geochim. Cosmochim. Ac., 60, 4037–4057, https://doi.org/10.1016/S00167037(96)002311, 1996.
Revels, J., Lubin, M., and Papamarkou, T.: ForwardMode Automatic Differentiation in Julia, arXiv [preprint], https://doi.org/10.48550/arXiv.1607.07892, 26 July 2016.
Schlitzer, R., Anderson, R. F., Dodas, E. M., Lohan, M., Geibert, W., Tagliabue, A., Bowie, A., Jeandel, C., Maldonado, M. T., Landing, W. M., Cockwell, D., Abadie, C., Abouchami, W., Achterberg, E. P., Agather, A., AguliarIslas, A., van Aken, H. M., Andersen, M., Archer, C., Auro, M., de Baar, H. J., Baars, O., Baker, A. R., Bakker, K., Basak, C., Baskaran, M., Bates, N. R., Bauch, D., van Beek, P., Behrens, M. K., Black, E., Bluhm, K., Bopp, L., Bouman, H., Bowman, K., Bown, J., Boyd, P., Boye, M., Boyle, E. A., Branellec, P., Bridgestock, L., Brissebrat, G., Browning, T., Bruland, K. W., Brumsack, H.J., Brzezinski, M., Buck, C. S., Buck, K. N., Buesseler, K., Bull, A., Butler, E., Cai, P., Mor, P. C., Cardinal, D., Carlson, C., Carrasco, G., Casacuberta, N., Casciotti, K. L., Castrillejo, M., Chamizo, E., Chance, R., Charette, M. A., Chaves, J. E., Cheng, H., Chever, F., Christl, M., Church, T. M., Closset, I., Colman, A., Conway, T. M., Cossa, D., Croot, P., Cullen, J. T., Cutter, G. A., Daniels, C., Dehairs, F., Deng, F., Dieu, H. T., Duggan, B., Dulaquais, G., Dumousseaud, C., EchegoyenSanz, Y., Edwards, R. L., Ellwood, M., Fahrbach, E., Fitzsimmons, J. N., Russell Flegal, A., Fleisher, M. Q., van de Flierdt, T., Frank, M., Friedrich, J., Fripiat, F., Fröllje, H., Galer, S. J. G., Gamo, T., Ganeshram, R. S., GarciaOrellana, J., GarciaSolsona, E., GaultRingold, M., George, E., Gerringa, L. J. A., Gilbert, M., Godoy, J. M., Goldstein, S. L., Gonzalez, S. R., Grissom, K., Hammerschmidt, C., Hartman, A., Hassler, C. S., Hathorne, E. C., Hatta, M., Hawco, N., Hayes, C. T., Heimbürger, L.E., Helgoe, J., Heller, M., Henderson, G. M., Henderson, P. B., van Heuven, S., Ho, P., Horner, T. J., Hsieh, Y.T., Huang, K.F., Humphreys, M. P., Isshiki, K., Jacquot, J. E., Janssen, D. J., Jenkins, W. J., John, S., Jones, E. M., Jones, J. L., Kadko, D. C., Kayser, R., Kenna, T. C., Khondoker, R., Kim, T., Kipp, L., Klar, J. K., Klunder, M., Kretschmer, S., Kumamoto, Y., Laan, P., Labatut, M., Lacan, F., Lam, P. J., Lambelet, M., Lamborg, C. H., Le Moigne, F. A. C., Le Roy, E., Lechtenfeld, O. J., Lee, J.M., Lherminier, P., Little, S., LópezLora, M., Lu, Y., Masque, P., Mawji, E., Mcclain, C. R., Measures, C., Mehic, S., Barraqueta, J.L. M., van der Merwe, P., Middag, R., Mieruch, S., Milne, A., Minami, T., Moffett, J. W., Moncoiffe, G., Moore, W. S., Morris, P. J., Morton, P. L., Nakaguchi, Y., Nakayama, N., Niedermiller, J., Nishioka, J., Nishiuchi, A., Noble, A., Obata, H., Ober, S., Ohnemus, D. C., van Ooijen, J., O'Sullivan, J., Owens, S., Pahnke, K., Paul, M., Pavia, F., Pena, L. D., Peters, B., Planchon, F., Planquette, H., Pradoux, C., Puigcorbé, V., Quay, P., Queroue, F., Radic, A., Rauschenberg, S., Rehkämper, M., Rember, R., Remenyi, T., Resing, J. A., Rickli, J., Rigaud, S., Rijkenberg, M. J. A., Rintoul, S., Robinson, L.F., RocaMartí, M., Rodellas, V., Roeske, T., Rolison, J. M., Rosenberg, M., Roshan, S., Rutgers van der Loeff, M. M., Ryabenko, E., Saito, M. A., Salt, L. A., Sanial, V., Sarthou, G., Schallenberg, C., Schauer, U., Scher, H., Schlosser, C., Schnetger, B., Scott, P., Sedwick, P. N., Semiletov, I., Shelley, R., Sherrell, R. M., Shiller, A. M., Sigman, D. M., Singh, S. K., Slagter, H. A., Slater, E., Smethie, W. M., Snaith, H., Sohrin, Y., Sohst, B., Sonke, J. E., Speich, S., Steinfeldt, R., Stewart, G., Stichel, T., Stirling, C. H., Stutsman, J., Swarr, G. J., Swift, J. H., Thomas, A., Thorne, K., Till, C. P., Till, R., Townsend, A. T., Townsend, E., Tuerena, R., Twining, B. S., Vance, D., Velazquez, S., Venchiarutti, C., VillaAlfageme, M., Vivancos, S. M., Voelker, A. H. L., Wake, B., Warner, M. J., Watson, R., van Weerlee, E., Alexandra Weigand, M., Weinstein, Y., Weiss, D., Wisotzki, A., Woodward, E. M. S., Wu, J., Wu, Y., Wuttig, K., Wyatt, N., Xiang, Y., Xie, R. C., Xue, Z., Yoshikawa, H., Zhang, J., Zhang, P., Zhao, Y., Zheng, L., Zheng, X.Y., Zieringer, M., Zimmer, L. A., Ziveri, P., Zunino, P., and Zurbrick, C.: The GEOTRACES Intermediate Data Product 2017, Chem. Geol., 493, 210–223, https://doi.org/10.1016/j.chemgeo.2018.05.040, 2018.
SCOR Working Group: GEOTRACES – An international study of the global marine biogeochemical cycles of trace elements and their isotopes, Geochemistry, 67, 85–131, https://doi.org/10.1016/j.chemer.2007.02.001, 2007.
Soetaert, K. and Meysman, F.: Reactive transport in aquatic ecosystems: Rapid model prototyping in the open source software R, Environ. Model. Softw., 32, 49–60, https://doi.org/10.1016/j.envsoft.2011.08.011, 2012.
Soetaert, K., Herman, P. M. J., and Middelburg, J. J.: A model of early diagenetic processes from the shelf to abyssal depths, Geochim. Cosmochim. Ac., 60, 1019–1040, https://doi.org/10.1016/00167037(96)000130, 1996.
Sridhar, A., Tissaoui, Y., Marras, S., Shen, Z., Kawczynski, C., Byrne, S., Pamnany, K., Waruszewski, M., Gibson, T. H., Kozdon, J. E., Churavy, V., Wilcox, L. C., Giraldo, F. X., and Schneider, T.: Largeeddy simulations with ClimateMachine v0.2.0: a new opensource code for atmospheric simulations on GPUs and CPUs, Geosci. Model Dev., 15, 6259–6284, https://doi.org/10.5194/gmd1562592022, 2022.
Stoppels, H.: ILU for SparseMatrixCSC, GitHub [code], https://github.com/haampie/IncompleteLU.jl (last access: 26 July 2023), 2022.
Sulpis, O., Humphreys, M. P., Wilhelmus, M. M., Carroll, D., Berelson, W. M., Menemenlis, D., Middelburg, J. J., and Adkins, J. F.: RADIv1: a nonsteadystate early diagenetic model for ocean sediments in Julia and MATLAB/GNU Octave, Geosci. Model Dev., 15, 2105–2131, https://doi.org/10.5194/gmd1521052022, 2022.
Tagliabue, A., Bowie, A. R., Boyd, P. W., Buck, K. N., Johnson, K. S., and Saito, M. A.: The integral role of iron in ocean biogeochemistry, Nature, 543, 51–59, https://doi.org/10.1038/nature21058, 2017.
ten Thije Boonkkamp, J. H. M. and Anthonissen, M. J. H.: Extension of the Complete Flux Scheme to TimeDependent Conservation Laws, in: Numerical Mathematics and Advanced Applications 2009, Berlin, Heidelberg, 865–873, https://doi.org/10.1007/9783642117954_93, 2010.
ten Thije Boonkkamp, J. H. M. and Anthonissen, M. J. H.: The Finite VolumeComplete Flux Scheme for AdvectionDiffusionReaction Equations, J. Sci. Comput., 46, 47–70, https://doi.org/10.1007/s1091501093888, 2011.
Tossell, J. A.: Calculating the partitioning of the isotopes of Mo between oxidic and sulfidic species in aqueous solution, Geochim. Cosmochim. Ac., 69, 2981–2993, https://doi.org/10.1016/j.gca.2005.01.016, 2005.
Vance, D., Archer, C., Little, S. H., Köbberich, M., and de Souza, G. F.: The oceanic cycles of the transition metals and their isotopes, Ac. Geochim., 36, 359–362, https://doi.org/10.1007/s1163101701626, 2017.
Wang, Y. and Van Cappellen, P.: A multicomponent reactive transport model of early diagenesis: Application to redox cycling in coastal marine sediments, Geochim. Cosmochim. Ac., 60, 2993–3014, https://doi.org/10.1016/00167037(96)001408, 1996.
Werder, M.: Parameters.jl, GitHub [code], https://github.com/mauro3/Parameters.jl (last access: 26 July 2023), 2022a.
Werder, M.: UnPack.jl, GitHub [code], https://github.com/mauro3/UnPack.jl (last access: 26 July 2023), 2022b.
Zeebe, R. E. and WolfGladrow, D. A.: CO_{2} in Seawater: Equilibrium, Kinetics, Isotopes, Elsevier, Amsterdam, 346 pp., ISBN 9780444509468, 2001.
Zheng, Y., Anderson, R. F., van Geen, A., and Kuwabara, J.: Authigenic molybdenum formation in marine sediments: a link to pore water sulfide in the Santa Barbara Basin, Geochim. Cosmochim. Ac., 64, 4165–4178, https://doi.org/10.1016/S00167037(00)004956, 2000.
 Abstract
 Introduction
 Model framework
 Grid, transport and boundary conditions
 Biogeochemical reactions
 Julia code files
 Numerical solver
 Model simulation and output
 Case studies
 Future developments
 Code and data availability
 Competing interests
 Acknowledgements
 Disclaimer
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Model framework
 Grid, transport and boundary conditions
 Biogeochemical reactions
 Julia code files
 Numerical solver
 Model simulation and output
 Case studies
 Future developments
 Code and data availability
 Competing interests
 Acknowledgements
 Disclaimer
 Financial support
 Review statement
 References