Articles | Volume 29, issue 2
Research article
 | Highlight paper
11 May 2022
Research article | Highlight paper |  | 11 May 2022

Lévy noise versus Gaussian-noise-induced transitions in the Ghil–Sellers energy balance model

Valerio Lucarini, Larissa Serdukova, and Georgios Margazoglou

We study the impact of applying stochastic forcing to the Ghil–Sellers energy balance climate model in the form of a fluctuating solar irradiance. Through numerical simulations, we explore the noise-induced transitions between the competing warm and snowball climate states. We consider multiplicative stochastic forcing driven by Gaussian and α-stable Lévy – α(0,2) – noise laws, examine the statistics of transition times, and estimate the most probable transition paths. While the Gaussian noise case – used here as a reference – has been carefully studied in a plethora of investigations on metastable systems, much less is known about the Lévy case, both in terms of mathematical theory and heuristics, especially in the case of high- and infinite-dimensional systems. In the weak noise limit, the expected residence time in each metastable state scales in a fundamentally different way in the Gaussian vs. Lévy noise case with respect to the intensity of the noise. In the former case, the classical Kramers-like exponential law is recovered. In the latter case, power laws are found, with the exponent equal to α, in apparent agreement with rigorous results obtained for additive noise in a related – yet different – reaction–diffusion equation and in simpler models. This can be better understood by treating the Lévy noise as a compound Poisson process. The transition paths are studied in a projection of the state space, and remarkable differences are observed between the two different types of noise. The snowball-to-warm and the warm-to-snowball most probable transition paths cross at the single unstable edge state on the basin boundary. In the case of Lévy noise, the most probable transition paths in the two directions are wholly separated, as transitions apparently take place via the closest basin boundary region to the outgoing attractor. This property can be better elucidated by considering singular perturbations to the solar irradiance.

1 Introduction

1.1 Multistability of the Earth's climate

The climate system comprises the following five interacting subdomains: the atmosphere, the hydrosphere (water in liquid form), the upper layer of the lithosphere, the cryosphere (water in solid form), and the biosphere (ecosystems and living organisms). The climate is driven by the inhomogeneous absorption of incoming solar radiation, which sets up nonequilibrium conditions. The system reaches an approximate steady state, where macroscopic fluxes of energy, momentum, and mass are present throughout its domain, and entropy is continuously generated and expelled into the outer space. The climate features variability on a vast range of spatial and temporal scales as a result of the interplay of forcing, dissipation, feedbacks, mixing, transport, chemical reactions, phase changes, and exchange processes between the subdomains (see Peixoto and Oort1992, Lucarini et al.2014a, Ghil2015, and Ghil and Lucarini2020).

In the late 1960s Budyko (1969) and Sellers (1969) independently proposed that, in the current astronomical and astrophysical configuration, the Earth could support two distinct climates, namely the present-day warm (W) state and a competing one characterised by global glaciation, usually referred to as the snowball (SB) state. Their analysis was performed using one-dimensional energy balance models (EBMs), which, despite their simplicity, were able to capture the essential physical mechanisms in action, i.e. the interplay between two key feedbacks. The Boltzmann feedback is associated with the fact that warmer bodies emit more radiation, and it is a negative, stabilising one. Instead, the instability of the system is due to the presence of the so-called ice–albedo feedback, whereby an increase in the ice-covered fraction of the surface leads to a further temperature reduction for the planet because ice efficiently reflects the incoming solar radiation. These mechanisms are active at all spatial scales, including the planetary one (see Budyko1969, and Sellers1969). Such pioneering investigations of the multistability of the Earth's climate were later extended by Ghil (1976) (see also the later analysis by Ghil and Childress1987), who provided a comprehensive mathematical framework for the problem, based on the study of the bifurcations of the system. The main control parameter defining the stability properties is the solar irradiance S. Below the critical value SW→SB, only the SB state is permitted, whereas, above the critical value SSB→W, only the W state is permitted. Such critical values, which determine the region of bistability, are defined by bifurcations that emerge when, roughly speaking, the strength of the positive, destabilising feedbacks becomes as strong as the negative, stabilising feedbacks. Many variants of the models proposed by Budyko and Sellers have been discussed in the literature, all featuring by and large rather similar qualitative and quantitative features (Ghil1981; North et al.1981; North1990; North and Stevens2006). Furthermore, these models have long been receiving a great deal of attention from the mathematical community regarding the possibility of proving the existence of solutions and evaluating their multiplicity (Hetzer1990; Díaz et al.1997; Kaper and Engler2013; Bensid and Díaz2019).

Only later were these predictions confirmed by actual data. Indeed, geological and paleomagnetic evidence suggests that, during the Neoproterozoic era, between 630×106 and 715×106 years ago, the Earth went, at least twice, into major long-lasting global glaciations that can be associated with the SB state (see Pierrehumbert et al.2011, and Hoffman et al.1998). Multicellular life emerged in our planet shortly after the final deglaciation from the last SB state (Gould1989). The robustness and importance of the competition between the Boltzmann feedback and the ice–albedo feedback in defining the global stability properties of the climate has been confirmed by investigations performed using higher complexity models (Lucarini et al.2010; Pierrehumbert et al.2011), including fully coupled climate models (Voigt and Marotzke2010). While the mechanisms described above are pretty robust, the concentration of greenhouse gases and the boundary conditions defined by the extent and position of the continents have an impact on the values of SW→SB and SSB→W, as well as on the properties of the competing states. The presence of multistability has a key importance in terms of determining habitability conditions for Earth-like exoplanets (see Lucarini et al.2013, and Linsenmeier et al.2015).

Additionally, several results indicate that the phase space of the climate system might well be more complex than the scenario of bistability described above. Various studies (Lewis et al.2007; Abbot et al.2011; Lucarini and Bódai2017; Margazoglou et al.2021) performed with highly nontrivial climate models report the possible existence of additional competing states, up to a total of five (Brunetti et al.2019; Ragon et al.2022). In Margazoglou et al. (2021), it is argued that, in fact, one can see the climate as a multistable system where multistability is realised at different hierarchical levels. As an example, the tipping points (Lenton et al.2008; Steffen et al.2018) that characterise the current (W) climate state can be seen as a manifestation of a hierarchically lower multistability with respect to the one defining the dichotomy between the W and SB states.

1.2 Transitions between competing metastable states: Gaussian vs. Lévy noise

Clearly, in the case of autonomous systems where the phase space is partitioned in more than one basin of attraction of the corresponding attractors and the basin boundaries, the asymptotic state of the system is determined by its initial conditions. Things change dramatically when one includes time-dependent forcing which allows for transitions between competing metastable states (Ashwin et al.2012). In particular, following the viewpoint originally proposed by Hasselmann (1976), whereby the fast variables of the climate system act as stochastic forcings for the slow ones (Imkeller and von Storch2001), the relevance of studying noise-induced transitions between competing states has become apparent (Hänggi1986; Freidlin and Wentzell1984). This viewpoint, where the noise is usually assumed to be Gaussian distributed, has provided very valuable insight on the multiscale nature of climatic time series (Saltzman2001) and is related to the discovery of phenomena like stochastic resonance (Benzi et al.1981; Nicolis1982).

Metastability is ubiquitous in nature, and advancing its understanding is a key challenge in complex system science at large (Feudel et al.2018). In general, the transitions between competing metastable states in stochastically perturbed multistable systems take place, in the weak noise limit, through special regions of the basin boundaries, which are named edge states. The edge states are saddles, and the trajectories initialised in the basin boundaries are attracted to them, but there is an extra direction of instability, so that a small perturbation sends an orbit towards one of the competing metastable states with a probability of one (Grebogi et al.1983; Ott2002; Kraut and Feudel2002; Skufca et al.2006; Vollmer et al.2009). In the case the edge state supports chaotic dynamics, we refer to it as the melancholia (M) state (Lucarini and Bódai2017). In previous papers, we have shown that it is possible to construct M states in high-dimensional climate models (Lucarini and Bódai2017) and to prove that the nonequilibrium quasi-potential formalism introduced by Graham (1987) and Graham et al. (1991) provides a powerful framework for explaining the population of each metastable state and the statistics of the noise-induced transitions. In the weak noise limit, edge states act as gateways for noise-induced transitions between the metastable states (Lucarini and Bódai2019, 2020; Margazoglou et al.2021); see also a recent study on a nontrivial metastable prey–predator model (Garain and Sarathi Mandal2022). The local minima and the saddles of the quasi-potential Φ, which generalises the classical energy landscape for non-gradient systems, correspond to competing metastable states and to edge states, respectively. In our investigation, the climate system is forced by adding a random – Gaussian-distributed – component to the solar irradiance, which impacts, in the form of multiplicative noise, only a small subset of the degrees of freedom of the system. We remark that such a choice of the stochastic forcing does not fully reflect the physical realism, as the variability of the solar irradiance has a more complex behaviour (Solanki et al.2013). Instead, noise acts as a tool for exploring the global stability properties of the system, and injecting noise as fluctuation of the solar irradiance has the merit of impacting the Lorenz energy cycle, thus effecting all degrees of freedom of the system (Lucarini and Bódai2020). See also the recent detailed mathematical analysis of the stochastically perturbed one-dimensional EBMs presented in Díaz and Díaz (2021).

A major limitation of this mathematical framework is the need to rigidly consider Gaussian noise laws, even if considerable freedom is left as to the choice of the spatial correlation properties of the noise. It seems natural to attempt a generalisation by considering the whole class of α-stable Lévy noise laws. Lévy processes (Applebaum2009; Duan2015), described in detail in Appendix A, are fundamentally characterised by the stability parameter α(0,2], where the α=2 case corresponds to the Gaussian case (which is, indeed, a special Lévy process). In what follows, when we discuss Lévy noise laws, we refer to α(0,2).

Note that α-stable Lévy processes have played an important role in geophysics as they have provided the starting point for defining the multiplicative cascades also referred to as universal multifractal. This framework has been proposed as way to analyse and simulate, at climate scales, the ubiquitous intermittency and heavy-tailed statistics of clouds (Schertzer and Lovejoy1988), rain reflectivity (Tessier et al.1993; Schertzer and Lovejoy1997), atmospheric turbulence (Schmitt et al.1996), and soil moisture (Millàn et al.2016). On longer timescales, multiplicative cascades have been used to interpret temperature records in the summit ice core Schmitt et al. (1995, see Lovejoy and Schertzer2013, for a summary of this viewpoint). Mathematicians, on the other hand, have defined a Lévy multiplicative chaos (Fan1997; Rhodes et al.2014) as a more mathematically tractable alternative to the universal multifractal. Finally, we remark that fractional Fokker–Planck equations have been proposed by Schertzer et al. (2001) to investigate the properties of nonlinear Langevin-type equations forced by an α-stable Lévy noise with the goal of analysing and simulating anomalous diffusion.

Following Ditlevsen (1999), it has become apparent that more general classes of α-stable Lévy noise laws might be useful for modelling noise-induced transitions in the climate system like Dansgaard–Oeschger events, which are sequences of periods of abrupt warming followed by slower cooling that occurred during the last glacial period (Barker et al.2011). The viewpoint by Ditlevsen (1999) was particularly effective in stimulating mathematical investigations into noise-induced escapes from attractors, where, as stochastic forcing, one chooses a Lévy, rather than Gaussian, noise (Imkeller and Pavlyukevich2006a, b; Chechkin et al.2007; Debussche et al.2013). Such analyses have clarified that a fundamental dichotomy exists with the classical Freidlin and Wentzell scenario mentioned above, even if phenomena like stochastic resonance can also be recovered in this case (Dybiec and Gudowska-Nowak2009; Kuhwald and Pavlyukevich2016). Whereas, in the Gaussian case, transitions between competing attractors occur as a result of the very unlikely combination of many steps all going in the right direction, in the Lévy case, transitions result from individual, very large and very rare jumps. Recently, Duan and collaborators have made fundamental progress in achieving a variational formulation of the Lévy noise-perturbed dynamical systems (Hu and Duan2020) and in developing corresponding methods for data assimilation (Gao et al.2016) and data analysis (Lu and Duan2020). In terms of applications, Lévy noise is becoming a more and more a popular concept and tool for studying and interpreting complex systems (Grigoriu and Samorodnitsky2003; Penland and Sardeshmukh2012; Zheng et al.2016; Wu et al.2017; Serdukova et al.2017; Cai et al.2017; Singla and Parthasarathy2020; Gottwald2021).

