<figure><img src="image_1.png" /></figure>

HAL
open science

# A punctuated equilibrium analysis of the climate evolution of cenozoic exhibits a hierarchy of abrupt transitions

Denis-Didier Rousseau, Witold Bagniewski, Valerio Lucarini

► To cite this version:

Denis-Didier Rousseau, Witold Bagniewski, Valerio Lucarini. A punctuated equilibrium analysis of the climate evolution of cenozoic exhibits a hierarchy of abrupt transitions. Scientific Reports, 2023, 13, 10.1038/s41598-023-38454-6. insu-04195490

HAL Id: insu-04195490

https://insu.hal.science/insu-04195490

Submitted on 4 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# A punctuated equilibrium analysis of the climate evolution of cenozoic exhibits a hierarchy of abrupt transitions

Denis-Didier Rousseau<fnref n="1" />,²,³✉, Witold Bagniewski⁴ & Valerio Lucarini⁵,⁶,⁷

The Earth's climate has experienced numerous critical transitions during its history, which have often been accompanied by massive and rapid changes in the biosphere. Such transitions are evidenced in various proxy records covering different timescales. The goal is then to identify, date, characterize, and rank past critical transitions in terms of importance, thus possibly yielding a more thorough perspective on climatic history. To illustrate such an approach, which is inspired by the punctuated equilibrium perspective on the theory of evolution, we have analyzed 2 key high-resolution datasets: the CENOGRID marine compilation (past 66 Myr), and North Atlantic U1308 record (past 3.3 Myr). By combining recurrence analysis of the individual time series with a multivariate representation of the system based on the theory of the quasi-potential, we identify the key abrupt transitions associated with major regime changes that separate various clusters of climate variability. This allows interpreting the time-evolution of the system as a trajectory taking place in a dynamical landscape, whose multiscale features describe a hierarchy of metastable states and associated tipping points.

Early evidence of abrupt transitions in Camp Century and Dye 3 Greenland ice cores<fnref n="1" />,² attracted a lot of attention from the paleoclimatic community before being well acknowledged and understood. These findings gave evidence of a sequence of previously unknown abrupt climatic variations. However, such transitions did not seem to agree with other marine and terrestrial records, which led to considerable debate in the field³⁻⁵. Over the course of several decades spent retrieving and studying more detailed paleorecords, the existence of these rapid climatic variations, known as Dansgaard-Oeschger events (DO), has become well accepted. Further support for these findings has recently been provided by the identification of additional abrupt transitions from the NGRIP ice core, which has been made possible by the increased temporal resolution of the record⁶. These additional events correspond to changes of either short duration or amplitude in the stable ¹⁸O over ¹⁶O isotopes ratio δ¹⁸O that visual or standard statistical inspections do not necessarily flag. We remind that the fractionation of oxygen isotopes can be used to reconstruct local temperature conditions.

The Earth climate has experienced numerous abrupt and critical transitions during its long history, well beyond the specific examples above⁷,⁸. As discussed in Rothman⁹ following Newell¹⁰, rapid variations change are more likely to lead to catastrophic consequences in the biosphere—as in the extreme case of mass extinction events—because it is hard for the evolutionary process to keep pace with shifting environmental conditions. Such transitions are often referred to as climatic tipping points (TPs), associated with possibly irreversible changes in the state of the system. The term TP was originally introduced in social sciences¹¹ and made popular more recently by Gladwell¹². The study of TPs has recently gained broad interest and perspective in Earth and Environmental sciences, especially with regard to the future of our societies under the present climate warming scenarios¹³⁻¹⁶. The term Tipping Elements (TE) was introduced in¹⁵,¹⁷, and subsequently adopted by other authors¹⁴,¹⁸ to characterize the specific components of the Earth System that are likely to experience a TP in the near future as a result of the ongoing climate crisis¹⁶,¹⁹⁻²¹. Recently, the concept of TP has been used to define,

<footnote n="1">¹Géosciences Montpellier, Université Montpellier, Montpellier, France. ²Institute of Physics-CSE, Division of Geochronology and Environmental Isotopes, Silesian University of Technology, Gliwice, Poland. ³Lamont Doherty Earth Observatory, Columbia University, Palisades, NY, USA. ⁴Ecole Normale Supérieure—Paris Sciences et Lettres, Laboratoire de Météorologie Dynamique, Paris, France. ⁵Department of Mathematics and Statistics, University of Reading, Reading, UK. ⁶Centre for the Mathematics of Planet Earth, University of Reading, Reading, UK. ⁷School of Systems Science, Beijing Normal University, Beijing, People's Republic of China. ⁸email: denis-didier.rousseau@umontpellier.fr</footnote>
in turn, rapid societal changes that might lead to potentially positive impacts in terms of addressing the climate crisis<sup>22,23</sup>

Here, we want to investigate critical transitions in the Earth climate history by looking at two key high-resolution datasets that show evidence of abrupt transitions. The first dataset is the CENOGRID benthic δ¹⁸O and δ¹³C record corresponding to the compilation of 14 marine records over the past 66 Myr, from the Cretaceous-Paleogene (K-Pg) extinction event till present<sup>24</sup>. The second dataset comprises the North Atlantic U1308 benthic δ¹⁸O, δ¹³C and δ¹⁸O bulk carbonate time series covering the past 3.3 Myr<sup>25</sup>.

While visual evidences of abrupt transitions have already been discussed for these datasets, we wish to identify key abrupt thresholds by applying the recurrence quantification analysis (RQA) to each individual univariate time series and supplementing it with the Kolmogorov–Smirnov (KS) test<sup>26</sup>, see<sup>27</sup> and discussion below. Then, the selected transitions are discussed in the context of the Earth climate history allowing the definition of dynamical succession of abrupt transitions. Such transitions are then interpreted taking into account the evolution of key climate factors such as CO₂ concentration, average global sea level, and depth of the carbonate compensation.

