# Diagenetic Self-Organization and Stochastic Resonance in a Model of Limestone-Marl Sequences

Authors: Ivan L’Heureux
Journal: Geofluids
 (2018)

---

## Abstract

Banded patterns in limestone-marl sequences (“rhythmites”) form widespread sediments typical of shallow marine environments. They are characterized by alternations of limestone-rich layers and softer calcareous-clayey material (marl) extending over hundreds of meters with a thickness of a few tens of meters. The banded sequences are usually thought to result from systematic variations in the external environment, but the pattern may be distorted by diagenetic nonlinear processes. Here, we present a reactive-transport model for the formation of banded patterns in such a system. The model exhibits interesting features typical of nonlinear dynamical systems: (i) the existence of self-organized oscillating patterns between a calcite-rich mode (“limestone”) and a calcite-poor one (“marl”) for fixed environmental conditions and (ii) bistability between these two modes. We then illustrate the phenomena of stochastic resonance, whereby the multistable system is driven by a small external periodic signal (the 100,000 years’ Milankovitch cycle comes to mind) that is too weak to generate oscillations between the states on its own. In the presence of random fluctuations, however, the system generates transitions between the calcite-rich and calcite-poor states in statistical synchrony with the external forcing. The signal-to-noise ratio exhibits many maxima as the noise strength is varied. Hence, this amplification effect is maximized for specific values of the noise strength.


---

## Body

## 1. Introduction

Rocks and minerals often exhibit rhythmic spatial variations in their chemical and/or physical properties over length scales varying from the micrometer to the kilometer [1–3]. Most of these rhythmic patterns are simple responses to systematic changes in the environmental conditions when the rock was formed (so-called extrinsic mechanisms [4]). For instance, varves in lake sediments directly reflect annual climate variations [5].On the other hand, in a nonlinear nonequilibrium geosystem, the presence of positive feedback may result in a pattern formation that is self-organized even in an external environment that does not change (intrinsic mechanisms) [2, 4]. Examples of such geochemical self-organization are found in the oscillatory zoning of solid solutions series [4, 6–9], Liesegang band formation [10–17], orbicular granites [7, 18], and cyclic layering in igneous bodies [19, 20].However, there exists another rhythmic pattern formation mechanism: stochastic resonance [21–25]. This class of phenomenon manifests itself in nonlinear systems that are multistable, whereby different sets of initial conditions lead to distinct asymptotically large-time solutions (so-called “attractors”) for the fixed parameter values. Assume that such a system is driven by a periodic external signal that is too weak on its own to induce transitions between the various attractors. However, if a random component is also included in the external signal, then transitions between distinct attractors become possible for sufficiently strong noise strength. Moreover, there exists a noise strength for which the response of the system at the frequency of the external signal is maximized. In other words, the presence of a random component amplifies the response of the system at the frequency of the external signal for some values of the noise strength.In this paper, we study a particular example of a geosystem exhibiting rhythmic patterns: calcareous rhythmites such as limestone-marl or limestone-shale sequences (Figure1). Calcareous rhythmites are widespread sediments typical of shallow marine environments. Examples are found in North America (Ordovician), Northern Europe (Silurian), Montana (Mississippian), and Central and Western Europe (Jurassic and Cretaceous) [26–30]. They are characterized by strong banded patterns alternating between weathering-resistant limestone-rich layers and less durable calcareous-clayey material. They form laterally extensive beds on the scale of hundreds of meters or even kilometers [28] and their thickness is typically a few tens of meters. The interlayer distance is usually of the order of several cm to dm [27]. The banded sequence is usually thought to be a result of external climatic fluctuations, such as those resulting from the Milankovitch cycle (extrinsic mechanism), but the pattern may be significantly distorted by diagenetic self-organization [28–31] or could even be entirely self-organized (intrinsic mechanism). This being said, the detailed pattern formation mechanism is still poorly understood [28].Figure 1
Limestone-marl rhythmites from the Vauréal formation (Baie de la Tour, Anticosti Island, Québec), dating from the Katian Stage (Upper Ordovician). The two walkers give an idea of the scale of the formation. Picture provided by André Desrochers.One of the classic models is due to Ricken [32] in which pressure dissolution of compacting interlayer calcite in the deep-burial environment provides the necessary components for the eventual cementation of the limestone layers. However, many limestone beds are not strongly compacted and there are strong pieces of evidence that cementation occurs early, in a shallower deposit. There are also pieces of evidence for dissolution of aragonite components and residual aragonite preserved in the limestone beds. Rather, Munnecke and Samtleben [26] proposed an early diagenetic model (Figure 2) in which immature aragonite-bearing sediment is preferentially dissolved through bacterial activity as the sediment is buried through an aragonite dissolution zone (ADZ) located at a fixed distance (typically a few decimeters) from the sediment-water interface. Calcite and terrigenous materials are not affected. Carbonate and calcium ions are then transported upwards and downwards by diffusion and potentially upwards by compaction-induced advection. Under propitious saturation conditions, the ions precipitate as calcite, thus cementing a portion of the sediment atop the ADZ. As this calcite-enriched section progresses through the ADZ, the aragonite dissolution rate decreases, production of carbonate and calcium ions is cut, and cementation ceases. As a result of burial, a new portion of immature fresh sediment passes through the ADZ and the cycle is ready to start again. Thus, through this diagenetic self-organization mechanism, a sequence of cemented calcite-rich layers (limestone) and carbonate-depleted sediment (calcareous marl) layers can be generated, even for constant external conditions.Figure 2
Schematic representation of Munnecke et al. [27] and Böhm et al.’s [28] conceptual model for the generation of self-organized rhythmites. Left diagram: as fresh sediment is buried, aragonite dissolution occurs in a bacterially active ADZ. Carbonate and calcium ions are transported by diffusion and advection and the sediment immediately above the ADZ is cemented by calcite (dominantly) and aragonite precipitation. Center diagram: as the sediment is buried, the cemented zone passes the ADZ with much-reduced aragonite dissolution rates. Right diagram: as the cemented zone gets buried below the ADZ, fresh sediment is made available again to strong aragonite dissolution, and the cycle repeats itself again. The result is a sequence of cemented sediments (limestone-rich) alternating with carbonate-depleted sections (marl).In Böhm et al. [28], a mathematical model implementing these ideas is proposed. The model is based on a cellular automaton approach in which a two-dimensional space is discretized, forming an array of cells of arbitrary size. The state of the system is reduced to three possible values: immature sediment, cemented sediment, and aragonite-depleted sediment. The dynamics of the system are defined by a set of rules occurring at every discrete time step (of value unity). The model was later modified to take differential compaction into consideration [29], but the reactive-transport dynamics of the dissolved components (calcium and carbonate) were not explicitly taken into account. Periodic self-organized solutions are generated in their models but they found that an external trigger was necessary to laterally synchronize the pattern. They also found that, in the presence of both self-organized diagenesis and an external fluctuating driving signal, the resulting pattern does not necessarily reflect the characteristics of the external signal: new frequencies may appear; others may disappear, while some frequencies are amplified and others reduced.In this paper, we propose a simple one-dimensional reactive-transport model that implements Munneke’s conceptual model. Our set of partial differential equations is presented in Section2 and describes the reactive transport of calcium and carbonate ions, the precipitation and dissolution of aragonite and calcite, and their advective transport, as well as the effects of compaction and chemical dissolution/precipitation on the sediment porosity. In Section 3, we present examples of numerical solutions. We establish the existence of diagenetic self-organization in the model, as well as multistability, that is, long-time asymptotic solutions that are characterized by qualitatively different compositions at the bottom of the sediment layer. This invokes the possibility of noise-induced transitions between the various states. In Section 4, we briefly review the salient features of stochastic resonance using the motion of a particle in a standard bistable potential. In Section 5, the diagenetic model is driven by a weak periodical signal (mimicking part of the Milankovitch cycle) superposed to random porosity variations at the sediment-water interface. Examples of noise-induced transitions between calcite-rich states (“limestone”) and calcite-poor states (“marl”) in stochastic synchronicity with the driving signal are presented. Some concluding comments are given in Section 6. Finally, a “List of Mathematical Symbols” Section is found.


## 2. Deterministic Diagenetic Model

In this section, we establish the basic conservation equations for fluid momentum and components mass, both solid (calcite, aragonite, and terrigenous materials) and dissolved (calcium and carbonate ions).


