Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 25, 2019;19:16]

Applied Mathematics and Computation xxx (xxxx) xxx

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations Daniel Tesfay a,b, Pingyuan Wei a, Yayun Zheng a,c, Jinqiao Duan d,∗, Jürgen Kurths e a

School of Mathematics and Statistics & Center for Mathematical Sciences, Huazhong University of Science and Technology, Wuhan 430074, China b Department of Mathematics, Mekelle University, P.O.Box 0231, Mekelle, Ethiopia c Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan 430074, China d Department of Applied Mathematics, Illinois Institute of Technology, Chicago, IL 60616, USA e Research Domain on Transdisciplinary Concepts and Methods, Potsdam Institute for Climate Impact Research, P.O. Box 601203, Potsdam 14412, Germany

a r t i c l e

i n f o

Article history: Received 31 March 2019 Revised 10 August 2019 Accepted 21 October 2019 Available online xxx Keywords: Maximal likely pathways Thermohaline circulation Lévy motion Stommel two compartment-model Fokker-Planck equation,

a b s t r a c t In this work we study the impact of non-Gaussian α -stable Lévy motion on transitions between metastable equilibrium states (or attractors) in a stochastic Stommel twocompartment model for the thermohaline circulation. By maximizing probability density of the solution process associated with a nonlocal Fokker-Planck equation, we compute maximal likely pathways and identify corresponding maximal likely stable equilibrium states. Our numerical results indicate that random fluctuations with small intensity induces a weakened thermohaline circulation when the Lévy noise stability index is from 0.1 to 0.7. Published by Elsevier Inc.

1. Introduction Nonlinear complex systems are in reality open to noisy perturbations. These noises, both Gaussian noise and nonGaussian noise, usually have important role on setting the dynamical behaviours of the systems. It is worth studying the influence of noise on global thermohaline circulation (THC) since it has a huge impact on the livelihood of human beings. This large-scale circulation in the ocean is decisive in the global climate change. The THC is driven by fluxes of heat and freshwater across the sea surface and posterior interior mixing of heat and salt [1]. Modeling the THC is a formidable challenge due to the lack of, for instance, universal equation of state which connects the density of water to temperature and salinity and the complicated form of the domain which is bounded by the edges of various continents. Taking the systemlevel approach, i.e., as a system the ocean is just a compartment filled with salt water where the circulation of the salt water solution is driven by the density (heat and salinity) difference, enables us make a headway. The oceanic conveyor belt (or THC) transports warm surface water, around 17 Sv in volume and 0.5-1.5 × 1015 W in heat on annual average, from low latitude (or equatorial) to high latitude (or polar) regime. The warm water from the low latitude regime cools down, gets denser and sinks on the high latitude regime. Then, the cold water at deeper levels, 2–3 km in depth, backs from the polar ∗

Corresponding author. E-mail addresses: [email protected] (D. Tesfay), [email protected] (P. Wei), [email protected] (Y. Zheng), [email protected] (J. Duan), [email protected] (J. Kurths). https://doi.org/10.1016/j.amc.2019.124868 0 096-30 03/Published by Elsevier Inc.

Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

ARTICLE IN PRESS

JID: AMC 2

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

to the equator regime as deep water-currents and surface water becomes more salty at the equator as water is removed by evaporation due to the greater heat. Evaporation makes the water heavier and thus less salty and as a consequence it upwells. To examine this circulation Stommel pioneered a simple deterministic two-compartment model [2]. Stommel’s idea has been further extended and improved through both conceptual and numerical models [3–5]. Actually compartment models are taken as the simplest setting to investigate stochastically forced systems. Climate records reveal incidence of abrupt climatic changes that might have some link with transition between stable equilibrium states of the THC [6]. Uncertain processes interrelated with random atmospheric fluctuations, challenges corresponding to unresolved scales and strange mechanisms worth consideration since they may lead to switching between the equilibrium states. It is customary to regard noisy fluctuations as Gaussian [7–10]. Velez-Belchi et al. [6] verified stochastic resonance with additive Brownian noise induces transitions between the different stable states of the THC. A pathwise analysis of slowly varying not too large a Brownian noise perturbation on compartment model for THC has shown the system spends significant amount of time in metastable equilibrium [11]. Multiple stable states and possible critical transitions between metastable states are discussed in [12]. All the aforementioned studies are restricted to Gaussian noise led perturbation. Noisy fluctuations, meanwhile, are demonstrated to be non-Gaussian in some complex systems such as spectral analysis of paleoclimatic data [13], metastability in climate systems [14], and bursty transition in gene expression [15]. The effect of Gaussian noise on THC has been well dealt with, while the influence of non-Gaussian stable Lévy noise on the dynamical behaviour of this system is scant. It is more expedient to model these random fluctuations which portray haphazard jumps, discontinuity, heavy tail distribution and bursting sample paths by a non-Gaussian (α -stable) Lévy motion. With respect to α -stable Lévy noise, Wang et al. [16] used numerical simulation of mean exit time to determine low variability salinity difference level domain. In this paper, we focus on the influence of non-Gaussian α -stable Lévy noise on stochastic compartment model for THC by investigating maximal likely attractors of the maximal likely pathways, i.e., evolution of the maximizer of probability density function in stochastic phase portraits as time increases. The subtlety of phase portraits for stochastic differential equations (SDE) can be supplemented by reconveying these geometrical entities in terms maximal likely pathways [17,18]. We will consider a simple model for THC as a scalar SDE in the following form:

dYt = f (Yt )dt + σ dLtα ,

Y (0 ) = Y0 ,

(1.1)

where f(Yt ) is a deterministic vector field, σ is the noise intensity, α is stability index and {Lt }, t ≥ 0 is a scalar non-Gaussian Lévy motion in Eqn. (1.1). The paper is structured as follows. Section 2 presents details of a stochastic Stommel 2-compartment model. Section 3 reviews maximal likely stable equilibrium states or maximal likely attractors. Section 4 discusses numerical results of computations on maximal likely pathways under different parameters, viz, stability index, noise intensity and (nondimensional) freshwater flux. 2. Model and method 2.1. Model A two-compartment model can be taken as the simplest setting to study the oceanic THC. In the model, the compartments represent the polar region and the equatorial region. Te and Se are the temperature and salinity in the equatorial compartment, and Tp and Sp are temperature and salinity of the polar compartment. Each compartment contains the same volume V of well mixed salt water solution. The temperature and salinity are uniform in each compartment but not necessarily the same in the compartments. The mass density of the salt water solution in each compartment depends on temperature and the degree of salinity. The compartments exchange thermal energy and salt water with one another via advective transport driven by the difference in densities between the two compartments. The exchange of freshwater in turn affects the salinity. Physically, this exchange occurs at the surface of equatorial and polar compartments. In the absence of advective transport, the temperature in the equatorial compartment relaxes to the local atmospheric temperature Tae = T0 + θ2 according to a time scale tr . Here θ is the average equator-to-pole atmospheric temperature difference and T0 is a reference temperature. For simplicity, we assume that the temperature in the polar compartment relaxes to its local atmospheric temperature Tap = T0 − θ2 according to the same time scale tr . Finally, we assume that in the absence of advective transport the flux 12 FS of freshwater changes salinity in the equatorial and polar compartments F

F

S S at constant rates 2H and − 2H , respectively. The depth of the model ocean is denoted by H. If FS is positive (negative), salinity increases (decreases) in the equatorial compartment while salinity decreases in the polar compartment with the rate proportional to FS S0 [19]. The fact that freshwater flux terms have opposite signs but are of equal magnitude is briefly explained as follows. Over the long term, evaporation is dominant over equatorial compartment whereas rainfall dominates in the polar compartment [7]. As a consequence of conservation of mass, the net of rainfall and evaporative flux for the total system is zero. A sketch of a variant of Stommel’s model is shown in Fig. 1. Following [7], we define the heat and salinity balances by the equations:



T˙e = −tr −1 Te −



T0 +

θ 2



−

1 Q (ρ )(Te − Tp ), 2

Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

3

Fig. 1. A variant of Stommel two-compartment model([20]).

 T˙p = −tr −1 Tp −

 T0 −

θ

 −

2

1 Q (ρ )(Tp − Te ), 2

1 FS S0 − Q (ρ )(Se − S p ), 2H 2 FS 1 S˙p = − S0 − Q (ρ )(S p − Se ). 2H 2 S˙e =

(2.1)

Transport of energy and saline solution between the two compartments is due to difference in their respective densities. Density of the saline solution is a nonlinear function of temperature and salinity. Since water expands and its density decreases with the increase in temperature; and water gets heavier as its salinity increases, its density increases with increasing salinity. We therefore use Taylor’s theorem to approximate the density in a compartment by the following linear equation of state:

ρ = 1 + βS (S − S0 ) − βT (T − T0 ), ρ0 where β T and β S are (positive) constant thermal expansion and haline contraction coefficients, respectively. The exchange of mass between the two compartments is modeled with the transition coefficient Q that is dependent only on the density gradient ρ . There are several models of the transition coefficient [7]. In this article we choose, for simplicity, a positive definite transport function which is independent of the direction of ρ . This is viable because in the closed advection of water between two compartments (Q), the quantity of equatorial water moving into the poles and vice versa remains unaltered in both directions of circulation. Thus,

1 q Q (ρ ) = + td V



ρ ρ0

2 ,

where td is diffusive mixing time scale between the two compartments that would occur when there is no density difference. Subtracting the equations in (2.1) and defining S ≡ Se − S p and T ≡ Te − Tp yields a coupled pair of equations for the time evolution of temperature and salinity difference between the compartments:

d T = −tr −1 (T − θ ) − Q (ρ )T , dt d S FS = S0 − Q (ρ )S. dt H If we define dimensionless variables x ≡ θT , system of equations.

(2.2) y≡

