Interacting tipping elements increase risk of climate domino effects under global warming

With progressing global warming, there is an increased risk that one or several tipping elements in the climate system might cross a critical threshold, resulting in severe consequences for the global climate, ecosystems and human societies. While the underlying processes are fairly well-understood, it is unclear how their interactions might impact the overall stability of the Earth's climate system. As of yet, this cannot be fully analysed with state-of-the-art Earth system models due to computational constraints as well as some missing and uncertain process representations of certain tipping elements. Here, we explicitly study the effects of known physical interactions among the Greenland and West Antarctic ice sheets, the Atlantic Meridional Overturning Circulation (AMOC) and the Amazon rainforest using a conceptual network approach. We analyse the risk of domino effects being triggered by each of the individual tipping elements under global warming in equilibrium experiments. In these experiments, we propagate the uncertainties in critical temperature thresholds, interaction strengths and interaction structure via large ensembles of simulations in a Monte Carlo approach. Overall, we find that the interactions tend to destabilise the network of tipping elements. Furthermore, our analysis reveals the qualitative role of each of the four tipping elements within the network, showing that the polar ice sheets on Greenland and West Antarctica are oftentimes the initiators of tipping cascades, while the AMOC acts as a mediator transmitting cascades. This indicates that the ice sheets, which are already at risk of transgressing their temperature thresholds within the Paris range of 1.5 to 2 ^∘C, are of particular importance for the stability of the climate system as a whole.

1 Introduction

1.1 Tipping elements in the climate system

The Earth system comprises a number of large-scale subsystems, the so-called tipping elements, that can undergo large and possibly irreversible changes in response to environmental or anthropogenic perturbations once a certain critical threshold in forcing is exceeded (Lenton et al., 2008). Once triggered, the actual tipping process might take several years up to millennia depending on the respective response times of the system (Hughes et al., 2013; Lenton et al., 2008). Among the tipping elements are cryosphere components such as the continental ice sheets on Greenland and Antarctica; biosphere components such as the Amazon rainforest, boreal forests and coral reefs; and large-scale atmospheric and oceanic circulation patterns such as monsoon systems and the Atlantic Meridional Overturning Circulation (AMOC). With continuing global warming, it becomes more likely that critical thresholds of some tipping elements might be exceeded, possibly within this century, triggering severe consequences for ecosystems, infrastructure and human societies. These critical thresholds can be quantified with respect to the global mean surface temperature (GMT), resulting in three clusters of tipping elements that are characterised by their critical temperature between 1 and 3, between 3 and 5, and above 5 ^∘C of warming compared with pre-industrial temperatures respectively (Schellnhuber et al., 2016). The most vulnerable cluster, which is already at risk between 1 and 3 ^∘C of warming, includes several cryosphere components, specifically mountain glaciers as well as the Greenland and West Antarctic ice sheets. Recent studies suggest potential early-warning indicators for these tipping elements, showing that some of them are approaching or might have already transgressed a critical threshold (Lenton et al., 2019; Caesar et al., 2018; Nobre et al., 2016; Favier et al., 2014).

1.2 Interactions between climate tipping elements

The tipping elements in the Earth system are not isolated systems; they interact on a global scale (Lenton et al., 2019; Kriegler et al., 2009). These interactions could have stabilising or destabilising effects, increasing or decreasing the probability of emerging tipping cascades, and it remains an important problem to understand how the interactions between the tipping elements affect the overall stability of the Earth system. Despite the considerable recent progress in global Earth system modelling, current state-of-the-art Earth system models cannot yet comprehensively simulate the non-linear behaviour and feedbacks between some of the tipping elements due to computational limitations (Wood et al., 2019). Furthermore, the interactions between tipping elements have only partially been described in a framework of more conceptual (but process-based) models, and our current understanding of the interaction structure of tipping elements is partly based on expert knowledge. For a subset of five tipping elements, an expert elicitation was conducted that synthesised a causal interaction structure and an estimation for the probability of tipping cascades to emerge (Kriegler et al., 2009). These studied tipping elements were the Greenland Ice Sheet, the West Antarctic Ice Sheet, the Atlantic Meridional Overturning Circulation (AMOC), the El Niño–Southern Oscillation (ENSO) and the Amazon rainforest (see Figs. 1 and S3). Although this network is not complete with respect to the physical interactions between the tipping elements and the actual set of tipping elements themselves (Wang and Hausfather, 2020; Lenton et al., 2019; Steffen et al., 2018), it presented a first step towards synthesising the positive and negative feedbacks between climate tipping elements. To the best of our knowledge, a systematic update of this assessment or a comparably comprehensive expert assessment has not been undertaken since Kriegler et al. (2009).

Figure 1 Interactions between climate tipping elements and their roles in tipping cascades. The Greenland Ice Sheet, the West Antarctic Ice Sheet, the Atlantic Meridional Overturning Circulation (AMOC) and the Amazon rainforest are depicted along with their main interactions (Kriegler et al., 2009). The links between the tipping elements are colour coded, where red arrows depict destabilising and blue arrows depict stabilising interactions. Where the direction is unclear, the link is marked in grey. A more thorough description of each of the tipping elements and the links can be found in Tables 1 and 2 as well as in Sect. 2. Where tipping cascades arise, the relative size of the dominoes illustrates how many ensemble members the respective climate component initiates tipping cascades in (red domino) or how many tipping cascades the respective climate component occurs in (blue domino). Standard deviations for these values are given in Fig. S1a and b. Generally, the polar ice sheets are found to more frequently take on the role of cascade initiators than the AMOC and Amazon rainforest.

Based on the network from this expert elicitation and a Boolean approach founded on graph grammars, an earlier study found that the strong positive–negative feedback loop between the Greenland Ice Sheet and the AMOC might act as a stabiliser to the Earth system (Gaucherel and Moron, 2017). Also, using the interaction network data of Kriegler et al. (2009), it has been shown that large economic damages due to tipping cascades could arise with respect to the social cost of carbon, using a stochastic and dynamic evaluation of tipping points in an integrated assessment model (Cai et al., 2016). Other studies also quantified the economic impacts of single climate tipping events and tipping interactions (Lemoine and Traeger, 2016; Cai et al., 2015). In the light of recent studies that hypothesise a considerable risk of current anthropogenic pressures triggering tipping cascades, up to a potential global cascade (towards a so-called “hothouse state” of the Earth system) (Lenton et al., 2019; Steffen et al., 2018), we here aim at developing a conceptual dynamic network model that can assess whether interactions of tipping elements have an overall stabilising or destabilising effect on the global climate state. As such, we view our approach as a hypotheses generator that produces qualitative scenarios (rather than exact quantifications or projections) that can then be further examined by more process-detailed Earth system models. In this way, the results of this study can lay the foundations and possibly guide towards a more detailed analysis with more complex models or data-based approaches.

1.3 Constraints from current observations and palaeoclimatic evidence

Observations over the past decades show that several tipping elements are already impacted by progressing global warming (Wang and Hausfather, 2020; Lenton et al., 2019; IPCC, 2014; Levermann et al., 2010). Ice loss from Greenland and West Antarctica has increased and accelerated over the past decades (Shepherd et al., 2018; Khan et al., 2014; Zwally et al., 2011). Recent studies suggest that the Amundsen Basin in West Antarctica might in fact have already crossed a tipping point (Favier et al., 2014; Rignot et al., 2014). The grounding lines of glaciers in this region are rapidly retreating, which could induce local marine ice sheet instabilities and eventually lead to the disintegration of the entire basin (Mercer, 1978; Weertman, 1974). Palaeoclimate records suggest that parts of Antarctica and larger parts of Greenland might already have experienced strong ice retreat in the past, especially during the Pliocene as well as during Marine Isotope Stage 5e and 11 (Dutton et al., 2015).

