Model of Early Diagenesis in the Upper Sediment with Adaptable complexity – MEDUSA (v. 2): a time-dependent biogeochemical sediment module for Earth System Models, process analysis and teaching

. MEDUSA is a time-dependent one-dimensional numerical model of coupled early diagenetic processes in the surface seaﬂoor (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) sea-ﬂoor sediment. In the vertical, the sediment is subdivided into two different zones. Solids (biogenic, mineral, etc.) raining down from the surface of the ocean are collected by the reactive mixed layer at the top. This is where chemical reactions take place. Solids are transported by bioturbation and advection, solutes by diffusion and bioirrigation. The classical coupled time-dependent early diagenesis equations (advection-diffusion reaction equations) are used to describe the evolutions of the 5 solid and solute components here. Solids that get transported deeper than the bottom boundary of the reactive mixed layer enter the second zone underneath, where reactions and mixing are neglected. Gradually as solid material gets transferred here from the overlying reactive layer, it is buried and preserved in a stack of layers that make up a synthetic sediment core. MEDUSA has been extensively modiﬁed since its ﬁrst release from 2007. The composition of the two phases, the processes (chemical reactions) and chemical equilibria between solutes are not ﬁxed any more, but get assembled from a set of XML 10 based description ﬁles that are processed by a code generator to produce the required Fortran code. 1D, 2D and 2D × 2D interfaces have been introduced to facilitate the coupling to common grid conﬁgurations and material compositions used in biogeochemical models. MEDUSA can also be run in parallel computing environments using the Message Passing Interface (MPI).

1 Introduction 1.1 Ocean-atmosphere :::::::::::::: Ocean-sediment exchange schemes: an overview Ocean biogeochemical cycle models call upon a variety of schemes of different complexity levels to represent ocean-sediment exchange fluxes.These can be classified into four major categories (Hülse et al., 2017): (1) reflective boundary conditions; (2) semi-reflective or conservative; (3) vertically-integrated dynamic models; (4) vertically resolved diagenetic models.These categories are similar but not identical to the levels in the classification of Soetaert et al. (2000): categories 3 and 4 respectively correspond to their level 3 and 4 descriptions; category 1 fits their level 2 description, while category 2 generalises this ::: the latter.
Reflective boundary conditions are the simplest of these: material reaching the model sea-floor (i.e., the deepest layer in the model water column) is unconditionally remineralised (oxidized, dissolved) there.Global mass conservation is obviously guaranteed with this approach, but the approach may be unrealistic in some places: calcite gets dissolved even if the sea-floor bathes in waters that are strongly supersaturated with respect to calcite or organic matter oxidized even if oxygen levels are low.This unrealistic behaviour can to some extent be alleviated with the semi-reflective or conservative scheme.Here, only part of the remineralisation products (nutrients, dissolved inorganic carbon, silica, etc.) of the solids reaching the sea-floor are returned to the bottom water; the remainder is returned to the surface, mimicking riverine input, and once again allowing for global mass conservation.The fraction remineralised can be made to vary in space and time and can also be different for different materials.Carbonate fractions remineralised can, e.g., be linked to the degree of saturation with respect to one carbonate mineral or another, and organic carbon fractions remineralised to the degree of oxygenation of the bottom waters.
Both schemes are attractive because of their convenient computational efficiencies.They do, however, not allow to take into account the complexities of the actual remineralisation pathways of the various biogenic components in the surface sediments, nor can they represent the temporary storage of such material in the surface sediment and delayed return of nutrients, dissolved inorganic carbon or silica to the ocean bottom waters.
The vertically-integrated dynamic category 3 encompasses ocean-sediment exchange schemes that explicitly include a single-box representation of the surface sediment.Mass balances of some, if not all, constituents of this single-layer sediment can be traced.Although termed "vertically-integrated" not all of the schemes that fall into this category can be traced back to some actual vertically resolved model that was vertically integrated.
Vertically resolved diagenetic models finally represent the most mechanistically oriented alternative to represent oceansediment exchange fluxes.Such models can take into account the complex interplay between various diagenetic processes (organic matter remineralisation, mineral dissolution or precipitation) and transport pathways (advection, bioturbation, solute diffusion in porewater, bioirrigation, etc.).They solve a set of coupled standard early diagenesis equations (Boudreau, 1997) for solid and dissolved component concentrations, generally in combination with law of mass action equations for chemical equilibria (e.g., for the carbonate system).
The actually required complexity of an adopted ocean-sediment exchange scheme depends of course on the application made.For short-term experiments (say a few decades to a few centuries) with high-resolution biogeochemical models carefully calibrated semi-reflective/conservative schemes are generally the only viable, but nevertheless perfectly acceptable option.
For long-term applications (simulation experiments exceeding several thousands of years, i.e., several ocean mixing cycles), vertically integrated or resolved schemes are required for realistic model responses to changing boundary conditions and forcings.
The surface sedimentary mixed-layer, where most of the processing of the deposited biogeochemically relevant material takes place only extends down to about 5 to 15 cm on global average (9.8 ± 4.5 cm according to Boudreau (1994)).As a result of the diagenetic processes in action there, strong concentration gradients are generated and sustained: the amplitude of the concentration differences observed over this depth interval may be comparable to those observed in the more than four orders of magnitude thicker overlying water column (3750 m on average) -see, e.g., the oxygen and pH profile data in section 3.3 below.
A realistic explicit representation of the surface sedimentary environment thus requires a vertical resolution of the same order of magnitude in terms of vertical layers or grid points than the complete water column above it.Typical vertically resolved early diagenesis present vertical resolutions of the order of ten to twenty layers (see Tab. 1).For comparison the water column in GENIE-1, which includes SEDGEM, (Ridgwell, 2007;Ridgwell and Hargreaves, 2007) has eight vertical layers, HAMOCC-2S (Heinze et al., 2009) has ten and the more recent HAMOCC 5.2 (Ilyina et al., 2013a) has forty layers.
Accordingly, sea-floor sediment modules (category 3 and 4 schemes) are not yet commonplace in global ocean biogeochemical models.Only four out of fifteen Earth System Models of Intermediate Complexity (EMICs) participating in the EMIC AR5 Intercomparison Project (Eby et al., 2013) are reported to have a sediment module included: Bern3D (Tschumi et al., 2011), DCESS-ESM (Shaffer et al., 2008), GENIE (Ridgwell and Hargreaves, 2007) and UVic 2.9 (Eby et al., 2009); one further participating EMIC, CLIMBER 2, is also routinely used with a sedimentary module included (e.g., Brovkin et al., 2012).Advanced high-resolution models generally call upon category 1 or 2 schemes for their sea-floor boundary condition, although there are exceptions.The HAMOCC (HAmburg Model of the Oceanic Carbon Cycle) family of models, whose origins reach back to Maier-Reimer (1984), actually has a long-standing history of explicitly taking sedimentary processes into account.HAMOCC 2 (Heinze et al., 1991) was the first one to get a fully coupled sediment module (Archer and Maier-Reimer, 1994), the oxic-only version of the calcite dissolution model of Archer (1991).Later, it received a purposely developed sediment module (Heinze et al., 1999;Heinze and Maier-Reimer, 1999).Archer et al. (2000) use HAMOCC 2 coupled to the much more complete diagenetic model MUDS (Archer et al., 2002), which considers the sequence of oxic, suboxic (via NO − 3 , FeOOH and MnO 2 reduction) and anoxic (via SO 2− 4 reduction) organic matter remineralisation pathways.Later developments of HAMOCC also included suboxic remineralisation pathways of organic matter in the standard sediment model of HAMOCC: in HAMOCC 5.1 (Maier-Reimer et al., 2005) denitrification was added, in HAMOCC 5.2 (Ilyina et al., 2013b) sulfate reduction.Gehlen et al. (2006) have coupled a slightly extended version of the sediment model of Heinze et al. (1999) to PISCES, the biogeochemical component of NEMO, with nitrate reduction and denitrification as an additional remineralisation pathway for organic matter.
The Gehlen et al. (2006) model was later also introduced as the sediment component into Bern3D (Tschumi et al., 2011).
Sea-floor sediments are not only relevant as processing " ::::::::: "processing : units" for biogenic material raining down from the surface euphotic layer, during which some part gets :::: parts ::: get remineralised (i.e., oxidised or dissolved), and some parts ::: the ::: rest : gets buried.Burial is, however, at first only temporary.Changes in the overlying boundary conditions (e.g., saturation conditions) may indeed lead to chemical erosion episodes where the surface sedimentary mixed-layer loses material faster than it is replenished by deposition from the water-column above.We are currently at the onset of such an episode: as ocean acidification due to the uptake of anthropogenic CO 2 progresses to the deep-ocean, the resulting change in the degree of saturation with respect to carbonate minerals is expected to enhance the dissolution of carbonates in the sea-floor surface sediments at depth so strongly that the dissolution rate will exceed the rate at which carbonate material gets deposited at the sediment water interface :::::::::::::::: (Archer et al., 1998).Previously buried carbonates will then return to the sedimentary mixed :::: layer as a result of the bioturbation activity, which tends to keep the surface mixed layer at a rather stable thickness, which seems to be controlled by the supply of organic matter (Boudreau, 1998).Archer (1996b) estimates that existing carbonates in surface sea-floor sediment can neutralize about 1600 GtC, considerably more than the ∼1000 GtC that may at most be emitted while still limiting global anthropogenic temperature change below 2 • C (e.g., Zickfeld et al., 2016), but much less than the estimated resources of fossil fuels of 8,543 -13,649 GtC (Bruckner et al., 2014, Tab. 7.2).
Finally it should not be forgotten that sea-floor sediments represent our most comprehensive source of information about past climate change and it is of course indispensable to understand how early diagenetic processes influence the sedimentary record, and it would be desirable to directly compare generated model (synthetic) sedimentary records to the observed records instead of thus opening new possibilities in terms of data assimilation.