### 2.1. Mass Conservation Equations

We use a reference frame with thex-spatial axis positively oriented towards the bottom with an origin fixed at the sea-sediment interface. The porous sediment is characterized by the volume fraction occupied by water (porosity ϕ) and a solid matrix made up of calcite (C), aragonite (A), and insoluble terrigenous materials (T). The advection velocities of the fluid phase and of the solid matrix are denoted as u and w, respectively. Within the aragonite microbially induced dissolution zone, A dissolves and releases calcium Ca+2 and carbonate ions CO3-2 which are transported by advection and diffusion before possibly precipitating in the sediment (and cementing it) as calcite or aragonite. The weight of sediment above a certain position contributes to its compaction. The volume fraction of water is also affected by the dissolution and precipitation reactions. These physicochemical factors affect the porosity which, in turn, modifies the hydraulic conductivity of the sediment and feedback on the advective transport of solutes.We define the dimensionless compositions of solid componentsCi = mi/ms as the mass mi of i={A,C,T} per total solid dry sediment mass ms. Then, by definition,(1)CT+CA+CC=1.The mass of component i per volume of dry sediment Vs is Ciρs, where ρs is the total dry sediment density:(2)ρs≡msVs=∑imi∑iVi=∑imi∑imi/ρi=∑iCi∑iCi/ρi=∑iCiρi-1.Here, Vi is the volume occupied by solid i and ρi is its density. The mass of component i per total sediment volume is Ciρs(1-ϕ). The diagenetic equation for the solid component i reads as follows:(3)∂∂tCiρs1-ϕ=-∂∂xCiρs1-ϕu+Ri,where u is the velocity of the solid sediment matrix (counted positive downwards) and Ri is the net reaction rate of component i reported per unit sediment volume.Dividing (3) by ρi, summing over all components i, and using (2) result in a continuity equation for the porosity:(4)∂∂t1-ϕ=-∂∂x1-ϕu+∑iRiρi.We now define the concentration of dissolved speciesk=Ca+2,CO3-2 as c^k, the number of moles of k per unit pore water volume (in M units). The corresponding diagenetic equation is(5)∂∂tϕc^k=-∂∂xϕc^kw+∂∂xϕDk∂c^k∂x+R^k,where w is the velocity of the pore water (counted positive downwards), Dk(ϕ) is the (tortuosity-corrected) diffusion coefficient of component k, and R^k is the net reaction rate (moles of component k reacted per unit time and sediment volume). We will use the following tortuosity factor, often used in diagenetic modeling [33]:(6)Dk=Dko1-log⁡ϕ2,where Dko is the diffusion of component k in seawater.Equation (5) can also be applied to water itself, with c^w≅ρw/μw, where ρw is the density of seawater (considered constant) and μw is its molar mass. We thus get another version of the continuity equation for porosity:(7)∂ϕ∂t=-∂∂xϕw+μwR^wρw.Adding (4) and (7) gives(8)-∂∂x1-ϕu-∂∂xϕw+μwR^wρw+∑iRiρi=0.The reactive term μwR^w/ρw describes the net change in water volume fraction as a result of precipitation/dissolution of the solid components and is nothing but -∑iRi/ρi. Thus,(9)∂∂x1-ϕu+∂∂xϕw=0and (4) may be rewritten as follows:(10)∂∂tϕ=-∂∂xϕw-∑iRiρi.It shows how porosity is modified by compaction and the reactive kinetics of the precipitation/dissolution of minerals. Equation (9) indicates that(11)1-ϕu+ϕw=At,where A(t) is a position-independent factor.Finally, using (10), the diagenetic equation (5) is simplified to(12)∂∂tc^k=-w∂∂xc^k+1ϕ∂∂xϕDk∂c^k∂x+1ϕR^k+c^k∑iRiρi.


### 2.2. Solving for the Velocities

The conservation of momentum takes the form of Darcy law. In our situation, this is(13)ϕw-u=Kϕ1-1gρw∂p∂x,where K is the hydraulic conductivity, p is pressure, and g is the acceleration of gravity. The left-hand side describes simply the flow rate of water with respect to the solid matrix. The conductivity is related to the porosity through empirical expressions, such as the Karman-Cozeny law, for instance [34, 35]:(14)K=d2gρw180ηϕ31-ϕ2,where d is the average grain diameter and η the viscosity of water. We will collect all the constant terms in a prefactor β and use instead(15)K=βϕ31-ϕ2Fϕ.Here, β is an empirical factor with units of velocity and F(ϕ) is a correction factor that becomes important for ϕ close to unity. After Hsu and Cheng [36], it is expected that(16)limϕ→1K=10βϕ21-ϕ.Thus, we choose [36](17)Fϕ=1-exp⁡-101-ϕϕ.In general [37], one can write the fluid pressure as the difference between the overburden stress (weight of sediment per unit area) and the effective stress σe in the sediment. Thus,(18)∂p∂x=gϕρw+1-ϕρs-∂σe∂x.In turn, the effective stress can be assumed [38] to be a function of the porosity. Its differential is related to the differential of the porosity through a linear response function H(ϕ) [39] (the inverse -1/H is called the pore space compressibility in Mello [40]):(19)dσe=Hϕdϕ.The Darcy flow (see (13)) then becomes with the help of (11)(20)ϕw-u=A-u=-Kϕ1-ϕρsρw-1+KϕHϕgρw∂ϕ∂x.The solid phase velocity at the water-sediment interface must match the sedimentation rate S. Denoting ϕ0 as the porosity at the boundary (water-sediment interface) and ρs0 as the dry sediment density at the boundary, we can eliminate A and evaluate the solid matrix advection field u as(21)u=S-Kϕ01-ϕ0ρs0ρw-1+Kϕ0Hϕ0gρw∂ϕ0∂x+Kϕ1-ϕρsρw-1-KϕHϕgρw∂ϕ∂x.Equation (11) then gives the fluid advection velocity w:(22)w=u+A-uϕ=S-Kϕ01-ϕ0ρs0ρw-1+Kϕ0Hϕ0gρw∂ϕ0∂x-Kϕ1-ϕ2ϕρsρw-1+Kϕ1-ϕϕHϕgρw∂ϕ∂x.It is not the object of this work to solve for the full elastic sediment problem. Instead, we will assume that, in the absence of reactions, the porosity field (denotedϕnr for “no reaction”) obeys an empirical relation [41–43] of the form(23)ϕnr=ϕ∞+ϕ0-ϕ∞exp⁡-bσe,where b and ϕ∞ are parameters. b is termed compaction coefficient [44] or sediment compressibility [45]. Then, the elastic response function is(24)H=dσedϕnr=-1bϕnr-ϕ∞.Substituting (22) in (10) shows that the pore compressibility generates a porosity diffusion-like term with a diffusion coefficient(25)Dϕ=-1-ϕKϕHgρw=βϕ31-ϕ1bgρwϕnr-ϕ∞F.As will be clear below, we will use the molecular diffusion coefficient of calcium to choose appropriate length and time scales. Estimates ofDϕ (with b varying from 1 to 100 (kPa)−1 [39]) show that the porosity diffusion coefficient is small compared with the molecular calcium diffusion coefficient. We will therefore neglect the contributions proportional to H in the advection terms (see (21) and (22)). On the other hand, we will keep the small diffusion term in the porosity transport equation (10) for the purpose of facilitating the stability of the numerical algorithm used to solve the system. It will be sufficient for our purpose to take this diffusion coefficient as a constant and estimate it for ϕ≅ϕnr≅ϕin, where ϕin is the (homogeneous) initial porosity.Summarizing, we have eight field variables: three compositionsci, two concentrations c^k, the porosity ϕ, and the velocities u and w. The equations describing their evolution are given by the five diagenetic equations (see (3) and (12)) and the continuity equation (10) with the advection velocities taken from (21) and (22). Instead of the diagenetic equation for the reactively inert component T, we will use definition (1) to obtain CT.


### 2.3. Reactions

