Articles | Volume 30, issue 3
Research article
21 Jul 2023
Research article |  | 21 Jul 2023

An approach for projecting the timing of abrupt winter Arctic sea ice loss

Camille Hankel and Eli Tziperman

Abrupt and irreversible winter Arctic sea ice loss may occur under anthropogenic warming due to the disappearance of a sea ice equilibrium at a threshold value of CO2, commonly referred to as a tipping point. Previous work has been unable to conclusively identify whether a tipping point in winter Arctic sea ice exists because fully coupled climate models are too computationally expensive to run to equilibrium for many CO2 values. Here, we explore the deviation of sea ice from its equilibrium state under realistic rates of CO2 increase to demonstrate for the first time how a few time-dependent CO2 experiments can be used to predict the existence and timing of sea ice tipping points without running the model to steady state. This study highlights the inefficacy of using a single experiment with slow-changing CO2 to discover changes in the sea ice steady state and provides a novel alternate method that can be developed for the identification of tipping points in realistic climate models.

1 Introduction

The Arctic is warming at a rate at least twice as fast as the global mean with profound consequences for its sea ice cover. Sea ice is already exhibiting rapid retreat with warming, especially in the summertime (Comiso and Parkinson2004; Nghiem et al.2007; Stroeve et al.2008; Notz and Stroeve2016; Stroeve and Notz2018), shortening the time that socioeconomic and ecological systems have to adapt. These concerns have motivated a large body of work dedicated to both observing present-day sea ice loss (Kwok and Untersteiner2011; Stroeve et al.2012; Lindsay and Schweiger2015; Lavergne et al.2019) and modeling sea ice to understand whether its projected loss is modulated by a threshold-like or “tipping point” behavior. Abrupt loss of Arctic sea ice could be driven by local positive feedback mechanisms (Curry et al.1995; Abbot and Tziperman2008; Abbot et al.2009; Kay et al.2012; Leibowicz et al.2012; Burt et al.2016; Feldl et al.2020; Hankel and Tziperman2021), remote feedback mechanisms that increase heat flux from the midlatitudes (Holland et al.2006; Park et al.2015), or by the natural threshold corresponding to the seawater freezing point (Bathiany et al.2016). If such an abrupt loss is caused by irreversible processes (typically, strong positive feedback mechanisms as opposed to the reversible mechanism of a freezing point threshold of Bathiany et al.2016), it is referred to here as a tipping point. A tipping point in the sense used here is a change in the number or stability of steady-state solutions (Ghil and Childress1987; Strogatz1994) as a function of CO2 and is also known as a bifurcation. We note that some of the climate literature uses the term tipping point in a more general sense of a relatively rapid change (e.g., Lenton2012). While most studies have concluded that there is no tipping point during the transition from perennial to seasonal ice cover (i.e., during the loss of summer sea ice), the existence of a tipping point during the loss of winter sea ice (transition to year-round ice-free conditions) continues to be debated in the literature (Eisenman2007, 2012; Eisenman and Wettlaufer2009; Notz2009). Wagner and Eisenman (2015) showed that a winter tipping point disappeared from a simple model of sea ice with no active atmosphere when a longitudinal dimension was added. On the other hand, other literature (e.g., Abbot and Tziperman2008; Hankel and Tziperman2021) has demonstrated the importance of atmospheric feedbacks, not included in the model of Wagner and Eisenman (2015), in inducing winter sea ice tipping point. Furthermore, three out of seven fully complex global climate models (GCMs) that lost their winter sea ice completely in the CMIP5 extended RCP8.5 scenario showed a very abrupt change in winter Arctic sea ice resembling a tipping point (Hezel et al.2014; Hankel and Tziperman2021). However, given the projected rapid changes to CO2 in the coming centuries and the slower response of the climate system, we do not expect future sea ice to be fully equilibrated to the CO2 forcing at a given time, making the standard steady-state tipping point analysis challenging. Thus, our first goal is to understand abrupt winter Arctic sea ice changes – which may or may not be due to tipping points – under rapidly changing CO2 forcing, where sea ice is not at equilibrium.

Tipping points imply a bistability (meaning that sea ice can take on different values for the same CO2 concentration), and hysteresis – an irreversible loss of sea ice even if CO2 is later reduced. Bistability (and therefore tipping points) can be tested for by running model simulations to steady state at many different CO2 values, which is computationally inefficient in expensive, state-of-the-art GCMs. GCM studies, therefore, tend to use a single experiment with very gradual CO2 increases and decreases (Li et al.2013) or even a faster CO2 change (Ridley et al.2012; Armour et al.2011) and look for hysteresis in sea ice that would imply the existence of a tipping point. These studies implicitly assume that such a run should approximate the behavior of the steady state at different CO2 concentrations. However, Li et al. (2013) further integrated two apparently bistable points and found that they equilibrated to the same value of winter sea ice: there was no “true” bistability at these two CO2 concentrations, the sea ice was simply out of equilibrium with the CO2 forcing. This calls into question the current use of time-changing CO2 runs to study the bifurcation structure of sea ice.

In light of the difficulties in using climate model runs with time-changing CO2 (hereafter “transient runs”), the first goal of this work is to understand the relationship between these transient runs and the steady-state value of sea ice in systems with and without bifurcations (since the existence of a bifurcation in winter sea ice remains unknown), and the second goal is to develop a new efficient method for the identification of tipping points from transient runs. Theoretical work in dynamical systems (Haberman1979; Mandel and Erneux1987; Baer et al.1989; Tredicce et al.2004) and studies related to bistability in the Atlantic Meridional Overturning Circulation (AMOC; Kim et al.2021; An et al.2021) have examined systems with tipping points when the forcing parameter (CO2 in our case) changes in time at a finite rate. They found that as the forcing parameter passes the bifurcation point, the system continues to follow the old equilibrium solution for some time before it rapidly transitions to the new one. Specifically, Kim et al. (2021) and An et al. (2021) find that the width of the hysteresis loop of AMOC is altered by the rate of forcing changes – this phenomenon is referred to as “rate-dependent hysteresis”. This rate dependence occurs in their case in a system that also has bistability and hysteresis in the equilibrium state. This type of analysis has, to our knowledge, not yet been applied in the context of winter sea ice loss under time-changing CO2 concentrations nor compared in systems with and without a bifurcation (that is, with and without an equilibrium hysteresis).

