The development of robust Early Warning Signals (EWSs) is necessary to quantify the risk of crossing tipping points in the present-day climate change. Classically, EWSs are statistical measures based on time series of climate state variables, without exploiting their spatial distribution. However, spatial information is crucial to identify the starting location of a transition process and can be directly inferred by satellite observations. By using complex networks constructed from several climate variables on the numerical grid of climate simulations, we seek for network properties that can serve as EWSs when approaching a state transition. We show that network indicators such as the normalized degree, the average length distance, and the betweenness centrality are capable of detecting tipping points at the global scale, as obtained by the MIT general circulation model in a coupled-aquaplanet configuration for CO concentration-driven simulations. The applicability of such indicators as EWSs is assessed and compared to traditional methods. We also analyze the ability of climate networks to identify nonlinear dynamical patterns.
Present-day anthropogenic CO forcing causes the Earth climate to vary in an extreme manner, which can lead to abrupt changes at a local or even global scale. Such abrupt changes in the climate are known as tipping processes. Detecting and anticipating them is crucial in our efforts to limit the consequences of climate change. Several techniques have been developed using time series analysis and have given rise to early warning signals (EWSs) of tipping processes. However, these methods do not consider the spatial evolution of climate change. To combine both temporal and spatial information, we consider the time evolution of climate networks. The idea is to grid the Earth’s surface with networks, thus capturing the spatial propagation of change. We test this method on a global-scale tipping using climate simulations. By including the information on the spatial distribution, these networks are able to serve as better EWS compared to classical methods and to reveal the regions where the tipping develops.
I. INTRODUCTION
Crossing tipping thresholds in the present-day climate change is raising increasing concern1,2 as the global mean temperature is getting closer to exceeding the Paris agreement range. Tipping elements have been identified with impacts at the regional and global scale.3–6 The interactions among these elements can likely induce tipping cascades with the risk of accelerating irreversible changes in the climate system.7–9 In this context, it becomes urgent to develop robust early warning signals (EWSs) for detecting the approach of tipping points.
In the past, a variety of methods has been explored to obtain EWS. The classical approach is based on the bifurcation theory and the phenomenon of “critical slowing down,”10,11 where changes in the variance, the lag-1 autocorrelation or high-order momenta of the system’s time series are used as EWS.12,13 However, several other techniques have been proposed in the literature such as estimating fluctuations decay,14 frequency change15 or flickering16 in noisy time series, and basin entropy17 near Hopf bifurcations. Moreover, physics-based indicators have also been tested for particular tipping elements, such as the freshwater transport for the tipping of the Atlantic Meridional Overturning Circulation (AMOC).18,19
In general, the above mentioned indicators provide information on the temporal development of a tipping process, without exploiting its spatial structure. However, critical slowing down arises when a system becomes increasingly sensitive to external perturbations and not necessarily near a catastrophic shift.20 This shows the need of including the information from the spatial processes during tipping events to have robust EWS and reduce the number of false positives.21 Spatial methods include the edge detectors in spatiotemporal data22 and indicators based on climate networks.23 These methods become particularly relevant in the current era of huge amounts of data coming from satellite observations with high spatiotemporal coverage.24 Network indicators have been successfully applied for studying climate subsystems, like the AMOC tipping25 and El Niño Southern Oscillation.26,27
Here, we investigate how network indicators behave as EWS during a tipping point at the planetary scale. Tipping between climatic attractors is simulated with the MIT general circulation model in a coupled-aquaplanet configuration28 by increasing or decreasing the atmospheric CO content at different rates. Using the Python package pyUnicorn,29 we construct networks based on linear and nonlinear spatial correlations of time series at each point of the numerical grid. We test several indicators that characterize the network topology on the local, mesoscopic, and global scales.30 We apply them to various networks, constructed from different climatic state variables, spanning the whole planet or only latitudinal strips with several time window sizes. The aim of this detailed analysis is to understand the ability of network indicators to predict the occurrence of global tipping points and to unveil nonlinear backbone structures in the climate system.
II. METHOD
In this section, we will briefly describe the main characteristics of the climate networks and the simulations data set we are going to use. To generate the climate networks and calculate the network indicators, we use the pyUnicorn package in python,29 as detailed in Appendix A.
A. Climate networks
A network is a set of points that interact with each other during a certain period of time. In this work, we will consider only undirected networks, i.e., without considering the direction of interrelationships.
The vertices (nodes) of a climate network are the spatial grid points of an underlying global climate data set, which can be derived from either observations or simulations. Specifically, we take as vertices the points on a longitude–latitude grid. Area-weighted indicators are then used to avoid the heterogeneous spatial distribution of vertices to bias the network topological proprieties29–31 (see Sec. II C and Appendix A).
In this work, we will consider two statistical ways of generating a complex network: the Pearson correlation coefficient (PCC) and the mutual information (MI).30 The reason for this choice is that PCC is commonly used in the literature but reveals only linear correlations between time series. In contrast, MI has the ability to reveal non-linear correlations, as shown in Donges et al.30
If two time series are heavily correlated linearly, the value will be large and will be close to . This is also true for some nonlinear correlations. However, cases where is large and small can unveil some non-linear features of the climate dynamics.30 Therefore, by comparing and , it is possible to identify a part of the contribution of non-linear dynamics in a complex network.
B. Networks indicators
All these indicators use different ranges of information. The degree centrality uses the local information from the vertices that are connected between each other. Whereas the local/global clustering coefficients provide information on a mesoscopic scale since they depend on the information from their connected neighbors. Finally, the average path length and the betweenness centrality use global topological information as they consider the shortest path among different connected pairs, which can span the whole network.30
C. Climate simulations
Using the MIT general circulation model (MITgcm),36–39 we consider the same coupled-aquaplanet configuration analyzed in Brunetti et al.,28 Ragon et al.,40 Brunetti and Ragon,41 i.e., a planet entirely covered by a 3000 m-deep ocean where full dynamical couplings between ocean and atmosphere are included. The CS32 cubed-sphere configuration (with grid points) was used, corresponding to a horizontal resolution of . For the networks, we interpolated this CS grid into a 2 longitude–latitude grid. Thus, the lattice is regular with nodes. In this configuration, five climatic attractors (which are stable over a multi-millennial time scale) have been found,28 denoted as snowball (S), waterbelt (WB), cold (C), warm (W), and hot (H) states as shown in the bifurcation diagram in Fig. 1. Here, we consider the atmospheric CO content as a forcing parameter. We have completed the bifurcation diagram in terms of the equilibrium surface air temperature (SAT), partially shown in the supplementary material of Brunetti and Ragon.41 This kind of diagrams allows one to know where the bifurcation-induced tipping points occur, the range of the stable branches and the regions of multistability, revealing the core dynamical structure of the climate system.
Bifurcation diagram for the coupled aquaplanet with the five steady states, the transient state and the two transitions C–H and C–S highlighted. Dashed lines correspond to qualitative sketches of unstable branches.
Bifurcation diagram for the coupled aquaplanet with the five steady states, the transient state and the two transitions C–H and C–S highlighted. Dashed lines correspond to qualitative sketches of unstable branches.
We investigate here the transitions at the global scale from one attractor to the other. We are going to focus on two types of transition: one from the cold state to the hot state (C–H) and one from the cold state to the snowball state (C–S), with intermediate transitions to a transient state around C (which is stable only over a timescale of hundred years, as verified in Brunetti and Ragon41) and to the waterbelt state, as shown in Fig. 1. The reverse transition from the hot state to the cold state is not possible in our configuration because of the hysteresis path between these steady states. Indeed, from the bifurcation diagram in Fig. 1, the climate system on the hot-state branch can only tip toward the snowball state by gradually reducing the forcing. Therefore, we choose to look at a transition starting from the same cold attractor as the upward transition, but going to colder states. With climate change we are expecting an upward transition (in terms of global temperature); nevertheless, it is important to verify which is the response of climate networks to both upward and downward transitions. Moreover, the presence of an intermediate transient state between the cold and waterbelt states enables us to estimate the network sensitivity.
Values of the start t1 and the end t2 of the tipping interval for different forcing rates in the C–H transition.
. | 0.1 ppm/yr . | 2 ppm/yr . | ||
---|---|---|---|---|
C–H . | . | . | . | . |
variable . | t1 (yr) . | t2 (yr) . | t1 (yr) . | t2 (yr) . |
SAT | 2400 | 3400 | 200 | 700 |
CC | 2300 | 3300 | 150 | 700 |
Humidity | 2400 | 3300 | 200 | 600 |
. | 0.1 ppm/yr . | 2 ppm/yr . | ||
---|---|---|---|---|
C–H . | . | . | . | . |
variable . | t1 (yr) . | t2 (yr) . | t1 (yr) . | t2 (yr) . |
SAT | 2400 | 3400 | 200 | 700 |
CC | 2300 | 3300 | 150 | 700 |
Humidity | 2400 | 3300 | 200 | 600 |
Values of the start t1 and the end t2 of the tipping interval for different forcing rates in the C–S transition.
. | 0.1 ppm/yr . | 0.4 ppm/yr . | ||
---|---|---|---|---|
C–S . | . | . | . | . |
variable . | t1 (yr) . | t2 (yr) . | t1 (yr) . | t2 (yr) . |
SAT | 300 | 1750 | 150 | NA |
CC | 300 | 1750 | 150 | NA |
Humidity | 300 | 1650 | 150 | NA |
. | 0.1 ppm/yr . | 0.4 ppm/yr . | ||
---|---|---|---|---|
C–S . | . | . | . | . |
variable . | t1 (yr) . | t2 (yr) . | t1 (yr) . | t2 (yr) . |
SAT | 300 | 1750 | 150 | NA |
CC | 300 | 1750 | 150 | NA |
Humidity | 300 | 1650 | 150 | NA |
Evolution of the edge density against the threshold value using (a) Pearson’s correlation coefficient (PCC) and (b) mutual information (MI) on the attractor (on the cold state at yr during the C–H transition shown in Fig. 3) and during tipping (at yr).
Evolution of the edge density against the threshold value using (a) Pearson’s correlation coefficient (PCC) and (b) mutual information (MI) on the attractor (on the cold state at yr during the C–H transition shown in Fig. 3) and during tipping (at yr).
Tipping is induced by varying the atmospheric CO content. We have performed several simulations using different forcing rates, ranging from 0.1 (which is sufficiently low to guarantee that the system relaxes toward the attractor at each step) to 2 ppm/yr (which is the present-day rate).
We consider three state variables with yearly frequency that provide different information on the climate dynamics. SAT is treated in detail in the main part of this article, while specific humidity (SH) and cloud cover (CC) can be found in Appendixes A– C. SAT and SH are prognostic variables in MITgcm, directly linked to the energy and water cycles, respectively, whereas the CC fraction is a diagnostic variable obtained from temperature, pressure, and specific humidity.42
D. Determination of the tipping time interval
To fully describe the transition, we define three main intervals for the temporal evolution of a climatic state variable: on the initial attractor, during the tipping process and on the final attractor. They are marked by the start and the end of tipping and can be precisely defined from the simulations by the following two criteria:41 (i) the ocean surface radiation imbalance is lower than 0.5 W m in the absolute value; (ii) under forcing conditions, the system does not shift monotonically in time toward other values. When one of these criteria is not satisfied, the system enters in the tipping (or transition) interval. The precise definition of the tipping interval is important regarding our analysis on detecting tipping points and early warning signals’ robustness. The above analysis is repeated for each state variable in each simulation, and the values of and are reported in Table I for the C–H transition and in Table II for the C–S transition.
Contrary to C–H, C–S shows the presence of an intermediate transient attractor, which has been analyzed in detail in Brunetti and Ragon.41 In addition, the ocean surface imbalance displays large fluctuations until around 1000 yr. This strong difference in the evolution profile compared to the C–H transition allows us to test the robustness of the networks and their capability of detecting fluctuations in the tipping process.
Temporal evolution of the mean surface air temperature during the transition from the cold to the hot state (C–H, left) and the cold to the snowball state (C–S, right), passing through the cold to transient (C–T), transient to waterbelt (T–WB), and waterbelt to snowball (WB–S) tipping intervals, when the atmospheric CO content increases or decreases at a rate of 0.1 ppm/yr. For selected time slices on the attractors and during tipping, the map of the normalized degree centrality for the PCC network is shown in the bottom panels.
Temporal evolution of the mean surface air temperature during the transition from the cold to the hot state (C–H, left) and the cold to the snowball state (C–S, right), passing through the cold to transient (C–T), transient to waterbelt (T–WB), and waterbelt to snowball (WB–S) tipping intervals, when the atmospheric CO content increases or decreases at a rate of 0.1 ppm/yr. For selected time slices on the attractors and during tipping, the map of the normalized degree centrality for the PCC network is shown in the bottom panels.
Temporal evolution of four PCC network indicators constructed with three different window sizes (WS) for the C–H transition.
Temporal evolution of four PCC network indicators constructed with three different window sizes (WS) for the C–H transition.
E. Threshold value for edge definition
The choice of the threshold above which pairs of vertices are considered to be connected is a critical step in the network construction. For PCC, we tested threshold values from 0.1 to 0.98 with increments of 0.02. For MI, we tested values from 0.005 to 1.5 with increments of 0.005. In both cases, we performed the test once while the system is on the attractor and once during the tipping process.
As shown in Fig. 2(a), during tipping, there is a saturating effect for PCC if we choose a value lower than 0.4 (red thick curve). Therefore, if we take a threshold below this value, the evolution of the edge density of the PCC is saturated, and the network dynamics cannot be properly studied during tipping. Therefore, the threshold in PCC needs to be carefully chosen. To ensure that our analysis is not affected by the saturating effect while still considering as many edges as possible, we use a threshold value of 0.5 for PCC. In contrast, there is no such saturating effect in MI, as shown in Fig. 2(b).
If we want to compare the two methods of network construction, it is more relevant to use the same edge density and not the same thresholds as they are not commensurate.26 The edges that are not present in the PCC adjacency matrix but are found in the MI matrix contribute to the non-linear dynamics of the complex network.30
We have also examined what happens when the edge density is fixed instead of the threshold. In this case, even though the overall behavior during the tipping process is similar to what is observed with a fixed threshold, it is difficult to track the evolution of all the spatial regions because only the strongest and most dominant edges exist and, therefore, other weaker correlations are neglected.
III. RESULTS
As detailed in the Method section, we consider how climate networks behave during an upward transition going from the cold to the hot state (C–H) and a downward transition going from the cold state to the snowball (C–S). We focus on networks generated using PCC, except in Subsection III E where networks generated with MI are introduced to search for their potential non-linear behavior.
A. General behavior
The evolution of SAT during the C–H transition is shown in the left panel of Fig. 3. The initial state is on the cold attractor at an atmospheric CO content of 320 ppm. By increasing it by 0.1 ppm per year, the system evolves along the stable branch of the cold attractor until the tipping starts. At this point, SAT increases much faster, finally arriving at the hot attractor. During this evolution, the climate network evolves accordingly. As an example, the bottom panels of Fig. 3 show the behavior of the normalized degree centrality in a PCC network. We can see that during the tipping process, the number of spatial correlations increases,21 as observed in phase transitions of thermodynamic physical systems.23 This is predicted by the slaving principle,43 where at a bifurcation all the dynamical modes relax and adapt to the slowly evolving one, thus giving rise to spatial correlations that extend to the whole system.
Before and after the tipping process, the network describes the internal variability of the attractor dynamics. In the cold state, the main correlations occur in the polar regions, while in the hot state, the patterns are more uniform and also involve the mid-latitude regions, in agreement with the principal component analysis reported in Brunetti et al.28
In the following, we will characterize in detail the ability of climate networks on the basis on both their spatial and temporal information to describe this type of global-scale tipping.
B. Sensitivity tests
As a preliminary step, we need to set some free parameters and evaluate their effects on the network dynamics. In particular, we consider the impact of changing the window size, as well as the discretization and intensity of the forcing rate. Unless otherwise specified, the tests are made using the C–H transition with a forcing rate of 0.1 ppm/yr and SAT as a climate state variable.
1. Time window size
Testing this free parameter is crucial for the reliability of the network outcomes analysis. Varying the window size affects the amount of information from the time series and, therefore, the significance of the network indicators.44 The task is to find a compromise between maximizing the amount of information at each time step on one side, and a high temporal resolution on the other side, with a sufficient number of points in the temporal evolution of the network in order to provide a proper temporal characterization of the tipping process.
The window size needs to be larger than the forcing frequency and smaller than the duration of the tipping process. In the C–H transition, the tipping process lasts nearly 1000 yr, therefore, if we want to track the network behavior during the transition, we need to generate a new network at least ten times. Thus, we test the following window sizes: 20, 50, and 100 yr. In our case, we are dealing with temporal series with yearly frequency. An analogous argument holds for time series with monthly, daily, or hourly frequency, although in this case, the daily and seasonal periodicity have to be properly taken into account when calculating the temporal correlations.
In Fig. 4, we see that the network indicators, averaged over the spatial domain, show a similar behavior for the three investigated window sizes. For all the indicators, the window size of 100 yr gives a stronger signal than the others. The window size of 20 yr gives a noisy evolution of the network behavior, sometimes without a clear signal when crossing the tipping point. This behavior can be explained by information theory, since short time series provide less information than longer ones and have higher uncertainty.44 Therefore, it is better to work with longer time series. However, in the search for early warning signals, we need to have enough points to detect a change of network behavior. The compromise is to choose a window size of 50 yr which is sufficient to generate relevant statistics and is 5% of the tipping transition. This is the value that we will use from now on.
2. Discretization of the forcing rate
There are two types of network stress: internal stress and environmental stress.45 The internal stress, which depends on the change of the number of vertices, does not occur in our climate networks by construction. On the other hand, the environmental stress acts on the spatial patterns of the network. As said before, the transition from one state to another is induced by an increase or decrease of the external forcing, in our case, the atmospheric CO content. Thus, we use this forcing as environmental stress, by considering the same average rate of 0.1 ppm/yr but with different temporal discretizations, 0.1 ppm every year, 2 ppm every 20 yr, and 4 ppm every 40 yr, starting and ending at the same CO concentration values. We focus the analysis around tipping.
To understand how the environmental stress acts on the network, we first need to verify if it affects the state variable that is extracted from the climate model. In Fig. 5, we compare the temporal evolution of SAT zonal averages for coarse and fine discretizations of the forcing. There are indeed localized significant differences between the cases, showing that the climate model is sensitive to the change of the forcing step. Now, we need to understand how this sensitivity affects spatial correlations. Figure 6 shows the evolution of the four considered network indicators, defined in Sec. II B. It turns out that the overall behavior of these indicators is similar and quite robust with respect to the change of the forcing step. The small discrepancies in the indicators are just unveiling the network sensitivity to the fact that in the case of coarse discretization the forcing is stronger, so that the system is farther away from the equilibrium than with fine discretization.
Temporal evolution of the zonal average of SAT for an increase of CO concentration of (a) 2 ppm each 20 yr; (b) 0.1 ppm per year. (c) Difference between (a) and (b).
Temporal evolution of the zonal average of SAT for an increase of CO concentration of (a) 2 ppm each 20 yr; (b) 0.1 ppm per year. (c) Difference between (a) and (b).
Temporal evolution of four network indicators in simulations with different discretization of the forcing rate: 0.1 ppm per year, 2 ppm each 20 yr, and 4 ppm each 40 yr.
Temporal evolution of four network indicators in simulations with different discretization of the forcing rate: 0.1 ppm per year, 2 ppm each 20 yr, and 4 ppm each 40 yr.
From now on, we will use yearly forcing steps.
3. Forcing rates
For the C–H transition, we compare the two following forcing rates: 0.1 and 2 ppm/yr. The first is sufficiently low to guarantee that the system follows the stable branch of the cold state until a bifurcation-induced tipping occurs, while the second corresponds to the present-day forcing rate and can give rise to rate-induced tipping46–49 in the coupled-aquaplanet configuration.41 Indeed, we verified that on the cold attractor (first 200 yr) the mean global SAT is C without forcing, C for the low forcing rate and C for the high forcing rate. The first two values are the same within the confidence interval, showing that the system has time to relax back to the attractor under the low forcing rate.50
As shown in Figs. 7(a) and 7(b), the SAT evolution is radically different. Starting from the same initial condition, the attractor lasts 200 yr for the 2 ppm/yr forcing rate, whereas it lasts 2400 yr for the 0.1 ppm/yr case. This has of course a direct consequence on the network indicators, which are not constant in the cold-attractor region, as shown in panels (e) and (f). In a general manner, we see that the behavior of the networks is similar and the indicators increase/decrease during the tipping process. While the latter starts for a CO content of 560 ppm when the rate is 0.1 ppm/yr, it is delayed to around 720 ppm for the 2 ppm/yr rate, thus overshooting the bifurcation-induced tipping that can be seen in Fig. 1. Hence, even if the tipping mechanism can be different in the two cases (bifurcation-induced vs rate-induced tipping for low or high forcing rates, respectively), the forcing rate does not have a major impact on the network, which in both cases reacts to the tipping process (see bottom row of Fig. 7).
Comparison of the network behavior for two forcing rates: 0.1 ppm/yr [panels (a), (c), and (d)] and of 2 ppm/yr [panels (b), (e), and (f)]. The shaded region corresponds to the tipping interval from the cold state to the hot state.
Comparison of the network behavior for two forcing rates: 0.1 ppm/yr [panels (a), (c), and (d)] and of 2 ppm/yr [panels (b), (e), and (f)]. The shaded region corresponds to the tipping interval from the cold state to the hot state.
C. Global and zonal networks
After the sensitivity tests, we describe now the behavior of the networks constructed over the entire numerical grid and over limited zonal areas during the C–H and C–S transitions.
The aquaplanet is divided into five latitudinal strips (covering the entire longitudinal range), namely, the north and south polar regions [with latitude ranging in and , respectively], the mid-latitude regions [ and ( ], and the equatorial region [ ], so that the zonal behavior of the networks can be investigated. The aim is to investigate if networks can identify dominant dynamical regions during the stable phase on an attractor and during the tipping process. The networks are generated independently for the full aquaplanet and for each zonal region. The network generation is thus performed by giving as input variables to pyUnicorn their values in a specific region (note that there is no significant difference if we generate the network globally and then divide it in regions).
As the aquaplanet is symmetric between the Northern and the Southern Hemispheres, we expect the network to behave in a symmetrical way with respect to the Equator. This can be verified by comparing the two mid-latitudinal and the two polar regions. There is indeed a symmetry between the indicators for the two hemispheres, as shown in Appendix B. In addition to unveiling the robustness of the networks, comparing the evolution of network indicators for the corresponding zones of both hemispheres bears information on how the internal variability of the simulations translates onto the networks.
1. Cold to hot (C–H) transition
As shown in Fig. 8 (main panels), the evolution of the network indicators of the C–H transition displays a strong and unique signal during tipping. For the normalized degree centrality, the global clustering coefficient and the average length distance, a strong increase close to and during the crossing of the tipping process is observed. This is consistent with the previously mentioned slaving principle near a bifurcation point, with spatial correlations occurring over the whole system.21,43 The behavior of the average length distance is expected since, as the system transitions, the vertices tend to correlate over longer distances. Conversely, the betweenness centrality is strongly decreasing and is approaching zero when the system tips. This is also expected, as the vertices correlate over longer distances so that the normalized degree centrality approaches 1. As all pairs of points tend to be connected, the possibility of alternative paths vanish.
Temporal evolution of the four network indicators for SAT in the C–H transition at a rate of 0.1 ppm/yr, at global (main plot, black curve) and zonal scales (blue, green and red curves correspond, respectively, to the polar, mid-latitude and equatorial regions, as indicated on the aquaplanet). The time range shaded in pink corresponds to the transition.
Temporal evolution of the four network indicators for SAT in the C–H transition at a rate of 0.1 ppm/yr, at global (main plot, black curve) and zonal scales (blue, green and red curves correspond, respectively, to the polar, mid-latitude and equatorial regions, as indicated on the aquaplanet). The time range shaded in pink corresponds to the transition.
The normalized degree centrality and the global clustering coefficient show higher values in the polar regions on the cold attractor, and they increase in all the zonal regions during the tipping process. The equatorial region shows an earlier signal of transition. This is likely due to the main dynamical feature in hot-climate attractors,28 where high surface temperatures and the consequent large water vapor-holding capacity of air can initiate deep convection cells, with latent heat released by the condensation of large amounts of moist, ascending air.51 The earlier signal of transition in the equatorial region is indeed also observed when specific humidity is used as state variable to construct the networks, as shown in Appendix C. The same strong signature in the equatorial region is observed in the average length distance. Finally, the intensity of the internal variability on the cold attractor of the betweenness centrality is smaller at mid-latitudes.
The evolution of the network indicators in the zonal regions unveils the importance of the local dynamics in climate. Some regions respond earlier to tipping than others and show different degrees of internal variability, depending on the attractor.
2. Cold to snowball (C–S) transition
This transition is shown in Fig. 3, right panels. Starting from the cold attractor and under a reduction of the CO content at a rate of 0.1 ppm/yr, the sea ice undergoes two main episodes of localized disruption and reconstruction at the ice edge at mid-latitudes, as can be seen from the simulation results (not shown). The mean SAT responds to these episodes through the ice-albedo feedback on a timescale of some decades.28 Afterwards, the climate system initiates the longest tipping from the cold to a transient state (denoted as C–T) at around C. After a shift to the waterbelt (T–WB), the climate system finally stabilizes to the snowball state (WB–S), where the ocean surface is completely covered by ice.
As shown in Fig. 9, the networks indicators react to the fluctuations before the C–T tipping, showing their high sensitivity to local spatial features. In a stronger manner, they react during C–T, T–WB, and WB–S tipping. The fact that the networks are sensitive to a transition that is not as clear as C–H unveils their reliability. This is especially relevant in the case of the Earth climate, where we do not necessarily expect direct transitions as C–H because of the presence of several interacting tipping elements. Moreover, our results show that networks can characterize both upward and downward transitions.
At the beginning, the betweenness centrality shows lower values in the polar region with respect to the other regions, while the other three network indicators show larger values, showing that active dynamical features are dominant at high latitudes where the ice area is increasing. Also the equatorial region is showing a strong signal during the main transition from the cold to the transient state (C–T) as it was the case for the C–H simulation. However, contrary to the C–H transition, the behavior of the equatorial region in the C–T transition is driven by the formation of the ice while it extends toward low latitudes. Therefore, we will consider the ice extent as the area where EWS are estimated in Sec. III D 2.
D. Early warning signals (EWSs)
We seek for statistical network indicators that are able to anticipate the transition by giving a clear and non-ambiguous signal. For this, we are going to use our regional study as it has unveiled that some regions are dynamically active and dominant during the transitions, thus more sensitive than the others. For SAT, the dominant region during the transition appears to be the equatorial region and the polar region for the cloud cover fraction, as shown in Appendix C. Therefore, the idea is to use the temporal evolution of the network indicators in the dominant region for each of the climate state variables, and compare it to traditional EWS methods10,11,52,53 such as the variance, the autocorrelation at lag-1 and the kurtosis.
These three classical indicators are estimated in the same region as the network ones. However, one has to keep in mind that the information of the regional dominance was extracted from the behavior of the complex network and not from a classical time series analysis. Here, we perform the EWS analysis for the two transitions, C–H and C–S. We are going to detail the results only for SAT, while the analysis made for the cloud cover fraction and the relative humidity provide similar results (not shown).
1. Cold to hot (C–H) transition
Figure 10 shows differences in the behavior of the candidate EWS indicators. On the plots, the 1 , 2 , and 3 intervals are labeled as multiples of the standard deviation of the distribution of each indicator on the cold attractor. A point is identified as an EWS when it overcomes the 3 interval. This strict limitation gives a strong confidence in the signal and overcomes the small fluctuating behavior of the networks due to the internal variability of the climate simulations. However, the choice of the 3 does not critically influences the number of EWS in our case.
Comparison of six EWS indicators in the equatorial region, during the C–H transition at a rate of 0.1 ppm/yr. (a)–(c). are network based and (d)–(f). are time series based.
Comparison of six EWS indicators in the equatorial region, during the C–H transition at a rate of 0.1 ppm/yr. (a)–(c). are network based and (d)–(f). are time series based.
For the three network-based methods [panels (a)–(c)], we see that the EWS occurs at 1850 yr for the normalized degree centrality, the average length distance and the betweenness centrality, whereas the kurtosis gives a signal of EWS at 2550 years [panel (d)]. The autocorrelation at lag-1 [panel (e)], which should increase near a tipping point,52 gives a signal at 2650 yr, with three contradictory low-level signals labeled by the red arrows. Also, we notice that the 3 interval takes 75% of the panel (e), due to the high internal variability on the cold attractor, which produces negative spikes exceeding in the autocorrelation at lag-1, challenging the normality assumption in that case. Finally, the variance decreases on the attractor and becomes an EWS at 2700 yr.
The three network indicators give comparable results, even if the ones using global topological information, such as the average length distance and the betweenness centrality, turn out to be more noisy than the normalized degree centrality, which is a local indicator. In summary, networks-based indicators, which take into account both spatial and temporal features of the dynamical process, perform better as EWS than the traditional ones. We also note that kurtosis and variance drift substantially while the system is still well on the attractor, which artificially widens their standard deviations. Conversely, the network indicators are remarkably stable over the attractor, which facilitates the definition of their typical behavior and eases the detection of the sudden deviations at the approach of tipping.
2. Cold to snowball (C–S) transition
We perform a similar analysis for the C–S transition. In this case, we focus on the ice formation, since we saw in Sec. III 1 2 that this is the main dynamical process that influences the overall dynamics. Therefore, the region where we construct the networks is the ice extent area, which evolves in time. The traditional indicators are also estimated in this region.
As shown in Fig. 11, the network-based indicators and the kurtosis to a certain extent detect at 300 yr the ice-disruption episodes before the main tipping (C–T). As already remarked in the case of the C–H transition (Sec. III D 1), the strong fluctuations in the autocorrelation at lag-1 do not enable us to see the EWS. These fluctuations are an artifact due to an extreme temporal stability of the zonally averaged SAT, within C. Therefore, the autocorrelation at lag-1 is mostly measuring random noise and is irrelevant. In contrast, during tipping, the fluctuations of the SAT increase to 0.1 C and become physically meaningful, so that the autocorrelation becomes relevant and its value stabilizes tending to 1.
Comparison of six EWS in the ice-extent region during the C–S transition at a rate of 0.1 ppm/yr, where (a) and (b) are network based and (d)–(f) are time series based.
Comparison of six EWS in the ice-extent region during the C–S transition at a rate of 0.1 ppm/yr, where (a) and (b) are network based and (d)–(f) are time series based.
Finally, the standard deviation shows a signal at 750 yr during C–T tipping, with strong oscillations within the 3 interval during the abrupt episodes of ice disruption. However, the fast succession of stable and tipping intervals in the C–S transition prevents an EWS to display clearly for each tipping process.
E. Non-linear climate dynamics
An additional advantage of climate networks is that they allow to find where some nonlinear dynamical patterns develop in the climate system. This can be understood by identifying the non-linear edges in the climate network. To do so, we generate a network with PCC to obtain the linear edges, and with MI to obtain the non-linear edges. Then, to perform the difference between these two networks in order to keep the non-linear part only, we use the method described in Donges et al.30 We compare the outcomes at two different time slices on the cold attractor and during the tipping process of the C–H transition.
A necessary condition that arises to do this comparison is to have the same edge density when generating the networks with PCC and MI.30 Thus, as in all the above analysis, we set the PCC threshold to during tipping. At 2650 yr [see the evolution plot in Fig. 3(a)], this corresponds to an edge density of 0.05. The same edge density is obtained in the MI network with a threshold of . Since we showed in Sec. II E that the PCC threshold during tipping is not saturated at 0.5, we keep an edge density of 0.05 for all the generated networks at both time slices, in the cold attractor (at yr) and during tipping (at yr), as summarized in Table III.
Threshold and edge density values for the cold attractor (50 yr) and tipping (2650 yr).
. | Cold attractor . | Tipping . | ||
---|---|---|---|---|
. | PCC . | MI . | PCC . | MI . |
Threshold | 0.10 | 0.03 | 0.5 | 0.29 |
Edge density | 0.05 | 0.05 | 0.05 | 0.05 |
. | Cold attractor . | Tipping . | ||
---|---|---|---|---|
. | PCC . | MI . | PCC . | MI . |
Threshold | 0.10 | 0.03 | 0.5 | 0.29 |
Edge density | 0.05 | 0.05 | 0.05 | 0.05 |
In Fig. 12, we show the spatial distribution of the proportion of edges in each vertex for the PCC and MI networks, respectively, as well as their difference, at the two time slices. We see that on the cold attractor (panels in the left column), the strongest linear correlations occur in the polar region, while non-linear ones are quite homogeneously occurring over all the planet. In contrast, during the tipping process (panels in the right column), the strongest linear correlations are due to the propagation of Kelvin waves in the equatorial region, while the patterns in the non-linear correlations correspond to the propagation of atmospheric Rossby waves, which transmit circulation anomalies in the surface temperature from the tropical region into the midlatitudes in each hemisphere.51,54,55 These linear (in blue) and non-linear (in red) features are clearly evident in the difference map [panel (f)]. The behavior of climate networks helps to evidence physical teleconnections at large scales, which become dominant during the tipping process.
Spatial distribution of the proportion of edges in each vertex. The first column corresponds to the cold attractor region ( yr) and the second column to tipping ( yr), with (a) and (b) PCC, (c) and (d) MI and (e) and (f) the difference between MI and PCC. The color bar indicates the proportion of edges for each vertex, with blue (negative) for a PCC dominance and red (positive) for a MI dominance.
Spatial distribution of the proportion of edges in each vertex. The first column corresponds to the cold attractor region ( yr) and the second column to tipping ( yr), with (a) and (b) PCC, (c) and (d) MI and (e) and (f) the difference between MI and PCC. The color bar indicates the proportion of edges for each vertex, with blue (negative) for a PCC dominance and red (positive) for a MI dominance.
IV. DISCUSSION AND CONCLUSIONS
In summary, using coupled-aquaplanet simulations showing a global-scale tipping from a climatic attractor to another, we have verified that network indicators have much better properties as early warning signals (EWSs) of tipping processes than traditional ones. This is due to the fact that they encode both spatial and temporal information.
Moreover, they show high sensitivity to local dynamical features, like sea ice formation or disruption, which enables one to use network indicators also for tipping occurring at the regional level. This is relevant, since some tipping elements in the climate system,5 like the Amazon rainforest dieback or the collapse of Greenland and West Antarctic ice sheets, can first have an impact on the subcontinental scale before affecting the functioning of the whole climate at the global scale through tipping cascade .9
We have tested that the ability of network indicators as EWS persists in systems forced at high rates with an overshoot of the bifurcation-induced tipping point. Beside rate-induced tipping, other tipping mechanisms can develop in the climate system, such as noise- or shock-induced tipping.41,56 These occur when the internal variability of the attractor dynamics reaches a critical threshold close to the basin boundary between alternative attractors, and can be triggered, for example, by perturbations in the volcanic activity57 or in the biological pump.58,59 Finding characteristic signatures in the behavior of network indicators and topology during tipping processes of different nature requires further investigation.
Finally, an additional advantage of complex networks is that alternative ways of network construction permit to isolate a part of the non-linear dynamical components. In the case of coupled-aquaplanets, we have seen that the main nonlinear feature on the interannual time scale is the propagation of Rossby waves from the Equator to mid-latitudes.55 In general, this type of analysis can be used to characterize the backbone structure of each attractor and the main teleconnections taking place on each attractor and during the tipping process. Using time series with daily or monthly frequencies, climatic phenomena with different characteristic timescales can be investigated using, for example, the so-called ordinal analysis.60,61 Moreover, the presence of non-linearity in the time series can be tested using surrogate datasets.62
Considering the huge current production of satellite data, methods that exploit both their spatial and temporal coverages, like those for network constructions and indicators, become of particular interest. We plan to apply such methods to present-day climate simulations participating in the Tipping Point Modelling Intercomparison Project (TIPMIP), which is presently defining its protocols and analysis procedures.
ACKNOWLEDGMENTS
We thank Cristina Masoller for useful discussions. We acknowledge the financial support from the Swiss National Science Foundation (Sinergia Project No. CRSII5_213539).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Laure Moinat: Conceptualization (supporting); Formal analysis (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Jérôme Kasparian: Methodology (equal); Supervision (equal); Writing – review & editing (equal). Maura Brunetti: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study and python scripts are available from the corresponding author upon reasonable request and is openly available on a Yareta Repository at the University of Geneva.63 Data were generated by the MIT general circulation model (version c67f) and by the python package pyUnicorn, which are openly available at https://github.com/MITgcm and https://www.pik-potsdam.de/∼donges/pyunicorn, Refs. 29 and 36, respectively.
APPENDIX A: THE pyUnicorn PACKAGE
The pyUnicorn package in python29 has the ability to combine the non-linear time series analysis and the complex network methods applied to climate. The package is open source. It can be applied to other fields, and therefore, we are not going to detail all the classes but only the ones that are relevant regarding our analysis. Its main functionality is to import several time series, to generate the complex network with the possibility to use different statistical methods for the threshold, and to give in output several files with the different analyses that have been performed.
Its inner structure is mainly dependent on the Numpy, Scipy, and igraph python libraries. The basic network generation is performed using the igraph library. More costly computation is done using C, C++, and FORTRAN compilers in Cython.64
This package is divided into five big classes: core, functnet, climate, timeseries, and utils. We are going to focus on the core and climate classes as we are only using these ones. The core class is the basic building of the network as it allows one to import data in the NetCDF format using the Data class. The climate class is built on GeoNetwork and Data, allowing one to build networks using different statistical constraints, such as climate.TsonisClimateNetwork and climate.MutualInfoClimateNetwork. The statistical constraints can be either fixed by climate.setthreshold or climate.setlinkdensity. All the functions that use the climate class take into consideration the heterogeneous spatial distribution of the vertices when the option node_weight_type = ‘surface’ is used, as in our case.
APPENDIX B: HEMISPHERICAL SYMMETRY
The symmetric behavior of the network indicators in the Northern and Southern Hemispheres of the aquaplanet can be observed by comparing the two mid-latitudinal and the two polar regions, as shown in Fig. 13. The internal variability on the attractor is translated to the small fluctuations of the network indicators, unveiling their robustness.
Comparison of the evolution of the network indicators for SAT in the two hemispheres in the case of the C–H transition at a rate of 0.1 ppm/yr: (a) mid-latitudinal regions and (b) polar regions.
Comparison of the evolution of the network indicators for SAT in the two hemispheres in the case of the C–H transition at a rate of 0.1 ppm/yr: (a) mid-latitudinal regions and (b) polar regions.
APPENDIX C: GLOBAL AND ZONAL NETWORKS FOR CLOUD COVER AND SPECIFIC HUMIDITY
We report here the behavior of networks constructed from two state variables complementary to the SAT detailed in the main text: cloud cover and specific humidity. The behavior of the same four network indicators as in the main text during the C–H and C–S transitions is shown in Figs. 14–17, which constitute direct counterparts of Figs. 8 and 9. The indicators constructed using cloud cover and specific humidity show a similar behavior to those based on SAT that we have discussed in the main text, in both the global and zonal regions. Thus, all these variables, which are accessible from satellite data, can be equally used to construct climate networks.
Temporal evolution of the four network indicators for specific humidity in the C–H transition at a rate of 0.1 ppm/1yr, at global and zonal scales.
Temporal evolution of the four network indicators for specific humidity in the C–H transition at a rate of 0.1 ppm/1yr, at global and zonal scales.
Temporal evolution of the four network indicators for cloud cover in the C–H transition at a rate of 0.1 ppm/1yr, at global and zonal scales.
Temporal evolution of the four network indicators for cloud cover in the C–H transition at a rate of 0.1 ppm/1yr, at global and zonal scales.
Temporal evolution of the four network indicators for specific humidity in the C–S transition at a rate of 0.1 ppm/yr, at global and zonal scales.
Temporal evolution of the four network indicators for specific humidity in the C–S transition at a rate of 0.1 ppm/yr, at global and zonal scales.
Temporal evolution of the four network indicators for cloud cover in the C–S transition at a rate of 0.1 ppm/yr, at global and zonal scales.
Temporal evolution of the four network indicators for cloud cover in the C–S transition at a rate of 0.1 ppm/yr, at global and zonal scales.