S β θ ( βT )

τ≡

t td

, we recast (2.2) as the following coupled dimensionless

s

dx = (−β (x − 1 ) − x[1 + μ2 (x − y )2 ] )dτ , dy = (F − y[1 + μ2 (x − y )2 ] )dτ .

(2.3)

Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

ARTICLE IN PRESS

JID: AMC 4

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

Fig. 2. (Color online) The potential V as a function of the salinity difference y.

Where β =

td tr

measures temperature restoring tensility, μ2 =

qtd (βT θ )2 V

is strength of the buoyancy-driven convection β S t

between the two compartments relative to the diffusive mixing. The parameter F = βS θ0Hd FS represents dimensionless freshT water forcing. The physical meaning of β  1 is that mixing between the compartments temperatures is sluggish compared to the rapid parity of each compartment’s temperature with its local temperature forcing. Salinity leads to non-linearity which causes the existence of multiple equilibria and thresholds in the THC. Due to the dominance of the β term in Eq. (2.3) the temperature difference x changes very little, and we can thus take x as x = 1 + O (β −1 ). This leaves us with the first order differential equation in y(t) (where τ is replaced by the usual notation of time t for convenience):

dy = (F − y[1 + μ2 (1 − y )2 ] )dt.

(2.4) 1+μ2

2 μ2

μ2

If F = F¯ is constant, Eq. (2.4) can be written as dy = −V  (y )dt where V (y ) = −F¯ y + 2 y2 − 3 y3 + 4 y4 . The strength of the transport between compartments is governed by density difference |ρ| = ρ0 |βS S − βT T | which in dimensionless form is |y − x|. In certain parameter ranges Eq. (2.4) has two stable equilibrium states separated by an unstable equilibrium. In one of the stable equilibrium states, the salinity difference is small and corresponds physically to relative large northward heat transport. This state is usually referred to as the on-state. In the other equilibrium state the salinity difference is large and corresponds to weak (or even reversed) circulation. This state is called off-state. Depending on the initial value of the salinity difference, the trajectory of y may lie either in the domain of attraction of the on-state or in the domain of attraction of the off-state. Using Cessi’s valuation of parameter values [7,19] F¯ = 1.1, μ2 = 6.2, V(y) is a double-welled potential with two stable minima at y− = 0.24 and y+ = 1.07 and an unstable maximum at yu = 0.69 as shown in Fig. 2. Numerical simulations of the mathematical model show that the properties of the solutions are sensitive to changes in values of the parameters, especially to changes in values of the freshwater flux [21,22]. F¯ determines stable states and transition between them. A multiplicity of models has displayed that the competition between thermal and saline forcing results in multiple equilibria. Clearly, the deterministic system has two stable steady states and the random driving moves the temperature and the salinity away from initial equilibrium. So a perturbation is certainly needed to induce a transition between both states. These perturbations are provided by background noise such as the freshwater flux. The fluctuations in the freshwater flux or in the external salinity inputs lead F¯ to vary following the difference between salinity increment rate at the equator and corresponding decline at north pole. Fluctuations in precipitation, evaporation or ice sheet melting in the high latitudes could possibly be taken as the required perturbations that move the THC over a tipping point from on(off)-state to off(on)-state. Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

JID: AMC

ARTICLE IN PRESS

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

5

Fig. 3. (Color online) Sketch of bifurcation diagram for the salinity difference y, as a function of freshwater forcing strength F¯ . The upper and lower branch blue lines mark the two stable steady state solutions, On and off; the red dotted line marks the unstable solution connecting the two stable ones. The range of F¯ -values between F¯ ≈ 0.9556 and F¯ ≈ 1.2963 is the multiple equilibria regime, where the stable states coexist under the same initial conditions. Changing F¯ further, the deterministic Stommel two-compartment model (2.4) with μ2 = 6.2 and F¯ as control parameter jumps from one solution to the other at the two bifurcation points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