In this simple model, we use four reactions.(i) Dissolution of aragonite (reaction DA), A → Ca+2 + CO3-2, enhanced by anaerobic oxidation of methane in the Munneke model [28]. The reaction rate is modeled as a first-order one with respect to aragonite:(26)RDA=1-ϕk2CAρsΩDA,where k2 is another rate coefficient (units of inverse time, 1/a). The term CAρs is the mass of aragonite per solid volume and multiplying by the factor 1-ϕ reports the rate per total sediment volume. The undersaturation factor ΩDA is(27)ΩDA=1-c^Cac^CO3KAmθx,where KA is the solubility of aragonite and m is a reaction order. It is understood that this reaction occurs only when the system is undersaturated, 1-c^Cac^CO3/KA>0. Finally, the factor θ(x) is a characteristic function that is zero when x is outside the aragonite dissolution zone (ADZ) and one otherwise. The ADZ is characterized by the position of its top edge (xd) and its thickness hd.(ii) Precipitation of aragonite (reaction PA): Ca+2 + CO3-2→ A. For simplicity, we neglect nucleation. The corresponding growth rate per unit sediment volume can be written as a first-order reaction:(28)RPA=1-ϕk1CAρsΩPA,where k1 is another rate coefficient (units of inverse time, 1/a). The oversaturation factor ΩPA is(29)ΩPA=c^Cac^CO3KA-1m′,where m′ is a reaction order (different from m in general). It is understood that this reaction is only effective when the system is oversaturated, c^Cac^CO3/KA-1>0.(iii) Dissolution of calcite (reaction DC), C → Ca+2 + CO3-2. The reaction rate is modeled as a first-order one with respect to calcite:(30)RDC=1-ϕk4CCρsΩDC,where k4 is a rate coefficient (units of inverse time, 1/a). The undersaturation factor ΩDC is(31)ΩDC=1-c^Cac^CO3KCn,where KC is the solubility of calcite and n is a reaction order. It is understood that this reaction takes place only when the system is undersaturated, 1-c^Cac^CO3/KC>0.(iv) Precipitation of calcite (reaction PC): Ca+2 + CO3-2→ C. Here,(32)RPC=1-ϕk3CCρsΩPC,where(33)ΩPC=c^Cac^CO3KC-1n′.k3 is a rate coefficient (in 1/a) and n′ is a reaction order (different from n in general). It is understood that the reaction occurs only when the argument of ΩPC is positive.The net reaction rates appearing in (3) and (5) are therefore(34)RC=RPC-RDC;RA=RPA-RDA;R^Ca=R^CO3=RDA+RDC-RPA-RPCμA.The molar mass of aragonite (or calcite) μA appears since, by definition, R^k corresponds to a change in a number of moles per unit sediment volume per unit time.


### 2.4. Boundary and Initial Conditions

The diagenetic equations for the solid components are first-order partial differential equations. We will fix the composition of the young sediment and the porosity at the water-sediment interface. We will also fix the concentrations of dissolved calcium and carbonate at that interface. On the other hand, one supplementary boundary condition is needed for the dissolved components and the porosity. At the bottom boundaryx=L, we will assume that the system has no diffusive flux. Mathematically,(35)Ci0,t=Ci0;c^k0,t=c^k0;ϕ0,t=ϕ0;∂c^kx,t∂xL=0;∂ϕx,t∂xL=0.The initial conditions are chosen as spatially homogeneous constants:(36)Cix,0=Ci,in;c^kx,t=c^k,in;ϕx,0=ϕin,x≠0.


### 2.5. Constant Density Approximation

The densities of the three solid components are fairly close to each other (2.950, 2.710, and about 2.8 g/cm3 for A, C, and T, resp. [27]). As a first approximation, we will neglect the variations in these densities. The pattern formation mechanism discussed here does not crucially depend on this approximation and the resulting simplification of the diagenetic equations is significant. We thus set ρs equal to a constant defined by the initial composition of the sediment and let ρA≅ρC≅ρT≅ρs. Combining (3) and (4) gives(37)∂∂tCi=-u∂∂xCi+1ρs1-ϕRi-Ci∑jRj.We can verify that summing these equations over i leads to an identity, as it should.Finally, (12) for the dissolved species becomes(38)∂∂tc^k=-w∂∂xc^k+1ϕ∂∂xϕDk∂c^k∂x+1ϕR^k+c^k∑iRiρs.


### 2.6. Scaling

It is convenient to express the dynamics of the system in reduced (dimensionless) form. We scale the dissolved species concentrations withKC and the advection velocities with a typical sedimentation rate S. We introduce time scale T∗ and length scale X∗ defined as follows:(39)T∗=DCaoS2;X∗=DCaoS.The evolution equations then take the following scaled forms in terms of the dimensionless variables x′=x/X∗, t′=t/T∗, c^k′=c^k/KC; U=u/S, and W=w/S (and dropping the ′):(40)∂∂tCA=-U∂∂xCA-Da1-CACAΩDA-ν1ΩPA+λCACCΩPC-ν2ΩDC;(41)∂∂tCC=-U∂∂xCC+Daλ1-CCCCΩPC-ν2ΩDC+CACCΩDA-ν1ΩPA;(42)∂∂tc^k=-W∂∂xc^k+1ϕ∂∂xϕdk∂c^k∂x+Da1-ϕϕδ-c^kCAΩDA-ν1ΩPA-λCCΩPC-ν2ΩDC;(43)∂∂tϕ=-∂∂xWϕ+dϕ∂2ϕ∂x2+Da1-ϕCAΩDA-ν1ΩPA-λCCΩPC-ν2ΩDC.Here, the dimensionless parameters are defined as follows:(44)Da=k2T∗=k2DCaoS2;λ=k3k2;ν1=k1k2;ν2=k4k3;dk=DkDCao;dϕ=DϕDCao;δ=ρsμAKC.Da can be interpreted as a Damköhler number. The under- and oversaturation factors become(45)ΩPA=c^Cac^CO3KCKA-1m;ΩDA=1-c^Cac^CO3KCKAm′θx;ΩPC=c^Cac^CO3-1n;ΩDC=1-c^Cac^CO3n′.Here, θ(x) is a characteristic function indicating that aragonite dissolution is activated in the ADZ only. It is defined as 1 when xd/X∗<x<(xd+hd)/X∗ and zero otherwise. Finally, taking into account the arguments of Section 2.2, we neglect terms of order dϕ in the advection terms (except for the small porosity diffusion term). The scaled velocities in (40)–(43) are(46)U=1-Kϕ0S1-ϕ0ρs0ρw-1+KϕS1-ϕρsρw-1;(47)W=1-Kϕ0S1-ϕ0ρs0ρw-1-KϕS1-ϕ2ϕρsρw-1.


### 2.7. Parameter Values

As far as possible, we will adopt conditions corresponding to warm shallow waters (temperature of 25°C, salinity of 35‰). Beside the parameters defining the boundary conditions, we need reasonable values for the material parameters appearing in the model. Some of them are well constrained (densities, molar masses). The diffusion coefficients Dko in water at various temperatures are available [33]. The solubilities KA and KC and the kinetic parameters m′ (aragonite dissolution) and n (calcite precipitation) are taken from Mackenzie and Andersson [46]. For simplicity, we have taken m=m′ and n′=n.The geometrical parameters of the system are its lengthL, the position xd below the water-sediment interface of the ADZ, and its thickness hd. The value of xd is compatible with the value mentioned in Böhm et al. [28]. The sedimentation rate S at the interface is of the order of 0.1–0.01 cm/a [47].The inverse of the rate coefficients1/k1≈1/k2 and 1/k3≈1/k4 is chosen to be smaller than the characteristic age of the system L/S, so that the reactions have many opportunities to occur during the advective transport of a sediment parcel across the system. Values of the rates between 0.1 a−1 and 5 a−1 are comparable to those used in a recent paleoclimate model [48].The parameterβ in the constitutive relation defining the hydraulic conductivity is not so well constrained. Values between 0.01 to 1.0 cm/a have been successfully used in modeling steady state compaction in sea sediments [39] and we adopt similar values here. For instance, a value of β = 0.1 cm/a corresponds to a permeability of about 10−13 cm2 for ϕ = 0.7. In Jourabchi et al. [39], values of the sediment compressibility b between about 1 and 100 (kPa)−1 have also been used. Its specific value is not crucial here (as long as the porosity diffusion coefficient is small compared to molecular diffusion). With the choiceb = 5 (kPa)−1, the porosity field is free from the artificial oscillations often seen in the numerical solution of convective systems in the presence of shock effects (due here to the discrete dissolution zone). For the same reason, the specific value of ϕ∞ is not crucial. It has been chosen to be a small quantity (0.01).In the simulations, the parametersk1, k3, β, S, xd, hd, and L were varied, as well as the water-sediment interface conditions (Ci0,c^k0,ϕ0) and the initial (uniform) sediment composition. We report here only results pertaining to two base (typical) scenarios with varying initial and boundary porosities. The parameters for the base scenarios are given in Table 1. To fix the idea, an initial sediment composition of 30% calcite, 60% aragonite, and 10% terrigenous materials was used, with initial calcium and carbonate concentrations both equal to 0.5 (in units KC).Table 1
Parameter values and initial conditions adopted in two base scenarios.
Parameters Scenario A Scenario B Unit μ A 100.09 100.09 g/mol μ W 18.45 18.45 g/mol ρ A 2.950 2.950 g/cm3 ρ C 2.710 2.710 g/cm3 ρ T 2.8 2.8 g/cm3 ρ W 1.023 1.023 g/cm3 D C a o 131.9 131.9 cm2/a D C O 3 o 272.6 272.6 cm2/a β 0.1 0.01 cm/a b 5 10 (kPa)−1 K A 10−6.19 10−6.19 M2 K c 10−6.37 10−6.37 M2 k 1 = k 2 1.0 0.01 a−1 k 3 = k 4 0.1 0.001 a−1 m = m ′ 2.48 2.48 - n = n ′ 2.80 2.80 - x d 50 50 cm h d 100 100 cm L 500 500 cm S 0.1 0.01 cm/a C A 0 = C A , i n 0.6 0.6 - C C 0 = C C , i n 0.3 0.3 - c ^ C a 0 = c ^ C a , i n 0.326 0.326 mM c ^ C O 3 0 = c ^ C O 3 , i n 0.326 0.326 mM


