Microbial Alkalinity Production and Silicate Alteration in Methane Charged Marine Sediments: Implications for Porewater Chemistry and Diagenetic Carbonate Formation

A numerical reaction-transport model was developed to simulate the effects of microbial activity and mineral reactions on the composition of porewater in a 230-m-thick Pleistocene interval drilled in the Peru-Chile Trench (Ocean Drilling Program, Site 1230). This site has porewater profiles similar to those along many continental margins, where intense methanogenesis occurs and alkalinity surpasses 100 mmol/L. Simulations show that microbial sulphate reduction, anaerobic oxidation of methane, and ammonium release from organic matter degradation only account for parts of total alkalinity, and excess CO2 produced during methanogenesis leads to acidification of porewater. Additional alkalinity is produced by slow alteration of primary aluminosilicate minerals to kaolinite and SiO2. Overall, alkalinity production in the methanogenic zone is sufficient to prevent dissolution of carbonate minerals; indeed, it contributes to the formation of cemented carbonate layers at a supersaturation front near the sulphate-methane transition zone. Within the methanogenic zone, carbonate formation is largely inhibited by cation diffusion but occurs rapidly if cations are transported into the zone via fluid conduits, such as faults. The simulation presented here provides fundamental insight into the diagenetic effects of the deep biosphere and may also be applicable for the long-term prediction of the stability and safety of deep CO2 storage reservoirs.