In order to analyze how the hysteresis of sea ice under time-changing forcing relates to the steady-state behavior of sea ice, we run a simple physics-based model of sea ice (Eisenman2007), configured in three different scenarios: with a large CO2 range of bistability, a small range of bistability, and no bistability in the equilibrium. These three scenarios span the range of possible behaviors of winter sea ice in state-of-the-art climate models. Each case is run with different rates of CO2 increase (ramping rates). We use results from this model and from an even simpler standard 1D dynamical system to demonstrate that the convergence of the transient behavior (under time-changing forcing) to the equilibrium behavior is very slow as a function of the ramping rate of CO2. In other words, even climate model runs with very slow-changing CO2 forcing may simulate sea ice that is considerably out of equilibrium near the period of abrupt sea ice loss. Finally, we propose a novel approach for uncovering the underlying equilibrium behavior – and thus the existence and location of tipping points – in comprehensive models where it is computationally infeasible to simulate steady-state conditions for many CO2 values. Such a method is important given the model-dependent nature of winter sea ice tipping points discussed above; uncovering the existence of sea ice tipping points in GCMs, which are the most realistic representation of Arctic-wide sea ice behavior that we have, is the next step toward understanding whether such tipping points exist in the real climate system. Our goal has some parallels to that of Gregory et al. (2004), who used un-equilibrated GCM runs to deduce the equilibrium climate sensitivity when fully equilibrated runs were computationally infeasible.

As mentioned above, some GCMs exhibit an abrupt change in winter sea ice that may be a tipping point, and others do not (Hezel et al.2014; Hankel and Tziperman2021). The reasons likely involve numerous differences in parameters and parameterizations. It is not obvious how to modify parameters in a single GCM to display all of these different behaviors. Therefore, we choose to use an idealized model of sea ice where we can directly produce different bifurcation behaviors to address our second goal and answer the question: is it possible to identify the CO2 at which tipping points occur without running the model to a steady state for many CO2 values? Answering such a question in a simple model is an obvious prerequisite to tackling the problem of identifying climate bistability in noisy, high-dimensional GCMs. In order to perform this analysis for each of the three scenarios mentioned above, we modify the strength of the albedo feedback via the choice of surface albedo parameters. The albedo values used here to generate the three scenarios are not meant to reflect realistic albedo values but rather allow us to represent in a single model the range of sea ice equilibria behaviors that may exist in different GCMs. We, therefore, follow in the footsteps of previous studies (e.g., Eisenman2007) that have also changed parameters outside of their physically relevant regime in order to understand summer sea ice bifurcation behavior; here we follow the same approach to understand when a winter sea ice bifurcation can be detected without running an expensive climate model to steady state.

2 Methods

2.1 Sea ice model

The sea ice model used follows Eisenman (2007) almost exactly, and its key features are depicted schematically in Fig. 1. The model contains four state variables: sea ice effective thickness (V, which is volume divided by the area of the model grid box), sea ice area (A), sea ice surface temperature (Ti), and mixed-layer temperature (Tml) for a single box representing the entire Arctic. Subsequent versions of this sea ice model have been used in Eisenman and Wettlaufer (2009), Eisenman (2012), and Wagner and Eisenman (2015). Those versions are derived from the model used here, making a few further modest simplifications (using a hyperbolic tangent function for surface albedo, assuming the ice surface temperature is in a steady state, combining all prognostic variables into one, i.e., enthalpy) that do not affect the qualitative behavior of the model (i.e., the nature of summer and winter sea ice bifurcations). We choose to implement the earlier model because it explicitly represents the key physical variables of ice volume, area, ocean temperature, and ice temperature as prognostic variables – as opposed to combining them all into a single enthalpy – and thus provides more transparency and interpretability. We, therefore, do not expect our results to change if we use any of the later model versions.

Figure 1Schematic showing some of the key features of the Eisenman (2007) model. Its four prognostic variables are ice volume, ice area, ice surface temperature, and ocean mixed-layer temperature. The full model equations can be found in the Supplement.


In the model, the atmosphere is assumed to be in radiative equilibrium with the surface, and the model is forced with a seasonal cycle of insolation; of poleward atmospheric heat transport from the midlatitudes; and of local optical thickness of the atmosphere, which represents cloudiness. Sea ice growth and loss are primarily determined by the heat budget at the bottom of the ice and are therefore set by the balance between ocean–ice heat exchanges and heat loss through the ice to the atmosphere. When conditions for surface melting are met (when the ice surface temperature is zero and net fluxes on the ice are positive), all surface heating goes into melting ice and the surface albedo of the ice is set to the melt pond albedo. The ocean temperature is affected by shortwave and longwave fluxes in the fraction of the box that is ice-free and by ice–ocean heat exchanges. When the ocean temperature reaches zero, all additional cooling goes into ice production, while the ocean temperature remains constant. The full equations of the sea ice model can be found in the original paper (Eisenman2007) and in the Supplement; here, we highlight a few minor ways in which our implementation differs. First, for simplicity, we do not model leads, which in the original model were represented by capping the ice fraction at 0.95 rather than 1. Second, we use an approximation to the seasonal cycle of insolation (Hartmann2015) using a latitude of 75 N. The atmospheric albedo is set to 0.425 to produce the same magnitude of the seasonal cycle as in the original model of Eisenman (2007).

2.2 Setup of simulations

In our transient-forcing scenarios (described below), we vary CO2 in time which affects the prescribed near-surface atmospheric midlatitude temperature (Tmidlat) and the atmospheric optical depth (N; see the Supplement). Specifically, we increase the annual mean of Tmidlat by 3 C per CO2 doubling and N by a ΔN that corresponds to 3.7 W m−2 per doubling. All model parameters are as in Eisenman (2007) except as mentioned below.