### 2.8. Numerical Approach

Starting with an initial uniform distribution for the five variablesCA, CC, c^Ca, c^CO3, and ϕ, system (40)–(43) is solved numerically on a regular grid of N=200 spatial steps Δx. The time is discretized in units of time steps Δt. The algorithm used to solve the advective-diffusion equations (see (42) and (43)) is a semi-implicit Fiadeiro-Veronis scheme [49]. This second-order scheme interpolates smoothly between the usual Crank-Nicholson scheme (used for pure diffusive transport) and the upstream scheme (used for pure advective problems) according to the value of the grid Péclet numbers WΔx/2dk or WΔx/2dϕ. The nonlinear terms in the transport coefficients and in the reactive terms are estimated by the advanced projection method [50], whereby these terms are evaluated at each intermediate time t+Δt/2 by Taylor expansion about t. The advective equations (39) and (40) are solved using an explicit upstream scheme. When Δt is small enough (5 × 10−6), the algorithm is stable, convergent, and relatively fast. With a 2.53 GHz Intel® Core™ i5 CPU, it takes about 50 seconds to generate a solution for a dimensionless time t=1.


## 2.1. Mass Conservation Equations

We use a reference frame with thex-spatial axis positively oriented towards the bottom with an origin fixed at the sea-sediment interface. The porous sediment is characterized by the volume fraction occupied by water (porosity ϕ) and a solid matrix made up of calcite (C), aragonite (A), and insoluble terrigenous materials (T). The advection velocities of the fluid phase and of the solid matrix are denoted as u and w, respectively. Within the aragonite microbially induced dissolution zone, A dissolves and releases calcium Ca+2 and carbonate ions CO3-2 which are transported by advection and diffusion before possibly precipitating in the sediment (and cementing it) as calcite or aragonite. The weight of sediment above a certain position contributes to its compaction. The volume fraction of water is also affected by the dissolution and precipitation reactions. These physicochemical factors affect the porosity which, in turn, modifies the hydraulic conductivity of the sediment and feedback on the advective transport of solutes.We define the dimensionless compositions of solid componentsCi = mi/ms as the mass mi of i={A,C,T} per total solid dry sediment mass ms. Then, by definition,(1)CT+CA+CC=1.The mass of component i per volume of dry sediment Vs is Ciρs, where ρs is the total dry sediment density:(2)ρs≡msVs=∑imi∑iVi=∑imi∑imi/ρi=∑iCi∑iCi/ρi=∑iCiρi-1.Here, Vi is the volume occupied by solid i and ρi is its density. The mass of component i per total sediment volume is Ciρs(1-ϕ). The diagenetic equation for the solid component i reads as follows:(3)∂∂tCiρs1-ϕ=-∂∂xCiρs1-ϕu+Ri,where u is the velocity of the solid sediment matrix (counted positive downwards) and Ri is the net reaction rate of component i reported per unit sediment volume.Dividing (3) by ρi, summing over all components i, and using (2) result in a continuity equation for the porosity:(4)∂∂t1-ϕ=-∂∂x1-ϕu+∑iRiρi.We now define the concentration of dissolved speciesk=Ca+2,CO3-2 as c^k, the number of moles of k per unit pore water volume (in M units). The corresponding diagenetic equation is(5)∂∂tϕc^k=-∂∂xϕc^kw+∂∂xϕDk∂c^k∂x+R^k,where w is the velocity of the pore water (counted positive downwards), Dk(ϕ) is the (tortuosity-corrected) diffusion coefficient of component k, and R^k is the net reaction rate (moles of component k reacted per unit time and sediment volume). We will use the following tortuosity factor, often used in diagenetic modeling [33]:(6)Dk=Dko1-log⁡ϕ2,where Dko is the diffusion of component k in seawater.Equation (5) can also be applied to water itself, with c^w≅ρw/μw, where ρw is the density of seawater (considered constant) and μw is its molar mass. We thus get another version of the continuity equation for porosity:(7)∂ϕ∂t=-∂∂xϕw+μwR^wρw.Adding (4) and (7) gives(8)-∂∂x1-ϕu-∂∂xϕw+μwR^wρw+∑iRiρi=0.The reactive term μwR^w/ρw describes the net change in water volume fraction as a result of precipitation/dissolution of the solid components and is nothing but -∑iRi/ρi. Thus,(9)∂∂x1-ϕu+∂∂xϕw=0and (4) may be rewritten as follows:(10)∂∂tϕ=-∂∂xϕw-∑iRiρi.It shows how porosity is modified by compaction and the reactive kinetics of the precipitation/dissolution of minerals. Equation (9) indicates that(11)1-ϕu+ϕw=At,where A(t) is a position-independent factor.Finally, using (10), the diagenetic equation (5) is simplified to(12)∂∂tc^k=-w∂∂xc^k+1ϕ∂∂xϕDk∂c^k∂x+1ϕR^k+c^k∑iRiρi.


## 2.2. Solving for the Velocities