At such a tipping point a sudden and (possibly) irreversible abrupt change occurs. The fold bifurcation is the simplest bifurcation in which an equilibrium of a nonlinear dynamical system can lose its stability. (See the bifurcation diagram in Fig. 3). Critical thresholds for these transitions correspond to catastrophic bifurcations. As we approach a catastrophic bifurcation, from which a nonlinear system like the oceanic THC will experience a finite jump to an alternative state, the current attractor is located within a shrinking base of attraction. Increasing the freshwater forcing reduces the polar-equator density difference which determines the overturning rate, while the equator-polar salt advection counterbalances this by the overturning circulation. A critical transition occurs when F¯ passes a bifurcation point. Since there are two control parameters, F¯ and μ2 , it is desirable to consider beyond codimension-one bifurcation events. The parameter space is now represented jointly by the (F¯ , μ2 )-plane. A greater perspective to the fold bifurcation is shown in Fig. 4 in terms of the cusp catastrophe surface. This cusp codimension-two phenomenon is observed because we have independent control of both F¯ and μ2 . Here, we can consider fold curves whose projections separate the control plane into regimes displaying either one or three coexisting solutions. These fold curves coalesce and vanish at the cusp. Clearly, from the shape of the catastrophe surface, we can see that the bifurcation diagrams change qualitatively from well-defined double well potentials for large values of μ2 to less welldefined double wells for μ2 near the cusp, and finally to a single well potential for small values of μ2 . In addition to the catastrophe surface, a sequence of 3 bifurcation diagrams, corresponding to 3 vertical slices of the catastrophe surface, would help display these cases. When μ2 is slightly greater than its value at the cusp of the catastrophe surface, we get a bifurcation diagram that looks like the one shown in Fig. 5(a). The dynamics of the process Y(t) (see Fig. 5(b)) will be very sensitive to noise in this case, but it offers the possibility of randomly hopping back and forth between the wells which could be explored and compared with the results obtained from the stochastic resonance studies. In the presence of high level noise, a jump between points in the cusp bifurcation surface can be understood as purely noise induced jump. Jumps in salinity are said to be abrupt transitions. Abrupt environmental phenomena in the oceanic circulation such as large bursts of freshwater from the North Atlantic due to massive iceberg releases can not be modeled by stochastic continuous model. Introducing a jump process into Eq. (2.4) makes the model more close to reality. In [13], it was established that the fast time scale perturbation includes Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

JID: AMC 6

ARTICLE IN PRESS

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

Fig. 4. (Color online) The cusp catastrophe surface of the deterministic Stommel two-compartment model in (2.4) with F¯ and μ2 as control parameters.

Fig. 5. (Color online)(a) A bifurcation diagram near to the cusp where F¯ ≈ 0.925, μ2 = 3.4. (b) A double shallow well potential with local minima around y = 0.5 and y = 0.9

a component with an α -stable distribution where moments of order less than α exist and tails behaved to be fat hinting that extreme events have maximal likelihood of occurrence. α -stable noises might occur in complex systems with many different time scales where dynamics becomes strongly recurrent. In the oceanic circulation, freshwater flux may be paramt eterized as a sum of time independent mean component, that is F¯ , and a noisy fluctuating process dL . As a result, we come dt out with the SDE:

dYt = (F¯ − Yt [1 + μ2 (1 − Yt )2 ] )dt + σ dLt .

(2.5)

Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

7

2.2. α -Stable lévy process A stable distribution Sα (μ, β , γ ) is the distribution for a stable random variable [17], where the parameters α ∈ (0, 2 ), μ ∈ [0, ∞ ), β ∈ [−1, 1] and γ ∈ (−∞, ∞ ) respectively are indices of stability, scale, skewness and shift. An α -stable scalar Lévy motion Ltα , α ∈ (0, 2 ) is a non-Gaussian stochastic process satisfying the conditions:[17,23–25]

(i) Lα = 0, almost surely; 0 (ii) Independent increments: For each n ∈ N and each 0 ≤ t1 < t2 < · · · < tn−1 < tn < ∞, the random variables Ltα − 2 Ltα , · · · , Ltαn − Ltα are independent; n−1

1

α α (iii) Stationary increments: Ltα − Lα s and Lt−s have the same distribution Sα ((t − s ) , 0, 0 ); (iv) Stochastically continuous sample paths: Sample paths are continuous in probability. For all δ > 0, all s ≥ 0; P(|Ltα − Lα s | ) > δ ) → 0 as t → s. Jump process with a specific size in the generating triplet (0, 0, ν α ) of the random variable Sα can be defined by Ltα = α Lt − Ltα− < ∞, t ≥ 0 where Ltα− is the left limit of the Lévy motion Ltα in Rn at any time t. The jump measure of Ltα is 1

να (dz ) = Cα |z|−(1+α ) dz, where Cα =