It has also been shown that the AMOC experienced a significant slowdown since the mid 20th century (Caesar et al., 2018), which has led to the weakest AMOC state in centuries (Caesar et al., 2021). Models from Phase 5 of the Coupled Model Intercomparison Project (CMIP5) have shown that weakening of the AMOC is currently largely caused by changing surface buoyancy fluxes (Levang and Schmitt, 2020). However, in the future, the overturning strength of the AMOC might also be adversely impacted by increased freshwater forcing of the North Atlantic ocean by meltwater influx from Greenland (Bakker et al., 2016; Böning et al., 2016). An AMOC slowdown also likely occurred during the last deglaciation in the Heinrich event 1 and Younger Dryas cold periods, as proxies from sea surface and air temperatures as well as climate model simulations suggest (Ritz et al., 2013).

The Amazon rainforest is not only directly impacted by anthropogenic climate change, including the increased risk of extensive drought events or heat waves (Marengo and Espinoza, 2015; Brando et al., 2014), but also by deforestation and fire (Thonicke et al., 2020; Malhi et al., 2009). This increases the likelihood that parts of it will shift from a rainforest to a savannah state, for instance through diminished moisture recycling (Staal et al., 2018; Zemp et al., 2017). It is suspected that the Amazon rainforest could be close to a critical extent of deforestation which might, along with global warming, suffice to initiate such a critical transition (Nobre et al., 2016). This could put 30 %–50 % of rainforest ecosystems at risk of shifting the rainforest to tropical savannah or dry forests (Nobre et al., 2016). From a local to regional point of view, the potential for critical transitions in the rainforest is further examined by more recent studies (Staal et al., 2020; Ciemer et al., 2019).

1.4 Structure of this work

In Sect. 2, we provide an overview of the biogeophysical processes governing the individual dynamics and interactions of the four tipping elements considered here as well as how these are represented in our conceptual network model. We also describe the construction of the large-scale Monte Carlo ensemble which enables us to propagate the parameter uncertainties inherent in the modelled tipping elements and their interactions. In Sect. 3, we explore how the critical threshold temperature ranges of the tipping elements change with increasing overall interaction strength. We also show which tipping elements initiate and transmit tipping cascades, revealing the characteristic roles of the tipping elements in the Earth system. Moreover, we discuss the distinct nature of ENSO as a potential tipping element and present results of a robustness analysis including this additional tipping element in our network model. Section 4 summarises the results and discusses the limitations of our approach. It also outlines possible further lines of research concerning tipping element interactions and the risks of emerging tipping cascades with more process-detailed models.

Figure 2 Schematic overview of the generalised tipping element and time series of a tipping cascade. (a) Exemplary bifurcation diagram of a tipping element with two stable regimes: the lower state indicates the stable baseline regime, and the upper state indicates the stable transitioned regime. For the Greenland Ice Sheet, for instance, these correspond to its pre-industrial, almost completely ice-covered state (stable baseline regime) and an almost ice-free state (stable transitioned regime), as can be expected in the long term for higher warming scenarios (Robinson et al., 2012). There are two ways that a tipping element can transgress its critical threshold (unstable manifold) and move into the transitioned state: an increase in the global mean surface temperature or via interactions with other climate components. In both cases, the tipping element converges to the stable transitioned regime indicated by the red hollow circles. (b) Exemplary time series showing a tipping cascade of two elements. Here, Greenland transgresses its critical temperature (T_{limit,Greenland}) first, i.e. would become ice-free. Through its interaction with the West Antarctic Ice Sheet, the West Antarctic Ice Sheet then transgresses the unstable manifold in the vertical direction (following the path of the red upward-directed arrow in panel a). This example is based on a scenario with global mean surface temperature increase of 1.6 ^∘C above pre-industrial levels and an interaction strength d=0.16 (see also Fig. 3).

2 Methods

In the following, we present our dynamic network approach for modelling tipping interactions and cascades in the Earth system. In Sect. 2.1, we motivate the use of a stylised equation to represent climate tipping elements in a conceptual manner. This equation exhibits a double-fold bifurcation (see Fig. 2)

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}=\left[-x_i^{\mathrm{3}}+x_i+c_i\right]{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_i}}.\end{array}$$

Here, x_i indicates the state of a certain tipping element, c_i is the critical parameter and τ_i the typical tipping timescale with i = (Greenland Ice Sheet, West Antarctic Ice Sheet, AMOC, Amazon rainforest). This approach has already been used frequently for qualitatively describing tipping dynamics in different applications and network types and has been applied to systems in climate, ecology, economics and political science (Klose et al., 2020; Krönke et al., 2020a; Wunderling et al., 2020a; Dekker et al., 2018; Brummitt et al., 2015; Abraham et al., 1991).

To describe the tipping elements' interactions, we extend Eq. (1) by a linear coupling term (Klose et al., 2020; Krönke et al., 2020a; Brummitt et al., 2015) to yield

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}=\left[\stackrel{\mathrm{Individual}\mathrm{dynamics}\mathrm{term}}{\overbrace{-x_i^{\mathrm{3}}+x_i+c_i}}+\stackrel{\mathrm{Coupling}\mathrm{term}}{\overbrace{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{\begin{array}{c}j\\ j\ne i\end{array}}{d}_{ij}\left(x_j+\mathrm{1}\right)}}\right]{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_i}},\end{array}$$

and we describe the physical interpretation of these interactions between the tipping elements in Sect. 2.2. While the first term (individual dynamics term) determines the dynamical properties of each individual tipping element, the second term (coupling term) describes the effects of interactions between tipping elements. If the prefactors in front of the cubic and the linear term are unity as in Eq. (2) and the additive coupling term is neglected (d_ij=0 for all i, j), the critical threshold values where qualitative state changes occur are $c_i^{\mathrm{1},\mathrm{2}}=\pm \sqrt{\mathrm{4}/\mathrm{27}}$ (Klose et al., 2020). The system described by this differential equation is bistable for values of the critical parameter between c₁ and c₂ and can here be separated into a transitioned and a baseline state, where $x_i=-\mathrm{1}$ denotes the baseline state and $x_i=+\mathrm{1}$ the completely transitioned one (see Fig. 2).

Building on these model equations, we describe the fully parameterised model and its parameters as it is used in this study in Sect. 2.3. Specifications of how tipping cascades are evaluated and timescales are chosen can be found in Sect. 2.4 and 2.5. Lastly, our large-scale Monte Carlo ensemble approach for the propagation of parameter and interaction network uncertainties is described in Sect. 2.6.

2.1 From conceptual to process-detailed models of climate tipping elements

In the conceptual network model investigated in this study, the main dynamics of each of the tipping elements are condensed to a non-linear differential equation with two stable states representing the current (baseline) state and a possible transitioned state capturing the qualitative dynamics of generalised tipping elements (see Eq. 1). This serves as a stylised representation of the Greenland Ice Sheet, the West Antarctic Ice Sheet, the AMOC and the Amazon rainforest. Here, we focus on these four out of a larger range of tipping elements in the cryosphere, biosphere, and oceanic and atmospheric circulation patterns that have been suggested in the literature (Schellnhuber et al., 2016; Scheffer et al., 2009; Lenton et al., 2008). In this study, we do not consider possible “back-tipping” (or hysteresis behaviour) of climate tipping elements, as the forcing represented by global mean surface temperature anomalies is only increased, never decreased, in our experiments. It is clear that the representation of a complex climate tipping element with all its interacting processes as well as positive and negative feedbacks in a stylised cusp bifurcation model is a strong simplification. In the following, we elaborate on why such a cusp bifurcation structure (Eq. 1) can nonetheless be assumed to capture the overall stability behaviour of these four tipping elements (Bathiany et al., 2016) before we introduce more mathematical details of our dynamical systems approach in Sect. 2.3.