The existence of TPs is intimately related to the multistability properties of the climate system, which have long been recognised in different contexts; see e.g.<sup>28–31</sup> and discussion in<sup>32,33</sup>. The multistability of the climate system comes the presence of more than one possible climate states for a given set of boundary conditions<sup>34</sup>. While earlier analyses have mostly evidenced the possibility of bistable behaviour, multistability can indeed include multiple competing states<sup>35–38</sup>. Recently, it has been proposed that the metastability properties of the climate system can be understood by interpreting the climate evolution as a diffusion process taking place in an effective dynamical landscape<sup>37,39,40</sup>defined by the the Graham's quasi-potential<sup>41</sup>. The local minima of the quasi-potential indicate the competing metastable states, with the transitions between such states occurring preferentially through the saddles of the quasi-potential. Such a viewpoint mirrors earlier proposals for interpreting biological evolution, namely the Waddington's epigenetic landscape<sup>42–46</sup>, see also relevant literature associated with synthetic evolution models like Tangled Nature<sup>49,50</sup>, and foresees the climatic transitions associated with the TPs as a manifestation of a dynamics characterized by punctuated equilibria<sup>47–50</sup>, The interplay between periods of stasis and rare transitions between competing metastable states seem to reflect the fact that climate fluctuations behave according to the dominance of stabilizing vs destabilizing feedbacks when considering short vs long time scales, respectively<sup>51</sup>. We will see that in most cases the RQA applied to the CENOGRID δ¹⁸O dataset flags TPs that disagree in terms of dating from those derived from the δ¹³C dataset. This is unsurprising because the two proxy data are sensitive to vastly different climatic processes. Nonetheless, the candidate TPs from both datasets come hand in hand with saddles of the bidimensional quasi-potential estimated from the bivariate time series. This indicates consistency between separate ways of detecting TPs. Additionally, TPs featuring faster characteristic time scales are associated with smaller-scale decoration of the quasi-potential, in agreement with what conjectured in<sup>37</sup>. The analysis of the benthic δ¹⁸O and δ¹³C suggests that the evolution of the climate in Cenozoic is characterized by a hierarchy of TPs due to an underlying multiscale quasi-potential.

## Results

**Detecting critical transitions of the past 66 Myr–3 Myr history of the earth climate.** The augmented KS test of the benthic δ¹⁸O record of the past 66 Ma identifies six major abrupt transitions corresponding to two major warming events at about 58 Ma and 56 Ma, followed by four major coolings at 47 Ma, 34 Ma, 14 Ma and 2.8 Ma respectively (Fig. 1A). The competing metastable states associated with these transitions feature rather long temporal persistence (a long time-window of 1–4 Ma is used, see Suppl. Mat.). These events are classical ones described from the literature<sup>52</sup>, where the first two transitions led to warmer conditions, while the latter four led to colder conditions. The same transitions are identified by employing the recurrence plot (RP) and recurrence rate (RR) analyses<sup>27,53</sup>, which also identify three more events, occurring at around 63 Ma, 40 Ma and 9.7 Ma, see Fig. 1C, Suppl. Tab. 1. We have chronologically labeled these TP<sub>O</sub>1 to TP<sub>O</sub>9, where the lower index refers to the used proxy data. As shown in Fig. 1B, TP<sub>O</sub>6 separates the climate variability in two separate macroclusters prior and after 34 Ma, the well-known Eocene–Oligocene Transition (EOT)<sup>54</sup>, which is a key step in the Cenozoic climate history and is associated with a major extinction event<sup>25,56</sup>.

The older macrocluster corresponds to very warm climate conditions and, using Marwan et al.'s nomenclature<sup>53</sup>, a disrupted variability. The average global temperature was estimated to be between 8 °C and 16 °C above the present day one, with no apparent presence of any major continental ice bodies<sup>24</sup> (Fig. 1A). Five major transitions are found within this period. They include the 2 major abrupt warmings at 58 Ma (TP<sub>O</sub>2) and 56 Ma (TP<sub>O</sub>3). These two transitions are thresholds towards much warmer oceanic deep water characterizing the first late Paleocene-Eocene hyperthermal<sup>57,58</sup> and Paleocene-Eocene Thermal Maximum (PETM)<sup>59,60</sup> respectively. The third one at about 47 Ma (TP<sub>O</sub>3) represents a transition towards cooler deep waters and named the Early-Middle Eocene cooling<sup>58</sup>. The last one at about 40 Ma (TP<sub>O</sub>4) inaugurates the period of continuous cooling due to the decrease of CO₂ concentrations that eventually leads to the TP<sub>O</sub>6 event<sup>61</sup>, which is associated with the build-up of the Antarctica ice sheets. TP<sub>O</sub>6 separates two fundamentally different modes of operation of the climate system, corresponding, as discussed later in Sect. 4, to two separate minima of the quasi—potential. From 34 Ma to present day the records feature more positive values of benthic δ¹⁸O associated with prevailing colder climate conditions. After TP<sub>O</sub>6, the climate featured mostly stationary conditions with a slight warming until the Middle Miocene Climate Transition (TP<sub>O</sub>7) that occurred around 14 Ma<sup>62,63</sup>. The last major transition (TP<sub>O</sub>9) occurred around 2.8 Ma, leading to the Pleistocene and the onset of glaciations in the Northern Hemisphere.

What discussed above provides a more systematic and robust counterpart of the analysis of TPs presented in<sup>24</sup>. In general, we prioritize the information coming from the δ¹⁸O record because of its very strong link to the Earth's temperature. Nonetheless, the analysis of the Benthic δ¹³C record performed along the same lines as above provides separate pieces of information on the critical transitions of the Earth's Climate. Benthic δ¹³C values
<figure><img src="image_4.png" /><figcaption>**Figure 1.** KS test and Recurrence Quantification Analysis (RQA) of CENOGRID benthic $\delta^{18}$O. (a) Time series in Ma BP with difference of the reconstructed and present Mean Global Temperature in pink). KS test identifying abrupt transitions towards warmer conditions in red and cooler or colder conditions in blue; (b) Recurrence plot (RP) with identification of the main two clusters prior and after 34 Ma. The main abrupt transitions identified are highlighted by red circles, and (c) Recurrence rate (RR). The pink crosses and vertical green lines indicate the abrupt transitions (TP) detected by the RQA. CENOGRID benthic $\delta^{18}$O data are from Westerhold et al. 24.</figcaption></figure>

characterize deep-water ventilation with high $\delta^{13}$C values in regions close to deep-water formation area. The KS analysis performed over a time window of 1–4 Myr individuates 14 TPs, with the RP suggesting an additional one,
located at around 34 Ma and associated with the EOT. We refer to these 15 TPs associated with the δ¹³C record as TPc 1 to 15; see Suppl. Figure 1, Suppl. Tab. 1. The interval 56.15 Ma–7.15 Ma is well characterized, showing some subclusters distributed around 34 Ma. The 56.15 Ma date groups δ¹³C values above 1‰, at the base of the record, while 7.15 Ma gathers the negative δ¹³C values which mainly occurs at the top of the record. The two periods are characterized by very different climatic conditions. The earlier climate regime features more input of carbon in the ocean while the later climate, instead, is characterized by higher presence of carbon in the atmosphere.