We configure the model in three different scenarios that yield a wide CO2 range of bistability in winter sea ice (Scenario 1), a small range of bistability in winter sea ice (Scenario 2), and no bistability in winter sea ice (Scenario 3). We do so by modifying the strength of the ice–albedo feedback by changing the albedos of bare ice (αi), melt ponds (αmp), and ocean (αo), as listed in Table S1 in the Supplement.

Figure 2Hysteresis runs (time-changing forcing) and equilibrium runs (fixed forcing) for average March sea ice effective thickness (sea ice volume divided by area of the grid cell; panels a, c, and d) and the simple ODE from Eq. 1 (b, d, and f). The first row corresponds to Scenario 1 (wide bistability), the second row to Scenario 2 (narrow bistability), and the third to Scenario 3 (no bistability). Blue lines indicate simulations with increasing forcing (CO2 or β), while red lines indicate simulations with decreasing forcing. Dashed and dotted black lines indicate the steady-state values of sea ice or the ODE variable x. These two black lines are different when the two initial conditions evolve to two different steady states. The legends indicate the different ramping rates (represented by darker colors for faster rates), which are in units of years per CO2 doubling in the case of the sea ice model. The green arrows demonstrate the direction of evolving sea ice effective thickness during the hysteresis experiments.


In each of the three scenarios, we tune the model (by adjusting the mean and amplitude of the atmospheric optical depth) to roughly match the observed seasonal cycle of ice thickness under pre-industrial CO2 ( 2.5–3.7 mEisenman2007). We then run each scenario with multiple CO2 ramping rates (expressed in “years per doubling”) with an initial stabilization period (fixed pre-industrial CO2), a period of exponentially increasing CO2 concentration (which corresponds to linearly increasing radiative forcing), another period of stabilization at the maximum CO2, a period of decreasing CO2, and a final period of stabilization at the minimum CO2 value (see Fig. S2 in the Supplement). Scenarios 2 and 3 are ramped to higher final CO2 values than Scenario 1 so that they lose all their sea ice. We also directly calculate the steady-state behavior of the sea ice (as done in the original study) by running many simulations with fixed CO2 values until the seasonal cycle of all the variables stabilizes. Because we expect multiple equilibria (which could be ice-free, seasonal ice, or perennial ice) at some CO2 values in Scenarios 1 and 2, we run these steady-state simulations starting with both a cold (ice-covered) and a warm (ice-free) initial condition in order to find these different steady states. In the ice-free initial condition runs, the ice–albedo feedback will still play an important role if the temperature cools sufficiently for ice to develop. At CO2 values for which the sea ice is bistable, the ice-free initial condition evolves to a perennially ice-free steady state, and the ice-covered initial condition evolves to a seasonally ice-covered steady state (seen by the dotted and dashed lines, respectively, in Fig. 2a and c).

2.3 Cubic ODE

The main points we are trying to make about the transient versus equilibrium behavior of winter sea ice near a tipping point are not unique to the problem of winter sea ice, and in order to demonstrate this, we use the simplest mathematical model that can display tipping points, following other studies that have also used such simple dynamical systems (Ditlevsen and Johnsen2010; Bathiany et al.2018; Ritchie et al.2021; Boers2021). The cubic ODE used, while much simpler than the sea ice model above, has some of the key characteristics of the sea ice system (it is a non-autonomous system due to the time-dependent forcing and has saddle–node bifurcations), which allows for direct comparison between the two models. The ODE equation,

(1) d x d t = - x 3 + δ x + β ( t ) , β ( t ) = β 0 + μ t ,

contains a time-changing forcing parameter, β(t) mimicking the effects of CO2 in the sea ice model. We consider this differential equation in three scenarios, paralleling those used with the sea ice model: in Scenario 1, δ=5 leading to a wide region of bistability; in Scenario 2, δ=1 leading to a narrow region of bistability; and finally, in Scenario 3, δ=0 leading to a mono-stable system. The different values of δ, therefore, produce the same three scenarios that were achieved in the sea ice model by modifying the strength of the ice–albedo feedback. We mimic the hysteresis experiments of the sea ice model with a sequence of ramping up and ramping down (using different ramping rates, μ) with values of β ranging from 10 to 10 to sweep the parameter space that contains the bifurcations. We calculate the steady states with fixed values of β (μ=0), starting with both a positive and a negative initial condition of x to yield two stable solutions when these exist.

We want to calculate the upper and lower CO2 values of the hysteresis region in runs with time-changing (i.e., transient) CO2 forcing. We do so by calculating the CO2 value at which the March sea ice area drops below a critical threshold (50 % ice coverage; results are insensitive to the specific value used) during increasing and decreasing CO2 integrations: we denote these CO2 values CO2i and CO2d, respectively (see Fig. S9). The difference between CO2i and CO2d is referred to below as the “hysteresis width” of the rate-dependent hysteresis whether an equilibrium hysteresis exists or not; this width approaches the width of bistability at very slow ramping rates.

2.4 A new method for predicting the CO2 of the sea ice tipping point

One of our main goals (see the Introduction) is to efficiently estimate the equilibrium behavior of sea ice, including the location of tipping points, without running the model to a steady state for many CO2 values. This would show that such estimation could be calculated for GCMs where tipping points cannot be detected using steady-state runs due to their computational cost. In order to estimate the values of CO2i and CO2d that would have occurred for an infinitely slow ramping rate (in other words, the range of CO2 for which there is bistability) using only the transient runs, we fit a polynomial of the form f(x)=mxc+b to CO2i and CO2d as functions of the ramping rate x. Because c is negative, the fitted parameter b represents the prediction of CO2i and CO2d at infinitely slow ramping rates, i.e., in the steady state. We also calculate the uncertainty in the fitted parameter b by block bootstrapping to account for autocorrelation; see the Supplement. Other fits to CO2i and CO2d as a function of ramping rates, such as an exponential function f(x)=a+bexp(-cx), could in principle be used, although we found that fit to be less good in our case.

3 Results

