Interactive comment on “ Influence of extreme events modeled by Lévy flight on global thermohaline circulation stability

Abstract. How will extreme events due to human activities and climate change affect the oceanic thermohaline circulation is a key concern in climate predictions. The stability of the thermohaline circulation with respect to extreme events, such as fresh-water oscillations, greenhouse gas accumulations and collapse of the Atlantic meridional overturning circulation, is examined using a conceptual stochastic Stommel two-compartment model. The extreme fluctuations are modeled by symmetric α-stable Lévy motions whose pathways are cádlág functions with at most a countable number of jumps. The mean first passage time, escape probability and stochastic basin of attraction are used to perform the stability analysis of on (off) equilibrium states. Our results argue that for model with weak fresh-water forcing strength, the greatest threat to the stability of the on-state represents noise with low jumps and higher frequency that can be seen as civilization-induced greenhouse gas accumulation. On the other hand, the off-state stability is more vulnerable to the agitations with moderate jumps and frequencies which can be interpreted as wind- driven circulations towards higher latitudes. Under the repercussion of stochastic noise, on to off transitions are more expected in the model if the fresh-water influx is strong. Moreover, transitions from one metastable state to another are equiprobable when the fresh-water input induces a symmetric potential well.



Introduction
Natural and civilization catalyzed fluctuations in climate have significant impact on the ocean and ocean circulation pattern 15 variations greatly affect climate (Chapman and Shackleton, 1999;Clark et. al., 2002). The global thermohaline circulation (THC), known as great ocean conveyor as well, consists of ocean tides and drifts and loops that have extensive ramifications on climate. THC is basically an outcome of the interplay of fresh-water with thermal energy along with the ocean-atmosphere interface and inside the ocean competition of temperature and salinity (Rahmstorf, 2003). This enormous oceanic process has a significant contribution in maintaining the equilibrium of Earth's energy framework by redistributing thermal energy of the 20 order 1.2×10 15 W northwards in the Atlantic ocean. A large proportion of the meridional overturning circulation (MOC) is usu-ally categorized as THC because MOC takes the lion's share in this heat penetration to the north pole (Ganachaud and Wunsch, 2001;Trenberth and Solomon, 1994). The warm and saltier surface water on interannual and even longer time-scales gets freshened and loses heat to the cold atmosphere. Subsequently, the water descends slightly to the bottom of the Atlantic since it gets denser than the underneath water. The cooled water eventually returns southward as deep current and the warm temperature 25 around the equatorial belt opts for upwelling. This thermal forcing is deterred by the surplus of rainfall and water discharge over evaporation at the low latitudes which instigates a fresh-water flux opposing the descend in and around the north pole.
The global THC is a combination of the floating of deepwater currents around the equator and in the southern oceans, the horizontal currents, and the descending and forming of deepwater in high latitudes. Climate reconstructions indicate that conveyor belt has indispensable contribution in the climate system transitions from cold to warm or from warm to cold climate 30 states (Rahmstorf, 1995;Bond et.al., 1997;Grootes and Stuiver, 1997). Studying the stability of this oceanic conveyor belt by analyzing the influence of internal and external agitations on its dynamical behavior is increasingly pressing presently.
To study how the THC transports properties latitudinally, a conceptual deterministic two-compartment model was forwarded by Stommel (Stommel, 1961). Shreds of evidence from sea observation and model simulation show that the strength of the thermohaline flow is sensitive to the surface fresh-water flux fluctuation (Jackson and Wood, 2017;Caesar et.al., 2018).  petition between Thermal versus saline forcing competition can lead to a multiple equilibria regime if the relaxation time-scales for the temperature and the salinity are distinct. The thermohaline flow system is bistable, one with strong circulation (analogous with the present set up), and the second state with a very weak flow, when the salinity difference is forced by a prescribed fresh-water flux (Marotzke, 2000). The multistability of THC is also verified by results obtained from different numerical models (Broecker, 1987). The deterministic compartment model has been further extended to include noisy thermal and saline 40 forcing oscillations (Huang et. al., 1992;Rahmstorf, 1996;Djikstra, 2005).
Forcing the global general ocean circulation model with some particularly large stochastic fresh-water fluctuations is found to trigger pulsation of transport from one stable configuration to the other (Mikolajewicz and Maier-Reimer, 1990). The Gaussian noise perturbed THC has been under extensive study. For instance, it was shown in (Vélez-Belchí et. al., 2001) that an increment of 5% of the fresh-water forcing in the THC could stimulate transitions between a high and low salinity difference 45 metastable states. In a time-dependent compartment model for THC with Brownian motion and moderate noise intensity, hysteresis does not adiabatically follow stationary distribution (Bergund and Gentz, 2002). Meanwhile, the noise forcing climate comprises of a non-Gaussian α-stable Lévy noise component (Fuhrer et. al., 1993;Ditlevsen, 1999). The occurrence more than a dozen of additional Dansgaard-Oeschger (D-O) events that took place during the last glacial period could not have been reproduced by using the Brownian process. Other extreme events such as stock market crashes in the financial industry, strokes 50 and seizures, and earthquakes are not continuous perturbation processes. The jumps in those events could better be modeled by Lévy flights (Kuhwald and Pavlyukevich, 2016).
Paleoclimatic data indicate the coincidence of transitions from strong to weak or from weak to strong thermohaline circulation states with the occurrence of extreme climatic variations (Vélez-Belchí et. al., 2001). We analyze the stability of metastable states of the THC model against extreme events modeled by stable Lévy motion by calculating three quantities, namely, mean 55 first passage (exit) time, first passage (escape) probability and stochastic basin of attraction that carry the dynamical information of the model. In the present form of the THC, analyzing the intensity and mechanism of the forcing schemes that could trigger such transitions and studying the stability of strong THC and weak THC equilibrium states is of fundamental importance.
Particularly, we will study the effect of extreme events on the scalar stochastic THC model by measuring the stability of equilibrium states of the salinity difference process Y t for various values of (nondimensional) fresh-water forcing and non-Gaussianity parameter α. In Eq. (1), V is a potential function (details are given in Section 2).
The paper is structured as follows. In Section 2, we discuss the simplified conceptual stochastic Stommel two-compartment model for thermohaline circulation. A brief introduction of the stability analysis measures is provided in Section 3. Stability analysis of the stochastic thermohaline circulation system is given and results obtained are presented in Section 4. Our findings 65 are summarized in Section 5.