The conservation of momentum takes the form of Darcy law. In our situation, this is(13)ϕw-u=Kϕ1-1gρw∂p∂x,where K is the hydraulic conductivity, p is pressure, and g is the acceleration of gravity. The left-hand side describes simply the flow rate of water with respect to the solid matrix. The conductivity is related to the porosity through empirical expressions, such as the Karman-Cozeny law, for instance [34, 35]:(14)K=d2gρw180ηϕ31-ϕ2,where d is the average grain diameter and η the viscosity of water. We will collect all the constant terms in a prefactor β and use instead(15)K=βϕ31-ϕ2Fϕ.Here, β is an empirical factor with units of velocity and F(ϕ) is a correction factor that becomes important for ϕ close to unity. After Hsu and Cheng [36], it is expected that(16)limϕ→1K=10βϕ21-ϕ.Thus, we choose [36](17)Fϕ=1-exp⁡-101-ϕϕ.In general [37], one can write the fluid pressure as the difference between the overburden stress (weight of sediment per unit area) and the effective stress σe in the sediment. Thus,(18)∂p∂x=gϕρw+1-ϕρs-∂σe∂x.In turn, the effective stress can be assumed [38] to be a function of the porosity. Its differential is related to the differential of the porosity through a linear response function H(ϕ) [39] (the inverse -1/H is called the pore space compressibility in Mello [40]):(19)dσe=Hϕdϕ.The Darcy flow (see (13)) then becomes with the help of (11)(20)ϕw-u=A-u=-Kϕ1-ϕρsρw-1+KϕHϕgρw∂ϕ∂x.The solid phase velocity at the water-sediment interface must match the sedimentation rate S. Denoting ϕ0 as the porosity at the boundary (water-sediment interface) and ρs0 as the dry sediment density at the boundary, we can eliminate A and evaluate the solid matrix advection field u as(21)u=S-Kϕ01-ϕ0ρs0ρw-1+Kϕ0Hϕ0gρw∂ϕ0∂x+Kϕ1-ϕρsρw-1-KϕHϕgρw∂ϕ∂x.Equation (11) then gives the fluid advection velocity w:(22)w=u+A-uϕ=S-Kϕ01-ϕ0ρs0ρw-1+Kϕ0Hϕ0gρw∂ϕ0∂x-Kϕ1-ϕ2ϕρsρw-1+Kϕ1-ϕϕHϕgρw∂ϕ∂x.It is not the object of this work to solve for the full elastic sediment problem. Instead, we will assume that, in the absence of reactions, the porosity field (denotedϕnr for “no reaction”) obeys an empirical relation [41–43] of the form(23)ϕnr=ϕ∞+ϕ0-ϕ∞exp⁡-bσe,where b and ϕ∞ are parameters. b is termed compaction coefficient [44] or sediment compressibility [45]. Then, the elastic response function is(24)H=dσedϕnr=-1bϕnr-ϕ∞.Substituting (22) in (10) shows that the pore compressibility generates a porosity diffusion-like term with a diffusion coefficient(25)Dϕ=-1-ϕKϕHgρw=βϕ31-ϕ1bgρwϕnr-ϕ∞F.As will be clear below, we will use the molecular diffusion coefficient of calcium to choose appropriate length and time scales. Estimates ofDϕ (with b varying from 1 to 100 (kPa)−1 [39]) show that the porosity diffusion coefficient is small compared with the molecular calcium diffusion coefficient. We will therefore neglect the contributions proportional to H in the advection terms (see (21) and (22)). On the other hand, we will keep the small diffusion term in the porosity transport equation (10) for the purpose of facilitating the stability of the numerical algorithm used to solve the system. It will be sufficient for our purpose to take this diffusion coefficient as a constant and estimate it for ϕ≅ϕnr≅ϕin, where ϕin is the (homogeneous) initial porosity.Summarizing, we have eight field variables: three compositionsci, two concentrations c^k, the porosity ϕ, and the velocities u and w. The equations describing their evolution are given by the five diagenetic equations (see (3) and (12)) and the continuity equation (10) with the advection velocities taken from (21) and (22). Instead of the diagenetic equation for the reactively inert component T, we will use definition (1) to obtain CT.


## 2.3. Reactions

In this simple model, we use four reactions.(i) Dissolution of aragonite (reaction DA), A → Ca+2 + CO3-2, enhanced by anaerobic oxidation of methane in the Munneke model [28]. The reaction rate is modeled as a first-order one with respect to aragonite:(26)RDA=1-ϕk2CAρsΩDA,where k2 is another rate coefficient (units of inverse time, 1/a). The term CAρs is the mass of aragonite per solid volume and multiplying by the factor 1-ϕ reports the rate per total sediment volume. The undersaturation factor ΩDA is(27)ΩDA=1-c^Cac^CO3KAmθx,where KA is the solubility of aragonite and m is a reaction order. It is understood that this reaction occurs only when the system is undersaturated, 1-c^Cac^CO3/KA>0. Finally, the factor θ(x) is a characteristic function that is zero when x is outside the aragonite dissolution zone (ADZ) and one otherwise. The ADZ is characterized by the position of its top edge (xd) and its thickness hd.(ii) Precipitation of aragonite (reaction PA): Ca+2 + CO3-2→ A. For simplicity, we neglect nucleation. The corresponding growth rate per unit sediment volume can be written as a first-order reaction:(28)RPA=1-ϕk1CAρsΩPA,where k1 is another rate coefficient (units of inverse time, 1/a). The oversaturation factor ΩPA is(29)ΩPA=c^Cac^CO3KA-1m′,where m′ is a reaction order (different from m in general). It is understood that this reaction is only effective when the system is oversaturated, c^Cac^CO3/KA-1>0.(iii) Dissolution of calcite (reaction DC), C → Ca+2 + CO3-2. The reaction rate is modeled as a first-order one with respect to calcite:(30)RDC=1-ϕk4CCρsΩDC,where k4 is a rate coefficient (units of inverse time, 1/a). The undersaturation factor ΩDC is(31)ΩDC=1-c^Cac^CO3KCn,where KC is the solubility of calcite and n is a reaction order. It is understood that this reaction takes place only when the system is undersaturated, 1-c^Cac^CO3/KC>0.(iv) Precipitation of calcite (reaction PC): Ca+2 + CO3-2→ C. Here,(32)RPC=1-ϕk3CCρsΩPC,where(33)ΩPC=c^Cac^CO3KC-1n′.k3 is a rate coefficient (in 1/a) and n′ is a reaction order (different from n in general). It is understood that the reaction occurs only when the argument of ΩPC is positive.The net reaction rates appearing in (3) and (5) are therefore(34)RC=RPC-RDC;RA=RPA-RDA;R^Ca=R^CO3=RDA+RDC-RPA-RPCμA.The molar mass of aragonite (or calcite) μA appears since, by definition, R^k corresponds to a change in a number of moles per unit sediment volume per unit time.


## 2.4. Boundary and Initial Conditions

The diagenetic equations for the solid components are first-order partial differential equations. We will fix the composition of the young sediment and the porosity at the water-sediment interface. We will also fix the concentrations of dissolved calcium and carbonate at that interface. On the other hand, one supplementary boundary condition is needed for the dissolved components and the porosity. At the bottom boundaryx=L, we will assume that the system has no diffusive flux. Mathematically,(35)Ci0,t=Ci0;c^k0,t=c^k0;ϕ0,t=ϕ0;∂c^kx,t∂xL=0;∂ϕx,t∂xL=0.The initial conditions are chosen as spatially homogeneous constants:(36)Cix,0=Ci,in;c^kx,t=c^k,in;ϕx,0=ϕin,x≠0.


## 2.5. Constant Density Approximation

The densities of the three solid components are fairly close to each other (2.950, 2.710, and about 2.8 g/cm3 for A, C, and T, resp. [27]). As a first approximation, we will neglect the variations in these densities. The pattern formation mechanism discussed here does not crucially depend on this approximation and the resulting simplification of the diagenetic equations is significant. We thus set ρs equal to a constant defined by the initial composition of the sediment and let ρA≅ρC≅ρT≅ρs. Combining (3) and (4) gives(37)∂∂tCi=-u∂∂xCi+1ρs1-ϕRi-Ci∑jRj.We can verify that summing these equations over i leads to an identity, as it should.Finally, (12) for the dissolved species becomes(38)∂∂tc^k=-w∂∂xc^k+1ϕ∂∂xϕDk∂c^k∂x+1ϕR^k+c^k∑iRiρs.


## 2.6. Scaling

It is convenient to express the dynamics of the system in reduced (dimensionless) form. We scale the dissolved species concentrations withKC and the advection velocities with a typical sedimentation rate S. We introduce time scale T∗ and length scale X∗ defined as follows:(39)T∗=DCaoS2;X∗=DCaoS.The evolution equations then take the following scaled forms in terms of the dimensionless variables x′=x/X∗, t′=t/T∗, c^k′=c^k/KC; U=u/S, and W=w/S (and dropping the ′):(40)∂∂tCA=-U∂∂xCA-Da1-CACAΩDA-ν1ΩPA+λCACCΩPC-ν2ΩDC;(41)∂∂tCC=-U∂∂xCC+Daλ1-CCCCΩPC-ν2ΩDC+CACCΩDA-ν1ΩPA;(42)∂∂tc^k=-W∂∂xc^k+1ϕ∂∂xϕdk∂c^k∂x+Da1-ϕϕδ-c^kCAΩDA-ν1ΩPA-λCCΩPC-ν2ΩDC;(43)∂∂tϕ=-∂∂xWϕ+dϕ∂2ϕ∂x2+Da1-ϕCAΩDA-ν1ΩPA-λCCΩPC-ν2ΩDC.Here, the dimensionless parameters are defined as follows:(44)Da=k2T∗=k2DCaoS2;λ=k3k2;ν1=k1k2;ν2=k4k3;dk=DkDCao;dϕ=DϕDCao;δ=ρsμAKC.Da can be interpreted as a Damköhler number. The under- and oversaturation factors become(45)ΩPA=c^Cac^CO3KCKA-1m;ΩDA=1-c^Cac^CO3KCKAm′θx;ΩPC=c^Cac^CO3-1n;ΩDC=1-c^Cac^CO3n′.Here, θ(x) is a characteristic function indicating that aragonite dissolution is activated in the ADZ only. It is defined as 1 when xd/X∗<x<(xd+hd)/X∗ and zero otherwise. Finally, taking into account the arguments of Section 2.2, we neglect terms of order dϕ in the advection terms (except for the small porosity diffusion term). The scaled velocities in (40)–(43) are(46)U=1-Kϕ0S1-ϕ0ρs0ρw-1+KϕS1-ϕρsρw-1;(47)W=1-Kϕ0S1-ϕ0ρs0ρw-1-KϕS1-ϕ2ϕρsρw-1.