In the following three subsections, we discuss the behavior of the sea ice model and the cubic ODE under time-changing forcing, the relationship of the transient and equilibrium behaviors, and a method that we propose for inferring the existence and location of tipping points from the transient behavior. Equilibrium hysteresis refers here to the path-dependent solution of a variable due to bistability and a bifurcation in the steady state (in other words, the loop traced by the steady-state solutions). The term rate-dependent hysteresis (An et al.2021; Manoli et al.2020) describes hysteresis loops that appear in time-changing forcing runs (rather than in the steady state) and that depend on the rate of forcing change. In our analysis rate-dependent hysteresis applies to both systems with and without equilibrium hysteresis: it refers to any differences in the results for increasing vs. decreasing CO2 simulations of sea ice that are altered by the rate of CO2 change.

3.1 Transient response of Arctic winter sea ice to time-changing CO2

Our goal in this section is to understand the relationship of winter sea ice forced with time-changing CO2 to its equilibrium state, both in cases with and without a sea ice tipping point. In Fig. 2a, c, and e, we plot the results of running all three scenarios (wide range of bistability, Scenario 1; narrow range of bistability, Scenario 2; and no bistability, Scenario 3) under time-changing (transient) and fixed CO2 values. In all scenarios, the experiments run with time-changing CO2 exhibit rate-dependent hysteresis; the hysteresis width (lower horizontal gray bar in Fig. 2a) is larger for faster ramping rates (Fig. 2a, c, and e). For Scenarios 1 and 2, which have a region of bistability and equilibrium hysteresis (upper gray bar in Fig. 2a), this corresponds to a widening from the equilibrium hysteresis (that would exist even with infinitely slow ramping rates), while in Scenario 3, this hysteresis occurs only in transient simulations and is due to the inertia in the system (the sea ice cannot respond instantaneously to forcing changes). In Scenarios 1 and 2, whose equilibrium solutions (dashed and dotted black lines in Fig. 2) have a tipping point and therefore an infinite gradient of sea ice thickness vs. CO2, the faster ramping rates also lead to more gradual (and finite) gradient of sea ice thickness vs. CO2.

The rate-dependent hysteresis loops across all scenarios at fast enough ramping rates (loops composed of the darkest blue and darkest red in Fig. 2a, c, and e) are qualitatively similar in shape, despite their different underlying steady-state structures. This similarity indicates that from a single hysteresis run with time-changing CO2, we cannot discern whether the underlying Arctic winter sea ice equilibrium behavior has a region of bistability or not nor how wide the region of true bistability is. In particular, a single hysteresis loop found from a time-changing forcing simulation would always overestimate the width of bistability if it was assumed to represent a quasi-steady state. This result demonstrates that the apparent sea ice hysteresis loop found by Li et al. (2013) could be due to a system without an equilibrium hysteresis, as they suggest, or due to a system with a narrower equilibrium hysteresis than the one implied by their transient simulation.

We now discuss the behavior of the simple cubic ODE (Eq. 1) under similarly time-changing forcing. Previous work in the dynamical systems literature (e.g., Haberman1979; Mandel and Erneux1987; Baer et al.1989; Breban et al.2003; Tredicce et al.2004; Kaszás et al.2019) has examined a variety of simple systems to understand the nature of bifurcations in the presence of a time-changing (“drifting” or “transient”) forcing parameter. In the climate literature, too (e.g., Ditlevsen and Johnsen2010; Bathiany et al.2018; Ritchie et al.2021; Boers2021), idealized dynamical systems similar to our Eq. (1) have been used to understand the predictability of tipping points in the presence of noise, and the ability to recover from such tipping points (“overshoot” scenarios). These works, as well as the AMOC study of An et al. (2021), found that a system with a bifurcation that is run with a time-changing forcing parameter can follow a given equilibrium value beyond the bifurcation value of the forcing parameter before undergoing the tipping point transition to the new equilibrium value. This is consistent with the out-of-equilibrium behaviors we find for sea ice in Scenarios 1 and 2. To our knowledge, the simple ODE used here has not yet been analyzed with our specific goal in mind: to compare the shape of rate-dependent hysteresis loops in generic dynamical systems both with and without bifurcations and to address the question of whether the equilibrium behavior can be inferred from the rate-dependent behavior of such systems.

To address these two goals, we configure Eq. (1) analogously to the sea ice model in three scenarios with wide bistability (Scenario 1), narrow bistability (Scenario 2), and no bistability (Scenario 3) and force it with a time-changing forcing parameter. In Fig. 2b, d, and f, we see that the three scenarios with similar dynamics (but different equilibrium structures) all display rate-dependent hysteresis, similar to the result from the sea ice model. Specifically, even when there is only one stable equilibrium solution in both models (Scenario 3, panels e and f of Fig. 2), there is still a narrow region of rate-dependent hysteresis. Thus, we find that the inability to tell if rate-dependent hysteresis in Arctic winter sea ice is accompanied by an underlying equilibrium hysteresis appears to be a generic feature of dynamical systems, which helps explain the challenges of interpreting the results of Li et al. (2013).

Mathematically, this 1D system is fundamentally different from the sea ice model because it is not periodically forced. We show in the Supplement that adding a sinusoidal forcing term to the ODE does not qualitatively change our results.

3.2 Slow convergence of the rate-dependent hysteresis to the equilibrium behavior

Our next objective is to demonstrate that it would require expensive runs in a GCM to approach the equilibrium behavior of sea ice using slower and slower-changing CO2 runs (hysteresis experiments). As we saw in Fig. 2, the rate of loss of sea ice with increasing CO2 is infinite (dashed and dotted black lines) in Scenarios 1 and 2 at the tipping points. On the other hand, the gradient of sea ice thickness with respect to CO2 is more gradual and finite under time-changing forcing (blue and red curves) but steepens as the ramping rate of CO2 decreases. We now quantify the rate of this steepening by examining the maximum gradient of sea ice loss during each transient simulation as a function of ramping rate (inverse of the years per doubling of CO2).