**Discussion of the detected critical transitions.** We next want to analyze TP₀ 1–9 in relation to different reconstructed paleoclimatic data, namely the global mean sea level (GMSL), the Pacific carbonate compensation depth (CCD), and CO₂ concentration (Fig. 2). Using benthic foraminifera δ¹⁸O and Mg/Ca records from high-resolution Pacific cores which were not included in the CENOGRID compilation, Miller et al.⁵⁸ have reconstructed variations in the GMSL over the past 66 Ma. By measuring the carbonate content in Pacific sediment cores and applying transfer functions, Pálike et al.⁶⁴ have generated a detailed Cenozoic record of the Pacific CCD, which denotes the depth below which carbonates dissolve. Finally, by compiling estimates from various proxies including foraminifera δ¹³C, boron isotopes, stomata, paleosols, Beerling & Royer⁶⁵ have produced a comprehensive Cenozoic record of the CO₂ concentration. The signature of the TP₀₆ is evident in the three records, corresponding to an abrupt decrease of the GMSL by about 70 m and of the CCD by around 1000 m. Additionally, TP₀₆ marks the start of a progressive decrease in CO₂ levels, from approximately 750 ppm to values of the order of 280 ppm.

Along the lines of Westerhold et al.²⁴, one can identify four competing states, "Warmhouse" (66 Ma-TP₀₁ and TP₀₃-TP₀₅) and "Hothouse" (TP₀₁-TP₀₃) climates in the earlier, warmer period, followed by the "Coolhouse" (TP₀₅-TP₀₆) and "Icehouse" climates (TP₀₆ to present); see Fig. 1. The first two states alternated in a warm-hot-warm sequence under extremely high CO₂ concentrations⁶⁵ as compared to those measured over the past 800 Kyr in the Antarctic ice cores (Fig. 2C). Before 34 Ma, one finds substantially larger values for GMSL, CCD depth, and CO₂ concentration whose average values are: +38 ± 15 m, 4600 ± 150 m, and 630 ± 300 ppmv respectively (Suppl. Tab. 2). Note that the CO₂ concentration is much higher than the levels observed over the past 800 Kyr in the Antarctic ice cores, which represent the reference states for the potential scenarios of climate warming outlined by the IPCC⁶⁷. Conversely, from 34 Ma until present, the Earth experienced much lower CO₂ concentrations and GMSL, thus generating the classical climate trend towards the recent ice-age conditions⁸,²⁴,⁵² (Fig. 2A,C). Indeed, the last 34 Myr show average values in GMSL, CO₂ concentration, and CCD depth of –3.5 ± 13 m, 330 ± 160 ppmv, and 3500 ± 400 m, respectively (Suppl. Tab. 2), which are much lower than in the older interval. This second set of means has underestimated values since it is based on data from the Miller et al.⁵⁸ dataset, which ends at 0.9 Ma,

These GMSL, CCD and CO₂ reconstructions show key transitions that fit with the CENOGRID thresholds deduced from the KS and RR analysis of the benthic δ¹⁸O as they signal an increase or a decrease in the global mean sea level corresponding roughly to warming or cooling episodes of the Earth history or strong variations in the concentration of atmospheric CO₂. The variations observed in CO₂, CCD and GMSL at the 9 identified TPs from the benthic δ¹⁸O record indicate heterogeneous characteristics prior to the TP6 major threshold (Suppl. Tab. 3). On the contrary, more homogenous features are noticed in the three climate proxies after TP₀₆, translating the occurrence of major reorganizations in the climate system, which become interesting to test at shorter timescales; see discussion below.

## Quasi-potential landscape and critical transitions

As discussed in detail in the Methods section, taking inspiration from the use of the Waddington epigenetic landscape to describe evolution⁴²⁻⁴⁶ and from the theory of punctuated equilibrium⁴⁹,⁵⁰, it has been proposed in³⁷,³⁹,⁴⁰ to study the global stability properties of the climate system by introducing an effective quasi-potential⁴¹, which generalizes the classical energy landscape framework often used to study metastable stochastic system. The quasi-potential formalism allows one to study general non-equilibrium system and to associate local maxima of the probability distribution function (pdf) of the system stable states. Additionally, saddles of the pdf are associated with Melancholia (M) states³⁴, which are unstable states living in the boundary between different basins of attraction. In the weak-noise limit, noise-induced transitions between competing stable states are expected to go through such M states. What we have done here is to construct the empirical bivariate pdf of climate system in the projected (δ¹³C, δ¹⁸O) and check whether TP₀ 1–9 indicated in Fig. 1 correspond, at least approximately, to saddles of the pdf. We find—see Fig. 3A—that, indeed, this is the case for TP₀s 1, 2, 4, 5, 6, 8, 9, whereas no agreement is found for TP₀s 3 and 7, which seem to take place in regions where the density is very small and very large, respectively. A very similar version of the bivariate pdf shown here had already been used in Westerhold et al.²⁴ to identify the Icehouse, the Coolhouse, the Warmhouse, and the Hothouse states, as groupings of nearby peaks corresponding to qualitatively similar climates.

A major improvement we propose here is to identify the transitions between such states, as well as less obvious transitions, by looking at the saddle points. TP₀₆ clearly emerges as the most important transitions, as it basically breaks the pdf into two separate parts. Note that the transitions shown in Fig. 1 are associated, because of the choice of a long time-window, with events that occur over long periods and that lead to persistent changes in the state of the system.

By combining the recurrence analysis with the quasi-potential formalism we can extract further relevant information. It is natural to ask ourselves what happens if, instead, we consider a catalogue of transitions for the δ¹⁸O record that are detected by considering shorter time windows (0.25 Ma–1 Ma) in the KS procedure. One finds—see Fig. 3B–11 of such transitions (see Suppl. Mat.). Once we report such fast TP₀s (FTP₀s) into the empirical bivariate pdf of the climate system in the projected (δ¹³C, δ¹⁸O) space, we find that they correspond to finer and smaller structures of the pdf as opposed to the case of the TPs. Hence, events that are associated with
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Variation through time of three main climate factors and comparison with the identified abrupt transitions (TPO) in the CENOGRID benthic δ¹⁸O. (a) Global Mean Sea Level in meters from Miller et al.⁵⁸. Identification of particular warm and of glaciation events. The Laurentide, GIW-WAIS and Ice free lines are from Miller et al.⁵⁸; (b) Carbonate Compensation Depth (CCD) in meters from Pálike et al.⁶⁴. The purple circles identify the TPs on this record; (c) Estimate of the CO₂ concentration in parts per million volume (ppmv) from Beerling & Royer⁶⁵. The Antarctica glaciation threshold at 750 ppmv and the NH glaciation threshold at 280 ppmv lines respectively are from DeConto et al.⁶⁶.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Probability density of the climate system in the projected CENOGRID benthic $\delta^{18}$O and $\delta^{13}$C space. (a) Chronologically ordered $TP_{OS}$ (diamonds) selected according to the KS methodology for $\delta^{18}$O with time window 1–4 My are shown. The two extra $TP_{OS}$ found via RP are indicated with a *. (b) Chronologically ordered fast $TP_{OS}$ (FTP$_{OS}$, diamonds) selected according to the KS methodology for $\delta^{18}$O with time window 0.25–1 Myr are shown. (c) Same as a), but for the $\delta^{13}$C record. The approximate timing of the TPs is indicated (rounded to .01 My). The 5 Ky-long portions of trajectories before and after each TP are also plotted.</figcaption></figure>