2.1.1 AMOC

Early conceptual models introduced in the 1960s showed that the AMOC could exhibit a cusp-like behaviour, using simplified box models based on the so-called salt–advection feedback (Stommel, 1961; Cessi, 1994). Many extensions and updates to this well-known box model approach have been put forward, each confirming the potential multi-stability of the AMOC (e.g. Wood et al., 2019). More complex Earth system models including a fully implicit ocean model (Huisman et al., 2010), Earth system models of intermediate complexity (EMICs, e.g. CLIMBER; Rahmstorf et al., 2005) and an atmosphere–ocean general circulation model (AOGCM, e.g. the FAMOUS model; Hawkins et al., 2011) have shown hysteresis behaviour which is qualitatively similar to Eq. (1). Furthermore, palaeoclimatic evidence suggests a bistability of the AMOC: in palaeoclimate records, Dansgaard–Oeschger events (see e.g. Crucifix, 2012) have been associated with large reorganisations of the AMOC (Ditlevsen et al., 2005; Timmermann et al., 2003; Ganopolski and Rahmstorf, 2002), where ice core data link the events to sea surface temperature increases in the North Atlantic. Even though there are considerable uncertainties, estimates from the literature suggest that the level of global warming sufficient for tipping the AMOC is between 3.5 and 6.0 ^∘C (Schellnhuber et al., 2016; Lenton, 2012; Levermann et al., 2012; Lenton et al., 2008), with the risk of crossing a critical threshold considerably increasing beyond 4 ^∘C above pre-industrial temperature levels (Kriegler et al., 2009).

2.1.2 Greenland Ice Sheet

Previous studies have shown that a double fold-like bifurcation structure for the ice sheets can arise from the melt–elevation feedback (Levermann and Winkelmann, 2016) as well as from the marine ice sheet instability and other positive feedback mechanisms (e.g. DeConto and Pollard, 2016; Schoof, 2007). In particular, dynamic ice sheet model simulations have identified irreversible ice loss once a critical temperature threshold is crossed (Toniazzo et al., 2004), leading to multiple stable states and hysteresis behaviour for the Greenland Ice Sheet (Robinson et al., 2012; Ridley et al., 2010). In Robinson et al. (2012), the critical temperature range for an irreversible disintegration of the Greenland Ice Sheet has been estimated between 0.8 and 3.2 ^∘C of warming above pre-industrial global mean surface temperature levels. Palaeoclimate evidence further suggests that there have been substantial, potentially self-sustained retreats of the Greenland Ice Sheet in the past. It has, for instance, been simulated that the Greenland Ice Sheet can disintegrate if warmer ocean conditions from the Pliocene are applied to an initially glaciated Greenland (Koenig et al., 2014). Further, Greenland was nearly ice-free for extended interglacial periods during the Pleistocene (Schaefer et al., 2016). Sea level reconstructions further suggest that large parts of Greenland could have been ice-free during Marine Isotope Stage 11 and the Pliocene (Dutton et al., 2015).

2.1.3 West Antarctic Ice Sheet

Compared with the case of the Greenland Ice Sheet, different processes make the West Antarctic Ice Sheet susceptible to tipping dynamics. As large parts of West Antarctica are grounded in marine basins, changes in the ocean are key in driving the evolution of the ice sheet. The marine ice sheet instability can trigger self-sustained ice loss where the ice sheet is resting below sea level on retrograde sloping bedrock (Weertman, 1974; Schoof, 2007). This destabilising mechanism is possibly already underway in the Amundsen Sea region (Favier et al., 2014; Joughin et al., 2014; Joughin and Alley, 2011). Once triggered, a single local perturbation via increased sub-shelf melting in the Amundsen region could lead to wide-spread retreat of the West Antarctic Ice Sheet (Feldmann and Levermann, 2015). Further, a recent study has shown strong hysteresis behaviour for the whole Antarctic Ice Sheet, identifying two major thresholds that lead to a destabilisation of West Antarctica around 2 ^∘C of global warming and large parts of East Antarctica between 6 and 9 ^∘C of global warming (Garbe et al., 2020). It is likely that the West Antarctic Ice Sheet has experienced brief but dramatic retreats during the past 5 Myr (Pollard and DeConto, 2009). Prior collapses have been suggested from deep-sea-core isotopes and sea level records (Gasson, 2016; Dutton et al., 2015; Pollard and DeConto, 2005).

2.1.4 Amazon rainforest

Conceptual models of the Amazon have identified multi-stability between rainforest, savannah and treeless states, leading to hysteresis (Staal et al., 2016, 2015; Van Nes et al., 2014). This hysteresis has been found to be shaped by local-scale tipping points of the Amazon rainforest, and its resilience might be diminished under climate change until the end of the 21st century (Staal et al., 2020). More complex dynamic vegetation models also found alternative stable states of the Amazon ecosystem (Oyama and Nobre, 2003) and suggest that rainforest dieback might be possible due to drying of the Amazon Basin under future climate change scenarios (Nobre et al., 2016; Cox et al., 2004, 2000). Observational data further support the potential for multi-stability of the Amazon rainforest (Ciemer et al., 2019; Hirota et al., 2011; Staver et al., 2011). While it remains an open question whether the Amazon has a single system-wide tipping point, the projected increase in droughts and fires (Malhi et al., 2009; Cox et al., 2008) is likely to impact the forest cover on a local to regional scale, which might spread to other parts of the region via moisture-recycling feedbacks (Zemp et al., 2017, 2014; Aragão, 2012). It is important to note that in contrast to the ice sheets and ocean circulation, the rainforest is able to adapt to changing climate conditions to a certain extent (Sakschewski et al., 2016). However, this adaptive capacity might still be outpaced if climate change progresses too rapidly (Wunderling et al., 2020c). A dieback of the Amazon rainforest has been found under a business-as-usual emissions scenario (Cox et al., 2004), which would be equivalent to a global warming of more than 3 ^∘C above pre-industrial levels until 2100 (≈3.5–4.5 ^∘C (see also Schellnhuber et al., 2016)), mainly due to more persistent El Niño conditions (Betts et al., 2004).

2.2 Physical interpretation of tipping element interactions

Based on these conceptual models as well as building on first coupled experiments with a discrete state Boolean model (Gaucherel and Moron, 2017) and economic impact studies (Cai et al., 2016; Lemoine and Traeger, 2016; Cai et al., 2015), we describe the interactions of the four tipping elements in a network approach using a set of linearly coupled, topologically equivalent differential equations (Kuznetsov, 2004). In the following, we go through the different main interactions of the four tipping elements considered here and expand on the underlying physical processes. Overall, the additional literature supports and refines the results from an early expert elicitation (Kriegler et al., 2009).

2.2.1 Greenland Ice Sheet → AMOC

Increasing freshwater input from enhanced melting of the Greenland Ice Sheet can lead to a weakening of the AMOC, as supported by palaeoclimate evidence as well as modelling studies (Driesschaert et al., 2007; Jungclaus et al., 2006; Rahmstorf et al., 2005). Palaeoclimatic records further suggest that the AMOC could exist in multiple stable states, based on observed temperature changes associated with meltwater influx into the North Atlantic (Blunier and Brook, 2001; Dansgaard et al., 1993). Therefore, it is likely that a tipping of the Greenland Ice Sheet would lead to a destabilisation of the AMOC (see Fig. 1).