Figure 3Maximum gradient of sea ice effective thickness with respect to CO2 in (a) and the maximum gradient of x with respect to the forcing parameter β in (b) during transient simulations. For the sea ice model (a) the data points from the 18 different runs are shown as faded points, with a superimposed line of best fit. For the cubic ODE (b) the maximum gradient lines corresponding to increasing and decreasing forcing time series are identical due to the symmetry around β=0 seen in Fig. 1b, d, and f.


In Fig. 3a, we plot the maximum gradient of March sea ice thickness with respect to CO2 during each hysteresis experiment, as a function of the CO2 ramping rate. In Scenarios 1 and 2 (wide and narrow bistability, respectively), the maximum gradient gets greater as the ramping rate is slower (Fig. 3a, negative slopes of solid and dashed lines), consistent with Fig. 2 (e.g., steepening from dark-blue to light-blue curves in Fig. 2a and b). In particular, the gradient approximately follows a negative power law as a function of ramping rate on both warming and cooling time series. In Scenario 3, the maximum gradient is nearly insensitive to the ramping rate (relatively flat dash-dotted lines seen in Fig. 3a). In Fig. 3b, we see a similar result for the simple ODE, as seen by the shallowing of the power law from Scenarios 1 to 3 (though here the slope in Scenario 3 is clearly nonzero). Notably, in the cubic ODE the power law in the case with the largest region of bistability (Scenario 1) is approximately given by max(dx/dβ)μ-1, where μ again is the ramping rate. The Supplement further explains the above convergence rate of µ−1.

A dependence on the reciprocal of the ramping rate in the case of wide bistability suggests that running a climate model with twice as gradual CO2 ramping leads to a less than a factor of 2 increase in the gradient max(dV/dCO2). This is an important result because this implies that the distance between the CO2 at the simulated transient tipping point and the CO2 of the true (equilibrium) tipping point (which we want to estimate) also only reduces by a factor of 2 when the ramping rate is reduced by a factor of 2. A greater power law slope (e.g., a slope of −2) would imply a much faster convergence to the equilibrium location of the tipping point. Thus, using more and more gradual ramping experiments may be an inefficient way to approach the equilibrium behavior of this physical system, suggesting the need for a more efficient approach, discussed next.

3.3 Predicting the steady-state behavior of sea ice using only transient runs

Our main novel result, presented next, is a method for finding the CO2 concentration at which a bifurcation (if any) occurs in the equilibrium using computationally feasible transient model runs instead of fixed-forcing steady-state runs. We are interested in this CO2 concentration because it determines the threshold beyond which significant sea ice loss is practically irreversible (Ritchie et al.2021). In our simple, inexpensive model, we can test the estimates of the bistability and associated tipping points derived from transient model runs against the known true tipping points and equilibrium structure that are found from fixed-forcing runs (see Methods). When used in a GCM, our method would provide a prediction for the existence and location of tipping points when the equilibrium value of sea ice is actually unknown. Thus, this section is a proof of concept that our new method can accurately determine whether observed rate-dependent hysteresis is caused by lag around a system with no bistability or tipping points or is caused by a rate-dependent widening of an equilibrium hysteresis loop in a system with tipping points.

Figure 4Estimating the equilibrium tipping point value from the rate-dependent hysteresis runs. In (a), the scatter points show the CO2 value of the right and left edges of the rate-dependent hysteresis (CO2i and CO2d, located along increasing (blue) and decreasing (red) CO2 time series, respectively) for different ramping rates. The dashed lines show the curve that is fitted to the scatter points, and the shaded blue and red bands show ± 2σ around the predicted values of CO2i and CO2d at infinitely slow ramping rates. The blue and red x symbols show the true equilibrium values of CO2i and CO2d (calculated from the fixed-CO2 runs starting with cold and warm initial conditions, respectively). In (b), we analyze the accuracy of this prediction as we use fewer transient runs. For the three scenarios, we show the result of sequentially excluding the most gradual ramping simulations from the curve-fitting process used for predictions. The dots and the corresponding bars represent the predicted equilibrium values of CO2i and CO2d, with ±2σ around the prediction, and dots moving away from the true value with larger error bars correspond to excluding more and more runs from the calculation.


In Fig. 4a, we plot a measure of the upper and lower CO2 values that correspond to the rightmost and leftmost edges of the rate-dependent hysteresis (by calculating the CO2 at which the March sea ice area crosses a critical threshold; see Methods and Fig. S9). We plot this measure for the warming (increasing greenhouse concentration) trajectories in blue (CO2i) and for the cooling (decreasing greenhouse) trajectories in red (CO2d), as a function of the ramping rate for all three scenarios. As expected, as the ramping rate gets slower CO2i and CO2d approach the CO2 values corresponding to the edges of the equilibrium hysteresis and the location of the true tipping points in the case of Scenarios 1 and 2 (denoted by the × symbols). In Scenario 3, CO2i and CO2d approach the same value (the rate-dependent hysteresis width approaches zero) because there is no bistability in the steady state.

Finally, we demonstrate that fitting a curve to the edges of the rate-dependent hysteresis (CO2i and CO2d) as a function of the ramping rate can be used to predict CO2i and CO2d at infinitely slow ramping rates (i.e., the edges of the equilibrium hysteresis). This would allow us to estimate the CO2 value corresponding to a bifurcation in the equilibrium behavior without running a model to a steady state. In Fig. 4a, we plot CO2i and CO2d and the curves that fit them (see Methods) as functions of the ramping rate and the predicted values of CO2i and CO2d at infinitely slow ramping rates with a 95 % confidence interval range shaded around them. We perform this fitting and estimation process using all the ramping experiments (18 different ramping rates total, as shown in Fig. 4a). We then repeat the fit using fewer and fewer experiments to explore how the uncertainty in predicted values of CO2i and CO2d increases as we move to only using a few fast ramping experiments that are more feasible when using full-complexity climate models. Figure 4b shows a summary of these analyses.