1+α α √ ( 2 ) ) 21−α π (1− α 2

is stability constant. (For details see [17] sec. 7.2.2.)

Brownian motion, Bt , is a symmetric 2-stable process. Bt has independent and stationary increments, continuous sample paths almost surely and has normal distribution N (0, t ). 2.3. Method We will consider a simple model for THC as a scalar SDE with additive noise, i.e. the intensity of the noise is independent of the state of the THC, in the form

dYt = f (Yt )dt + σ dLtα ,

Y ( 0 ) = y0 ,

(2.6)

where f (Yt ) = F¯ − Yt [1 + μ2 (1 − Yt )2 ] is a deterministic vector field, σ is the noise intensity and {Ltα }, t ≥ 0, is a symmetric scalar α -stable non-Gaussian Lévy motion. The generator A of the solution process Yt for the SDE (2.6) with triplet (0, 0, ν α ) is

Aϕ (y, t ) = f (y )ϕ  (y, t ) +



R1 \{0}

= f (y )ϕ  (y, t ) + σ α

[ϕ (y + σ z, t ) − ϕ (y, t )]να (dz )



R1 \{0}

[ϕ (y + z, t ) − ϕ (y, t )]να (dz ),

(2.7)

where ϕ (y ) ∈ C (R1 ) in Eqn (2.7). The Fokker-Planck equation of (2.6) in terms of the probability density function p(y, t) for the solution process Yt given initial condition Y0 = y0 is [17]:

p(y, 0 ) = δ (y − y0 ),

pt (y, t ) = A∗ p(y, t ),

In Eq. (2.8), δ represents the dirac function and by solving



R1 \{0}

Aϕ (y )υ (y )dy =



R1 \{0}

A∗ ,

(2.8) the adjoint operator of A in the Hilbert space

L 2 ( R1 ),

can be computed

ϕ (y )A∗ υ (y )dy,

for ϕ , υ in the domains of definition for the operators A and A∗ , respectively,

A∗ υ ( y ) = σ α



R1 \{0}

[υ (y + z ) − υ (y )]να (dz ).

Consequently, we obtain the following nonlocal Fokker-Planck equation

pt = −( f (y ) p(y, t ))y + σ α



R1 \{0}

[ p(y + z ) − p(y )]να (dz ).

(2.9)

In the limiting case the stability index α = 2, the SDE (2.6) reduces to the Itô-Stochastic differential equation

dYt = f (Yt )dt + σ dBt ,

Y ( 0 ) = y0

(2.10)

where Bt is Brownian motion. The Fokker-Planck equation for (2.10) is

pt = −( f (y ) p(y, t ))y +

1 2 σ p(y, t )yy , 2

p(y, 0 ) = δ (y − y0 ).

(2.11)

Numerical finite difference method as in [26] and standard difference method are used in simulating Eqs. (2.9) and (2.11) respectively. Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

ARTICLE IN PRESS

JID: AMC 8

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx Table 1 Parameters of the stochastic Stommel compartment model [19]. Parameter

Meaning

Value

Unit

tr H td ta q V

temperature relaxation time scale mean ocean depth diffusion time scale advective time scale transport coefficient ocean volume thermal expansion coefficient haline contraction coefficient reference salinity meridional temperature difference

25 4500 180 29 1.92 × 1012 300 × 4.5 × 8, 250 10−4 7.6 × 10−4 35 25

days m years years m3 s−1 km3 K −1 g(kg)−1 K

βT βS S0

θ

3. Maximal likely trajectories Phase portraits provide geometric pictures for lower-dimensional deterministic dynamical systems. For low-dimensional stochastic dynamical systems such as Eqs. (2.5) and (2.6), analogous geometric pictures can be generated by plotting trajectories of the maximizer ym (t) of the probability density function p(y, t|y0 ) of the solution process. These trajectories represent the maximum likelihood path of the process y(t). At a given time instant t, the maximizer of ym (t) of the probability density function p(y, t|y0 ) yields the most probable location of this trajectory at time t. The deterministic entity ym (t) follows the uppermost crest of the surface in the (y, t, p)-space leaving behind trace out of the so called the maximal likely trajectory starting at y0 , as time goes on. A state that attracts all nearby states is called a maximal likely stable equilibrium (MLSE) state. MLSE states are dependent on non-Gaussianity index α , noise intensity σ and freshwater flux F¯ . The most probable phase portrait [17,18], an extension of phase portraits to envision stochastic dynamics, is the state space with representative maximal likelihood pathways including stable equilibrium states. Most probable phase portraits, like the phase portraits are geometric entities. They depict sample pathways and stable equilibrium states that are maximal likely. We examine the number and location of MLSE salinity difference states as a parameter varies. We define bifurcation time as the time between the change in number of maximal likely equilibrium states. It is a time scale for the birth of a new most probable stable equilibrium state. 4. Result and discussion In this section, we exhibit the number and value of MLSE states for the stochastic THC model (2.6). For simplicity, we consider four maximal likely trajectories with initial points y0 = 0.2, 0.4, 0.8, 1.5, to both sides of the equilibria in the deterministic dynamical system. Due to the lengthy computation process time is capped to T = 50. We then examine maximal likelihood trajectories when system parameters change followed by a brief discussion on general behavior of the system solution. 4.1. Result As Lévy noise parameters vary: The deterministic dynamical system is bistable in some range of freshwater flux parameter that contains F¯ = 1.1. A perturbation with very small noise intensity on the deterministic system invokes an interesting phenomenon as the Lévy noise index varies. For σ = 0.001, in line with the existence of a transient state the interval of alpha stability can be partitioned in two. (i) 0.1α < 0.7: In this interval of stability index, we have low salinity difference transient MLSE state. The transient state terminates since maximal likelihood salinity difference trajectories in the stochastic THC model merge by jumping to maximal likely metastable equilibrium state as time counts on. Bifurcation time is a time interval for the emergence of a new MLSE state. This time scale predicts transition from low salinity difference maximal likelihood state (or strong THC) to high salinity difference maximal likelihood state (or weak THC). Jump from low salinity difference comes out at different bifurcation times as index of stability to high salinity difference maximal likely equilibrium state in the interval varies. As can be seen from Fig. 6(a), bifurcation time attains its minimum at t ≈ 15.4683 when α ≈ 0.15. This duration should be multiplied by the diffusion time, taken here to be td = 180 years as in Table (1). Corresponding dimensional value of this minimum time is therefore about 2784 years. As shown in Fig. 7(a), bifurcation time decreases first, but increases after it attains minimum as α increases in this interval. In Fig. 6(b) bifurcation time is at t ≈ 22.5279 (about 4055 years) when α = 0.3. The transient maximum likely equilibrium state lasts for relatively short time interval for small α before it ends jumping up to the maximal likely metastable equilibrium state. Next, we turn our focus on discussing why this jump takes place. An α -stable process moves predominantly by big jumps as α nears to 0, and small jumps override as α is approaches 2. The jumps that occur in the maximal likelihood trajectories Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

JID: AMC

ARTICLE IN PRESS

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

9

Fig. 6. (Color online) MLSE salinity difference states 0.24 and 1.07 for F¯ = 1.1, σ = 0.001 and (a) α = 0.15. (b) α = 0.3. (c) α = 1. (d) α = 2.

Fig. 7. (Color online) The stochastic model experiences jumps from on to off state when F¯ = 1.1. (a) The system is monostable when α < 0.7, and the minimum bifurcation time is t ≈ 15.4683 when α = 0.15 and σ = 0.001. (b) For σ = 0.1 bifurcation time is t ≈ 7.8275 when α = 0.5.

due to the Lévy noise perturbation in this interval are seldom but big in size, as a result the salinity difference trajectories originating from a small neighbourhood of the deep well have a substantial probability of bypassing the barrier wall to land in the shallow well. It can also be clearly seen that the bifurcation time increases or equivalently transient MLSE state lasts longer, with the increase in α bound to the interval under consideration. Since y is dimensionless salinity difference, it is proportional to the salinity difference between the equatorial and polar regions, higher y values imply reduced density difference.Thus, higher values of y may consequently lead to diminished THC. This implies noise with such a small noise intensity and some range of values of stability indices may invoke weakened THC. Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

JID: AMC 10

ARTICLE IN PRESS

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

Fig. 8. (Color online) MLSE salinity difference states for σ = 0.01, α = 0.05 and (a) F¯ = 0.9. (b) F¯ = 1. (c) F¯ = 1.15. (d) F¯ = 1.4.

For the purpose of comparison, the influence of Gaussian noise with the same noise intensity is shown in Fig. 6(d). We notice that the stochastic system (2.10) is bistable. The reason for transition in the Lévy noise case may be attributed to the jump in the evolution of salinity difference pathways. (ii) 0.7α < 2: Our numerical results display the stochastic Stommel two-compartment model for THC in the SDE (2.6) has two MLSE states. The stable equilibrium states in the absence of noise match with these MLSE states after perturbation. The stochastic system in this interval of stability index is, therefore, not affected by the noise. In Fig. 6(c), for instance, the two most likely stable salinity difference are at 0.24 and 1.07. An implication of this phenomenon in the perturbed THC model is the possibility that there is neither a shift in potential wells from the original location nor a state transition. Intensity of the noise is not large enough to stir a slight change in the MLSE states and the jumps occurring in the maximal likely pathways are lower than the height of the barrier wall between the deep potential well and the shallow potential well (see Fig. 2). Numerical experiments reveal the range of stability index where maximum likely transient equilibrium state prevails gets wider as noise intensity increases from 0.001 to 0.1 as shown in Fig. 7. In Fig. 7(a), the system is monostable when stability index is below 0.7 whereas the system is bistable when stability index exceeds 1.9 in Fig. 7(b). Mean while, the bifurcation time is affected adversely by this noise intensity increment. The minimum bifurcation time decreases from 15.4683 (about 2784 years) in Fig. 7(a) to 7.8275 (about 1409 years) in Fig. 7(b). Notice that the stability index also increases from 0.15 to 0.5. This result hints that the transition from weak to strong THC may be enhanced as the Lévy noise parameters increase. As value of F¯ varies: F¯ is of particular interest as it monitors the strength of the freshwater forcing and determines the stable states and transitions between stable states of THC. In the absence of noise, the model is bistable in the interval 0.96  F¯  1.3, otherwise there is only one stable state: on-state for F¯  0.96 or off-state when F¯  1.3. It is desirable to examine how this stochastic THC model is influenced as freshwater flux parameter varies. We consider very small values of the parameters α = 0.05 and σ = 0.01 and the same initial conditions. When the value of the freshwater flux parameter varies while the system is under the influence of non-Gaussian Lévy noise, our numerical results show that, for F¯ = 1.1 there are two MLSE states at 0.24 and 1.07. It is worth noting that MLSE states of the stochastic system under the given freshwater flux parameter, stability constant and noise intensity remain the same as that of the deterministic system. When F¯ = 0.9, as in Fig. 8(a), the stochastic system has only one MLSE state at 0.17. Similarly, at 0.185 for F¯ = 0.95. Actually in these, the stochastic system is monostable and these maximal likely low salinity difference states correspond to strong THC. Existence of transient maximal likely equilibrium state is an interesting event that takes place when F¯ = 1, and F¯ = 1.15. In Fig. 8(b), the transient maximum likely high salinity difference equilibrium state 1 merges to the maximal likely low salinity difference stable equilibrium state 0.205 at the bifurcation time t = 34.9 (about 6282 years) when F¯ = 1. For F¯ = 1.15, the transient maximal likely equilibrium state Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

JID: AMC

ARTICLE IN PRESS

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

11

Fig. 9. (Color online) PDF’s for α = 0.1, 0.5, 1, 1.5. For small values of α the system is more influenced by large jumps. Pathways with initial points in the domain of attraction of y− have higher probability to transition from On-state to off-state. While when the values of α increase the perturbation is influenced by small jumps. Thus, these pathways have higher probability of staying in the potential well where they originated.

is at 0.27 and the maximal likely metastable equilibrium state is at 1.09. The jump is from high to low salinity difference state and it takes place at bifurcation time t = 21.8606 (about 3934 years) as in Fig. 8(c). When there is a transient state, numerical simulation results show bifurcation time increases and the transient maximal likely stable equilibrium state is low salinity difference state or strong THC which finally merges to high salinity difference state or weak THC when F¯ < 1.1, while if F¯ > 1.1, bifurcation time decreases and the transition is reversed. The system has maximum likely high salinity difference equilibrium state (or weak THC) 1.17 if F¯ = 1.4 as shown in 8(d). Abrupt climate changes may have some degree of relationship with sudden fluctuations of the THC due to vast freshwater flux caused by break up and melting of iceberg and glaciers. This conveyor belt transfers huge amount of heat from low to high latitudes regulating the global climate by redistributing heat across the planet. Weak THC may have far reaching consequences in water sources, precipitation for farming and well-being of low latitude ecosystems. On the other hand, it may lead to cooler ocean temperature in certain pockets of high latitudes and warmer temperature in other pockets with potential to extra sea level rise along coasts and ascends sea temperature adversely affecting the oceanic ecology. 4.2. Discussion In order to acquire an overall qualitative understanding of the statistical behavior of the solution of the stochastic differential equation (2.6) as parameters vary, we calculate a solution’s θ -moment numerically. For every t > 0, E|Ltα |θ of α -stable Lévy process {Ltα } is finite or infinite according to 0 < θ < α or θ ≥ α , respectively [25]. In Table 2, α is fixed would mean so is p(y, t|y0 ) and E|Ym(t ) |θ decreases with decrease in θ . For small α , p(y, t|y0 ) is also small. This indicates that probability of salinity difference actual pathways from an initial point y0 inside a domain of attraction of the on metastable state to transition to the off metastable state is relatively high. If α is large, on the other hand, small jumps have more influence on the system. As a consequence, transition is less probable. Figure 9 shows PDF of an ensemble of salinity difference Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

ARTICLE IN PRESS

JID: AMC 12

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

Fig. 10. (Color online) Around 50 simulations of system (2.6) for F¯ = 1.1, σ = 0.001 and (a) α = 0.15, with starting from the point y0 = 0.2. (b) α = 0.3, with starting from the point y0 = 0.5. (c) α = 1, with starting from the point y0 = 0.8. (d) α = 1.5, with starting from the point y0 = 1.5. Table 2 The θ -moment for the solutions of the stochastic differential Eq. (2.6) E|Ym(t ) |θ is finite if 0 < θ < α only.

α \E|Yt |θ \θ

0.05

0.25

0.5

0.75

1

1.25

0.1 0.5 1 1.5

5.1820 × 10−13 0.0012 0.0064 0.0074

– 0.0011 0.0055 0.0062

– – 0.0047 0.0051

– – 0.0041 0.0043

– – – 0.0037

– – – 0.0033

pathways starting from points near to y0 in the domain of attraction of the on metastable state (See Fig. 10). In agreement with the result in Table 2, transition is more likely for small α . Acknowledgments The authors would like to thank Xiujun Chen, Hui Wang, Xiaoli Chen and Ziying He for helpful discussions. In addition, we thank the anonymous referee for his or her detailed review which led to further improvement of paper. This work was partially supported by grants NSFC, China 11531006, 11771449, 11801192 and NSF, China 1620449. References [1] [2] [3] [4] [5] [6] [7] [8] [9]

H. Rahmstorf, Thermohaline circulation: the current climate, Nature 421 (699) (2003). H. Stommel, Thermohaline convection with two stable regimes of flow, Tellus 13 (1961) 224–230. H.A. Djiksta, W. Weijer, Stability of the global ocean circulation: basic bifurcation diagrams, J. Phys. Oceanogr. 35 (2005) 933–948. R.X. Huang, J.R. Luyten, H.M. Stommel, Multiple equilibrium states in combined thermal and saline circulation, J. Phys. Oceanogr. 22 (1992) 231–246. S. Rahmstorf, On the freshwater forcing and transport of the atlantic thermohaline circulation, Clim. Dyn. 12 (1996) 799–811. P. Vélez-Belchí, A. Alvarez, P. Colet, J. Tintoré, Stochastic resonance in the thermohaline circulation, Geophys. Res. Lett. 28 (2001) 2053–2056. P. Cessi, A simple box model of stochastically-forced thermohaline flow, J. Phys. Oceanogr. 24 (1994) 1911–1920. S.M. Griffies, E. Tziperman, A linear thermohaline oscillator driven by stochastic atmospheric forcing, J. Clim. 8 (1995) 2440–2453. G. Rong, L. Quan, Y. Yuangen, D. Haiyou, M. Chengzhang, J. Ya, Y. Ming, Noise decomposition principle in a coherent feed-forward transcriptional regulatory loop, Front. Physiol. 7 (600) (2016), doi:10.3389/fphys.2016.0 060 0.

Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868

JID: AMC

ARTICLE IN PRESS

[m3Gsc;November 25, 2019;19:16]

D. Tesfay, P. Wei and Y. Zheng et al. / Applied Mathematics and Computation xxx (xxxx) xxx

13

[10] L. Yonskai, Y. Ming, Z. Xiufen, The linear imterplay of intrinsic and extrinsic noises ensures a high accuracy of a cell fate selection in budding yeast, Sci. Rep. 4 (2014) 5764. [11] N. Bergund, B. Gentz, Metastability in simple climate models: pathwise analysis of slowly driven Langevin equations, Stoch. Dyn. 2 (2002) 327–356. [12] P. Good, J. Bamber, K. Halladay, A. Harper, L. Jackson, G. Kay, B. Kruijt, J. Lowe, O. Phillips, J. Ridley, M. Srokosz, C. Turlay, P. Williamson, Recent progress in understanding climate thresholds: ice sheets, the atlantic meridional overturning circulation, tropical forests and response to ocean acidification, Prog. Phys. Geogr. 42 (2018) 24–60. [13] P.D. Ditlevsen, Observation of α -stable noise induced millennial climate changes from an ice-core record, Geophys. Res. Lett. 26 (1999) 1441–1444. [14] L. Serdukova, Y. Zheng, J. Duan, J. Kurths, Metastability for discontinuous dynamical systems under Lévy noise: case study on amazonian vegetation, Sci. Rep. 7 (9336) (2017), doi:10.1038/s41598- 017- 07686- 8. [15] K. Niraj, S. Abhudai, V.K. Raul, Transcriptional bursting in gene expression: analytical results for general stochastic models, PLoS Comput. Biol. 11 (2015) e1004292. [16] X. Wang, J. Wen, J. Li, J. Duan, Impact of α -stable Lévy noise on the stommel model for the thermohaline circulation, Discret. Contin. Dyn. Syst. 17 (2012) 1575–1584. [17] J. Duan, An Introduction to Stochastic Dynamics, Cambridge University Press, Cambridge University Press, New York, USA, 2015. [18] Z. Cheng, J. Duan, L. Wang, Most probable dynamics of some nonlinear systems under noisy fluctuations, Commun. Nonlinear Sci. Numer. Simul. 30 (2016) 108–114. [19] H.A. Dijkstra, Nonlinear Climate Dynamics, Utrecht University, Cambridge University Press, New York, USA, 2013. [20] H.A. Dijkstra, Dansgaard-oeschger events[PDF], retrieved from https://www.whoi.edu/fileserver.do?id=243550&pt=10&p=124956, 2015. [21] J.M. Gregory, O.A. Saenko, A.J. Weaver, The role of atlantic freshwater balance on the hysteresis of the meridional overturning circulation, Clim. Dyn. 21 (2003) 707–717. [22] A. Timmermann, H. Gildor, M. Schulz, E. Tziperman, Coherent resonant millennial-scale climate oscillations triggered by massive meltwater pulses, J. Clim. 21 (2003) 707–717. [23] D. Applebaum, Lévy processes and stocahstic calculus, Cambridge Studies in Advanced Mathematics, 2nd, Cambridge University Press, Cambridge, 2004. [24] G. Samorodnitsky, M.S. Taqqu, Stable Non-Gaussian Random Processes, Chapman and Hall, New York, USA, 1994. [25] K. Sato, Lévy Processes and Infinitely Divisible Distributions, Cambridge University Press, USA, 1999. [26] T. Gao, J. Duan, X. Li, Fokker-planck equation for stochastic dynamical systems with symmetric Lévy motions, Appl. Math. Comput. 278 (2016) 1–20.

Please cite this article as: D. Tesfay, P. Wei and Y. Zheng et al., Transitions between metastable states in a simplified model for the thermohaline circulation under random fluctuations, Applied Mathematics and Computation, https://doi.org/10. 1016/j.amc.2019.124868