2.2.2 AMOC → Greenland Ice Sheet

Conversely, if the AMOC weakens, leading to a decline in its northward surface heat transport, Greenland might experience cooler temperatures (e.g. Jackson et al., 2015; Timmermann et al., 2007; Stouffer et al., 2006), which would have a stabilising effect on the ice sheet. With version 3 of the Hadley Centre Global Environment Model (HadGEM3), it has been shown that temperatures in Europe could drop by several degrees if the AMOC collapses, regionally up to 8 ^∘C (Jackson et al., 2015). A cooling trend in sea surface temperatures (SSTs) over the subpolar gyre, as a result of a weakening AMOC, has been confirmed by recent reanalysis and observation data (Caesar et al., 2018; Jackson et al., 2016; Frajka‐Williams, 2015; Robson et al., 2014). This “fingerprint” translates a reduction in overturning strength of 1.7 Sv per century to 0.44 K SST cooling per century (Caesar et al., 2018). AMOC regime shifts between weaker and stronger overturning strength during the last glacial period have been associated with large regional temperature changes in Greenland, for example during Dansgaard–Oeschger or Heinrich events (Barker and Knorr, 2016). Moreover, there is palaeoclimatic evidence from 3.6 Myr ago that a weaker North Atlantic current as part of the AMOC fostered Arctic sea ice growth which might have preceded continental glaciation in the Northern Hemisphere at that time (Karas et al., 2020). Based on these findings, we assume that a weakening of the AMOC would have a stabilising effect on the Greenland Ice Sheet (see Fig. 1).

2.2.3 West Antarctic Ice Sheet → AMOC

It remains unclear whether increased ice loss from the West Antarctic Ice Sheet has a stabilising or destabilising effect on the AMOC (see Fig. 1). Swingedouw et al. (2009) identified different processes based on freshwater hosing experiments into the Southern Ocean, which could be associated with a melting West Antarctic Ice Sheet (Swingedouw et al., 2009). Using the LOVECLIM1.1 EMIC, the authors found both enhancing and weakening effects on the AMOC strength. First, deep-water adjustments are observed. This means that an increase in the North Atlantic Deep Water formation is observed in response to a decrease in Antarctic Bottom Water production due to the release of freshwater in the Southern Ocean. This mechanism has been termed the so-called “bipolar ocean see-saw”. Second, salinity anomalies in the Southern Ocean are distributed to the North Atlantic, which dampens the North Atlantic Deep Water formation (compare to Seidov et al., 2005). Third, the North Atlantic Deep Water formation is enhanced by a strengthening of Southern hemispheric winds in response to a Southern hemispheric cooling. The reason for the stronger winds is the greater meridional temperature gradient between a cooler Antarctic region (due to the hosing experiment) and the Equator. This effect has been termed the “Drake Passage effect” (Toggweiler and Samuels, 1995).

Overall, the first and the third mechanism tend to strengthen the AMOC, whereas the second process would rather lead to a weakening of the AMOC. The specific timescales and relative strengths of these mechanisms are as of yet unclear (Swingedouw et al., 2009). In a coupled ocean–atmosphere model, a slight weakening of the AMOC was detected for a freshwater input of 1.0 Sv in the Southern Ocean over 100 years (Seidov et al., 2005). However, other studies suggest that the overturning strength of the AMOC remains at a high level if influenced by freshwater input from the West Antarctic Ice Sheet due to the effects from the bipolar ocean see-saw by decreasing Antarctic Bottom Water formation as described above (Swingedouw et al., 2008).

2.2.4 AMOC → West Antarctic Ice Sheet

The interaction from the AMOC to the West Antarctic Ice Sheet is destabilising (see Fig. 1). If the AMOC shut down, sea surface temperature anomalies could appear as the northward heat transport is diminished significantly. This could then lead to a warmer South and colder North, as observed in modelling studies (Weijer et al., 2019; Timmermann et al., 2007; Stouffer et al., 2006; Vellinga and Wood, 2002). A model intercomparison study for EMICs and AOGCMs found a sharp decrease in surface air temperatures over the Northern Hemisphere, whereas a slight increase over the Southern Hemisphere and around the Antarctic Ice Sheet was observed (Stouffer et al., 2006). In their study (Stouffer et al., 2006), a forcing of 1.0 Sv was applied to the northern part of the North Atlantic Ocean. Therefore, we set this link as destabilising in the interaction network model (see Fig. 1).

2.2.5 Greenland Ice Sheet ↔ West Antarctic Ice Sheet

The direct interaction between the Greenland and the West Antarctic ice sheets via sea level changes can be regarded as mutually destabilising, although with different magnitudes (see Fig. 1). It is a well-known phenomenon from tidal changes that grounding lines of ice sheets are varying (e.g. Sayag and Worster, 2013). Therefore, the Greenland Ice Sheet and the West Antarctic Ice Sheet could influence each other by sea level rise if one or the other cryospheric element would melt. Gravitational as well as elastic and rotational impacts would then enhance the sea level rise if one of the huge ice sheets melted first, as then only the other ice sheets would exert strong gravitational forces on ocean waters (Kopp et al., 2010; Mitrovica et al., 2009). The impact of this effect would be larger if Greenland became ice-free earlier than West Antarctica, as many marine-terminating ice shelves are located in West Antarctica, but the interaction is destabilising in both directions (see Fig. 1).

2.2.6 AMOC → Amazon rainforest

Lastly, the interaction between the AMOC and the Amazon rainforest is set as unclear (see Fig. 1). It is suspected that the Intertropical Convergence Zone (ITCZ) would be shifted southward if the AMOC collapsed. This could cause large changes in seasonal precipitation on a local scale and could, as such, have strong impacts on the Amazon rainforest (Jackson et al., 2015; Parsons, 2014). In the ESM2M Earth system model, it has been found that a strongly suppressed AMOC, through a 1.0 Sv freshwater forcing, leads to drying over many regions of the Amazon rainforest (Parsons, 2014). However, some regions would receive more rainfall than before. On a seasonal level, the wet season precipitation is diminished strongly, whereas the dry season precipitation is significantly increased (Jackson et al., 2015; Parsons, 2014). This could have consequences for the current vegetation that is adapted to this seasonal precipitation, especially in places where the seasonality is strong. However, overall, it remains unclear whether a tipped AMOC would have a reducing or increasing influence on the precipitation in South America. Instead, it might differ from location to location and is set as unclear in our study (see Fig. 1).

2.3 Dynamic network model of interacting tipping elements

