Basin stability and limit cycles in a conceptual model for climate tipping cascades
Nico Wunderling, Maximilian Gelbrecht, Ricarda Winkelmann, Jürgen Kurths, Jonathan Donges
BBasin stability and limit cycles in a conceptual modelfor climate tipping cascades
Nico Wunderling, , , * Maximilian Gelbrecht, , Ricarda Winkelmann , ,J ¨urgen Kurths , , & Jonathan F. Donges , * Earth System Analysis, Potsdam Institute for Climate Impact Research (PIK),Member of the Leibniz Association, 14473 Potsdam, Germany Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam, Germany Department of Physics, Humboldt University of Berlin, 12489 Berlin, Germany Complexity Science, Potsdam Institute for Climate Impact Research (PIK),Member of the Leibniz Association, 14473 Potsdam, Germany Lobachevsky State University Nizhny Novgorod, Nizhny Novgorod, Russia, Stockholm Resilience Centre, Stockholm University, Stockholm, SE-10691, Sweden * Correspondences should be addressed to [email protected] or [email protected] a r X i v : . [ phy s i c s . a o - ph ] N ov ipping elements in the climate system are large-scale subregions of the Earththat might possess threshold behavior under global warming with large po-tential impacts on human societies. Here, we study a subset of five tipping ele-ments and their interactions in a conceptual and easily extendable framework:the Greenland and West Antarctic Ice Sheets, the Atlantic Meridional Over-turning Circulation (AMOC), the El-Ni ˜no Southern Oscillation (ENSO) andthe Amazon rainforest. In this nonlinear and multistable system, we performa basin stability analysis to detect its stable states and their associated Earthsystem resilience. By combining these two methodologies with a large-scaleMonte Carlo approach, we are able to propagate the many uncertainties asso-ciated with the critical temperature thresholds and the interaction strengths ofthe tipping elements. Using this approach, we perform a system-wide and com-prehensive robustness analysis with more than 3.5 billion ensemble members.Further, we investigate dynamic regimes where some of the states lose stabil-ity and oscillations appear using a newly developed basin bifurcation analysismethodology. Our results reveal that the state of four or five tipped elementshas the largest basin volume for large levels of global warming beyond 4 ˝ Cabove pre-industrial climate conditions, representing a highly undesired statewhere a majority of the tipping elements reside in the transitioned regime.For lower levels of warming, states including disintegrated ice sheets on WestAntarctica and Greenland have higher basin volume than other state config-urations. Therefore in our model, we find that the large ice sheets are of par-ticular importance for Earth system resilience. We also detect the emergenceof limit cycles for 0.6% of all ensemble members at rare parameter combina-tions. Such limit cycle oscillations mainly occur between the Greenland Ice heet and AMOC (86%), due to their negative feedback coupling. These limitcycles point to possibly dangerous internal modes of variability in the climatesystem that could have played a role in paleoclimatic dynamics such as thoseunfolding during the Pleistocene ice age cycles. Introduction
During the last decades, the field of tipping elements has become a major point of interest incomplex systems and network science [1, 2]. They have been used in the description of variousfields such as in financial markets, technological progress, ecology or in climate science [e.g.,3, 4, 5, 6]. Tipping elements can interact across scales in space and time [7] which couldpotentially lead to catastrophic domino effects [8] or, for instance, lead to a hothouse cliamtestate in the case of climate tipping elements [9].In the climate system, tipping elements are subregions of the Earth system that can exhibitthreshold behavior, where a small forcing perturbation can be sufficient to invoke a strong non-linear response of the system that can qualitatively change the state of the whole region orsystem due to internal, self-enforcing feedbacks [6]. Climate tipping elements comprise sys-tems from the cryosphere (e.g. Greenland, Antarctic Ice Sheet, Permafrost), the biosphere (e.g.Amazon rainforest, coral reefs) and large-scale circulation systems (e.g. Monsoon systems, At-lantic Meridional Overturning Circulation) [6, 10]. Their potential tipping to alternative stateswould be associated with severe impacts on the biosphere and threaten human societies [11].It has been suggested that several climate tipping elements are at risk or on the way of trans-gressing into an undesired state even at global warming levels below the 2.0 ˝ C goal of theParis Agreement [10, 11, 12]. Among others, tipping elements that already show warning sig-nals of degradation at present times [11, 12] are: the West Antarctic Ice Sheet where partsin the Amundsen Bay (Pine Island & Thwaites region) are suspected to have been destabi-lized [13, 14, 15], the AMOC which experienced a major slowdown of 15% from 1950 tonow [16], the Amazon rainforest which might approach a tipping point due to climate changeand deforestation [17]. Critical deforestation ratios might lie between 20 to 40%, where cur-4ent deforestation is reaching 20% [17, 18]. Furthermore, the Greenland Ice Sheet loses massat an accelerating pace [19, 20] and the frequency of major El-Ni˜no events are suggested toincrease twofold and strong ENSO effects will occur more often as global warming contin-ues [21, 22]. However, others highlight that large uncertainties are related to future changes ofENSO and whether major El-Ni˜no events will become more frequent or intense under globalwarming [23, 24].Furthermore, contradicting a common misunderstanding, tipping elements do not necessarilytip immediately after the crossing of their tipping point, but their tipping time trajectory mighttake very long and appear smooth [25]. For instance for the large ice sheets, the disintegrationtime scale could be on the order of several centuries up to millennia as has been suggested bymodeling studies [26, 27, 28].For most of the tipping elements, there is a critical temperature range at which they are suspectedto leave their current safe state separating the climate tipping elements into three groups [10].The first group comprises elements that might transgress their state within the limits of the ParisAgreement (Paris, 2015) of 2 ˝ C above pre-industrial and with that, these are the most vulner-able climate tipping elements with respect to global warming. This group contains mainlycryosphere elements (Arctic summer sea ice, West Antarctic Ice Sheet, Greenland Ice Sheetand Alpine glaciers) as well as the Coral reefs that are likely to be lost even when global warm-ing is restricted to 2 ˝ C above pre-industrial. The second group might tip at temperatures above3 ˝ C (for instance the Amazon rainforest, AMOC or ENSO) and the most resilient group onlyat temperatures around 5 ˝ C above pre-industrial or higher (e.g., parts of the Antarctic ice sheet,permafrost or Arctic winter sea ice).However, the tipping elements in the climate system are not independent of each other, butconnected [29, 11] and the knowledge about the exact interaction structure is sparse and par-5ially based on experts that, for instance, suggested an interaction structure, including sign andstrength, for a subset of five tipping elements: The Greenland Ice Sheet, the West AntarcticIce Sheet, the AMOC, the Amazon rainforest and the ENSO [29]. Behind each connection be-tween two tipping elements within this subset, there is a physical process or set of processes (seeTab. A.1). For instance the impact of the Greenland Ice Sheet on the AMOC due to freshwaterinput from melting ice slows down the AMOC on the one hand and a weakening AMOC on theother hand cools latitudes in the northern hemisphere. Note that this subset network of tipping isneither complete in the number and selection of tipping elements, nor is it comprehensive in thepossible connection pathways and their potential strength between the tipping elements. Thereare also earlier investigations on tipping points [30] and the interaction of tipping points [31, 32]in the context of economic damage and the social cost of carbon using further developed ver-sions of the integrated assessment model DICE [33, 34]. Here, Cai et al. (2016) [31] explicitlybase their findings on the interactions of tipping elements from the expert elicitation in Kriegleret al. (2009) [29].Following the elaborations above, in this work we aim at investigating the resilience of variousattractors for interacting climate tipping elements and we want to elucidate the role that differenttipping cascades have in that regard. The approach put forward here can easily be adapted tomore tipping elements and further interaction structures once they are more comprehensivelyunderstood [see also 35].We explore the stability landscape and the dynamics of a subset of five tipping elements rep-resented by normal form fold bifurcations based on known interactions across scales in timeand space between these tipping elements (Fig. 1) [29]. These five tipping elements are theGreenland and West Antarctic Ice Sheets, the AMOC, the ENSO and the Amazon rainfor-est [29]. We introduce the model of interacting tipping elements in section. We use the con-6ept of basin stability [36] in order to determine the basin sizes of various attractors of thismultistable system (see section 2). Based on a large-scale Monte Carlo ensemble, this method-ology gives an estimate how stable and resilient various attractors are. It has been appliedto many dynamic systems before such as power grids, neuronal models and further nonlinearsystems [37, 38, 39, 40, 41]. Furthermore, especially for larger coupling strengths, Hopf bifur-cations can occur, thus invoking oscillatory limit-cycle solutions of the model. For the detectionand quantification of these types of limit cycle attractors, we apply a newly developed bifurca-tion algorithm that is able to identify different dynamical properties in complex systems: theMonte Carlo Basin Bifurcation analysis (MCBB) [42].In the following, we first introduce the methodological approach of this work: the model ofinteracting tipping elements (section 2.2), the basin stability approach (section 2.3), the con-struction of the large-scale Monte Carlo ensemble (section 2.4) and the Monte Carlo Basin Bi-furcation analysis (section 2.5). Then, we evaluate the basin volume in our model (section 3.1)and quantify the occurrence of limit cycles in our model (section 3.2). Lastly, the results withrespect to the climate system are discussed and summarized in sections 4 and 5.
For the purpose of investigating the dynamical properties of the subset of five tipping elements,we further developed a conceptual network approach that is fully dynamic and captures themain nonlinear dynamical properties of tipping elements [35, 44, 45]. The actual physicalprocesses behind the tipping elements are not explicitly modeled to maintain an accessible andcontrollable structure. The modeling of complex systems using conceptual approaches is apopular tool and has been successfully applied to, among others, ecology, social systems orepidemiology [5, 46, 47]. 7 reenland Ice Sheet (disintegration)
West Antarctic Ice Sheet (disintegration)
Amazonrainforest (dieback)
El-Niño Southern Oscillation (permanentness) ++ + ? + + GMT2.0°C GMTGMTGMTGMT4.8°C3.2°C5.3°C4.0°C
Greenland Ice SheetAMOCWest Antarctic IceSh.ENSOAmazon rainforest S t a t e o f t i pp i n g e l e m e n t A B
Figure 1.
A) World map with connections shown for five tipping elements where the inter-action structure is known from an expert elicitation [29, 43]. Each link represents a physicalmechanism and has a certain strength (see supplementary Tab. A.1). A positive arrow impliesan effect that drives the tipping element closer to its tipping point, a negative arrow drives thetipping element away from its tipping point and a question mark stands for an unclear direction.B) (Fold) Bifurcation diagram of each of the tipping elements without coupling. On the x-axis,the average global warming required for tipping is shown in degrees above pre-industrial globalmean temperature levels for the respective tipping element.8 .1 Tipping elements and interactions
Given that, despite major advances, current EMICs (Earth system models of intermediate com-plexity) and GCMs (global circulation models) are not yet able to fully represent the nonlin-ear behavior of some Earth system components together with their interactions, but physicsbased models and equations as well as paleo climate observations suggest the existence of suchproperties for many tipping elements as for instance the Greenland and (West) Antarctic IceSheet [26, 28, 48, 49], the AMOC [50, 51, 52, 53, 54], the Amazon rainforest [55, 56, 57, 58, 59]or the ENSO [60, 61], the conceptual approach chosen here demonstrates an option how tomodel interactions between tipping elements. Thus, we put this model forward as a first steptowards a more process-detailed assessment of tipping elements and their interactions. Thisalso emphasizes that future research could focus on developing more complex, emulator- orEMIC/GCM-like models of tipping elements to investigate their nonlinear interplay such as hasrecently been developed for the Antarctic Ice Sheet [49].While we have described why we use a conceptual description of the main dynamics of thetipping elements directly above, we outline the physical mechanisms of there interactions here-after. Although some interactions between the tipping elements are better understood and evalu-ated than others, we describe the physical mechanisms of all of them in the following separatedinto destabilizing ( ` ), stabilizing ( ´ ) and unclear links (?) (see Fig. 1). This description aimsto provide a basic physical understanding, but cannot resolve the problem of how strong exactlyeach of these interactions is. Therefore, the interaction structure is kept as described later inEq. 5 with the multiplicative interaction strength factor d .1. Destabilizing interactions:I) Greenland Ice Sheet Ñ AMOC: When the Greenland Ice Sheet starts to melt, it has a9iminishing influence on the overturning strength of the AMOC due to freshwater inputinto the North Atlantic. This has been observed in modeling studies [16, 51, 62, 63] andin observations [64].II) AMOC Ñ West Antarctic Ice Sheet: When the AMOC collapses, sea surface temper-ature anomalies arise due to the collapse of the northward heat transport of the AMOC.This results in a cold north and a warm south of the equator as shown by modeling stud-ies [65, 66, 67, 68].III) Greenland Ice Sheet Ñ West Antarctic Ice Sheet (and vice versa): The shift of ground-ing lines due to changing sea level is a well-known phenomenon from tidal changes [e.g.69]. Thus, if the sea level rises due to global warming, the floating ice shelves couldbe lifted which is likely to result in grounding line retreat. Furthermore, gravitationalchanges as well as elastic and rotational effects might then amplify the sea level changeif one of the large ice sheets disintegrates first because the gravitational attraction thenonly emanates from the other, remaining ice sheet [70, 71]. The effect would be strongerif Greenland melts first since the West Antarctic Ice Sheet has more marine terminatingglaciers and ice shelves.IV) AMOC Ñ ENSO: There are two opposing effects that have proposed, which describehow the AMOC might influence the ENSO: (i) It has been suggested that oceanic Kelvinwaves originate from a colder North Atlantic and travel southward. Then, in westernAfrica, Rossby waves would be emitted towards the north and the south, which are thentranslated back into Kelvin waves that travel into the Pacific ocean. This effect woulddeepen the Pacific thermocline and weaken the amplitude of ENSO [60]. (ii) With aweaker AMOC, the northern tropical Atlantic would in turn become cooler and northerlytrade winds would be intensified over the northeastern tropical Pacific. It has been ar-gued that this could result in a southward displacement of the Pacific ITCZ leading to a10ea surface temperature anomaly [72]. At the same time, it is argued that Rossby wavesare sent into the northeast tropical Pacific [73]. This would intesify ENSO due to windstress interaction from AMOC. Overall, it is believed that mechanism (ii) is stronger thanmechanism (i) [60]. Furthermore, it can be extracted from complex Earth system modelsthat a decrease in AMOC intensity indeed strengthens the variability of ENSO [61, 74].V) ENSO Ñ Amazon rainforest: Literature studies suggest that droughts related to cli-mate variabilities such as the El-Ni ˜no Southern Oscillation can affect the stability of theAmazon rainforest [75, 76, 77]. Using an EMIC, it has been found that a permanentEl-Ni˜no state would endanger substantial portions of the Amazon basin due to a reor-ganization and reduction of water access in the South American tropics via teleconnec-tions [78].VI) ENSO Ñ West Antarctic Ice Sheet: The interaction between ENSO and the WestAntarctic Ice Sheet is one of the least certain interactions as has already been stated inKriegler et al (2009) [29]. Nevertheless, there are hints for warming oceanic effects fromEl-Ni˜no in the Amundsen and Ross Sea region, while La Ni˜na would cool this region.At the same time, atmospheric effects could have an opposite effect, which would offsetthe oceanic effect [79]. Besides that, it has been found with satellite observations that iceshelves gain height, but yet lose mass during El-Ni˜no events in the Amundsen and RossSea region [80]. While the primary driver of melt in West Antarctica is the warm oceanwater below ice shelves, an extended period of surface melting has been observed duringJanuary 2016, which is likely promoted by the strong El-Ni ˜no event in this year [81].Since it is expected that the frequency of major El-Ni ˜no events will increase during cli-mate change [21], we set this interaction positive (see Tab. A.1 and Fig. 1).2. Stabilizing interactions: 11) AMOC Ñ Greenland Ice Sheet: For a decreasing overturning strength of the AMOC,the northern hemisphere is cooled since the heat transport towards the North Atlanticwould be weakened. This has been observed in modeling studies [16, 65, 66, 82].II) ENSO Ñ AMOC: Using reanalysis data, evidence has been found that the transport ofwater vapor out of the tropical Atlantic is enhanced [83]. Comparing La Ni ˜na and El-Ni˜noconditions, it was found in this study that El-Ni˜no conditions lead to a stronger northernAMOC on a multi-decadal timescale. However, another study questions this finding anddoes not find a strong impact on the deepwater formation from AMOC [84]. Therefore,this interaction is less well established from literature and therefore considered of lowstrength, but with a negative sign (see Fig. 1 and Tab. A.1).3. Unclear interaction direction:I) AMOC Ñ Amazon rainforest: When the AMOC shuts down, the intertropical conver-gence zone (ITCZ) is likely dislocated southward, leading to large changes in seasonalprecipitation on a local to regional degree. This might then impact parts of the Amazonrainforest [68, 82, 85]. Still, it is unknown as to whether this interaction is positive ornegative and might differ from region to region. Therefore, this link is set as unclear (seeFig. 1).II) West Antarctic Ice Sheet Ñ AMOC: A literature study using a coupled ocean-atmospheremodel found a decrease in the AMOC for high freshwater inputs from the West AntarcticIce Sheet [86]. However, another study detected a stabilization of the AMOC if influencedby freshwater input from West Antarctica. This is ascribed to the effects from the bipolarocean seesaw due to decreasing Antarctic Bottom Water formation [87]. With an EMIC,is has been found from using freshwater input experiments into the Southern Ocean thatdifferent processes could enhance or slow down the AMOC [88]: (i) The deep water ad-12ustments via the bipolar ocean seesaw tend to intensify the NADW formation. (ii) TheNADW is strengthened by southern hemispheric wind increase representing an ocean-atmosphere interaction. (iii) Salinity anomalies from the Southern Ocean are distributedto the North Atlantic weakening the NADW [compare to 86]. Overall, the processes (i)and (iii) strengthen the AMOC and process (ii) weakens it. However, the exact time scaleand efficiencies of these processes have been rated unknown as of yet [88].III) Amazon rainforest Ñ ENSO: Under a dieback of the Amazon rainforest, the moisturesupply to the atmosphere will significantly change, also since the atmospheric moisturerecycling feedback over the Amazon basin would break down [89, 90, 91]. However, it isunclear whether and to which extent this would then impact ENSO.
In our conceptual model, we divide the dynamics x i of the considered tipping elements i intotheir individual dynamics f i p x i q and a direct interaction term g i p (cid:126)x q ” g i p x q . This yields τ i x i “ f i p x i q ` g i p x q , (1)where τ i is the typical time that passes when a tipping element undergoes a critical transitionfrom one state to another. We model the individual dynamics of each of the tipping elementswith the general tipping approach (CUSP equation [8, 92]) f i p x i q “ ´ a i x i ` b i x i ` c i a i , b i , c i P R , (2)where a i ą and b i ą . Assuming additive separability of the interactions between the tippingelements and linear interactions, the interaction term g i p x q becomes g i p x q “ ÿ j g ij p x i , x j q linear “ interactions ÿ j A ij x j . (3)13ere, A ij is the interaction structure and strength, which is set to zero if there is no connectionbetween the tipping elements i and j . Altogether, Eq. 1 becomes τ i x i “ ´ a i x i ` b i x i ` c i ` ÿ j A ij x j . (4)Each tipping element x i following this equation possesses two fold bifurcations at ˘ a p a i q{p b i q and has already been investigated in theoretical works on tipping cascades [92], but also in var-ious contexts where nonlinear behavior is important as for instance in policy, environmentalissues, economy or climate [8, 93]. For these equations exist a framework that allows to investi-gate tipping cascades on larger networks with regard to their interaction structure in the networkas well as microstructures that are decisive for finding emergent tipping cascades [44, 45].In our model, we specify the interaction structure and strength term A ij by setting it equal to amultiplicative factor d times the actual link strength s ij between each pair of tipping elements.Therefore, A ij “ d { ¨ s ij . The link strength values s ij are taken from the expert elicitation [29].The factor { is used for normalization reasons since then d P r , s . If we now additionallyset a “ , b “ and c i “ ´a { { T limit, i ¯ ¨ ∆ GMT, the tipping elements are described by thefollowing nonlinear, ordinary differential equation (all parameters of Eq. 5 are explained in theTabs. A.1 and A.2) dx i dt “ »—– ´ x i ` x i ` a { T limit, i ¨ ∆ GMT ` d ¨ ÿ jj ‰ i s ij p x j ` q fiffifl τ i . (5)Here, x i is the state of the tipping element (see Fig. 1B) and i stands for the considered tippingelements i “ Greenland Ice Sheet , West Antarctic Ice Sheet , AMOC , ENSO , Amazon rainforest.We choose these five tipping elements since their interaction structure is known from an expertelicitation [29]. The increase of the global mean temperature above pre-industrial is denoted by ∆ GMT, T limit, i is the critical temperature threshold of the respective element. The last term is14he coupling term, where d is a general multiplicator that determines the strength of the inter-action term in comparison to the other, individual dynamics terms. The parameter d is variedbetween 0, meaning no interactions, and 1, where the interactions become as important as theindividual dynamics. Following this, one might tend to assume that the individual dynamicsof the tipping element influences the tipping element more than the interaction effect. Thismight make smaller coupling parameters more realistic than higher ones. In Eq. 5, s ij is thelink strength that is based on the expert elicitation [29] and τ i is a typical timescale at which acertain tipping element transgresses its state.This typical tipping time scale ranges from decades for the Amazon rainforest to several millen-nia for the large ice sheets (see Appendix Tab. A.2). Then, our system of differential equationsis integrated forward in time using scipy.odeint [94] until more than 20 times the Greenland IceSheet’s typical transition time scale has passed. This is equal to 100,000 years simulation time.This is the time when equilibrium is reached in the simulations. However, we are not intendingto compute an exact time scale for tipping or tipping cascades here, but we are rather interestedin the system’s attractors and their stability properties. This is why we denote model years inarbitrary units instead of giving an exact time, also since this would be beyond the scope of thisconceptual model (see Fig. 1A). Note that we adapted the link from ENSO to AMOC from un-certain to negative compared to the original results of the expert elicitation on tipping elementinteractions [29] since there is only a dampening process known in literature [43].There are considerable uncertainties associated with this approach, especially with the criticaltemperature at which a certain tipping element transgresses its state T limit, i as well as in thestrength of the interactions s ij . The uncertainties of these two parameters are shown in theAppendix Tabs. A.1 and A.2. Thus, with Eq. 5, we model tipping events and cascades undercertain conditions of global warming (GMT) and the interaction strength ( d ).15 .3 Basin stability We are interested in the stability properties of different attractors within the state space. Anappropriate tool to investigate the stability landscape of such states is the so-called basin sta-bility [36, 39]. Basin stability is a nonlinear stability measure for the resilience of an attrac-tor to disturbances. Where traditional measures such as the computation of Lyapunov expo-nents or Master stability functions rely on linear approximations in reaction to small perturba-tions [95, 96], basin stability approaches can also consider large perturbations. Such perturba-tions can occur in Earth system components such as the large ocean circulations or the Amazonrainforest [36, 97]. The basin stability is an established algorithm focusing on the stabilitylandscape of the entire phase space, while other nonlinear stability measures such as surviv-ability [41], stability threshold [98], constrained basin stability [99] and topology of sustain-able management [100] approaches focus on the stability of parts of the state space or desiredregimes in it. Therefore, basin stability computations are a first step that aims to quantify thestability of different attracting states, but do not aim to study potential desired regimes as wouldbe required for the other mentioned methods. The concept of basin stability has been appliedto many multistable systems. Examples comprise the Amazon rainforest [36], the stability innetworks of power grids [37], neuronal models [38] and further nonlinear systems such as incoupled network systems [39], oscillators [101, 102] or chimera states [103]. While basin sta-bility can widely be applied, it has its limitations, for instance in cases where basins becometoo peculiar, e.g. for riddled basins with holes [40]. Since this is not the case in our model ofinteracting climate tipping elements, we utilize basin stability in this work.An attractor A is defined as the minimal compact invariant set A Ď X , where X is the entirestate space [104]. B p A q Ď X is the basin of attraction of A which comprises all states fromwhich the system converges to A . The basin stability or the basin volume V p A q is then quan-tified as the probability that a system will return to a certain attractor A after a perturbation16 B p A q : “ ż R B p A q dµ P r , s , (6)where B p A q is 1 in case x P B p A q and 0 otherwise. µ is a measure on the state space X that encodes the relevance of a certain perturbation and our knowledge about the system. Theestimation of the integral in Eq. 6 can be difficult, but in our system it can be assumed thatthe estimation of the basin volume can be estimated via a Monte Carlo ensemble. The totalvolume of a basin of attraction is then measured as the fraction of simulations with randomlychosen initial conditions that end up in that certain attractor N p A q over the total number ofinitial conditions N p Ω q V p A q “ P p A q “ N p A q{ N p Ω q . (7)Here, P p A q is the probability that a random initial condition ends up in the basin of attractor A . To assign the basin volume V p A q with the probability P p A q , it is required that the space ofinitial conditions is covered well and uniformly. Therefore in this work, it is necessary to extendthe classical concept of basin stability since it is not only required to sample the space of initialconditions sufficiently well, but also to sample over the uncertainties in the model parametersthemselves (see Tabs. A.1 and A.2). Thus, we need to set up a very large-scale Monte Carloensemble of several billion ensemble members whose construction details can be directly foundbelow. In order to apply the concept of basin stability in a meaningful way, the state space must be cov-ered well enough. However, in this application, the parameters of the models have uncertaintiesthemselves in the critical temperature thresholds and the interaction strength and structure. Thismeans, we need a way of covering the many uncertainties in these various parameters as wellas the state space itself. Therefore, it is necessary and useful to combine basin stability with a17arge scale Monte Carlo sample that covers an adequate extent of the phase space and parameterspace. This is what explain in the hereafter.The basic Monte Carlo ensemble without the extension for basin stability is set up as follows:for each pair of global mean temperatures (GMT) and interaction strengths d, there is a sam-ple of size 100 constructed with initial conditions from the uncertainty range in T limit, i and s ij using a latin hypercube algorithm [105] (see Tabs. A.1 and A.2). Latin hypercube samplingis an extension to the usual random sampling and is used to improve the space coverage ofinitial conditions. Therefore, the space of initial conditions is separated into its dimensions,i.e., the number of different initial parameters (here 17, see Tabs. A.1 and A.2). Then, it issecured that only one sample occurs in each axis hyperplane (compare to the N-rooks problemin mathematics). We apply this sampling procedure for each of the 27 different network setupsthat arise from the permutation (positive, negative, zero) of the three uncertain links (Amazonrainforest Ñ ENSO, AMOC Ñ Amazon and West Antarctic Ice Sheet Ñ AMOC, see Fig. 1A).This then leads to 2700 samples. These 2700 samples are computed for each global mean tem-perature increase up to 8 ˝ C above pre-industrial which can be reached in business as usualscenarios RCP8.5 extended from 2100 to 2500 [10] in steps of 0.1 ˝ C and coupling constant dbetween 0.0 and 1.0 in steps of 0.02 accounting for 864.000 simulation runs.The extension of the Monte Carlo ensemble, integrating basin stability is detailed below: thebasin stability of the system for each of these 864.000 samples is computed by permuting theinitial state of each of the five tipping elements within its limit, i.e., between the untipped(x= ´ . ) and the tipped state (x= ` . ). The state variables of the five tipping elements resultin a five dimensional state vector 18ec basin stability, i “ t x Greenland Ice Sheet p q , x AMOC p q , x West Antarctic Ice Sheet p q ,x ENSO p q , x Amazon rainforest p qu , where vec basin stability, i P r´ . ` . s for each tipping element. However, vec basin stability, i cannotbe permuted in a completely random way, but each of its five dimensions needs to be permutedin an independent way since there is a strong nonlinearity at state equal zero for each of thefive dimensions. Of course, in principle if there would be infinite computational resources, wewould not need to take this nonlinearity into account, but would be able to increase the size ofthe Monte Carlo ensemble even further. But since this is not the case, we need to “manually”account for this important nonlinear property. This means that the sign of each state must beequally probable, i.e.: P sign “ P p´ , ´ , ´ , ´ , ´q “ P p` , ´ , ´ , ´ , ´q “ P p´ , ` , ´ , ´ , ´q ““ . . . “ P p` , ` , ` , ´ , `q “ P p` , ` , ` , ` , ´q “ P p` , ` , ` , ` , `q ““ { “ . . (8)This can be achieved when random starting conditions are drawn from each of the 32 combina-tions of P sign . Hence, for each of the 32 combinations, we chose 10 different initial conditionsending up with 320 different settings. For the 320 randomly chosen perturbations (i.e., the ini-tial conditions of the tipping elements), we again used a latin hypercube algorithm [105]. Thatmeans it fulfils the condition that each of the 32 different possible signs of the initial conditionsin their five-dimensional subspace (one dimension for each tipping element) is covered equallyoften.Altogether, we employ a very large ensemble of simulations to compute the basin stability of19 total “ ¨ . “ . . . « . ¨ samples. How the final state can depend onthe initial conditions is shown exemplary for three timelines in Fig. 2A-C. s t a t e [ a . u . ] Tipped regimeTransition regimeUntipped regime
Model time [a.u.]
Greenland West Antarctica AMOC ENSO Amazon rainforest
A B C
Figure 2.
A-C) Timelines at GMT=2.5 ˝ C above pre-industrial levels and d “ . withdifferent initial conditions (ICs) as used to probe basin stability. A) IC “ r´ , ´ , ´ , ´ , ´ s (today’s initial conditions, no element is tipped at t “ ), B) IC “ r´ . , ´ . , . , ´ . , . s ,C) IC “ r . , . , ´ . , . , . s . Coupling nonlinear ODEs as in the model described here, invokes the possibility of furthertypes of bifurcations besides fold bifurcations. Here, we utilize Monte Carlo Basin BifurcationAnalysis [42] to uncover system attractors and estimate their basins of attraction finding Hopf-Bifurcations and thus oscillating solutions converging to limit cycle attractors. MCBB is anovel, numerical approach to analyze multistable systems, quantify and track their asymptoticstates in terms of their basins of attraction by utilizing random sampling and clustering methods.Since MCBB is based on Monte Carlo ensembles and we are interested in a quantitative measureof interesting dynamical properties (here occurring limit cycles, i.e., Hopf-bifurcations), it is awell suited method for our purposes. It has also been applied to other nonlinear systems suchas the Dodds-Watts model, the Kuramoto model or Stuart-Landau oscillators [42].MCBB aims to find classes of attractors that collectively share the largest basins of attractionsof the system. Similar attractors, at different parameter values, have to share similar valuesof invariant measures ρ and the difference of theses measures has to smoothly vanish if the20arameter difference goes to zero. If this is the case they are regarded as being part of thesame class of attractors. N tr trajectories of the system, here 140 000, with randomized initialconditions and parameters are integrated. In order to identify the different classes of attractors,suitable statistics S i are measured on every system dimension for every trajectory, here, themean, variance and the Kullbach-Leibler divergence to a normal distribution. Hence, for everystatistic i , S i is a N tr ˆ matrix. A distance matrix of each trajectory to each other is computedfrom these statics with D ij “ ÿ k w k ÿ l | S k,il ´ S k,jl | ` w | p p i q ´ p p j q | (9)where p p i q is the control parameter used to generate the i -th trajectory and w i are free parametersof the method, here r
1; 0 .
5; 0 .
5; 1 s which is the default recommendation for these parameters.This distance matrix is used as an input for a density-based clustering algorithm such as DB-SCAN which can find if this notion of continuity between different trajectories exists and thuseach cluster corresponds to a different class of attractors. For further details on MCBB, referto [42]. When applying this to the conceptual model for climate tipping points, not only thedifferent possible states of tipped elements are found, but also different classes of oscillatingstates induced by Hopf-Bifurcations are found. For the MCBB analysis, the parameter uncer-tainties were varied randomly within the same bounds as for the previously described basincomputations. The initial conditions of five tipping elements were chosen to all start at ´ ,i.e., not tipped for the results presented in the main text, and at random between ´ and ` forthe results presented in the appendix. The computations are performed with the Julia libraryMCBB.jl. 21 Results
We compute the basin stability of each potential state that could be governed by the network offive tipping elements.. The present day state could be considered as some kind of safe state forthe Earth system when all five tipping elements are in a negative state. On the other hand therecould be a state where all five tipping elements reside in the positive, tipped state. In betweenthere are intermediate scenarios, where some tipping elements already crossed their thresholdsand others did not. In Fig. 3, we show the average basin stability for each of these six possiblesituations, i.e., with zero, one, two, three, four and five tipped elements. In this experiment, weperturbed the initial conditions of all tipping elements at the same time. The fraction of initialconditions that end up in the respective basin are plotted as the color.In general, we observe that the size of the basin of attraction for higher global warming levelsbecomes larger for a higher number of tipped elements as would be expected. For high levelsof warming, the basin of five tipped elements dominates.For increasing interaction strength, the volume of the basins with three or less tipped elementsdecreases (Fig. 3A-D). Contrasting this, the basin volume with four tipped elements increaseswith increasing interaction strength, while the basin for five tipped elements first increases andthen decreases again (Fig. 3E, F). the last issue is due to the strong negative feedback loopbetween the Greenland Ice Sheet and the AMOC. In such cases of high coupling, the AMOCtips, but safeguards the Greenland Ice Sheet which reaches the untipped regime for global meantemperature increases above 4 ˝ C and interaction strengths above 0.5. This poses a hypothet-ical scenario which would only be realistic if the interaction strength between Greenland andAMOC is very high, but this behavior has also been observed in experiments of tipping cascadesearlier [35]. 22
Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e A B C DE F
Figure 3.
Basin volume for each number of tipped elements in dependence of interactionstrength and global mean temperature increase. A-F) Basin volume of 0 to 5 tipped elements.23or instance, in the basin volume plot of zero, one or two tipped elements, the number of statesthat equilibrate in this state is very small (Fig. 3A, B and C). For temperature increases above2 ˝ C the associated basin volume is close to zero for all interaction strengths. At the same time,the size of the basin decreases for higher coupling strengths.The uncertainties of the basin volumes are quantified as standard deviation in the appendix(Fig. B.1). We find that uncertainties generally increase for a higher amount of tipped elementsas well as for higher interaction strengths. The standard deviation is highest for small tem-perature increases and high coupling strengths since here, the attractors depend on the initialconditions in terms of the critical temperature thresholds and initial coupling constants (seeTab. A.1). The basin of four and five tipped elements show a regime of increased standard de-viation for temperatures around 2-5 ˝ C above pre-industrial and interaction strength parametersof more than 0.2. This is probably due to the fact that in this regime the state of the GreenlandIce Sheet has a large variation because of its strong negative feedback loop to AMOC. Thus,whether this element tips, also depends a lot on the explicit initial conditions of the state as wellas on parameters (see Tabs. A.1 and A.2). Outside and around this regime, the uncertainty issmaller since either Greenland is not tipped with high certainty for lower temperature increases(below 2 ˝ C) or tipped with high certainty at higher temperature increases (above 5 ˝ C).There exists a narrow range of global mean temperature increases when single tipping elementscan transgress their state without triggering a tipping cascade. This range is mostly locatedbelow 1 ˝ C above pre-industrial for low coupling strength and well below 1 ˝ C for higher in-teraction strengths (see Fig. 3B and Fig. B.2). If we separate this response into the respectivesingular tipping elements, we can see that above an interaction strength of 0.2-0.4, the Green-land Ice Sheet and ENSO cannot tip without causing a cascade due to their strong interactionslinks to AMOC or the Amazon rainforest, respectively (for more details see Appendix B).24dditionally, we investigate some important intermediate states in more detail, where some el-ements are in the tipped regime, while others are not. It was found that several tipping cascadesof size two and three are more frequent than others, for instance a tipping cascade betweenthe Greenland and the West Antarctic Ice Sheet is more likely than, for instance, a cascade be-tween the AMOC and the Amazon rainforest [35]. Thus, we investigate the basin volume thatcorresponds to such cascades.We find that the ice sheets appear to be of particular importance for the stability of the Earthsystem in our model because they have a high basin stability in both, when exactly two andexactly three elements are tipped (Figs. 4 and 5). Although a potential disintegration of the icesheets can take several centuries up to millennia, states including tipped ice sheets seem to bemore stable than states without tipped ice sheets. This is also consistent with the earlier resultthat the large ice sheets are the initiators of many cascades in the studied model [35].In case exactly two elements are tipped (Fig. 4), the basin of the Greenland and the West Antarc-tic Ice Sheet is the only one which has increased basin volume for low interaction strength andglobal warming levels of 1-3 ˝ C above pre-industrial. This would represent a scenario in which,both, the Greenland Ice Sheet as well as the West Antarctic Ice Sheet are triggered and becomeice free on long time scales without a tipping of the AMOC. This could for example be the casewhen global warming is higher than necessary to safeguard the large ice sheets, but low enoughsuch that the time of their disintegration is slow enough such that the freshwater input into theAMOC does not stop their functioning.In parallel, if exactly three elements are tipped, the combinations that include the Greenlandand the West Antarctic Ice Sheet have a higher basin stability at low interaction strength. Here,global warming levels are up to 4 ˝ C above pre-industrial (Fig. 5).We compare the basin volumes of these scenarios, where exactly two or three tipping elementsare in the tipped regime and the large ice sheets are among these tipped elements (see Fig. 6).25
Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-AMOC tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-WAIS tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-ENSO tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] AMOC-WAIS tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] AMOC-ENSO tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] AMOC-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] WAIS-ENSO tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] WAIS-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] ENSO-AMAZ tipped B a s i n v o l u m e A BC DE FG HI J
Figure 4.
Basin stability for exactly two tipped elements for all pairs of tipping elements (panelA to J). Most of the basin volumes are very small, but the joint basin of the Greenland andthe West Antarctic Ice Sheet for low interaction strength and the joint basin of the AMOC andthe West Antarctic Ice Sheet for high interaction strength and low temperature is increased.The abbreviations in the title stand for: Greenland Ice Sheet (GIS), West Antarctic Ice Sheet(WAIS), Atlantic Meridional Overturning Circulation (AMOC), El-Ni˜no Southern Oscillation(ENSO) and Amazon rainforest (AMAZ). 26
Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-AMOC-WAIS tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-AMOC-ENSO tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-AMOC-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-WAIS-ENSO tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-WAIS-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS-ENSO-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] AMOC-WAIS-ENSO tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] AMOC-WAIS-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] THC-ENSO-AMAZ tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] WAIS-ENSO-AMAZ tipped B a s i n v o l u m e B AC DE FG HI J
Figure 5.
Basin stability for exactly three tipped elements (panel A to J). The volume of basinsis large, where the Greenland and the West Antarctic Ice Sheet are included, for low interactionstrength and global warming levels of up to 4 ˝ C above pre-industrial.27e observe that the basin volume is highest between 1-4 ˝ C above pre-industrial levels foran interaction strength of 0.1. We find that the basin volume is largest at intermediate inter-action strengths d (mainly below 0.2) for a global mean temperature increase of 2 ˝ C abovepre-industrial levels. We also reveal that the basin volume for two tipped ice sheets (red curve)is lower than for exactly three tipped elements including the two ice sheets (other curves). Sincemany basin volumes of exactly two or three tipped elements are very close to zero (see Figs. 4and 5) and the basin volumes including tipped ice sheets are different from zero, this empha-sizes again that the ice sheets could be of special interest for the resilience of the Earth systemwith respect to tipping dynamics.Furthermore, some basin volumes are increased for low to intermediate levels of global warmingand high interaction strengths (above « . ). It is likely that such scenarios are less realisticsince, either such a low increase of the global mean temperature is improbable, or such highinteraction strength would pose the unlikely scenario that interactions are as important as theindividual dynamics of the tipping elements. This would be the case when the interactionstrength d approaches . (see Figs. 4, 5 and compare to Eq. 5).28 Global mean temperature [°C] B a s i n v o l u m e Interaction strength d = 0.10
GIS-WAIS tippedGIS-AMOC-WAIS tippedGIS-WAIS-ENSO tippedGIS-WAIS-AMAZ tipped 0.0 0.2 0.4 0.6 0.8 1.0
Interaction strength d [a.u.] B a s i n v o l u m e Global mean temperature = 2.0°C
GIS-WAIS tippedGIS-AMOC-WAIS tippedGIS-WAIS-ENSO tippedGIS-WAIS-AMAZ tipped A B Figure 6.
Basin volume for selected tipping elements residing in the tipped and untipped regimefor A) an interaction strength of d “ . and B) a global mean temperature increase of 2.0 ˝ Cabove pre-industrial levels. The shaded colors represent the standard deviation arising fromthe 27 different network setups due to the structural uncertainties of the network of tippingelements (compare Sect. 2.4 and see Fig. 1A). The standard deviation increases with increasinginteraction strength d in panel B) since the equilibrium state depends on the initial conditionsof the respective ensemble member. This equilibrium state is influenced significantly morefor higher interaction strengths, but due to stabilizing, destabilizing and unclear links in thenetwork of tipping elements, the change of the equilibrium states fluctuates more and showsa higher uncertainty. Therefore, the standard deviation increases. The same is valid for panelA) at a region where (parts) the critical temperature ranges T limit, i are located for many of thetipping elements, that is mainly between 2-4 ˝ C above pre-industrial.29 .2 Oscillatory states
Furthermore, from the basin stability results we aim to separate off limit cycle attractors in thestate space. The results from MCBB (Monte Carlo Basin Bifurcation [42]) identify the param-eter regimes where Hopf Bifurcations occur and the tipping elements start to show Kadyrovoscillations. Such Kadyrov oscillations have already been found in the early literature on dy-namical systems of the CUSP type [92]. As shown in Fig. 7 for initial conditions at ´ for alltipping elements, this is most prominently the case for large interaction strengths and mediumtemperature increase values. Here, about every tenth solution is oscillating. This is due to thefact that uncertainties are largest in these regimes. For smaller interaction strength values, limitcycles can still occur but are much rarer with an occurrence at about 1% of all solutions. Of allthese limit cycle oscillations almost all (95%) have a significant amplitude (Standard deviation ą .000.020.040.060.080.100.120.14 Interaction strength d [a.u.] Δ G l o b a l m ea n t e m p e r a t u r e [ ° C ] r e l . b a s i n s i z e o s c . s t a t e s A Δ Global mean temperature [°C] B r e l . b a s i n s i z e o s c . s t a t e s Interaction strength d [a.u.] C GISAMOCWAISENSOAMAZmodel time [a.u.] example trajectory GIS-AMOC oscillation D other osc.GIS - AMOC other osc.GIS - AMOC Figure 7.
Oscillating states in case the initial conditions of all tipping elements are r´ , ´ , ´ , ´ , ´ s . For random initial conditions, see Fig. C.1. A) Occurrence of oscillatingstates with respect to global warming and interaction strength parameter d. B) Dependence oflimit cycles and their main type on the temperature increase at high interaction strength (dashedmagenta line in panel A). C) Dependence of limit cycles and their main type on the couplingstrength at a temperature increase of 2 ˝ C above pre-industrial (dashed blue line in panel A). D)Example time series for a limit cycle of AMOC and GIS.31
Discussion
We find that the only dominating stable state in the long term, for large temperature increasesaround and above 4.0 ˝ C above pre-industrial levels, is the one with four or five tipped elements.Our results emphasize that the ice sheets could be of special importance for the stability ofthe climate system regarding their increased basin volume in case more than one element istipped. Based on the known interactions from Kriegler et al. [29] this makes sense, since theinteractions between the ice sheets, especially from Greenland to West Antarctica, are strongdue to potentially rising sea level that might cause grounding line retreat [69].Of course, the ice sheets interact with global modes of ocean variability like the AMOC andreduce its overturning strength, but in our model these interactions are not sufficient to tip theAMOC over in many cases. These states with disintegrated ice sheets are especially relevantexhibiting a high basin volume for intermediate climate warming scenarios consistent with theclimate target of the Paris Agreement that aims at limiting global warming to well below 2 ˝ Cabove pre-industrial levels [108]. Limit cycle oscillations between the tipping of some elementshave been detected at some rare parameter configurations, mainly between the Greenland IceSheet and the AMOC. Although it remains unclear whether such (Kadyrov) oscillations haveoccurred in the climate system, they point to possible relevant internal modes of variability in theclimate system. In principle such limit cycle behavior could have played a role in paleo climatedynamics such as in the Pleistocene ice age cycles [106, 107]. Further, the individual dynamicsare not the sole determinant of the final state of the tipping elements since the network effectscan cause additional tipping events. Through this network interaction, it is therefore possiblethat cascades of tipping events emerge, even before the actual critical temperature threshold forsome of the tipping elements is reached [35]. 32
Conclusion
In this work, we study a conceptual model of five climate tipping elements based on a systemof coupled, nonlinear differential equations. We investigate the stability of different dynamicalregimes with respect to its stable states applying the concept of basin stability using a verylarge-scale Monte Carlo simulation of more than 3.5 billion ensemble members. Following thatapproach, we are able to propagate the numerous uncertainties thoroughly which are associatedwith the critical temperature thresholds and interaction strengths. With a Monte Carlo basinbifurcation analysis tool, we detected oscillatory states within our system.We observe that the largest basin volume is that of the basin, where all five tipping elementsare in the transgressed state, especially for large levels of global warming. We also detect thatthe ice sheets are of special importance for the stability of states, where the large cryospherecomponents reside in the transgressed state, while the other tipping elements do not. We alsodetect Hopf-bifurcations for few parameter configurations (0.6%), mainly taking place betweenthe Greenland Ice Sheet and the AMOC (86%).Our complex dynamical networks approach strongly simplifies the nature of tipping elementsas well as their interaction structure. However, it can serve to integrate simplified conceptsof tipping elements until coupled, process-based models are developed that can resolve therespective nonlinearities in the Earth system in more detail since current state-of-the-art Earthsystem models cannot yet model all these nonlinearities due to a lack of comprehensive process-understanding and computational constraints. It is further important to note that some studieshave hypothesized that major changes in ENSO are possible [60, 61] based on conceptual mod-els [109, 110], but however, whether this is evidence for a permanent and potentially evenirreversibly tipped ENSO remains uncertain and debated. Surely, ENSO exerts strong feed-33acks onto the climate system that will increase if major El-Ni˜no events become more frequent,for instance through strong drying trends over Amazonia. Furthermore, in earlier research wefound that the main results of our model remain robust under the omission of ENSO such thatwe decided to investigate the more complex case and included ENSO here, even though the useof Eq. 5 is only a topologically equivalent dynamical equation [for more details see 35]. Whilesome literature studies present ENSO among the list of potential tipping elements [6, 10, 29], itstill remains uncertain whether ENSO is a tipping element in a strict sense.Overall, our network approach can easily be adapted to further tipping elements as soon astheir interaction structure would be understood. It is also possible to probe the effect of dif-ferent structural interaction hypotheses to further tipping elements within the scope of an un-certainty analysis, as has already been performed here for three interaction links. Further, theresults of our study motivate that it could be worthwhile to look into the dynamics in moredetail using process-detailed Earth system models. Especially the role of the large sheets in thestability landscape and oscillations between climate system components could be of interest.Even though, there is some knowledge about the interaction structure present in literature (seeSect. 2.1), a new expert elicitation might be worthwhile because the knowledge about the in-teractions between the tipping elements has surely widened since the original expert elicitationfrom Kriegler et al. (2009) [29]. 34 eferences [1] Scheffer M 2009
Critical transitions in nature and society vol 16 (Princeton UniversityPress)[2] Watts D J 2002
Proc. Natl. Acad. Sci. Nature
J Prod. Innovat. Manag. Nature
Proc. Natl. Acad. Sci.
Science
J. Roy. Soc. Interface et al. Proc. Natl. Acad. Sci.
Nat. Clim. Change Nature
Earth Syst. Dynam. Disc.
Science
Nat. Clim. Change The Cryosphere Disc.
Nature
Sci. Adv. eaba2949[18] Nobre C A, Sampaio G, Borma L S, Castilla-Rubio J C, Silva J S and Cardoso M 2016 Proc. Natl. Acad. Sci. et al.
Contribution of working group I to the fifth assessmentreport of the intergovernmental panel on climate change [20] Zwally H J, Li J, Brenner A C, Beckley M, Cornejo H G, DiMARZIO J, Giovinetto M B,Neumann T A, Robbins J, Saba J L et al.
J. Glaciol. et al. Nat. Clim. Change Proc.Natl. Acad. Sci.
Nat. Clim.Change et al. Nat. Geosci. Trends Ecol. Evol. The Cryosphere Sci. Adv. e1500589[28] Robinson A, Calov R and Ganopolski A 2012 Nat. Clim. Change Proc. Natl. Acad.Sci.
Proc. Natl. Acad. Sci.
Nat. Clim. Change Nat. Clim. Change A question of balance: Weighing the options on global warmingpolicies (Yale University Press)[34] Nordhaus W 2014
J. Assoc. Environ. Resour. Econ. Earth Syst. Dynam. Disc.(Preprint doi: 10.5194/esd-2020-18) [36] Menck P J, Heitzig J, Marwan N and Kurths J 2013
Nat. Phys. New J. Phys. Sci. Rep. Phys. Rev. E New J. Phys. Sci. Rep. New J. Phys. Trends Ecol. Evol. Chaos Phys. Rev. E
Science
Sci. Rep. Clim. Dynam. Nature
Tellus et al. Geophys. Res. Lett. [52] Hawkins E, Smith R S, Allison L C, Gregory J M, Woollings T J, Pohlmann H andDe Cuevas B 2011 Geophys. Res. Lett. [53] Mecking J, Drijfhout S, Jackson L and Graham T 2016 Clim. Dynam. Clim. Dynam. Geophys. Res. Lett. [56] Cox P M, Betts R, Collins M, Harris P P, Huntingford C and Jones C 2004 Theor. Appl.Climatol. Science
Glob. Change Biol. Ecol. Complex. J. Climate Earth Syst. Dynam. Geophys. Res. Lett. [63] Driesschaert E, Fichefet T, Goosse H, Huybrechts P, Janssens I, Mouchet A, MunhovenG, Brovkin V and Weber S 2007 Geophys. Res. Lett. [64] Robson J, Hodson D, Hawkins E and Sutton R 2014 Nat. Geosci. et al. J. Climate et al. J. Climate Climatic Change J. Geophys.Res.-Oceans
Geophys. Res. Lett. Science
ClimaticChange
J. Climate J. Climate Geophys. Res.Lett. [75] Holmgren M, Hirota M, Van Nes E H and Scheffer M 2013 Nat. Clim. Change et al. Front. Ecol. Environ. Philos. Trans. R. Soc. B
Earth Syst. Dynam. [79] Bertler N, Naish T, Mayewski P and Barrett P 2006 Adv. Geosci. Nat.Geosci. et al. Nat. Communs. Clim. Dynam. Geophys. Res. Lett. Journal of climate Geophys. Res. Lett. Global Planet. Change Geophys. Res. Lett. [88] Swingedouw D, Fichefet T, Goosse H and Loutre M F 2009 Clim. Dynam.
Nature
Sci. Rep. Nature Communications Int. J. Bifurc. Chaos Roy. Soc. Open Sci. et al. Nat. Methods Phys. Rev. E Phys. Rev. Lett. Nonlinear physical oceanography: a dynamical systems approachto the large scale ocean circulation and El Nino vol 28 (Springer Science & BusinessMedia)[98] Klinshov V V, Nekorkin V I and Kurths J 2015
New J. Phys. Phys. Rev. E Earth Syst. Dynam. Sci. Rep. Phys. Rev. E Sci. Rep. Commun. Math. Phys. https://pythonhosted.org/pyDOE/index.html accessed: 2020-06-08[106] Ditlevsen P, Mitsui T and Crucifix M 2020 Clim. Dynam. Philos. Trans. Royal Soc. A https://treaties.un.org/pages/ViewDetails.aspx?src=TREATY&mtdsg_no=XXVII-7-d&chapter=27&clang=_en [109] Timmermann A, Jin F F and Abshagen J 2003
J. Atmos. Sci. Mon. Weather Rev.
Geosci. Model Dev. cknowledgements This work has been carried out within the framework of PIK’s FutureLab on Earth Resiliencein the Anthropocene. N.W., M.G. and R.W. acknowledge the financial support by the IRTG1740/TRP 2015/50122-0 project funded by DFG and FAPESP. N.W. is grateful for a scholar-ship from the Studienstiftung des deutschen Volkes. J.F.D. is grateful for financial support bythe Stordalen Foundation via the Planetary Boundary Research Network (PB.net), the EarthLeague’s EarthDoc program and the European Research Council Advanced Grant project ERA(Earth Resilience in the Anthropocene, ERC-2016-ADG-743080). We are thankful for financialsupport by the Leibniz Association (project DominoES). The authors gratefully acknowledgethe European Regional Development Fund (ERDF), the German Federal Ministry of Educa-tion and Research and the Land Brandenburg for supporting this project by providing resourceson the high performance computer system at the Potsdam Institute for Climate Impact Research.
Author contributions
CRediT (Contributor Roles Taxonomy) statement:
N.W.:
Conceptualization, Formal analy-sis, Investigation (Basin stability, Monte Carlo), Visualization (Introduction, Methods, Basinstability), Writing - Original Draft, Writing - Review & Editing.
M.G.:
Formal analysis, Inves-tigation (Monte Carlo Basin Bifurcation), Visualization (Oscillatory states), Writing - Review& Editing.
R.W.:
Conceptualization, Writing - Review & Editing, Supervision, Funding ac-quisition.
J.K.:
Writing - Review & Editing, Funding acquisition.
J.F.D.:
Conceptualization,Writing - Review & Editing, Supervision, Funding acquisition.
Code and data availability
The data that support the findings of this study are available from the corresponding authorupon reasonable request. The code for the Monte Carlo ensemble construction and the con-44eptual Earth system that support the findings of this study are freely (3-clause BSD license)available on github under the following doi: 10.5281/zenodo.4153102. The algorithm onthe Monte Carlo Basin Bifurcation (MCBB) is available directly from the GitHub repositoryhttps://github.com/ma-ximilian-gelbrecht/MCBB.jl/. For the use of MCBB, please also con-fer [42].
Note on color maps
This paper makes use of the conceptually uniform colormaps developed by [111].45 ppendixA Parameter uncertainties
In the following tables (Tabs. A.1 and A.2), we list the critical temperatures T limit, i for therespective tipping element and the interactions between them together with their uncertainties. Interaction Link strength range s ij (a.u.) Process Greenland Ñ AMOC r` ` s Freshwater inflowAMOC Ñ Greenland r´ ´ s AMOC breakdown, Greenland coolingGreenland Ñ West Antarctica r` ` s Grounding line retreatENSO Ñ Amazon rainforest r` ` s Drying over AmazoniaENSO Ñ West Antarctica r` ` s Warming of Ross and Amundsen seasAMOC Ñ Amazon rainforest r˘ ˘ s Changes in hydrological cycleWest Antarctica Ñ AMOC r˘ ˘ s Increase in meridional salinity gradient ( ´ ),Fast advection of freshwater anomalyto North Atlantic ( ` )AMOC Ñ ENSO r` ` s Cooling of North-East tropical Pacific with thermo-cline shoaling and weakening of annual cycle in EEPWest Antarctica Ñ Greenland r` ` s Grounding line retreatENSO Ñ AMOC r´ ´ s Enhanced water vapor transport to PacificAMOC Ñ West Antarctica r` ` . s Heat accumulation in Southern OceanAmazon rainforest Ñ ENSO r˘ ˘ . s Changes in tropical moisture supply
Table A.1.
Each interaction in the network of Fig. 1 has a specific link strength range and aspecific physical process that is connected to the respective interaction. The link strength rangesare computed from literature values [43, 29] such that they can be used in Eq. 5. For a more indepth description please be referred to Wunderling et al. (2020) [35].
B More basin stability results
Here, we show the standard deviation of the basin volume for 0 to 5 tipped elements (Fig. B.1)and the basin volume for one specific tipped element (Fig. B.2) to complement the results fromFig. 3. 46
Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e A BC DE F
Figure B.1.
Standard deviation over the 27 different network realizations of the basin volumenormalized to one for each number of tipped elements (panels A-F) in dependence of interactionstrength and global mean temperature increase. Mean values can be found in Fig. 3.47
Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] GIS tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] AMOC tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] WAIS tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] ENSO tipped B a s i n v o l u m e Global mean temperature [°C] I n t e r a c t i o n s t r e n g t h d [ a . u . ] AMAZ tipped B a s i n v o l u m e A BC DE F
Figure B.2.
Basin stability for single elements showing that the basin volume for the AMOC,WAIS (West Antarctic Ice Sheet) and AMAZ (Amazon rainforest) are qualitatively similar sinceit is possible that only this particular element tips. However, for ENSO and the GIS (Green-land Ice Sheet) this is not the case. It can onyl very rarely happen that these elements tip onthemselves for high interaction strength even at low temperature increases since both of thempossess a very strong link to another element that they would draw along into the tipped state.For ENSO, this is the Amazon rainforest and for the Greenland Ice Sheet it is the AMOC. Notethat the color bar is different for panel A) to improve visibility.48 ipping element T limit, i [ ˝ C] τ i [a.u.] Greenland 0.8 – 3.2 4900West Antarctica 0.8 – 5.5 2400AMOC 3.5 – 6.0 300ENSO 3.5 – 7.0 300Amazon rainforest 3.5 – 4.5 50
Table A.2.
Critical temperature range T limit, i of the five tipping elements as taken from theliterature [10], see also Eq. 5. The typical tipping time scale τ i is given in model years (inarbitrary units) since it is beyond the scope of this model to make predictions about the exacttipping times. However, certain differences in tipping times as used here can be decisive whethera tipping event occurs or not. For more information see Wunderling et al. (2020) [35]. C Oscillatory regimes for random initial conditions
Here, we show the results of a Monte Carlo Basin Bifurcation analysis for random initial con-ditions (Fig. C.1). We find that limit cycles occur more frequently when the initial conditionsare randomly shuffled. 49 .0 0.2 0.4 0.6 0.8