RADIv1: a non-steady-state early diagenetic model for ocean sediments in Julia and MATLAB/GNU Octave

We introduce a time-dependent, one-dimensional model of early diagenesis that we term RADI, an acronym accounting for the main processes included in the model: chemical Reactions, Advection, molecular and bio-Diffusion, and 15 bio-Irrigation. RADI is targeted for study of deep-sea sediments, in particular those containing calcium carbonates (CaCO3). RADI combines CaCO3 dissolution driven by organic matter degradation with a diffusive boundary layer and integrates stateof-the-art parameterizations of CaCO3 dissolution kinetics in seawater, thus serving as a link between mechanistic surfacereaction modelling and global-scale biogeochemical models. RADI also includes CaCO3 precipitation, providing a continuum between CaCO3 dissolution and precipitation. RADI integrates components rather than individual chemical species for 20 accessibility and is straightforward to compare against measurements. RADI is the first diagenetic model implemented in Julia, a high-performance programming language that is free and open source, and it is also available in MATLAB/GNU Octave. Here, we first describe the scientific background behind RADI and its implementations. Then, we evaluate its performance in three selected locations and explore other potential applications, such as the influence of tides and seasonality on early diagenesis in the deep ocean. RADI is a powerful tool to study the time-transient and steady-state response of the sedimentary 25 system to environmental perturbation, such as deep-sea mining, deoxygenation or acidification events.


Introduction
The seafloor, which covers ~70% of the surface of the planet and modulates the transfer of materials and energy from the biosphere to the geosphere, remains for the vast majority unexplored. Today, this rich, poorly understood ecosystem is threatened locally by deep-sea mining activities (e.g. ploughing of the seabed), because it contains abundant valuable minerals 35 and metals essential for the energy transition (Thompson et al., 2018). The deep ocean is also being perturbed globally by climate change, including seawater acidification caused by the uptake of ~10 billion tons of anthropogenic carbon dioxide (CO2) into the ocean each year (Perez et al., 2018;Gruber et al., 2019), roughly a quarter of our total annual emissions (Friedlingstein et al., 2020). In this context, it is important to improve our understanding of the seafloor's response to environmental change. 40 Accumulation of sinking biogenic aggregates and lithogenic particles at the seafloor provides reactive material that regulates the chemical composition of sediment porewaters. Whereas biogenic particles typically sink through the water column at rates from a few meters to hundreds of meters per day (Riley et al., 2012), the same particles accumulate in sediments much more slowly, typically a few centimeters per thousand years (Jahnke, 1996). The residence time of solid particles in the top centimeter of sediments is therefore very long (a few hundred or thousand years) compared to their residence time in the 45 water column (a few weeks). Additionally, while solutes are dispersed by advection in the water column, molecular diffusion dominates in porewaters, which is slower. The long residence time of reactive solid material in surface sediments, coupled with the slow diffusive transport of dissolved species, can lead to large gradients in chemical composition between sediment porewaters and the overlying seawater, inducing solute fluxes between the two (Hammond et al., 1996). Thus, the top few millimeters of the seafloor play a significant role in many major marine biogeochemical cycles. 50 The overall rate of biogeochemical reactions is determined by the slowest, 'rate-limiting' step, which can be (i) transport to or from the reaction site or (ii) the reaction kinetics of the particle at the mineral-water interface. At the seafloor, the rate-limiting step for many biogeochemical reactions is solute transport via molecular diffusion through the sediment porewaters or through the diffusive boundary layer (DBL). The DBL is a thin film of water extending up to a few millimeters above the sediment-water interface in which molecular diffusion is the dominant mode of solute transport. The presence of a 55 DBL above the sediment-water interface ( Fig. 1) has been reported by several investigators (Morse, 1974;Archer et al., 1989b;Gundersen and Jørgensen, 1990;Santschi et al., 1991;Glud et al., 1994) and its thickness depends on the composition and roughness of the substrate, as well as on the flow speed of the overlying seawater (Chriss and Caldwell, 1982;Dade, 1993;Røy et al., 2002;Han et al., 2018). Diffusive fluxes of solutes across the sediment-water interface are driven by concentration gradients between the overlying seawater and the sediment column being considered. If most of the concentration gradient for 60 a given solute occurs within the porewaters, rather than within the DBL, then the diffusive flux of this solute is termed 'internal' or 'sediment-side controlled' (Boudreau and Guinasso, 1982). Conversely, if the majority of the concentration gradient for a given solute is within the DBL, the chemical flux across the sediment-water interface is termed 'external' or 'water-side transport-controlled'. In practice, the chemical exchange of most solutes is controlled by a combination of both regimes termed RADI uses the same set of reactive-transport partial differential equations as implemented in CANDI (Boudreau, 130 1996b), i.e., for each solute component, and for each solid component, where v is the concentration of a given component, t is time, φ is sediment porosity and φs is the solid-volume fraction, d is the 135 effective molecular diffusion coefficient and b is the bioturbation coefficient, z is depth, u is the porewater burial velocity and w is the solid burial velocity, α is the irrigation coefficient, vw is the concentration of a solute in the bottom waters and ΣR is the net production rate from all biogeochemical reactions for a given component. Each of these terms will be described in detail later in the Model description section. These partial differential equations are solved numerically using the method of lines described in Boudreau (1996b). Instead of searching for steady-state solutions directly, RADI computes the 140 concentrations of a set of solids and solutes at each depth and time steps following a time vector set by the user. The user determines the simulation time depending on the objectives, e.g., multimillennial to predict a steady state, or a few days to study the response of the sedimentary system to high frequency cyclic phenomena such as tides. For initial conditions, the user can choose between predefined uniform values (e.g., set all concentrations to zero) or a set of saved concentrations (e.g., from a previous simulation that has reached steady state). T is the total simulation time, dt is the temporal resolution, i.e., the interval 145 between each timestep, and t refers to the array of modelled timepoints. All time units are in years (a). The interface between the surface sediment and overlying seawater, conventionally set at a sediment depth z = 0, represents the top layer of RADI's vertical axis (Fig. 1). The bottom layer of the model is at a sediment depth Z. Between these limits, n layers are present, each being separated by a constant vertical gap dz. Depth units are in meters. The values assigned to dz and dt depend on the nature of the problem and on the kinetics of the chemical reactions. In the present study, all cases use dz = 2 mm and dt = 1/128000 150 a, i.e., ~4 minutes. If a lower dz is used, dt needs to be lowered as well to preserve numerical stability. In general, the ratio dz/dt should be kept below a value of 256 m/a. If dz is divided by two, dt needs to be divided by two as well, and the speed at which RADI runs will be reduced by a factor of four.  RADI operates on a static, user-defined porosity profile. Sediment porosity, φz in Fig. 1, refers to the porewater 180 volume fraction in the sediment (dimensionless) and typically decreases exponentially with sediment depth due to steady-state compaction. The sediment porosity profile is parametrized following Boudreau (1996b) as: where ∞ is the porosity at great depth, 0 is the porosity at the sediment-water interface, and β is an attenuation coefficient expressed in m -1 . A typical deep-sea sediment porosity profile is shown in Fig. 1. Here the measured porosity profile at station 185 7 of cruise NBP98-2 (Sayles et al., 2001) is fit using ∞ = 0.87, 0 = 0.915, and β = 33 m -1 . The solid volume fraction ( , dimensionless) is defined as: where R(vt,z) quantifies the rate of change of vt,z due to chemical reactions, A(vt,z) quantifies the rate of change of vt,z due to advection, D(vt,z) quantifies the rate of change of vt,z due to molecular and bio-diffusion, and I(vt,z) quantifies the rate of change 205 of vt,z due to bio-irrigation. In general, only the subscript zs are explicitly written out in this document, for variables and parameters that vary with depth. The ts are implicit but excluded for clarity.