MEDUSA: from version 1 to version 2
The first version of the Model of Early Diagenesis in the Upper Sediment,1 MEDUSA -hereafter MEDUSA-v1 -was described in Munhoven (2007).It is a time-dependent vertically resolved biogeochemical model of the early diagenesis processes in the seafloor ::::::: sea-floor : sediment.MEDUSA-v1 included clay, calcite, aragonite and organic matter as solid components, and CO 2 , HCO − 3 , CO 2− 3 and O 2 as porewater solutes.Besides that configuration, two others (unpublished) had been developed: one which furthermore included opal and dissolved silica and one which also included the 13 C isotopic signatures of all carbonbearing components.
Right from the beginning, MEDUSA had been developed as a sediment module for the diverse ocean biogeochemical models used in our research group, ranging from box to three-dimensional models, the latter with diverse grid configurations and also various sets of chemical tracers.Furthermore, our research interests required a model that could be used in studies dealing with time scales ranging from tens of years to hundreds of thousands of years.Accordingly, the model requirements were laid out as follows: (1) the model code should be customizable to accommodate different chemical compositions; (2) the model should offer the possibility to be coupled to strongly different host model grid layouts; (3) it should be possible to run the model with variable time-steps; (4) the model must be able to cope with chemical erosion, i.e., be able to recover previously buried material from deeper layers and to return it to the chemically reactive mixed-layer.
The customization options offered by MEDUSA-v1 had to be selected with pre-processor directives in the code.Extending the capabilities of the model on the basis of that mechanism had become more and more cumbersome and difficult to manage with time.The code was therefore revised in depth and only the parts related to the transport terms in the equations and the equation system solver-the framework system-were kept.The rest of the code was from now on built on purpose for each application with a code configuration and generation tool that would produce and assemble the parts related to the components, processes (reactions) and chemical equilibria required.A code generator was developed to read in the required information, such as chemical and physical properties of components, chemical reactions describing diagenetic processes and chemical equilibria from a series of description files.These description files use a format based upon the eXtensible Markup Language (XML) syntax (W3C, 2008).Organising the information in an XML tree offers attractive flexibility: such a tree can be easily extended for later developments and it is possible to access any particular information wherever it is located in a file.XML thus offers a high degree of compatibility between subsequent versions of the configurator, which can always extract the relevant information as long as the required mark-up tags remain present.Above all, XML files remain mostly human-readable and the possibility to insert comments makes it possible to ensure the traceability of the stored information.
As the complete tool was meant to require only a Fortran compiler to be built, a simple library, called µXML, for reading and processing basic ASCII encoded XML files in Fortran 95 was developed as a prerequisite.