In this subsection, we describe the details of the employed dynamic network model, the foundations of which are given by Eqs. (1) and (2). The critical parameter c_i of tipping element i is modelled as a function of the global mean surface temperature, i.e. $c_i=\sqrt{\frac{\mathrm{4}}{\mathrm{27}}}\cdot \frac{\mathrm{\Delta }\mathrm{GMT}}{T_{\mathrm{limit},i}}$, where T_limit,i is the critical temperature, and ΔGMT is the increase in the global mean surface temperature above pre-industrial levels. Note that ΔGMT denotes the global mean temperature at the surface, which should not be confused with a volume average temperature. This parameterisation implies that a state change is initiated as soon as the increase in GMT exceeds the critical temperature ($\frac{\mathrm{\Delta }\mathrm{GMT}}{T_{\mathrm{limit},i}}>\mathrm{1}$). In Table 1, the limits ΔT_limit,i are noted between which the critical temperature T_limit,i is uniformly drawn at random (see Sect. 2.6). In addition, we model the physical interactions between the tipping elements as a linear coupling (first-order approach). The coupling term $\frac{\mathrm{1}}{\mathrm{2}}\sum _j{d}_{ij}(x_j+\mathrm{1})$ consists of a sum of linear couplings to other elements x_j with $d_{ij}=d\cdot s_{ij}/\mathrm{5}$. It is necessary to add +1 to x_j such that the direction (sign) of coupling is only determined by d_ij and not by the state x_j. Thus, Eq. (2) becomes

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}=\left[-x_i^{\mathrm{3}}+x_i+\sqrt{{\displaystyle \frac{\mathrm{4}}{\mathrm{27}}}}\cdot {\displaystyle \frac{\mathrm{\Delta }\mathrm{GMT}}{T_{\mathrm{limit},i}}}+d\cdot \sum _{\begin{array}{c}j\\ j\ne i\end{array}}{\displaystyle \frac{s_{ij}}{\mathrm{10}}}\left(x_j+\mathrm{1}\right)\right]{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_i}}.\end{array}$$

Here, d is the overall interaction strength parameter that we vary in our simulations, and s_ij is the link strength based on the expert elicitation (Kriegler et al., 2009) (see Table 2 and Sect. 2.6). The prefactor $\frac{\mathrm{1}}{\mathrm{10}}$ sets the coupling term of Eq. 3 to the same scale as the individual dynamics term by normalising s_ij when d is varied between 0.0 and 1.0. The geophysical processes behind the interactions between the tipping elements are listed in Table 2 and are described and referenced in Sect. 2.2.

Table 1 Nodes in the modelled network of interacting tipping elements. For each tipping element in the network (see Fig. 1), a range of critical temperatures ΔT_limit is known from a review of the literature (Schellnhuber et al., 2016). Within this temperature range, the tipping element is likely to undergo a qualitative state transition.

In this network of tipping elements, very strong interactions exist, as detailed above. For each tipping element, there are two potential reasons for a state transition, either through the increase in GMT or through the coupling to other tipping elements (Fig. 2a).

The overall interaction strength d is described as a dimensionless parameter (see Eq. 3) that is varied over a wide range in our simulations, i.e. for d∈[0; 1], to account for the uncertainties in the actual physical interaction strength between the tipping elements. This way a range of different scenarios can be investigated. An interaction strength of zero implies no coupling between the elements such that only the individual dynamics remain. When the interaction strength reaches high values of around one, the coupling term is of the same order of magnitude as the individual dynamics term. In principle, more complex and data- or model-based interaction terms could be developed. However, while some interactions (e.g. between Greenland Ice Sheet and AMOC) have been established with EMICs such as CLIMBER-2 and LOVECLIM as well as global circulation models (GCMs) (Wood et al., 2019; Sterl et al., 2008; Driesschaert et al., 2007; Jungclaus et al., 2006; Rahmstorf et al., 2005), other interactions are less well-understood, potentially leading to biased coupling strengths (see also Sect. 2.2). Due to the sparsity of data concerning tipping interactions in the past, it remains challenging to extract the interaction parameters from palaeoclimatic evidence. Therefore, we attempt to include the full uncertainty ranges concerning the different model parameters and interaction strengths. To this end, we run large ensembles of simulations over long timescales. This is important as the disintegration of the ice sheets, for instance, would play out over thousands of years (Winkelmann et al., 2015; Robinson et al., 2012). Due to computational constraints, studying such an ensemble of millennial-scale simulations is typically not feasible with more complex Earth system models.

We propagate the considerable uncertainties linked to the parameters of the tipping elements and their interactions with a large-scale Monte Carlo approach (see Sect. 2.6).

2.4 Parameterisation of the tipping elements' intrinsic timescales

The four tipping elements in the coupled system of differential equations form a so-called “fast–slow system” (Kuehn, 2011), describing a dynamical system with slowly varying parameters compared with fast changing states x_i. We include the typical transition times τ_i from the baseline to the transitioned state in Eq. (3) based on values from the literature (Lenton et al., 2008; Robinson et al., 2012; Winkelmann et al., 2015), setting the tipping timescales for the Greenland Ice Sheet, West Antarctic Ice Sheet, AMOC and the Amazon rainforest to 4900, 2400, 300 and 50 years for a reference warming of 4 ^∘C above pre-industrial GMT respectively. The tipping timescale is calibrated at this reference temperature in the case of vanishing interaction between the elements. After calibration, the tipping time is allowed to scale freely with changes in the GMT and the interaction strength d.

We integrate all model simulations to equilibrium, such that the simulation time is at least 20 times larger than the longest assumed tipping timescale of 4900 years. As the actual absolute tipping times derived from our model simulations are difficult to interpret, our results should not be taken as a projection of how long potential tipping cascades would take to unfold. Rather, following our conceptual approach, we are interested in the relative differences (not the absolute values) between the typical tipping times, as they can be decisive as to whether a cascade emerges or not. Therefore, the figures below show the time in arbitrary units (a.u.; see Figs. 2 and 3).

Figure 3 Time series of tipping cascades. Exemplary time series of states for each of the four investigated tipping elements, simulated here until equilibrium is reached. For comparability reasons, the parameter settings for the time series are the same (exact parameters can be found in Table S1), and all time series are computed for ΔGMT increases of 1.4, 1.5, 1.6 and 1.7 ^∘C above pre-industrial (columns). Couplings are constant for each row. Tipping cascades as shown here are defined as the number of transitioned elements at a fixed interaction strength and ΔGMT compared to the simulation with a slightly higher ΔGMT (ΔGMT increase by 0.1 ^∘C) but same interaction strength. If, between these two simulations, some of the tipping elements alter their equilibrium state, a tipping cascade of the respective size occurred and is counted as such. (a) A single tipping event for an interaction strength of 0.08. Tipping occurs at 1.6 ^∘C. (b) A tipping cascade with two elements for an interaction strength of 0.16. (c) A tipping cascade with three elements for an interaction strength of 0.24. For other initial conditions, interaction strengths and global mean surface temperatures (ΔGMT) tipping cascades with four elements can occur too. Additionally, we marked the baseline and the transitioned regime as grey hatched areas. Between the hatched areas, the state is not stable, and a critical state transition occurs. In the lower grey area, the element is said to be in the baseline regime, and in the upper grey region, the element is said to be in the transitioned regime.

Table 2 Interaction links in the network of tipping elements. For each link in the network of Fig. 1, there is a strength and a sign for each interaction of the tipping elements. The sign indicates if the interaction between the tipping elements is increasing or decreasing the danger of tipping cascades. Following Kriegler et al. (2009), the strength s_ij gives an estimate in terms of increased or decreased probability of cascading transitions (Kriegler et al., 2009). For example, if Greenland transgresses its threshold, the probability that the AMOC does as well is increased by a factor of 10 (see entry for Greenland → AMOC). A random number between +1 and $s_{\mathrm{Greenland}\to \mathrm{AMOC}}=+\mathrm{10}$ is then drawn for our simulations and used for s_ij in Eq. (3). Conversely, the probability that Greenland transgresses its threshold if the AMOC is in the transitioned state is decreased by a factor of $\frac{\mathrm{1}}{\mathrm{10}}$. A random number between −1 and $s_{\mathrm{AMOC}\to \mathrm{Greenland}}=-\mathrm{10}$ is then drawn. The main physical processes that connect pairs of tipping elements are described in this table and in Sect. 2.2. The link strengths are grouped into strong, intermediate and weak links. Note that in the expert elicitation (Kriegler et al., 2009), there has been an estimation of the maximum increase or decrease in the tipping probability in the case that the element that starts the interaction is already in the transitioned state. For example, the link between Greenland and AMOC is given as [1; 10] in Kriegler et al. (2009) and is modelled here as a randomly drawn variable between 1 and 10 for s_ij. An example of an unclear coupling would be the link between West Antarctica and the AMOC which is given as [0.3; 3] in Kriegler et al. (2009) and we translate into an s_ij between −3 and 3. In general, the values are drawn between 1 and the respective maximum value s_ij if the interaction between i and j is positive or between −1 and the negative maximum value s_ij if the interaction between i and j is negative.