faster time scales are associated with smaller jumps between secondary maxima in the pdf belong to a hierarchically lower rung than those occurring over longer time scales. This seems to support the proposal made in³⁷.

As additional step, we repeat the same analysis leading to Fig. 3A by using, instead, the $\delta^{13}$C record, which features, as mentioned above, a total of 15 $TP_{CS}$. Figure 3C shows clearly that the $TP_{CS}$ are, for the most part, saddles that had not been flagged by the $TP_{OS}$. We have evidence that the same saddle is crossed more than once in a back-and-forth fashion few millions of years apart (e.g. $TP_{C4}$ & $TP_{C5}$; $TP_{C6}$ & $TP_{C7}$; $TP_{C12}$ & $TP_{C13}$), which supports the dynamical interpretation discussed here. Comparing Fig. 3A and C, one discovers that the
same saddle is responsible for TP<sub>O</sub>8 and for TP<sub>C</sub>15 even if the dating is different. Similarly, the same saddle is responsible for TP<sub>O</sub>5 and TP<sub>C</sub>6 & TP<sub>C</sub>7. Two key climatic features that appear as TPs in both records: the PETM is captured by TP<sub>O</sub>3 and TP<sub>C</sub>3 in both records, while TP<sub>O</sub>6 and TP<sub>C</sub>8 represent the EOT.

## The recent past

The past 3.3 Myr record from North Atlantic core U1308 can be considered as a blow-up of the CENOGRID dataset (Fig. 4). As previously mentioned, the last 3.3 Myr have been defined as an Icehouse climate state, with the appearance, development, and variations of the NHIS²⁴, whilst the Antarctic ice sheets had already mostly reached their maximal expansion. The variations in the deep-water temperature, as expressed by the benthic δ¹⁸O, are interpreted as an indicator of the continental ice volume with clear interglacial-glacial successions⁶⁸⁻⁷⁰. The Icehouse state is characterized by a change of the interplay between benthic δ¹³C and δ¹⁸O, which corresponds

