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.
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 (Figure 1). 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].
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.
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 Section 2 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 the -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 and , respectively. Within the aragonite microbially induced dissolution zone, A dissolves and releases calcium Ca+2 and carbonate ions 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 components = as the mass of per total solid dry sediment mass . Then, by definition,The mass of component per volume of dry sediment is , where is the total dry sediment density:Here, is the volume occupied by solid and is its density. The mass of component per total sediment volume is . The diagenetic equation for the solid component reads as follows:where is the velocity of the solid sediment matrix (counted positive downwards) and is the net reaction rate of component reported per unit sediment volume.
Dividing (3) by , summing over all components , and using (2) result in a continuity equation for the porosity:
We now define the concentration of dissolved species as , the number of moles of per unit pore water volume (in M units). The corresponding diagenetic equation iswhere is the velocity of the pore water (counted positive downwards), is the (tortuosity-corrected) diffusion coefficient of component , and is the net reaction rate (moles of component reacted per unit time and sediment volume). We will use the following tortuosity factor, often used in diagenetic modeling [33]:where is the diffusion of component in seawater.
Equation (5) can also be applied to water itself, with , where is the density of seawater (considered constant) and is its molar mass. We thus get another version of the continuity equation for porosity:Adding (4) and (7) givesThe reactive term describes the net change in water volume fraction as a result of precipitation/dissolution of the solid components and is nothing but . Thus,and (4) may be rewritten as follows:It shows how porosity is modified by compaction and the reactive kinetics of the precipitation/dissolution of minerals. Equation (9) indicates thatwhere is a position-independent factor.
Finally, using (10), the diagenetic equation (5) is simplified to
2.2. Solving for the Velocities
The conservation of momentum takes the form of Darcy law. In our situation, this iswhere is the hydraulic conductivity, is pressure, and 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]:where is the average grain diameter and the viscosity of water. We will collect all the constant terms in a prefactor and use insteadHere, β is an empirical factor with units of velocity and is a correction factor that becomes important for close to unity. After Hsu and Cheng [36], it is expected thatThus, we choose [36]
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 in the sediment. Thus,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 [39] (the inverse is called the pore space compressibility in Mello [40]):The Darcy flow (see (13)) then becomes with the help of (11)The solid phase velocity at the water-sediment interface must match the sedimentation rate . Denoting as the porosity at the boundary (water-sediment interface) and as the dry sediment density at the boundary, we can eliminate and evaluate the solid matrix advection field asEquation (11) then gives the fluid advection velocity :
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 for “no reaction”) obeys an empirical relation [41–43] of the formwhere and are parameters. is termed compaction coefficient [44] or sediment compressibility [45]. Then, the elastic response function is
Substituting (22) in (10) shows that the pore compressibility generates a porosity diffusion-like term with a diffusion coefficient
As will be clear below, we will use the molecular diffusion coefficient of calcium to choose appropriate length and time scales. Estimates of (with 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 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 , where is the (homogeneous) initial porosity.
Summarizing, we have eight field variables: three compositions , two concentrations , the porosity ϕ, and the velocities and . 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 , we will use definition (1) to obtain .
2.3. Reactions
In this simple model, we use four reactions.
(i) Dissolution of aragonite (reaction DA), A → Ca+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:where is another rate coefficient (units of inverse time, ). The term is the mass of aragonite per solid volume and multiplying by the factor reports the rate per total sediment volume. The undersaturation factor iswhere is the solubility of aragonite and is a reaction order. It is understood that this reaction occurs only when the system is undersaturated, . Finally, the factor is a characteristic function that is zero when is outside the aragonite dissolution zone (ADZ) and one otherwise. The ADZ is characterized by the position of its top edge and its thickness .
(ii) Precipitation of aragonite (reaction PA): Ca+2 + → A. For simplicity, we neglect nucleation. The corresponding growth rate per unit sediment volume can be written as a first-order reaction:where is another rate coefficient (units of inverse time, ). The oversaturation factor iswhere is a reaction order (different from in general). It is understood that this reaction is only effective when the system is oversaturated, .
(iii) Dissolution of calcite (reaction DC), C → Ca+2 + . The reaction rate is modeled as a first-order one with respect to calcite:where is a rate coefficient (units of inverse time, ). The undersaturation factor iswhere is the solubility of calcite and is a reaction order. It is understood that this reaction takes place only when the system is undersaturated, .
(iv) Precipitation of calcite (reaction PC): Ca+2 + → C. Here,where is a rate coefficient (in ) and is a reaction order (different from in general). It is understood that the reaction occurs only when the argument of is positive.
The net reaction rates appearing in (3) and (5) are thereforeThe molar mass of aragonite (or calcite) appears since, by definition, 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 boundary , we will assume that the system has no diffusive flux. Mathematically,The initial conditions are chosen as spatially homogeneous constants:
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 equal to a constant defined by the initial composition of the sediment and let . Combining (3) and (4) givesWe can verify that summing these equations over leads to an identity, as it should.
Finally, (12) for the dissolved species becomes
2.6. Scaling
It is convenient to express the dynamics of the system in reduced (dimensionless) form. We scale the dissolved species concentrations with and the advection velocities with a typical sedimentation rate . We introduce time scale and length scale defined as follows:The evolution equations then take the following scaled forms in terms of the dimensionless variables , , ; , and (and dropping the ′):Here, the dimensionless parameters are defined as follows:Da can be interpreted as a Damköhler number. The under- and oversaturation factors becomeHere, is a characteristic function indicating that aragonite dissolution is activated in the ADZ only. It is defined as 1 when and zero otherwise. Finally, taking into account the arguments of Section 2.2, we neglect terms of order in the advection terms (except for the small porosity diffusion term). The scaled velocities in (40)–(43) are
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 in water at various temperatures are available [33]. The solubilities and and the kinetic parameters (aragonite dissolution) and (calcite precipitation) are taken from Mackenzie and Andersson [46]. For simplicity, we have taken and .
The geometrical parameters of the system are its length , the position below the water-sediment interface of the ADZ, and its thickness . The value of is compatible with the value mentioned in Böhm et al. [28]. The sedimentation rate at the interface is of the order of 0.1–0.01 cm/a [47].
The inverse of the rate coefficients and is chosen to be smaller than the characteristic age of the system , 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 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 choice b = 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 parameters , , , , , , and were varied, as well as the water-sediment interface conditions 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 ).
2.8. Numerical Approach
Starting with an initial uniform distribution for the five variables , , , , and ϕ, system (40)–(43) is solved numerically on a regular grid of spatial steps . The time is discretized in units of time steps . 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 or . 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 by Taylor expansion about . The advective equations (39) and (40) are solved using an explicit upstream scheme. When 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 .
3. Deterministic Solutions
Figure 3 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).
(a) a
(b) a
(c) a
(d) a
(e) Steady state
One 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, . 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 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. Figure 4 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 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)).
(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, Figure 5 shows the evolution of the sediment composition at the bottom of the system using two different initial porosities 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 defines the basin of attraction of attractor . For fixed , as the porosity at the top boundary is varied, the boundaries of the basins of attraction shift. Consequently, different steady states are also obtained when the parameter changes. Figure 6(a) illustrates the domain of existence of each attractor in the (, ) space for scenario A.
(a) ;
(b) = 0.55; = 0.72
One can understand the nature of the steady states as follows. A state 0-C (as illustrated in detail in Figure 3) corresponds to a situation where the advection velocity 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 . On the other hand, although the sediment at 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 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), 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 plays an important role in generating pulsating solutions. As is seen from (47), the absolute value of is higher for larger values of the boundary porosity and of the overall porosity, in agreement with Figure 6(a).
Figure 7 shows a phase diagram in (, ) 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.
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 system moving in a bistable symmetric potential :Here, is the weak amplitude of the external periodic signal, ω is its frequency, and 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 σ:where the symbol refers to averaging with respect to the noise realizations and is the Dirac-delta function. In the undriven deterministic case (), it is seen that the system is indeed bistable: two stable steady states coexist at and an unstable fixed point is present at . The potential barrier height separating the two basins of attraction and is . 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 .
Figure 8 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 , and we varied the noise intensity . It is seen that for small noise intensity (Figure 8(a)), transitions between and 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 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: [24].
(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 Figure 8(f) for σ = 0.05. The distribution is basically a decaying function modulated by periodical peaks centered at , where is the period of the driving signal and 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:where is the constant background porosity, α is the amplitude of the small periodic driving signal, and ω is its frequency. The stochastic term 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:The OU process describes a stochastic variable 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: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 can be exactly related to that at time [57]:where is a random number sampled from a normal Gaussian distribution of variance one and average zero. Given an initial value , the numerical implementation of the OU process is straightforward with the application of (53) at every time step . We chose to be sampled from the equilibrium distribution, satisfying (52): .
Figure 9 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 and = 0.6. The other parameters were α = 0.03, ω = 2π/105 rd/a (corresponding to a driving frequency = 13.19 in reduced unit), and τ = 0.1 and σ was varied. Deterministically (, 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)).
(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 frequency 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.
List of Mathematical Symbols
| : | Space-independent parameter in (11) (cm/a) |
| : | Sediment compressibility ((kPa)−1) |
| : | Mass of component per total dry sediment mass (-) |
| : | Initial homogeneous values of (-) |
| : | Values of at the water-sediment interface (-) |
| : | Molar concentration of dissolved species (M) |
| : | Initial homogeneous values of (M) |
| : | Values of at the water-sediment interface (M) |
| : | Average grain diameter (cm) |
| Da: | Damköhler number (-) |
| : | (-) |
| : | (-) |
| : | Molecular diffusion coefficient of ion in the sediment (cm2/a) |
| : | Molecular diffusion coefficient of ion in water (cm2/a) |
| : | Porosity diffusion coefficient (cm2/a) |
| : | Correction factor in the hydraulic conductivity (see (17)) (-) |
| : | Acceleration of gravity (cm2/s) |
| : | Response function (kPa) |
| : | Thickness of the aragonite dissolution zone (cm) |
| : | Hydraulic conductivity (cm/a) |
| : | Solubility of component (M2) |
| : | Rate coefficients (a−1) |
| : | Length of the system (cm) |
| , : | Reaction order for aragonite precipitation and dissolution (-) |
| : | Mass of component in the sediment (g) |
| : | Total mass of dry sediment (g) |
| : | Number of discrete spatial steps (-) |
| , : | Reaction order for calcite precipitation and dissolution (-) |
| : | Fluid pressure (kPa) |
| : | Reaction rate for solid component (g/a-) |
| : | Reaction rate for ion (mol/a-) |
| : | Reaction rates for precipitation and dissolution (g/a-) |
| : | Reaction rate for water (mol/a-) |
| : | Sedimentation rate (cm/a) |
| : | Time (a) |
| : | Time scale (a) |
| : | Solid matrix velocity (cm/a) |
| : | Scaled solid matrix velocity (-) |
| : | Volume of solid component (cm3) |
| : | Total volume of solid components (cm3) |
| : | Pore water velocity (cm/a) |
| : | Scaled pore water velocity (-) |
| : | Position from the water-sediment interface (cm) |
| : | Position of the top edge of the aragonite dissolution zone (cm) |
| : | Position scale (cm) |
| : | Dynamical variable in standard stochastic resonance (-) |
| : | Amplitude of the periodic driving signal (scaled) (-) |
| : | Constant in the hydraulic conductivity (see (15)) (cm/a) |
| : | Random number sampled from a normal Gaussian (-) |
| : | |
| : | Time step (-) |
| : | Spatial step (-) |
| : | Random variable in the Ornstein-Uhlenbeck process (-) |
| : | Viscosity of water (Pa-s) |
| : | Characteristic function for the dissolution zone (-) |
| : | Ratio of rate coefficients in (40)–(42) (-) |
| : | Molar mass of water (g/mol) |
| : | Molar mass of calcite or aragonite (g/mol) |
| : | Ratio of rate coefficients in (40)–(42) (-) |
| : | Random white noise process (-) |
| : | Density of component (g/cm3) |
| : | Total solid density (g/cm3) |
| : | Initial value of the total solid density (g/cm3) |
| : | Density of water (g/cm3) |
| : | Noise intensity (scaled unit) (-) |
| : | Effective stress (kPa) |
| : | Scaled correlation time of the Ornstein-Uhlenbeck process (-) |
| : | Porosity (-) |
| : | Initial homogeneous porosity (-) |
| : | Stationary unreactive system porosity in (24) (-) |
| : | Porosity at the water-sediment interface (-) |
| : | Parameter in (23) (-) |
| : | Bistable potential in standard stochastic resonance (-) |
| : | Frequency of the periodic driving signal (scaled) (-) |
| : | Oversaturation or undersaturation factors (-). |
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was funded by a grant from the Natural Sciences and Engineering Research Council of Canada. The author acknowledges André Desrochers for providing the picture in Figure 1.