Vertical partitioning of the sediment column
The complete sediment column is subdivided into three (or four) different vertically stacked parts (called realms), as illustrated on Fig. 1: (1) REACLAY, the top-most part extending downwards from the sediment top at the sediment-water interface and where the chemical reactions are taken into consideration; (2) TRANLAY, the transition layer of changing thickness just underneath, acting as a temporary storage to connect REACLAY to the underlying (3) CORELAY, a stack of sedimentary layers representing the deep sediment, i. e., the sediment core.Additionally an optional diffusive boundary layer (DBL-not to scale on Fig. 1) acting as a diffusive barrier to the sediment-water exchange of solutes can be included on top of the REACLAY realm.
REACLAY includes the bioturbated sedimentary mixed-layer, where most of the reactions relevant for early diagenesis take place (organic matter remineralisation, carbonate dissolution etc.).

Equations in the DBL and REACLAY realms
In the REACLAY realm, MEDUSA solves the standard time-dependent diagenetic equation (e.g., Berner, 1980;Boudreau, 1997), which can be written for each sediment component i (solid or solute) in generic form as 200 In this equation, t is time and z depth below the sediment-water interface ::: SWI : (positive downwards -see Fig. 1).Ĉi denotes the concentration of i in moles for solutes and in kg for solids per unit volume of total sediment (solids plus porewater).Ĵi is the local transport (advection and diffusion), per unit surface area of total sediment.Ŝi = Ri + ri + Qi represents the net source-minus-sink balance for constituent i per unit volume of total sediment, where Ri is the net reaction rate, equal to the difference between production and destruction (or decay) rates, ri is the net fast reaction rate, that is going to be filtered out of 205 the equations by equilibrium considerations and Qi the non-local transport (considered only for solutes).The total sediment concentrations Ĉi are related to the more directly accessible phase-specific concentrations C s i (for solids) and C f i (for solutes) by Ĉi = ϕ s C s i and Ĉi = ϕ f C f i , respectively.ϕ s and ϕ f denote the volume fraction of bulk solids and of porewater in the total sediment, linked to porosity ϕ by ϕ s (z) = 1 − ϕ(z) and ϕ f (z) = ϕ(z).The porosity profile ϕ(z) is prescribed but may be different for each column in multi-column set-ups.
In the DBL (if any), only equations for solutes are considered and porosity is set to 1. Solids are supposed to rain through the DBL and directly enter REACLAY at its top surface.

Transport
Solids are transported by advection throughout the sediment column and subject to bioturbation in the surface mixed layer.
Bioirrigation provides a non-local transport mode for solutes.In MEDUSA, the source-sink approach (Boudreau, 1984) is used to quantify the effect of bioirrigation: Here α is the bioirrigation "constant", which may be depth-dependent, and C oc i the concentration of solute i in the irrigation channels, set equal to the solute's concentration in the seawater overlying the sediment.

Chemical Reactions and Equilibria
The set of chemical reactions and equilibria to consider is completely dependent upon the given application, i.e., on the chemical composition of the sediments required, the diagenetic processes to consider (e.g., organic matter remineralisation, possibly following several pathways, carbonate dissolution, ) and the equilibria between components of solute systems (e.g., carbonate, phosphate, borate systems) to take into account.MEDUSA does not include a standard composition and reaction/equilibrium network but must be configured to fit the complexity requirements of a given application: including one or more classes of organic matter (solid or dissolved), one or more types of carbonates, one or more organic matter degradation processes, etc.
The chemical interconversion reactions represented by the ri terms in the source-minus-sink term Ŝi are orders of magnitudes faster than all other reactions.They are supposed to evolve in quasi-equilibrium.The ri are therefore eliminated from the partial differential equation system by considering appropriate linear combinations of selected equations and by including the thermodynamic equilibrium equations in the system of equations.The partial differential equation system is thus converted into a differential algebraic equation system.The subroutines required to evaluate the source and sink terms related to chemical reactions and to convert the complete system to a differential algebraic system are generated by the companion MEDUSA COnfiguration and COde GENeration tool (MEDUSACOCOGEN) described in section 2.4 below.
The static volume conservation equation is used to replace one of the solids' evolution equations: Here I s denotes the inventory of solid components considered in the model configuration and ϑ i the partial specific volume of solid i.We suppose that the densities ρ i of individual solid components are constant and independent of each other.In this case ϑ i = 1/ρ i and the ϑ i are also constant.However, although the density of each solid constituent is constant, the average density of the solid phase may vary both in space and time as chemical reactions proceed and modify the sediment composition and influence the advection rate profiles.To ensure compatibility with other early diagenesis models that assume that the solid phase has a constant density (both in space and time), all reactive (non inert) solid components can optionally be declared as volumeless.