Oceanic thermohaline circulation model
Global thermohaline circulation is an oceanographic phenomena that refers to the movement of ocean waters across both hemispheres and is responsible for the heat transfer and redistribution, acting as a regulator of the global climate. The schematic functioning of THC is shown in Fig. 1. The main engine of this circulation is the difference in density between ocean currents 70 ∆ρ, which is determined by the salinity S e , S p and the temperature T e , T p of the water and can be represented by where β T = 0.17 × 10 −3 o C −1 is the thermal expansion coefficient and β S = 0.75 × 10 −3 psu −1 is the haline contraction coefficients, respectively. The surface ocean waters in the subtropical regions due to intense evaporation F s /2 have high salinity S e , however the high water temperature T e maintains the low density and prevent surface waters sinking.

75
In high latitude areas, the formation of dense water is mainly associated with lower temperatures T p and increased salinity S p due to the formation of ice. Thus, in the polar regions, the increase in surface water density causes it to sink and displace deep water. In this way, the origin of the thermohaline circulation is a vertical flow of surface water 1 2 Q(∆ρ), diving to an intermediate depth or close to the bottom, depending on the density of that water. The systems of superficial and deep circulation of the oceans are interconnected. The continuation is a horizontal flow: the recently sunk waters repel in the horizontal direction 80 the deep waters that occupied this place. In this way, the cold, dense waters sink and slowly flow towards the equator. Thermal energy and salinity balances can be defined by the system of differential equations (Cessi, 1994) (the dots represent derivative  (Cessi, 1994). Each compartment represents the waters of the equatorial and polar oceans with the same volumes V = 300 × 4.5 × 8, 250 km 3 and temperature Te/Tp and salinity Se/Sp characteristic of each of them. The other system parameters are the mean ocean depth H = 4500 m, the exchange mass function Q, the density gradient ∆ρ, the freshwater flux Fs, the equatorial atmospheric temperature Tae, the polar atmospheric temperature Tap, the reference temperature T0 = 5 o C, the reference salinity S0 = 35 psu, the meridional temperature difference θ = 25 K and the temperature relaxation time scale tr = 25 days. with respect to time): where  The time evolution of temperature ∆T ≡ T e − T p and salinity difference ∆S ≡ S e − S p between the compartments are obtained by subtracting the conservation equations (3), respectively.
The original system (4) with the substitutions x ≡ ∆T θ , y ≡ ∆Sβs θβT , τ ≡ t t d is reduced to the dimensionless system of evolution equations (Cessi, 1994;Tesfay et. al., 2020) where β the temperature restoration tensility, µ 2 the buoyancy-driven convection strength and F dimensionless fresh-water forcing are defined as The dynamical system (5) can be further simplified, since the diffusion time scale t d is much larger than the temperaturerestoring time scale t r . Thus, taking the approximation x = 1 + O(β −1 ), we get the first order differential equation in y(t) 105 where τ is replaced by the usual notation of time t for convenience. Considering that F is constant, Eq. (7) can be written as The variation in the fresh-water forcing strength F , Fig 2 ( is stable, called off -state and the other is unstable one. In the off -state y is large and corresponds to weak (or even reversed) circulation. Only stable off -state takes place in the THC systems with F > 1.2963, Fig 2 (c). With global warming, by the end of the twenty-first century, the mean surface air temperature due to harmful human activities leading to the accumulation of 115 greenhouse gases in the atmosphere will increase by 2−6 o C (Chapman and Shackleton, 1999;Ganachaud and Wunsch, 2001;Rahmstorf, 2003). This large-scale warming of the climate can cause an increase in frequency or/and intensity of extreme events on the time-scale of decades, such as high-latitude precipitations, the Greenland ice sheet deglaciation, the freshwater outflows to the oceans that certainly affect the dynamic aspects of thermohaline circulation (Clark et. al., 2002;Jackson and Wood, 2017;Gregory et. al., 2017). An increase in the mass of freshwater in the global ocean reduces its density and thereby complicates its 120 deep immersion, which can slow down THC and even "switch it off" if the parameters cross the tipping threshold (Stommel, 1961). Such extreme and fast events cannot be modeled by deterministic models, since they do not take into account the uncertainty, unpredictability and the likelihood of their occurrence. The Brownian motion predicts the behavior of such random fluctuation with very low accuracy, since it has continuous sample paths and normally distributed increments. It was proved (Ditlevsen, 1999;Marotzke, 2000) that the Lévy process, characterized by heavy-tailed distribution and discontinuous cadlag  paths, with a certain precision simulates and predicts these rare events. Consider now that the fresh-water flux can be written as sum of a stochastic component dL α t with a parameter F which is independent of time; then (7) generalizes into the Itô equation, Here L α t is a symmetric α-stable Lévy process, with α ∈ (0, 2), defined on the probability space (Ω, F, P ) (See Methods).

130
The drift term V ′ (Y t ) satisfies a Lipschitz condition with jump measures, then the solution of the SDE (9) exists and is unique (Applebaum, 2004). We focus our attention on the three THC models with different values of the parameter F ; a weak F =1, fresh-water forcing F =1.126 which induces a symmetric potential function and a strong F =1.28. These values represent different geometries of a double-well potential function as in Fig. 2 (d), (e) and (f).

135
This section summarizes definitions and the main properties of the stability measures used in the present analysis. Also a brief summary to the type of stochastic process chosen to model extreme events is given.

A symmetric α-stable scalar Lévy motion
Climate extreme events show random fluctuations having sample pathways with intermittent jumps and heavy tails. Natural and more appropriate candidate for modeling such a non-Gaussian process is an α-stable Lévy motion. L α t with 0 < α < 2 is 140 a stochastic process that satisfies the following properties: a) L α 0 = 0, almost sure; b) L α t has independent increments; c) stationary increments L α t − L α s and L α t−s have the same symmetric α-stable distribution, i.e. S α ((t − s) 1 α , 0, 0); 145 d) stochastically continuous sample paths, i.e., for every s > 0, L α t → L α s in probability, as t → s. The probability density function for L α t is defined by where f α is the probability density function for the standard symmetric α-stable random variable Y ∼ S α (1, 0, 0) (for more details see (Applebaum, 2004;Duan, 2015). The generating triplet of L α t is (0, 0, ν α ), with the jump measure, i.e. the expected 150 value of the number of jumps of size dz during the unit time, is defined as: where c α is the intensity constant. When α ∈ (0, 1), the α-stable Lévy motion has finite variation, otherwise, when α ∈ [1, 2) it is unbounded.

Mean first exit time 155
It is defined as the first exit time from a deterministic domain D ⊂ R 1 of attraction of y (on/of f ) as follows: and the mean exit time or the mean residence time of the process in the on(off)-state domain is denoted as u(y) Eτ (ω, y) ≥ 0.
It has been proven (Duan, 2015) that the mean exit time of the stochastic system (9) for an orbit starting at y ∈ D, satisfies the following nonlocal partial differential equation with an external boundary condition where A is the generator defined as Here D c is the complement set of D in R 1 . More over, the generator can be interpreted as Au = lim t→0 Eu(yt)−u t , for every 165 u ∈ C 2 (R 1 ).

Escape probability
The likelihood that the salinity difference process Y t exits firstly from the on(off)-state domain D by landing in the set U ∈ D c belonging to the off(on)-state domain is represented by 170 and solves the following integro-differential equation with the Balayage-Dirichlet boundary condition 0, y ∈ D c \U.

175
The evolution of the probability density function p(y, t) for the solution paths Y t of the SDE (9) is governed by a Fokker-Planck with the initial condition Y 0 = y 0 . Where A * is the adjoint operator of A (14) in the Hilbert space L 2 (R 1 ) that has the explicit form in the case of a symmetric α-stable Lévy motion. Thus, the equation (17) is specified as where f (y) = −V ′ is a vector field of the SDE (9).
Sufficient conditions for the existence and regularity of the probability density p(y, t) in the Lévy processes driven SDE (9) is established under Hörmander's condition by using the Malliavin calculus with jump. For more details, See (Chen et.al., 2015;Zhang, 2014) and the references therein.

Stochastic basin of attraction
The SBA is an important theoretical and practical tool that helps to describe a metastable behavior of a system. Stochastic basin of attraction (SBA) quantifies the stability of a metastable state in a dynamical system with noise perturbation in terms of size of the basin depending on the escape probability (Serdukova et. al., 2016). SBA is the collection of initial conditions of solution processes that have low (high) probability of exit (return) from (to) a neighborhood of an attractor. This geometric 190 tool is applicable to models with small noise and noise that is a function of an order parameter.
By Definition (Serdukova et. al., 2016): SBA of the attractor K with the open deterministic domain of attraction D is the set The stability analysis of the off and on-states in the THC model is based on the three main measures: escape probability, mean first exit time and stochastic basin of attraction which are frequently used in behaviour prediction of the stochastic dynamical system trajectories and whose main properties are summarized in Section 3. In this Section we discuss the results obtained and interpret them from the climatological point of view. The potential function V (y) (8) for the weak fresh-water input, as shown in Fig. 3 (a), is asymmetric and has the deepest on-state (y = 0.2) the widest 0.8 stability basin. The length of the 200 deterministic basin of attraction for off -state (y = 1) is 0.26 units smaller than that for on-state. This indicates that under the Lévy perturbations the transition from off -state to on-state is more likely. In fact, the off -state basin under the noise with α = 1 decreases by 6.75 times, while the on-state basin reduces only by 1.98. The greatest threat to the stability of on-state represents noise with low jumps of higher frequency, i.e. α = 1.5. Extreme events of short time-scale (decades or less) such as human-induced greenhouse gas emission are more likely to destabilize the on-state. As a response to increase in greenhouse 205 gas concentration, atmospheric water cycle boost is expected. This leads to a greater excess of precipitation over evaporation  in high latitudes, freshening and reducing the density of surface water, and hence its tendency to sink and weaken the THC.
The off -state stability, on the other hand, is more vulnerable to the perturbations with moderate jumps and frequencies, i.e α = 1. Centennial-scale wind driven circulation towards higher latitudes causes a strong salting of the tropical North Atlantic and subsequently enhances the recovery of THC.

210
The analysis of the trajectories' first mean exit time (12) and (13) from the deterministic basin Fig. 4(a) confirms the stability results based on the SBA.
Actually, the trajectories remain longer in the vicinity of on-state than in the off -state basin as the stability index (that is, α) increases. The relationship between the residence time and stability index can be understood as the extreme event modeled by Lévy motion with small jumps and high frequency (α = 1.5) contributes to faster escape compared to the events characterized 215 by moderate or big jumps that correspond to moderate or low probability (α = 1 or 0.5).
For strong fresh-water forcing the potential function V (y) is also asymmetric, but now with the reversed depth: off -state y = 1.14 here is deeper than on-state y = 0.37. In this model Fig. 3 (b), off -state is more stable than on-state and this implies that under the influence of stochastic noise, the transitions from on to off are more expected.
Comparing the dimension of the effects that noise with different α parameter values causes in the global thermohaline 220 circulation, it should be noted that the extreme events modeled by Lévy motion with α = 0.5 does not reduce the stability of the off -state significantly, inducing only 0.114 unit cutback in the basin width. In the case of weak THC stable state, large fresh-water outburst due to massive ice melting is rare but high impact event that would rather favour a total shutdown. The The symmetry in the potential V (y) is presented for F = 1.126 (see Fig. 3 (c)) generating the equality of the lengths of the on(off)-states stability basins. Also, under the force of stochastic noise, the transition from one state to another is equiprobable.
The greatest destabilizer of the states is the perturbation with moderate jumps and probability noise, alternatively the events 230 characterized by low probability but high jumps bring less transient impact between the states. In Fig 4 (b), the mean exit time of the orbits from the neighborhoods of on(off)-states is the same and manifests the highest values for the orbits carried by Lévy noise with α = 0.5. This shows as well that such type of perturbations has minimal influence on the stability of the states.
In Fig. 5 the D I and D C II SBA criteria (Serdukova et. al., 2016) based on the escape probabilities (15) and (16) for the stable on(off)-states of the THC model with different fresh-water forcing F are shown: (a)-(d) for weak fresh-water input, (e)-(h) for 235 symmetric potential well inducing fresh-water input and (i)-(l) for strong forcing. The first criteria defines the set D I of the initial salinity difference y 0 that originate the trajectories with the m < 0.3 probability of escape from the deterministic basin of on-state Fig. 5 (a), (e), (i) and off -state (d), (h), (l). The trajectories with probability M > 0.8 of return to the interval D I have their initial points y 0 included by the second criteria in the on-state's SBA Fig. 5 (b), (f), (j) (off -state's SBA (c), (g), (k)).
The geometry of the probability density functions (17) and (18) of the salinity difference process Y t Fig. 6, once again confirms 240 the previous conclusions about the stability of the on(off)-states. Therefore, the process is most likely to remain in the state that represents the deepest valley in the potential function V (y), this is in on-state for weak forcing and in off -state for strong forcing. For symmetrical potential well inducing forcing, the orbits with equal probability visit on or off -states. The transition between the states is more likely if the THC system undergoes to Lévy perturbations with low jumps and high probability noise as shown in Fig. 6 (c), since the difference between the probabilities of staying in each of the states is smaller. Therefore, for probabilities decreases with the increase in α parameter (which means a decrease in height and an increase in the frequency of jumps), resulting in a difference of 1.81 times for α = 1 Fig. 6 (b) and 1.42 times for α = 1.5 Fig. 6 (c).