The contribution by Gottwald (2021) is especially worth recapitulating because of its methodological clarity. There, the idea is, following Ditlevsen (1999), to provide a conceptual deterministic climate model able to generate a Lévy-noise-like signal to describe, at least qualitatively, abrupt climate changes similar to Dansgaard–Oeschger events. A key building block is the idea proposed in Thompson et al. (2017) that a Lévy noise can be produced by integrating the so-called correlated additive and multiplicative (CAM) noise processes, which are defined by starting from standard Gaussian processes. The other key ingredient is to consider the atmosphere as the fast component in the multiscale model and deduce, using homogenisation theory (Pavliotis and Stuart2008; Gottwald and Melbourne2013), that its influence on the slower climate components can be closely represented as a Gaussian forcing. Finally, the temperature signal is cast as the integral of a CAM process.

We remark that Gaussian and Lévy noise can be associated with stochastic forcings of a fundamentally different nature. One might think of Gaussian noise as being associated to the impact of very rapid unresolved scales of motion on the resolved ones Pavliotis and Stuart (2008). Instead, one might interpret α-stable Lévy noise as describing, succinctly, the impact of what in the insurance sector are called acts of God (e.g. an asteroid hitting the Earth, a massive volcanic eruption, or the sudden collapse of the West Antarctic ice sheet).

1.3 Outline of the paper and main results

We consider here the Ghil–Sellers Earth's EBM (Ghil1976), a diffusive, one-dimensional energy balance system, governed by a nonlinear reaction–diffusion parabolic partial differential equation. We stochastically perturb the system by adding random fluctuations to the solar irradiance; therefore, the noise is introduced in a multiplicative form. We study the transitions between the two competing metastable climate states and carry out a comparison of the effect of considering Lévy vs. Gaussian noise laws of weak intensity ε.

The main challenges of the problem are (a) the fact that we are considering dynamical processes occurring in infinite dimensions (Doering1987; Duan and Wang2014; Alharbi2021) and (b) the consideration of multiplicative Lévy noise laws (Peszat and Zabczyk2007; Debussche et al.2013). We characterise noise-induced transitions between the competing climate basins and quantify the effect of noise parameters on them by estimating the statistics of escape times and empirically constructing mean transition pathways called instantons.

The results obtained confirm that, in the weak noise limit ε→0, the mean residence time in each metastable state driven by Gaussian vs. Lévy noise has a fundamentally different dependence on ε. Indeed, as expected, in the Gaussian case, the residence time grows exponentially with ε−2, which is, thus, in basic agreement with the well-known Kramers (1940) law and the previous studies performed on climate models (Lucarini and Bódai2019, 2020). Instead, in the case of α-stable noise laws, the residence time increases with εα. We perform simulations for α={0.5,1.0,1.5}. The obtained scaling can be explained by effectively treating the Lévy noise as a compound Poisson process, and this is in agreement with what is found for low-dimensional dynamics (Imkeller and Pavlyukevich2006a, b) and for the infinite dimensional stochastic Chafee–Infante reaction–diffusion equation (Debussche et al.2013) in the case of additive noise. This might indicate that such scaling laws are more general than what has been so far assumed.

Furthermore, we find clear confirmation that, in the case of Gaussian noise in the weak noise limit, the escape from either attractor's basin takes place through the edge state. Indeed, the most probable paths for both thawing and freezing processes meet at the edge state and have distinct instantonic and relaxation sections. In turn, for Lévy noise in the weak noise limit, the escapes from a given basin of attraction occur through the boundary region closest to the outgoing attractor. Hence, the paths are very different from the Gaussian case (especially so for the freezing transition) and, somewhat surprisingly, are identical regardless of the value of α considered. These properties can be better understood by studying the impact of including singular perturbations to the value of the solar irradiance.

The rest of the paper is organised as follows. In Sect. 2, we present the Ghil–Sellers EBM and summarise its most important dynamical aspects and the steady-state solutions and their stability. The stochastic partial differential equation obtained by randomly perturbing the solar irradiance in the EBM is given in Sect. 3, where we also clarify the mathematical meaning of the solution of the stochastic partial differential equation. Section 3 also introduces the mean residence time and most probable transition path between the competing climate states. The numerical methods are also briefly presented. In Sect. 4, we discuss our main results. In Sect. 5, we present our conclusions and perspectives for future investigations. Finally, Appendix A presents a succinct description of α-stable Lévy processes, Appendix B sketches the derivation of the scaling laws for mean residence times presented in Debussche et al. (2013), Appendix C explores the behaviour and dynamics of singular Lévy perturbations of different duration, and Appendix D presents a tabular summary of the statistics of the problem.

2 The Ghil–Sellers energy balance climate model

The Ghil–Sellers EBM (Ghil1976) is described by a one-dimensional nonlinear, parabolic, reaction–diffusion partial differential equation (PDE) involving the surface temperature T field and the transformed space variable x=2ϕ/π[-1,1], where ϕ[-π/2,π/2] is the latitude. The model describes the processes of energy input, output, and diffusion across the domain and can be written as follows:

(1) C ( x ) T t = D I ( x , T , T x , T x x ) + D II ( x , T ) - D III ( T ) ,

where C(x) is the effective heat capacity, and T=T(x,t) has boundary and initial conditions, as follows:

(2) T x ( - 1 , t ) = T x ( 1 , t ) = 0 , T ( x , 0 ) = T 0 ( x ) .

The equation does not depend explicitly on the time t. The subscripts t and x refer to partial differentiation. The first term – DI – on the right-hand side of Eq. (1) can be written as follows:

(3) D I ( x , T , T x , T x x ) = 4 π 2 cos ( π x / 2 ) [ cos ( π x / 2 ) K ( x , T ) T x ] x ,

and describes the convergence of meridional heat transport performed by the geophysical fluids. The function K(x,T) is a combined diffusion coefficient, expressed as follows:


The empirical functions k1(x) and k2(x) are eddy diffusivities for sensible and latent heat, respectively, and g(T) is associated with the Clausius–Clapeyron relation, which describes the relationship between temperature and saturation water vapour content of the atmosphere.

The second term – DII – on the right-hand side of Eq. (1) describes the energy input associated with the absorption of incoming solar radiation and can be written as follows:

(6) D II ( x , T ) = μ Q ( x ) 1 - α a ( x , T ) ,

where 𝒬(x) is the incoming solar radiation, and αa(x,T) is the surface reflectivity (albedo), which is expressed as follows:

(7) α a ( x , T ) = b ( x ) - c 1 T m + min T - c 2 z ( x ) - T m , 0 c ,

where the subscript {}c denotes a cutoff for a generic quantity h, defined as follows:

(8) h c = h min h h min , h h min < h < h max , h max h max h .

The term c2z(x) in Eq. (7) indicates the difference between the sea level and surface level temperatures, and b(x) is a temperature independent empirical function of the albedo. The parameterisation given in Eqs. (7)–(8) encodes the positive ice–albedo feedback. The relative intensity of the solar radiation in the model can be controlled by the parameter μ.

The last term – DIII – on the right-hand side of Eq. (1) describes the energy loss to space by outgoing thermal planetary radiation and is responsible for the negative Boltzmann feedback. It can be written as follows:

(9) D III ( T ) = σ T 4 1 - m tanh ( c 3 T 6 ) ,

where σ is the Stefan–Boltzmann constant, and the emissivity coefficient is expressed as 1−mtanh (c3T6). Such a term describes, in a simple yet effective way, the greenhouse effect by reducing infrared radiation losses. The values of the empirical functions C(x),Q(x),b(x),z(x),k1(x),k2(x) at discrete latitudes and empirical parameters c1,c2,c3,c4,c5,σ,m,Tm are taken from Ghil (1976), as modified in Bódai et al. (2015). The choice of the empirical functions and parameters are extensively discussed in Ghil (1976). Of course, one might reasonably wonder about the robustness of our modelling strategy. Indeed, a plethora of EBMs analogous to the one described here have been presented in the literature, where slightly different parameterisations for the diffusion operator, for the albedo, and for the greenhouse effect are introduced. Such models are in fundamental agreement, both in terms of physical (Ghil1981; North et al.1981; North1990; North and Stevens2006) and mathematical properties (Hetzer1990; Díaz et al.1997; Kaper and Engler2013; Bensid and Díaz2019).

Figure 1(a) Stationary solutions TW(x), TSB(x), and TM(x) in Kelvins (K) of the zonally averaged energy balance model Eq. (1). (b) Bifurcation diagram of the average temperature T as a function of control parameter μ.


In this study, we consider μ=1.05. For this value of μ, two stable asymptotic states – the W and the SB states – co-exist (see Fig. 1b). Indeed, a codimension of one manifold separates the basins of attraction of the W and SB states. We refer to DW (DSB) as the basin of attraction of the W (SB state). We refer to B as the basin boundary, which includes a single edge state M. Therefore, the system has three stationary solutions TW(x), TSB(x), and TM(x) for the W, SB, and M state, respectively, as shown in Fig. 1a. In Ghil (1976), the three stationary solutions were obtained by equating Tt to 0, and it was shown, through linear stability analysis, that the stationary solutions TW and TSB are stable, while TM is unstable. In Bódai et al. (2015) the unstable solution TM was constructed using a modified version of the edge-tracking algorithm (Skufca et al.2006).

Following previous studies (Bódai et al.2015; Lucarini and Bódai2019, 2020; Margazoglou et al.2021), when visualising our results, we apply a coarse graining to the phase space of the model. In what follows, we perform a projection on the plane spanned by the spatially averaged temperature T and the averaged Equator minus the poles' temperature difference ΔT, which is defined as follows:

(10)T=[T(x,t)]01,(11)ΔT=[T(x,t)]01/3-[T(x,t)]1/31, where(12)[T(x,t)]xlxh=xlxhcos(πx/2)T(x,t)dxxlxhcos(πx/2)dx.

Such a representation allows for a minimal yet still physically relevant description of the system. Indeed, changes in the energy budget of the system (warming versus cooling) are, to a first approximation, related to variations in T, while the large-scale energy transport performed by the geophysical fluids is controlled by ΔT. The boundary between high and low latitude in Eq. (11) is established at x=±1/3, i.e. at 30 N/S. Additionally, in some visualisations, we consider, as a third coordinate, the fraction of the surface with a below-freezing temperature (therefore, we expect 1 for global glaciation and 0 for no ice). We refer to this variable as I, and it is an attempt to extract an observable that resembles the sea ice percentage of the Earth's surface. Thus, the stationary solutions TW(x), TSB(x), and TM(x), in terms of ΔT and T, correspond to ΔTW= 16 K, ΔTSB= 8.3 K, ΔTM= 17.5 K, TW= 297.7 K, TSB= 235.1 K, TM= 258 K, IW=0.2, ISB=1, and IM=1.

3 Background and methodology

3.1 Stochastic energy balance model

In order to analyse the influence of random perturbations on the deterministic dynamics of the climate model described in Sect. 2, we perturb the relative intensity μ of the solar irradiance by including a symmetric α-stable Lévy process and rewrite Eq. (1) in the form of the following stochastic partial differential equation (SPDE):

(13) C ( x ) T t = D I x , T , T x , T x x + D II ( x , T ) 1 + ε / μ L ˙ α ( t ) - D III ( T ) ,

where the boundary and initial conditions defined by Eq. (2) apply to the stochastic temperature field 𝒯. Here the parameter ε>0 controls the noise intensity, and (Lα(t)t≥0) is a symmetric α-stable process defined in Appendix A. We consider symmetric processes because we want to have a simple mathematical model allowing for transitions in both the SB→W and the W direction →SB direction. Instead, a strongly skewed process would have made it very hard to explore the full phase space because a lack of symmetry would invariably favour one of the two transitions. As mentioned before, we refer to the Lévy case if the stability parameter α(0,2), so that we consider a jump process. We recall that the jumps become more frequent and less intense as α increases.