## 2.7. Parameter Values

As far as possible, we will adopt conditions corresponding to warm shallow waters (temperature of 25°C, salinity of 35‰). Beside the parameters defining the boundary conditions, we need reasonable values for the material parameters appearing in the model. Some of them are well constrained (densities, molar masses). The diffusion coefficients Dko in water at various temperatures are available [33]. The solubilities KA and KC and the kinetic parameters m′ (aragonite dissolution) and n (calcite precipitation) are taken from Mackenzie and Andersson [46]. For simplicity, we have taken m=m′ and n′=n.The geometrical parameters of the system are its lengthL, the position xd below the water-sediment interface of the ADZ, and its thickness hd. The value of xd is compatible with the value mentioned in Böhm et al. [28]. The sedimentation rate S at the interface is of the order of 0.1–0.01 cm/a [47].The inverse of the rate coefficients1/k1≈1/k2 and 1/k3≈1/k4 is chosen to be smaller than the characteristic age of the system L/S, so that the reactions have many opportunities to occur during the advective transport of a sediment parcel across the system. Values of the rates between 0.1 a−1 and 5 a−1 are comparable to those used in a recent paleoclimate model [48].The parameterβ in the constitutive relation defining the hydraulic conductivity is not so well constrained. Values between 0.01 to 1.0 cm/a have been successfully used in modeling steady state compaction in sea sediments [39] and we adopt similar values here. For instance, a value of β = 0.1 cm/a corresponds to a permeability of about 10−13 cm2 for ϕ = 0.7. In Jourabchi et al. [39], values of the sediment compressibility b between about 1 and 100 (kPa)−1 have also been used. Its specific value is not crucial here (as long as the porosity diffusion coefficient is small compared to molecular diffusion). With the choiceb = 5 (kPa)−1, the porosity field is free from the artificial oscillations often seen in the numerical solution of convective systems in the presence of shock effects (due here to the discrete dissolution zone). For the same reason, the specific value of ϕ∞ is not crucial. It has been chosen to be a small quantity (0.01).In the simulations, the parametersk1, k3, β, S, xd, hd, and L were varied, as well as the water-sediment interface conditions (Ci0,c^k0,ϕ0) and the initial (uniform) sediment composition. We report here only results pertaining to two base (typical) scenarios with varying initial and boundary porosities. The parameters for the base scenarios are given in Table 1. To fix the idea, an initial sediment composition of 30% calcite, 60% aragonite, and 10% terrigenous materials was used, with initial calcium and carbonate concentrations both equal to 0.5 (in units KC).Table 1
Parameter values and initial conditions adopted in two base scenarios.
Parameters Scenario A Scenario B Unit μ A 100.09 100.09 g/mol μ W 18.45 18.45 g/mol ρ A 2.950 2.950 g/cm3 ρ C 2.710 2.710 g/cm3 ρ T 2.8 2.8 g/cm3 ρ W 1.023 1.023 g/cm3 D C a o 131.9 131.9 cm2/a D C O 3 o 272.6 272.6 cm2/a β 0.1 0.01 cm/a b 5 10 (kPa)−1 K A 10−6.19 10−6.19 M2 K c 10−6.37 10−6.37 M2 k 1 = k 2 1.0 0.01 a−1 k 3 = k 4 0.1 0.001 a−1 m = m ′ 2.48 2.48 - n = n ′ 2.80 2.80 - x d 50 50 cm h d 100 100 cm L 500 500 cm S 0.1 0.01 cm/a C A 0 = C A , i n 0.6 0.6 - C C 0 = C C , i n 0.3 0.3 - c ^ C a 0 = c ^ C a , i n 0.326 0.326 mM c ^ C O 3 0 = c ^ C O 3 , i n 0.326 0.326 mM


## 2.8. Numerical Approach

Starting with an initial uniform distribution for the five variablesCA, CC, c^Ca, c^CO3, and ϕ, system (40)–(43) is solved numerically on a regular grid of N=200 spatial steps Δx. The time is discretized in units of time steps Δt. The algorithm used to solve the advective-diffusion equations (see (42) and (43)) is a semi-implicit Fiadeiro-Veronis scheme [49]. This second-order scheme interpolates smoothly between the usual Crank-Nicholson scheme (used for pure diffusive transport) and the upstream scheme (used for pure advective problems) according to the value of the grid Péclet numbers WΔx/2dk or WΔx/2dϕ. The nonlinear terms in the transport coefficients and in the reactive terms are estimated by the advanced projection method [50], whereby these terms are evaluated at each intermediate time t+Δt/2 by Taylor expansion about t. The advective equations (39) and (40) are solved using an explicit upstream scheme. When Δt is small enough (5 × 10−6), the algorithm is stable, convergent, and relatively fast. With a 2.53 GHz Intel® Core™ i5 CPU, it takes about 50 seconds to generate a solution for a dimensionless time t=1.


## 3. Deterministic Solutions

