Basin stability and limit cycles in a conceptual model for climate tipping cascades

Tipping elements in the climate system are large-scale subregions of the Earth that might possess threshold behavior under global warming with large potential impacts on human societies. Here, we study a subset of five tipping elements and their interactions in a conceptual and easily extendable framework: the Greenland and West Antarctic Ice Sheets, the Atlantic Meridional Overturning Circulation (AMOC), the El-Nino Southern Oscillation (ENSO) and the Amazon rainforest. In this nonlinear and multistable system, we perform a basin stability analysis to detect its stable states and their associated Earth system resilience. Using this approach, we perform a system-wide and comprehensive robustness analysis with more than 3.5 billion ensemble members. Further, we investigate dynamic regimes where some of the states lose stability and oscillations appear using a newly developed basin bifurcation analysis methodology. Our results reveal that the state of four or five tipped elements has the largest basin volume for large levels of global warming beyond 4 {\deg}C above pre-industrial climate conditions. For lower levels of warming, states including disintegrated ice sheets on West Antarctica and Greenland have higher basin volume than other state configurations. Therefore in our model, we find that the large ice sheets are of particular importance for Earth system resilience. We also detect the emergence of limit cycles for 0.6% of all ensemble members at rare parameter combinations. Such limit cycle oscillations mainly occur between the Greenland Ice Sheet and AMOC (86%), due to their negative feedback coupling. These limit cycles point to possibly dangerous internal modes of variability in the climate system that could have played a role in paleoclimatic dynamics such as those unfolding during the Pleistocene ice age cycles.


Introduction
During the last decades, the field of tipping elements has become a major point of interest in complex systems and network science [1,2]. They have been used in the description of various fields 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 could potentially lead to catastrophic domino effects [8] or, for instance, lead to a hothouse climate state in the case of climate tipping elements [9]. In the climate system, tipping elements are subregions of the Earth system that can exhibit threshold 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 or system due to internal, self-enforcing feedbacks [6]. Climate tipping elements comprise systems from the cryosphere (e.g. Greenland, Antarctic Ice Sheet, Permafrost), the biosphere (e.g. Amazon rainforest (AMAZ), coral reefs) and large-scale circulation systems (e.g. Monsoon systems, AMOC) [6,10]. Their potential tipping to alternative states would 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 transgressing into an undesired state even at global warming levels below the 2.0 • C goal of the Paris agreement [10][11][12]. Among others, tipping elements that already show warning signals of degradation at present times [11,12] are: the West Antarctic Ice Sheets (WAIS) where parts in the Amundsen Bay (Pine Island & Thwaites region) are suspected to have been destabilized [13][14][15], the AMOC which experienced a major slowdown of 15% from 1950 to now [16], the AMAZ which might approach a tipping point due to climate change and deforestation [17]. Critical deforestation ratios might lie between 20% to 40%, where current deforestation is reaching 20% [17,18]. Furthermore, the GIS loses mass at an accelerating pace [19,20] and the frequency of major El-Niño events are suggested to increase twofold and strong El-Nino Southern Oscillation (ENSO) effects will occur more often as global warming continues [21,22]. However, others highlight that large uncertainties are related to future changes of ENSO and whether major El-Niño events will become more frequent or intense under global warming [23,24].
Furthermore, contradicting a common misunderstanding, tipping elements do not necessarily tip immediately after the crossing of their tipping point, but their tipping time trajectory might take very long and appear smooth [25]. For instance for the large ice sheets, the disintegration time scale could be on the order of several centuries up to millennia as has been suggested by modeling studies [26][27][28].
For most of the tipping elements, there is a critical temperature range at which they are suspected to 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 Paris agreement (Paris, 2015) of 2 • C above pre-industrial and with that, these are the most vulnerable climate tipping elements with respect to global warming. This group contains mainly cryosphere elements (Arctic summer sea ice, WAIS, GIS and Alpine glaciers) as well as the Coral reefs that are likely to be lost even when global warming is restricted to 2 • C above pre-industrial. The second group might tip at temperatures above 3 • C (for instance the AMAZ, AMOC or ENSO) and the most resilient group only at 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, but connected [11,29] and the knowledge about the exact interaction structure is sparse and partially based on experts that, for instance, suggested an interaction structure, including sign and strength, for a subset of five tipping elements: the GIS, the WAIS, the AMOC, the AMAZ and the ENSO [29]. Behind each connection between two tipping elements within this subset, there is a physical process or set of processes (see table A.1). For instance the impact of the GIS on the AMOC due to freshwater input from melting ice slows down the AMOC on the one hand and a weakening AMOC on the other hand cools latitudes in the northern hemisphere. Note that this subset network of tipping is neither complete in the number and selection of tipping elements, nor is it comprehensive in the possible connection pathways and their potential strength between the tipping elements. There are 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 versions of the integrated assessment model DICE [33,34]. Here, Cai et al (2016) [31] explicitly base their findings on the interactions of tipping elements from the expert elicitation in Kriegler et al (2009) [29].
Following the elaborations above, in this work we aim at investigating the resilience of various attractors for interacting climate tipping elements and we want to elucidate the role that different tipping cascades have in that regard. The approach put forward here can easily be adapted to more tipping elements and further interaction structures once they are more comprehensively understood (see also [35]).
We explore the stability landscape and the dynamics of a subset of five tipping elements represented by normal form fold bifurcations based on known interactions across scales in time and space between these tipping elements (figure 1) [29]. These five tipping elements are the GIS and WAIS, the AMOC, the ENSO and the AMAZ [29]. We introduce the model of interacting tipping elements in section. We use the concept of basin stability [36] in order to determine the basin sizes of various attractors of this multistable system (see section 2). Based on a large-scale Monte Carlo ensemble, this methodology gives an estimate how stable and resilient various attractors are. It has been applied to many dynamic systems before such as power grids,  [29,43]. Each link represents a physical mechanism and has a certain strength (see supplementary table A.1). A positive arrow implies an effect that drives the tipping element closer to its tipping point, a negative arrow drives the tipping 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 global mean temperature levels for the respective tipping element. neuronal models and further nonlinear systems [37][38][39][40][41]. Furthermore, especially for larger coupling strengths, Hopf bifurcations can occur, thus invoking oscillatory limit-cycle solutions of the model. For the detection and quantification of these types of limit cycle attractors, we apply a newly developed bifurcation algorithm that is able to identify different dynamical properties in complex systems: the Monte Carlo Basin Bifurcation analysis (MCBB) [42].
In the following, we first introduce the methodological approach of this work: the model of interacting tipping elements (section 2.2), the basin stability approach (section 2.3), the construction of the large-scale Monte Carlo ensemble (section 2.4) and the Monte Carlo Basin Bifurcation 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 with respect to the climate system are discussed and summarized in sections 4 and 5.