Boundary conditions
The differential equation systems describing the evolutions of the compositions of the DBL and REACLAY and the TRANLAY realms have to be completed by boundary conditions connecting them to the ovelying :::::::: overlying : seawater, the underlying CORELAY :::::::: TRANLAY and between each other.Solute concentrations at the top are derived from those in the overlying ocean At the bottom of the REACLAY realm, flux continuity is adopted for both solutes and solids.For solutes, this requires the concentration gradient to reduce to zero (Neumann boundary condition) since u(z) ≡ 0 is adopted here.For solids, a variety of effective boundary conditions arise -and the types may change in time -depending on whether biodiffusion has vanished there or not and whether the sediment is burying (w + B > 0, where w + B denotes the advection rate on the outer side of REACLAY's bottom interface) or eroding (wz + B < 0 :::::: w + B < 0) solid material at the bottom of REACLAY.

TRANLAY and CORELAY
TRANLAY is a TRANsition LAYer that collects the solids leaving REACLAY through the bottom (see Fig. 1).As soon as its thickness exceeds a given threshold value (by default 1 cm) at the end of a time step, one or more new sediment core layers are formed and subtracted from TRANLAY to be transferred to CORELAY, which is managed as a last-in-first-out stack of sediment layers.
In general, material is only preserved in TRANLAY and CORELAY; it is assumed that no chemical reactions take place there.
However, we have to make one exception to this rule: reactions that are part of a radioactive decay chain are still taken into account in these two realms to avoid physically unrealistic results.Radioactive material that would have left REACLAY and be returned there later during a chemical erosion event would contribute to create unrealistically young concentrations in REACLAY if radioactive decay would have been temporarily suspended for a more or less extended time.
The bottom point of the grid that covers REACLAY is always a node.The nature of the topmost point depends on whether a DBL is included or not: if no DBL is included, the topmost interface of REACLAY-the sediment-water interface, SWI-is ::::::: SWI-is located at a grid node; if a DBL is included, the interface at the top of REACLAY is mapped onto a grid vertex, defined by a virtual node located above the top of REACLAY and the topmost interior node of REACLAY.Similarly, the bottom of the DBL is located at a vertex of the DBL grid, defined by a virtual node located below the bottom of the DBL and the lowest node inside the DBL (which is most often the top of the DBL).Detailed information about the grid generation can be found in the technical report ' : "MEDUSA Technical Reference' : " in the Supplement.
In multi-column setups-this would be the most common usage for coupling to biogeochemical cycle models-every sediment column in MEDUSA must have the same number of grid-points, but the spacing and extent of each of these may be different.

Spatial discretisation
::::::: Discrete ::::::::: equations The discrete version of the evolution equation for a component i in a cell j represented by the node situated at z j writes Here, ::::: where ∆t = t n − t n−1 :: is ::: the :::::: implicit :::: time :::: step : and h j is the distance that separates :::::: between : the vertices z j− 1 2 = 1 2 (z j + z j−1 ) and z j+ 1 2 = 1 2 (z j+1 +z j ) ::: that ::::::: delimit :: the :::: cell : j.Equation ( 6) is slightly modified at the bottom-most node where it relates to a half-cell only and one may chose to formally express the mass-balance equations for that half-cell at the representative node (which is actually the bottom node of the grid, resp.), or at some intermediate point between that node and the delimiting vertex.In the absence of a DBL, a similar procedure is adopted at the top-most node.
The numerical schemes adopted in MEDUSA have been selected with the physical meaningfulness of the results in mind.Accordingly, positiveness of the calculated concentration evolutions was deemed indispensable.The discretisation of the advective part of the local transport term in the equations thus requires an upwinding approach.One may chose between a first-order full upwind and a second order exponential fitting scheme, known elsewhere as the Allen-Southwell-Il'in or the Scharfetter-Gummel scheme (Hundsdorfer and Verwer, 2003).It is closely related to the scheme of Fiadeiro and Veronis (1977): on regularly spaced grids both schemes lead to identical discrete forms of the equations.The exponential fitting scheme is, however, better suitable for the flux-conservative finite volume approach on irregularly spaced grids adopted in MEDUSA as it allows for exact mass conservation.
Unlike steady-state models, where the solids' advection rate is always oriented downwards relative to the sediment-water interface, MEDUSA has to be able to cope with solids' advection rates that may have any orientation and that may even change their orientation with time.Both upwinding schemes automatically handle this complication.

Code Organisation
The MEDUSA common framework includes the subroutines to assemble the equation system and its Jacobian, and to solve the equation system, modules to make available fundamental data (physical constants, unit conversion parameters, etc.), modules to hold the forcing data and the intermediate results.It also provides the core management system for multiple sediment columns, which can be processed in a sequential or a parallel fashion.For parallel processing, Message Passing Interface (MPI) calls are included and can be activated by a pre-processor switch.Three different MPI interfaces are provided: 1D, 2D and 2D×2D, respectively for a sequential linear distribution of the sediment columns, a two-dimensional ordering (typically longitude-latitude) and a hierarchically ordered two-dimensional array of two-dimensional arrays of sediment columns.
In multi-column setups the chemical composition (solids and porewater solutes) must be the same in all the columns.It is nevertheless possible for each organic matter type (solid or solute) to have different C:N:P:O:H ratios in each column.These individual ratios must, however, remain constant with time.
The framework system must be completed with the specific parts required for a particular application to build a working instance of MEDUSA.This includes first of all the requested composition in terms of solids and solutes, the reaction network of the diagenetic processes to consider and the chemical equilibria between components of solute systems.MEDUSA must also be aware of the material characteristics such as densities, molar compositions, and some thermodynamic properties for solid components, or diffusion coefficients for solutes.Rate laws for the different processes under consideration need to be specified.
The information that is required for producing the Fortran code is collected in a series of XML files that define the composition of the sediment and describe the components, processes and equilibria to be considered.These are processed by the MEDUSA COnfigurator and COde GENerator, MEDUSACOCOGEN, to generate as diverse things as modules providing index parameters to address single components by meaningful names, subroutines to calculate molar masses of the components, subroutines to evaluate reaction rates of all the components and the corresponding derivatives with respect to the relevant component concentrations, to modules to assure the I/O to NETCDF files.The complete diagenesis model is bundled into an object library (libmedusa.a) to be linked with the application (host model).

MEDUSACOCOGEN: the MEDUSA COnfigurator and COde GENerator
The : "Reference Guide to the Configuration and Code Generation Tool MEDUSACOCOGEN : " : in the Supplement provides an exhaustive description of the procedure to follow to build a working MEDUSA application.The formats of the required files are described in full detail with commented examples.The library of rate law functions and equilibrium relationships is also presented in detail.Further information can be found in the example applications provided in the code.Here, we present only a general overview about the functionality of the code generator.

Main Building List
The main building list provides the names of the description files of the solids, solutes and solute systems to consider in a particular model configuration, as well as the process and the equilibrium description file names.

Sediment components: solids, solutes and solute systems
Components (solutes or solids) can be of several types and be part of different classes.We distinguish three types of components: ignored (default), normal or parameterized (for solutes only -see below).Evolution equations are only generated for normal components.
Solids actually encompass all the characteristics of solid phase components: their concentrations, their age or production time, their isotopic signatures, . . .A solid's description file provides information about its physical and chemical properties, such as intrinsic density, alkalinity content per mole (both mandatory), chemical composition, molar mass (optional).A solid can be part of one of four classes: basic solid (default), (particulate) organic matter, solid colour or solid production time.The basic solid class includes all physical solids but organic matter, be they reactive or not.For numerical stability reasons, it is mandatory to include at least one inert solid in the model sediment components, to be flagged as mud.There is a dedicated class for organic matter, offering special functionality.Chemical composition is mandatory for this class and can be set in terms of their C:N:P ratios, from which the actual composition is then derived by CH 2 O, NH As mentioned above, there is a special type of solutes: parametrized.For parametrized solutes, no evolution equations are generated.Their description files must therefore include a code snippet to calculate their abundance or to derive their value from specific boundary conditions.Typical examples of parametrized constituents are the calcium concentration, which can be derived from salinity, or the saturation concentration of CO 2− 3 with respect to calcite, which can be calculated from the degree of saturation at the boundary or from the solubility product.The description files of normal solutes must include a code snippet to calculate its diffusion coefficient in free seawater.Solutes' descriptions must also include their specific alkalinity content per mole.We only distinguish between two classes of solutes: basic solute (default) and (dissolved) organic matter.Just like particulate organic matter, dissolved inorganic :::::: organic : matter must be characterised by its chemical composition in terms of its C:N:P(:O:H) ratios.For basic solutes, the chemical composition is optional.
Finally, a solute system is a set of solutes that MEDUSA considers as a total sum in the equilibrium calculations.Typical examples of solute systems are DIC (dissolved inorganic carbon, composed of CO 2 , HCO − 3 and CO 2− 3 ) or the borate system (B(OH) 3 and B(OH) − 4 ).Solutes that are part of a solute system are considered to be in local chemical equilibrium with each other.Solute system description files simply list the description files of the solute components that make them up.The code generator internally compiles Total Alkalinity as a solute system, based upon the specific alkalinity contents declared in the solute description files included in the model configuration.

Processes and equilibria
Description files for processes include a representation of the underlying chemical reaction that translates its effect on the various model constituents and specify the rate law to apply.MEDUSACOCOGEN currently provides twenty-one different rate law formulations (actually thirty as several of them have a few variants).Similarly, equilibrium description files must include a representation of the chemical equilibrium and specify the law of mass action to use for it.Expressions for laws of mass action (equilibrium relationships) require less variety than process rate laws: the four provided library routines cover the most common cases.
Rate laws and laws of mass action are provided in specially formatted Fortran 95 modules (so-called MODLIB files).The source code of a MODLIB file is preceded by a header (protected by a conditional inclusion pre-compiler directive) that provides meta-data helping to identify and classify the different parameters required.The module itself defines (1) a derived type structure that encapsulates the parameter values and relevant component index references and (2) a subroutine to evaluate the rate-law expression, resp.the equilibrium relationship, for given concentration and parameter values, and the derivatives with respect to the concentrations of the model components.MODLIB files for laws of mass action must further include a subroutine to set the equilibrium constant (currently derived from the boundary conditions) and another one to evaluate the scale to applied to the equilibrium relationship.The current collection can be easily extended by adding MODLIB files to the library (for details about how to do this, please refer to the MEDUSACOCOGEN reference guide in the Supplement).
Chemical equilibria are always taken into account in the DBL and the REACLAY realms.Processes are usually considered only in the REACLAY realm.However, processes involving only solutes can also be considered in the DBL.In addition, processes that represent radioactive decay of solid trace elements (i.e., of elements whose volume is considered negligible and that do therefore not impinge on the advection rate profile) should be declared to apply also in the TRANLAY and CORELAY realms (where reactions are normally stalled).This way, adequate corrections can be applied in case a sediment column becomes subject to chemical erosion and previously buried material gets remixed into the REACLAY realm.Without such corrections, the material returned to TRANLAY or REACLAY would appear too young.

Test Case Applications
Three applications have been selected to illustrate the functionality of MEDUSA: (1) a replication of the ALL simulation experiment from Munhoven (2007), supplemented with an extended configuration that attaches an age tracer to calcite; (2) a coupling simulator where the boundary conditions that would normally be provided by a biogeochemical model that MEDUSA would be coupled to are read from a file; (3) a MEDUSA configuration with complex composition and reaction network as considered in state-of-the-art early diagenesis models used for the analysis of site observational data.:: All ::: of ::: the :::: test :::: case :::::::::: applications ::: use ::::: α ≡ 0 ::: (no :::::::::::: bioirrigation).

MEDMBM-PT -Coupled simulation experiment with chemical erosion and resolved sedimentary records
MEDUSA produces truly resolved and, via its CART-based age/production time control system, consistently dated synthetic sedimentary records.In order to illustrate the potential of this combination of features, the ALL experiment from Munhoven ( 2007) is repeated here with a MEDUSA configuration equivalent to the one used in that study, extended to attach age control information to calcite particles.

Application description
Munhoven (2007) used MBM (Multi-Box-Model) coupled to MEDUSA-v1 to explore the implications of the rain ratio changes proposed to explain the glacial-interglacial atmospheric CO 2 variations (see, e.g., Archer and Maier-Reimer, 1994)   represent the evolutions of the carbonate compensation depth (CCD, defined here as the depth where the carbonate dissolution rate is equal to the carbonate deposition rate) and of the calcite saturation horizon (CSH, i.e., the depth at which saturation with respect to calcite is reached).The two distributions are broadly similar but present noticeable differences.There is first of all a systematic time lag of about 3 to 8 kyr, best discernible in the regions where %CaCO 3 exceeds 70%.This time lag is due to the fact that the age tracer tracks the average age of the surface sediment at burial, which is different from zero at all times.
Another important difference is the alteration of the actual evolution (Fig. 2a) in the sedimentary record (Fig. 2b).The most striking differences are visible during times of shoaling CCD at depths where %CaCO 3 is lower than 70%.As shown by the long-dashed line on Fig. 2b, which traces the evolution of the limit between the magenta and the black zones from Fig. 2a (the 0% calcite line), shifted by 5 kyr towards the greater ages to account for the average calcite burial age, up to 20 kyr of calcite history have been deleted from the record between 4700 and 6700 m d.b.s.l.due to chemical erosion.
MEDUSA always includes at least internally a TRANLAY-CORELAY stack of sediment layers underneath the REACLAY realm in order to handle possible chemical erosion events.These layers can a priori be only approximately dated from the recorded time of burial for each layer, i.e., the time when a given layer is separated from TRANLAY and transferred to CORELAY.With the CART based production time information attached to a solid constituent, this shortcoming can be overcome.

COUPSIM -Coupling simulator
MEDUSA has already been coupled to several ocean biogeochemistry and Earth System Models.First results have been published for the coupling to the Community Earth System Model, CESM (Kurahashi-Nakamura et al., 2020); other coupling projects are well advanced (Moreira Martinez et al., 2016;Völker et al., 2020).
To illustrate the procedure of coupling MEDUSA to a biogeochemical model and to assess the computational cost of a typical sediment module for a real three-dimensional ocean biogeochemistry model, a coupling simulator, COUPSIM, was developed where the boundary conditions that would normally be provided by the host biogeochemical model are read from a file.

Application description
We use results obtained with the coupled Biogeochemistry/Ecosystem/Circulation model, BEC (Moore et al., 2004),3 that were made publicly available (Moore et al., 2005)  The resulting global distributions of the solid fractions of calcite, total organic carbon and opal (resp.denoted %Calcite, %TOC and %Opal hereafter) were then evaluated against observational data (Seiter et al., 2004) on the basis of their respective standard deviations and correlation coefficients.The results of these experiments are summarised in the Taylor diagram (Taylor, 2001) on Fig. 3.For that diagram each one of the three distributions was normalized with respect to the standard deviation of its respective observational counterpart in order to be able to report all the results on a common scale.The peculiar distributions of the different characteristic points on that diagram indicate that there is an structural incompatibility between the data and the possible model results.Calcite points cluster around a correlation coefficient of 0.38 and a standard deviation of 1.3, despite a large range of values (5-1000 %day −1 ) used for the dissolution rate constant (k c in Table 4) in the experiments.The opal and organic carbon points align on two beams with correlation coefficient values between 0.23 and 0.27, and between 0.35 and 0.45, respectively, each one with a large range of standard-deviations.For the organic matter dissolution rate laws values between 0.005 and 0.032 yr −1 were adopted independently for the two rate constants (k ox and k nr in Table 4); for opal dissolution, rate constant values between 0.04 and 0.07 yr −1 , together with asymptotic concentrations ranging between 500 and 700 µmol/L were used (resp.k o and C os in Table 4).
Among all the experiments carried out, the one that offered the best compromise in terms of standard deviations coming as close as possible to the observations (i.e., to 1 after the normalization) and maximizing the correlation was selected as the best- ).: For opal, a comparatively high value of 0.05 yr −1 (mol m −3 ) −1 had to be adopted in order to avoid widespread opal-dominated sea-floor sediments.For comparison, the rate constant of 30 yr −1 of Boudreau (1990) translates to 0.0034 yr −1 (mol m −3 ) −1 for the rate-law formulation adopted here.The asymptotic ("saturation") concentration for opal dissolution had to be set to 700 µM, which ranges at the cold-water end (≤ 1 • C of the experimentally derived values of Dixit et al. (2001)).
The resulting average surface model sediment composition at steady state is compared to the target data of Seiter et al.
( the maximum in the equatorial East Pacific.This latter is, however, too narrow and similarly to the calcite maxima too sharply delimited. Among the structural incompatibilities, the band of calcite-rich sediments along the entire rim of the Pacific Ocean in the Northern Hemisphere stands out.These are not seen in the data for the simple reason that the CSH is actually shallower than 1000 m almost everywhere along this band, except in the Sea of Japan (or East Sea) (Yool et al., 2001).In the BEC results, this whole band is supersaturated with respect to calcite.No parameter combination can override this supersaturation and significantly reduce the amount of carbonate preserved.The reflective boundary condition actually contributes to this unrealistic supersaturation, possibly even causes it: as all of the calcite that reaches the deepest ocean cells there unconditionally dissolves the degree of super-saturation of the bottom waters is artificially increased.In the coupled CESM-MEDUSA model experiment of Kurahashi-Nakamura et al. (2020), this carbonate-rich band is not produced.Another important feature of the model %Calcite is the absence of carbonates at the seafloor ::::::: sea-floor in the Atlantic north of the equator and south of 30 • N, again in contradiction with the data.This results from very low CaCO 3 deposition rates of the order of 1-3 g/m 2 /yr, combined with high lithogenic deposition rates of the order of 15 g/m 2 /yr in the BEC results, which caps %Calcite at about 5-15% a priori even in the absence of dissolution.The extended areas of carbonate-rich sea-floor sediments along the North Pacific rim and of strongly diluted sediments in the central Atlantic contribute to the poor correlation coefficient, whatever the model parameter combinations chosen.
In general, the model sediment compositions show far more pronounced contrasts than the observed ones.In all three distributions, the more diffuse features of the global distributions seen in the data are missing.Intermediate values are widespread in the observed distributions, especially for organic carbon where there are only a few isolated spots at the upper end of the depicted range.The model produces in most regions surface sediments that are either poor or extremely rich in calcite.There are only few areas with calcite fractions between 20 and 80 %.The same holds for organic carbon and opal, albeit to a lesser degree for the latter which presents more extended areas with intermediate abundances.The sharp contrasts for %Calcite can partly be explained by to the vertical grid resolution of the BEC version used by Moore et al. (2004).The CSH in the South Pacific is situated at a depth of about 2.5-3.3 km and at about 3-3.7 km in the Indian Ocean (Yool et al., 2001).At these depths, the vertical resolution of BEC is about 340-450 m.As the transition zone from carbonaterich to carbonate-poor sediment depths is about 500 m thick there, it is clear that this transition zone cannot be satisfactorily option would be to call upon sub-grid scale depth profiles and to attach several (two to four) sediment columns to each model sea-floor grid element.This would of course increase the computational burden of the sediment part of the model, but as can be deduced from the execution times reported and discussed below, these would not represent a major hindrance, given the efficiency of the numerical procedures adopted in MEDUSA.
The %Opal distribution is essentially correlated to the deposition flux rates of opal, which themselves present a highly contrasted distribution.%TOC appears to depend on the organic deposition flux rate and the bottom water oxygen distribution.
The preservation tongue in the equatorial East Pacific results from a well delimited deposition rate maximum, in combination with mid to low oxygen concentration rates ::::::::::: concentrations : (20-160 µmol/L).The high organic deposition rates in the Southern Ocean and along the west coasts of Africa do, on the other hand, not lead to high %TOC as these regions are much better oxygenated (160-250 µmol/L).

Execution times
The execution times for the best-fit experiment for different computing environments (serial, parallel MPI with the MPI-1D and MPI-2D interfaces) are reported in Table 5. MPI-2DT2D execution times are not reported as they are always within a few percent of those for MPI-2D.The results discussed above were for steady state, i.e., one infinitely long time step.To assess the behaviour can be simulated in MEDUSA by calling upon the volumeless solids option at compile time.With this option, only the main inert solid is considered to have a finite density.All the others are considered to be tracers without significant volume, but only mass (i.e., they have zero partial specific volume and infinite density).Both approaches are mathematically equivalent and both have an important drawback: they allow to transport physically unrealistic amounts of non-inert material (e. g., calcite mass fractions exceeding 100%), as we will see below.The model sediment columns were assumed to extend down to 82 cm as in the original study; we used a REACLAY grid with 321 nodes, without a DBL.Porosity and bioturbation profiles were prescribed, using the information provided by Jourabchi et al. (2008).Boundary and forcing conditions (solute concentrations, calcite saturation state at the SWI, solids' burial rate) were also taken from Jourabchi et al. (2008).For our purpose, the solids' burial rate was converted to an equivalent flux of inert material (dubbed clay) across the SWI.

Adjustment procedures
The adjustment was carried out in two stages, for each of the thirteen stations.At the first stage, the organic matter deposition rate and the degradation rate constants were adjusted in order to reproduce the measured oxygen concentration profiles; calcite dissolution was ignored.For the second stage, these two parameters were held fixed at the values obtained during the first stage.
This time, the calcite deposition rate and the calcite dissolution rate constant were adjusted in order to reproduce the measured pH profiles.We only considered calcite dissolution rate orders 1, 2 and 4.5, disregarding the relatively uncommon order 0.5.
In contrast to Jourabchi et al. (2008), we did not adopt static initial values for the calcite dissolution rate constant at the second step.We tested random perturbations to these in order : to : allow a deeper exploration of the results space.In a few instances, this allowed us to significantly improve on the results of Jourabchi et al. (2008).Fits were carried out with MPFIT (Markwardt, 2009) 5 with the GNU Data Language (GDL, v. 0.9.6).The minimisation procedure implemented in MPFIT is based upon the Levenberg-Marquardt algorithm.For the pH profile fitting, calculations were based upon H + concentrations derived from the pH data.

Results and discussion
O 2 data are comparatively straightforward to fit.The model profiles obtained here are essentially identical to those of Jourabchi et al. (2008) and we are therefore not discussing them here.The detailed results can be found in the jeasim_vl.odsspreadsheet file in the work/jeasim directory of the code and data archive provided in the Supplement, on the sheet named jeasim_adj0; the resulting O 2 profile fits are shown on Figs.S1 and S2 in the : "Additional Results : " in the Supplement.In terms of the relative misfits, the quality of the fits obtained here is even in most instances slightly superior to those of Jourabchi et al. (2008), except for site 2, where ours is 37% worse.Three others (sites 19, 40 and 59) have similar relative misfits (within ±10% of those of Jourabchi et al. (2008); all the rest present between about 20 and 80% lower relative misfits.
pH profile adjustments, on the other hand, were far less straightforward.We therefore proceeded in several steps (up to three).First, the fitting procedure was carried out with only the pH profile data as a constraint.As mentioned above, the way the application is designed (with the volumeless solids option), physically unrealistic mass fractions exceeding 100% cannot be    no calcite dissolution; full -calcite dissolution kinetics of order 4.5; dotted -order 2; short dashes -linear.  .pH profiles with data-constrained surface sediment %calcite.Lines are as follows: long dashed -no calcite dissolution; fullcalcite dissolution kinetics of order 4.5; dotted -order 2; short dashes -order 1 (linear).Please notice that the n = 4.5 fit at site 2 and all the fits for sites 19 and 20 require calcite mass fractions in the deposition fluxes exceeding 100% (between 140 and 730%).
We concur with Jourabchi et al. (2008) in concluding that there are probably other constituents, reactions and processes not considered in our models that may significantly impinge on the proton exchange reactions and thus influence the pH profiles in a way that can thus not be reproduced.However, there are additional assumptions that must also be considered to explain the model-data discrepancies.The empirical relationship used to derive the burial rate from the sea-floor depth is possibly unreliable for usage at isolated sites.At some stations, a model configuration with a DBL might be required.This is almost certainly the case at site 120, for which Wenzhöfer et al. (2001) have made out a DBL thickness of 725 ± 25 µm.6

Conclusions
In its version 2 MEDUSA, the Model of Early Diagenesis in the Upper Sediment with Adaptable complexity, offers a flexible approach to develop tailored sediment modules suitable for coupling to ocean biogeochemical cycle models, but also for process analysis and teaching.MEDUSA is based upon the standard time-dependent diagenetic equation (e.g., Berner, 1980;Boudreau, 1997) with the vertical as the only spatial dimension.The chemical composition of the solid and solute phases can be chosen to fit the requirements of the scientific question to address.The network of diagenetic reactions (redox processes, mineral dissolution, . . . ) and chemical equilibria in the porewaters can be made exactly as complex as required.The Application Programming Interfaces (APIs) allow the coupling to the most commonly encountered geographical grid layouts in biogeochemical models (unstructured, two-dimensional, or hierarchically ordered two-dimensional arrays of two-dimensional areas).
They have furthermore been designed so that existing model configurations can be easily extended (or simplified).The adopted numerical procedures allow for large range of time steps, making MEDUSA suitable for long simulation experiments with wellresolved representations of the sea-floor (several thousands of sediment columns).Parallel processing for multi-column set-ups is supported via Message Passing Interface (MPI) instructions that can be optionally activated.
Three test case applications have been selected and presented in this model description paper.The first one revisits the study of Munhoven (2007) which used the predecessor version MEDUSA-v1.A model configuration with exactly the same functionality as MEDUSA-v1 was set up and extended by the age monitoring facilities offered by MEDUSA.These latter are based upon the Constituent-oriented Age and Residence time Theory, CART (Deleersnijder et al., 2001;Delhez et al., 1999), which has a well-established record in the analysis of geophysical fluid flow models.CART can be easily transposed to the standard equations describing the early diagenesis of seafloor ::::::: sea-floor : sediments.The newly generated model code accurately reproduces the original results.The added age-control information provides insight into the relationship between the actual time-dependent evolution of the sea-floor sedimentary mixed-layer composition and the resulting sedimentary record.
For the second application, a coupling simulator was developed to illustrate the opportunities and to assess the computational requirements of using a fully coupled vertically resolved early diagenesis model as the ocean-sediment exchange scheme in a model of ocean biogeochemical cycles.MEDUSA offers attractive flexibility requiring reasonable computational overburden.

Figure 1 .
Figure 1.Partitioning of the sediment column in MEDUSA: an optional diffusive boundary layer (DBL) on top of the main part of the model sediment where diagenetic reactions and advective-diffusive transport take place (REACLAY), the transition layer (TRANLAY) and the core represented by the stack of layers (CORELAY).The bottom of the bioturbation zone may coincide with the bottom of REACLAY.If the optional DBL is omitted, zW = zT.See text for further details.

Figure 2 .
Figure 2. Average surface sediment CaCO3 content for the ALL experiment from Munhoven (2007).(a) Actual surface sediment %CaCO3 as a function of time.The grey line traces the evolution of the calcite saturation horizon (CSH), the white line the carbonate compensation depth (CCD) where the dissolution rate equals the deposition rate.(b) Synthetic sediment core record from cores "drilled" at 100 m sea-floor depth intervals from 2050 to 6850 m d.b.s.l. in the low-latitude Indo-Pacific box, as a function of age of the layers.The long dashed line indicates the limit between the black and magenta zones in panel (a), shifted by 5 kyr to the right.

Figure 3 .
Figure 3. Taylor diagram of the results of the parameter sensitivity experiments carried out with COUPSIM: %Calcite (blue circles); %TOC (green diamonds); %Opal (blue ::: red).The full symbols indicate the combination found for the best-fit experiment.Maps for the surface sediment compositions of this best-fit experiment are shown on Fig. 4.
) on Fig. 4. Calcite-rich sediments are produced in the Indian and South Pacific; the calcite-rich sediments along the Atlantic mid-ocean ridge are only reproduced in the South Atlantic and at mid-latitudes in the North Atlantic.The sediments that are richest in TOC are located in the equatorial East Pacific.The opal belt in the Southern Ocean stands out, as well as resolved.There are several options to alleviate this shortcoming.The most obvious one is of course to increase the vertical resolution of the host model.The CESM version used by Kurahashi-Nakamura et al. (2020) has sixty vertical levels.Another

Figure 5 .
Figure5.pH profiles, unconstrained with respect to the calcite mass fraction in the deposition flux and in the surface solids, for stations that produced physically impossible fits (calcite mass fractions exceeding 100% either in the mass deposition flux or the solids concentration profile).Lines are as follows: long dashed -no calcite dissolution; full -calcite dissolution kinetics of order 4.5; dotted -order 2; short dashes -linear.

Figure 6 .
Figure 6.pH profiles for experiments with calcite mass fractions in the deposition flux limited to 90%.Lines are as follows: long dashed - Figure7.pH profiles with data-constrained surface sediment %calcite.Lines are as follows: long dashed -no calcite dissolution; full - between age tracers carried by different components) as the discrete representation of the underlying evolution equations of the age tracer and the age-carrying component have exactly the same structure and thus suffer from the same numerical dispersion.The original theory has been reformulated here in terms of production time instead of age.Production time is easier to handle than age in the evolution equations used in MEDUSA.That time remains constant in the absence of chemical reactions and mixing whereas age will continue to evolve with time and thus require a sustained virtual advection or reaction, even if the material carrying the age tracer gets transferred to the CORELAY realm (where mixing and reactions are ignored).The amended theory is detailed in the technical report Deleersnijder et al. (2001)cks, or completely in terms of the C:N:P:O:H composition.The solid's colour class can be used for immaterial(volume-less)properties of solids, such as classical colour tracers or isotopic properties.Each component of this class is linked to another solid from which it inherits physical and chemical properties.For components in the solid 's production time class, the code generator produces adequate equations for age (concentration) tracers.These equations follow the Constituent-oriented Age and Residence time Theory, CART, developed byDelhez et al. (1999)andDeleersnijder et al. (2001).The CART approach provides a means to avoid numerical differentiation in the calculated evolution of the age tracer and the age-carrying component (and also : "Early Diagenesis in Sediments.A one-dimensional model formulationwhich is included : " in the Supplement.

Table 5 .
Execution times of the COUPSIM best-fit experiment.One unit of CPU time corresponds here to 9:14.73 min.The executing platform had an Ubuntu 16.04.6LTS operating system (64-bit kernel 4.4.0-185-generic)running on a 1.90 GHz Intel® Core™ i7-8650U CPU; the codes were compiled with GFORTRAN 5.4.0, using optimisation level -O0.Experiments carried out with the MPI interfaces used OPENMPI v. 1.10.2 and were run on as many processes as indicated in brackets.