Figure3 illustrates the numerical solution obtained for scenario A in the case where both the initial porosity and the porosity at the sea-sediment boundary are equal to 0.5 and 0.6, respectively. Figures 3(a)–3(d) show the aragonite and calcite composition, the dissolved species concentrations, and the porosity fields at four different times, whereas Figure 3(e) shows the resulting final profiles. Notwithstanding the presence of temporary propagating waves in the mineral composition profiles, the final profile is time-independent and exhibits only a single calcite peak without banding patterns. The mineral composition at the bottom boundary in the steady state has 76.4% calcite and 0% aragonite (and 23.4% terrigenous materials). We will characterize this steady state by a qualitative assessment of the sediment mineral composition at its bottom boundary. Thus, the case illustrated here is denoted by 0-C (for no aragonite and nonzero calcite).Figure 3
Numerical solution for scenario A withϕin = 0.5 and ϕ0 = 0.6. Aragonite and calcite composition, dissolved species concentrations and porosity fields at four different times (a–d), and final steady state profiles (e).
(a)
t = 1685 a (b)
t = 3370 a (c)
t = 5055 a (d)
t = 6749 a (e)
Steady stateOne can understand the general features of the profiles as follows. The dissolved species concentrations reach a steady state profile relatively early in the simulation. Porosity decreases slightly with position as a result of compaction. Recall that, whereas calcite dissolution can occur wherever the system is locally undersaturated, dissolution of aragonite occurs only in the ADZ,xd<x<xd+hd. The calcium and carbonate ions released induce a local oversaturation with respect to calcite and a large calcite peak appears. The solid matrix advection velocity U remains positive throughout the simulation. Through downward advection, the bottom edge of the aragonite depletion zone is shifted out of the system, so that, eventually, no aragonite is left in the bottom part of the system. Meanwhile, the calcite peak splits into two smaller peaks. The topmost peak reaches a steady state, whereby the reaction kinetics balances the advective transport of the components, whereas the bottom peak is advected out of the system. Nevertheless, some residual calcite remains at the bottom of the system, hence a 0-C steady state.For scenario A, there exists a combination of initial and boundary porosities for which the large-time solution is oscillatory. Thus, it appears that this simple reactive-transport implementation of Munneke’s conceptual model is propitious to diagenetic self-organized periodic solutions. Figure4 illustrates two such situations. Here, the aragonite depletion zone is advected downwards and the bottom of the system is free from aragonite. On the other hand, the top edge of the aragonite dissolution zone releases carbonate and calcium ions, generating a broad precipitating calcite peak and a local porosity decrease as cementation processes (the term proportional to ΩPC in (43)). As a result of the porosity decrease, the pore water velocity decreases in absolute value (see (15) and (47)), thus decreasing the upstream advective transport of dissolved carbonate and calcium. Eventually, the calcite peak becomes undersaturated and starts to dissolve, increasing the porosity and increasing the advective transport of dissolved species until the cycle starts over again. More complex self-organized oscillatory solutions (with three calcite peaks per cycle) are also possible (Figure 4(b)).Figure 4
Numerical solutions for scenario A illustrating self-organized oscillations. The calcite composition at the bottom of the system (x=L) is plotted as a function of time. (a) ϕin = ϕ0 = 0.8. Oscillation period = 17,290 a. (b) Complex oscillation obtained with ϕin = 0.8 and ϕ0 = 0.74. Average time between troughs = 26,288 a.
(a) (b)Moreover, the model exhibits another interesting nontrivial signature of nonlinear behaviour: multistability. This occurs when distinct steady state solutions are obtained for different initial values of the variables, leaving all parameters unchanged otherwise. To illustrate this feature for scenario A, Figure5 shows the evolution of the sediment composition at the bottom of the system using two different initial porosities ϕin and leading to two different steady states (or attractors) as time is large: a 0-0 state (Figure 5(a)), a state without aragonite or calcite, and a state A-C with nonzero aragonite and calcite (Figure 5(b)). The set of initial values evolving towards an attractor of type i defines the basin of attraction of attractor i. For fixed ϕin, as the porosity at the top boundary ϕ0 is varied, the boundaries of the basins of attraction shift. Consequently, different steady states are also obtained when the parameter ϕ0 changes. Figure 6(a) illustrates the domain of existence of each attractor in the (ϕin, ϕ0) space for scenario A.Figure 5
Numerical solutions for scenario A illustrating bistability. All variables are taken at the bottom of the system (x=L) and are plotted as a function of time. (a) ϕin = 0.8 and ϕ0 = 0.72. A final steady state 0-0 is obtained. (b) ϕin = 0.55 and ϕ0 = 0.72. A final steady state A-C is obtained.
(a)
ϕ i n = 0.8; ϕ0=0.72 (b)
ϕ i n = 0.55; ϕ0 = 0.72Figure 6
Phase diagram in (ϕin, ϕ0) parameter space for scenario A. The corresponding attractors are denoted by their steady state characteristic (0-0, A-C, 0-C) or their oscillatory behaviour (Osc).One can understand the nature of the steady states as follows. A state 0-C (as illustrated in detail in Figure3) corresponds to a situation where the advection velocity U remains positive and the dissolution of aragonite is complete. As mentioned above, the bottommost edge of the dissolution-induced aragonite depletion zone is then advected outside the system so that no aragonite is left at the bottom of the sediment, at x=L. On the other hand, although the sediment at x=L remains weakly undersaturated with respect to calcite, its dissolution is compensated by its advective transport. The adjacent 0-0 state differs from the 0-C state in that the undersaturation with respect to calcite is larger, so that all the calcite at x=L dissolves. The existence of aragonite in the A-C state is due to the following mechanism: if the porosity at the boundary is sufficiently larger than the porosity at some position, compaction will result in an upwards drift of the solid matrix (in (46), U becomes negative). If this occurs in the neighborhood of the aragonite dissolution zone, the depletion layer will shift towards the sediment surface, thus leaving a residual aragonite concentration at the bottom of the system close to its initial value. Finally, the oscillatory regime occurs for large initial and boundary porosities. Indeed, as explained above, the pore water advection term W plays an important role in generating pulsating solutions. As is seen from (47), the absolute value of W is higher for larger values of the boundary porosity and of the overall porosity, in agreement with Figure 6(a).Figure7 shows a phase diagram in (ϕin, ϕ0) space for scenario B (Table 1). Scenario B is characterized by overall slower kinetics: the rate coefficients of dissolution/precipitation and the sedimentation rate are smaller, and the sediment has a finer texture (smaller β). In this case, no self-oscillatory solutions are obtained but the system exhibits multistability between many steady states. This scenario was chosen to illustrate the presence of stochastic resonance in the limestone-marl system. But first, a review of the main features of stochastic resonance is presented for the sake of completeness.Figure 7
Phase diagram in (ϕin, ϕ0) parameter space for scenario B. The corresponding attractors are denoted by their steady state characteristic (0-0, A-C, 0-C or A-0).


## 4. Stochastic Resonance: A Primer

In the presence of an external stochastic signal, nonlinear dynamical systems may exhibit interesting nontrivial phenomena, such as noise-induced transitions [25, 51, 52], stochastic coherence [25], or stochastic resonance [24, 25]. It is the latter class of dynamical behaviour that we explore in our model. For the sake of completeness, only the salient features of the phenomenon of stochastic resonance are described. More details are found in the classic reviews of Jung [23] and Gammaitoni et al. [24].Three conditions are necessary to obtain stochastic resonance: (i) the deterministic (noiseless) dynamical system must exhibit multistability, whereby more than one attractor coexists, (ii) the system is driven by a weak deterministic periodic signal, and (iii) the external signal also includes a stochastic contribution. In the absence of an external signal, one can think of the dynamical system as moving in a multiwell potential. It evolves to one attractor or another, depending on its initial condition. In the presence of the weak external periodic signal, the perturbations to the system are too small to cause the system to transit from one potential well to another. The signal just oscillates slightly around one attractor as the potential barriers separating the basins of attraction fluctuate slightly. However, in the presence of noise, the fluctuations may induce the system to cross the potential barriers. If the noise parameters are chosen appropriately, the barrier-crossing dynamics can occur in stochastic synchrony with the weak periodic signal: these noise-induced transitions occur more easily as the potential barrier is smaller, and this happens once every cycle of the periodic signal. The result consists in a series of large transitions from one attractor to another, with a higher probability of this occurring once per cycle, even though the external periodic signal alone (without noise) is too weak to induce these transitions.To illustrate this qualitative description, one can use the standard case of a single-variable dynamical systemy(t) moving in a bistable symmetric potential Φ(y):(48)dydt=-dΦdy+αsin⁡ωt+ξt;Φy=-y22+y44.Here, α is the weak amplitude of the external periodic signal, ω is its frequency, and ξ(t) is an external stochastic driving term. In the simple standard case, this noise term is chosen to be white [51], with zero average and a time correlation characterized by noise intensity σ:(49)ξtξt′=2σδt-t′,where the symbol  refers to averaging with respect to the noise realizations and δ(t) is the Dirac-delta function. In the undriven deterministic case (α=σ=0), it is seen that the system is indeed bistable: two stable steady states coexist at y=±1 and an unstable fixed point is present at y=0. The potential barrier height separating the two basins of attraction y>0 and y<0 is ΔΦ=Φ(0)-Φ(±1)=1/4. In the presence of the periodic forcing signal, the potential barrier height fluctuates periodically but, in order to preserve the double well structure of the potential, α must be chosen smaller than 2/33.Figure8 illustrates the numerical solution of (48) and (49) using a stochastic first-order Euler algorithm [53]. We used α = 0.15 and ω = 2π/200, an initial condition y(0)=1, and we varied the noise intensity σ. It is seen that for small noise intensity (Figure 8(a)), transitions between y~1 and y~-1 occur rarely: the noise is not sufficiently strong to induce common transition between the steady states. On the other hand, when the noise is too strong (Figure 8(c)), transitions occur often but no coherent pattern is obtained in the transition frequency. For intermediate values of the noise intensity (Figure 8(b)), transitions occur in stochastic synchrony with the external periodic signal. An example of power spectrum is illustrated in Figure 8(d) for the case corresponding to Figure 8(a) (σ = 0.05). The signal-to-noise ratio SNR is then defined as the ratio of the power spectrum at the driving frequency f=ω/2π to the value of the background power spectrum at that same frequency [25]. Then, one of the classic ways to characterize quantitatively the stochastic resonance phenomenon is to plot SNR as a function of the noise intensity. Such a plot is illustrated in Figure 8(e) and exhibits the signature of stochastic resonance phenomena: a local maximum in the signal-to-noise ratio for a nonzero noise intensity. Qualitatively, the presence of such a maximum is expected when the noise-induced average escape rate from a potential well (the “Kramers” rate) matches twice the frequency of the external signal: exp⁡(-ΔΦ/s)≅22πf [24].Figure 8
(a–c) Time series illustrating the evolution of a variabley in a double well potential, subjected to a weak periodic signal and to white noise fluctuations for a noise intensity of (a) σ = 0.05, (b) σ = 0.1, and (c) σ = 0.15. (d) Power spectrum of (a). (e) Signal-to-noise ratio as a function of the noise intensity showing a maximum at σ~0.11. (f) Residence time distribution for the region y>0 for σ = 0.05.
(a)
σ = 0.05 (b)
σ = 0.1 (c)
σ = 0.15 (d)
σ = 0.05 (e) (f)Another feature of stochastic resonance is shown by plotting the residence time distribution. This is a histogram of the time between transitions. One example is illustrated in Figure8(f) for σ = 0.05. The distribution is basically a decaying function modulated by periodical peaks centered at kP/2, where P is the period of the driving signal and k is an odd integer. This is easy to understand: suppose that a transition from one well to the other just occurred, most likely as the potential well is the lowest. Then, after a waiting time of half a period, the potential barrier will likely reach its minimum again and transitions back to the first well will be favored. If the system misses the first time, it will have to wait a full period before favorable conditions occur again, and so on.