Reactions
In RADI, biogeochemical reactions operate on solutes and solids throughout the entire sediment column, including the very top and bottom layers. R(vz) is the net rate at which vz is being consumed (negative R) or produced (positive R) by these reactions. Biogeochemical reactions in RADI (Table 3) are grouped into three categories: (i) organic matter degradation, 220 (ii) oxidation of reduced metabolites (organic matter degradation byproducts), and (iii) dissolution or precipitation of calcium carbonate minerals. RADI has been designed for early diagenesis in deep-sea sediments, so formation and re-oxidation of metal sulfide minerals are not considered.

Organic matter degradation
Organic carbon deposited on the seafloor originates mainly from primary production in the upper ocean or on land 225 and, to a lesser extent, from the ocean interior via chemoautotrophy. Despite the differences in origin, detrital organic matter found in marine sediments typically has the same composition: ~60% proteins, ~20% lipids, ~20% carbohydrates, and a fraction of other compounds (Hedges et al., 2002;Burdige, 2007;Middelburg, 2019). Here, the stoichiometry of organic matter is represented by the coefficients c (for carbon), n (for nitrogen), and p (for phosphorus). By default, the c:p ratio is set to 106:1 and the n:p ratio set to 16:1, following the Redfield ratio that describes the average composition of phytoplankton biomass 230 (Redfield, 1958), but these values can easily be adjusted. In RADI, c/p is denoted RC, n/p is denoted RN, and p/p is denoted RP, which is unity. Organic matter is also simplified here as an elementary carbohydrate (CH2O). In reality, loss of H and O during biosynthesis of proteins, lipids, and polysaccharides occurs (Anderson, 1995;Hedges et al., 2002;Middelburg, 2019), which results in an effective molar ratio of O2 consumed to C degraded of ~1.2 during aerobic respiration (Anderson and Sarmiento, 1994) instead of 1 as assumed here (Table 3). 235 Observations show that some organic compounds are preferentially degraded and become selectively depleted (Cowie and Hedges, 1994;Lee et al., 2000). As a result, the bulk reactivity of organic matter decreases with increasing age (Middelburg, 1989). Degradation of organic matter deposited at the seafloor typically follows a sequential utilization of available oxidants, O2, NO3 -, MnO2, Fe(OH)3, and SO4 2-, followed by methanogenesis (Froelich et al., 1979;Berner, 1980;Arndt et al., 2013). All organic matter degradation reactions implemented in RADI are shown in Table 3. 240 To account for the decrease in organic matter degradation rate with sediment depth, we separate organic matter into fractions of different reactivity, and we assign a rate constant to each of the degradable fractions. Following Jørgensen (1978), Westrich and Berner (1984), and Soetaert et al. (1996b), three different classes of organic matter are considered: refractory, slow-decay, and fast-decay. The refractory organic matter class is not reactive during the timescales considered here. The fastand slow-decay organic matter fractions each have a depth-dependent, oxidant-independent reactivity. The overall organic 245 matter degradation rate decreases with depth because the quantity of organic matter and the relative proportions of fast and slow-decay materials decline with depth. Organic matter is degraded following the sequential utilization of available oxidants.
The oxidant limitation is represented by a Michaelis-Menten-type (also termed 'Monod') function, in which each oxidant has an associated half-saturation constant (Koxidant in mol m -3 ) that symbolizes the oxidant concentration at which the process proceeds at half its maximal speed (Soetaert et al., 1996b). The presence of some oxidants may also inhibit other metabolic 250 pathways; this is represented by an inhibition constant (Koxidant′ in mol m -3 ) that is specific to each oxidant. These limiting and inhibitory functions have been widely used (Boudreau, 1996b;Van Cappellen and Wang, 1996;Soetaert et al., 1996b;Couture et al., 2010) and they allow a single equation to be used for each component across the entire model depth range, and also permit some overlap between the different pathways, as observed in nature (Froelich et al., 1979). In RADI, the overall degradation of fast-or slow-decay organic carbon occurs at a rate: 255 where POC is the rate constant for the degradation of a given type of organic carbon (fast-or slow-decay) expressed in a -1 , [ ] is the concentration of organic carbon (fast-or slow-decay) in sediments, and . is the sum of the fractions of organic carbon degraded by each oxidant (dimensionless, always equal to one), given by: (8f).
Half-saturation and inhibition constants for each oxidant used in RADI are given in Table 4. The degradation rate constant of organic carbon, POC , is computed as a function of the flux of organic carbon reaching the seafloor and is sedimentdepth dependent (Archer et al., 2002): 270 where FPOC is the total flux of organic carbon reaching the seafloor (i.e., fast, slow, and refractory), in mol m -2 a -1 . The numbers 1.3 × 10 -4 and 1.5 × 10 -1 have been tuned to best fit observations of both a Southern Ocean and North Atlantic Station, see the 275 Model evaluation section 3.

Oxidation of organic-matter degradation by-products
Organic matter degradation reactions primarily change oxidants (e.g., O2, NO3 -, MnO2, Fe(OH)3, SO4 2-) into their reduced forms (e.g., H2O, N2, Mn 2+ , Fe 2+ , H2S; Table 3). If oxygen is introduced into the system, or if the reduced metabolites diffuse upwards in oxic porewaters, then these reduced byproducts are converted back into their oxidized form and the energy 280 contained in them becomes available to the microbial community, though these energetics are not considered in RADI.
Here, four redox reactions involving organic-matter degradation byproducts are implemented ( Table 3): oxidation of each of Fe 2+ , Mn 2+ , ΣH2S and ΣNH3. These four reactions consume porewater total alkalinity (TAlk) but do not alter porewater ΣCO2 (Table 3), thus locally acidifying porewaters. Here, we use the TAlk definition of Dickson (1981), in which TAlk is defined as "the number of moles of hydrogen ion equivalent to the excess of proton acceptors (bases formed from weak acids 285 with a dissociation constant K ≤ 10 -4.5 and zero ionic strength) over proton donors (acids with K > 10 -4.5 ) in one kilogram of sample". This scheme should be sufficient for all open-ocean applications but may not be suitable for coastal and anoxic environments with extensive metal-sulfide mineral turnover, which require a more complete set of redox reactions such as that from the CANDI model of Boudreau (1996b). Additional components and reactions can easily be implemented in future versions (see Future Developments section 5). The rate constants for these four redox reactions are taken from Boudreau 290 (1996b) and reported in Table 4.