2.5 Modelling protocol and evaluation of tipping cascades

In our network model, if the critical temperature threshold of a tipping element is surpassed, it transgresses into the transitioned state (Fig. 2a) and can potentially increase the likelihood of further tipping events via its interactions: for instance, the increased freshwater influx from a disintegration of the Greenland Ice Sheet can induce a weakening or even collapse of the AMOC (Fig. 2b). In our model simulations, we consider increases in the global mean surface temperature from 0 up to 8 ^∘C above the pre-industrial average, which could be reached in worst-case scenarios such as the extended Representative Concentration Pathway 8.5 (RCP8.5) by the year 2500 (Schellnhuber et al., 2016; IPCC, 2014).

For each tipping element, we start from the baseline (non-tipped) state (where x_i is negative). Global warming or interactions with the other parts of the climate system can then cause the element to tip into the transitioned state (see Fig. 2). When the critical parameter reaches $\sqrt{\frac{\mathrm{4}}{\mathrm{27}}}$ from below (i.e. when ΔGMT reaches T_limit,i), the stable baseline state x_i reaches $-\frac{\mathrm{1}}{\sqrt{\mathrm{3}}}$ in the case of an autonomous tipping element. Therefore, the threshold for the baseline state is defined as $x_i^{-}=-\frac{\mathrm{1}}{\sqrt{\mathrm{3}}}$. If the critical parameter increases above $\sqrt{\frac{\mathrm{4}}{\mathrm{27}}}$, the state x_i is larger than $x_i^{-}$, the stability of the lower stable state is lost and a state transition towards the upper stable $x_i^{+}$ occurs. As for the lower stable state $x_i^{-}$, the stable transitioned state is defined for states $x_i>x_i^{+}=+\frac{\mathrm{1}}{\sqrt{\mathrm{3}}}$.

We identify and define tipping cascades at a fixed interaction strength d and GMT as the number of additionally tipped elements in equilibrium (as defined above) after an incremental GMT increase of 0.1 ^∘C. The tipping element with a critical temperature threshold closest to the GMT at this point is counted as the initiator of the cascade. All tipping elements that appear in a particular cascade are counted as an occurring tipping element in that tipping cascade.

With increasing global mean surface temperature and interaction strength, generally more tipping cascades occur (Fig. 3). However, the size, the timing and the occurrence of cascades can also depend critically on the specific initial conditions (Wunderling et al., 2020b), which are not varied in the experiments presented here. In an exemplary simulation, we show how a global mean surface temperature increase from 1.5 to 1.6 ^∘C triggers the Greenland Ice Sheet to transition to an ice-free state in one realisation of our Monte Carlo ensemble at low interaction strength (Fig. 3a). For larger interactions strengths, the West Antarctic Ice Sheet as well as AMOC might then also tip as part of a tipping cascade that was initiated by the Greenland Ice Sheet in this case (Fig. 3b, c). The initial conditions and parameters for the specific example in Fig. 3 can be found in Table S1 in the Supplement.

2.6 Monte Carlo sampling and propagation of uncertainties

As the strength of interactions between the tipping elements is highly uncertain, a dimensionless interaction strength is varied over a wide range in our network approach to cover a multitude of possible scenarios. To cope with the uncertainties in the critical threshold temperatures and in the link strengths between pairs of tipping elements (see Eq. 3 and Tables 1 and 2), we set up a Monte Carlo ensemble with approximately 3.7 million members in total.

This Monte Carlo ensemble is generated as follows: for each combination of global mean surface temperature ΔGMT and overall interaction strength d, we create 100 realisations from a continuous uniform distribution of randomly drawn parameter sets for critical threshold temperatures T_limit,i and interaction link strengths s_ij based on the uncertainty ranges given above (see Tables 1 and 2). As our model has 11 parameters with uncertainties (4 critical threshold temperature parameters and 7 interaction link strength parameters), we use a Latin hypercube sampling to construct a set of parameters for each ensemble members such that the multidimensional space of sampled parameters is more successfully covered than with a usual random sample generation (Baudin, 2013).

We also sample all nine different interaction network structures that arise when we permute all possibilities (negative, zero, positive) arising from the two unclear links between the AMOC and the Amazon rainforest, and between West Antarctica and the AMOC (see Table 2 and Fig. 1). For each of these nine network structures, we compute the same 100 starting conditions that we received from our Latin hypercube sampling. Thus, in total, we compute 900 samples for each GMT (0.0–8.0 ^∘C, step width of 0.1 ^∘C) and interaction strength (0.0–1.0, step width of 0.02) combination, resulting in a large ensemble of 3.7 million members overall.

Our approach is conservative in the sense that there are several destabilising interactions that are not considered here (Lenton et al., 2019; Steffen et al., 2018). Further, by sampling uncertain parameters from a uniform distribution, we are treating lower and higher threshold temperatures as well as strong and weak link interactions equally, potentially resulting in a more balanced ensemble. Additional knowledge about the critical threshold temperatures and interaction link strengths would considerably improve our analysis.

3 Results

3.1 Shift in effective critical threshold temperatures due to interactions

Owing to the interactions between the tipping elements, their respective critical temperatures (previously identified for each element individually, see Fig. 4a) are effectively shifted to lower values (except for Greenland, see Fig. 4b and c). For West Antarctica and the AMOC, we find a sharp decline for interaction strengths up to 0.2 and an approximately constant critical temperature range afterwards. The effective critical temperature for the Amazon is only marginally reduced due to the interactions within the network, as it is only influenced by the AMOC via an unclear link.

Figure 4 Shift of critical temperature ranges due to interactions. (a) Critical global mean surface temperatures for each of the four investigated tipping elements without taking interactions into account (as reproduced from the literature; Schellnhuber et al., 2016). Grey bars indicate the standard deviation arising when drawing from a random uniform distribution between the respective upper and lower temperature limits. These bars correspond to the critical temperature ranges in case of zero interaction strength in panels (b) and (c). Change in critical temperature ranges with increasing interaction strength for the Greenland Ice Sheet and West Antarctic Ice Sheet (b) and the Atlantic Meridional Overturning Circulation (AMOC) and the Amazon rainforest (c). The standard deviation of the critical temperatures for each tipping element within the Monte Carlo ensemble is given as respective coloured shading.

In particular, the ensemble average of the critical temperature at an interaction strength of d=1.0 is lowered by about 1.2 ^∘C (≈40 %) for the West Antarctic Ice Sheet, 2.75 ^∘C (≈55 %) for the AMOC and 0.5 ^∘C (≈10 %) for the Amazon rainforest respectively (see Fig. S2). This is due to the predominantly positive links between these tipping elements (see Fig. 1).