The predicted values of CO2i and CO2d are remarkably accurate for all scenarios (points approaching the red and blue x symbols in Fig. 4b), even when excluding several of the slower ramping experiments. This is an important test because when this method is applied to a GCM, one would only have a smaller number of faster-ramping experiments due to computational limitations. The uncertainties (indicated by the shaded blue and red bars around the points) in the predictions grow when excluding more experiments from the curve fitting process but still remain very low, especially for Scenarios 1 and 2. In predicting CO2d for Scenario 3, the uncertainties are a bit higher because the functional form of our fit does not represent this case as well as the others, leading to serial correlation in the residuals. The structure in the residuals can be used to guide the choice of the functional form used to fit such model output in future applications. This same method and functional form can also successfully predict the equilibrium structure of our simple ODE (Eq. 1), with even smaller uncertainties in the prediction when using very few ramping experiments (see Fig. S11). Finally, we can use the difference between the distributions CO2i and CO2d to calculate the probability that bistability – and thus a tipping point – exists (see the Supplement). Another very similar approach using only the difference between CO2i and CO2d (i.e., the hysteresis width) as a function of the ramping rate is also shown in Fig. S10.

Overall, these results demonstrate the potential for using several shorter runs with time-changing CO2 forcing to efficiently estimate the CO2 value of the tipping points and predict the existence of bistability in GCMs where equilibrium runs or long, slow-ramping hysteresis runs are computationally infeasible.

4 Discussion

We have shown that a single climate model hysteresis run with time-changing (transient) forcing cannot be used to conclusively estimate the true location of Arctic winter sea ice tipping points, the range of bistability in the steady state, and even the existence of bistability at all. We demonstrated that the transient sea ice responses under time-changing CO2 reflect the generic behavior of a nonlinear dynamical system (e.g., our Eq. 1): specifically, we showed that systems with and without bistability can also produce qualitatively indistinguishable rate-dependent hysteresis behavior. We also find that very long model runs are needed to identify whether the system approaches a bifurcation (Fig. 3) and at what CO2 this occurs. We showed that even in runs with a very slow-changing CO2, the system can be surprisingly far from the equilibrium as it undergoes a tipping point, consistent with the work of Li et al. (2013). In addition, even with a very slow-ramping experiment, one would always have to perform additional expensive fixed-forcing experiments (as done by Li et al.2013) to confirm that the experiment was indeed in quasi-equilibrium. Instead, we propose a novel method that uses a few fast-ramping experiments to efficiently predict the true range of bistability and provide uncertainty estimates on this prediction.

We demonstrated that the method we propose can accurately predict the steady-state behavior of sea ice in a simple model; now we discuss applying this method to a GCM. First, we note that while we use a highly idealized model of sea ice in this study, the method developed deals with identifying bistability in complex systems with unknown equilibrium structures more generally. This means that the framework should be applicable to other models (including GCMs), since moving from fast- to slower-ramping rates allows convergence to the equilibrium behavior. It could also be used in the context of vastly different climate problems, for example, in identifying the abrupt transitions to a moist greenhouse (Popp et al.2016), runaway greenhouse (Goldblatt et al.2013), or snowball Earth state (Hyde et al.2000). The functional form used to fit the transient runs, as well as the level of certainty achieved from a given number of experiments, would likely depend on the given model and climate problem analyzed. Possible challenges in finding the functional best fit to the transient runs might mirror those of Gregory et al. (2004), who encountered difficulties when trying to fit a line to un-equilibrated GCM runs with a different goal of deducing the equilibrium climate sensitivity. We suggest that a careful examination of the residuals from a given fit can help guide the choice of functional form.

The generality of the method also highlights another advantage: the same set of ramping experiments in a GCM could be used to analyze all suspected tipping elements in the Earth's climate system simultaneously. The main challenge we anticipate in applying this method to GCMs comes from the significant stochastic variability and multiple timescales of forcings that may render the calculated width of the rate-dependent hysteresis more uncertain in a GCM. Nonetheless, using multiple runs to estimate the width of the bistability of a given climate variable and providing a quantified uncertainty in such a prediction should offer a potential improvement over using a single hysteresis experiment.

We can estimate the efficiency of the proposed approach over more standard ones when applied in a GCM. Taking the experimental setup of Li et al. (2013) as a guide, we can assume that a slow-ramping experiment to four times CO2 requires a 2000-year ramp-up and ramp-down with at minimum a 2500-year equilibration period after each ramp (though they actually allowed the model to equilibrate for nearly 6000 years). Within the 500 ppm width of the rate-dependent hysteresis found by Li et al. (2013), 10 fixed-forcing experiments 2500 years long would be needed to test for bistability and estimate the tipping point location at a relatively crude accuracy of 100 ppm. This leads to a total of 34 000 simulation years. On the other hand, if we used our proposed approach, we could run three ramping experiments with fast to intermediate rates of 100, 200, and 400 years to quadruple CO2. We would run only one experiment to complete equilibration after ramp-up (2500 years) and run the others only until they lost their sea ice, using the ice-free steady-state run to conduct the three ramp-downs. This yields a total of approximately 6400 simulation years and computational savings of over a factor of 5. Using only three ramping experiments is sufficient to get an estimate of the equilibrium hysteresis width and location, but the uncertainty in the estimate could still be high.

Finally, our results indicate that rate-dependent hysteresis and irreversibility of Arctic winter sea ice are expected to be relevant for realistic rates of CO2 increase. While rate-dependent hysteresis has been explored in other climate contexts (e.g., AMOC; Kim et al.2021; An et al.2021), previous work on Arctic winter sea ice has typically sought to identify equilibrium hysteresis in sea ice because it would imply irreversibility of sea ice loss, generally ignoring the out-of-equilibrium behavior of sea ice under rapid CO2 changes. The SSP585 scenario in CMIP6 corresponds to a ramping rate of approximately 60 years per CO2 doubling; this is a rate at which sea ice in our idealized model already exhibits rate-dependent hysteresis and significant deviation from its steady state (see Figs. 2 and S2). Since we identify rate-dependent hysteresis in sea ice here in all scenarios, even without a deep ocean and subsequent recalcitrant warming (Held et al.2010), we expect rate-dependent hysteresis to be even more pronounced in GCMs and in the real climate when such long-timescale components are included. We, therefore, conclude that on policy-relevant timescales the significant irreversibility of winter Arctic sea ice involved in rate-dependent hysteresis is likely to occur in the real climate system due to the expected lagged response regardless of whether an actual bifurcation (tipping point) in the equilibrium exists.