CaCO3 dissolution and precipitation
RADI includes two CaCO3 polymorphs: low-Mg calcite and aragonite, but more could be added in future versions, 310 e.g., high-Mg calcite and/or vaterite. Calcite and aragonite both have different dissolution kinetics, in which their dissolution rates increase as the undersaturation state of seawater with respect to calcite (1 -Ωca,z) or aragonite (1 -Ωar,z) increases (Keir, 1980;Walter and Morse, 1985;Sulpis et al., 2017;Dong et al., 2019;Naviaux et al., 2019b). Here, Ωz is the sediment-depthdependent saturation state of seawater with respect to calcite or aragonite, computed as [Ca 2+ ]z •[CO3 2-]z / K * sp, where K * sp is the stoichiometric solubility constant of calcite or aragonite at in situ temperature, pressure, and salinity, as given in Mucci 315 (1983) and Millero (1995). At each time step, Ωz is computed using porewater [Ca 2+ ]z and [CO3 2-]z from the previous time step, the latter being calculated as a function of TAlk and the proton concentration [H + ]. At each model time step, the total hydrogen ion concentration [H + ] is computed from TAlk and ∑CO2 using a single Newton-Raphson iteration from the previous time step (Humphreys et al., 2021):  (Uppström, 1974), plus equilibrium constants for silicic acid (Sillén et al., 1964) and phosphoric acid (Yao and Millero, 1995). Its derivative is computed following the approach of CO2SYS, see Humphreys et al. (2021). The carbonate ion concentration is then computed as: 325 where K1 * and K2 * are the first and second dissociation constants for carbonic acid, respectively, taken from Lueker et al.

(13).
In these expressions, the dissolution rate constant ( , in a -1 ) and the reaction order ( , unitless) are mineral-specific.
The dissolution rate constants implicitly account for each mineral's specific surface area. Similar formulations have previously been implemented to describe calcite dissolution rates (e.g. Archer et al. 2002) but, in most cases, with a high reaction order 335 and with a tuned rate constant independent of solution chemistry (Fig. 2). Such discretizations are convenient but lack a mechanistic description of the controls on calcite dissolution in seawater (Adkins et al., 2021).
The latest advances using isotope-labelling approaches to study carbonate dissolution kinetics show abrupt changes in dissolution mechanism depending on solution saturation state with either calcite or aragonite (Subhas et al., 2017;Dong et al., 2019;Naviaux et al., 2019a;Naviaux et al., 2019b). Close to equilibrium, dissolution occurs primarily at sites on the crystal 340 surfaces that are most exposed to the solution, e.g., steps and kinks. Further away from equilibrium, dissolution etch pits are activated at surface sites associated with defects and impurity atoms. Far away from equilibrium, there is enough free energy for dissolution etch pits to occur anywhere on the mineral surface, without the aid of crystal defects (Adkins et al., 2021).
However, at temperatures most relevant to the deep oceans, ~5 °C or less, the defect-assisted dissolution mechanism is skipped (Naviaux et al., 2019b) and only the step-edge retreat (close to equilibrium) and homogeneous etch-pit formation (far away 345 from equilibrium) dissolution regimes remain (Naviaux et al., 2019b) (Fig. 2). For both aragonite and calcite, while homogeneous etch-pit formation is indeed associated with a high-order dependency on the solution saturation state, step-edge retreat dissolution rates dominating near equilibrium show very little dependence on seawater saturation (Dong et al., 2019;Naviaux et al., 2019a). This could have significant consequences for the predicted carbonate dissolution rate near equilibrium: saturation-state independent step-edge retreat dissolution will always be predicted to be faster close to equilibrium than 350 dissolution associated with a high reaction order because a high reaction order forces the dissolution rate to converge to zero as the solution gets closer to equilibrium (Fig. 2). The rate constants are tuned to best fit the observations in the two stations presented in the Model Evaluation sections 3. We use kdiss. ca. = 6.3 x 10 -3 a -1 for 0.828 < Ωca < 1, kdiss. ca. = 20 a -1 for Ωca ≤ 0.828, kdiss. ar. = 3.8 x 10 -3 a -1 for 0.835 < Ωar < 1, and 370 kdiss. ar. = 4.2 x 10 -2 a -1 for Ωar ≤ 0.835. Both calcite and aragonite dissolution rate constants are lower than the values reported in the original publications. We suspect that i) the reactive surface area of grains in sediments is much smaller than their specific surface area measured using adsorption isotherms via the BET method and ii) unaccounted dissolution inhibitors are present in sediments, such as dissolved organic carbon (Naviaux et al., 2019a). A comparison of the steady-state [CO3 2-] and [Calcite] porewater profiles predicted by RADI using the tuned rate constants implemented in RADIv1 and the original rate 375 constants is shown in Fig. S1.
Calcite precipitation is also included in the model and its rate (solid blue line in Fig. 2b) is parameterized with the following function: . ., where kprec. ca is the precipitation rate constant set to 0.4 mol m -3 a -1 and η is equal to 1.76. The precipitation reaction order is 380 taken from Zuddas and Mucci (1998), corrected for a seawater-like ionic strength of 0.7 mol kg -1 . The precipitation / dissolution rate continuum implemented in RADI (see Fig. 2) is very different from what a classic model with only calcite dissolution following high-reaction order kinetics would display. For comparison, the dissolution rate of calcite using a dissolution rate constant kdiss of 100 % per day and a reaction order η of 4.5, as implemented in most diagenetic models, including Archer (1991), is shown in Fig. 2a. Mechanistic interpretations of the "kinks" in the dissolution rate profiles and of a non-zero 385 dissolution rate near equilibrium still require more research but the implications of these features for our understanding of marine CaCO3 cycles can be explored with the present model.