## 5. Stochastic Resonance Applied to the Limestone-Marl Diagenetic Model

In this section, we apply the concept of stochastic resonance to our multistable model: we investigate the effect of a weak periodic external signal driving the model, superposed to a source of random fluctuations. The weak external periodic signal could represent a fluctuating climate proxy, such as the Milankovich cycle, which is believed to be instrumental in inducing ice ages [54, 55]. This cycle is characterized by variations in the earth’s orbital parameters, such as the precession of the earth’s axis (about 20 ka), changes in its obliquity (about 40 ka), and variations in the earth’s eccentricity (about 100 ka). As a result of relatively large variations in solar insolation induced by the first two cycles, it is not surprising that the associated climate response is strong. However, based on simple energy budget arguments, the climate response to variations in eccentricity is expected to be weak [56].Not surprisingly, if the precession and obliquity forcing is strong enough, a periodic small-scale sedimentary sequence is generated. But, in order to investigate the potential impact of the stochastic resonance mechanism on rhythmite formation, we rather consider the weakest orbital forcing term in our model: the eccentricity signal at 100 ka. Notwithstanding the fact that postdepositional compaction is not modeled in our system, a pattern generated at that period would correspond to large-scale variations in the thickness of layer bundles [31].Concretely, we assume that the small eccentricity driving signal translates in weak periodical variations in the sediment deposition conditions. For simplicity, we keep all model parameters constant as well as the chemical composition of the incoming sediment but vary its structural properties, that is, its porosity. We thus modify our diagenetic model by using the following time-dependent porosity boundary condition:(50)ϕ0t=ϕ0+αsin⁡ωt+εt,where ϕ0 is the constant background porosity, α is the amplitude of the small periodic driving signal, and ω is its frequency. The stochastic term ε(t) corresponds to the noisy fluctuations that couple the system to the environment. We choose an Ornstein-Uhlenbeck (OU) [51] colored process to characterize the dynamics of the noise. It is more realistic than white noise and is described by two parameters, a noise intensity parameter σ and a noise correlation time τ, as follows:(51)dεdt=-ετ+ξtτ;ξt=0;ξtξt′=2σδt-t′.The OU process describes a stochastic variable ε(t) that includes a linear damping term acting over a time scale τ and that is itself driven by white noise with intensity σ. The long-time correlation of the OU process can easily be obtained:(52)εtεt′≅στexp⁡-t-t′τ.It is seen that the correlation decays over a time scale τ and, in the limit where τ approaches zero, the OU process approaches a white noise process of intensity σ. The solution of the linear equation (51) at time t+Δt can be exactly related to that at time t [57]:(53)εt+Δt=εtexp⁡-Δtτ+γtστ1-exp⁡-2Δtτ1/2,where γ(t) is a random number sampled from a normal Gaussian distribution of variance one and average zero. Given an initial value ε(0), the numerical implementation of the OU process is straightforward with the application of (53) at every time step Δt. We chose ε(0) to be sampled from the equilibrium distribution, satisfying (52): ε(0)=γ(0)σ/τ.Figure9 shows various features of the numerical solution of the limestone-marl system (see (40)–(43)) with a stochastically driven porosity boundary condition as in (50) and (51). The case corresponding to scenario 2 was adopted (Table 1 and Figure 7) with ϕ0=0.68 and ϕin = 0.6. The other parameters were α = 0.03, ω = 2π/105 rd/a (corresponding to a driving frequency f = 13.19 in reduced unit), and τ = 0.1 and σ was varied. Deterministically (σ=0, Figure 9(a)), the forcing amplitude α is sufficiently small for the system to oscillate slightly about the corresponding 0-C steady state without transiting to other steady states. But for a nonzero noise intensity, noise-induced transitions to other steady states (0-0 or A-0) occur (Figures 9(b)–9(d)). For intermediate values of the noise intensity (Figure 9(c)), transitions occur in stochastic synchrony with the external periodic signal. A power spectrum of the calcite composition is shown in Figure 9(e) for σ = 0.001, as well as a 0-C state residence time distribution for σ = 0.00125 (Figure 9(f)). The latter graph is qualitatively similar to the case of stochastic resonance in a double well potential (Figure 8(f)).Figure 9
(a) Dimensionless time series illustrating the deterministic evolution of the calcite compositionCc and the terrigenous composition CT at the bottom of the sediment when the system is subjected to weak periodic driving of period 105 a. The parameter values are those of case 2. The transients are omitted. (b–d) Calcite composition Cc when the system is also subjected to colored noise fluctuations for a noise intensity of σ = 0.0001 (b), σ = 0.001 (c), and σ = 0.005 (d). (e) Power spectrum of (c). (f) 0-C residence time distribution for σ = 0.00125. (g) Signal-to-noise ratio as a function of the noise intensity.
(a)
σ = 0 (b)
σ = 0.1 × 10−3 (c)
σ = 1 × 10−3 (d)
σ = 5 × 10−3 (e)
σ = 0.001 (f)  (g)Again, we can define a signal-to-noise ratio SNR as the ratio of the power spectrum at the driving frequencyf to the value of the background power spectrum at that frequency. Figure 9(g) shows a plot of the SNR for various noise intensities σ. The presence of peaks in this graph is a signature of stochastic resonance. In fact, multipeak stochastic resonance has been observed in various systems. In a chain of Fermi-Pasta-Ulam oscillators driven at the boundary by noise [58] and a weak periodic signal, multipeaks in the SNR plot were linked to the coexistence of various stable and metastable stationary regimes. Multipeaks were also seen in an array of coupled monostable nonlinear oscillators [59] and in a quartic bistable potential subjected to two colored noises [60]. In the latter case, a transition from one-peak to two-peaks stochastic resonance was observed as the correlation time of one of the colored noise increased. In our case, the presence of colored noise (with a fairly high correlation time) and of multistability in the deterministic system could both play a role. But further analysis is beyond the framework of this presentation.


## 6. Conclusion

In this contribution, we propose a simple reactive-transport model for the formation of limestone-marl sequences, implementing the conceptual ideas of Munnecke et al. [27] and Böhm et al. [28] in a mathematical framework. The reactions involved include calcite dissolution and precipitation and aragonite precipitation. Aragonite dissolution occurs through bacteriological activity in a discrete zone located at a fixed distance below the sediment surface. Solute diffusion in a porous environment and sediment compaction are considered. Self-organized oscillating solutions were obtained for some parameter values. The system typically exhibits multistability between various steady states or limit cycle solutions.Coupling the system to a weak external periodic signal superposed to a random contribution presents the necessary ingredients for stochastic resonance to occur in the system. For illustration purpose, we chose a weak 100,000 years Milankovitch signal that affects the physical properties of the incoming sediment (porosity). We found that, for a noise strength sufficiently high, switches between a calcite-rich state and a terrigenous-rich one are obtained, in statistical synchrony with the weak external signal.The analysis presented here is only preliminary. Other parameter values need to be explored; the coupling between the system and its environment does not need to involve only the surface sediment porosity and other types of stochastic resonance patterns could be obtained, particularly in the regime where the system is bistable between a steady state and a limit cycle. But these findings open the way to a more thorough understanding of the nontrivial effects of noise on the formation of limestone-marl sequences influenced by a weak external signal, such as the Milankovitch 100,000 years’ cycle.



---

*Source: 4968315-2018-02-28.xml*