Code availability

An implementation of the Eisenman (2007) sea ice model in Python used for this study can be found on Zenodo at (Hankel2022).

Data availability

No data sets were used in this article.


The supplement related to this article is available online at:

Author contributions

CH and ET designed the research project and prepared the paper together; CH implemented the model and conducted the experiments.

Competing interests

The contact author has declared that neither of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to thank Ian Eisenman for his helpful input during the project and for the guidance in using his sea ice model. We thank the anonymous reviewers for their constructive feedback. Eli Tziperman thanks the Weizmann Institute for its hospitality during parts of this work.

Financial support

This research has been supported by the National Science Foundation (grant no. AGS-1924538).

Review statement

This paper was edited by Stefano Pierini and reviewed by three anonymous referees.


Abbot, D. S. and Tziperman, E.: Sea ice, high-latitude convection, and equable climates, Geophys. Res. Lett., 35, L03702,, 2008. a, b

Abbot, D. S., Walker, C., and Tziperman, E.: Can a convective cloud feedback help to eliminate winter sea ice at high CO2 concentrations?, J. Climate, 22, 5719–5731,, 2009. a

An, S.-I., Kim, H.-J., and Kim, S.-K.: Rate-Dependent Hysteresis of the Atlantic Meridional Overturning Circulation System and Its Asymmetric Loop, Geophys. Res. Lett., 48, e2020GL090132,, 2021. a, b, c, d, e

Armour, K., Eisenman, I., Blanchard-Wrigglesworth, E., McCusker, K., and Bitz, C.: The reversibility of sea ice loss in a state-of-the-art climate model, Geophys. Res. Lett., 38, L16705,, 2011. a

Baer, S. M., Erneux, T., and Rinzel, J.: The slow passage through a Hopf bifurcation: delay, memory effects, and resonance, SIAM J. Appl. Math., 49, 55–71, 1989. a, b

Bathiany, S., Notz, D., Mauritsen, T., Raedel, G., and Brovkin, V.: On the potential for abrupt Arctic winter sea ice loss, J. Climate, 29, 2703–2719, 2016. a, b

Bathiany, S., Scheffer, M., Van Nes, E., Williamson, M., and Lenton, T.: Abrupt climate change in an oscillating world, Sci. Rep.-UK, 8, 1–12, 2018. a, b

Boers, N.: Observation-based early-warning signals for a collapse of the Atlantic Meridional Overturning Circulation, Nat. Clim. Change, 11, 680–688, 2021. a, b

Breban, R., Nusse, H. E., and Ott, E.: Lack of predictability in dynamical systems with drift: scaling of indeterminate saddle-node bifurcations, Phys. Lett. A, 319, 79–84, 2003. a

Burt, M. A., Randall, D. A., and Branson, M. D.: Dark warming, J. Climate, 29, 705–719, 2016. a

Comiso, J. C. and Parkinson, C. L.: Satellite observed changes in the Arctic, Phys. Today, 8, 38–44,, 2004. a

Curry, J. A., Schramm, J. L., and Ebert, E. E.: Sea ice–albedo climate feedback mechanism, J. Climate, 8, 240–247, 1995. a

Ditlevsen, P. D. and Johnsen, S. J.: Tipping points: Early warning and wishful thinking, Geophys. Res. Lett., 37, L19703,, 2010. a, b

Eisenman, I.: Arctic catastrophes in an idealized sea ice model, 2006 Program of Studies: Ice (Geophysical Fluid Dynamics Program), 133–161, 2007. a, b, c, d, e, f, g, h, i

Eisenman, I.: Factors controlling the bifurcation structure of sea ice retreat, J. Geophys. Res.-Atmos., 117, D01111,, 2012. a, b

Eisenman, I. and Wettlaufer, J. S.: Nonlinear threshold behavior during the loss of Arctic sea ice, P. Natl. Acad. Sci. USA, 106, 28–32, 2009. a, b

Feldl, N., Po-Chedley, S., Singh, H. K., Hay, S., and Kushner, P. J.: Sea ice and atmospheric circulation shape the high-latitude lapse rate feedback, npj Climate and Atmospheric Science, 3, 1–9, 2020. a

Ghil, M. and Childress, S.: Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory and Climate Dynamics, Springer-Verlag, New York, ISBN 978-0-387-96475-1, 1987. a

Goldblatt, C., Robinson, T. D., Zahnle, K. J., and Crisp, D.: Low simulated radiation limit for runaway greenhouse climates, Nat. Geosci., 6, 661–667, 2013. a

Gregory, J., Ingram, W., Palmer, M., Jones, G., Stott, P., Thorpe, R., Lowe, J., Johns, T., and Williams, K.: A new method for diagnosing radiative forcing and climate sensitivity, Geophys. Res. Lett., 31, L03205,, 2004. a, b

Haberman, R.: Slowly varying jump and transition phenomena associated with algebraic bifurcation problems, SIAM J. Appl. Math., 37, 69–106, 1979. a, b

Hankel, C.: camillehankel/sea_ice_thermo_0d: 0D-Sea-Ice-Model, Zenodo [code],, 2022. a

Hankel, C. and Tziperman, E.: The Role of Atmospheric Feedbacks in Abrupt Winter Arctic Sea Ice Loss in Future Warming Scenarios, J. Climate, 34, 4435–4447, 2021. a, b, c, d

Hartmann, D. L.: Global physical climatology, vol. 103, Chapter 2, Newnes, Elsevier, ISBN 978-0-12-328531-7, 2015. a

Held, I. M., Winton, M., Takahashi, K., Delworth, T., Zeng, F., and Vallis, G. K.: Probing the fast and slow components of global warming by returning abruptly to preindustrial forcing, J. Climate, 23, 2418–2427, 2010. a