In contrast, the critical temperature range for the Greenland Ice Sheet tends to be raised due to the interaction with the other tipping elements, accompanied by a significant increase in overall uncertainty. This can be explained by the strong negative feedback loop between Greenland and the AMOC that is embedded in the assumed interaction network (see Table 2, see also Gaucherel and Moron, 2017). On the one hand, enhanced meltwater influx into the North Atlantic might dampen the AMOC (positive interaction link); on the other hand, a weakened overturning circulation would lead to a net-cooling effect around Greenland (negative interaction link). Thus, the state of Greenland strongly depends on the specific parameter values in critical threshold temperature and interaction link strength of the respective Monte Carlo ensemble members.

Overall, the interactions are more likely to lead to a destabilisation within the network of climate tipping elements with the exception of the Greenland Ice Sheet.

Table 3 Share of tipping events in ensemble simulations. For different maximum values of the interaction strength d (first column), the share of ensemble simulations is shown that have a tipping event or cascade (third column) within the Paris limit until the global mean surface temperature increase reaches 2.0 ^∘C above pre-industrial. This means that 61 % of all ensemble members contain a tipping event or cascade, whereas 39 % do not (second column) if all interaction strengths until 1.0 are considered (see Fig. 5a and b). Overall, the fraction of tipping events stays the same and does not decrease for lower maximum interaction strengths. However, the distribution of tipping events and cascade sizes changes, i.e. the number of large cascades decreases with lower maximal interaction strength. This is shown in the split last column that displays the share of cascades with one, two, three or four elements.

Figure 5 Tipping cascades for all interaction strengths between 0.0 and 1.0. (a, b) For global warming up to 2.0 ^∘C above pre-industrial levels, the coloured shading illustrates the fraction of model representations in the Monte Carlo ensemble without tipping events (grey), with a single tipping event (purple), and with cascades including two (red), three (dark orange) or four (light orange) elements. (c–e) Occurrence of tipping cascades with two, three or four elements as a function of global mean surface temperature increase. The counts are normalised to the highest value of the most frequent tipping cascade (in cascades with two elements). (f) Dominant cascades with two and three elements for temperature increases from 0 to 8 ^∘C above pre-industrial. Other cascades are not shown, as their relative occurrence is comparatively much smaller. The standard deviation represents the difference between the possible ensemble realisations of the interaction network (see Sect. 2.3). Hence, it tends to be larger for cascades where unclear interaction links are involved, e.g. for the AMOC–Amazon rainforest cascade (compare Fig. 1 and Table 2).

3.2 Risk of emerging tipping cascades

Tipping cascades occur when two or more tipping elements transgress their critical thresholds for a given temperature level (see Sect. 2.5). We evaluate the associated risk as the share of ensemble simulations in which such tipping cascades are detected. For global warming up to 2.0 ^∘C, tipping occurs in 61 % of all simulations (Fig. 5a). This comprises the tipping of individual elements (22 %) as well as cascades including two elements (21 %), three elements (15 %) and four elements (3 %; see Fig. 5b). As the coupling between the tipping elements is highly uncertain, we introduce an upper limit to the maximum interaction strength and vary it from 0.0 to 1.0 (see Table 3). The highest value of 1.0 implies that the interaction between the elements is as important as the non-linear threshold behaviour of an individual element (see Eq. 3). For lower values, the interaction plays a less dominant role. We find that the occurrence of tipping events does not depend significantly on the maximum interaction strength; however, the cascade size decreases for lower values.

Tipping cascades are first induced at warming levels around 1 ^∘C above pre-industrial GMT, where the lower bound of the critical temperature range for the Greenland Ice Sheet is exceeded. The bulk of tipping cascades, however, is found between a 1 and 3 ^∘C GMT increase. This is true for all cascade sizes (see Fig. 5c–e). For temperatures above a 3 ^∘C GMT increase, cascades occur less frequently as most of the tipping elements already transgress their individual threshold before this temperature is reached. The most prevalent tipping cascades with two or three elements, as simulated in our network approach, consist of cascading transitions between the ice sheets and/or the AMOC, summing up to 80 % of all tipping cascades with two or three elements (Fig. 5f).

3.3 Different roles of tipping elements

For each of the four tipping elements, we systematically assess their role within the network model, generally distinguishing between initiators (triggering a cascade), followers (last element in a tipping chain) and mediators (elements in-between).

We find that in up to 65 % of all ensemble simulations, the Greenland Ice Sheet triggers tipping cascades. At the same time, it occurs as frequently in cascades as the other tipping elements (around 29 % of all cases, see Fig. 1). Thus, we call Greenland a dominant initiator of cascades. Following this argument for Greenland, the West Antarctic Ice Sheet is both an initiator and mediator of cascades, as it occurs often in cascades (31 %) and, likewise, often acts as the initiator (23 %). Although the frequency of occurrence in cascades is very similar for the AMOC as for the two large ice sheets, it is a dominant mediator of cascades as it does not initiate many cascades (13 %). Lastly, the Amazon rainforest is a pure follower in cascades because it is only influenced directly by the AMOC and cannot influence any other tipping element itself in our model due to the given interaction network structure (see Fig. 1). The reason why the ice sheets often act as initiators of tipping cascades in our model is likely because their critical threshold ranges tend to be lower than for the other tipping elements (see Fig. 4a). Many cascades are then passed on to other tipping elements, especially the AMOC. Thus, the role of the AMOC as the main mediator of cascades can be understood from a topological point of view as the AMOC is the most central network element with many connections to the other tipping elements. As such, the AMOC connects the two hemispheres and can be influenced by both the Greenland Ice Sheet and the (West) Antarctic Ice Sheet, as is also suggested by the literature (Wood et al., 2019; Ivanovic et al., 2018; Hu et al., 2013; Swingedouw et al., 2009; Rahmstorf et al., 2005).

3.4 Structural robustness and sensitivity analysis including ENSO

While many tipping elements (including the ice sheets, the AMOC and the Amazon rainforest) to a first approximation exhibit a transition between two or more alternative stable states, often described by the paradigmatic double-fold bifurcation (Scheffer et al., 2009; Lenton et al., 2008) as discussed above, tipping of the El Niño–Southern Oscillation (ENSO) could rather imply a transition from irregular oscillatory occurrences to a more permanent state of strong El Niño conditions (Dekker et al., 2018; Lenton et al., 2008; Kriegler et al., 2009). In coupled experiments for the AMOC and ENSO with conceptual models, it was found that a tipping AMOC can lead to a Hopf bifurcation in ENSO (Dekker et al., 2018; Timmermann et al., 2005). Overall, changes in the frequency of major El Niño events seem likely, also based on intermediate complexity and conceptual models (Dekker et al., 2018; Timmermann et al., 2005), but whether this poses the possibility of a permanent El Niño state remains debated. A more frequent occurrence of El Niño events could have strong impacts on global ecosystems up to a potential dieback of the Amazon rainforest (Duque-Villegas et al., 2019).

While some studies have emphasised the uncertainty about future ENSO changes (Kim et al., 2014; Collins et al., 2010), another study found that the frequency of El Niño events could increase twofold in climate change scenarios in simulations of the CMIP3 and CMIP5 climate model ensembles as well as in perturbed physics experiments (Cai et al., 2014). Also, some ENSO characteristics appear to respond robustly to global warming (Kim et al., 2014; Power et al., 2013; Santoso et al., 2013), such as an intensification of ENSO-driven drying in the western Pacific and rainfall increases in the central and eastern equatorial Pacific due to non-linear responses to surface warming (Power et al., 2013). Moreover, from an observational point of view, it was found that the global warming trend since the early 1990s has enhanced the Atlantic capacitor effect which might lead to more favourable conditions for major El Niño events on a biennial rhythm (Wang et al., 2017). Palaeoclimate evidence from the Pliocene (5.3–2.6 Myr before present) with atmospheric CO₂ levels comparable to today's conditions suggests that there may have been permanent El Niño conditions during that epoch (Fedorov et al., 2006; Ravelo et al., 2006; Wara et al., 2005). However, it must be noted that the Pliocene was different in terms of the continental configuration compared with today. Particularly, the Panama gateway was open for at least part of the Pliocene, resulting in tropical interactions between Atlantic and Pacific ocean waters (Haug and Tiedemann, 1998).