We define L˙(t)=Q(x)[1-αa(x,T)]L˙α(t), as the generalised derivative of a stochastic process in a suitably defined functional space. Equation (13) features multiplicative noise. The research interest on this type of SPDE (Doering1987; Peszat and Zabczyk2007; Duan and Wang2014; Alharbi2021) is mainly focused on defining weak, strong, mild, and martingale solutions, in specifying under which conditions these solutions exist and are unique, and in constructing numerical approximation schemes for the solutions (Davie and Gaines2000; Cialenco et al.2012; Burrage and Lythe2014; Jentzen and Kloeden2009; Kloeden and Shott2001), among other aspects.

First, let us define the concept of a mild solution in this context. Let (Ω,F,P) be a given complete probability space and H(,,) a separable Hilbert space with a norm and inner product ,. Equation (13) can be rewritten in the more general form, as follows:

(14) T t = A ( x ) [ E ( x , T ) T x ] x + F ( x , T ) + ε G ( x , T ) L ˙ α ( t ) , T x ( - 1 , t ) = T x ( 1 , t ) = 0 , T ( x , 0 ) = T 0 ( x ) ,

where A,E,F,G are Lipschitz functions defined on [-1,1]×H and G(x,T)L˙α(t)=L˙(t). Under certain assumptions (Yagi2010), the problem (Eq. 14) is formulated as a Cauchy problem, for which a local mild solution, a progressively measurable process 𝒯(t), for all t[0,tF] and T0H has the following integral representation:

(15) T ( t ) = Ψ ( t ) T 0 + 0 t Ψ ( t - s ) Υ ( T ( s ) ) d s + ε 0 t Ψ ( t - s ) G ( T ( s ) ) d β + ε 0 t Ψ ( t - s ) G ( T ( s ) ) d γ ,

where the dependence on x is kept implicit, and β (γ) is the Poisson random measure (compensated Poisson random measure) defined through Lévy–Itô decomposition. Instead, Ψ(t) with t0 is the evolution operator having the generalised semigroup property for the family of sector operators with the bounded inverses, and Υ(T)=T+F(x,T), 𝒯∈H is a nonlinear operator, which we assume to be Lipschitz continuous. Following the abstract theory presented in Yagi (2010), under certain structural assumptions for the operators Ψ and Υ and for the functional space, one can prove that the solution Eq. (15) is the unique local mild solution of Eq. (14).

As mentioned above, things are radically different for the special case α=2, which corresponds to Gaussian noise. In this case, we revisit Eq. (14), and we define L˙α=2(t)=W˙(t), where (W(t)t≥0) is a Wiener process. We then define W˙(t)=G(x,T)W˙(t) as the generalised derivative of a Wiener process in a suitably defined functional space.

3.2 Noise-induced transitions: mean escape times

By incorporating stochastic forcing into the system, its long-time dynamics change significantly, allowing transitions between the competing basins. This dynamical behaviour is called metastability and is graphically captured by Fig. 2, where, in Fig. 2a and b, a typical spatiotemporal evolution of the temperature field is shown for the stability parameters α=0.5 and α=1.5, respectively. Instead, in Fig. 2c and d, the temporal evolution of the global temperature T and of the averaged Equator and poles' temperature difference Δ𝒯 (as defined in Eqs. 1011) is correspondingly shown. In what follows, we investigate the time statistics and the paths of the transitions between such basins.

In a complete probability space (Ω,F,P) we define the first exit time τx of a cádlág mild solution T(;x) of Eq. (13), starting at the xDW/SB domain of a W/SB climate stable state as follows:

(16) τ x ( ω ) = inf t > 0 | T t ( ω , x ) D W / SB , ω Ω , x H .

The mean escape time is then expressed by 𝔼[τx(ω)]. In the case of the infinite dimensional multistable reaction–diffusion system described by Chafee–Infante equation under the influence of additive infinite-dimensional α-stable Lévy noise – α(0,2) – it was shown (Debussche et al.2013) that, in the weak noise limit, ε→0 the mean escape time from one of the competing basins of attraction increases as εα. In such a limit, the jump diffusion system reduces the Markov chain to a finite state, with values in the set of stable states. Details of this method are given in Appendix B. Similar results have been obtained for bistable one-dimensional stochastic differential equations (SDEs; Imkeller and Pavlyukevich2006a, b). The basic reason behind this result is that, in order to study the transitions between the competing basins of attraction, one can treat the Lévy noise as a compound Poisson process, where jumps arrive randomly, according to a Poisson process, and the size of the jumps x is given by a stochastic process that obeys a specified probability distribution. For a symmetric α-stable Lévy process, such a distribution asymptotically decreases as |x|-1-α, as discussed in Appendix A. Let us assume that positive values of x bring the state of the system closer to the basin boundary (as in the case of positive fluctuations of the solar irradiance when studying escapes from the SB state). Assuming a simple geometry for the basin boundary, we find that a transition takes place when an event larger than a critical value xcrit>0 is realised. The probability of such an event scales with xcrit-α. A similar argument applies when considering transitions triggered by negative fluctuations of the stochastic variable. Small-size events, which occur frequently and correspond to the non-occurrence of jumps, do not actually play any relevant role in determining the transitions, while they are responsible for the variability within each basin of attraction.

We now consider the case α=2. While the corresponding finite dimensional problem is thoroughly documented in the literature (Freidlin and Wentzell1984) and has been applied in a similar context by some of the authors (Lucarini and Bódai2019, 2020; Ghil and Lucarini2020; Margazoglou et al.2021), the treatment of infinite dimensional SDEs driven by an infinite dimensional Wiener process via the Freidlin–Wentzell theory requires further extension. In the present context, we refer to Budhiraja and Dupuis (2000) and Budhiraja et al. (2008) and references therein, where the problem of an infinite dimensional reaction–diffusion equation driven by an infinite dimensional Wiener process has been addressed.

We assume that steady-state conditions and ergodicity are met, and we also assume that the analysing system is bistable and a unique edge state is present at the basin boundary, as in the case studied here. In the case of Gaussian noise, transitions between the competing basins of attraction are not determined by a single event as in the 0<α<2 case but, instead, occur as a result of very unlikely combinations of subsequent realisations of the stochastic variable acting as a forcing. In the weak noise limit, the transitions occur according to the least unlikely (yet very unlikely) chain of events, whose probability is described using a large deviation law (Varadhan et al.1985). One finds that the mean escape time from either basin of attraction decreases exponentially with increasing noise intensity ε and is given by a generalised Kramers' law, as follows:

(17) E [ τ W / SB ( ε ) ] exp 2 Δ Φ W M / SB M ( T ) ε 2 ,

where ΔΦWM=ΦM(T)-ΦW(T) is the height of the quasi-potential barrier in the W attractor; correspondingly, ΔΦSBM(T)=ΦM(T)-ΦSB(T) is the height of the quasi-potential barrier in the SB attractor, and Φ is the Graham quasi-potential mentioned above (Graham1987; Graham et al.1991).

3.3 Noise-induced transitions: most probable transition paths

In the weak noise limit, the most probable path to escape an attractor is defined by a class of trajectories named instantons (Grafke et al.2015, 2017; Bouchet et al.2016; Grafke and Vanden-Eijnden2019) or maximum likelihood escape paths (Lu and Duan2020; Dai et al.2020; Hu and Duan2020; Zheng et al.2020). However, note that different noise laws result into possibly radically different instantonic trajectories (Dai et al.2020; Zheng et al.2020).

In our case, the theory indicates that, if the stochastic forcing is Gaussian, under a rather general hypothesis, the instanton will connect the attractor W/SB with the edge state M, which then acts as gateway for noise-induced transitions. Once the quasi-potential barrier is overcome, a free-fall relaxation trajectory links M with the competing attractor SB/W. For equilibrium systems, (e.g. for gradient flows), where a detailed balance is achieved, the relaxation and instantonic trajectories within the same basin of attraction are identical. On the contrary, for non-equilibrium systems, the relaxation and instantonic trajectories will differ and will only meet at the attractor. (See a detailed discussion of this aspect and of the dynamical interpretation of the quasi-potential Φ in Lucarini and Bódai2020, and Margazoglou et al.2021). Instead, if the noise is of Lévy type, the theory formulated for simpler equations suggests that the instanton will connect the attractor with a region on the basin boundary that is the nearest, in the phase space, to the attractor, as the concept of the quasi-potential is immaterial (Imkeller and Pavlyukevich2006a, b).

In general, the maximum likelihood transition trajectory 𝒯M(t) can be defined (Zheng et al.2020; Lu and Duan2020) as a set of system states at each time moment t[0,tf] that maximises the conditional probability density function p(.|.;.) of the passage from the origin stable state ϕW/SB to the destination stable state ϕSB/W and is expressed as follows:

(18) T M ( t ) = arg max x p ( T ( t ) = x | T ( 0 ) = x 0 ; T ( t f ) = x f ) = p T ( t f ) = x f | T ( t ) = x p T ( t ) = x | T ( 0 ) = x 0 p ( T ( t f ) = x f | T ( 0 ) = x 0 ) ,

where x0 (xf) belongs to the basin of attraction DW/SB (DSB/W), and p(.|.) is the probability density function evolving according to the Fokker–Planck equation (Risken1996). This method is applicable either if efficient numerical algorithms are available to solve the Fokker–Planck equation associated to the studied stochastically driven system, or, empirically, when considering a large ensemble of simulations. Note that this is not an asymptotic approach, i.e. it does not require a weak noise limit ε→0 and is applicable for systems with either Gaussian or non-Gaussian noise. Yet, in the weak noise limit, the definition (Eq. 18) leads to constructing the optimal transition paths described above.

In the following section, for practical purposes, we construct such optimal transition path in the coarse-grained 2D phase space (T,ΔT) and 3D phase space (T,ΔT,I) of the variables defined in Sect. 2 by averaging the ensemble of transitions connecting the two competing states in the weak noise limit.

Figure 2The metastable behaviour of the solution path of a stochastic energy balance model (Eq. 13), for ε=0.04, T0= 300 K, and t(0,300) years, for (a) α=0.5, (b) α=1.5, and the (c) temperature average T and (d) temperature contrast at low and high latitudes. Δ𝒯, for ε=0.01, α=0.5,andα=1.5. Red, green, or blue dash-dotted lines portray the stationary climate states of TW/TM/TSB, respectively.


3.4 Numerical methods

We solve Eq. (13) through the MATLAB pdepe function, which is well suited for solving 1D parabolic and elliptic PDEs. We discretise the 1D space with a regular grid of 201 grid points, following Bódai et al. (2015).

The time span of integration t[0,Tf], varies for different cases, with Tf(105,15×105) years, with a time stepping of 1 year. Each year, we consider a different value for the relative solar irradiance by extracting a random number Zj (see Eq. 19). To simulate the stochastic noise term εLα(t), which is added in the parameter μ in Eq. (13), we use the recursive algorithm from Duan (2015). The process values Lα(t1),,Lα(tN) at each moment tj, j∈ℕ are obtained via the following:

(19) L α ( t j ) = L α ( t j - 1 ) + ( t j - t j - 1 ) 1 α Z j , j = 1 , , N ,

where the second term is an independent increment, and Zj are the independent standard symmetric α-stable random numbers generated by an algorithm in Weron and Weron (1995, see also Grafke et al.2015, for a detailed explanation of the steps above). For illustrative reasons, some sample solutions of Eq. (13) for different values of α are shown in Fig. 2a and b.

For the numerical simulations discussed below, we consider α=(0.5,1.0,1.5,2) and ε(0.0001,0.3). We select ε in such a way that the noise intensity is strong enough to induce at least an order of 10 transition, given our constraints in the time length of the simulations, and weak enough that we are not far from the weak noise limit, where the scaling laws discussed above apply and transitions paths are well organised. Our simulations are performed by taking the Itô interpretation for the stochastic equations.