INTRODUCTION
Microbial activity below the seafloor (in the deep biosphere) affects global cycling of carbon. In sediments on many continental margins, microbes convert buried organic matter through a series of reactions to methane (CH 4 ) and dissolved inorganic carbon (DIC), both of which can return back to the water column through advection or diffusion (e.g., Dickens, 2003;Krumins et al., 2013). In some areas, high rates of organic carbon decomposition lead to oversaturation of porewater with respect to different carbonate minerals (e.g., Baker and Burns, 1985;Moore et al., 2004;Meister et al., 2006;Meister et al., 2007;Meister, 2015;Wehrmann et al., 2016). Precipitation of diagenetic carbonates can even occur within carbonate-free ocean margin sediment sequences (e.g., Pisciotto and Mahoney, 1981;Kelts and McKenzie, 1984;Baker and Burns, 1985). Importantly, these carbonates add to the total worldwide carbon burial flux, and variations in their accumulation over time may have contributed to past changes in global carbon cycling, such as during times of widespread anoxia (Schrag et al., 2013;Sun and Turchyn, 2014). Nevertheless, gaps remain in the understanding of organic carbon decomposition and diagenetic carbonate formation, particularly in sediment sequences characterized by strong methanogenesis.
While production of DIC drives diagenetic carbonate precipitation, a simultaneous increase in the pH buffering capacity and total alkalinity must happen to increase the activity of dissolved CO 3 2− and, hence, the saturation state with respect to various carbonate phases. The total alkalinity (titration alkalinity) has been defined as "excess of proton acceptors over donors with respect to carbonic acid": TA [ (Dickson, 1981;Middelburg et al., 2020). Several microbially mediated reactions increase alkalinity, including iron-reduction (Raiswell and Fisher, 2000;Wehrmann et al., 2009), organoclastic sulphate reduction (OSR), and anaerobic oxidation of methane (AOM; Moore et al., 2004;Ussler and Paull, 2008;Meister, 2013). In general, these processes progress with depth below the seafloor, according to their redox potential and free energy yield (e.g., Froelich et al., 1979). They also affect alkalinity differently. Notably, OSR produces two moles of DIC and two moles of alkalinity per mole of sulphate consumed, with a 1:1 ratio of TA:DIC (Eq. 1): By contrast, AOM generates more alkalinity relative to DIC (TA: DIC 2:1; Eq. 2): In the absence of any terminal electron acceptor, subseafloor fermentation of organic matter can still occur, releasing CH 4 , organic acids, and CO 2 to porewater. A dominant overall process is "methanogenesis", which proceeds mainly through autotrophic reduction of CO 2 with H 2 as an electron donor (Conrad, 1999;Meister and Reyes, 2019): Crucially, the overall pathway results in excess CO 2 , but without producing further alkalinity. A current problem in our collective knowledge of the deep biosphere and the role of methanogenesis for subseafloor biogeochemical processes revolves around the fate of CO 2 and extreme porewater alkalinity. Production of CO 2 (Eq. 3) should lead to acidification and undersaturation of porewater with respect to carbonate minerals. This is because CO 2 dissolves in water and generates protons, but does not change alkalinity (Eq. 4): However, porewaters along continental margins and within zones of significant methanogenesis consistently exhibit very high alkalinity, often exceeding 50 mmol/L ( Figure 1) and sometimes surpassing 160 mmol/L (Wefer et al., 1998). Furthermore, the high alkalinity mostly reflects dissolved HCO 3 − concentrations (except near depths of HS − production), the pH usually exceeds 7, and early diagenetic carbonate (dolomite) often has positive δ 13 C values, indicative of precipitation under methanogenenic conditions (Claypool and Kaplan, 1974;Kelts and McKenzie, 1984). As side-stepped in various works (e.g., Chatterjee et al., 2011;Meister et al., 2011), there is a "lost proton problem". How can sufficient alkalinity, as HCO 3 − , occur in such sub-seafloor environments to the point of driving significant amounts of carbonate precipitation?
Several additional sources of alkalinity in marine sediments have been discussed in previous studies (e.g., Chatterjee et al., 2011;Meister et al., 2011). Given a nominal C/N ratio >7:1, organic matter degradation generally produces ammonia in most continental margin sediments (e.g., Arndt et al., 2013). Due to uptake of a proton, each mole of NH 4 + (from NH 3 ) also produces one mole of alkalinity: While NH 4 + accumulates in the porewater, it may adsorb to clay mineral surfaces, in exchange for Ca 2+ and Mg 2+ , which are readily released into the porewater Ockert et al., 2014;Mavromatis et al., 2015). Dissolved organic species, which almost certainly form in methanogenic systems, also need consideration (Smith, 2005).
Several authors (e.g., Aloisi et al., 2004;Wallmann et al., 2008;Meister et al., 2011;Archer et al., 2012;Solomon et al., 2014) have suggested that alteration of silicate minerals, particularly clay minerals, occurs within methanogenic zones and may produce sufficient alkalinity to drive diagenetic carbonate precipitation. At low pH and low temperatures, silicates may alter to secondary minerals by loss of oxides of alkali and earth alkali metals, which react to hydroxides upon protonation, thereby buffering acidification of dissolved CO 2 (Amiotte Suchet et al., 2003;Colbourn et al., 2015;Eq. 6): Often silicate alteration reactions also show a loss of SiO 2 , which upon water uptake reacts to H 4 SiO 4 and which, hence, does not contribute to total alkalinity production. As a consequence of Eq. 6, alkalinity also can be expressed by the balance of conservative ions. This has been shown by Wolf-Gladrow et al. (2007), who derived an explicit conservative expression of total alkalinity: TA ec [Na + ] + 2 [Mg 2+ ] + 2 [Ca 2+ ] TSO 4 ] (here only including the components relevant for the considered system). The advantage of this expression is that all parameters are conservative and can be directly measured. TA ec also provides a clearer understanding of how mineral reactions affect the porewater alkalinity.
Both ion exchange with NH 4 + and silicate mineral alteration may enrich porewaters in dissolved Mg 2+ , Na + , and K + Wallmann et al., 2008), while Ca 2+ typically becomes depleted due to carbonate precipitation. But which scenario is correct, and why and how extremely high alkalinity arises in zones of high microbial activity, remains unclear. Understanding the production of alkalinity and its effect on carbonate equilibrium is not only essential for understanding diagenetic carbonate formation and the natural carbon cycle, but also links to long-term storage of CO 2 in geological rock reservoirs (Bickle et al., 2007;Kasina et al., 2014).
In this study, we examine data from Ocean Drilling Program (ODP) Site 1230 on the Peru Margin (D'Hondt et al., 2003;Donohue et al., 2006). This site was drilled to examine microbiological and porewater expressions in an exemplary modern-day methanogenic system. Although Site 1230 has one of the most complete porewater datasets in ocean drilling history, the available information has yet to be placed into context. We developed a numerical reaction-transport model to simulate alkalinity production and its effect on the carbonate equilibrium as a result of microbial activity and mineral reactions in natural methanogenic zones over millions of years. Rates of microbial metabolic reactions (OSR and methanogenesis) are linked to an overall organic matter degradation rate. The rate of organic matter degradation (R TOC ; TOC total organic carbon) was determined using a presumed decay function and determining the decay parameters by fitting the simulated porewater profiles to the measured porewater data. Also rates of ion exchange on clay mineral surfaces and rates of dissolution/precipitation of the minerals detected by X-ray diffraction analysis were determined in this way. We then assessed how the different processes affect total alkalinity and carbonate saturation. This simulation aims at clarifying the longstanding enigma of how diagenetic carbonates form in deep methanogenic zones below the seafloor.

STUDY SITE AND MEASURED POREWATER PROFILES
Intense coastal upwelling and high primary productivity characterize surface waters along the Peruvian continental margin. In 2002, ODP Leg 201 drilled and cored Site 1230 on the lower slope of the Peru margin (9°6.7525′ S / 80°35.0100′ W; Supplementary Figure S1A) at 5,086 m water depth, and within 100 m from ODP Site 685. The site comprises five holes (A-E), from which cores collected sediment to 278 m below seafloor (mbsf). Except for a few short pressure cores to collect gas, sediment recovery over the upper 216 m was accomplished using an advanced hydraulic piston core (APC) tool while recovery below this depth was accomplished mostly using an extended core barrel (XCB) tool.
The overall sequence has two main lithological units, which match those recovered at Site 685 (Supplementary Figure S1B; Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 756591 D' Hondt et al., 2003). The upper 230 m (Unit I) consist of Pleistocene-Holocene diatom ooze with variable mud content, traces of biogenic carbonate (nannofossils and foraminifera) and commonly pyrite. Total organic carbon content ranges between 2 and 3.5 wt%, with a nearly constant C/N ratio of about 9 (Prokopenko et al., 2006). By contrast, total inorganic carbon (TIC) content of bulk sediment mostly lies below 0.2 wt%, but reaches up to 1 wt% in the middle part of the unit (Meister et al., 2005). A 10-m-scale cyclicity in chromaticity (D'Hondt et al., 2003) spans the upper 118 m and reflects alternating domains enriched in diatom-ooze or dark-green clay, presumably linked to glacial and interglacial intervals. A change in chromaticity coincides with an interval of lower diatom content between 118 and 148 mbsf. The clay-rich diatom ooze shows signs of increasing consolidation and fissility towards the base of Unit I (D'Hondt et al., 2003).
A sharp boundary at 216 mbsf delineates Unit I and Unit II, and marks a downward shift to strongly consolidated Miocene sediments. The boundary likely represents a décollement surface within the accretionary prism (D'Hondt et al., 2003;Matmon and Bekins, 2006), as it is characterized by a 20-cm-thick fault breccia cemented by syntectonic dolomite (Meister et al., 2011) and a ∼4 Ma gap in age (Suess 1988).
Porewater profiles show major changes in fluid composition with depth, including a steep and near linear negative sulphate gradient with a thin sulphate-methane transition zone (SMT) at approximately 9 mbsf. Sulphide reaches ca. 10 mmol/L at the SMT and decreases to zero at 20 mbsf. Below the SMT, a methanogenic zone prevails throughout the sequence. Extremely high concentrations of up to 300 mmol/L methane were reconstructed by Spivack et al. (2006), and gas hydrates were detected at 82 and 148 mbsf (D'Hondt et al., 2003). Measured DIC and alkalinity both reach concentrations of more than 150 mmol/L near 150 mbsf (D'Hondt et al., 2003). Below 150 mbsf, alkalinity, ammonia, and Mg 2+ concentrations decrease downwards towards the Unit I-II boundary (D'Hondt et al., 2003;Donohue et al., 2006). A sudden change in porewater chemistry occurs at this depth (D'Hondt et al., 2003), as indicated by more radiogenic Sr isotope values; likely this is due to fluid transported along the fault (Kastner et al., 1990;Meister et al., 2011).

Mineral Analysis
Sediment samples from ODP Site 1230 were selected from undisturbed core sections recovered from 6 to 254 mbsf. For bulk mineralogical analysis, powdered sediment samples were analysed by X-ray diffraction (XRD) using a Panalytical X'Pert PRO diffractometer (Cu-Kα radiation, 40 kV, 40 mA, step size 0.0167, 5 s per step).
For clay mineral separations, six sediment samples were slightly crushed by hand to pieces of approximately 3 mm. Organic matter was removed by oxidation with diluted H 2 O 2 , and the samples were further disaggregated with a 400 W ultrasonic probe for 3 min. The <2 μm grain-size fraction was separated from bulk sediment by sedimentation in an Atterberg-cylinder for 24 h and 33 min (Engelhardt et al., 1974;Gaucher et al., 2004) and dried at 50°C. The homogenized samples were then saturated with K + and Mg 2+ ions and oriented samples were prepared by dispersing 10 mg of clay in 1 ml of water, pipetting the suspensions onto glass slides, and drying at room temperature. The oriented, K-and Mg-saturated samples were analysed in air-dried state and after vapour solvation with either ethylene glycol (K-samples) or glycerol (Mgsamples) at 60°C for 24 h to identify expandable clay minerals like smectite and vermiculite. Additional K-saturated samples were heated to 550°C to destroy kaolinite and expandable clay minerals (Eslinger and Pevear, 1988;Jasmund and Lagaly, 1993;Moore and Reynolds, 1997). Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 756591

Reaction Transport Model
Concentrations of dissolved species were simulated with a one-dimensional reaction-transport model ( Figure 2) according to Fick's second law of diffusion in the form given by Boudreau (1997): where C is the concentration of a solute (mmol/L), z is the depth in metres below seafloor (mbsf), dt is the time interval (years), and ϕ is porosity. The first term on the righthand side is a downward advection term, simulating burial with the sedimentation rate ω, whereby a steady-state compaction (dϕ/ dt 0; ϕ·ω is constant) was assumed. The downward decreasing porosity was calculated as: where ϕ 0 and ϕ ∞ are constant porosities at sediment depth z 0 and ∞, respectively, and b is the porosity e-folding distance (the distance over which (ϕ 0 -ϕ ∞ ) decreases by a factor of e). The second term on the righthand side of Eq. 7 accounts for diffusion, where D is the diffusion constant (m 2 /s; see below). The tortuosity τ was calculated from porosity according to Boudreau (1997): The third term on the righthand side of Eq. 7 represents any source or sink (or combination thereof), where R is the rate of production or consumption of the solute (mmol L −1 a −1 ). For dissimilatory metabolic reactions (sulphate reduction and methanogenesis, according to Eqs 1, 3), the source and sink terms are stoichiometrically related to the degradation rate of TOC. The degradation of TOC was calculated using the reactive continuum (RC) model of Boudreau and Ruddick (1991): where TOC(t) is the TOC at sediment age t, TOC 0 is the initial TOC upon sedimentation, and a RC and ν are fitting parameters in the reactive continuum model. The source/sink s TOC was calculated from the derivative of organic matter decay after time: where ρ s is the density of the dry sediment and M C is the molecular weight of carbon. Also, TOC degradation is linked stoichiometrically to the release of ammonium via a prescribed C/N ratio. For sulphate reduction and AOM, a Monod kinetic term was applied (Treude et al., 2003;Arndt et al., 2006;Arndt et al., 2009;Knab et al., 2008): where K s and K AOM are the half saturation constants for sulphate reduction and AOM, respectively, and k AOM is the first-order rate constant for methane. Most sulphide produced by sulphate reduction is consumed by Fe-sulphide formation at the sulphide-iron front below the SMT. To calculate the sink of sulphide, FeS-precipitation was stoichiometrically coupled with reductive dissolution of Fe, with a Monod-term to limit reaction at very small sulphide concentrations. The rate constant was varied to fit the measured sulphide profile. Due to the high pressure and low temperature at >5,000 m water depth, free gas phase should not reside within the porewater over the upper few hundred metres below the seafloor. However, dissolved CH 4 concentrations can surpass those for gas hydrate saturation. The saturation conditions were calculated according to an approach given by Tishchenko et al. (2005), based on Duan et al. (1992). Gas hydrate formation and dissociation was included in the transport model as a source/sink-term, with a rate depending on the oversaturation of dissolved CH 4 . Burial of the solid phase (given in mmol equivalents per litre of porewater) was calculated using the sedimentation rate.

Cation Exchange on Mineral Surfaces
For the calculation of cation exchange on clay minerals, we assumed that adsorbed ions are in equilibrium with surrounding pore fluids. This is justified, because radiotracer experiments (Ockert et al., 2014) show that Ca 2+ -NH 4 + exchange occurs very rapidly, over a few hours. The conditional exchange coefficients (K Ex ) were found by solving a partition function with respect to ammonium for each cation (here exemplary for Ca 2+ ): and by assuming the sum of all mole fractions (X) equals to one (Boatman and Murray, 1982). The cation concentrations were fitted to the porewater concentrations by varying the exchange coefficients. Rates of adsorption and desorption were calculated from the depth gradients of each adsorbed ion times the burial velocity and were included as source/sink terms in the reactiontransport function.

Speciation and Mineral Reactions
Measured concentrations of dissolved ions in porewater of marine sediment generally do not account for speciation. For example, Mg concentrations determined via inductively coupled plasma atomic emission spectrometry, such as done on ODP Leg 201, inextricably include those from Mg 2+ , MgCl + , and MgSO 4 (0) (e.g., Katz and Ben-Yaakov, 1980). For the aqueous speciation and calculation of mineral saturation states the program Phreeqc (Version 3.5.0; Parkhurst and Appelo, 2013) was used. The reaction-transport model was linked to the PhreeqRM module (Parkhurst and Wissmeier, 2015). From the initial ion concentrations (here assumed as seawater concentrations), an initial (hypothetical) solution containing up to 90 species was calculated ( Figure 2). The speciation was then renewed at each time step. Charge balance during transport was maintained, using the correction factor for Coulombic effects (Boudreau et al., 2004) and leaving chloride concentration variable. Charge balance also Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 756591 was maintained for reactions. The program determines ion activity product (IAP) and the solubility product (K SP ). The saturation indices (SI log IAPlog K SP ) of calcite, dolomite, and all silicate phases detected by XRD were calculated for in-situ temperature and pressure conditions. Rates of mineral precipitation or dissolution were calculated as k·SI (cf. Lasaga, 1998;Morse et al., 2007), whereby the rate constants k were assumed to exponentially decrease with depth, and the decay parameters were found by fitting the porewater profiles to measured data (Supplementary Material). Sources and sinks of ions were calculated in stoichiometric proportion to the amounts of mineral dissolved or precipitated. Since Al in porewater is extremely depleted, it was not possible with the present model to simulate reaction rates with respect to Al species using a reasonable time step. For this reason, Al was kept constant at a meaningful concentration, i.e., at an intermediate concentration at which kaolinite remains always supersaturated while other silicate minerals remain undersaturated. Mineral dissolution and precipitation were then stoichiometrically coupled so that no Al was added to or removed from the solution.
In Eq. 14 the activity coefficient γ i is calculated for ion "i" as a function of ionic charge z, the ion radius r, the ionic strength I of the solution, and the sum of interactions with each ion k in the solution, with the interaction coefficient ϵ between ions "i" and "k" and the molar concentration c k .

Parameterization and Boundary Conditions
The reaction-transport equation was solved using the finite differences method with an explicit-implicit scheme. The nonlinear terms, including different source and sink terms, were solved explicitly. The simulations were run until a steady state was reached. All parameters, their symbols, values (or ranges), units, and references are listed in Table 1. The porewater concentration data used for fitting the model are mostly ODP Leg 201 shipboard data (D'Hondt et al., 2003; supplemented with shore-based measurements of Ca, K, Mg, and Na by Donohue et al. (2006). For the initial solution, major ion concentrations of seawater (Table 1) were used, and the calculated solution was also set at the seafloor as an upper domain boundary condition. Elevated DIC concentrations may occur in deep-sea bottom water due to accumulation of CO 2 from aerobic respiratory production (Zeebe, 2007), which accordingly affects the pH of the boundary condition. The lower domain boundary was defined at 230 mbsf, where the concentrations are assumed as constant. This stipulation is necessary because diverse species show a gradient across the lower boundary due to the influence of fluid advection along the fault at 230 mbsf.
A minimal long-term sedimentation rate of 0.1 m/ka was assumed based on the thickness of the Pleistocene-Holocene interval of 230 m according to the Leg 112 age model (Suess, 1988). The porosity function was fitted to the measured porosity data. Diffusion coefficients were calculated for P, T, and salinity at each depth, using the functions and constants given in Boudreau (1997). No temperature and pressure derivatives of the diffusion coefficients were considered, since these effects were determined to be negligible.
We used the rate constants of AOM given in Arndt et al. (2006) and references therein and the half-saturation constant for organoclastic sulphate reduction as reported by Tarpgaard et al. (2017). Initial TOC content and parameters a and ν were found by fitting the concentration profiles of metabolites, such as sulphate and ammonium to the measured data and by fitting the SMT to 9 mbsf (Supplementary Material). In particular, varying a and ν also affects the curvature of the sulphate profile . Production of NH 4 + was stoichiometrically linked to the rate of TOC-degradation, assuming a C/N ratio of 9 (Prokopenko et al., 2006; also consistent with ratios given in; Burdige and Komada, 2013).
Cation concentrations were fitted to the measured porewater data by varying the adsorption constants, using a CEC of 100 meq/100 g based on Ockert et al. (2014). The conditional  For speciation and calculation of saturation states we used the parameters in the database for the specific ion interaction theory, given in the Phreeqc package (sit.dat database). The minerals were selected from the database according to the mineral assemblage detected by XRD and according to their stoichiometric composition commonly observed in marine sediments (e.g., Marinoni et al., 2008;Park et al., 2019; Table 2). For thermodynamic calculations in-situ pressure and temperature at each depth were used, with a linear temperature gradient between the mudline-temperature of 1.7 and 11.2°C at 278 mbsf at Site 1230 (D'Hondt et al., 2003).
Looking at a depth plot of the Mg + glycerol saturated clay fractions ( Figure 3B), it is obvious that only the uppermost two samples contain large amounts of smectite (18 Å, 4.9°2θ). The deeper samples show only small smectite peaks, but additionally contain another expandable clay mineral, vermiculite (14 Å, 6.3°2 θ). Illite, chlorite, and kaolinite are present in samples of all depths.

Biogeochemical Reactions and Effects on Alkalinity
The basic parameters used for the geochemical simulation, temperature, pressure, sedimentation rate, porosity, TOC content, and C/N ratio were fitted as good as possible to the measured data ( Figures 4A-D; Supplementary Figure S3). Modelled porewater profiles, shown in Figure 5, reach a steady state after <1 Ma, even though 2.3 Ma are needed for burial of the gas hydrates below 230 m at a sedimentation rate of 0.1 m/ka. Therefore, the time needed to reach a steady state is well within the time frame of deposition of the Pleistocene interval, throughout which the sediment composition does not fundamentally change. The geochemistry within the underlying Miocene interval is decoupled by fluid flow along the décollement, essentially setting the boundary conditions for the Pleistocene evolution of the porewater profiles.
The measured sulphate profile is well reproduced by the model ( Figure 5A), with an SMT at 9 mbsf. However, the rather linear gradient of the measured sulphate profile is difficult to reconcile with sulphate reduction rates of more than 1,000 pmol/cm 3 d typically measured near the surface using radiotracer experiments (e.g., Parkes et al., 2005), which are about three orders of magnitude higher than the ones used in the present simulation (cf. also several magnitudes lower values modelled by Wang et al., 2008). Such high rates of organoclastic sulphate reduction should result in a more curved sulphate profile  and a lower contribution of methanogenic activity. Our simulation clearly shows that, assuming a C/N ratio of 9.85 and a TOC 0 of 3.5 wt% and the organic carbon degradation parameters a 25,000 a and ν 0.33, not only the solid phase parameters, but also the distribution of sulphate and methane could be well reproduced (Supplementary Figure S3). Methane falls within the range of concentrations determined with the argon method  down to a depth of 60 mbsf ( Figure 5A), where gas hydrate saturation is reached. Below this depth, a large scatter in the methane data probably reflects a heterogeneous distribution of gas hydrates, which contrasts with the simulation showing gas hydrate contents increasing continuously with burial and microbial production of methane. Besides sulphate and methane, ammonium concentration provides a reliable indicator for dissimilatory microbial activity (Heini et al., 2015), and the ammonium profile is well reproduced by the model in consistency with the measured C/N ratios ( Figure 5B). It is important to notice Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 756591 that the ammonium curve strongly depends on the parameters used for the organic carbon degradation function. Generally, a more reactive organic matter pool, with a smaller value for "a", leads to a steeper increase of ammonium in the upper part and a lower increase in the lower part of the profile. Hence, ammonium distribution can be used as a further constraint to fine-tune the organic matter degradation function. As a result of the constrained metabolic activity, bicarbonate shows a steep increase above the SMT, while dissolved CO 2 reaches its highest concentrations only below 100 mbsf ( Figure 5C). The high CO 2 concentrations relate to the pH drop to about 6 in the methanogenic zone (short-dashed line in Figure 5D). The simulated DIC increases with depth to more than 250 mmol/L around 150 mbsf ( Figure 5E), which is considerably larger than measured data, which reach only ca. 150 mmol/L. This is expected, since a significant portion of DIC was probably lost during core recovery. In-situ loss of CO 2 (e.g., Paull et al., 1996) due to partitioning of CO 2 into escaping methane gas bubbles can be excluded, as no gas phase is present at the high pressure at 5,000 m water depth. Sampling loss of CO 2 would also explain the offset of the simulated pH from the measured pH. Indeed, the measured pH could be largely restored by re-equilibrating the simulated porewater solutions with 2 L of headspace per litre of porewater at atmospheric pressure (dashed lines in Figures 5D,E). This means that CO 2 escape would have mainly occurred by equilibration with a closed headspace during core recovery or sample storage (in the core liner or in a sample container), whereas an equilibration with the open atmosphere did not occur.
The top 10 mbsf are characterized by a strong increase in porewater alkalinity, which is largely the result of sulphate reduction and AOM. However, below the SMT (∼9 mbsf) methanogenesis would only cause a small increase in alkalinity, due to the release of ammonium. Phosphate concentrations only reach up to 0.5 mmol/L (D'Hondt et al., 2003), so that the effect on alkalinity would be small. Assuming microbial activity alone, alkalinity only reaches ∼100 mmol/L between 100 and 150 mbsf ( Figure 5E, dotted line), which is lower than in the measured profile. Thus, the model scenario that includes microbial activity without ion exchange or mineral reactions reproduces the concentration profiles of the main metabolites relatively well, but it cannot explain the concentrations of some of the cations (mainly Mg 2+ and K + ) and the full extent of alkalinity in the porewater.

Mineral Reactions and Their Effect on Porewater Chemistry
In the microbial activity scenario discussed above several conservative ions that count towards total alkalinity are not reproduced correctly (Figures 5F,G; dotted lines). Most prominently, the increase in Mg 2+ significantly contributes to total alkalinity (TA ec ). The increase in Mg 2+ may be a result of release in exchange for NH 4 + . It is well known that especially smectites have a high ion exchange capacity due to their large surface area and also due to their ability to expand (e.g., Boatman and Murray, 1982). Therefore, we tested the effect of ion exchange, assuming a relatively high CEC of 100 meq/100 mg, which is typical for smectites (Ockert et al., 2014), and an extremely high content of exchanger of 50% in the sediment (approx. 10% smectite occurs in the uppermost two samples), however this did not show any significant effect on  Meister et al., 2005); (D) measured C/N ratios from Site 1230 (Prokopenko et al., 2006). The solid lines indicate the values assumed for the model.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 756591 major cations Ca 2+ , Mg 2+ , Na + , K + , and NH 4 + in the porewater. This may seem surprising, because such a high content of an exchanger provides an ion exchange capacity on the order of hundreds of meq per liter of porewater solution. However, the minimal effect can be explained by the low sedimentation rate of 0.1 m/ka resulting in a very slow exchange of ions. Therefore, another explanation must be found for the discrepancy in simulated and measured alkalinity and conservative ion concentrations.
The alternative effect that could modify the ion content of the porewater would be recrystallization (Wallmann et al., 2008, and further references above). Unlike during ion exchange, the minerals are structurally decomposed and/or modified during this process, as commonly observed during the degradation of sheet silicates by weathering. Under diffusion-limited conditions in sediments and with a high rock/water ratio, concentrations would increase upon dissolution of some minerals until equilibrium is reached. If several different mineral phases are present, as in the modelled case, some minerals dissolve, others precipitate, allowing for larger mass transfer to occur (e.g., Michalopoulos and Aller, 1995). Effectively, the water acts as a medium, and mineral reactions are driven by the thermodynamic stabilities of the minerals under prevailing P/T conditions. Analysed samples contain various amounts of vermiculite, chlorite, illite, smectite, and kaolinite, but there is no significant change in the peak pattern with depth indicative of ongoing phyllosilicate alteration. Also, the depletion of smectite below 70 mbsf must result from changing sedimentation, because smectite remains supersaturated and would rather form than dissolve throughout the profile. The coexistence of more phases than allowed by the Gibbs phase rule is indicated by the results from XRD ( Figures 3A,B), so that equilibrium of the solution with all mineral phases cannot be reached. Instead, they are probably limited by reaction kinetics.
The model results show that increasing acidification leads to a minor decrease in the saturation indices of the silicate phases k-feldspar, albite, smectite, and illite, (Figures 6A,B), while they remain supersaturated at the SMT. In contrast, vermiculite and chlorite are strongly undersaturated below the SMT ( Figure 6C). Porewater ( Figure 6D) is most supersaturated with respect to kaolinite throughout the section and only slightly undersaturated at the sediment-water interface. Even if the simulated saturation indices are somewhat arbitrary, because the Al concentrations in porewater are assumed, the response to pH and ionic compositions seems reasonable. In particular, the suggested dissolution and precipitation reactions based on the saturation states correspond to known mineral alteration reactions: vermiculite and chlorite commonly weather to illite, and illite weathers to smectite (Chamley, 1989). According to the classical Goldich series, silicate-rich minerals are more resistant to chemical weathering at low temperature (Goldich, 1938;Kowalewski and Rimstidt, 2003). Most insoluble are Al 2 O 3 octahedra, and accordingly the solubility decreases with increasing Al-content (Chamley, 1989). Saturation states in Figure 6 are consistent with this rule, showing Al-poor minerals chlorite and vermiculite to be more undersaturated and more susceptible to acidification of porewater. In contrast, Al-rich minerals, feldspar and kaolinite are supersaturated. Kaolinite precipitation is then limited by Al available from the dissolving clay minerals. As a result of alteration of Al-poor to Alrich minerals, alkali and Earth alkali metal ions are preferentially released to the solution, along with some silica. During dissolution and precipitation, charge balance is maintained as the cations dissolve from the silicates in oxide form and the oxide is readily protonated upon dissolution, yielding free ions and H 2 O (Eq. 6), hence, contributing to the excess of conservative cations and increasing TA ec .
By allowing the undersaturated silicates to dissolve and supersaturated silicates to precipitate, while adjusting the kinetic constants of different mineral phases, it was possible to improve the fit of porewater profiles to the measured data, in particular the increase of Mg 2+ and K + (Figures 5F,G; solid lines). It was found that the simulated porewater profiles best fit to the measured data if the reaction rate decreases exponentially with depth (Supplementary Figure S4). The dissolution of minerals other than chlorite and vermiculite had no significant effect on porewater chemistry, as their departure from equilibrium was minimal (Supplementary Figure S5). While the dissolution of chlorite only delivers Mg 2+ , K-vermiculite also provides sufficient K + to reproduce the K + profile. Also, a minor amount of Na + may originate from mineral reactions, although albite, the Na-endmember of plagioclase, is not undersaturated and is, thus, unlikely a significant source of Na. After inclusion of additional conservative cations released from silicate alteration, the simulated alkalinity almost entirely matches the measured alkalinity ( Figure 5E, solid line).
Alternatively, also volcanic glass may react with porewater and provide in particular K + and Na + . Volcanic glass in ash layers is fairly reactive and may undergo the generalized reaction (Scholz et al., 2013): (K, Na) 2 O; MgO; CaO; ·Al 2 O 3 ; ·SiO 2 + H 4 SiO 4 + 6CO 2 + 3H 2 O → Al 2 Si 2 O 5 (OH) 4 + K + + Na + + Mg 2+ + Ca 2+ + 6HCO 3 − Although this has been suggested as a source of alkalinity for diagenetic carbonate formation elsewhere (e.g., Wehrmann et al., 2016), no major ash layers were reported from the sedimentary record at Site 1230. Interlayers of expandable smectites could furnish a further source of Na + . For example, a strong increase in Na + with depth, which is decoupled from the chlorinity, occurs at Namibia margin Site 1084 and was interpreted as being derived from clay-rich sediment with up to 30% smectite in the clay fraction (Kastanja et al., 2006). However, as discussed above, ion exchange alone did not significantly affect the porewater chemistry.
Below 150 mbsf, full speciation using the measured concentrations and alkalinity would result in a charge imbalance. Since the model conserves charges, it cannot be forced to adopt this unbalance. Chloride was arbitrarily chosen to compensate the excess of negative charges, so that charge balance is maintained, which results in a decrease in Cl − concentration. Despite this uncertainty near the bottom of the domain, the high alkalinity between 50 and 200 mbsf is real and can only be reached by including silicate alteration.
In addition to these mineral reactions, opal needs discussion as the sediment contains abundant diatoms (D'Hondt et al., 2003), as also indicated by an elevated baseline in the XRD-patterns. The measured concentrations of dissolved silica are clearly below the saturation of opal-A, but the simulations follow the data fairly well with opal-C/T as an equilibrium phase ( Figures 5H, 6D). Indeed, it is surprising that the dissolution of opal-A, which is available in large amounts, does not dominate silica concentrations in the porewater. But it is known that opal-A dissolution can be inhibited, as diatoms are covered by organic matter (Van Cappellen, 1996). Opal-C/T as the thermodynamically more stable phase then should form, according to Ostwald's step rule (e.g., Meister et al., 2014). The curve in the silica concentration profile near the seafloor is due to the fact that silica is strongly undersaturated in seawater. The slight increase with depth is explained by the increasing solubility with temperature.
Overall, alteration of silicate minerals contributes to the very high total alkalinity (∼150 mmol/L) measured at ODP Site 1230. This additional alkalinity effect significantly buffers acidification of the porewater by the production of CO 2 during microbial methanogenesis. In contrast, the production of alkalinity by anaerobic metabolisms (organoclastic sulphate reduction and AOM) and the dissimilatory release of ammonium alone would not be sufficient to explain the measured alkalinity. High (>80 mmol/L) alkalinity concentrations characterize porewaters recovered from within the uppermost few hundred metres below the seafloor at many drill sites along continental margins of the world. Away from the Peru Margin, examples include DSDP Site 262 (Timor Sea, Indian Ocean;Cook, 1974), ODP Sites 994, 995, and 997 (Blake Ridge, Atlantic Ocean;Paull et al., 1996), Site 1019 (California Margin, Pacific; Lyle et al., 1997), Site 1082 (Namibian Margin, Atlantic; Murray and Wigley, 1998) and IODP Sites U1426 and U1427 (Sea of Japan, west Pacific;Tada et al., 2015; see also the examples in Figure 1). All these sites have several commonalities: high accumulation rates of clay and organic carbon, a shallow SMT, and within underlying porewaters, high methane, elevated NH 4 + and K + , and low Ca 2+ concentrations. However, unlike Site 1230, some of them have low Mg 2+ concentrations in porewater and lack gas hydrate. We suspect that microbial methane production and consequent mineral reactions are ubiquitous processes on continental margins around the world, but the rates and details vary depending on parameters embedded in our modelling. These would include pressure and temperature, but also the rate and composition of sedimentary inputs, from organic carbon to aluminosilicates. Not all minerals react equally. In particular, the Mg-rich clay minerals, vermiculite and chlorite dominating in areas of physical weathering under cold and temperate conditions (Chamley, 1989), and to some extent probably also volcanic ash, provide a strong capacity to buffer the pH by increasing alkalinity in the porewater. Future modelling efforts at other sites should lead to further testing and refinement of ideas presented here.

Factors Controlling Carbonate Precipitation
Having established the effects of biogeochemical activity and silicate alteration on the alkalinity and DIC content, we can now assess the factors controlling the saturation state of carbonates ( Figure 6E). Calcite and dolomite are only slightly supersaturated or undersaturated in the bottom water, which is partially due to high amounts of CO 2 produced by aerobic respiration in the water column, below the oxygen minimum zone at a water depth of 5,000 m (e.g., Zeebe, 2007). In the top few metres of the sediment, pH generally converges to values near 7 as a result of sulphate reduction (Soetaert et al., 2007;Meister, 2013;Meister, 2014), whereby additional alkalinity from AOM results in a sharp terrace at the SMT at 9 mbsf (solid line in Figure 5D), and saturation indices of calcite and dolomite reach maxima ( Figure 6E). Ca 2+ concentration decreases from the surface to the SMT and also Mg 2+ shows a kink at this depth, suggesting recent or ongoing precipitation of calcite and dolomite. The kinks can be largely reproduced by the model if dolomite precipitation is included. Indeed, a hard-lithified dolomite from 6.5 mbsf shows δ 13 C values more negative than −30‰, which is a clear indication of an AOM-induced dolomite (Meister et al., 2007;Meister et al., 2011;Meister et al., 2019a).
In the methanogenic zone below 9 mbsf, DIC increase to 250 mmol/L would cause a pH decrease to <6, but taking into account alkalinity production from mineral reactions more moderate pH values above 6 are reached ( Figure 5D). The saturation indices of carbonates remain near saturation, which is due to the buffering effect of the alkalinity reaching 150 mmol/L near 100 mbsf. Hence, alkalinity production due to NH 4 + release and silicate alteration largely prevents the dissolution of carbonates in the methanogenic zone, as it would be expected based on CO 2 production alone. Carbonates may even become slightly Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 756591 supersaturated due to this effect, but their precipitation would be hampered by the slow supply of Ca 2+ due to the long diffusion distances from the seafloor. Calcite was detected by XRD in the bulk sediments in the top 150 mbsf, with higher amounts between 50 and 150 mbsf ( Figure 3A). This carbonate is clearly part of the sediment, whereas an ex-situ precipitation due to degassing during core recovery can be excluded, as from porewater with 4 mmol/L Ca 2+ only ca. 1 g CaCO 3 per kg sediment could precipitate, which would be far below the detection limit of the XRD analysis. Besides low-Mg calcite also traces of Mg-calcite were detected. While the calcite may be diagenetic, precipitation of significant amounts of carbonate is probably not possible at present core depth, due to lack of Ca supply. Thus, the disseminated calcite detected between 50 and 150 mbsf formed most likely in the past, when the according sediment was located sufficiently near to the sediment surface to allow for a sufficient supply of Ca 2+ from seawater, but not as shallow as during the formation of the lithified layer of dolomite at 6.5 mbsf.
The formation of hard lithified diagenetic carbonates vs. fine disseminated carbonate is a longstanding problem (cf. Garrison et al., 1984). At Site 1230 no dolomite layers occur between 6.5 and 220 mbsf throughout the zone of abundant methane production. Having demonstrated that porewaters are not carbonate-undersaturated in the methanogenic zone, it cannot be argued that past dolomite layers have been dissolved in the methanogenic zone. Apparently, they have never formed, despite the intense methane production that should have maintained a pronounced SMT. However, a possible explanation for the lack of dolomite layers could be that initially methane was not available for intense AOM at the SMT (scenario in Figure 7A). Due to the great water depth, a gas hydrate zone established and expanded as the Pleistocene interval increased in thickness. At this stage, disseminated calcite or dolomite could have formed ( Figure 7B). Only when the thickness of the methanogenic zone became sufficiently large to induce an SMT at shallow depth, a dolomite layer could have formed at 6.5 mbsf ( Figure 7C). The SMT may have migrated a few metres upward and downward, possibly due to variation of bottom water temperature or variation of TOC and sedimentation rate, which are not resolved by the simulation.
While carbonate formation in relation with AOM has been discussed, it is still not clear how carbonates (in particular dolomite) with a strongly positive carbon isotope signature may form. A 20-cm-thick, hard-lithified dolomite breccia was recovered from 230 mbsf at Site 1230, exhibiting a δ 13 C of up to 15‰ (Meister et al., 2011). Inorganic carbon in methanogenic zones typically shows positive δ 13 C values, as DIC becomes enriched and CH 4 depleted in 13 C as a result of fractionation within the methanogenic pathways (see Meister and Reyes, 2019 for details). At Site 1230, Meister et al. (2011) showed the influence of fluid with more radiogenic 87 Sr/ 86 Sr ratios transported along an active fault zone at 230 mbsf, providing Ca 2+ and inducing precipitation of a dolomitic fault breccia. As Ca 2+ was delivered into the midst of a zone with high alkalinity, large amounts of dolomite could have precipitated ( Figure 7C). The formation of a dolomite breccia showing a syntectonic cement structure with compromise boundaries (Meister et al., 2011) also suggests a rapid precipitation. However, these conditions are specific to the tectonic situation in an accretionary prism, where fluid is driven upwards, but they still do not explain how hard lithified dolomite layers form in methanogenic zones elsewhere.
More extensive dolomite layers showing a methanogenic δ 13 C signature are found, for example, in the Miocene Monterey Fm. of FIGURE 7 | Conceptual model to explain the formation of diagenetic carbonates in relation to the evolution of the methanogenic zone: (A) During Pleistocene sedimentation over the Miocene, dissimilatory reactions lead to the onset of a methanogenic zone. (B) Within the expanding methanogenic zone alkalinity builds up, but insignificant carbonates (disseminated calcite or dolomite; blue shaded area) are precipitated as the diffusion distance for Ca 2+ and Mg 2+ are too long. (C) The methanogenic zone has now further expanded and a gas hydrate zone has established. The gradients near the surface are steep, with a shallow SMT driving carbonate precipitation (lithified dolomites; blue bars). Also at depth, fluid flow along the décollement provides Ca 2+ , inducing carbonate precipitation. The two types of dolomite show contrasting δ 13 C values, derived from dissimilatory sulphate reduction and methanogenesis, respectively.
California (e.g., Murata et al., 1969;Burns and Baker, 1987;Mozley and Burns, 1993;Loyd et al., 2012) and in many organic carbon rich marine sediments drilled by ODP (e.g., Kelts and McKenzie, 1984;Thornburg and Suess, 1992;Rodriguez et al., 2000). As these sites are partially located on the shelf, they may often be affected by CH 4 loss under dynamic conditions (Arning et al., 2012;Contreras et al., 2013), and gas escape must be taken into account, whereby partitioning of CO 2 in methane bubbles could lead to escape of CO 2 and thereby inducing carbonate formation via pH increase (Meister and Reyes, 2019;Meister et al., 2019b). However, none of these effects would be expected in the deep sea, where free gas is not present at low temperatures. In order to simulate carbonate diagenesis on the shelf, gas transport will have to be included in future modelling studies.

CONCLUSION
The Neogene sediment sequence at ODP Site 1230 in the Peru-Chile Trench contains porewaters with extreme alkalinity, reaching ∼150 mmol/L at ca. 100 mbsf. Full-speciation reaction-transport modelling shows that these concentrations result from several reactions. Much of the alkalinity is produced by organoclastic sulphate reduction and AOM in the upper part of the profile but also through release of ammonium from organic matter degradation throughout the sequence. The alteration of silicate minerals, in particular vermiculite and chlorite to kaolinite and opal-C/T, contributes substantial amounts of Mg 2+ and K + and further increases alkalinity in the methanogenic zone. In combination, microbial processes and clay mineral alteration produce sufficient alkalinity within the methanogenic zone to buffer acidification caused by increased DIC (up to 250 mmol/ L) from dissimilation of organic matter and, thus, to prevent undersaturation and dissolution of carbonates. Within the methanogenic zone, diagenetic carbonate formation is largely calcium-limited, due to the long diffusion distances, but dolomite beds readily form if Ca 2+ is supplied, such as at the SMT or along a fault zone. The SMT dolomites generally show a negative δ 13 C signature, whereas dolomites forming within the methanogenic zone show a positive δ 13 C signature, but in the absence of a gas phase, a shift of the dolomitization front due to CO 2 degassing does not occur within the gas hydrate stability zone.
While our simulation provides insight into carbonate diagenesis in deep methanogenic zones, it also would be applicable to human-made CO 2 storage reservoirs. Our study shows that reservoirs rich in specific clay minerals, vermiculite and chlorite, would have an extremely high capacity to trap CO 2 as bicarbonate, although the kinetics of these reactions are rather slow. In any case, having a quantitative model at hand that can realistically simulate carbonate diagenesis in marine sediments will be essential to understand the role of sub-surface fluidmicrobe-mineral interactions in the global carbon cycle.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
PM: Design of the study, supervision, developing the model, interpretion, writing the manuscript; GH: Modelling with Phreeqc, wrote early version; EP: Advising numerical modelling, provided comments to the manuscript; SG: X-ray diffraction, clay mineral analysis, comments to the manuscript; GD: Shipboard and shore-based porewater analysis, comments to the manuscript; CB: Code writing in C++; BL: Code writing, code testing, advice for mathematical solutions, comments to the manuscript.