Advection
The solid burial velocity at the sediment-water interface, w0 in m a -1 , is given by: where is the flux of a solid species at the sediment-water interface (mol m -2 a -1 ), is the molar mass of that solid (g mol -1 ), and is its solid density (g m -3 ). The solid and porewater burial velocity at greater depth are assumed to be equal, and are computed as: Thus, the porewater burial velocity, u, at all depths is: 395 and the solid burial velocity, w: Depth profiles of u and w are shown in Fig. 1, computed from the solid fluxes at station 7 of cruise NBP98-2 (Sayles et al., 2001), seen in Model Evaluation Section 3.2. In Fig. 1, the sharp porosity decline in the top centimeters of the sediments causes 400 the solid fraction at ~5-cm depth to be roughly 50% higher than just below the interface. This leads to a solid burial velocity decrease of about the same magnitude (Fig. 1).
Advection is implemented following Boudreau (1996b), where the advection rate (Az, in mol m -3 a -1 ) for solutes is given by: where dz is the effective diffusion coefficient for a given solute at a given depth (in m 2 a -1 ), d° is the "free-solution" molecular diffusion coefficient for a given solute (in m 2 a -1 ), and θz is the depth-dependent tortuosity (unitless) defined in the Diffusion section 2.4. For solids, a more sophisticated weighted-difference scheme is employed (Fiadeiro and Veronis, 1977;Boudreau, 1996b): where bz is the depth-dependent bioturbation coefficient (m 2 a -1 ) and The parameter Peh is half of the 'cell Peclet number', which expresses the influence of advection relative to bioturbation across 415 a distance separating two points of the grid, centered at the depth z. If bioturbation dominates (Peh << 1), e.g., near the sedimentwater interface, σz tends toward zero and a centered-difference discretization is implemented. If advection dominates (Peh >> 1), e.g., deeper in sediments, σz tends toward unity and backward-difference discretization prevails, see Eq. (20). This differencing scheme, originally developed by Fiadeiro and Veronis (1977), maintains stability in the entire sediment column (Boudreau, 1996b). 420

Diffusion
The diffusion flux of any species depends on its effective diffusion coefficient, dz(v), which varies with depth within the sediment.
For each solute, free-solution diffusion coefficients, denoted dz°(v), were computed at in-situ temperatures (Li and Gregory, 1974;Boudreau, 1997;Schulz, 2006). For solute variables representing several individual species (e.g., ΣPO4, ΣCO2), 425 the diffusion coefficient of the dominant species was considered. Given the high proportion of HCO3relative to CO3 2and CO2 (aq) in seawater and porewaters (see Fig. S2), the diffusion coefficient of HCO3was adopted for both TAlk and ΣCO2.
However, this approach may not be suited for sedimentary environments in which pH is lower than 7, because a greater proportion of dissolved inorganic species would then be under the form of carbonic acid, or CO2 (aq), which has a higher diffusion coefficient than HCO3 - (Fig. S2). Free-solution diffusion coefficients, their temperature dependencies and their 430 sources are reported in Table 5. The diffusion of solutes in the porewaters is slower than in an equivalent volume of water as a result of the physical barriers caused by the presence of solid grains in a sediment. To correct for this effect, we follow Boudreau (1996b) and compute the effective diffusion coefficient for a given solute as: where so-called tortuosity (θ) is defined as (Boudreau, 1996a): 435 = √1 − 2 ln( ) (24). increasing flux of total organic carbon reaching the seafloor ( POC ) and attenuates in low-oxygen conditions. This pattern was also observed by Smith et al. (1997) and Smith and Rabouille (2002). As in Archer et al. (2002), we couple both bioturbation and irrigation to the incoming carbon deposition flux (Fig. 3) rather than water depth or sediment accumulation rate (Boudreau, 1994;Middelburg et al., 1997;Soetaert et al., 1996c), although all these quantities are related to each other. From an ecological perspective, more carbon to the seafloor represents more food available to benthic communities, hence more biological 460 transport. Linking bioturbation activity to carbon deposition flux also allows for a direct coupling with Earth system models  (Boudreau, 1997;Schulz, 2006) dz°(Fe 2+ ) 0.001076 + 0.000466 × Tw (Li and Gregory, 1974) dz°(Mn 2+ ) 0.009625 + 0.000481 × Tw (Li and Gregory, 1974) simulating carbon sinking fluxes in the ocean. Following Archer et al. (2002), we express the surficial bioturbation mixing rate (b0, in m 2 a -1 ) as: b 0 = (2.32 • 10 −6 )( POC • 10 2 ) 0.85 where POC is expressed in mol m -2 a -1 . The bioturbation mixing rate at all depths (bz, in m 2 a -1 ) is: 465 where the characteristic depth = 8 cm, following Archer et al. (2002), and [O 2 ] is the oxygen concentration in the bottom waters. This depth-dependent bioturbation mixing rate is common to all solids and its depth distribution is shown in Fig. 3 as a function of in situ [O 2 ] and POC . The effective diffusion coefficient for solids is then set as: . 470 (Bio)diffusion is implemented in RADI following the centered difference discretization scheme from Boudreau (1996b). At sediment depth z, where 0 < z < Z, for both solutes and solids: where dz(v) is the relevant effective diffusion coefficient.

Irrigation 480
The mixing of solutes caused by burrow flushing or ventilation occurs through an ensemble of processes collectively termed irrigation. Macroscopic burrows are often present in the seafloor sediment, with a complex three-dimensional structure and filled with oxygenated water that is ventilated for aerobic respiration. In a one-dimensional framework, this causes apparent internal sources or sinks of porewater solutes at particular depths (Boudreau, 1984;Emerson et al., 1984;Aller, 2001).
Mathematically, this is parameterized as a non-local exchange function, i.e., a first-order kinetic reaction: 485 and the irrigation coefficient at all depths as: where the characteristic depth is 5 cm (Archer et al., 2002). The depth-distribution of the irrigation coefficient is shown in

Boundary conditions
Modeling of advection and diffusion processes requires appropriate boundary conditions in the layers above and below (zdz and z + dz, respectively). Effective values of each variable immediately adjacent to the modelled depth domain are calculated following Boudreau (1996b) and used to compute the effects of advection and diffusion in the top and bottom layers using the same equations as within the sediment itself. 500 At the sediment-water interface, RADI enables prescribed solid fluxes and a diffusive boundary layer control for solutes. Following Boudreau (1996b), we calculate advection and diffusion at z = 0 for solutes and solids as: and v (-dz) respectively. Here, θ is the tortuosity, δ is the boundary layer thickness (expressed in m, see Fig. 1), and vw is the solute concentration above the diffusive boundary layer, i.e., in the bottom waters. At the sediment depth z = Z, v(Z+dz) falls outside the depth range of the model. The bottom boundary condition demands that concentration gradient disappear, which can be translated by the following for both solutes and solids: . 510 This 'no-flux' bottom boundary condition should be appropriate here because we set Z so that all action occurs at shallower depth. However, if anaerobic methane oxidation or subsurface weathering are included in future versions, a 'constant' flux boundary conditions might need to be included.