We remark that, when we consider Lévy noise, it does happen that, for some years, the solar irradiance has negative values. Of course these conditions bear no physical relevance, and are a necessary result of considering unbounded noise. Nonetheless, we have allowed for this to occur in our simulations in order to be able to stick to the desired mathematical framework. We remind the reader that this study does not aim at capturing, with any high degree of realism, the description of the actual evolution of climate. At any rate, as can be understood from the discussion below in Sect. 4.2.2 and from what is reported in Appendix C, were we to consider longer-lasting (e.g. 2 years vs. 1 year) fluctuations of the solar irradiance, a satisfactory exploration of the transitions between the competing W and SB states would be possible with a greatly reduced occurrence of such unphysical events, and the basic reason for this being the presence of a larger factor (tj-tj-1)1α in Eq. (19).

4 Results and discussion

In what follows, we aim at addressing three main questions: (1) what are the temporal statistics of the SB→W and W→SB transitions? (2) What are the typical transition pathways? (3) What are the fundamental differences between transitions caused by Gaussian vs. Lévy noise? A summary of the results of the numerical simulations is given in Table D1 in Appendix D, including the sample size, i.e. number of transitions, point estimates for mean escape time, and their 0.95 confidence intervals for exits from both the W and SB basins. See the data availability section for information on how to access the supplement (Lucarini et al.2022), which contains the raw data produced in this study and some illustrative animations portraying noise-induced transitions between the two competing metastable states.

4.1 Mean escape time

Our analysis confirms that there is a fundamental dichotomy in the statistics of mean escape times between Lévy noise and Gaussian noise-induced transitions.

Figure 3a shows the dependence of the mean escape time from either attractor on ε and α for the Lévy case. The red circles (blue squares) correspond to escapes from the W (SB) basin (see Lucarini et al.2022, for additional details). The scaling ε-α presented in Eq. (B6) is shown by the dotted black line for each value of α. We also portray the best power law fit of the mean residence time with respect to ε for each value of α; the confidence intervals of the exponent are shown in Table 1. Our empirical results seem to indicate, at least in this case, an agreement with the εα scaling presented and discussed earlier in the paper. This points at the possibility that the εα scaling might apply in more general conditions than what has been, as of yet, rigorously proven, and this is specifically the case when multiplicative Lévy noise is considered. The stochastically perturbed trajectories forced by Lévy noise consist of jumps, and the probability of occurrence of a high jump, which can trigger the escape from the reference basin of attraction, is polynomially small in noise intensity ε.

The Gaussian case – where no jumps are present – is portrayed in Fig. 3b. We show, in semi-logarithmic scale, the mean residence time versus 1/ε2. We perform a successful linear fit of the logarithm of the mean residence time in either attractor versus 1/ε2, and using Eq. (17), we obtain an estimate of the local quasi-potential barrier ΔΦW/SBM, which is half of the slope of the corresponding straight lines of the linear fit (see the last column of Table 1). We conclude that, for μ=1.05, the local minimum of Φ corresponding to the W state is deeper than the one corresponding to the SB state.

Figure 3Estimates of the mean escape time 𝔼(τ) (in years) from the W (blue circles) and SB (red squares) states, as a function of the noise intensity ε. (a) Lévy noise for α=0.5, 1.0, and 1.5, with the dotted line being the corresponding prediction from Eq. (B6), while the straight green (SB) and dashed black (W) lines are the fittings of Eq. (B6) of the relevant dataset. (b) Gaussian noise, with straight green (SB) and dashed black (W) lines being the fit of Eq. (17), is shown.


Table 1Estimates of the exponent α via the fitting of Eq. (B6) for the Lévy case (three first columns) and of the energy barrier ΔΦW/SBM via the fitting of Eq. (17) for the Gaussian case (last column). In parenthesis, the estimated error of the last digit is shown.

Download Print Version | Download XLSX

4.2 Escape paths for the noise-induced transitions

We now explore the geometry of the transition paths associated with the metastable behaviour of the system. We first discuss the case of Gaussian noise because it is indeed more familiar and more extensively studied.

4.2.1 Gaussian noise

We estimate the transition paths by averaging among the escape plus relaxation trajectories using the run performed with the weakest noise (see Table D1). We first perform our analysis in the 2D-projected state space defined by (T,ΔT). We prescribe two small circular-shaped regions enclosing the two deterministic attractors and search the time series of the portions of the whole trajectory that leave one of these regions and reach the other one. This creates two subsets of our full dataset from which we build a 2D histogram for each of the SB→W and W→SB transitions in the projected space. We then estimate the most probable transition paths by finding, for each bin value of T, the peak of the histogram in the Δ𝒯 direction. The distributions are very peaked, and almost indistinguishable estimates for the instantonic and relaxation trajectories are obtained when computing the average of Δ𝒯 according to the 2D histogram conditional on the value of T.

Figure 4(a) Invariant measure in the 2D-projected state space defined by T,ΔT. The coloured points indicate the deterministic attractors of the SB (blue), M (green), and W (red) states, and the blue (red) line is the stochastically averaged transition paths for the W→SB (SB→W) transitions. Dashed (solid) lines are the relaxation (instantonic) trajectories. The arrows show the direction of transitions. The top left inset shows that the dark blue (red) contours portray the ensembles of the transition paths between W→SB (SB→W). Here, the system is driven by Gaussian noise with ε=0.14. (b) The invariant measure and most probable transition paths (W→SB in blue and SB→W in red) in the 3D-projected state space are defined by T,ΔT,I. The darker brown shading indicates a higher probability density for the corresponding isosurface. The 2D projections on the T,I, and (Δ𝒯,ℐ) planes are shown. The location of the M state is given by a pink square. Here, the system is driven by Gaussian noise, with ε=0.16.


In the background of Fig. 4a, we show the empirical estimate of the invariant measure in the 2D-projected state space defined by (T,ΔT). Additionally, we indicate the position of the deterministic attractors, where the blue (red) circle corresponds to the SB (W) state and of the M state (green square). In the inset of Fig. 4a, we present the ensemble of W→SB (SB→W) transitions as deep blue (red) contours. The most probable transition paths are shown in blue for the W→SB and in red for the SB→W. The instantonic portion of the blue (red) line is the one connecting the W (SB) attractor to the M state and is portrayed as a solid line, while the relaxation portion, connecting the M state with the SB (W) attractor, is portrayed as a dashed line. Within each basin of attraction, the instantonic and relaxation trajectories do not coincide, and, instead, only meet at the corresponding attractor and at the M state. This is particularly clear for the W state. The presence of such a loop, proving the existence of non-vanishing probability currents and the breakdown of detailed balance, is a signature of non-equilibrium dynamics, which was also observed in Margazoglou et al. (2021) and has, instead, gone undetected in Lucarini and Bódai (2019) and Lucarini and Bódai (2020). See Lucarini et al. (2022) for some illustrative simulations of the transitions.

Let us provide some physical interpretation of how the transitions occur. Looking at the SB→W most probable path, the escape includes a simultaneous increase in T and Δ𝒯. In practice, a SB→W transition takes place when, starting at the SB state, one has a (rare) sequence of positive anomalies in the fluctuating solar irradiance μ̃, i.e. μ̃>μ. While the planet is warming globally, the Equator is warming faster than the poles, resulting in a positive rate ΔT˙>0, because it receives, in relative and absolute terms, more incoming solar radiation. Considering that the Equator also in the SB state is warmer than the poles, the melting of the ice conducive to the transition occurs first at the Equator, with a subsequent decrease in the albedo in low latitudes. Once the system crosses the M state, and supposing that persistent μ̃<μ do not appear at this stage, the system will relax towards the W state. The relaxation includes a consistent global warming of the planet but with a change in sign in the rate of ΔT˙, and a subsequent decrease in Δ𝒯, implying that as soon as the temperature at the Equator has risen enough, the poles will then warm at a faster pace because the ice–albedo effect kicks in.

The global freezing of the planet associated with the W→SB transition is qualitatively similar but not identical to the reverse SB→W process. Notice a considerable overlap of the transition paths ensembles in both basins of attraction, shown as red and blue contours in the inset of Fig. 4a. This implies the presence of less extreme non-equilibrium conditions compared to what was observed in Margazoglou et al. (2021), where the W→SB and SB→W transitions occurred through fundamentally different paths (see the discussion therein, especially regarding the role of the hydrological cycle).

Figure 4b presents the optimal transition paths W→SB and SB→W in a three-dimensional projection, where we add, as a third coordinate, the variable , which indicates the fraction of the surface that has subfreezing temperatures (𝒯< 273.15 K). On the sides of the figure, two two-dimensional projections on the (𝒯,Δℐ) and on the (Δ𝒯,ℐ) planes are shown. Here, darker brown shadings indicate the higher density of points and the red and blue dots sample the highest probability for the SB→W and W→SB transitions paths, respectively. One could argue that the presence of an intersection between the SB→W and W→SB highest probability transition paths in Fig. 4a could have been a simple effect of 2D projection. Instead, we see here that the SB→W and W→SB most probable transition paths also cross in the 3D projection in a well-defined region, which indeed corresponds to the M state (pink square).

4.2.2 Lévy noise

There is scarcity of rigorous mathematical results regarding the weak noise limit of the transition paths between competing states in metastable stochastic systems forced by multiplicative Lévy noise. Indeed, the derivation of analytical results for this type of system largely remains an open problem. Recently, for stochastic partial differential equations with additive Lévy and Gaussian noise, the Onsager–Machlup action functional has been derived in Hu and Duan (2020), leading to a precise formulation of the most probable transition paths. Hence, we do not have solid mathematical results to interpret what we describe below, where, instead, we need to use heuristic arguments. As far as we know, this is the first attempt to estimate the most probable transition pathway between the metastable states in infinite stochastic systems with multiplicative pure Lévy process.

Figure 5A two-dimensional projection of the invariant measure on T,ΔT for different choices of α for Lévy noise. (a) α=0.5 and ε=0.0001, (b) α=1 and ε=0.004, and (c) α=1.5 and ε=0.01. The blue (red) line corresponds to the W→SB (SB→W) most probable transition path, and the arrows show the direction of the transitions. The coloured points indicate the deterministic attractors of the SB (blue), M (green), and W (red) states. In the inset of the left top corner plot in panels (a–c), the dark blue (red) contours are the ensembles of the transition paths between W→SB (SB→W).


A striking feature in Fig. 5 is that the invariant measure and the structure of the most probable transition paths (SB→W and W→SB), in the weak noise limit, are fundamentally different between the Lévy case and the Gaussian one. The invariant measure is highly peaked (dark red in the colour scheme) in a small region around the deterministic attractors, as, most typically, the Lévy noise fluctuations of μ̃ are very small. Additionally, the most probable transition paths depend very weakly on the chosen value for the stability parameter α. This suggests that the geometry of most probable path of transitions does not depend on the frequency and height of the Lévy diffusion jumps but rather on the qualitative fact that we are considering jump processes. Note that each panel of Fig. 5 is computed using data coming from the weakest noise considered for the corresponding value of α (see Table D1).

The W→SB most probable transition path is characterised by the simultaneous decrease in both T and Δ𝒯. This implies that the jump leads to a rapid and direct freezing of the whole planet. The stochastically averaged path crosses the basin boundary far from the M state. The most probable SB→W transition follows, instead, a path that is somewhat similar to the one found for the Gaussian case. We then argue that the closest region in the basin boundary to SB attractor is not too far from the M state. Further insight on difference between the Gaussian and Lévy case can be found by looking at the animations included in Lucarini et al. (2022).

Figure 6Comparison between the supercritical and subcritical paths that are constricted, using singular perturbations and the ensemble of trajectories corresponding to Lévy noise-induced transitions. (a) SB→W transitions. (b) W→SB transitions. Singular perturbations are where the subcritical paths are depicted as thin solid and dashed lines. The supercritical paths, leading to transition, are depicted as thick solid and dashed lines. See the text for further details. Here, α=1.0 and ϵ=0.004.


Figure 7Same as Fig. 6, respectively, using the same 3D projection, lines, and colour coding as in Fig. 4. The subcritical and supercritical paths are depicted as thin and thick solid lines, respectively. The corresponding dashed lines shown in Fig. 6 are not reported here, for simplicity. Here, α=1.0 and ϵ=0.004.


4.2.3 Lévy noise – singular perturbations