Methods
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 the main nonlinear dynamical properties of tipping elements [35,44,45]. The actual physical processes behind the tipping elements are not explicitly modeled to maintain an accessible and controllable structure. The modeling of complex systems using conceptual approaches is a popular tool and has been successfully applied to, among others, ecology, social systems or epidemiology [5,46,47].

Tipping elements and interactions
Given that, despite major advances, current EMICs (Earth system models of intermediate complexity) and GCMs (global circulation models) are not yet able to fully represent the nonlinear behavior of some Earth system components together with their interactions, but physics based models and equations as well as paleo climate observations suggest the existence of such properties for many tipping elements as for instance the GIS and (West) Antarctic ice sheet [26,28,48,49], the AMOC [50][51][52][53][54], the AMAZ [55][56][57][58][59] or the ENSO [60,61], the conceptual approach chosen here demonstrates an option how to model interactions between tipping elements. Thus, we put this model forward as a first step towards a more process-detailed assessment of tipping elements and their interactions. This also emphasizes that future research could focus on developing more complex, emulator-or EMIC/GCM-like models of tipping elements to investigate their nonlinear interplay such as has recently been developed for the Antarctic ice sheet [49].
While we have described why we use a conceptual description of the main dynamics of the tipping elements directly above, we outline the physical mechanisms of there interactions hereafter. Although some interactions between the tipping elements are better understood and evaluated than others, we describe the physical mechanisms of all of them in the following separated into destabilizing (+), stabilizing (−) and unclear links (?) (see figure 1). This description aims to provide a basic physical understanding, but cannot resolve the problem of how strong exactly each of these interactions is. Therefore, the interaction structure is kept as described later in equation (5) with the multiplicative interaction strength factor d.
(a) Destabilizing interactions: 1. GIS → AMOC: when the GIS starts to melt, it has a diminishing influence on the overturning strength of the AMOC due to freshwater input into the North Atlantic. This has been observed in modeling studies [16,51,62,63] and in observations [64]. 2. AMOC → WAIS: when the AMOC collapses, sea surface temperature 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 studies [65][66][67][68]. 3. GIS → WAIS (and vice versa): the shift of grounding 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 could be lifted which is likely to result in grounding line retreat. Furthermore, gravitational changes as well as elastic and rotational effects might then amplify the sea level change if one of the large ice sheets disintegrates first because the gravitational attraction then only emanates from the other, remaining ice sheet [70,71]. The effect would be stronger if Greenland melts first since the WAIS has more marine terminating glaciers and ice shelves. 4. AMOC → ENSO: there are two opposing effects that have proposed, which describe how the AMOC might influence the ENSO: (i) it has been suggested that oceanic Kelvin waves originate from a colder North Atlantic and travel southward. Then, in western Africa, Rossby waves would be emitted towards the north and the south, which are then translated back into Kelvin waves that travel into the Pacific ocean. This effect would deepen the Pacific thermocline and weaken the amplitude of ENSO [60]. (ii) With a weaker AMOC, the northern tropical Atlantic would in turn become cooler and northerly trade winds would be intensified over the northeastern tropical Pacific. It has been argued that this could result in a southward displacement of the Pacific ITCZ leading to a sea surface temperature anomaly [72]. At the same time, it is argued that Rossby waves are sent into the northeast tropical Pacific [73]. This would intensify ENSO due to wind stress interaction from AMOC. Overall, it is believed that mechanism (ii) is stronger than mechanism (i) [60]. Furthermore, it can be extracted from complex Earth system models that a decrease in AMOC intensity indeed strengthens the variability of ENSO [61,74]. 5. ENSO → AMAZ: literature studies suggest that droughts related to climate variabilities such as the ENSO can affect the stability of the AMAZ [75][76][77]. Using an EMIC, it has been found that a permanent El-Niño state would endanger substantial portions of the Amazon basin due to a reorganization and reduction of water access in the south American tropics via teleconnections [78]. 6. ENSO → WAIS: the interaction between ENSO and the WAIS is one of the least certain interactions as has already been stated in Kriegler et al (2009) [29]. Nevertheless, there are hints for warming oceanic effects from El-Niño in the Amundsen and Ross Sea region, while La Niña would cool this region. At the same time, atmospheric effects could have an opposite effect, which would offset the oceanic effect [79]. Besides that, it has been found with satellite observations that ice shelves gain height, but yet lose mass during El-Niño events in the Amundsen and Ross sea region [80]. While the primary driver of melt in west Antarctica is the warm ocean water below ice shelves, an extended period of surface melting has been observed during January 2016, which is likely promoted by the strong El-Niño event in this year [81]. Since it is expected that the frequency of major El-Niño events will increase during climate change [21], we set this interaction positive (see table A.1 and figure 1). (b) Stabilizing interactions: 1. AMOC → GIS: for a decreasing overturning strength of the AMOC, the northern hemisphere is cooled since the heat transport towards the north Atlantic would be weakened. This has been observed in modeling studies [16,65,66,82]. 2. ENSO → AMOC: using reanalysis data, evidence has been found that the transport of water vapor out of the tropical Atlantic is enhanced [83]. Comparing La Niña and El-Niño conditions, it was found in this study that El-Niño conditions lead to a stronger northern AMOC on a multi-decadal timescale. However, another study questions this finding and does 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 low strength, but with a negative sign (see figure 1 and table A.1). (c) Unclear interaction direction: 1. AMOC → AMAZ: when the AMOC shuts down, the intertropical convergence zone (ITCZ) is likely dislocated southward, leading to large changes in seasonal precipitation on a local to regional degree. This might then impact parts of the AMAZ [68,82,85]. Still, it is unknown as to whether this interaction is positive or negative and might differ from region to region. Therefore, this link is set as unclear (see figure 1). 2. WAIS → AMOC: a literature study using a coupled ocean-atmosphere model found a decrease in the AMOC for high freshwater inputs from the WAIS [86]. However, another study detected a stabilization of the AMOC if influenced by freshwater input from west Antarctica. This is ascribed to the effects from the bipolar ocean 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 that different processes could enhance or slow down the AMOC [88]: (i) the deep water adjustments via the bipolar ocean seesaw tend to intensify the NADW formation. (ii) The NADW is strengthened by southern hemispheric wind increase representing an ocean-atmosphere interaction. (iii) Salinity anomalies from the southern ocean are distributed to 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 scale and efficiencies of these processes have been rated unknown as of yet [88]. 3. AMAZ → ENSO: under a dieback of the AMAZ, the moisture supply to the atmosphere will significantly change, also since the atmospheric moisture recycling feedback over the Amazon basin would break down [89][90][91]. However, it is unclear whether and to which extent this would then impact ENSO.

