{'fraction_non_alphanumeric': 0.033116264294790344, 'fraction_numerical': 0.017947903430749684, 'mean_word_length': 4.631708407871199, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Attractors and bifurcation diagrams in complex climate models

The climate is a complex non-equilibrium dynamical system that relaxes toward a steady state under the continuous input of solar radiation and dissipative mechanisms. The steady state is not necessarily unique. A useful tool to describe the possible steady states under different forcing is the bifurcation diagram, that reveals the regions of multi-stability, the position of tipping points, and the range of stability of each steady state. However, its construction is highly time consuming in climate models with a dynamical deep ocean, interactive ice sheets or carbon cycle, where the relaxation time becomes larger than thousand years. Using a coupled setup of MITgcm, we test two techniques with complementary advantages. The first is based on the introduction of random fluctuations in the forcing and permits to explore a wide part of phase space. The second reconstructs the stable branches and is more precise in finding the position of tipping points.


Introduction
The climate is a dynamical complex system that is fuelled by the incoming solar radiation and reaches a steady state under the effect of dissipation over a multitude of temporal and spatial scales (Ghil & Lucarini, 2020). Under a given forcing (such as astronomical Milankovitch cycles, or the increasing atmospheric CO 2 content due to volcanism or present-day anthropogenic emissions) the system can be driven out of the steady state. In such conditions, the system can reach critical thresholds (or tipping points) where its properties abruptly change, often in an irreversible way. Such behaviour can be illustrated by the so-called bifurcation diagram (BD), where a state variable (for example, the mean surface air temperature) is plotted as a function of the driving force, as schematically shown in Figure 1 (second row), together with the corresponding potential curve (first row) for the different types of tipping mechanism (Ashwin et al., 2012;Vanselow et al., 2022): 1) B-tipping, that occurs when the deterministic system dynamics reach a bifurcation; 2) N-tipping, when the internal variability increases up to exceed the height of a critical barrier separating two basins of attraction, so that the system can access another dynamical solution; 3) S-tipping, related to sudden shocks that induce the passage from one state to another, as happened when a huge asteroid presumably caused the Cretaceous-Tertiary extinction 66 million years ago; 4) R-tipping, when the rate of change of the forcing or some internal parameter exceed a critical value, allowing the system to cross, for example, a moving basin boundary (note, however, that Rtipping can be induced by different mechanisms and does not necessarily requires multistability (Vanselow et al., 2022) BDs are easily obtained using energy balance models (Budyko, 1969;Sellers, 1969;Ghil, 1976), the simplest in the hierarchy of climate models. In intermediate complexity models (Lucarini & Bódai, 2017) or low resolution general circulation models (Rose, 2015;Zhu & Rose, 2022), BDs can be still constructed with a reasonable computational cost. As the complexity of the model increases, the amount of CPU time needed to perform series of simulations that explore a huge range of driving force and initial conditions toward steady states becomes prohibitive. Indeed, the standard method requires convergence toward the steady state (or attractor) to obtain the BD, that means to perform simulations over time scales of the order of several times the relaxation time of the included climatic components (McGuffie & Henderson-Sellers, 2005), which can be 10 3 years for the deep ocean, 10 4 years for the carbon equilibration between atmosphere and ocean (Zhu & Rose, 2022), and even higher for dynamical ice sheets.
Using a general circulation model (thus, at the top of the hierarchy) in a simplified configuration (i. e. a planet entirely covered by the ocean), Ragon et al. (2022) obtained the bifurcation diagram using the standard method over time scales of thousand years, thus excluding nonlinear feedbacks on longer time scales such as those induced by ice sheets and vegetation. The MIT general circulation model  (MITgcm) was used in its coupled atmosphere-oceansea ice configuration with 2.8 • horizontal resolution, 15 levels in an ocean depth of 3 km and 5 pressure layers in the atmospheric column. The atmospheric module is based on SPEEDY (Molteni, 2003) that provides a rather realistic representation of the radiative scheme despite the coarse vertical resolution, with the advantage of requiring less computer resources than state-of-the-art atmospheric modules. On the other hand, the ocean dynamics is very accurate and this is crucial to include nonlinear feedbacks on millennial time scales. Using this setup, it is possible to run nearly 200 years in 1 day on clusters at the University of Geneva. The same MITgcm setup, with an additional land module, has been succesfully applied to investigate the ocean dynamics during the Jurassic (Brunetti et al., 2015), the present-day climate (Brunetti & Vérard, 2018), and the climatic oscillations in the Early Triassic.
Here we compare the standard technique to construct the BD with two additional methods which require lower computational costs. We will describe such techniques in Section 2, the resulting BDs in Section 3 and we draw our conclusions on pros and cons for each method, and future developments in Section 4.

Construction of Bifurcation Diagrams
We first describe the standard technique and the corresponding BD. Since, in this case, convergence to the dynamical attractors is required over time scales comparable to the longest feedback process in the setup (in our case, until deep ocean equilibrium is reached), this BD is used for comparisons with those obtained with less computationally expensive techniques that will be described in the following subsections.

Standard Technique
The standard method to construct the BD is based on the theory of dynamical systems. The phase space can be divided in basins of attraction, i. e the minimal invariant closed sets attracting an open set of initial conditions as time goes to infinity (Strogatz, 1994) for given values of the internal parameters (such viscosity, diffusion coefficients, albedo of different surfaces, . . . ) and external forcing. Note that under a constant forcing, the system is ergodic (Tél et al., 2020), i. e sufficiently long temporal averages correspond to large ensemble averages.
Thus, the first step is to consider a huge ensemble of initial conditions spanning a large part of phase space. Starting from such initial conditions, the system then evolves in time until the attractors are reached. Statistically steady-state conditions are realized within each attractor in the climate system when its mean annual surface energy balance F s , i. e. the sum of sensible, latent, net solar and longwave radiation fluxes at the surface, becomes nearly zero. In general, simulation runs over several times n the relaxation time t relax are needed to guarantee convergence on the attractor (Drótos et al., 2017). t relax depends on the nonlinear feedbacks implemented in the numerical simulation. In an aquaplanet, there are no ice sheets (characterised by a time scale of the order of 10 5 yr), and vegetation (10 2 -10 3 yr). An active carbon cycle between ocean and atmosphere (10 4 yr) is also excluded, thus the deep-ocean dynamics is the process with the largest relaxation time of the order of 10 3 yr. In our setup, it turns out to be between 500 and 2000 yr, depending on the attractor (see Figure 2b in Brunetti et al. (2019)).
Five steady climates have been found (Brunetti et al., 2019) under the same forcing, represented by the same amount of incoming solar radiation (S 0 = 342 W m −2 ) and atmospheric CO 2 content (fixed at 326 ppm): snowball (where ice covers the entire surface), waterbelt (where an ocean belt survives near the equator), cold state (with an ice cap extending to 43 degree latitude), warm state (with an ice cap comparable to the present one, up to 60 degree) and hot state (a planet without ice). Simulations are stopped when the surface energy imbalance F s becomes lower than 0.2 W m −2 in absolute value in Brunetti et al. (2019), corresponding to an ocean temperature drift dT o /dt = F s /(c p ρh) (Marshall and Plumb 2008, p. 229) lower than 0.058 • C per century, with c p = 4000 J K −1 kg −1 being the specific heat capacity, ρ = 1023 kg m −3 the seawater density, and h = 3000 m the ocean depth. Indeed, under such conditions, there is essentially no drift in the annual averages of the climatic observables.
The second step is to slightly vary the forcing, starting from each attractor, and let the system relax again (Drótos et al., 2017). In this way, stable branches of the BD can be constructed. In practise, in the aquaplanet simulations, Ragon et al. (2022) performed the following: the simulations are initialized from the five attractors found in Brunetti et al. (2019) at S 0 = 342 W m −2 with slightly different values of the incoming solar radiation and same CO 2 content until convergence. This procedure is repeated until a tip-ping point is reached, that is where a sudden shift of the state variable to a different attractor is observed.
Since the attractors are complex dynamical objects living in a high dimensional space (Brunetti et al., 2019;Falasca & Bracco, 2022), the projection of their invariant measure on a given state variable, as the surface air temperature T , is arbitrary (Faranda et al., 2019;Tél et al., 2020). It turns out that different quantities are maximised in each attractor, as found in Ragon et al. (2022), i. e. total precipitation and surface temperature in hot state, intensity of the Lorenz energy cycle in warm state, heat transport in cold state, available potential and kinetic energies in waterbelt, while snowball minimises all the above quantities. Thus, each climatic state would be better represented by a different projection. However, in the present study, the projection is performed, as commonly done in the literature, in terms of the surface air temperature which spans an interval of more than 70 • C from the snowball to the hot state, thus differentiating well all the attractors.
The BD obtained in such a way is shown in Figure 1 of Ragon et al. (2022), where it can be seen that the range of stability of the warm state is very small in comparison to the others, while the snowball is stable over all the range of forcing that has been explored, from 334 to 350 W m −2 . The positions of B-tipping are also evident, bracketing either end of the stable ranges for warm and cold states, and the cold end of waterbelt and hot state. Note that the warm state may become unstable when carbon exchanges between atmosphere and ocean are included, as shown in Zhu and Rose (2022) using a setup with horizontal resolution of 3.75 • .
Such BD has been obtained using a very computationally expensive technique. The total simulation time t tot can be estimated as t tot ∼ n(N +M ) t relax , where n ∼ 10 to reach convergence, N ∼ 40 is the number of initial conditions in the first step, M ∼ 40 is the number of points on the stable branches of the bifurcation diagram (see Figure 1 in Ragon et al. (2022)), giving approximately t tot ∼ 800 t relax ∼ 8 · 10 5 yr. Of course, the effective time can be optimized by launching N simulations in parallel with different initial conditions, and by considering separately each stable branch. Note, however, that the second term in t tot cannot be reduced further, since the initial conditions at one forcing value are the final state of the previous one.

Method I: Random Fluctuations in the Forcing
This method is based on the idea described in Lucarini and Bódai (2020); Margazoglou et al. (2021) of studying noise induced transitions between basins of attraction. In order to explore the entire phase space of the system, and not only the single basin of attraction allowed in nonlinear deterministic dynamics, random fluctuations of the incoming solar radiation are introduced around a given value (S 0 = 342 W/m 2 in our setup) at discrete times. Even if noise directly affects only a small fraction of degrees of freedom, it propagates to all, since scales are interconnected as well as dynamical and thermodynamical fields in the coupled climate model (Lucarini & Bódai, 2020).
If η is a random number taken from a normal distribution with zero mean and standard deviation σ, a new value of the incoming solar radiation S is prescribed as S = S 0 (1+ η) at regular temporal intervals ∆t 1 . We have tested several values of the standard deviation, ranging from σ = 0.01 to 0.1, and of the temporal interval, from ∆t 1 = 1 yr to 100 yr. The smaller ∆t 1 , the more uniformly the phase space is filled, since the system cannot relax toward the attractor, and thus the reconstruction of the position of the attractors becomes unfeasible. On the contrary, too large values of ∆t 1 are computational expensive. On the other hand, the standard deviation σ regulates the range of forcing values being explored and the portion of phase space than can be accessed, a weaker noise corresponding to rarer transitions between the basins of attraction (Margazoglou et al., 2021).
In our MITgcm setup, the values that reproduce better the BD obtained with the standard method in a reasonable amount of time are σ = 0.025 and ∆t 1 = 10 yr. Such values are model-and setup-dependent. For example, in the Planet Simulator (PlaSim), an open-source intermediate complexity climate model with slab ocean (Fraedrich et al., 2005), the relaxation time is governed by upper ocean processes and is of the order of several decades for 5.6 • of horizontal resolution and 10 atmospheric levels (Margazoglou et al., 2021). In this case, the reconstruction of the invariant measure is obtained using σ = 0.18 and ∆t 1 = 1 yr (Margazoglou et al., 2021). We find that a good reconstruction occurs after a total simulation time of the order of t tot ∼ 10 4 yr. Using different seeds and/or initial conditions for the random variable, such total time can be easily partitioned over different set of runs.

Method II: Reconstruction of Stable Branches
The first step coincides with that in the standard technique: the attractors at a given forcing are needed as starting points and are obtained by scanning a large ensemble of initial conditions (Brunetti et al., 2019). However, the second step differs. We require the change in forcing to be small enough so that the invariant measure on the attractor remains nearly the same (Tél et al., 2020), in particular its mean and standard deviation. In other words, while the forcing is changed by ∆S, the corresponding variation on the attractor is small enough that the surface energy imbalance remains nearly zero. This excludes for the presence of R-tipping and guarantees quasi-ergodicity of the system. In practise, starting from each attractor found during the first step, the forcing is changed by +∆S (−∆S) each ∆t 2 = N 2 yr in order to determine the upper (lower) stable branch in the bifurcation diagram, until a tipping point is reached.
The tipping point is attained when one of the three following criteria is satisfied: (i) the standard deviation within the N 2 points for each forcing value becomes larger than the internal variability on the attractor (an early warning of critical slowing down of the dynamics (Scheffer et al., 2011)); (ii) the N 2 points turn out to be ordered in time, pointing from the original attractor toward a new one; (iii) the surface energy imbalance F s becomes much larger than zero (in our setup, larger than 0.5 W m −2 in absolute value, which is slightly larger than the threshold 0.2 W m −2 used to characterise the convergence to an attractor). When one of these criteria is satisfied, the system is no more on the initial dynamical attractor and is approaching the unstable branch where a shift toward a new basin of attraction takes place.
We have tested different values of ∆S from 0.1 to 0.5 W m −2 and N 2 from 10 to 100. Like in the previous case, the choice is model dependent. We have checked that in our setup good agreement with the bifurcation diagram from the standard technique (that is very accurate albeit highly time consuming) is obtained when ∆S = 0.25 W m −2 and N 2 ≥ 20. With less points, the criterium based on the standard deviation is not applicable (low statistics). With larger ∆S the position of B-tipping is less accurate and the requirement of remaining on the same original attractor cannot be satisfied, while too small values require large CPU time. The total simulation time can be estimated as t tot = nN t relax + M 2 ∆t 2 , where M 2 is the number of points in the forcing along all the branches. The second term is smaller than the analogous in the standard model since ∆t 2 t relax .

Results
We take advantage of the BD obtained with the standard (computational expensive) technique to test the other two methods described in Sections 2.2 and 2.3, and understand under which conditions they can be applied.

Method I
We construct the bifurcation diagram by plotting the normalised 2-dimensional histogram (invariant measure) in terms of the global annual surface air temperature and the forcing, as shown in Fig. 2 using a standard deviation of the normal distribution equal to σ = 0.025 and time interval ∆t 1 = 10 yr. Other BDs obtained for different values of σ and ∆t 1 can be found in Supporting Information (Fig. S1). As can be seen, the overall structure of the phase plane corresponds quite well to that in the standard BD (see Fig. 1 in Ragon et al. (2022)). The main attractors (hot, cold, waterbelt and snowball) can be easily recognised and correspond to regions of high density of points. However, the uncertainty in the exact position of the attractors is large, so that, for example, the warm state, with a short stable branch, cannot be distinguished from the hot climate and, in general, it is not possible to precisely infer the edges of the stable branches.
Moreover, a region with increased density different from the ones that correspond to the five attractors appears at T ∼ −20 • C that may correspond to an additional steady state. This can be checked by performing simulations without noise starting from initial conditions in such region of high density, and let the system relaxes towards a steady state under fixed forcing, i. e. until the surface energy imbalance becomes negligible. It turns out that the feature at T ∼ −20 • C is only transient, as shown in Supporting Information ( Figure S2): the system remains near such value of temperature for nearly 100 yr with a small surface imbalance but eventually is attracted to the waterbelt.

Method II
Starting from the five attractors found in Brunetti et al. (2019) at S 0 = 342 W/m 2 , the stable branches are derived by changing the forcing by ∆S each N 2 years. The resulting BD is shown in Fig. 3 using ∆S = 0.1 W m −2 and N 2 = 100 (20 for waterbelt upper branch). As can be seen, the five stable branches can be precisely recovered. If ∆S is too large, the edges of the stable branches are not correctly reproduced, as shown in Supporting Information ( Figure S3). Moreover, by applying the criteria listed in Section 2.3 to determine the position of B-tipping, we check that they correspond well to those found with the standard technique (see Figure 3a and its enlargements (panels b, c, d, e, f), where the red arrows correspond to B-tipping and colors correspond to the time ordering of the N 2 points). In Figure 4, two examples of evolution (hot to colder climates and cold to snowball) show when the system abandons the attractor at the point where the standard deviation of the temperature becomes larger than the internal variability on the attractor and/ or |F s | > 0.5 W m −2 (following the criteria (i) and (iii) in Section 2.3, respectively).
Another interesting feature in Figures 3a,c is the transition from the cold state to colder climates. When the system looses stability on the cold branch, it is attracted first by transient structures, like the one at T = −20 • C. Then, instead of relaxing to the waterbelt, it is directly attracted to the snowball. This can be understood by remembering that the climate state lives on a high dimensional space that we are arbitrarily projecting to a single state variable (the surface air temperature). In a higher dimensional space, the two climatic trajectories (the transition from cold to snowball and the waterbelt stable branch) would not cross each other. This behavior and the fact that an analogous crossing occurs in the transition from hot to snowball (see Figure 3a) show that the waterbelt climate is not well described by a single state variable.
Finally, the dependence of the invariant measure can be investigated, as done here, in terms of the incoming solar radiation S for a fixed atmospheric CO 2 content. The opposite can also be considered: different CO 2 for fixed S (Zhu & Rose, 2022). An exam- ple is shown in Supporting Information ( Figure S4) where, however, the MITgcm setup does not include the feedback of the atmospheric carbon with the ocean.

Summary and Conclusions
We have proposed two methods to construct BDs in general circulation models (GCM), where the number and timescale of nonlinear feedback mechanisms make the computational costs of the standard method prohibitive. Such diagrams store crucial information about the nonlinear structure of the climate system, like the regions of multistability, the kind of tipping, the amplitude of climatic oscillations or the intensity of forcing necessary to give rise to a climatic shift.
In the first method, random fluctuations of the forcing allow to explore the entire phase space of the climate system (Margazoglou et al., 2021). The choice of the parameters for such method, i. e. the standard deviation of the Gaussian fluctuations and the temporal interval after which the forcing randomly changes its value, is model dependent. However, once the parameters fixed, such method gives an overall picture of the phase space projected on a given state variable, revealing the number of attractors and their position. The disadvantage is that the attractors with small stable branches may be blurred and the position of B-tippings is not accurate.
The second method reconstructs the stable branches of the attractors by changing the forcing by a small amount that guarantees the system to remain on the same attractor until a B-tipping is reached. The signature of a shift is given by increased variance of the internal variability, non-null surface energy imbalance or an ordered temporal sequence of points toward another climatic state. The advantage of such a method is that the reconstruction of B-tipping is quite accurate with much lower computational costs with respect to the standard method.
The proposed methods have been described by projecting the high-dimensional space of the attractors to a single state variable (global annual mean of surface air temperature). Of course, they can be applied using other projections, and averages over different temporal and spatial scales. This is particularly useful in order to develop early warning indicators that are based not only on temporal series but also on spatial information (Nijp et al., 2019). From the BD one can infer the forcing that induces a tipping and perform a detailed analysis in space to identify where the main changes take place (gridded data in general circulation models and remote sensing permit such kind of spatial studies). Spatial early warning indicators do not need a long temporal record to get a meaningful signal, thus playing a crucial role to identify tipping mechanisms from datasets with irregular or infrequent temporal resolution (Génin et al., 2018).
Apart from the fundamental and practical interest of obtaining the BD for our presentday climate, there are many aspects that are worth to be analysed. An open question is, for example, whether the attractors and the BD are model dependent. The robustness of the results presented here should be tested against other climatic models. For the comparison, analogous numerical setups should be chosen at the level of horizontal/vertical resolutions and included nonlinear feedback mechanisms. Important information on the description of the nonlinear processes and their interplay can be gained from such comparison of BDs produced by different models, that can be used to improve algorithms and to correct biases. This is why we suggest including the comparison of BDs in coarseresolution GCMs with simplified configurations in the Tipping Point Model Intercomparison Project (TipMIP) beside that of Earth System Models in the present-day configuration.

Open Research
The data that support the findings of this study were generated by the MIT general circulation model that is openly available on GitHub (http://mitgcm.org/, https:// github.com/MITgcm/MITgcm, version c67f).