Based on what is discussed in Sect. 3.2, we expect that the transitions occur through the nearest region to the outgoing attractor in the basin boundary. We now try to clarify the properties of the most probable escape paths in the Lévy noise case by considering an additional set of simulations, taking inspiration from the edge-tracking algorithm (Skufca et al.2006). The idea is to exploit the fact that large jumps drive the transitions in the Lévy noise case. Starting from the deterministic SB state, we apply, in Eq. (6), singular perturbations of the form μμ+κδ(t) and bracket the critical value κcritSBW, leading to a transition to the W state as κcritSBW[κcritSBW,s,κcritSBW,u], where the simulation performed with the supercritical (subcritical) value of κ=κcritSBW,u (κ=κcritSBW,s) asymptotically reaches the competing (comes back to the initial) steady state. We obtain κcritSBW,s1.149  y and κcritSBW,u1.1492  y. Starting from the W state, we follow a similar procedure and find κcritWSB[κcritWSB,u,κcritWSB,s]. We obtain κcritWSB,s-1.3458  y and κcritWSB,u-1.346  y. In both cases, the value of κ of the supercritical and subcritical paths differs by δκ≈0.0002  y.

The projections on the 2D phase space spanned by (𝒯,Δ𝒯) of the supercritical and subcritical paths corresponding to the SB→W (W→SB) transition are shown in Fig. 6a (Fig. 6b), using the thick and thin black dashed lines, respectively. The basin boundary is indicated by a cyan (magenta) line for the SB→W (W→SB) transition. The steps to estimate the basin boundary are presented in Appendix C. Note that we are using, as background, the invariant measure and the subset of transitions referring to Lévy noise simulations performed using α=1.0. Nonetheless, what is discussed below would apply equally well had we chosen to consider as background, instead, data coming from the simulation performed with α=0.5 or α=1.5. Indeed, for the link we propose between transitions due to the Lévy noise and the case of singularly perturbed trajectories, what matters depends on the discontinuous nature of the Lévy noise.

In Fig. 6a (Fig. 6b) the supercritical and subcritical paths are superimposed on the ensemble of the trajectories of the SB→W (W→SB) transitions due to Lévy noise. The lines are better visible in the insets. By construction, after the perturbation is applied, the supercritical and subcritical orbits are close to the basin boundary. Hence, they are attracted towards the M state before being repelled towards the final asymptotic state. For comparison, we also portray, for both the SB→W and the W→SB case, an additional pair of supercritical and subcritical paths that are constructed using values of κ that differ by δκ=0.3  y, which is a much larger difference than the one mentioned previously (for the dashed lines). The paths depicted as thick (thin) solid lines cross (do not cross) the basin boundary. When looking at the W→SB transitions due to Lévy noise, we understand that, if the perturbation sends the orbit near the basin boundary, then the subsequent evolution of the system follows the supercritical paths. Of course, the Lévy perturbation often overshoots the basin boundary. In this case, after the transition, the orbit is not necessarily attracted towards the M state, whereas it converges directly to the final SB state. The signature of the attracting influence of the M state persists in the stochastically averaged transition trajectory. Note the bending towards higher values of Δ𝒯 before the eventual convergence to the SB state. In the case of the SB→W transitions, a similarly good correspondence between the supercritical path and the stochastically averaged transition trajectory is found.

Further support to the viewpoint proposed here is given by Fig. 7a and b, which are constructed along the lines of Fig. 4b and portray the supercritical and subcritical paths and the stochastically averaged transition trajectory realised via Lévy noise, with α=1.0 for the SB→W and W→SB case, respectively.

The transitions shown in Figs. 6 and 7 have been obtained by considering a discrete approximation of Dirac's δ, where the forcing acts at a constant value for 1  y. Specifically, Dirac's δ(t) is approximated as Δτ(t), where Δτ(t)=1/τ, if 0<t<τ, and vanishes elsewhere. The results are virtually unchanged if one considers τ<1  y because the main dynamical processes of the re-equilibration of the system act on longer timescales. The effect of the negative feedbacks of the system starts to become apparent when considering slower perturbations, lasting 2 or more years. Indeed, the resilience of the system to transitions is reduced when, ceteris paribus, faster perturbations are considered (see the analyses in this direction dealing with stability of the large-scale ocean circulation; Stocker and Schmittner1997; Lucarini et al.2005, 2007; Alkhayuon et al.2019). Nonetheless, also in this case, the agreement with the results presented in Figs. 6 and 7 is considerable. We remark that considering longer-lasting perturbations allows one to observe W→SB transitions without using, at any time, (unphysical) negative values for the solar irradiance. This is reassuring in terms of the robustness and of the physical sense of our results. Further details on the impact of considering different values of τ are reported in Appendix C.

5 Conclusions

It is a well-known that, as a result of the competition between the Boltzmann stabilising feedback and the ice–albedo destabilising feedback, under current astronomical and astrophysical conditions, the climate system is multistable, as at least two competing and distinct climates are present, i.e. the W and the SB. More recent investigations indicate that the partition of the phase space of the climate system might be more complex, as more than two asymptotic states might be present, some of them, possibly, associated with small basins of attraction.

For deterministic multistable systems, the asymptotic state of an orbit depends uniquely on the initial condition, and, specifically, on which basin of attraction it belongs to. The presence of stochastic forcing allows for transitions to occur between competing basins, thus giving rise to the phenomenon of metastability. Gaussian noise as a source of stochastic perturbations has been widely studied by the scientific community in recent years and provided very fruitful insight into the multiscale nature of the climatic time series. However, it has become apparent that more general classes of α-stable Lévy noise laws might also be suitable for modelling the observed climatic phenomena. In this regard, it is important to achieve a deeper understanding of the possible noise-induced transitions between competing stable climate states under α-stable Lévy perturbations and compare them with the Gaussian case.

As a starting point in this direction, we have studied the influence of different noise laws on the metastability properties of the randomly forced Ghil–Sellers EBM, which is governed by a nonlinear, parabolic, reaction–diffusion PDE. In the deterministic version of the model, we have three steady-state solutions, i.e. two stable, attractive climate states and one unstable saddle, corresponding to the edge state. The stable states correspond to the well-known W and SB climates. There is a fundamental dichotomy in the properties of the noise-induced transitions determined by whether we consider a stochastic forcing of intensity ε with a Gaussian versus an α-stable Lévy noise law. Note that, instead, the spatial structure of the noise is unchanged. This indicates that the phenomenology associated with the metastable behaviour depends critically on the choice of the noise law. Not many studies have investigated, numerically or through mathematical theory, the properties of transitions in metastable systems driven by multiplicative Lévy noise, as done here.

First, in the weak noise limit ε→0, the mean residence times inside either competing basin of attraction for diffusions driven by Gaussian vs. Lévy noise have a fundamentally different dependence on ε. Our results show that the logarithm of the mean residence time for Gaussian diffusions scales with ε−2, while, instead, a much weaker dependence is found for the Lévy case. Indeed, we find that the mean residence time is proportional to εα, where α is the stability parameter of the noise law. This result is in agreement with what has been proven in some special cases for additive Lévy noise and might indicate that these scaling properties are more general than usually assumed. We propose a simple argument based on approximating the Lévy noise as a composed Poisson process to support the applicability of the result in general circumstances, but, clearly, detailed mathematical investigations in this direction are needed.

Second, the results obtained for the most probable transition paths confirm that, in the weak noise limit, escapes from basins of attraction driven by Gaussian noise take place through the edge state. Additionally, instantonic and relaxation portions within each basin of attraction are clearly distinct, indicating nonequilibrium conditions that are, yet, qualitatively similar. In turn, Lévy diffusions leave the basin through the boundary region closest to the outgoing attractor, which seems to be the vicinity of the edge state when the thawing transition is considered. The freezing transition, instead, proceeds along a path that is fundamentally different. Finally, the most probable transition paths for the Lévy case appear to depend very weakly on the value of the stability parameter α, but seem, instead, to be determined by the nature of the Lévy noise of being a jump process. Indeed, we suggest that these properties can be better understood by considering that, to a first approximation, the transitions due to the Lévy diffusion correspond to supercritical paths associated with Dirac's δ-like singular perturbations to the solar irradiance. This viewpoint seems of general relevance to other problems where Lévy noise is responsible for exciting transitions between competing metastable states.

Our findings provide strong evidence that choosing noise laws other than Gaussian leads to fundamental changes in the metastability properties of a system, both in terms of statistics of the transitions between competing basins of attraction and most probable paths for such transitions. Leaving the door open for general noise laws might be relevant, both for interpreting observational data and for performing modelling exercises for the climate system and complex systems in general.

Let us give an example of the impact of making a wrong assumption on the nature of the acting stochastic forcing. Were we to naively interpret one of the panels in Fig. 5 as resulting from the dynamics of a dynamical system perturbed by Gaussian noise, then we would have to conclude that the unperturbed deterministic system possesses at least two edge states on the basin boundary separating the competing basins of attraction (see Margazoglou et al.2021, for a case where this situation applies). Hence, we would infer fundamentally wrong properties on the geometry of the phase space. Additionally, we would infer fundamentally different properties for the drift term.

Recent developments in data-driven methods based on the formalism of the Kramers–Moyal equation allow one to test accurately whether data are compatible with the hypothesis that stochasticity in the dynamics enters as a result of Gaussian noise or more general form of random forcing (Rydin Gorjão et al.2021b; Li and Duan2022). Indeed, we point the reader to the recent contribution by Rydin Gorjão et al. (2021a), which shows that the analysis of proxy climatic datasets indicates the need to go beyond Langevin equation-based modelling, as they discover that it is necessary to treat noise as the sum of continuous and discontinuous processes. This indicates the need to consider, in future modelling exercises, the possibility of investigating the properties of metastable systems where the stochastic forcing comes as the result of simultaneous Gaussian and α-stable Lévy noise perturbations.

Appendix A: Stochastic perturbations of Lévy type

In this section, we provide a summary of the basic properties of a symmetric α-stable Lévy process in a Hilbert space in which the solutions to SPDE (Eq. 13) are defined. We also repeat some properties in the n space that are more familiar to a wide audience of readers. It is pertinent to refer to the distribution law of Lévy increments, its characteristic function, the Lévy–Itô decomposition, and the Lévy jump measure for a deeper study of the metastable behaviour of the stochastic climate system (Eq. 13). Let (Ω,F,P) be a given complete probability space and H(,,) a separable Hilbert space with the norm and inner product ,. A stochastic process (Lα(t)t≥0) is a symmetric α-stable Lévy process in H if it satisfies the following:

  1. Lα(0)=0, almost surely.

  2. Independent increments – for any n∈ℕ and 0t1<t2<<tn-1<tn, the vector is as follows:

    (A1) ( L α ( t 1 ) - L α ( t 0 ) , , L α ( t n ) - L α ( t n - 1 ) ) ,

    where there is a family of independent random vectors in H.

  3. Stationary increments – for 0 l<t random vectors Lα(t)−Lα(l) and Lα(tl) have the same law 𝔏(.) in H, as follows:

    (A2) L ( L α ( t ) - L α ( l ) ) = L ( L α ( t - l ) ) .

    This law in n is a symmetric α-stable distribution L(.)=Sα((t-l)1α,0,0), i.e. zero skewness and shift parameters, with a stability parameter α(0,2] and a scaling parameter (t-l)1α. The stable distribution by the generalised central limit theorem (Schertzer and Lovejoy1997) is a limit in the distribution as n→∞ of the normalised sum Yn=1Bni=1n(Xi-Mn) of n independent identically distributed random variables Xi, with a common probability distribution function F(x), which does not necessarily have to have moments of both the first and second order. A necessary and sufficient condition for this is as follows (Keller and Kuske2000; Burnecki et al.2015):

    (A3) F ( x ) = [ c 1 + r 1 ( x ) ] | x | - α , x < 0 , = 1 - [ c 2 + r 2 ( x ) ] x - α , x > 0 ,

    with 0<α2, c1, and c2 positive constants, r1(x)→0 as x- and r2(x)→0 as x+. When this condition holds and α=2, then we can set Bn=h(n), where h(n) satisfies h2=nln h, and the stable distribution is just the Gaussian law.

  4. Sample paths are continuous in probability, i.e. for any t 0 and η>0, as follows:

    (A4) lim l t P ( L α ( t ) - L α ( l ) > η ) = 0 .

    For α(0,2) the symmetric α-stable Lévy process in n has a characteristic function of the following form:

    (A5) E e i u , L α ( t ) = e - C ( α ) t | | u | | α , u R n , t 0 ,

    where C(α)=π-1/2Γ((1+α)/2)Γ(n/2)Γ((n+α)/2), and Γ(.) is the Gamma function. In the case where α=2, we set C(2)=1/2 and (A5) becomes the characteristic function of a standard Brownian motion. However, the Brownian motion cannot be seen as a weak limit of α-stable Lévy process because of the divergence C(α)→∞ as α→2. The properties of the sample paths of Lα(t) are, in fact, quite different for α=2 and α<2. First, the α-stable Lévy process is a discontinuous, pure jump process, while the Brownian motion has continuous paths. Second, the Brownian motion has moments of all orders, whereas E|Lα(t)|γ< if γ<α. It can also be proved that the tails of Lα(t) are heavy, i.e. P(Lα(t)>u)u-α,u , which is quite the opposite of the exponentially light Gaussian tails. Moreover, for α(0,1), the path variation in Lα(t) is bounded on finite time intervals and unbounded for α[1,2).