<figure><img src="image_7.png" /><figcaption>**Figure 4.** RQA of U1308 benthic δ¹⁸O. (a) Time series in ka BP; (b) RP; and (c) RR. Pink crosses and green lines as in Fig. 1. TPo9-10 and RTPo1-4 abrupt transitions are identified from the RR. U1308 benthic δ¹⁸O data are from Hodell & Channell²⁵.</figcaption></figure>
to a new relationship between the carbon cycle and climate⁷¹. Indeed, one finds a very strong correlation between the two records (Pearson's coefficient being approximately -0.6). The correlation mainly results from the fact that the time series have approximately a common quasi-periodic behavior due to amplified response to the astronomical forcing as dictated by the Milankovich theory. Note that the presence of such almost regular resonant oscillations makes the use of the quasi-potential framework not particularly useful for describing the dynamics of the system, so that we will not pursue this approach for the analysis of the Quaternary records.

The KS augmented test and the RQ of the benthic $\delta^{18}$O agree in identifying six abrupt transitions dated at approximately 2.93 Ma, at 2.52 Ma, at 1.51 Ma, at 1.25 Ma, at 0.61 Ma and at 0.35 Ma (Fig. 4, Suppl. Tab. 1). They characterize the dynamics of North Hemisphere ice sheets (elevation and spatial expansion) and agree with the classical transitions characterizing the Marine Isotope Stages (MIS) as already observed in numerous records covering the same interval (see Rousseau et al.⁷² and references therein).

The first two transitions broadly match the previously discussed CENOGRID's $TP_{O9}$ associated with the onset of the Pleistocene (note that the interval between the two is much smaller than the resolution needed to separate to $TP_{O8}$) and are followed by four more recent TPs associated with the benthic $\delta^{18}$O record ($RTP_{O8}$). There is clear evidence of the Mid-Pleistocene (MPT) critical transition, between 1.25 Ma and 0.8 Ma, during which a shift occurred from climate cycles dominated by a 40-kyr periodicity (due to obliquity) to 100-Kyr periodicity (due to eccentricity) dominated ones⁷³⁻⁷⁷. The 1.25 Ma date is particularly significant, since it is followed by an Increase in the amplitude of glacial-interglacial fluctuations.

A complementary RQA of the $\delta^{18}$O bulk carbonate record from U1308, which characterizes episodes of iceberg calving into the North Atlantic Ocean IRD released into the North Atlantic Ocean²⁵, and therefore illustrates the dynamics of the Northern Hemisphere ice sheets (NHIS), yields similar dates to those obtained for the benthic $\delta^{18}$O record (see Suppl. Figure 3, Suppl. Tab. 1)⁷². Indeed, one finds abrupt transitions at 2.75 Ma, at 1.5 Ma, at 1.25 Ma, at 0.9 Ma and 0.65 Ma. Finally, as opposed to the CENOGRID $\delta^{13}$C RP, U1308 $\delta^{13}$C RP shows a drifting pattern similar to that of benthic $\delta^{18}$O, with only 2 key transitions at 2.52 Ma and 0.48 Ma (see Suppl. Figure 3). Note that the 0.48 Ma transition does not have any equivalent in the benthic $\delta^{18}$O records.

## Discussion

Studying the same CENOGRID dataset, Boettner et al.⁷⁸ identified 9 geological transitions at respectively 62.1 Ma, 55.9 Ma, 33.9 Ma, 23.2 Ma, 13.8 Ma, 10.8 Ma and 7.6 Ma. Four of them are indeed identical to those determined in the present study: 62.1 Ma, 55.9 Ma, 33.9 Ma and 13.8 Ma, the first two being preceded by a significant early warning signal⁷⁸, which is instead absent in the case for the EOT key transition.

Based on the results of both the RQA and the KS test of the $\delta^{18}$O and $\delta^{13}$C time series considered in this study, and the bivariate analysis performed using the framework of the quasi-potential theory, we propose a succession of critical transitions as described in Fig. 5. The critical transitions $TP_{O1}$ to $TP_{O9}$ shaped the Earth climate towards the onset and development of the Southern ice sheets and the later build-up of the NHIS. $TP_{O9}$ is followed by four more recent RTPs during the Quaternary, which steered the evolution of the ice sheets and of the climate as a whole until the present day. The climatic evolution during the Cenozoic until about 3 Ma seems to conform to a punctuated equilibrium framework, where the $TP_{O8}$s are associated with rapid transitions between rather different metastable modes of operation of the climate system. In particular, a key step that separates two rather diverse sets of climatic states occurred around 34 Ma at the EOT ($TP_{O6}$). Without the major drop in GMSL, in $CO_2$ concentrations, and in CCD, the Earth climate could have been different. However, after $TP_{O6}$, the Earth climate entered new dynamical regimes marked by much lower $CO_2$ concentrations, a lower GMSL, and a lower CCD. The remodeling of oceanic basins and mountain uplifts changed the marine and atmospheric circulations patterns, which played a crucial role in initiating and shaping the development of the NHIS.

Interestingly, the analysis of the $\delta^{13}$C time series identifies a different set of critical transitions, with the exception of the PETM ($TP_{O2}$ and $TP_{C3}$) and the EOT ($TP_{O6}$ and $TP_{C8}$), which are identified for both analyzed proxies. Looking into the bivariate pdf in the projected ($\delta^{13}$C, $\delta^{18}$O) space allows a better understanding of the nature and the origin of the TPs separately detected by studying the recurrence properties of the univariate time series. Indeed, we are in most cases able to associate both the $TP_{O8}$s and the $TP_{C8}$s to transitions across saddles of the effective quasi-potential. Clearly, some critical transitions might be more easily detectable when examining one time series rather than the other, because different proxies might be more sensitive to the active climatic processes, yet the approach taken here allows placing all TPs within a common ground. Additionally, TPs associated with faster processes and occurring between slower TPs correspond to transitions across smaller-scale saddles of the quasi-potential, hence revealing an analogy between the hierarchy of TPs and the multiscale nature of the quasi-potential.

More recently, variations in the extent and volume of the NHIS have contributed to the occurrence of the millennial variability marked by the Bond cycles, which have been best described during the last climate cycle. However, the onset of these cycles has been proposed to date back to 0.9 Ma⁷². Human activity is now rapidly pushing the Earth system towards the limits of its safe operating space associated with the occurrence of TPs; see also Rothman⁹ for a geochemical perspective focusing on the ongoing perturbation imposed on the marine carbon cycle. This concern is supported by actual observations impacting numerous tipping elements (see Lenton¹⁴) through very drastic tipping cascades¹⁶. This paves the way for a possible upcoming major transition, which might lead us to climate conditions that are fundamentally different from what has been observed in the recent or more distant past, and, at the very least, could bring us into a climate state with a much reduced or absent NHIS³⁸. This potential major transition, leading to de facto irreversible changes for the climate and the biosphere, could become a boundary between the Cenozoic icehouse and a new, warmer and radically different climate state compared to Pleistocene conditions.
<figure><img src="image_8.png" /></figure>

**Figure 5.** Evolution of the Earth Climate history among 2 different tipping dynamical landscapes and proposal for a potential third one. The first dynamical landscape, in light red, corresponds to the Hot-Warm House time interval. The second one, in light blue, represents the Cold-Ice house time interval. The third one, in light green, highlights the potential new dynamical landscape represented by the Anthropocene time interval. The different abrupt transitions identified in the present study are reported as TPos or RTPos to differentiate the major tipping points from the critical transitions characterizing transitions of lighter significance in the climate history. Various plate tectonic and ice sheet events are indicated and supported by maps of plate movements and North and South Hemisphere ice sheets. The Antarctica maps are from Pollard & DeConto<sup>79</sup>, Northern Hemisphere ice sheet maps are from Batchelor et al.<sup>80</sup>. The paleogeographic maps have been generated using the Ocean Drilling Stratigraphic Network (ODSN) plate tectonic reconstruction service: <https://www.odsn.de/odsn/services/paleo map/paleomap.html>. The red arrows on the tectonic maps indicate the key events that correspond to abrupt transition.

## Methods

**Recurrence plots and Kolmogorov-Smirnov test.** The Kolmogorov-Smirnov (KS) test is a robust method for accurately detecting discontinuities in a particular time series and is therefore a very precise way for determining the timing of abrupt transitions. Our method, described by Bagniewski et al.<sup>26</sup> is modified from the two-sample KS test and has been successfully applied to various geological time series<sup>24,25,69</sup>. The method uses the KS statistic, $D_{KS}$, to compare two sample distributions taken before and after a potential transition point within a sliding window. If the samples do not belong to the same continuous distribution, i.e., $D_{KS}$ is greater than a predefined threshold, a transition point is identified.

This classic two-sample KS test is augmented by additional criteria to refine the results and pinpoint significant abrupt transitions indicative of a true climatic shift. This involves discarding transitions below a rate-of-change threshold since smaller changes in the time series might be due to data errors or short-term variability within intervals shorter than the proxy record's sampling resolution. Furthermore, as the frequency with which the KS test detects transitions largely depends on the chosen window length, $D_{KS}$ is calculated for different window length, within a range that corresponds to the desired time scale at which a given paleorecord is to be investigated. The analysis starts by identifying transitions using the longest window, which has the largest sample size and thus carries greater statistical significance. Subsequently, the method incorporates transitions detected using shorter windows to capture transitions occurring on shorter time scales. For more details see Bagniewski et al.<sup>26</sup>, Rousseau et al.<sup>72</sup>.

To gain further insight into the climate story revealed by proxy records, we performed an analysis based on recurrence plots (RPs), first introduced by Eckmann<sup>81</sup> in the study of dynamical systems and later popularized in
the climate sciences by Marwan et al. ⁵³,⁸². The RP for a time series {$x_i$: $i = 1, ..., N$} is constructed as a square matrix in a Cartesian plane with one copy {$x_i$} of the series on the abscissa and another copy {$x_j$} on the ordinate, with both axes representing time. A dot is entered into a position (i, j) of the matrix when $x_j$ is sufficiently close to $x_i$. For the details—such as how “sufficiently close” is determined—we refer to Marwan et al. ⁵⁰.

The square matrix of dots that is the visual result of RP exhibits a characteristic pattern of vertical and horizontal lines, indicating recurrences. These lines sometimes form clusters that represent specific periodic patterns. Eckmann⁸¹ distinguished between large-scale *typology* and small-scale *texture* in the interpretation of RPs. The most interesting typologies in RP applications are associated with recurrent patterns that are not strictly periodic and are thus challenging to detect by purely spectral approaches to time series analysis. RPs offer an important advantage by enabling a visual identification of nonlinear relationships and dynamic patterns that characterize high-dimensional systems subject to time-dependent forcing, such as the climate system. However, in periodically forced systems, the recurrence structure is dominated by repetitive, periodic patterns rather. This may pose challenges when attempting to identify meaningful transitions or regime shifts in climate at time scales strongly affected by orbital forcing.

Marwan et al.⁸² extensively discussed the various measures used to objectively quantify the typologies observed in RPs. Collectively known as recurrence quantification analysis (RQA), these measures include the Recurrence Rate (RR), which represents the probability of a specific state recurring within a given time interval. The RR is calculated by determining the density of recurrence points along the diagonal of the RP within a sliding window. Since low RR values correspond to an unstable behavior of the system, the minima of the RR may be used to identify abrupt transitions. Here we follow the approach by Bagniewski et al.²⁶ and select local minima based on their prominence.

**Quasi-potential, Melancholia states, and critical transitions.** Traditionally, tipping points are schematically represented as being associated with the bifurcation occurring for a system described by a one-dimensional effective potential when a change in the value of a certain parameter leads to a change in the number of stable equilibriums. Hence, conditions describing the nearing of a tipping point can be related to the presence of slower decay of correlations (critical slowing down). This viewpoint, while attractive, suffers from many mathematical issues due to the fact that the true dynamics of the system occurs in a possibly very high dimensional space. Tantet et al.⁸³,⁸⁴ and Santos Gutierrez & Lucarini⁸⁵ have introduced a mathematically rigorous framework for the occurrence of tipping points that clarifies the link between rate of decay of correlations, sensitivity of the system to perturbations, and robustness of the unperturbed dynamics.

Here, we wish to take a different angle on the problem. Instead of focusing on the individual tipping points, we attempt to capture the global stability properties of the system. We take inspiration from the application of the Waddington epigenetic landscape to describe morphological evolution⁴²⁻⁴⁶ and from the theory of punctuated equilibrium⁴⁹,⁵⁰, which associates periods of stasis (characterized by relatively stable morphology) with the emergence of new species through abrupt changes (named cladogenesis) from a previous species. Lucarini & Bodai³⁹,⁴⁰ proposed describing the global properties of the climate system using the formalism of quasi-potential⁴¹. Roughly speaking, assuming that the system lives in $\mathbb{R}^N$, and its dynamics is described by a stochastic differential equation the probability that its state is within the volume $d\vec{x}$ around the point $\vec{x} \in \mathbb{R}^N$ is given by $P(\vec{x}, d\vec{x}) = \rho(\vec{x})d\vec{x}$ where $\rho(\vec{x}) \approx e^{-\frac{2\Phi(\vec{x})}{\varepsilon^2}}$ is the probability distribution function (pdf) and $\Phi(\vec{x})$ is the quasipotential.

The function $\Phi(\vec{x})$ depends in a nontrivial way on the drift term and noise law defining the stochastic differential equation. This setting generalizes the classical energy landscape and applies to a fairly large class of stochastic dynamical systems. One can see the dynamics of the system as being driven towards lower values $\Phi(\vec{x})$ (plus an extra rotational effect that is typical of non-equilibrium systems), while the stochastic forcing noise makes sure that the system is erratically pushed around. Hence, the minima of $\Phi(\vec{x})$ correspond to local maxima of the pdf, and the saddles (which coincide for both $\Phi(\vec{x})$ and $\rho(\vec{x})$) coincide with the Melancholia (M) states³⁴. Such M states are unstable states of the system that live at the boundary between basins of attraction and are the gateways for the noise-induced transitions between competing stable states. Margazoglou et al.³⁷ applied this method to the investigation of the metastability properties of an intermediate complexity climate model and suggested that the presence of decorations of the quasi-potential at different scales could be interpreted as being associated with a hierarchy of tipping points. Indeed, passing near M states is intimately associated with the occurrence of critical transitions. Hence, the construction of the quasi-potential $\Phi(\vec{x})$ can be seen as the structural counterpart of the investigation of the time-evolution of the system and of its critical transitions. Zhou et al.⁸⁶ give a complete overview of different methods applied to perform such an analysis.

**Data and materials availability**

All data generated by the present study from the main text or the supplementary materials will be submitted to PANGAEA data repository. U1308 marine data are available at https://doi.org/10.1594/PANGAEA.871937 (Hodell and Channell, 2016b). CENOGRID data are available at https://doi.org/10.1594/PANGAEA.917503 (Westerhold, 2020).

Received: 8 December 2022; Accepted: 8 July 2023

Published online: 12 July 2023
## References

1. Dansgaard, W., Johnsen, S. J., Moller, J. & Langway, C. C. One thousand centuries of climatic record from camp century on the Greenland ice sheet. *Science* **166**, 377–381 (1969).

2. Dansgaard, W. et al. A new Greenland deep ice core. *Science* **218**, 1273–1277 (1982).

3. Broecker, W. Climatic change—are we on the brink of a pronounced global warming. *Science* **189**, 460–463 (1975).

4. Broecker, W. S. et al. Can the Greenland climatic jumps be identified in records from ocean and land? *Quat. Res.* **30**, 1–6 (1988).

5. Broecker, W. S. & Denton, G. The role of ocean-atmosphere reorganizations in glacial cycles. *Geochim. Cosmochim. Acta* **53**, 2465–2501 (1989).

6. Rasmussen, S. O. et al. A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: Refining and extending the INTIMATE event stratigraphy. *Quat. Sci. Rev.* **106**, 14–28 (2014).

7. Hoffman, P. F. et al. Snowball Earth climate dynamics and Cryogenian geology-geobiology. *Sci. Adv.* **3**, (2017).

8. Scotese, C. R., Song, H., Mills, B. J. W. & van der Meer, D. G. Phanerozoic paleotemperatures: The earth's changing climate during the last 540 million years. *Earth-Sci. Rev.* **215**, (2021).

9. Rothman, D. H. Thresholds of catastrophe in the Earth system. *Sci. Adv.* **3**, e1700906 (2017).

10. Newell, N. D. Crises in the history of life. *Sci. Am.* **208**, 76–95 (1963).

11. Grodzins, M. Metropolitan segregation. *Sci. Am.* **197**, 33– (1957).

12. Gladwell, M. *The tipping point: how little things can make a big difference*. (Back Bay Books, 2002).

13. Lenton, T. et al. Tipping elements in the Earth's climate system. *Proc. Nat. Acad. Sci. U. S. A.* **105**, 1786–1793 (2008).

14. Lenton, T. Environmental Tipping Points. In: *Annual review of Environment and Resources* (eds. Gadgil, A. & Liverman, D.) vol. 38 1–29 (2013).

15. Lenton, T. et al. Climate tipping points—too risky to bet against. *Nature* **575**, 592–595 (2019).

16. Brovkin, V. et al. Past abrupt changes, tipping points and cascading impacts in the Earth system. *Nat. Geosci.* **14**, 550–558 (2021).

17. Lenton, T. & Williams, H. On the origin of planetary-scale tipping points. *Trends Ecol. Evol.* **28**, 380–382 (2013).

18. Steffen, W. et al. Trajectories of the earth system in the anthropocene. *Proc. Nat. Acad. Sci. U. S. A.* **115**, 8252–8259 (2018).

19. Kriegler, E., Hall, J., Held, H., Dawson, R. & Schellnhuber, H. Imprecise probability assessment of tipping points in the climate system. *Proc. Nat. Acad. Sci. U. S. A.* **106**, 5041–5046 (2009).

20. Wunderling, N., Donges, J., Kurths, J. & Winkelmann, R. Interacting tipping elements increase risk of climate domino effects under global warming. *Earth Syst. Dyn.* **12**, 601–619 (2021).

21. Boers, N. Observation-based early-warning signals for a collapse of the atlantic meridional overturning circulation (vol 11, pg 680, 2021). *Nat. Clim. Change* **11**, 1001–1001 (2021).

22. Otto, I. et al. Social tipping dynamics for stabilizing Earth's climate by 2050. *Proc. Nat. Acad. Sci. U. S. A.* **117**, 2354–2365 (2020).

23. Lenton, T. et al. Operationalising positive tipping points towards global sustainability. *Glob. Sustain.* **5**, (2022).

24. Westerhold, T. et al. An astronomically dated record of Earth's climate and its predictability over the last 66 million years. *Science* **369**, 1383–1387 (2020).

25. Hodell, D. A. & Channell, J. E. T. Mode transitions in Northern Hemisphere glaciation: Co-evolution of millennial and orbital variability in Quaternary climate. *Clim. Past* **12**, 1805–1828 (2016).

26. Bagniewski, W., Ghil, M. & Rousseau, D. D. Automatic detection of abrupt transitions in paleoclimate records. *Chaos* **31**, 113129 (2021).

27. Bagniewski, W., Rousseau, D.-D. & Ghil, M. The PaleoJump database for abrupt transitions in past climates. *Sci. Rep.* **13**, 4472 (2023).

28. Stommel, H. Thermohaline convection with 2 stable regimes of flow. *Tellus* **13**, 224–230 (1961).

29. Sellers, W. D. A global climatic model based on the energy balance of the earth-atmosphere system. *J. Appl. Meteorol. Climatol.* **8**, 392–400 (1969).

30. Budyko, M. Effect of solar radiation variations on climate of earth. *Tellus* **21**, 611 (1969).

31. Ghil, M. Climate stability for a sellers-type model. *J. Atmospheric Sci.* **33**, 3–20 (1976).

32. Ghil, M. A century of nonlinearity in the geosciences. *Earth Space Sci.* **6**, 1007–1042 (2019).

33. Ghil, M. & Lucarini, V. The physics of climate variability and climate change. *Rev. Mod. Phys.* **92**, (2020).

34. Lucarini, V. & Bodai, T. Edge states in the climate system: Exploring global instabilities and critical transitions. *Nonlinearity* **30**, R32–R66 (2017).

35. Lewis, J., Weaver, A. & Eby, M. Snowball versus slushball Earth: Dynamic versus nondynamic sea ice? *J. Geophys. Res.-Oceans* **112**, (2007).

36. Abbot, D., Voigt, A. & Koll, D. The Jormungand global climate state and implications for Neoproterozoic glaciations. *J. Geophys. Res.-Atmospheres* **116**, (2011).

37. Margazoglou, G., Grafke, T., Laio, A. & Lucarini, V. Dynamical landscape and multistability of a climate model. *Proc. R. Soc. A-Math. Phys. Eng. Sci.* **477**, (2021).

38. Ragon, C. et al. Robustness of competing climatic states. *J. Clim.* **35**, 2769–2784 (2022).

39. Lucarini, V. & Bodai, T. Transitions across melancholia states in a climate model: Reconciling the deterministic and stochastic points of view. *Phys. Rev. Lett.* **122**, (2019).

40. Lucarini, V. & Bodai, T. Global stability properties of the climate: Melancholia states, invariant measures, and phase transitions. *Nonlinearity* **33**, R59–R92 (2020).

41. Graham, R. Macroscopic potentials, bifurcations and noise in dissipative systems. In: *Fluctuations and stochastic phenomena in condensed matter* 1–34 (Springer, 1987).

42. Waddington, C. Canalization of development and the inheritance of acquired characters. *Nature* **150**, 563–565 (1942).

43. Waddington, C. H. *The strategy of the genes*. (Allen & Unwin, 1957).

44. Goldberg, A. D., Allis, C. D. & Bernstein, E. Epigenetics: A landscape takes shape. *Cell* **128**, 635–638 (2007).

45. Baedke, J. The epigenetic landscape in the course of time: Conrad Hal Waddington's methodological impact on the life sciences. *Stud. Hist. Philos. Sci. Part C: Stud. Hist. Philos. Biol. Biomed. Sci.* **44**, 756–773 (2013).

46. Allen, M. Compelled by the diagram: Thinking through CH Waddington's epigenetic landscape. *Contemporaneity* **4**, 119 (2015).

47. Christensen, K., de Collobiano, S., Hall, M. & Jensen, H. Tangled nature: A model of evolutionary ecology. *J. Theor. Biol.* **216**, 73–84 (2002).

48. Jensen, H. J. Tangled Nature: A model of emergent structure and temporal mode among co-evolving agents. *Eur. J. Phys.* **40**, (2019).

49. Eldredge, N. & Gould, S. J. Punctuated equilibria: An alternative to Phyletic Gradualism. In: *Models in Paleobiology* (ed. Schopf, T. J. M.) 82–243 (Freeman Cooper, 1972).

50. Gould, S. J. & Eldredge, N. Punctuated equilibria: The tempo and mode of evolution reconsidered. *Paleobiology* **3**, 115–151 (1977).

51. Arnscheidt, C. W. & Rothman, D. H. Presence or absence of stabilizing Earth system feedbacks on different time scales. *Sci. Adv.* **8**, eadc 9241 (2022).

52. Zachos, J., Pagani, M., Sloan, L., Thomas, E. & Billups, K. Trends, rhythms, and aberrations in global climate 65 Ma to present. *Science* **292**, 686–693 (2001).
53. Marwan, N., Carmen Romano, M., Thiel, M. & Kurths, J. Recurrence plots for the analysis of complex systems. *Phys. Rep.* **438**, 237–329 (2007).

54. Zachos, J., Breza, J. & Wise, S. Early Oligocene ice-sheet expansion on Antractica—stable isotope and sedimentological evidence from Kerguelen Plateau. *Southern Indian-Ocean. Geology* **20**, 569–573 (1992).

55. Ivany, L., Patterson, W. & Lohmann, K. Cooler winters as a possible cause of mass extinctions at the eocene/oligocene boundary. *Nature* **407**, 887–890 (2000).

56. Sun, J. et al. Synchronous turnover of flora, fauna, and climate at the Eocene-Oligocene Boundary in Asia. *Sci. Rep.* **4**, (2014).

57. Thomas, E., Zachos, J. C. & Bralower, T. J. Deep-sea environments on a warm earth: latest Paleocene-early Eocene. in *Warm Climates in Earth History* (eds. Huber, B. T., Macleod, K. G. & Wing, S. L.) 132–160 (Cambridge University Press, 1999). https://doi.org/10.1017/CBO9780511564512.006.

58. Miller, K. et al. Cenozoic sea-level and cryospheric evolution from deep-sea geochemical and continental margin records. *Sci. Adv.* **6**, (2020).

59. Kennett, J. & Stott, L. Abrupt deep-sea warming, palaeoceanographic changes and benthic extinctions at the end of the Paleocene. *Nature* **353**, 225–229 (1991).

60. Dickens, G. R., Castillo, M. M. & Walker, J. C. A blast of gas in the latest Paleocene: Simulating first-order effects of massive dissociation of oceanic methane hydrate. *Geology* **25**(3), 259–262 (1997).

61. Lauretano, V. et al. Eocene to oligocene terrestrial southern hemisphere cooling caused by declining pCO₂. *Nat. Geosci.* **14**, 659–664 (2021).

62. Kennett, J. & Shackleton, N. Oxygen isotopic evidence for development of psychrosphere 38 Myr ago. *Nature* **260**, 513–515 (1976).

63. Shackleton, N. J. Oceanic carbon isotope constraints on oxygen and carbon dioxide in the cenozoic atmosphere. In: *The Carbon Cycle and Atmospheric CO₂: Natural Variations Archean to Present* 412–417 (American Geophysical Union (AGU), 1985). doi:https://doi.org/10.1029/GM032p0412.

64. Pálike, H. et al. A Cenozoic record of the equatorial Pacific carbonate compensation depth. *Nature* **488**, 609–614 (2012).

65. Beerling, D. & Royer, D. Convergent cenozoic CO₂ history. *Nat. Geosci.* **4**, 418–420 (2011).

66. DeConto, R. et al. Thresholds for cenozoic bipolar glaciation. *Nature* **455**, 652–U52 (2008).

67. IPCC. *Climate Change 2021: The physical science basis. Contribution of working group i to the sixth assessment report of the inter-governmental panel on climate change*. vol. In Press (Cambridge University Press, 2021).

68. Chappell, J. & Shackleton, N. J. Oxygen isotopes and sea-level. *Nature* **324**, 137–140 (1986).

69. Shackleton, N. J. The 100,000-year ice-age cycle identified and found to lag temperature, carbon dioxide, and orbital eccentricity. *Science* **289**, 1897–1902 (2000).

70. Elderfield, H. et al. Evolution of ocean temperature and ice volume through the mid-pleistocene climate transition. *Science* **337**, 704–709 (2012).

71. Turner, S. K. Pliocene switch in orbital-scale carbon cycle/climate dynamics. *Paleoceanography* **29**, 1256–1266 (2014).

72. Rousseau, D., Bagniewski, W. & Ghil, M. Abrupt climate changes and the astronomical theory: Are they related?. *Clim. Past* **18**, 249–271 (2022).

73. Shackleton, N. J. & Opdyke, N. D. Oxygen isotope and palaeomagnetic evidence for early Northern Hemisphere glaciation. *Nature* **270**, 216–223 (1977).

74. Pisias, N. G. & Moore, T. C. The evolution of pleistocene climate: A time-series approach. *Earth Planet. Sci. Lett.* **52**, 450–458 (1981).

75. Ruddiman, W. F., Raymo, M., Martinson, D. G., Clement, B. M. & Backman, J. Pleistocene evolution =: Northern Hemisphere ice sheets and North Atlantic Ocean. *Paleoceanography* **4**, 353–412 (1989).

76. Clark, P. U. & Pollard, D. Origin of the middle Pleistocene transition by ice sheet erosion of regolith. *Paleoceanography* **13**, 1–9 (1998).

77. Clark, P. U. et al. The middle Pleistocene transition: Characteristics, mechanisms, and implications for long-term changes in atmospheric pCO₂. *Quat. Sci. Rev.* **25**, 3150–3184 (2006).

78. Boettner, C., Klinghammer, G., Boers, N., Westerhold, T. & Marwan, N. Early-warning signals for Cenozoic climate transitions. *Quat. Sci. Rev.* **270**, (2021).

79. Pollard, D. & DeConto, R. Continuous simulations over the last 40 million years with a coupled Antarctic ice sheet-sediment model. *Palaeoclimatol., Palaeoecol., Palaeogeogr.* **537**, (2020).

80. Batchelor, C. L. et al. The configuration of northern hemisphere ice sheets through the Quaternary. *Nat. Commun.* **10**, (2019).

81. Eckmann, J. P., Kamphorst, S. O. & Ruelle, D. Recurrence plots of dynamical systems. *Europhys. Lett.* **4**, 973–977 (1987).

82. Marwan, N., Schinkel, S. & Kurths, J. Recurrence plots 25 years later—Gaining confidence in dynamical transitions. *EPL* **101**, (2013).

83. Tantet, A., Lucarini, V. & Dijkstra, H. Resonances in a chaotic attractor crisis of the Lorenz flow. *J. Stat. Phys.* **170**, 584–616 (2018).

84. Tantet, A., Lucarini, V., Lunkeit, F. & Dijkstra, H. Crisis of the chaotic attractor of a climate model: A transfer operator approach. *Nonlinearity* **31**, 2221–2251 (2018).

85. Santos Gutiérrez, M. & Lucarini, V. On some aspects of the response to stochastic and deterministic forcings. *J. Phys. A: Math. Theor.* **55**(42), 425002. https://doi.org/10.1088/1751-8121/ac90fd (2022).

86. Zhou, J. X., Aliyu, M. D. S., Aurell, E. & Huang, S. Quasi-potential landscape in complex multi-stable systems. *J. R. Soc. Interface* **9**, 3539–3553 (2012).

## Acknowledgements

We would like to thank our colleagues from Horizon 2020 TiPES project (grant No. 820970) and especially Michael Ghil for useful discussions about this study. This is LDEO contribution and TiPES contribution 192. VL acknowledges useful exchanges with T. Bodai, R. Boerner, R. Deeley, G. Margazoglou, C. Nesbitt, D. Rothman and L. Serdukova. We would like to thank the editor and reviewers' feedback, which has helped us to improve the clarity and quality of our work.

## Author contributions

Conceptualization: D.D.R.; Methodology, investigation and Visualization: D.D.R., V.L., W.B.; Writing—original draft: D.D.R.; Writing—review & editing: D.D.R., V.L., W.B.

## Funding

This research has been supported by the European Commission, Horizon 2020 Framework Programme (TiPES, grant no. 820970), by the Marie Cure ITN Critical Earth (Grant No. 956170) and from the EPSRC Project No. EP/T018178/1.
Competing interests

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-023-38454-6.

**Correspondence** and requests for materials should be addressed to D.-D.R.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023