Conclusion
Fluctuations in the THC patterns influence climate and civilization-induced and natural climatic changes also have impacts on 250 THC. The THC is bistable only for some domain of the nondimensional fresh-water parameter values. Extreme events such as greenhouse gas emissions, collapse of major ice sheets and global warming induce a transition between both states. It has, therefore, great importance to sufficiently scrutinize the stability of strong (weak) THC states against these extreme events. In this work, we have performed analysis of stability of the on-state and off -state in the stochastic two-compartment thermohaline circulation model by considering three different values of the fresh-water input control parameter.

255
Random noise agitations in geophysical complex dynamical systems have been regarded as continuous perturbation or Gaussian processes. However, the paleoclimatic data sets signify random fluctuations in swift climatic transitions have a non-Gaussian distribution with heavy tail and pathways that are cádlág functions with at most a countable number of jumps.
Therefore, it is more fitting applying non-Gaussian symmetric α-stable Lévy processes in modeling the influence of extreme events in the conceptual stochastic Stommel two-compartment model. Actually, these non-Gaussian processes are becoming 260 increasingly popular lately.
We applied three concepts, mean residence time, first passage probability and stochastic basin of attraction. Each of these quantities is helpful in understanding the stability of the strong THC (small salinity difference y) state and the weak THC (large salinity difference y) against stochastic perturbations (extreme events) and predicting transitions between these states.
The main conclusion of our work are the following: For the weak fresh-water flux, on to off -state transition is more probable 265 under extreme events modeled by Lévy noise perturbations with smaller jumps but with high frequency. This SBA analysis result is confirmed by calculating MET of sample paths. Greenhouse gas concentration contributes to a faster escape from the strong THC state (which is analogous to the current day setting). In the deterministic case, on-state has the widest basin of attraction as compared to both symmetric well potential inducing and strong fresh-water influxes. Our results for strong fresh-water input indicate extreme events of moderate jumps and frequencies facilitate the recovery of the strong THC, i.e., the 270 transition from off to on-state is boosted. Strong THC state attraction basin suffers its largest shrink to extreme events with a moderate jump and probability Lévy noise and noise-driven pathways linger in the attraction basin of the weak THC state. Both strong and weak THC state have equal basin stability and transitions between the states are equiprobable for the symmetric potential well inducing fresh-water forcing. The geometry of the probability density function also shows that pathways shuttle between both states. Extreme events characterized by medium jumps and frequency destabilize the states most, while events 275 with smaller jumps and high frequency bring less transient impact.
Extreme climate events can have severe aftermaths on the ocean. Warming of the globe due to the greenhouse gas accumulation may cause magnificent melting of polar ice caps and glaciers followed by fresh-water pulses of considerable volume followed by a collapse of the THC. Total shut down of the THC is a high risk rare event with catastrophic effects on human