Julia and MATLAB/GNU Octave implementations
We have implemented RADI both in Julia (Humphreys and Sulpis, 2021) and in MATLAB/GNU Octave (Sulpis et 515 al., 2021). Both implementations use similar nomenclature and provide identical results. Documentation for both is available from https://radi-model.github.io. The Julia implementation is available from https://github.com/RADI-model/Radi.jl and the MATLAB/GNU Octave implementation is available from https://github.com/RADI-model/Radi.m. Julia (https://julialang.org) is a high-level, high-performance, and cross-platform programming language that is free and open source (Bezanson et al., 2017). Its high-performance stems primarily from just-in-time (JIT) compilation of code 520 before execution, which has been built-in since its origin. RADI uses Julia's multiple-dispatch paradigm, a core feature of the language, which improves the readability of the code and reduces the scope for errors. Specifically, each modelled component of the sediment column is either a porewater solute or a solid. These components are initialized in the model as variables either of Solute or Solid type. Advection and diffusion are governed by different equations for porewater solutes than for solids but the same top-level functions (advect! and diffuse!) can be used within RADI to calculate the effects of these 525 processes for both component types; the multiple-dispatch paradigm ensures that the correct equations are automatically used on the basis of the type of the input variable. While the model has been designed to solve a single profile at a time, Julia's support for parallelized computation (across multiple processors) would also support efficient computations across a series or grid of vertical profiles.
As of version R2015b, MATLAB also features JIT compilation with a corresponding execution speed-up. However, 530 MATLAB is an expensive, proprietary software, which limits how widely it can be used. The MATLAB implementation also runs in GNU Octave (https://www.gnu.org/software/octave/), which is a free and open-source clone of MATLAB. However, GNU Octave executes more slowly than MATLAB for a variety of reasons, including a lack of JIT compilation.
For a model that necessarily includes long simulations with relatively short time-steps, computational speed is an important consideration. Our testing indicates that the Julia implementation runs ~3 times faster than the MATLAB (R2020a) 535 implementation and ~70 times faster than the GNU-Octave implementation.
Simultaneously developing the model in two languages allowed us to quickly identify and remedy bugs and typographical errors in both implementations. Each was coded independently, with equations and parameterizations writtenout from the original sources, thus avoiding code copy-and-paste errors. Frequent comparisons were made throughout this process to ensure that the results were consistent. For a typographical error to survive to the final models would therefore 540 require an identical mistake to have been made independently in both implementations. The risk of such errors is thus substantially reduced by our dual-language approach. Where errors were identified, in some cases they were subtle enough that they may otherwise not have been noticed, while still causing meaningful errors in final model results.

Model evaluation
To evaluate the performance of RADI, we used in-situ data obtained at three different locations and compared our 545 predictions to the measured porewater and sediment solid phase composition profiles. We used these comparisons to tune the CaCO3 dissolution and POC degradation rate constants; all other parameters were assigned a priori using values from the literature. Thus, we did not aim to reproduce observations as well as possible by tuning a wide selection of parameters. Instead, we evaluated whether a generic approach using measured deposition fluxes and bottom-water conditions could explain observations while tuning only the inorganic and organic reactivity constants. 550

North-western Atlantic Ocean
First, RADI was compared to the porewater and sediment composition measurements of station #9 described in Hales et al. (1994), located in the North-western Atlantic Ocean (24.33°N / 70.35°W) at a 5210-m depth. The bottom-water TAlk and ΣCO2 were 2342 μmol kg -1 and 2186 μmol kg -1 , respectively, bottom-water in situ temperature was 2.2°C, salinity was 34.9, and oxygen concentration was 266.6 μmol kg -1 (Hales et al., 1994). The computed bottom-water saturation state with 555 respect to calcite was 0.88. The only CaCO3 polymorph reaching the seafloor was assumed to be calcite. The calcite flux to the seafloor was set to 0.20 mol m -2 a -1 (20.02 g m -2 a -1 ) and the POC flux to 0.18 mol m -2 a -1 , which correspond to the low end of fluxes measured by sediment traps on the continental slope (Hales et al., 1994). The clay flux was set to a value of 26 g m -2 a -1 to fit the calcite sediment surface concentration measured by Hales et al. (1994). The porosity at the sediment-water interface was set to that measured by Sayles et al. (2001) in the Southern Pacific Ocean station, see Fig. 1. Following the 560 diffusive boundary layer distribution from Sulpis et al. (2018), δ at the station location was set to 938 μm. This value represents an annual-mean estimate derived using a number of assumptions, e.g., considering the sediment-water interface to be a horizontal surface and neglecting sediment roughness. A complete description of the environmental parameters for this North-Atlantic station, along with their sources, is available in Table S1.
RADI was run using the environmental conditions described above and the steady-state concentration profiles of O2, 565 ΣNO3, calcite, and POC were compared with observations. Complete methods for solute and solids measurements are described in Hales et al. (1994). Briefly, porewater O2 concentration was measured both in-situ using microelectrodes and onboard (along with ΣNO3) from the retrieved box core (Hales et al., 1994). The steady-state calcite, TAlk and ΣCO2 profiles were compared to those obtained from a RADI simulation with "traditional", 4.5-order calcite dissolution kinetics, see Fig. 2, all other variables being unchanged. 570 (Hales et al., 1994).

575
RADI predicts a porewater O2 concentration decreasing from the bottom-water value to zero at ~20 cm-depth (Fig.   4). In the top 2 cm, the RADI porewater O2 predictions near the surface are in good agreement with the in-situ microelectrode measurements. The RADI-predicted [O2] is lower than that measured on-board, but [ΣNO3] is well reproduced by RADI.
RADI predicts that organic matter respiration is mainly aerobic (see Table 3a) until about 20 cm-depth. Between 20 and 35-580 cm depth, ΣNO3 is the preferred oxidant for organic matter degradation (see Table 3b), which leads to a decrease in porewater [ΣNO3] in both RADI predictions and observations, followed by ΣSO4 deeper than 35-cm depth. The calcite profile is relatively well reproduced by RADI in the top 20 cm but the measured calcite concentration drop below 20-cm depth is not well reproduced. When "traditional" 4.5-order calcite dissolution kinetics are implemented, calcite concentrations are underestimated by about a third in the top 20 cm. The observed lower calcite concentrations below 20-cm depth may be 585 attributed to a lower calcite accumulation rate to the seafloor in the past, whereas the model considers accumulation of solids to be unchanged through time. Calcite concentration predicted by RADI does increase again below 25-cm depth due to porewater supersaturation (Ωca at 40 cm-depth is about 1.05), but this increase is too small to be noticed on the figure.

Southern Pacific Ocean
RADI was also compared with data collected at the station #7 mooring #3 described in Sayles et al. (2001) The bottom-water chemical composition was taken from the GLODAPv2 1x1° climatologies (Lauvset et al., 2016), 595 linearly interpolated over depth, latitude, and longitude to match the station location and seafloor depth. The bottom-water insitu temperature was 0.84°C, salinity was 34.696, oxygen concentration was 215.7 μmol kg -1 , and calculated saturation state with respect to calcite was 0.85. Solid fluxes at this station were measured by Sayles et al. (2001) using sediment traps collecting sinking particles between the months of November and December 1997. Their deepest sediment trap available was at a depth of 3257 m, i.e. 600 meters above the seafloor, which we assume to be representative of sinking fluxes to the seafloor, 600 although the loss of material after collection usually causes sediment traps to underestimate the real sinking fluxes (Buesseler et al., 2007). The only CaCO3 polymorph reaching the seafloor was assumed to be calcite and its flux was set to 0.25 mol m -2 a -1 (25.02 g m -2 a -1 ), rather than using the sediment-trap value, in order to fit the calcite sediment surface concentration measured by Sayles et al. (2001). This CaCO3 flux to the seafloor is slightly higher than the measured CaCO3 flux at 600 meters above the seafloor in mid-January 1997 (0.19 mol m -2 a -1 ; Sayles et al. 2001). The POC flux was set to 0.14 mol m -2 a -605 1 (4.62 g OM m -2 a -1 ), which is also slightly higher than the measured POC flux averaged between the months of November and December 600 meters above the seafloor (0.11 mol m -2 a -1 ). The clay flux, which we considered to be the total measured sediment flux minus the assumed POC and calcite fluxes, was 32 g m -2 a -1 . The porosity profile was tuned to best fit the porosity measurements at this station (Sayles et al., 2001, see Fig. 1). Finally, using the diffusive boundary layer world map computed in Sulpis et al. (2018) based on bottom current speeds at in-situ temperature and pressure measurements, the diffusive 610 boundary layer thickness (δ) at the station location was set to 715 μm. A complete description of the environmental parameters for this station, along with their sources, is available in Table S2. RADI was run using the environmental conditions described above to compare the steady-state concentration profiles of TAlk, ΣCO2, O2, ΣNO3, calcite, and POC to in-situ measurements. Methods for solutes and solids concentration analyses are described in Sayles et al. (2001). Briefly, TAlk, ΣCO2, and ΣNO3 were sampled in situ using the Woods Hole Interstitial 615 Marine Probe (Sayles, 1979) while O2 was sampled at a higher depth-resolution but in the ship laboratory using whole-core squeezing (Bender et al., 1987). We also compare the RADI steady-state concentration profiles with those obtained from a simulation with "traditional" 4.5-order calcite dissolution kinetics, all other variables being the same. RADI predicts porewater O2 concentrations that are slightly higher than observed (Fig. 5). Because RADI does not predict porewater O2 to go to zero until about the 30-cm depth, the dominant organic matter degradation pathway switches from mainly aerobic to ΣNO3 at about the 30-cm depth. Nevertheless, the RADI-predicted ΣNO3 profile is lower than observed 625 values, particularly toward the bottom of the resolved depth. The TAlk and ΣCO2 porewater profiles predicted by a RADI simulation using 4.5-order calcite dissolution kinetics are very similar to those using the new calcite dissolution kinetics scheme, but the predicted calcite concentration is lower.

Central equatorial Pacific Ocean
As a third case study to evaluate the performance of RADI, solute fluxes through the sediment-water interface 630 computed from model steady-state runs were compared to fluxes measured using benthic chambers. The comparison took place at station #W-2 described in Berelson et al. (1994) and Hammond et al. (1996), located in the Central equatorial Pacific Ocean (0°N / 139.9°W) at a depth of 4370 m. Bottom-water in-situ temperature was set to 1.40°C and salinity to 34.69 (Lauvset et al., 2016). Bottom-water oxygen concentration was set to 159.7 μmol kg -1 and the bottom-water saturation state with respect to calcite computed using the carbonate system solver within RADI was 0.78. For the purposes of this evaluation, the CaCO3 635 flux to the seafloor was assumed to be entirely calcite. The calcite flux was set to 0.22 mol m -2 a -1 , which represents 22.02 g of calcite m -2 a -1 , the POC flux was 0.20 mol m -2 a -1 , that is, 6.6 g of organic matter m -2 a -1 , and the clay flux was set to 2.0 g m -2 a -1 . The steady-state calcite content within the top cm was 61 dry wt %, in line with CaCO3 contents observed in this area (Archer, 1996;Hammond et al., 1996). The porosity profile was built using an attenuation coefficient β = 33 m -1 , 0 = 0.85, which is the measured surface porosity (Hammond et al., 1996), and ∞ = 0.74, which is the measured porosity at depth 640 (Berelson et al., 1994), see Eq.
(3). The DBL thickness, δ, was fixed to a value of 1 mm. A complete description of the environmental parameters for this station, along with their sources, is available in Table S3.
The diffusive fluxes for a given solute (Jv) between the sediment-water interface and the bottom waters occur as a response to the concentration gradient within the DBL and can be expressed by: where v0 and vw are solute concentrations at the sediment-water interface and in bottom waters, respectively. In this definition, a positive Jv indicates a solute release from the sediment porewaters to the bottom waters while a negative Jv represents a solute flux towards the sediment.

Discussion on model performance
These three model evaluation examples allowed us to determine a set of organic carbon degradation-rate constants, CaCO3 dissolution-rate constants and organic-carbon flux composition (fast-decay, slow-decay or refractory) that can best reproduce sediment and porewater measurements in all stations while keeping all other model parameters to values from the 660 literature. In each station, bottom-water composition was fixed from observations. Due to a lack of adequate data, solid deposition fluxes were tuned to best fit observed CaCO3 and POC contents in sediments, except in the North-western Atlantic station where the CaCO3 and POC fluxes were taken from measurements and in the Southern Pacific station where the clay flux was inferred from observations. The POC composition that allows to best fit porewater and sediment data in the three stations was as follows: 70% of fast-decay POC, 27% of slow-decay POC, and 3% of refractory POC. The tuned fast-and 665 slow-decay POC degradation rate constants are reported in Organic matter degradation section 2.2.1 and are of similar orders of magnitude as in most other models (Arndt et al., 2013). The tuned CaCO3 dissolution rate constants are reported in CaCO3 dissolution and precipitation section 2.2.3 and are two orders of magnitude lower than their laboratory-based values (Fig. S1), which we attribute to the presence of dissolution inhibitors (e.g., dissolved organic carbon, Naviaux et al., 2019a) or to the reactive surface area of natural grains in situ being lower than in laboratory experimental settings. Thus, with its current 670 settings, RADI should be able to accurately predict porewater chemistry and sediment composition in deep-sea environments, provided that the POC and CaCO3 deposition fluxes are known.
In the central equatorial Pacific, all RADI diffusive fluxes are within the uncertainty range of observations except the ΣNO3 flux, which is underestimated by RADI. The low ΣNO3 flux could be attributed, for example, to the presence of organic matter with a stoichiometry different than the Redfield ratio used in the current version of RADI or to errors in the nitrification 675 parameters.
In addition, we note than the choice of calcite dissolution kinetics implemented in RADI does not seem to have a large impact on TAlk and ΣCO2 porewater profiles, but rather affect the predicted calcite concentration. RADI's step-edge retreat dissolution regime and its low reaction order induce calcite dissolution rates near equilibrium that are orders of magnitude higher than what is predicted in a high-order rate law (Fig. 2). That may cause the predicted calcite concentration 680 to be higher when step-edge retreat dissolution is included in the kinetics than when a high-order rate law is implemented, because porewaters can get closer to equilibrium with respect to calcite more easily.

Potential model applications
In the following section, we continue to analyze the results obtained using the environmental conditions of the equatorial Pacific station #W-2 and present a few examples of relevant model applications. RADI can be used to study both 685 steady state and transient conditions but in the following subsections we focus on time-dependent problems, since transient diagenetic models are underrepresented in the literature.

Seasonal variability
At the seafloor, the fluxes of sinking material regulating the chemical composition of sediment porewaters are patchy in both space and time. Seafloor microbes and macrofauna respond quickly to pulses of organic matter delivery to the seafloor 690 (Smith et al., 1992), causing short-term variability of sediment oxygen consumption (Smith and Baldwin, 1984;Smith et al., 1994). In addition, both the POC and CaCO3 fluxes to the deep seafloor are strongly affected by seasonal flux variability (Billett et al., 1983;Smith and Baldwin, 1984;Lampitt, 1985;Lampitt et al., 1993;Lampitt et al., 2010). In the northeast Atlantic Ocean, at 3000-m depth, Lampitt et al. (2010) have shown that the summer POC and CaCO3 fluxes can be ~ 10 and ~ 4 times higher, respectively, than the winter-time minima. The seasonal coupling between organic matter and CaCO3 fluxes 695 to the seafloor and the state of upper-ocean ecosystem is the result of rapid vertical transport of these materials (Sayles et al., 1994). If the fluxes of reactive material reaching the seafloor is affected by seasons, early diagenesis could display a seasonal signal and this should be taken into account when interpreting sedimentary data (Martin and Bender, 1988;Sayles et al., 1994;Soetaert et al., 1996a). Here we use RADI to assess how the porewater chemistry and solid composition of deep-sea sediments may be impacted by seasonally varying fluxes. 700 Seasonally time-varying solid fluxes to the seafloor (F, in mol m -2 a -1 ) can be represented with the following sinusoidal function: where t is time in years, Δt is the time period separating two maxima (here set to 1 year), and ΔF is the amplitude. We assume that all CaCO3 settling at the seafloor is calcite and set the mean FCalcite to 0.22 mol m -2 a -1 , its amplitude change ΔFCalcite to 705 0.11 mol m -2 a -1 , the mean FPOC to 0.20 mol m -2 a -1 , and its amplitude change ΔFPOC to 0.10 mol m -2 a -1 . All other parameters correspond to values from the Central equatorial Pacific stations described in the Model evaluation section 3.3.
Seasonal variations in the calcite and POC fluxes reaching the seafloor are visible in sediment profiles of both solids and solutes (Fig. 7). Nonetheless, the amplitude of concentration changes separating an annual minimum from an annual maximum is very small, barely if (at all) detectable by observations. The annual amplitude is about 0.5 wt% for calcite, 0.07 710 wt% for POC, 2 μmol kg -1 for [O2], and 7 μmol kg -1 for [ΣCO2]. While concentrations at the sediment-water interface respond quickly to seasonal flux changes, there is a phase lag that increases with depth between the concentrations of solids and porewater solutes and the seasonally changing fluxes to the seafloor. Thus, it is possible that porewaters and solid particles at the top mm-thick sediment layers are never really at a steady state but always lagging behind seasonal changes, even if these are minimal. This is in agreement with earlier modeling (Martin and Bender, 1988) and observational (Sayles et al., 1994) 715 studies, indicating that biogeochemical reactions and bioturbation at the deep seafloor are too slow to show a discernible seasonal signal. However, this might not be the case for sites receiving more reactive organic matter (Soetaert et al., 1996a).

Tidal cycles
This section explores the applicability of RADI for studying the response of sediments to higher frequency phenomena such as tides. The DBL thickness is dependent on the overlying current speed (Levich, 1962;Santschi et al., 1983;Lorke et al., 2003): slower currents generate thicker DBLs whereas faster currents cause the DBL to thin (Larkum et al., 2003;725 Lorke et al., 2003;Higashino and Stefan, 2004). In the deep-sea, an important contributor of benthic current speeds are tidal forces, which makes tidal currents a potentially important contributor to biogeochemical exchanges across the sediment-water interface (Egbert and Erofeeva, 2002;Sulpis et al., 2019). If tidal current-speed fluctuations induce DBL thickness fluctuations, they may induce solute concentration fluctuations at the sediment-water interface, thus affecting early diagenesis. The strongest tidal currents occur during the transition from high to low tides. For semidiurnal tides, the time period separating a low from 730 a high tide is ~6 hours (Pugh, 1987). Setting the average DBL thickness to δ = 1mm and assuming that tides generate δ fluctuations with an amplitude Δδ = 0.5mm, the time-dependent δ can be expressed as: where t is time in years, and Δt is set to 1/1461 a ( ~6 hours). RADIv1 was run using the steady state solutes and solids depth profiles from the equatorial Pacific station #W-2 as initial conditions, with a DBL thickness fluctuating in response to tidal 735 currents computed using Eq. (37).

740
While none of the solids seem to respond to tidal velocity fluctuation, due to their slow accumulation rate, solutes show a clear response (Fig. 8). At the sediment-water interface, the simulated porewater The implications are potentially important for our interpretation of porewater microprofiles. pH (Archer et al., 1989b;Cai and Reimers, 1993;Zhao and Cai, 1999;Cai et al., 2000), pCO2 (Cai et al., 2000;Zhao and Cai, 1997), CO3 2- (de Beer et al., 2008;Han et al., 2014;Cai et al., 2016), O2 (Revsbech et al., 1980;Reimers, 1987;Archer et al., 1989a;Sosna et al., 2007), and even dissolved Fe, Mn or S(-II) (Brendel and Luther, 1995) microelectrodes have been developed during the past decades. 750 According to the results presented here, microprofiles, which capture instantaneous snapshots of porewater chemistry, should show appreciable differences depending on when they are carried out during a tidal cycle. That organic matter degradation rates inferred from oxygen microprofiles span a wide range (Archer et al., 1989a;Arndt et al., 2013;Wenzhöfer et al., 2016) may be, among other factors, due to the dependency on tidal and other ocean-bottom current fluctuations. To adequately capture O2 consumption rate in sediments, O2 fluxes should be measured and integrated over a period of time longer than a 755 tidal cycle (Berg et al., 2022).

Benthic chambers
RADI can also be used in the calibration of sensors and the optimization of sampling protocols and experimental designs. In the DBL, molecular diffusion is the dominant mode of solute transport and laboratory experiments of CaCO3 dissolution in seawater suggest that diffusion through the DBL is the rate-limiting step for CaCO3 dissolution at the seafloor 760 in the absence of organic matter respiration (Sulpis et al., 2017;Boudreau et al., 2020). Nevertheless, earlier assessments of in-situ CaCO3 dissolution at the sediment-water interface in the central equatorial Pacific indicated that DBL thickness does not impact overall dissolution rates (Berelson et al., 1994).
In their study, Berelson et al. (1994) deployed a set of free-sinking benthic chambers onto the seafloor. In each chamber, the portion of the chamber exposed above the sediment-water interface was sealed and isolated from external bottom 765 waters, and water samples were drawn during the incubation period. Each incubation lasted between 80 and 120 hours.
Chambers were stirred with a paddle at various rates to quantify the dependency of the measured diffusive fluxes across the sediment-water interface on the DBL thickness, which were calibrated via anhydrite dissolution to the 300-to-600 μm range.
Seeing no influence of the stirring rate on the measured diffusive fluxes, Berelson et al. (1994) discarded the hypothesis of fast, surficial carbonate dissolution and instead argued for slow, high-reaction-order calcite dissolution kinetics at the seafloor, 770 as subsequently implemented in most models. To better interpret the results from a benthic-chamber study such as that of Berelson et al. (1994), RADI can be used to predict the time-response of diffusive fluxes across the DBL following an instantaneous change of DBL thickness due, for instance, to a change in paddle stirring rate within a chamber. RADI was run using the steady-state solutes and solids depth profiles from the equatorial Pacific station #W-2 as initial conditions. In the initial run, the DBL thickness was set to 1 mm. We simulated the response of this sediment to an 775 instantaneous change in the DBL thickness, with one model run representing a situation where δ increases from 1 to 5 mm (e.g., a slow stirring rate) and one model run representing a δ drop from 1 mm to 200 μm (e.g., a fast stirring rate). Following a 5-fold increase in δ, diffusive fluxes of TAlk, ΣCO2, and O2 initially decrease by a factor ~5 but then increase back as the solute concentrations at the interface adapt to the new DBL (Fig. 9). Two hours after the δ increase, diffusive fluxes converge toward a new steady state. Following a 5-fold decrease in δ, diffusive fluxes immediately increase by the same magnitude, but 780 go back close to their initial value within an hour as the interfacial porewater concentrations adjust to the new DBL. These results suggest that the incubation periods of the Berelson et al. (1994) benthic chamber experiments were long enough to let porewaters adjust to the changes caused by the paddle stirrers and confirm that, in an organic-matter and CaCO3-rich sediment such as in the Equatorial Pacific, the influence of a DBL on diffusive fluxes across the sediment-water interface should indeed be limited. Additionally, these results confirm the observation by Berelson et al. (1994)

Additional applications
The time-dependent problems presented above focus on relatively short timescales (from minutes to months). A non-795 steady-state model such as RADI can also be used to project the sediment response to perturbations over longer periods of time. Examples include estimating the effect of negative emission technologies such as coastal enhanced weathering with olivine on early diagenesis Montserrat et al., 2017), of deep-sea mining (Haffert et al., 2020) or bottom trawling (Trimmer et al., 2005;van de Velde et al., 2018;De Borger et al., 2021), the impacts of a decadal bottomwater deoxygenation event such as in the Saint Lawrence estuary (Jutras et al., 2020), or the present anthropogenic CO2 800 transient. However, the current version of RADI cannot deal with long-term major dissolution (erosion) events, because the burial velocity calculation scheme currently implemented does not account for solid mass gain or loss within the sediment. To study the response of sediments to a global-scale ocean acidification event such as the Paleocene-Eocene Thermal Maximum (Zachos et al., 2005;Cui et al., 2011), a different burial velocity calculation scheme would have to be implemented, such as that adopted by Munhoven (2021). 805

