Transitions across Melancholia States in a Climate Model: Reconciling the Deterministic and Stochastic Points of View

The Earth is well-known to be, in the current astronomical configuration, in a regime where two asymptotic states can be realised. The warm state we live in is in competition with the ice-covered snowball state. The bistability exists as a result of the positive ice-albedo feedback. In a previous investigation performed on a intermediate complexity climate model we have identified the unstable climate states (Melancholia states) separating the co-existing climates, and studied their dynamical and geometrical properties. The Melancholia states are ice-covered up to the mid-latitudes and attract trajectories initialised on the basins boundary. In this paper, we study how stochastically perturbing the parameter controlling the intensity of the incoming solar radiation impacts the stability of the climate. We detect transitions between the warm and the snowball state and analyse in detail the properties of the noise-induced escapes from the corresponding basins of attraction. We determine the most probable paths for the transitions and find evidence that the Melancholia states act as gateways, similarly to saddle points in an energy landscape.

The Earth, for a vast range of parameters controlling its radiative budget, e.g. the intensity of the solar irradiance and the concentration of greenhouse gases, including the presentday astronomical configuration and atmospheric composition, supports two co-existing climates. One is the warm state we live in, and the other one is the snowball state, featuring global glaciation and extremely low surface temperatures [1,2]. Indeed, events of onset and decay of snowball conditions have taken place in the Neoproterozoic [3]. The bistability of the climate system comes from the competition between the positive ice-albedo feedback (ice reflects efficiently the solar radiation) and the negative Boltzmann feedback (a warmer surface emits more radiation), with the tippings point realised when the negative and positive feedbacks are equally strong, with ensuing loss of bistability. With simple models one can identify, within the bistability region, unstable solutions-Melancholia states -sitting in-between the two stable climates. Melancholia (M) states are, far from the tipping points, ice-covered up to the mid-latitudes. Small perturbations applied to trajectories initialised on the M states lead to the system falling into either asymptotic state [4][5][6]. Improving our understanding of the related critical transitions is a key challenge for geoscience and has strong implications in terms of planetary habitability [2,[7][8][9].
The goal of this letter is to explore, using a simplified yet Earth-like climate model, the phase space of the climate system by taking advantage of the rich dynamics resulting from adding stochastic perturbations, and, in particular, by focusing on noise-induced transitions between the warm (W) and snowball (SB) attractors and linking this with the global stability properties analysed in [9] using tools and ideas of highdimensional deterministic dynamical systems. The methodology proposed here is of general relevance for studying multistable systems [10] and, specifically, for studying in a novel way the properties of the Earth tipping elements [11].
Multistable systems are extensively investigated both in natural and social sciences [10] and they can be introduced as follows. We consider a smooth autonomous continuoustime dynamical system acting on a smooth finite-dimensional compact manifold M. We define x(t, x 0 )=S t (x 0 ) as an orbit at time t, where S t is the evolution operator, and x 0 the initial conditions at t = 0. We write the corresponding set of ordinary differential equations asẋ is a smooth vector field. The system is multistable if it possesses more than one asymptotic states, defined by the attractors Ω j , j = 1, . . . , J. The asymptotic state of the orbit is determined by its initial condition, and the phase space is partitioned between the basins of attraction B j of the attractors Ω j and the boundaries ∂B l , l = 1, . . . , L separating such basins. If the dynamics is determined by the energy landscape U (x), with F (x) = −∇U (x), the attractors are the local minima x j , j = 1, . . . , J of U (x), the basin boundaries ∂B l , l = 1, . . . , L are the mountain crests, which are smooth manifolds, each possessing a minimum energy saddle -a mountain pass. More generally,the basin boundaries can be strange geometrical objects with co-dimension smaller than one. Orbits initialized on the basin boundaries ∂B l , l = 1, . . . , L are attracted towards invariant saddles. Such saddles Π l , l = 1, . . . , L can feature chaotic dynamics [12][13][14][15]. The definition of the basin boundaries and of the saddles is key for understanding the global stability properties of the system and its global bifurcations.
In [9,16], we adapted the edge tracking algorithm presented in [17,18] and constructed in the bistable region the M states separating the two co-existing realisable W and SB climates. The critical transitions are associated to boundary crises [14] associated to collisions between the M state and one of the stable climates, and are flagged by a diverging linear response [19]. In [9] we constructed the M states for an intermediate complexity climate model with O(10 5 ) degrees of freedom. We showed that the M state has, in a range of values of the control parameter µ (ratio between the considered solar irradiance and the present-day value), chaotic dynamics, leading to weather variability and to a limited horizon of predictability; see the caption of Fig. 1. Since this instability is much faster than the climatic one due to the icealbedo feedback, the basin boundary is a fractal set with nearzero codimension, in agreement with results obtained in lowdimensional cases [12,20]. Near the basin boundary there is de facto no Lorenz' [21] predictability of the second kind 1 .
Here, building on [9], we study how a fluctuating solar irradiance can trigger transitions between the W and SB states, and investigate the typical paths of such transitions. The climate model is constructed by coupling the primitive equations atmospheric model PUMA [22] with the Ghil-Sellers energy balance model [6] (see also [23,24]), which describes succinctly oceanic heat transports. The ocean model describes effectively the ice-albedo feedback, and defines the slow manifold of the system. The coupling is realised by relaxing the atmospheric temperature to an adiabatic profile anchored to the ocean surface temperature, and by incorporating vertical heat fluxes. The ocean temperature T S (t, φ, λ), where φ is latitude and λ is longitude, evolves as follows: where S * is the present-day solar irradiance (the factor 4 comes from the Earth-Sun geometry [25]), and the heat capacity C and the geometrical factor I depend explicitly on φ. The albedo α depends on φ and, critically, on T S , with a rapid transition from high albedo for low values of T S (α max = 0.6) to low albedo for T S 260 K (α min = 0.2), which fuels the ice-albedo feedback. Finally, O is the outgoing radiation, , and increases with T S , (accounting for the Boltzmann feedback beside the greenhouse effect), D is a diffusion operator parametrizing the meridional heat transport, and χ describes the ocean-atmosphere heat exchange [9].
The stochastic perturbation modulates the solar irradiance through the factor (1 + σdW/dt), where σ controls the noise intensity noise, and dW is the increment of a Wiener process. The noise is multiplicative because dW/dt is multiplied by the factor 1 − α(φ, T S ) in Eq. 1. Adding a Gaussian random variable of variance σ 0 at each time step ∆t (1 hour), when numerically integrating the model, corresponds to having a relative fluctuation of the solar irradiance Noise-induced escapes from attractors have been intensely studied [26][27][28]. We frame our problem by considering a stochastic differential equation in Itô form written as In the case of nondegenerate additive noise and for a class of multiplicative noise laws, the Freidlin-Wentzell [29] theory and extensions thereof [20,30,31] show that in the weak-noise limit σ → 0 the invariant measure can be written as a large deviation law: where Z(x) is the pre-exponential factor and Φ(x) is the pseudo-potential 2 , which has local minima at the deterministic attractors Ω j , j = 1, . . . , J. Both the M states and the attractors can be chaotic: if so, Φ has constant value over each M state and each attractor, respectively [30,31]. The probability that an orbit with initial condition in B j does not escape from it over a time p decays as: whereτ σ is the expected escape time and where ∆Φ = Φ(Π l ) − Φ(Ω j ) is the pseudo-potential barrier height [20]; In general, one may need to add a correcting prefactor in Eq. 3 [20]. In the weak-noise limit, the transition paths follow the instantons, which are minimizers of the Freidlin-Wentzell action [27,28,32,33]. An instanton connects a point in Ω j to a point in Π l ; if these sets are not fixed points, the instanton is not unique, as the pseudo-potential is constant on Ω j and Π l . 3. An optimal algorithm for estimating ∆Φ is reported in [34].
The results given in Eqs. 2-3 apply in the case of multiplicative noise laws if one assumes that the noise correlation matrix is well-behaved, according to what discussed in [35,36]. We believe that in our system such a condition applies, essentially because the factor 1 − α(φ, T S ) is bounded in the phase space between 0.4 (α = α max = 0.6, ice-cover) and 0.8 (α = α min = 0.2, very warm conditions with absence of ice cover). In the phase space region near the SB attractor, we have that 1 − α(φ, T S ) ∼ 0.4 since the temperature T S is extremely low and the planet is fully glaciated, so that α(φ, T S ) is constant, with α(φ, T S ) ∼ α min . Near the W attractor, the properties of α(φ, T S ) are more complex, because only part of the planet is glaciated. If [X] the global average of the spatial field X, and Y is the long term average of Y , we have that, typically, ∂[ α ]/∂[ T S ] < 0, because a decrease in [ T S ] leads to moving the ice line equatorward, thus leading to higher average albedo. Then, near the W attractor, noise enhances the instability linked to the W → SB transition, and one expects that for finite noise the peak of the invariant measure is shifted to lower values of [T S ] with respect to the deterministic attractor.
The ratio of the variance of the noise in the W vs SB attractors can be estimated as where the typical albedo of the W (SB) attractor is α W ≈ 0.3 (α SB ≈ α max = 0.6). The two attractors have different microscopic (and macroscopic) temperatures.
We show now our results. We treat two cases inside the region of bistability depicted in Fig. 1, namely µ = 0.98 (close to the tipping point µ W →SB ) and µ = 1.0.
In the case of µ = 0.98, we consider noise intensities ranging from σ τ = 0.5% to σ τ = 1.4%, with τ = 100 years (y). For each value of σ τ , we initialise 50 orbits in the W basin of attraction and study the statistics of the escape times towards the SB attractor. When the transition takes place, we stop the integration. We observe (not shown) that for each value of σ τ the escape times are to a good approximation exponentially distributed, see Eq. 3. The expectation value of the transition timesτ σ is shown in Fig. 2. Indeed,τ σ obeys to a good approximation what shown in Eq. 3, so that the difference of the potential Φ is half of the slope of the straight line. For reference, we have that for σ 100y = 0.5% the average escape time is about 5.2 × 10 3 y. We can predict that the escape rate increases to about 1.2 × 10 7 y when σ 100y ∼ 0.3%.
We then look at the transition paths. Following [9,16], we choose to consider the reduced phase space spanned by [T S ] and by ∆T S , which is the difference between the spatial averages of T S in the latitudinal belts [0, 30 • N ] and [30 • N, 90 • N ], respectively. This reduced phase space provides a minimal yet physically informative viewpoint on the problem. Figure 3 depicts (σ τ = 1.0%), the transient twodimensional probability distribution function (pdf)ρ constructed using the above-described 50 simulations, where the statistics is collected only until the W → SB transition is realised. Note thatρ is not the invariant density of the system. The transitions typically take place along a very narrow band linking the W attractor and the M state 3 . We construct an estimate of the instanton associated to the W → SB transition by conditionally averaging the orbits according to the value of [T S ]. To a good approximation, the instanton connects the W attractor to the M state, and follows a path of decreasing probability. We remark that we do not find evidence of different paths for escape vs relaxation trajectories, which, instead, is a signature of non-equilibrium [37]. This can be explained by considering that, as discussed in [16], the ocean model evolves approximately in an energy landscape.
When µ = 0.98, the SB → W transitions are rather rare unless one considers relatively large values of σ. This is due to the much lower value of Φ at the SB attractor than at the W attractor (see Eq. 3), so that (see Eq. 2) the fraction of the invariant measure supported near the W attractor is extremely small. We focus next on the case of µ = 1.0, where the population is split more evenly between the W and SB attractors. We are then able to construct for each value of σ τ the invariant measure using a single orbit of the system, provided we can observe a sufficient number of transitions. We use σ 100y = 1.5%. Our results are shown in Fig. 4 for a trajectory lasting ≈ 6.0 × 10 4 y and characterised by 92 SB → W and W → SB transitions, whose average rates are consistent with an occupation of about 35% for the W basin of attraction, and of about 65% for the SB basin of attraction. The projection of the invariant measure on the ([T S ], ∆T S ) plane shows that the peaks of the pdfs are very close to the W and SB attractors (note the predicted slight shift for the W case, visible because the noise is stronger than in Fig. 3), and that the agreement further improves when considering the two marginal pdfs (top left and bottom right insets). We can construct both the W → SB and the SB → W instantons, whose starting and final points agree remarkably well with the attractors and the M state. We discover that the instantons follow a path of monotonic descent, following closely the crests of the pdf), with the minimum at the M state. We also run a simulation lasting ≈ 2.7 × 10 4 y using σ 100y = 1.8%. We obtain 73 SB → W and W → SB transitions. Comparing the statistics of this run with what is shown in Fig. 4, we find good agreement with the predictions of Eqs. 2-3 (not shown).
Concluding, in this letter we have studied the problem of noise-induced transitions between the W and the SB state of the climate system using a intermediate complexity stochastic model. The deterministic version of this model had been used to construct the M states for a vast range of values of the solar irradiance [9]. Including stochastic perturbations allows for exploring the phase space of the system. We have cho-sen to study the impact of fluctuations in the solar irradiance, which entails adding a multiplicative noise. In particular, the SB climate reflects more radiation because large ice-covered surfaces lead to higher albedo, so that the effective noise will be weaker than the one acting near the W attractor. We have explained why a large deviation theory-based mathematical framework is able to explain very satisfactorily our results. We show that one can construct the pseudo-potential defining the escape rate from the basins of attraction and the natural measure of the stochastically perturbed system. In future investigations we plan to extend the analysis performed here to values of µ covering the whole range of bistability, in order to understanding how the pseudo-potential depends on µ. Using a low-dimensional projection of the invariant measure, we show that the instantons connect the attractors with the M states, following a path of steepest descent. This gives a key connection between the stochastic and the deterministic points of view on study of the multistability of the climate.
A further link between the stochastically perturbed and deterministic system can be described as follows. If we consider values of µ just below the critical one µ W →SB ∼ 0.965 defining the W → SB tipping point, one can find long-lived transient chaotic trajectories. It is worth investigating whether such transients are associated to a invariant saddle emerging as a result of the boundary crisis. We have observed that the way such long-lived trajectories collapse to the SB state is very similar to the way, when µ = 0.98, stochastically perturbed orbits initialised in the W state basin of attraction perform the transition. In future studies we will explore such similarity looking into the possible presence of a ghost state [38].
The knowledge of M states is key to understanding tipping points. When an attractor and a M states are nearby, the system's response to perturbations diverges, as the correlations decay very slowly [19]. These findings extend and generalise the classical point of view given in [11]. What is proposed here, combined with the framework given in [39], could be key for understanding and predicting the criticalities in the trajectories of the Earth system [40], including those leading to very hot, non-habitable conditions [41]. A tipping element we will investigate along the lines of the present study is the Atlantic meridional overturning circulation [11].