Model
In our conceptual model, we divide the dynamics x i of the considered tipping elements i into their individual dynamics f i (x i ) and a direct interaction term g i ( x) ≡ g i (x). This yields where τ i is the typical time that passes when a tipping element undergoes a critical transition from one state to another. We model the individual dynamics of each of the tipping elements with the general tipping approach (CUSP equation [8,92]) where a i > 0 and b i > 0. Assuming additive separability of the interactions between the tipping elements and linear interactions, the interaction term g i (x) becomes Here, A ij is the interaction structure and strength, which is set to zero if there is no connection between the tipping elements i and j. Altogether, equation (1) becomes Each tipping element x i following this equation possesses two fold bifurcations at ± (4a 3 i )/(27b i ) and has already been investigated in theoretical works on tipping cascades [92], but also in various contexts where nonlinear behavior is important as for instance in policy, environmental issues, economy or climate [8,93]. For these equations exist a framework that allows to investigate tipping cascades on larger networks with regard to their interaction structure in the network as 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 a multiplicative factor d times the actual link strength s ij between each pair of tipping elements. Therefore, A ij = d/10 · s ij . The link strength values s ij are taken from the expert elicitation [29]. The factor 1/10 is used for normalization reasons since then d ∈ [0, 1]. If we now additionally set a = 1, b = 1 and c i = 4/27/T limit, i · ΔGMT, the tipping elements are described by the following nonlinear, ordinary differential equation (all parameters of equation (5) are explained in the tables A.1 and A.2) Here, x i is the state of the tipping element (see figure 1(B)) and i stands for the considered tipping elements i = GIS, WAIS, AMOC, ENSO, AMAZ. We choose these five tipping elements since their interaction structure is known from an expert elicitation [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 is the coupling term, where d is a general multiplicator that determines the strength of the interaction term in comparison to the other, individual dynamics terms. The parameter d is varied between 0, meaning no interactions, and 1, where the interactions become as important as the individual dynamics. Following this, one might tend to assume that the individual dynamics of the tipping element influences the tipping element more than the interaction effect. This might make smaller coupling parameters more realistic than higher ones. In equation (5), s ij is the link strength that is based on the expert elicitation [29] and τ i is a typical timescale at which a certain tipping element transgresses its state. This typical tipping time scale ranges from decades for the AMAZ to several millennia for the large ice sheets (see appendix table A.2). Then, our system of differential equations is integrated forward in time using scipy.odeint [94] until more than 20 times the GIS'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 intending to compute an exact time scale for tipping or tipping cascades here, but we are rather interested in the system's attractors and their stability properties. This is why we denote model years in arbitrary units instead of giving an exact time, also since this would be beyond the scope of this conceptual model (see figure 1(A)). Note that we adapted the link from ENSO to AMOC from uncertain to negative compared to the original results of the expert elicitation on tipping element interactions [29] since there is only a dampening process known in literature [43].
There are considerable uncertainties associated with this approach, especially with the critical temperature at which a certain tipping element transgresses its state T limit,i as well as in the strength of the interactions s ij . The uncertainties of these two parameters are shown in the appendix tables A.1 and A.2. Thus, with equation (5), we model tipping events and cascades under certain conditions of global warming (GMT) and the interaction strength (d).

Basin stability
We are interested in the stability properties of different attractors within the state space. An appropriate tool to investigate the stability landscape of such states is the so-called basin stability [36,39]. Basin stability is a nonlinear stability measure for the resilience of an attractor to disturbances. Where traditional measures such as the computation of Lyapunov exponents or master stability functions rely on linear approximations in reaction to small perturbations [95,96], basin stability approaches can also consider large perturbations. Such perturbations can occur in Earth system components such as the large ocean circulations or the AMAZ [36,97]. The basin stability is an established algorithm focusing on the stability landscape of the entire phase space, while other nonlinear stability measures such as survivability [41], stability threshold [98], constrained basin stability [99] and topology of sustainable management [100] approaches focus on the stability of parts of the state space or desired regimes in it. Therefore, basin stability computations are a first step that aims to quantify the stability of different attracting states, but do not aim to study potential desired regimes as would be required for the other mentioned methods. The concept of basin stability has been applied to many multistable systems. Examples comprise the AMAZ [36], the stability in networks of power grids [37], neuronal models [38] and further nonlinear systems such as in coupled network systems [39], oscillators [101,102] or chimera states [103]. While basin stability can widely be applied, it has its limitations, for instance in cases where basins become too peculiar, e.g. for riddled basins with holes [40]. Since this is not the case in our model of interacting 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 entire state space [104]. B(A) ⊆ X is the basin of attraction of A which comprises all states from which the system converges to A. The basin stability or the basin volume V(A) is then quantified as the probability that a system will return to a certain attractor A after a perturbation where 1 B(A) is 1 in case x ∈ B(A) 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. The estimation of the integral in equation (6) can be difficult, but in our system it can be assumed that the estimation of the basin volume can be estimated via a Monte Carlo ensemble. The total volume of a basin of attraction is then measured as the fraction of simulations with randomly chosen initial conditions that end up in that certain attractor N(A) over the total number of initial conditions N(Ω) Here, P(A) is the probability that a random initial condition ends up in the basin of attractor A. To assign the basin volume V(A) with the probability P(A), it is required that the space of initial conditions is covered well and uniformly. Therefore in this work, it is necessary to extend the classical concept of basin stability since it is not only required to sample the space of initial conditions sufficiently well, but also to sample over the uncertainties in the model parameters themselves (see tables A.1 and A.2). Thus, we need to set up a very large-scale Monte Carlo ensemble of several billion ensemble members whose construction details can be directly found below.

Monte Carlo ensemble to compute basin stability
In order to apply the concept of basin stability in a meaningful way, the state space must be covered well enough. However, in this application, the parameters of the models have uncertainties themselves in the critical temperature thresholds and the interaction strength and structure. This means, we need a way of covering the many uncertainties in these various parameters as well as the state space itself. Therefore, it is necessary and useful to combine basin stability with a large scale Monte Carlo sample that covers an adequate extent of the phase space and parameter space. 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 sample 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 tables A.1 and A.2). Latin hypercube sampling is an extension to the usual random sampling and is used to improve the space coverage of initial conditions. Therefore, the space of initial conditions is separated into its dimensions, i.e., the number of different initial parameters (here [17], see tables A.1 and A.2). Then, it is secured that only one sample occurs in each axis hyperplane (compare to the N-rooks problem in mathematics). We apply this sampling procedure for each of the 27 different network setups that arise from the permutation (positive, negative, zero) of the three uncertain links (AMAZ → ENSO, AMOC → Amazon and WAIS → AMOC, see figure 1(A)). This then leads to 2700 samples. These 2700 samples are computed for each global mean temperature increase up to 8 • C above pre-industrial which can be reached in business as usual scenarios RCP8.5 extended from 2100 to 2500 [10] in steps of 0.1 • C and coupling constant d between 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: the basin stability of the system for each of these 864.000 samples is computed by permuting the initial state of each of the five tipping elements within its limit, i.e., between the untipped (x = −1.0) and the tipped state (x = +1.0). The state variables of the five tipping elements result in a five dimensional state vector vec basin stability, i = {x Greenland ice sheet (0), x AMOC (0), x west Antarctic ice sheet (0), where vec basin stability, i ∈ [−1.0; +1.0] for each tipping element. However, vec basin stability,i cannot be permuted in a completely random way, but each of its five dimensions needs to be permuted in an independent way since there is a strong nonlinearity at state equal zero for each of the five dimensions. Of course, in principle if there would be infinite computational resources, we would not need to take this nonlinearity into account, but would be able to increase the size of the 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 be equally probable, i.e.: This can be achieved when random starting conditions are drawn from each of the 32 combinations of P sign . Hence, for each of the 32 combinations, we chose 10 different initial conditions ending up with 320 different settings. For the 320 randomly chosen perturbations (i.e., the initial conditions of the tipping elements), we again used a latin hypercube algorithm [105]. That means it fulfils the condition that each of the 32 different possible signs of the initial conditions in their five-dimensional subspace (one dimension for each tipping element) is covered equally often.
Altogether, we employ a very large ensemble of simulations to compute the basin stability of N total = 320 · 864.000 = 3.569.184.000 ≈ 3.6 · 10 9 samples. How the final state can depend on the initial conditions is shown exemplary for three timelines in figures 2(A)-(C).

Monte Carlo basin bifurcation analysis
Coupling nonlinear ODEs as in the model described here, invokes the possibility of further types of bifurcations besides fold bifurcations. Here, we utilize Monte Carlo basin bifurcation analysis [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 a novel, numerical approach to analyze multistable systems, quantify and track their asymptotic states 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 measure of interesting dynamical properties (here occurring limit cycles, i.e., Hopf-bifurcations), it is a well suited method for our purposes. It has also been applied to other nonlinear systems such as 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 attractions of the system. Similar attractors, at different parameter values, have to share similar values of invariant measures ρ and the difference of theses measures has to smoothly vanish if the parameter difference goes to zero. If this is the case they are regarded as being part of the same class of attractors. N tr trajectories of the system, here 140 000, with randomized initial conditions 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, the mean, variance and the Kullbach-Leibler divergence to a normal distribution. Hence, for every statistic i, S i is a N tr × 5 matrix. A distance matrix of each trajectory to each other is computed from these statics with where p (i) is the control parameter used to generate the ith trajectory and w i are free parameters of the method, here [1; 0.5; 0.5; 1] which is the default recommendation for these parameters. This distance matrix is used as an input for a density-based clustering algorithm such as DBSCAN which can find if this notion of continuity between different trajectories exists and thus each cluster corresponds to a different class of attractors. For further details on MCBB, refer to [42]. When applying this to the conceptual model for climate tipping points, not only the different possible states of tipped elements are found, but also different classes of oscillating states induced by Hopf-bifurcations are found. For the MCBB analysis, the parameter uncertainties were varied randomly within the same bounds as for the previously described basin computations. The initial conditions of five tipping elements were chosen to all start at −1, i.e., not tipped for the results presented in the main text, and at random between −1 and +1 for the results presented in the appendix. The computations are performed with the Julia library MCBB.jl.

Basin stability
We compute the basin stability of each potential state that could be governed by the network of five tipping elements. The present day state could be considered as some kind of safe state for the Earth system when all five tipping elements are in a negative state. On the other hand there could be a state where all five tipping elements reside in the positive, tipped state. In between there are intermediate scenarios, where some tipping elements already crossed their thresholds and others did not. In figure 3, we show the average basin stability for each of these six possible situations, i.e., with zero, one, two, three, four and five tipped elements. In this experiment, we perturbed the initial conditions of all tipping elements at the same time.
The fraction of initial conditions 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 levels becomes larger for a higher number of tipped elements as would be expected. For high levels of warming, the basin of five tipped elements dominates.
For increasing interaction strength, the volume of the basins with three or less tipped elements decreases (figures 3(A)-(D)). Contrasting this, the basin volume with four tipped elements increases with increasing interaction strength, while the basin for five tipped elements first increases and then decreases again (figures 3(E) and (F)). The last issue is due to the strong negative feedback loop between the GIS and the AMOC. In such cases of high coupling, the AMOC tips, but safeguards the GIS which reaches the untipped regime for global mean temperature increases above 4 • C and interaction strengths above 0.5. This poses a hypothetical scenario which would only be realistic if the interaction strength between Greenland and AMOC is very high, but this behavior has also been observed in experiments of tipping cascades earlier [35].
For instance, in the basin volume plot of zero, one or two tipped elements, the number of states that equilibrate in this state is very small (figures 3(A)-(C)). For temperature increases above 2 • 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 (figure B.1). We find that uncertainties generally increase for a higher amount of tipped elements as well as for higher interaction strengths. The standard deviation is highest for small temperature increases and high coupling strengths since here, the attractors depend on the initial conditions in terms of the critical temperature thresholds and initial coupling constants (see table A.1). The basin of four and five tipped elements show a regime of increased standard deviation for temperatures around 2 • C-5 • C above pre-industrial and interaction strength parameters of more than 0.2. This is probably due to the fact that in this regime the state of the GIS 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 well as on parameters (see tables A.1 and A.2). Outside and around this regime, the uncertainty is smaller 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 elements can transgress their state without triggering a tipping cascade. This range is mostly located below 1 • C above pre-industrial for low coupling strength and well below 1 • C for higher interaction strengths (see figures 3(B) and B.2). If we separate this response into the respective singular tipping elements, we can see that above an interaction strength of 0.2-0.4, the GIS and ENSO cannot tip without causing a cascade due to their strong interactions links to AMOC or the AMAZ, respectively (for more details see appendix B).
Additionally, we investigate some important intermediate states in more detail, where some elements are in the tipped regime, while others are not. It was found that several tipping cascades of size two and three are more frequent than others, for instance a tipping cascade between the GIS and the WAIS is more likely than, for instance, a cascade between the AMOC and the AMAZ [35]. Thus, we investigate the basin volume that corresponds to such cascades.
We find that the ice sheets appear to be of particular importance for the stability of the Earth system in our model because they have a high basin stability in both, when exactly two and exactly three elements are tipped (figures 4 and 5). Although a potential disintegration of the ice sheets can take several centuries up to millennia, states including tipped ice sheets seem to be more stable than states without tipped ice sheets. This is also consistent with the earlier result that the large ice sheets are the initiators of many cascades in the studied model [35].
In case exactly two elements are tipped (figure 4), the basin of the GIS and the WAIS is the only one which has increased basin volume for low interaction strength and global warming levels of 1 • C-3 • C above pre-industrial. This would represent a scenario in which, both, the GIS as well as the WAIS are triggered and become ice free on long time scales without a tipping of the AMOC. This could for example be the case when global warming is higher than necessary to safeguard the large ice sheets, but low enough such that the time of their disintegration is slow enough such that the freshwater input into the AMOC does not stop their functioning.
In parallel, if exactly three elements are tipped, the combinations that include the GIS and the WAIS have a higher basin stability at low interaction strength. Here, global warming levels are up to 4 • C above pre-industrial ( figure 5).
We compare the basin volumes of these scenarios, where exactly two or three tipping elements are in the tipped regime and the large ice sheets are among these tipped elements (see figure 6). We observe that the basin volume is highest between 1 • C-4 • C above pre-industrial levels for an interaction strength of 0.1. We find that the basin volume is largest at intermediate interaction strengths d (mainly below 0.2) for a global mean temperature increase of 2 • C above pre-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). Since many basin volumes of exactly two or three tipped elements are very close to zero (see figures 4 and 5) and the basin volumes including tipped ice sheets are different from zero, this emphasizes again that the ice sheets could be of special interest for the resilience of the Earth system with respect to tipping dynamics.
Furthermore, some basin volumes are increased for low to intermediate levels of global warming and high interaction strengths (above ≈ 0.5). It is likely that such scenarios are less realistic since, either such a low increase of the global mean temperature is improbable, or such high interaction strength would pose  (5)).