Future developments
One advantage of RADI is that it is easily tunable by the user: adding new components is straightforward as long as the chemical reactions are known. In future releases, we plan to add oxygen, carbon, and calcium isotopes as individual components in order to predict the diagenetic response of isotopic signals. Additionally, adsorption/desorption reactions on clay surfaces could be a critically important advance, especially regarding the prediction of sedimentary pH profiles (Meysman 810 et al., 2003), as RADI currently treats clay minerals as non-reactive.
The representation of organic matter in the current version of the model is oversimplified. All reactive organic matter in RADI is associated with a 'Redfield' stoichiometry but marine organic matter can considerably deviate from this ideal (Martiny et al., 2013;Teng et al., 2014). Finally, a model is only as good as its assumptions. RADI is targeted to study deep-sea, carbonate sediments. To be 815 used in coastal environments, additional biogeochemical reactions would be necessary, particularly those involving methane and iron sulfide. Close to the shore, sediments become more permeable and the assumption of molecular diffusion as the dominant mode of solute transport in porewaters does not hold. In very shallow environments that are subject to high wave energy, pressure-induced advection in the sediment porewaters also needs to be included (Huettel et al., 2014). Moreover, coastal sediments have typically lower pH than open-ocean sediments, which may render our assumption of both ΣCO2 and 820 TAlk diffusing with a fixed diffusion coefficient set to that of the HCO3ion inaccurate, see Fig. S2.

