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₂O₂, 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²⁺ 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 (Mg-samples) 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).

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): ∂ ϕ C ∂ t =   − ∂ ϕ ωC ∂ z + ∂ ∂ z ( ϕ D τ 2 ⋅ ∂ C ∂ z ) ± ∑ ​ ϕ R (7) 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: ϕ ( z ) = ϕ ∞ + ( ϕ 0 − ϕ ∞ ) e − z / b (8) where ϕ₀ and ϕ_∞ are constant porosities at sediment depth z = 0 and ∞, respectively, and b is the porosity e-folding distance (the distance over which (ϕ₀–ϕ_∞) 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²/s; see below). The tortuosity τ was calculated from porosity according to Boudreau (1997): τ 2 = 1 − ln ϕ 2 (9)

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⁻¹ a⁻¹). 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): TOC ( t ) = TOC 0 [ a RC ( a RC + t ) ] ν (10) where TOC(t) is the TOC at sediment age t, TOC₀ 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: s TOC = ∂ ( TOC ρ s   ( 1 − ϕ ) 100  M C   ϕ ) / ∂ t (11) 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): s AOM = k AOM [ CH 4 ] [ SO 4 2 − ] K S , AOM + [ SO 4 2 − ] (12) 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₄ 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₄. Burial of the solid phase (given in mmol equivalents per litre of porewater) was calculated using the sedimentation rate.

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²⁺-NH₄ ⁺ 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²⁺): K NH 4 − Ca = X NH 4 2 X Ca ⋅ [ Ca ] [ NH 4 ] 2 (13) 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 reaction-transport function.

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²⁺, MgCl⁺, and MgSO₄ ⁽⁰⁾ (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 was maintained for reactions. The program determines ion activity product (IAP) and the solubility product (K_SP).

The saturation indices (SI = log IAP – log 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.

Activity correction was calculated according to the Brønsted-Guggenheim-Satchard model (also known as specific ion interaction theory, SIT; Brønsted, 1921; Guggenheim and Turgeon, 1955), using the formalism presented by Grenthe and Plyasunov (1997; Eq. 14). log γ i = − 0.51 ⋅  z i 2   I 1 + 1.5  r i   I +   ∑ k ϵ i , k  c k     (14)

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.

The reaction-transport equation was solved using the finite differences method with an explicit-implicit scheme. The non-linear 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; 2004) 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₂ 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 (Meister et al., 2013). Production of NH₄ ⁺ 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 exchange coefficients only control the relative proportions in which the ions are adsorbed or desorbed, and they were adjusted to fit the data. Downward decreasing dissolution rate constants for individual clay minerals were found by fitting the data (Supplementary Material), while the precipitation of kaolinite was linked via the transfer of aluminium (see discussion for details).

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).

Bulk XRD-analyses (Figure 3A) show quartz (peak at 3.34 Å; 26.67° 2θ) as the most abundant phase. The second most abundant phase is feldspar, whereby albite/anorthite (3.19 Å; 27.92° 2θ) is more abundant than K-feldspar (3.24 Å, 27.53° 2θ) in all samples. Small amounts of mica or illite (peak-positions at 10 Å, 8.83° 2θ) are present in all samples. Calcite is detected in the samples from 6 mbsf and 57 to 144 mbsf. Two separate carbonate phases, a low-Mg calcite at 3.03 Å (29.4° 2θ) and a Mg calcite at 3.02 Å (29.5° 2θ), occur. A halite (NaCl) peak at 2.82 Å (31.7° 2θ) is usually present as the samples were not washed before analysis.

All analysed samples contain different amounts of smectite, vermiculite, chlorite, kaolinite, and illite (Figure 3B). Additionally, the clay fractions contain small amounts of quartz, feldspar, and calcite (only in one sample at 70.3 mbsf). Here, the identification of the clay minerals is explained for the sample from 21.8 mbsf (Supplementary Figure S2). Smectite was identified by a broad peak at 14.1 Å (6.3° 2θ) with Mg saturation which shifted to 11.8 Å (7.5° 2θ) with K saturation and collapsed to 9.9 Å (8.9° 2θ) after heating to 550°C. Saturation of the K-sample with ethylene glycol expanded smectite again to 16.8 Å (5.3° 2θ), the Mg-sample with glycerol to 18 Å (4.9° 2θ). Chlorite was identified by the peaks at 14.1, 7.07, 4.7 and 3.53 Å (6.3°, 12.5°, 18.8°, and 25.1° 2θ), which did not change position during treatments. Illite peaks at 9.9, 4.97 Å, and 3.33 Å (8.9°, 17.8°, and 26.7° 2θ) also kept their positions. Kaolinite peaks at 7.13 Å and 3.57 Å (12.4° and 24.9° 2θ) disappeared after heating to 550°C. Mg-saturated vermiculite was recognized based on a strong peak at 14 Å (6.3° 2θ), which shifted to 10 Å (8.83° 2θ) with K saturation, while saturation with Mg and glycerol (MgGly) did not change the 14 Å peak position (only present in traces in the sample from 21.8 mbsf, but more abundant in the other samples).

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.

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³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 (Meister et al., 2013) and a lower contribution of methanogenic activity. Our simulation clearly shows that, assuming a C/N ratio of 9.85 and a TOC₀ 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 (Spivack et al., 2006) 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 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₂ reaches its highest concentrations only below 100 mbsf (Figure 5C). The high CO₂ 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₂ (e.g., Paull et al., 1996) due to partitioning of CO₂ 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₂ 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₂ 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²⁺ and K⁺) and the full extent of alkalinity in the porewater.

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²⁺ significantly contributes to total alkalinity (TA_ec). The increase in Mg²⁺ may be a result of release in exchange for NH₄ ⁺ (Von Breymann et al., 1990). 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 major cations Ca²⁺, Mg²⁺, Na⁺, K⁺, and NH₄ ⁺ 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₂O₃ 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 Al-rich 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₂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²⁺ 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²⁺, 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 , N a ) 2 O ; M g O ; C a O ; ⋅ A l 2 O 3 ; ⋅ S i O 2 ] + H 4 S i O 4 + 6 C O 2 + 3 H 2 O → A l 2 S i 2 O 5 ( O H ) 4 + K + + N a + + M g 2 + + C a 2 + + 6 H C O 3 − (15)

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₂ 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₄ ⁺ and K⁺, and low Ca²⁺ concentrations. However, unlike Site 1230, some of them have low Mg²⁺ 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.

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₂ 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²⁺ concentration decreases from the surface to the SMT and also Mg²⁺ 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 δ¹³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₄ ⁺ release and silicate alteration largely prevents the dissolution of carbonates in the methanogenic zone, as it would be expected based on CO₂ production alone. Carbonates may even become slightly supersaturated due to this effect, but their precipitation would be hampered by the slow supply of Ca²⁺ 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²⁺ only ca. 1 g CaCO₃ 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²⁺ 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 δ¹³C of up to 15‰ (Meister et al., 2011). Inorganic carbon in methanogenic zones typically shows positive δ¹³C values, as DIC becomes enriched and CH₄ depleted in ¹³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 ⁸⁷Sr/⁸⁶Sr ratios transported along an active fault zone at 230 mbsf, providing Ca²⁺ and inducing precipitation of a dolomitic fault breccia. As Ca²⁺ 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 δ¹³C signature are found, for example, in the Miocene Monterey Fm. of 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₄ loss under dynamic conditions (Arning et al., 2012; Contreras et al., 2013), and gas escape must be taken into account, whereby partitioning of CO₂ in methane bubbles could lead to escape of CO₂ 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.