Oscillatory states
Furthermore, from the basin stability results we aim to separate off limit cycle attractors in the state space. The results from MCBB [42] identify the parameter regimes where Hopf Bifurcations occur and the tipping This equilibrium state is influenced significantly more for higher interaction strengths, but due to stabilizing, destabilizing and unclear links in the network of tipping elements, the change of the equilibrium states fluctuates more and shows a higher uncertainty. Therefore, the standard deviation increases. The same is valid for panel (A) at a region where (parts) the critical temperature ranges T limit,i are located for many of the tipping elements, that is mainly between 2 • C-4 • C above pre-industrial. elements start to show Kadyrov oscillations. Such Kadyrov oscillations have already been found in the early literature on dynamical systems of the CUSP type [92]. As shown in figure 7 for initial conditions at −1 for all tipping elements, this is most prominently the case for large interaction strengths and medium temperature increase values. Here, about every tenth solution is oscillating. This is due to the fact that uncertainties are largest in these regimes. For smaller interaction strength values, limit cycles can still occur but are much rarer with an occurrence at about 1% of all solutions. Of all these limit cycle oscillations almost all (95%) have a significant amplitude (standard deviation > 0.1) in at least one tipping element. The most common limit cycles are simultaneous oscillations of AMOC and GIS as shown in figure 7(D). They make up about 86% of all oscillating states found. The reason for this predominant oscillation is that there is a strong negative feedback loop between the GIS and the AMOC via freshwater input from Greenland that weakens the AMOC, while on the other side a weaker AMOC cools the northern hemisphere (see e.g. [16,29]). Still, whether such oscillations could indeed exist in the climate system remains speculative, but in principle there is evidence of oscillatory behavior in paleo data of the Earth system [106,107].