Although neither the incremental nor the marginal distributions of a Lévy process in general are representable by the elementary functions, the Lévy motion is completely determined by the Lévy–Khintchine formula, which specifies the characteristic function of the Lévy process.

If Lα(t) is a symmetric α-stable Lévy process in H, then, in the following:

  1. The characteristic function of the Lévy–Khintchine formula is as follows:


    where, in the following:

    (A6) ψ ( h ) = H e i h , y - 1 - i h , y 1 { 0 < y 1 } ν ( d y ) ,

    where 1S is the indicator function for a set S, taking 1 on S and 0, otherwise, and ν is a Borel measure (also called the Lévy jump measure) in H, for which H(1y2)ν(dy)< with 1y2=min{1,y2}. A Borel measure, in addition, can be defined as the expected value of the number of jumps of specified size Q in the unit time interval, i.e. ν(Q)=EN(1,Q)(ω), ω∈Ω.

  2. In the Lévy–Itô decomposition, for any sequence of positive radii rn→0 and On={yH|rn+1<||y||rn}, there exists a sequence of independent compensated compound Poisson processes (L¯n(t))t0, n0 in H, with jump measures νn(B)=ν(BOn) for B∈ℬ(H), the Borel σ-algebra in H, and n1, which satisfy , almost surely for all t0 as follows:


    If Lα(t) is a symmetric α-stable Lévy process in n with the generating triplet (0,0,να), then there exists an independent Poisson random measure N on R+×(Rn{0}) (quantifying the number of jumps of Lα(t)), such that, for each t0, the following applies:

    (A9) L α ( t ) = y < 1 y N ̃ ( t , d y ) + y 1 y N ( t , d y ) ,

    where Ñ(dt,dx)=N(dt,dx)-να(dx)dt is the compensated Poisson random measure, and να(dx) is the jump measure. The small y<1 (large y1) jumps are controlled by Ñ(t,dy) ( N(t,dy) ).

  3. Its Lévy jump measure ν is symmetric in the sense that ν(-Q)=ν(Q) for Q∈ℬ(H) and has the following specific geometry:

    (A10) ν ( Q ) = Q ν ( d y ) = Q d r r 1 + α σ ( d s ) ,

    where r=y and s=y/y and σ:B(B1(0))[0,) is an arbitrary finite Radon measure on the unit sphere of H. The jump measure for a symmetric α-stable Lévy motion Lα(t) in n is defined by the following:

    (A11) ν α ( d u ) = c ( n , α ) d u | | u | | n + α ,

    with the intensity constant c(n,α)=αΓ((n+α)/2)21-απn/2Γ(1-α/2), where Γ(.) is the Gamma function (see Duan2015, and Applebaum2009).

One can come to a more intuitive interpretation of the stability parameter α(0,2) variation. For smaller values of α, the process is characterised by higher jumps with a lower frequency. As α increases, jumps decrease in height, and the frequency of their occurrence increases.

Appendix B: Probabilistic theory for the Lévy noise-induced escape

We briefly recapitulate here the main ideas behind the proof given in Debussche et al. (2013) of how the mean residence time in the competing metastable states of stochastically perturbed Chafee–Infante reaction–diffusion PDE scales with the intensity ε of the additive L(t) α-stable Lévy noise that acts as stochastic forcing.

One proceeds by considering the decomposition in the driving Lévy process by regularly varying the jump measure ν into small ξε and large ηε jump components. Let ΔtL=L(t)-L(t-) denote the jump increment of L at time t0, and 1ερ for ε, ρ(0,1) the jump height threshold of L. The process ηε is a compound Poisson process consisting of all jumps of height ΔtL>ε-ρ, with the following intensity:

(B1) β ε = ν 1 ε ρ B 1 c ( 0 ) ε α ρ ,

and the jump probability measure outside the ball 1ερB1(0) is as follows:

(B2) ν 1 ε ρ B 1 c ( 0 ) / β ε ,

where B1(0) is a ball of unit radius in H centred at the origin. The occurrence time of a kth large jump is defined recursively by the following:

(B3) Z 0 = 0 , Z k = inf { t > Z k - 1 | Δ t L > ε - ρ } , k 1 .

The waiting times between successive ηtε jumps have an exponential distribution Zk-Zk-1Exp(βε).

Small jump processes ξε=L-ηε, due to the symmetry of Lévy measure ν, are a mean zero martingale in H with finite exponential moments. Probabilistic events causing small jumps in the stochastic solution of the system are not able to overcome the force of its deterministic stable state, and therefore, do not contribute to the exit from the basin of attraction. Formally, during the time between two large jumps tk=Zk-Zk-1, the solution of Eq. (13), following the deterministic path (Eq. 1), returns to a small vicinity of the stable equilibria ϕW/SB, as follows:

(B4) sup x D W / SB sup Z k - 1 t Z k T ( t ) - T ( t ) 0 for ε 0 .

When a first large jump occurs, the solution process moves to the neighbouring domain of attraction with the following probability:

(B5) P ϕ W / SB + ε Δ t i L D W / SB = P Δ t i L 1 ε D W / SB c - ϕ W / SB = ν 1 ε ( D W / SB ) c - ϕ W / SB 1 ε ρ B 1 c ( 0 ) ν ( 1 ε ρ B 1 c ( 0 ) ) ε α ( 1 - ρ ) .

This is the probability that, at time ti, there will be a jump increment ΔtiL that exceeds the distance between the attractor and its domain of attraction boundary, as expressed by the jump probability measure (Eq. B2). In the zero-noise limit, the mean residence time in a basin of attraction is given by the following:

(B6) E [ τ ( ε ) ] i = 1 E [ Z i ] P inf j : ϕ W / SB + ε Δ t j L D W / SB = i E [ t 1 ] P ϕ W / SB + ε Δ t 1 L D W / SB i = 1 i 1 - P ϕ W / SB + ε Δ t 1 L D W / SB i - 1 1 ε α ρ ε α ( 1 - ρ ) 1 ε α ( 1 - ρ ) 2 = 1 ε α ,

that is, by the sum of all the mean values of a large jump occurrence time times the probability that the minimum of a sample of size i of jump increments is sufficiently large to reach into the neighbouring domain of attraction. Thus, at the random time instant of large jumps, the solution process transitions, in an abrupt move, from one attractor to another. Such behaviour of the random dynamical system is known as a metastability.

In Debussche et al. (2013), it was proven that, in the timescale λ(ε)=ν1εB1c(0),ε>0 the metastable shifting of the diffusion process between neighbourhoods of the two attractors represents a continuous time Markov chain in a state space {ϕSB,ϕW}, with a transition rate matrix 𝔔, as follows:

(B7) Q = 1 μ B 1 c ( 0 ) - μ ( ( D SB - ϕ SB ) c ) μ ( ( D SB - ϕ SB ) c ) μ ( ( D W - ϕ W ) c ) - μ ( ( D W - ϕ W ) c ) ,

where μ(⋅) is the limit measure of ν.

Appendix C: Transitions induced by singular perturbations

In Sect. 4.2.3, and in particular Figs. 6 and 7, we have studied the effect of singular perturbations of a Lévy kick. The idea is that transitions in a system perturbed by Lévy noise are primarily driven by rare large jumps. By applying a singular perturbation of the form μμ+κδ(t) (where μ=μ0=1.05 throughout), we have been able to bracket the critical values κcritWSB,s allowing for transitions between the two attractors. The expression κδ(t) is approximated as κΔτ(t), where Δτ(t)=1/τ if 0<t<τ and vanishes elsewhere. In Sect. 4.2.3 the results are shown for τ=1 year.

Table C1Values of supercritical κcrit,u for SB→W (second column) and W→SB (third column) for different approximations Δτ(t) of Dirac's δ(t).

Download Print Version | Download XLSX

Figure C1Comparison between the supercritical paths constructed using singular perturbations of different duration τ, where, in the legend, m stands for month and y stands for years. (a) SB→W transitions. (b) W→SB transitions. Different values of τ correspond to different coloured lines. The basin boundaries are now depicted as thick black lines. Here, α=1.0 and ϵ=0.004.


We performed additional simulations to locate the supercritical and subcritical values of κ for τ= 1 month, 6 months, 2 years, and 4 years. The corresponding supercritical values of κcritSBW,u and κcritWSB,u are shown in Table C1. In Fig. C1, we plot the corresponding supercritical transition trajectories for the values of Table C1 at different durations. Notice that now we use coloured solid lines for the supercritical cases. To estimate the basin boundary, we record the final point of when the forcing was active, in the (𝒯,Δ𝒯) projected space, for each duration. This point for each case is particularly visible, in the insets of Fig. C1, as a rapid reflection of the trajectory, which then follows closely the basin boundary (depicted as a thick black line). The basin boundary we can explore through this procedure is then estimated by linking the points obtained when considering various values of τ. Notice that the estimated basin boundaries are slightly different when looking at the two SB→W and W→SB transitions, as the basin boundaries have folds than cannot be captured in the sampled 2-dimensional projection used in Fig. C1.

Finally, as stated earlier, from the third column of Table C1, we remark that, when considering forcings with durations of, for example, 2 years and longer, transitions from the W to the SB state can be achieved while retaining, at all times, a positive value for the solar irradiance because, while the forcing is active, its value is μ+κ/τ.

Appendix D: Estimates for the mean escape time

We report in Table D1 a summary of the statistics of the escape times from the W state and from the SB state for various choices of the noise law.

Table D1Estimates and 95 % confidence intervals for the mean escape time τ from the W state and from the SB state, for Lévy noise with (a) α=0.5, (b) α=1.0, (c) α=1.5, and (d) Gaussian noise. No indicates the average number of transitions occurring in 105 years of temporal evolution.

Download Print Version | Download XLSX

Data availability

All the data used to produce the figures contained in this paper are publicly available in the supplement (Lucarini et al.2022) through the data repository

Video supplement