Given the particular uncertainties regarding ENSO compared with the other tipping elements considered in our analysis, we excluded it and its interactions with the other tipping elements from the main analysis above. However, we performed a comprehensive structural robustness and sensitivity analysis including ENSO as a tipping element (see also Figs. S3–S8): for this purpose, we choose to represent ENSO in the same way as the other tipping elements, although the use of Eq. (1) is not entirely appropriate for ENSO. Rather, the potential tipping behaviour could be conceptualised by a Hopf bifurcation (i.e. a transition from a limit cycle leading to oscillating behaviour to a stable fixed point attractor) instead of a fold bifurcation (Dekker et al., 2018; Timmermann et al., 2003; Zebiak and Cane, 1987).

A typical transition time of 300 years is chosen, the critical temperature threshold lies between 3.5 and 7.0 ^∘C above pre-industrial levels (Schellnhuber et al., 2016) and our analysis is based on simulations of 11 million ensemble members arising from the 27 different network combinations from the three unclear links AMOC → Amazon rainforest, West Antarctica → AMOC and Amazon rainforest → ENSO (see Fig. S3). The interactions including ENSO are described in detail in the Supplement (see Table S2 and description there).

Our robustness analysis reveals that the roles of the tipping elements remain qualitatively the same: the ice sheets remain strong initiators of tipping cascades (in 40 % of cases for the Greenland Ice Sheet, and 28 % of cases for the West Antarctic Ice Sheet). The AMOC mainly acts as a mediator and only initiates 5 % of all cascades (see Fig. S3). In this extended network of tipping elements, ENSO tends to take on an intermediate role. As it is strongly coupled to the Amazon rainforest, it initiates many cascades including the Amazon rainforest, especially at temperature levels above 3 ^∘C (see Fig. S4). However, ENSO also mediates tipping cascades from the AMOC to the West Antarctic Ice Sheet or the Amazon rainforest. Generally, we also find that the interactions destabilise the overall network of tipping elements apart from the Greenland Ice Sheet (Figs. S5, S6). The change in the critical temperature range for the Amazon rainforest is larger and is shifted more towards lower temperature levels due to the influence of ENSO. Overall, the model results remain robust, also with respect to the occurrence and size of tipping cascades (see Fig. S7), suggesting a certain degree of structural stability of our analysis.

4 Discussion and conclusions

It has been shown previously that the four integral components of the Earth's climate system mainly considered here are at risk of transgressing into undesirable states when critical thresholds are crossed (Schellnhuber et al., 2016; Lenton et al., 2008). Over the past decades, significant changes have been observed for the polar ice sheets as well as for the Atlantic Meridional Overturning Circulation (AMOC) and the Amazon rainforest (Lenton et al., 2019). Should these climate tipping elements eventually cross their respective critical temperature thresholds, this may affect the stability of the entire climate system (Steffen et al., 2018).

In this study, we show that this risk increases significantly when considering interactions between these climate tipping elements and that these interactions tend to have an overall destabilising effect. Altogether, with the exception of the Greenland Ice Sheet, interactions effectively push the critical threshold temperatures to lower warming levels, thereby reducing the overall stability of the climate system. The domino-like interactions also foster cascading, non-linear responses. Under these circumstances, our model indicates that cascades are predominantly initiated by the polar ice sheets and mediated by the AMOC. Therefore, our results also imply that the negative feedback loop connecting the Greenland Ice Sheet and the AMOC might not be able to stabilise the climate system as a whole, a possibility that was raised in earlier work using a Boolean modelling approach (Gaucherel and Moron, 2017).

While our conceptual model evidently does not represent the full complexity of the climate system and is not intended to simulate the multitude of biogeophysical processes or to make predictions of any kind, it allows us to systematically assess the qualitative role of known interactions of some of the most critical components of the climate system. The large-scale Monte Carlo approach further enables us to systematically take into account and propagate the substantial uncertainties associated with the interaction strengths, interaction directions and the individual temperature thresholds. This comprehensive assessment indicates structurally robust results that allow qualitative conclusions, despite all of these uncertainties.

In our Monte Carlo approach employed for propagating parameter uncertainties, we assume that all parameters including critical threshold temperatures and interaction link strengths are statistically independent. However, this is likely not the case in the climate system where, for example, interaction link strengths associated with the AMOC to Greenland and West Antarctica would be expected to be correlated. Further analyses would have to consider the effects of such interdependencies.

Overall, this work could form the basis of a more detailed investigation using more process-detailed Earth system models that can represent the full dynamics of each tipping element and their interactions. Major advances have been made in developing coupled Earth system models; however, computational constraints have so far prohibited a detailed interaction analysis as is presented in this work.

In the future, these more complex climate models might be driven with advanced ensemble methods for representing and propagating various types of uncertainties in climate change simulations (Daron and Stainforth, 2013; Stainforth et al., 2007), which would comprise a significant step forward in the current debate on non-linear interacting processes in the realm of Earth system resilience. Some examples of relevant processes that could be investigated with more complex models are the following: first, the changing precipitation patterns over Amazonia due to a tipped AMOC, i.e. whether rainfall patterns will increase or decrease and whether this would be sufficient to induce a tipping cascade in (parts of) the Amazon rainforest. This would shed more light on the AMOC–Amazon rainforest interaction pair. Second, the influence of the disintegration of the West Antarctic Ice Sheet on the AMOC could be further studied by introducing freshwater input into the Southern Ocean surrounding the West Antarctic Ice Sheet similar to the hosing experiments that have been performed for the Greenland Ice Sheet (Wood et al., 2019; Hawkins et al., 2011; Rahmstorf et al., 2005). Here, some studies suggest that freshwater input into the Southern Ocean at a modest rate would not impact the AMOC as much as freshwater input into the North Atlantic (Ivanovic et al., 2018; Hu et al., 2013; Swingedouw et al., 2009), while higher melt rates could have more severe impacts on the AMOC (Swingedouw et al., 2009). With carefully calibrated coupled ice–ocean models, including dynamic ice sheets (e.g. Kreuzer et al., 2020), ice–ocean tipping cascades could be studied in more detail.

Further, the timescales for potential tipping dynamics need to be more rigorously explored in contrast to the conceptual approach used here. It is important to note that the transition of one tipping element has a delayed effect on the other elements, especially in the case of the comparatively slowly evolving ice sheets. Their temperature threshold is lower than for the other tipping elements considered here, and their disintegration would unfold over the course of centuries up to millennia (Winkelmann et al., 2015; Robinson et al., 2012; Lenton et al., 2008). Therefore, meltwater influx into the ocean and changes in sea level would affect the state of other tipping elements only after a significant amount of time. Therefore, our analysis of emerging tipping cascades needs to be understood in terms of committed impacts over long timescales due to anthropogenic interference with the climate system mainly in the 20th and 21st centuries, rather than short-term projections.

Finally, it appears worthwhile to perform an updated expert elicitation along the lines of Kriegler et al. (2009), where additional interactions, tipping elements and a better understanding of the interaction strengths would help to narrow down the space of possible scenarios and uncertainties that have been investigated here.