Discussion
We find that the only dominating stable state in the long term, for large temperature increases around 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 of the climate system regarding their increased basin volume in case more than one element is tipped. Based on the known interactions from Kriegler et al [29] this makes sense, since the interactions between the ice sheets, especially from Greenland to west Antarctica, are strong due 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 and reduce its overturning strength, but in our model these interactions are not sufficient to tip the AMOC over in many cases. These states with disintegrated ice sheets are especially relevant exhibiting a high basin volume for intermediate climate warming scenarios consistent with the climate target of the Paris agreement that aims at limiting global warming to well below 2 • C above pre-industrial levels [108]. Limit cycle oscillations between the tipping of some elements have been detected at some rare parameter configurations, mainly between the GIS and the AMOC. Although it remains unclear whether such (Kadyrov) oscillations have occurred in the climate system, they point to possible relevant internal modes of variability in the climate system. In principle such limit cycle behavior could have played a role in paleo climate dynamics such as in the Pleistocene ice age cycles [106,107]. Further, the individual dynamics are not the sole determinant of the final state of the tipping elements since the network effects can cause additional tipping events. Through this network interaction, it is therefore possible that cascades of tipping events emerge, even before the actual critical temperature threshold for some of the tipping elements is reached [35].

Conclusion
In this work, we study a conceptual model of five climate tipping elements based on a system of coupled, nonlinear differential equations. We investigate the stability of different dynamical regimes with respect to its stable states applying the concept of basin stability using a very large-scale Monte Carlo simulation of more than 3.5 billion ensemble members. Following that approach, we are able to propagate the numerous uncertainties thoroughly which are associated with the critical temperature thresholds and interaction strengths. With a Monte Carlo basin bifurcation 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 elements are in the transgressed state, especially for large levels of global warming. We also detect that the ice sheets are of special importance for the stability of states, where the large cryosphere components reside in the transgressed state, while the other tipping elements do not. We also detect Hopf-bifurcations for few parameter configurations (0.6%), mainly taking place between the GIS and the AMOC (86%).
Our complex dynamical networks approach strongly simplifies the nature of tipping elements as well as their interaction structure. However, it can serve to integrate simplified concepts of tipping elements until coupled, process-based models are developed that can resolve the respective nonlinearities in the Earth system in more detail since current state-of-the-art Earth system 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 studies have hypothesized that major changes in ENSO are possible [60,61] based on conceptual models [109,110], but however, whether this is evidence for a permanent and potentially even irreversibly tipped ENSO remains uncertain and debated. Surely, ENSO exerts strong feedbacks onto the climate system that will increase if major El-Niño events become more frequent, for instance through strong drying trends over Amazonia. Furthermore, in earlier research we found that the main results of our model remain robust under the omission of ENSO such that we decided to investigate the more complex case and included ENSO here, even though the use of equation (5) is only a topologically equivalent dynamical equation (for more details see [35]). While some literature studies present ENSO among the list of potential tipping elements [6,10,29], it still 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 as their interaction structure would be understood. It is also possible to probe the effect of different structural interaction hypotheses to further tipping elements within the scope of an uncertainty analysis, as has already been performed here for three interaction links. Further, the results of our study motivate that it could be worthwhile to look into the dynamics in more detail using process-detailed Earth system models. Especially the role of the large sheets in the stability landscape and oscillations between climate system components could be of interest. Even though, there is some knowledge about the interaction structure present in literature (see section 2.1), a new expert elicitation might be worthwhile because the knowledge about the interactions between the tipping elements has surely widened since the original expert elicitation from Kriegler et al (2009) [29].

Note on color maps
This paper makes use of the conceptually uniform colormaps developed by [111].  figure 1 has a specific link strength range and a specific physical process that is connected to the respective interaction. The link strength ranges are computed from literature values [29,43] such that they can be used in equation (5). For a more in depth description please be referred to   [35].  [10], see also equation (5). The typical tipping time scale τ i is given in model years (in arbitrary units) since it is beyond the scope of this model to make predictions about the exact tipping times. However, certain differences in tipping times as used here can be decisive whether a tipping event occurs or not. For more information see Wunderling et al (2020) [35].  Basin stability for single elements showing that the basin volume for the AMOC, WAIS and AMAZ are qualitatively similar since it is possible that only this particular element tips. However, for ENSO and the GIS this is not the case. It can only very rarely happen that these elements tip on themselves for high interaction strength even at low temperature increases since both of them possess a very strong link to another element that they would draw along into the tipped state. For ENSO, this is the AMAZ and for the GIS it is the AMOC. Note that the color bar is different for panel (A) to improve visibility. Oscillating states for random initial conditions. Limit cycles occur more often than for initial conditions at −1 for all tipping elements. On the other side, the limit cycle oscillation between the AMOC and the GIS is still dominating.