Illustrative animations portraying noise-induced transitions can be found in Lucarini et al. (2022 We, furthermore, uploaded the relevant material to YouTube at the following: (last access: 5 May 2022), which is the Lévy SB→W; (last access: 5 May 2022), which is the Lévy W→SB; (last access: 5 May 2022), which is the Gaussian SB→W; and (last access: 5 May 2022), which is the Gaussian W→SB.

Author contributions

VL conceptualized the paper, developed the methodology, validated the results, conducted the data analysis. and took the lead role in writing and revising the paper. LS and GM conceptualized the paper, developed the methodology and software, and helped with the writing and revising of the paper.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Nonlinear Processes in Geophysics. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Centennial issue on nonlinear geophysics: accomplishments of the past, challenges of the future”. It is not associated with a conference.


The authors wish to thank Peter Ashwin, Reyk Börner, Jesus I. Diaz, Jinqiao Duan, Michael, Ghil, Tobias Grafke, Serafim Kalliadasis, Alessandro Laio, Xue-Mei Li, and Greg Pavliotis, for the useful exchanges on various topics covered in this paper. Valerio Lucarini wishes to thank Myles Allen, for suggesting that we look into 3D projections of the phase space when studying transitions paths. The authors acknowledge the support provided by the EU Horizon 2020 project TiPES (grant no. 820970). Valerio Lucarini acknowledges the support provided by the EPSRC project (grant no. EP/T018178/1). This paper is dedicated to K. Hasselmann.

Financial support

This research has been supported by the Horizon 2020 (TiPES; grant no. 820970) and the Engineering and Physical Sciences Research Council (grant no. EP/T018178/1).

Review statement

This paper was edited by Daniel Schertzer and reviewed by two anonymous referees.


Abbot, D. S., Voigt, A., and Koll, D.: The Jormungand global climate state and implications for Neoproterozoic glaciations, J. Geophys. Res.-Atmos., 116, D18103,, 2011. a

Alharbi, R.: Nonlinear parabolic stochastic partial differential equation with application to finance, Doctoral thesis (PhD), University of Sussex, Brighton, (last access:5 May 2022), 2021. a, b

Alkhayuon, H., Ashwin, P., Jackson, L. C., Quinn, C., and Wood, R. A.: Basin bifurcations, oscillatory instability and rate-induced thresholds for Atlantic meridional overturning circulation in a global oceanic box model, P. Roy. Soc. A-Math. Phy., 475, 20190051,, 2019. a

Applebaum, D.: Lévy Processes and Stochastic Calculus, Cambridge Studies in Advanced Mathematics, 2nd edn., Cambridge University Press,, 2009. a, b

Ashwin, P., Wieczorek, S., Vitolo, R., and Cox, P.: Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system, Philos. T. R. Soc. A, 370, 1166–1184,, 2012. a

Barker, S., Knorr, G., Edwards, R. L., Parrenin, F., Putnam, A. E., Skinner, L. C., Wolff, E., and Ziegler, M.: 800,000 Years of Abrupt Climate Variability, Science, 334, 347–351,, 2011. a

Bensid, S. and Díaz, J. I.: On the exact number of monotone solutions of a simplified Budyko climate model and their different stability, Discrete Cont. Dyn.-B, 24, 1033–1047, 2019. a, b

Benzi, R., Sutera, A., and Vulpiani, A.: The mechanism of stochastic resonance, J. Phys. A-Math. Gen., 14, L453–L457,, 1981. a

Bódai, T., Lucarini, V., Lunkeit, F., and Boschi, R.: Global instability in the Ghil–Sellers model, Clim. Dynam., 44, 3361–3381, 2015. a, b, c, d

Bouchet, F., Gawedzki, K., and Nardini, C.: Perturbative calculation of quasi-potential in non-equilibrium diffusions: a mean-field example, J. Stat. Phys., 163, 1157–1210,, 2016. a

Brunetti, M., Kasparian, J., and Vérard, C.: Co-existing climate attractors in a coupled aquaplanet, Clim. Dynam., 53, 6293–6308,, 2019. a

Budhiraja, A. and Dupuis, P.: A variational representation for positive functionals of infinite dimensional Brownian motion, Probability and Mathematical Statistics–Wroclaw University, 20, 39–61, 2000. a

Budhiraja, A., Dupuis, P., and Maroulas, V.: Large deviations for infinite dimensional stochastic dynamical systems, Ann. Probab., 36, 1390–1420,, 2008. a

Budyko, M. I.: The effect of solar radiation variations on the climate of the Earth, Tellus, 21, 611–619,, 1969. a, b

Burnecki, K., Wylomanska, A., and Chechkin, A.: Discriminating between Light- And heavy-tailed distributions with limit theorem, PLoS ONE, 10, e0145604,, 2015. a

Burrage, K. and Lythe, G.: Accurate stationary densities with partitioned numerical methods for stochastic partial differential equations, Stochastics and Partial Differential Equations: Analysis and Computations, 2, 262–280,, 2014. a

Cai, R., Chen, X., Duan, J., Kurths, J., and Li, X.: Lévy noise-induced escape in an excitable system, J. Stat. Mech. Theory E., 2017, 063503,, 2017. a

Chechkin, A., Sliusarenko, O., Metzler, R., and Klafter, J.: Barrier crossing driven by Levy noise: Universality and the Role of Noise Intensity, Phys. Rev. E, 75, 041101,, 2007. a

Cialenco, I., Fasshauer, G. E., and Ye, Q.: Approximation of stochastic partial differential equations by a kernel-based collocation method, Int. J. Comput. Math., 89, 2543–2561,, 2012. a

Dai, M., Gao, T., Lu, Y., Zheng, Y., and Duan, J.: Detecting the maximum likelihood transition path from data of stochastic dynamical systems, Chaos, 30, 113124,, 2020. a, b

Davie, A. M. and Gaines, J. G.: Convergence of numerical schemes for the solution of parabolic stochastic partial differential equations, Math. Comput., 70, 121–134,, 2000. a

Debussche, A., Högele, M., and Imkeller, P.: The dynamics of nonlinear reaction-diffusion equations with small lévy noise, in: Lecture Notes in Mathematics, Springer, Berlin,, 2013. a, b, c, d, e, f, g

Díaz, G. and Díaz, J. I.: Stochastic energy balance climate models with Legendre weighted diffusion and an additive cylindrical Wiener process forcing, Discrete Cont. Dyn.-S.,, 2021. a

Díaz, J. I., Hernández, J., and Tello, L.: On the Multiplicity of Equilibrium Solutions to a Nonlinear Diffusion Equation on a Manifold Arising in Climatology, J. Math. Anal. Appl., 216, 593–613,, 1997. a, b

Ditlevsen, P. D.: Observation of α-stable noise induced millennial climate changes from an ice-core record, Geophys. Res. Lett., 26, 1441–1444,, 1999. a, b, c

Doering, C. R.: A stochastic partial differential equation with multiplicative noise, Phys. Lett. A, 122, 133–139,, 1987. a, b

Duan, J.: An introduction to stochastic dynamics, Cambridge University Press, New York, 2015. a, b, c

Duan, J. and Wang, W.: Effective Dynamics of Stochastic Partial Differential Equations, Elsevier, Boston,, 2014. a, b

Dybiec, B. and Gudowska-Nowak, E.: Lévy stable noise-induced transitions: stochastic resonance, resonant activation and dynamic hysteresis, J. Stat. Mech. Theory E., 2009, P05004,, 2009. a

Fan, A. H.: Sur les chaos de Lévy stables d'indice 0<α<1, Ann. Sci. Math. Québec, 1, 53–66, 1997. a

Feudel, U., Pisarchik, A. N., and Showalter, K.: Multistability and tipping: From mathematics and physics to climate and brain–Minireview and preface to the focus issue, Chaos, 28, 033501,, 2018. a

Freidlin, M. I. and Wentzell, A. D.: Random perturbations of dynamical systems, Springer, New York, 1984. a, b

Gao, T., Duan, J., Kan, X., and Cheng, Z.: Dynamical inference for transitions in stochastic systems with α-stable Lévy noise, J. Phys. A-Math. Theor., 49, 294002,, 2016. a

Garain, K. and Sarathi Mandal, P.: Stochastic sensitivity analysis and early warning signals of critical transitions in a tri-stable prey–predator system with noise, Chaos, 32, 033115,, 2022. a

Ghil, M.: Climate stability for a Sellers-type model, J. Atmos. Sci., 33, 3–20,<0003:CSFAST>2.0.CO;2, 1976. a, b, c, d, e, f

Ghil, M.: Energy-Balance Models: An Introduction, in: Climatic Variations and Variability: Facts and Theories: NATO Advanced Study Institute First Course of the International School of Climatology, Ettore Majorana Center for Scientific Culture, Erice, Italy, March 9–21, 1980, edited by: Berger, A., Springer Netherlands, Dordrecht, 461–481,, 1981. a, b

Ghil, M.: A mathematical theory of climate sensitivity or, How to deal with both anthropogenic forcing and natural variability?, in: Climate Change: Multidecadal and Beyond, edited by Chang, P. C., Ghil, M., Latif, M., and Wallace, J. M., World Scientific/Imperial College Press, 31–51, 2015. a

Ghil, M. and Childress, S.: Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics, Springer-Verlag, Berlin, 1987. a

Ghil, M. and Lucarini, V.: The physics of climate variability and climate change, Rev. Mod. Phys., 92, 035002,, 2020. a, b

Gottwald, G.: A model for Dansgaard-Oeschger events and millennial-scale abrupt climate change without external forcing, Clim. Dynam., 56, 227–243,, 2021. a, b

Gottwald, G. A. and Melbourne, I.: Homogenization for deterministic maps and multiplicative noise, P. Roy. Soc. A-Math. Phy., 469, 20130201,, 2013. a

Gould, S. J.: Wonderful Life: The Burgess shale and the Nature of History, W.W. Norton, New York, 1989. a

Grafke, T. and Vanden-Eijnden, E.: Numerical computation of rare events via large deviation theory, Chaos, 29, 063118,, 2019. a

Grafke, T., Grauer, R., and Schäfer, T.: The instanton method and its numerical implementation in fluid mechanics, J. Phys. A-Math. Theor., 48 333001,, 2015. a, b

Grafke, T., Schäfer, T., and Vanden-Eijnden, E.: Long Term Effects of Small Random Perturbations on Dynamical Systems: Theoretical and Computational Tools, in: Recent Progress and Modern Challenges in Applied Mathematics, Modeling and Computational Science, edited by: Melnik, R., Makarov, R., and Belair, J., Fields Institute Communications, Springer, New York, NY,, pp. 17–55, 2017. a

Graham, R.: Macroscopic potentials, bifurcations and noise in dissipative systems, in: Fluctuations and Stochastic Phenomena in Condensed Matter, edited by: Garrido, L., Springer Berlin Heidelberg, 1–34, ISBN',978-3-540-47401-2, 1987. a, b

Graham, R., Hamm, A., and Tél, T.: Nonequilibrium potentials for dynamical systems with fractal attractors or repellers, Phys. Rev. Lett., 66, 3089–3092,, 1991. a, b

Grebogi, C., Ott, E., and Yorke, J. A.: Fractal Basin Boundaries, Long-Lived Chaotic Transients, and Unstable-Unstable Pair Bifurcation, Phys. Rev. Lett., 50, 935–938,, 1983. a

Grigoriu, M. and Samorodnitsky, G.: Dynamic Systems Driven by Poisson/Lévy White Noise, in: IUTAM Symposium on Nonlinear Stochastic Dynamics, edited by: Namachchivaya, N. S. and Lin, Y. K., Springer Netherlands, Dordrecht, 319–330,, 2003. a

Hänggi, P.: Escape from a metastable state, J. Stat. Phys., 42, 105–148, 1986. a

Hasselmann, K.: Stochastic climate models, Part I. Theory, Tellus, 28, 473–485, 1976. a

Hetzer, G.: The structure of the principal component for semilinear diffusion equations from energy balance climate models, Houston J. Math., 16, 203–216, 1990. a, b

Hoffman, P. F., Kaufman, A. J., Halverson, G. P., and Schrag, D. P.: A Neoproterozoic Snowball Earth, Science, 281, 1342–1346,, 1998. a

Hu, J. and Duan, J.: Onsager-Machlup action functional for stochastic partial differential equations with Levy noise, arXiv [preprint],, 4 December 2020. a, b, c

Imkeller, P. and Pavlyukevich, I.: First exit times of SDEs driven by stable L'evy processes, Stoch. Proc. Appl., 116, 611–642,, 2006a. a, b, c, d

Imkeller, P. and Pavlyukevich, I.: Lévy flights: transitions and meta-stability, J. Phys. A-Math. Gen., 39, L237–L246,, 2006b. a, b, c, d

Imkeller, P. and von Storch, J. S.: Stochastic Climate Models, Birkhauser, Basel, 2001. a

Jentzen, A. and Kloeden, P. E.: The numerical approximation of stochastic partial differential equations, Milan J. Math., 77, 205–244,, 2009. a

Kaper, H. and Engler, H.: Mathematics and climate, SIAM, Philadelphia, 2013. a, b

Keller, J. and Kuske, R.: Rate of convergence to a stable law, SIAM J. Appl. Math., 61, 1308–1323,, 2000. a

Kloeden, P. E. and Shott, S.: Linear-implicit strong schemes for itô-galkerin approximations of stochastic PDES, Journal of Applied Mathematics and Stochastic Analysis, 14, 697341,, 2001. a

Kramers, H.: Brownian motion in a field of force and the diffusion model of chemical reactions, Physica, 7, 284–304,, 1940. a

Kraut, S. and Feudel, U.: Multistability, noise, and attractor hopping: The crucial role of chaotic saddles, Phys. Rev. E, 66, 015207,, 2002. a

Kuhwald, I. and Pavlyukevich, I.: Stochastic Resonance in Systems Driven by α-Stable Lévy Noise, International Conference on Vibration Problems 2015, Procedia Engineer., 144, 1307–1314,, 2016. a

Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping elements in the Earth's climate system, P. Nat. Acad. Sci. USA, 105, 1786–1793, 2008. a

Lewis, J. P., Weaver, A. J., and Eby, M.: Snowball versus slushball Earth: Dynamic versus nondynamic sea ice?, J. Geophys. Res., 112, C11014,, 2007. a

Li, Y. and Duan, J.: Extracting Governing Laws from Sample Path Data of Non-Gaussian Stochastic Dynamical Systems, J. Stat. Phys., 186, 30,, 2022. a

Linsenmeier, M., Pascale, S., and Lucarini, V.: Climate of Earth-like planets with high obliquity and eccentric orbits: Implications for habitability conditions, Planet. Space Sci., 105, 43–59,, 2015. a

Lovejoy, S. and Schertzer, D.: The Weather and Climate: Emergent Laws and Multifractal Cascades, Cambridge University Press, Cambridge, 2013. a

Lu, Y. and Duan, J.: Discovering transition phenomena from data of stochastic dynamical systems with Lévy noise, Chaos, 30, 093110,, 2020. a, b, c

Lucarini, V. and Bódai, T.: Edge states in the climate system: exploring global instabilities and critical transitions, Nonlinearity, 30, R32–R66,, 2017. a, b, c

Lucarini, V. and Bódai, T.: Transitions across Melancholia States in a Climate Model: Reconciling the Deterministic and Stochastic Points of View, Phys. Rev. Lett., 122, 158701,, 2019. a, b, c, d, e

Lucarini, V. and Bódai, T.: Global stability properties of the climate: Melancholia states, invariant measures, and phase transitions, Nonlinearity, 33, R59–R92,, 2020. a, b, c, d, e, f, g

Lucarini, V., Calmanti, S., and Artale, V.: Destabilization of the thermohaline circulation by transient changes in the hydrological cycle, Clim. Dynam., 24, 253–262,, 2005. a

Lucarini, V., Calmanti, S., and Artale, V.: Experimental mathematics: Dependence of the stability properties of a two-dimensional model of the Atlantic ocean circulation on the boundary conditions, Russ. J. Math. Phys., 14, 224–231,, 2007. a

Lucarini, V., Fraedrich, K., and Lunkeit, F.: Thermodynamic analysis of snowball Earth hysteresis experiment: Efficiency, entropy production and irreversibility, Q. J. Roy. Meteor. Soc., 136, 2–11,, 2010. a

Lucarini, V., Pascale, S., Boschi, R., Kirk, E., and Iro, N.: Habitability and Multistability in Earth-like Planets, Astron. Nachr., 334, 576–588,, 2013. a

Lucarini, V., Blender, R., Herbert, C., Ragone, F., Pascale, S., and Wouters, J.: Mathematical and physical ideas for climate science, Rev. Geophys., 52, 809–859,, 2014a. a

Lucarini, V., Serdukova, L., and Margazoglou, G.: Lévy-noise versus Gaussian-noise-induced Transitions in the Ghil-Sellers Energy Balance Model, figshare,, 2022. a, b, c, d, e, f

Margazoglou, G., Grafke, T., Laio, A., and Lucarini, V.: Dynamical landscape and multistability of a climate model, P. Roy. Soc. A-Math. Phy., 477, 20210019,, 2021. a, b, c, d, e, f, g, h, i

Millàn, H., Cumbrera, R., and Tarquis, A. M.: Multifractal and Levy-stable statistics of soil surface moisture distribution derived from 2D image analysis, Appl. Math. Model., 40, 2384–2395,, 2016. a

Nicolis, C.: Stochastic aspects of climatic transitions – response to a periodic forcing, Tellus, 34, 308–308,, 1982. a

North, G. and Stevens, M.: Energy-balance climate models, in: Frontiers of Climate Modeling, edited by: Kiehl, J. T. and Ramanathan, V., Cambridge University Press, Cambridge, 52–72,, 2006. a, b

North, G. R.: Multiple solutions in energy balance climate models, Global Planet. Change, 2, 225–235,, 1990. a, b

North, G. R., Cahalan, R. F., and Coakley Jr., J. A.: Energy balance climate models, Rev. Geophys., 19, 91–121,, 1981. a, b

Ott, E.: Chaos in dynamical systems, 2nd edn., Cambridge University Press, Cambridge,, 2002. a

Pavliotis, G. and Stuart, A.: Multiscale methods, Texts in applied mathematics, Springer, New York, NY, 2008. a, b

Peixoto, J. P. and Oort, A. H.: Physics of Climate, AIP Press, New York, New York, 1992. a

Penland, C. and Sardeshmukh, P. D.: Alternative interpretations of power-law distributions found in nature, Chaos, 22, 023119,, 2012. a

Peszat, S. and Zabczyk, J.: Stochastic Partial Differential Equations with Levy Noise: An Evolution Equation Approach, Cambridge University Press,, 2007. a, b

Pierrehumbert, R., Abbot, D., Voigt, A., and Koll, D.: Climate of the Neoproterozoic, Annu. Rev. Earth Pl. Sc., 39, 417–460,, 2011. a, b

Ragon, C., Lembo, V., Lucarini, V., Vérard, C., Kasparian, J., and Brunetti, M.: Robustness of Competing Climatic States, J. Climate, 35, 2769–2784,, 2022. a

Rhodes, R., Sohier, J., and Vargas, V.: Levy multiplicative chaos and star scale invariant random measures, Ann. Probab., 42, 689–724,, 2014. a

Risken, H.: The Fokker–Planck equation, Springer, Berlin, 1996. a

Rydin Gorjão, L., Riechers, K., Hassanibesheli, F., Witthaut, D., Lind, P. G., and Boers, N.: Changes in stability and jumps in Dansgaard–Oeschger events: a data analysis aided by the Kramers–Moyal equation, Earth Syst. Dynam. Discuss. [preprint],, in review, 2021a. a

Rydin Gorjão, L., Witthaut, D., Lehnertz, K., and Lind, P. G.: Arbitrary-Order Finite-Time Corrections for the Kramers-Moyal Operator, Entropy, 23, 517,, 2021b. a

Saltzman, B.: Dynamical Paleoclimatology: Generalized Theory of Global Climate Change, Academic Press New York, New York, 2001. a

Schertzer, D. and Lovejoy, S.: Multifractal simulations and analysis of clouds by multiplicative processes, Atmos. Res., 21, 337–361,, 1988. a

Schertzer, D. and Lovejoy, S.: Universal multifractals do exist!: Comments on “a statistical analysis of mesoscale rainfall as a random Cascade”, J. Appl. Meteorol., 36, 1296–1303,<1296:UMDECO>2.0.CO;2, 1997. a, b

Schertzer, D., Larchevêque, M., Duan, J., Yanovsky, V. V., and Lovejoy, S.: Fractional Fokker–Planck equation for nonlinear stochastic differential equations driven by non-Gaussian Lévy stable noises, J. Math. Phys., 42, 200,, 2001. a

Schmitt, F., Lovejoy, S., and Schertzer, D.: Multifractal analysis of the Greenland Ice-Core Project climate data, Geophys. Res. Lett., 22, 1689–1692,, 1995. a

Schmitt, F., Schertzer, D., Lovejoy, S., and Brunet, Y.: Multifractal temperature and flux of temperature variance in fully developed turbulence, Europhys. Lett., 34, 195–200,, 1996. a

Sellers, W. D.: A global climatic model based on the energy balance of the earth-atmosphere system, J. Appl. Meteorol., 8, 392–400, 1969. a, b

Serdukova, L., Zheng, Y., Duan, J., and Kurths, J.: Metastability for discontinuous dynamical systems under Lévy noise: Case study on Amazonian Vegetation, Sci. Rep.-UK, 7, 9336,, 2017. a

Singla, R. and Parthasarathy, H.: Quantum robots perturbed by Levy processes: Stochastic analysis and simulations, Commun. Nonlinear Sci., 83, 105142,, 2020. a

Skufca, J. D., Yorke, J. A., and Eckhardt, B.: Edge of Chaos in a Parallel Shear Flow, Phys. Rev. Lett., 96, 174101,, 2006. a, b, c

Solanki, S. K., Krivova, N. A., and Haigh, J. D.: Solar Irradiance Variability and Climate, Annu. Rev. Astron. Astr., 51, 311–351,, 2013. a

Steffen, W., Rockström, J., Richardson, K., Lenton, T. M., Folke, C., Liverman, D., Summerhayes, C. P., Barnosky, A. D., Cornell, S. E., Crucifix, M., Donges, J. F., Fetzer, I., Lade, S. J., Scheffer, M., Winkelmann, R., and Schellnhuber, H. J.: Trajectories of the Earth System in the Anthropocene, P. Nat. Acad. Sci. USA, 115, 8252–8259,, 2018. a

Stocker, T. F. and Schmittner, A.: Influence of CO2 emission rates on the stability of the thermohaline circulation, Nature, 388, 862–865,, 1997. a

Tessier, Y., Lovejoy, S., and Schertzer, D.: Universal multifractals: theory and observations for rain and clouds, J. Appl. Meteorol. Clim., 32, 223–250,<0223:UMTAOF>2.0.CO;2, 1993. a

Thompson, W. F., Kuske, R. A., and Monahan, A. H.: Reduced α-stable dynamics for multiple time scale systems forced with correlated additive and multiplicative Gaussian white noise, Chaos, 27, 113105,, 2017. a

Varadhan, S. R. S.: Large deviations and applications, Society for Industrial and Applied Mathematics Philadelphia, 75 pp.,, 1985. a

Voigt, A. and Marotzke, J.: The transition from the present-day climate to a modern Snowball Earth, Clim. Dynam., 35, 887–905,, 2010. a

Vollmer, J., Schneider, T. M., and Eckhardt, B.: Basin boundary, edge of chaos and edge state in a two-dimensional model, New J. Phys., 11, 013040,, 2009. a

Weron, A. and Weron, R.: Computer simulation of Levy alpha-stable variables and processes, Chaos – The Interplay Between Stochastic and Deterministic Behaviour, edited by: Garbaczewski, P., Wolf, M., and Weron, A., Springer Berlin Heidelberg, Berlin, Heidelberg, 379–392, ISBN 978-3-540-44722-1, 1995. a

Wu, J., Xu, Y., and Ma, J.: Lévy noise improves the electrical activity in a neuron under electromagnetic radiation, PloS one, 12, e0174330–e0174330,, 2017.  a

Yagi, A.: Dynamical Systems, in: Abstract Parabolic Evolution Equations and their Applications, Springer Monographs in Mathematics, Springer Berlin Heidelberg,, 2010. a, b

Zheng, Y., Serdukova, L., Duan, J., and Kurths, J.: Transitions in a genetic transcriptional regulatory system under Lévy motion, Sci. Rep.-UK, 6, 29274,, 2016. a

Zheng, Y., Yang, F., Duan, J., Sun, X., Fu, L., and Kurths, J.: The maximum likelihood climate change for global warming under the influence of greenhouse effect and Lévy noise, Chaos, 30, 013132,, 2020. a, b, c

Short summary
In most of the investigations on metastable systems, the stochastic forcing is modulated by Gaussian noise. Lévy noise laws, which describe jump processes, have recently received a lot of attention, but much less is known. We study stochastic versions of the Ghil–Sellers energy balance model, and we highlight the fundamental difference between how transitions are performed between the competing warm and snowball states, depending on whether Gaussian or Lévy noise acts as forcing.