Code availability
The current versions of RADI in both Julia and MATLAB/GNU Octave are freely available from GitHub (https://github.com/RADI-model) under the GNU General Public License v3. The exact version of the model used to produce the results used in this paper is archived on Zenodo (RADI.jl v0.3; https://doi.org/10.5281/zenodo.5005650; v1 will be released 825 after review), along with input data and scripts to run the model for all the simulations presented in this paper. RADI users should cite both this publication and the relevant Zenodo reference (Humphreys and Sulpis, 2021;Sulpis et al., 2021).

Data availability
Sediment and porewater composition, porosity, and solid fluxes data for the Southern Pacific Ocean station described in Sayles et al. (2001) are available at http://usjgofs.whoi.edu/jg/dir/jgofs/southern/nbp98_2/. Sediment and porewater composition for 830 the North-western Atlantic Ocean station described in Hales et al. (1994)  Conceptualization, methodology, validation, resources, writingreview and editing, supervision 840

Competing interests
The authors declare that they have no conflict of interest.

Acknowledgements
Thanks are due to Bernard P. Boudreau whose CANDI model (Boudreau, 1996b) was a large source of inspiration during the creation of the present RADI model, and to Daniel L. Johnson for fruitful discussions. We thank David Burdige and one 845 anonymous reviewer for their constructive feedback. We also thank Lukas van de Wiel for assistance with the Utrecht