Hezel, P. J., Fichefet, T., and Massonnet, F.: Modeled Arctic sea ice evolution through 2300 in CMIP5 extended RCPs, The Cryosphere, 8, 1195–1204,, 2014. a, b

Holland, M. M., Bitz, C. M., and Tremblay, B.: Future abrupt reductions in the summer Arctic sea ice, Geophys. Res. Lett., 33, L23503,, 2006. a

Hyde, W. T., Crowley, T. J., Baum, S. K., and Peltier, W. R.: Neoproterozoic 'snowball Earth' simulations with a coupled climate/ice-sheet model, Nature, 405, 425–429, 2000. a

Kaszás, B., Feudel, U., and Tél, T.: Tipping phenomena in typical dynamical systems subjected to parameter drift, Sci. Rep.-UK, 9, 8654,, 2019. a

Kay, J. E., Holland, M. M., Bitz, C. M., Blanchard-Wrigglesworth, E., Gettelman, A., Conley, A., and Bailey, D.: The influence of local feedbacks and northward heat transport on the equilibrium Arctic climate response to increased greenhouse gas forcing, J. Climate, 25, 5433–5450, 2012. a

Kim, H.-J., An, S.-I., Kim, S.-K., and Park, J.-H.: Feedback processes modulating the sensitivity of Atlantic thermohaline circulation to freshwater forcing timescales, J. Climate, 34, 5081–5092, 2021. a, b, c

Kwok, R. and Untersteiner, N.: The thinning of Arctic sea ice, Phys. Today, 64, 36–41, 2011. a

Lavergne, T., Sørensen, A. M., Kern, S., Tonboe, R., Notz, D., Aaboe, S., Bell, L., Dybkjær, G., Eastwood, S., Gabarro, C., Heygster, G., Killie, M. A., Brandt Kreiner, M., Lavelle, J., Saldo, R., Sandven, S., and Pedersen, L. T.: Version 2 of the EUMETSAT OSI SAF and ESA CCI sea-ice concentration climate data records, The Cryosphere, 13, 49–78,, 2019. a

Leibowicz, B. D., Abbot, D. S., Emanuel, K. A., and Tziperman, E.: Correlation between present-day model simulation of Arctic cloud radiative forcing and sea ice consistent with positive winter convective cloud feedback, J. Adv. Model. Earth Sy., 4, M07002,, 2012. a

Lenton, T. M.: Arctic climate tipping points, Ambio, 41, 10–22, 2012. a

Li, C., Notz, D., Tietsche, S., and Marotzke, J.: The transient versus the equilibrium response of sea ice to global warming, J. Climate, 26, 5624–5636, 2013. a, b, c, d, e, f, g, h

Lindsay, R. and Schweiger, A.: Arctic sea ice thickness loss determined using subsurface, aircraft, and satellite observations, The Cryosphere, 9, 269–283,, 2015.  a

Mandel, P. and Erneux, T.: The slow passage through a steady bifurcation: delay and memory effects, J. Stat. Phys., 48, 1059–1070, 1987. a, b

Manoli, G., Fatichi, S., Bou-Zeid, E., and Katul, G. G.: Seasonal hysteresis of surface urban heat islands, P. Natl. Acad. Sci. USA, 117, 7082–7089, 2020. a

Nghiem, S., Rigor, I., Perovich, D., Clemente-Colón, P., Weatherly, J., and Neumann, G.: Rapid reduction of Arctic perennial sea ice, Geophys. Res. Lett., 34, L19504,, 2007. a

Notz, D.: The future of ice sheets and sea ice: Between reversible retreat and unstoppable loss, P. Natl. Acad. Sci. USA, 106, 20590–20595, 2009. a

Notz, D. and Stroeve, J.: Observed Arctic sea-ice loss directly follows anthropogenic CO2 emission, Science, 354, 747–750, 2016. a

Park, D.-S. R., Lee, S., and Feldstein, S. B.: Attribution of the recent winter sea ice decline over the Atlantic sector of the Arctic Ocean, J. Climate, 28, 4027–4033, 2015. a

Popp, M., Schmidt, H., and Marotzke, J.: Transition to a moist greenhouse with CO2 and solar forcing, Nat. Commun., 7, 1–10, 2016. a

Ridley, J. K., Lowe, J. A., and Hewitt, H. T.: How reversible is sea ice loss?, The Cryosphere, 6, 193–198,, 2012. a

Ritchie, P. D., Clarke, J. J., Cox, P. M., and Huntingford, C.: Overshooting tipping point thresholds in a changing climate, Nature, 592, 517–523, 2021. a, b, c

Stroeve, J. and Notz, D.: Changing state of Arctic sea ice across all seasons, Environ. Res. Lett., 13, 103001,, 2018. a

Stroeve, J., Serreze, M., Drobot, S., Gearheard, S., Holland, M., Maslanik, J., Meier, W., and Scambos, T.: Arctic sea ice extent plummets in 2007, Eos T. Am. Geophys. Un., 89, 13–14, 2008. a

Stroeve, J. C., Serreze, M. C., Holland, M. M., Kay, J. E., Malanik, J., and Barrett, A. P.: The Arctic's rapidly shrinking sea ice cover: a research synthesis, Climatic Change, 110, 1005–1027, 2012. a

Strogatz, S.: Nonlinear dynamics and chaos, Westview Press, ISBN 0-201-54344-3, 1994. a

Tredicce, J. R., Lippi, G. L., Mandel, P., Charasse, B., Chevalier, A., and Picqué, B.: Critical slowing down at a bifurcation, Am. J. Phys., 72, 799–809, 2004. a, b

Wagner, T. J. and Eisenman, I.: How climate model complexity influences sea ice stability, J. Climate, 28, 3998–4014, 2015. a, b, c

Short summary
We present a novel, efficient method for identifying climate tipping point threshold values of CO2 beyond which rapid and irreversible changes occur. We use a simple model of Arctic sea ice to demonstrate the method’s efficacy and its potential for use in state-of-the-art global climate models that are too expensive to run for this purpose using current methods. The ability to detect tipping points will improve our preparedness for rapid changes that may occur under future climate change.