Hysteresis of the Greenland ice sheet from the Last Glacial Maximum to the future

The Greenland Ice Sheet (GrIS) has undergone accelerated ice-mass loss in recent decades and it is expected to be one of the main contributors to global sea-level rise in the coming century. Due to the existence of positive feedbacks governing its mass balance, it is thought to be a tipping element of the Earth system. Its stability has been studied under temperatures ranging from the present day to a global warming of +4 K, showing a threshold behavior leading to an ice-free state for warmer temperatures. However, its stability at lower temperatures has not been studied yet. Here we use the ice-sheet model Yelmo to obtain the stability diagram of the GrIS for the full range of glacial-interglacial temperatures, with regional summer air temperature anomalies relative to present extending from a climate representative of the Last Glacial Maximum (−12 K) to a warmer climate (+4 K). We find that the hysteresis persists in almost the entire studied range. Consistent with previous studies, a critical threshold is found between +1.2 and +1.8 K of regional summer air temperature anomaly, associated with atmospheric feedbacks that are represented by the coupled regional energy-moisture balance model REMBO. In addition, a second threshold is found for negative temperature anomalies, which is mainly driven by ocean warming that triggers the marine ice-sheet instability in the northeastern region of the glacial GrIS. The existence of this threshold is consistent with transient studies of the GrIS over the last glacial cycle.

1 Introduction

The Greenland ice sheet (GrIS) has undergone an accelerated ice-mass loss in the last decades, partly as a result of global warming (Meredith et al., 2019; Hanna et al., 2021; Otosaka et al., 2023). There is a close agreement between the different methods used to estimate its mass loss, yielding a contribution of 14 ± 2 mm to sea-level rise between 1992 and 2020 (Otosaka et al., 2023). As summarized in the Fifth (Bindoff et al., 2014) and Sixth (Eyring et al., 2021) Assessment Reports of the Intergovernmental Panel for Climate Change (IPCC AR5 and AR6, respectively), several studies attribute both the temperature increase and the sea-level rise to anthropogenic forcing (Bamber et al., 2019; Fox-Kemper et al., 2021; Marcos and Amores, 2014; Meredith et al., 2019; Slangen et al., 2014).

Ice-mass loss, once triggered, is amplified by at least two positive feedbacks: the melt-elevation and the melt-albedo feedback. As a consequence, the GrIS responds nonlinearly to an increase of temperature and shows a threshold behavior: if a temperature threshold is exceeded, the GrIS can reach a qualitatively different equilibrium state, with strongly reduced ice volume. The GrIS is therefore considered to be a tipping element (Lenton et al., 2008). Ice-sheet modelling studies furthermore suggest that the GrIS shows multistability and hysteresis with respect to the temperature forcing (Solgaard and Langen, 2012; Robinson et al., 2012; Höning et al., 2023). As a result, if the rise in temperatures caused by anthropogenic climate change continues in the future, the GrIS could reach a new, ice-free equilibrium state, and it may not be able to recover its current extent and volume unless the temperature is decreased well below pre-industrial values. Therefore, the ice volume reduction could be irreversible on very long timescales. However, this does not mean that the changes would necessarily happen abruptly in time: due to the slow response of the ice sheets to atmospheric changes, variations in volume take relatively long times to manifest (Alley et al., 2005). The ice-free equilibrium state works as an attractor to which the GrIS approaches in global warming scenarios. The final state of the ice sheet depends on the rate of global warming and the time above the critical threshold (Bochow et al., 2023).

Despite the widespread agreement concerning the threshold behavior of the GrIS, the precise values of the critical threshold for the future remain highly uncertain (Fox-Kemper et al., 2021). Simulations indicate a critical mean annual global temperature threshold ranging between 1.6 (Robinson et al., 2012) and 3.1 K (Gregory and Huybrechts, 2006) above pre-industrial values. In addition, the hysteresis behavior of the GrIS has only been studied under temperatures above present-day values (Solgaard and Langen, 2012; Robinson et al., 2012; Bochow et al., 2023; Höning et al., 2023). Ice-sheet volume hysteresis for lower temperatures has only been studied in the context of the complete Northern Hemisphere and as a function of the insolation (Calov and Ganopolski, 2005; Abe-Ouchi et al., 2013).

Here, our main goal is to investigate the GrIS critical thresholds and stability across a range of temperatures, from conditions representative of the Last Glacial Maximum (LGM; ca. 21 kyr), when summer atmospheric temperatures in central Greenland were approximately 12 K below present-day values (Buizert et al., 2018), to future globally warmer conditions, with a summer regional anomaly of +4 K. To this end, we calibrate the ice-sheet model in order to correctly represent the GrIS both at the LGM and the present day. We investigate whether the hysteresis persists in colder climates and whether additional tipping points exist. Finally, we characterize the change in equilibrium states at each transition, and we study the dominant mechanisms along the stability diagram.

The work is structured as follows: in Sect. 2 we describe the the ice-sheet model Yelmo coupled to the regional energy-moisture atmospheric model REMBO used for our simulations and the experimental setup; in Sect. 3 we present and analyze the results of the experiments; in Sect. 4 we discuss the main results and relate them to the existing literature. Finally, the main conclusions of the work are given in Sect. 5.

2 Methods

2.1 Model

We use the ice-sheet model Yelmo (Robinson et al., 2020) coupled with the regional energy-moisture balance atmospheric model REMBO (Robinson et al., 2010). Yelmo is a three-dimensional thermomechanical ice-sheet model that has been previously used to simulate the evolution of the Greenland (Bochow et al., 2023; Tabone et al., 2024), the North American (Moreno-Parada et al., 2023) and the Antarctic ice sheets (Blasco et al., 2021; Juarez-Martinez et al., 2024). To compute the ice velocities, Yelmo uses the depth-integrated-viscosity approximation (DIVA) that considers longitudinal, lateral, and vertical shear stresses (Goldberg, 2011). In DIVA, the effective viscosity and the horizontal velocity gradients are replaced by their depth-integrated values in the effective strain rate equation. This allows for a more efficient two-dimensional solution, which yields close results to the computationally much more expensive Stokes problem (Robinson et al., 2022). A more detailed description of Yelmo is found in Robinson et al. (2020).

The basal friction is calculated using the regularized Coulomb friction law (Joughin et al., 2019) as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{b}}=c_{\mathrm{b}}{\left({\displaystyle \frac{\left|{\mathbf{u}}_{\mathrm{b}}\right|}{\left|{\mathbf{u}}_{\mathrm{b}}\right|+u_{\mathrm{0}}}}\right)}^q{\displaystyle \frac{{\mathbf{u}}_{\mathrm{b}}}{\left|{\mathbf{u}}_{\mathrm{b}}\right|}},\end{array}$$

where u_b is the basal velocity, q the friction law exponent, u₀ the regularization term for the Coulomb law set at 100 m yr⁻¹, and c_b the basal yield stress, defined as:

$$\begin{array}{}\text{(2)}& c_{\mathrm{b}}=\tan \left(\mathit{\varphi }\right)N.\end{array}$$

Here, N is the effective pressure, calculated following Bueler and van Pelt (2015). The till friction angle, ϕ, is a linear function of bedrock elevation, with a lower limit of 1° when the bedrock is 700 m below sea level and an upper limit of 40° when the bedrock is 700 m above sea level. As in Tabone et al. (2024), we impose a scaling factor of 0.1 on c_b in the northeastern basin. This accounts for the softer properties of the bed and allows for a better representation of the North East Greenland Ice Stream (NEGIS).

To represent glacial isostatic adjustment (GIA), Yelmo is coupled to the regional isostasy model FastIsostasy (Swierczek-Jereczek et al., 2024). We employ here an elastic lithosphere and a viscous mantle (ELVA) with a spatially homogeneous upper-mantle viscosity of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{21}}\text{Pa}\text{s}$ in the first 88 km and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{21}}\text{Pa}\text{s}$ in the following 600 km, and gravitational interactions between the ice sheet, the solid Earth and the ocean.

The surface mass balance (SMB) is calculated by the regional energy-moisture balance (REMBO) model, which follows an ITM scheme, described below in Sect. 2.2. The geothermal heat flow is taken from Shapiro and Ritzwoller (2004) and is imposed as a boundary condition at a depth of 2 km into the bedrock. Submarine melting at the base of ice shelves is parametrized by a linear equation dependent on ocean temperature (Sect. 2.3). Calving follows the Von Mises stress criterion (Morlighem et al., 2016) and it is defined as a function of the effective stress at the ice front (Lipscomb et al., 2019). While other calving criteria exist that represent the calving phenomenon with greater complexity, the computational cost of applying these methods at the full ice-sheet scale makes it preferable to employ simpler criteria such as Von Mises. However, this law produces satisfactory results consistent with previous studies of the GrIS using Yelmo (Tabone et al., 2024) and has demonstrated good performance in reproducing observed calving rates across different glaciers (Morlighem et al., 2016; Choi et al., 2018; Goelzer et al., 2017). In addition, since the presence of ice is considered implausible in the deep ocean, calving is applied where the bedrock depth exceeds 1000 m.

2.2 Atmospheric forcing

REMBO is a 2D regional energy-moisture balance climate model that simulates daily precipitation and temperature fields for Greenland, as well as the surface mass balance (SMB) through a one-layer snowpack model. The temperature and humidity over the ocean around Greenland are imposed as boundary conditions, for which the climatological mean from years 1958–2001 of ERA-40 reanalysis (Uppala et al., 2005) is used to represent the present-day conditions. The atmosphere over Greenland is then simulated through the vertically integrated energy balance and moisture balance equations. Additionally, REMBO takes into account the orography of Greenland and the planetary albedo, including its seasonal changes, thus reproducing the associated albedo and elevation feedback, the latter with a lapse rate of 6.5 K km⁻¹ (like in Robinson et al., 2012). It does not simulate the atmospheric broader general circulation and has its own biases, particularly in precipitation. Despite relatively low resolution of REMBO (100 km), the surface boundary conditions – and therefore the melt-elevation and melt-albedo feedbacks – are computed at the higher ice-sheet model resolution of 16 km. As a result, the model is able to represent the key atmosphere-ice feedbacks relevant to this study at a relatively high resolution. It has been previously used to simulate both the past and the future evolution of the GrIS (Robinson et al., 2012, 2011; Calov et al., 2015; Bochow et al., 2023).

The SMB is defined as the difference between the accumulation and the runoff. REMBO follows the insolation-temperature melt (ITM) method, which was originally developed by Pellicciotti et al. (2005) and van den Berg et al. (2008). The surface melt rate M_s follows the linear equation:

$$\begin{array}{}\text{(3)}& M_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}{L}_{\mathrm{m}}}}[c+\mathit{\lambda }T+{\mathit{\tau }}_{\mathrm{a}}(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}}\left)S\right],\end{array}$$

where ρ_w is the water density; L_m is the latent heat of ice melting; c and λ are empirical parameters set to $-\mathrm{55}{\mathrm{W}\mathrm{m}}^{-\mathrm{2}}$ and $-\mathrm{10}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{K}^{-\mathrm{1}}$; τ_a is the transmissivity of the atmosphere; α_s the surface albedo; and S the insolation at the top of the atmosphere. Therefore, the dependence of melting on both absorbed insolation and temperature is taken into account. The surface albedo is parametrized as a function of the snow thickness and the type of ground surface (Robinson et al., 2010).

REMBO is bidirectionally coupled to the Yelmo ice-sheet model. The SMB and surface temperature obtained by REMBO force the ice-sheet evolution in Yelmo and the evolution of the topography calculated in Yelmo is an input for REMBO (Robinson et al., 2010).

To force REMBO, we impose monthly temperature anomalies at the boundaries. Here, as in Robinson et al. (2012), we assume that annual temperatures follow a sinusoidal cycle where ΔT_DJF=2ΔT_JJA, so that the annual temperature anomaly is $\mathrm{\Delta }{T}_{\mathrm{ann}}=(\mathrm{\Delta }{T}_{\mathrm{DJF}}+\mathrm{\Delta }{T}_{\mathrm{JJA}})/\mathrm{2}=\mathrm{1.5}\mathrm{\Delta }{T}_{\mathrm{JJA}}$. Since temperature during the summer months is more important for melting, we use ΔT_JJA as input.

While relative humidity r is held constant, specific humidity depends on temperature according to

$$\begin{array}{}\text{(4)}& Q=Q_{\mathrm{sat}}\left(T\right)\cdot r,\end{array}$$

where Q_sat(T) is the saturation specific humidity, which follows the Clausius–Clapeyron relation. The choice to maintain constant relative humidity rather than specific humidity in REMBO is based on the stronger temperature dependence of the latter. Total precipitation is then calculated as a function of specific humidity, thus accounting for temperature variations:

$$\begin{array}{}\text{(5)}& P=(\mathrm{1}+k|\mathrm{\nabla }{z}_{\mathrm{s}}\left|\right){\displaystyle \frac{Q}{\mathit{\tau }}},\end{array}$$

where ∇z_s is the surface elevation gradient, τ is the water turnover time in the atmosphere, and k is an empirical parameter.

REMBO overestimates precipitation in the southwestern and northern regions (Robinson et al., 2010). To address this systematic bias, we apply a spatially varying bias correction following Wetterhall et al. (2012) and based on P_MAR, the 1961–1990 mean monthly precipitation from Tedesco et al. (2023). Thus, the corrected precipitation is:

$$\begin{array}{}\text{(6)}& P_{\mathrm{corr}}=\mathrm{P}\cdot \mathit{\delta }P,\text{where}\mathit{\delta }P={\displaystyle \frac{P_{\mathrm{MAR}}}{P_{\mathrm{REMBO},\mathrm{ref}}}}\end{array}$$

where P is the precipitation without correction, and P_REMBO,ref is the precipitation simulated by REMBO for the present day without correction. Gaussian smoothing is first applied to both P_MAR and P_REMBO,ref so that we only correct large-scale features. After this, the correction factor δP is computed internally within REMBO at its native resolution of 100 km. To maintain a consistent field with the boundary conditions at each moment and to avoid introducing a bias from present-day conditions, δP is restricted to values between 0.6 and 1.4 (see Fig. A3). This correction serves to improve the agreement of the seasonality and overall magnitude of precipitation, particularly in the north.

2.3 Oceanic forcing

The ocean forcing of the GrIS is treated in very different manners throughout modelling studies. Some fully neglect it, taking into account only atmospheric forcings (Stone et al., 2013) while others parameterize sub-shelf melting as a function of the ice-shelf draft (i.e., the ice-thickness below the ocean surface; Bradley et al., 2018). Still, others do not use an explicit forcing but constrain their model with relative sea level (Simpson et al., 2009; Lecavalier et al., 2014). Finally, some include parameterizations of basal melting as a function of ocean temperatures (e.g. Tabone et al., 2018). Results show that ice-ocean interactions play a key role in the past GrIS evolution, notably at times and in areas where the ocean-ice contact zone is largest, such as the LGM, when the GrIS margins reached the continental shelf break. Even if present marine-terminating glaciers cover only a small part of Greenland in the present day, their advance and retreat significantly affects the ice geometry and dynamics in the continental zone, and this effect is amplified as the ice-ocean contact zone expands (Tabone et al., 2018).

As in Tabone et al. (2018), we formulate the basal melting rate at the grounding line B_gl using an anomaly approach:

$$\begin{array}{}\text{(7)}& B_{\mathrm{gl}}\left(t\right)=\mathit{\kappa }\mathrm{\Delta }{T}_{\mathrm{ocn}}\left(t\right)+B_{\mathrm{ref}},\end{array}$$

where κ=15 m i.e. yr⁻¹ °C⁻¹ is the heat flux between the ice and the ocean at the ice-ocean interface; ΔT_ocn(t) is the anomaly between the ocean temperature at time t and at present day; and B_ref represents the basal melting rate at the present day. Here, B_ref is set at 50 m i.e. yr⁻¹ following observations of basal melting at the grounding line (Wilson et al., 2017). All parameters in Eq. (7) are assumed to be spatially uniform along the marine boundaries of the GrIS. Using a single value is a coarse approximation; however, in the absence of a complete sub-shelf melt map for Greenland, constructing a reliable 2D spatial field would be subjected to large uncertainties. Therefore, for the idealized experiments conducted in this study, we consider this approximation to be appropriate.

On the other hand, the sub-shelf basal melting rate B_sh is lower than that at the grounding line:

$$\begin{array}{}\text{(8)}& B_{\mathrm{sh}}=\mathit{\gamma }{B}_{\mathrm{gl}},\end{array}$$

where γ=0.1, in accordance with observations of Greenland glaciers (Wilson et al., 2017).

ΔT_ocn is related to the annual atmospheric temperature anomaly ΔT_ann as follows:

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{T}_{\mathrm{ocn}}=f_{\mathrm{ocn}}\mathrm{\Delta }{T}_{\mathrm{ann}}.\end{array}$$

The scaling factor between oceanic and annual atmospheric temperature f_ocn is an empirical parameter that is highly uncertain. Here we calculate the mean ratio between oceanic and atmospheric temperatures over the ocean surrounding Greenland using PMIP3 model outputs (CCSM4, CNRM-CM5, FGOALS-g2, IPSL-CM5A-LR, MIROC-ESM, MPI-ESM-P, MRI-CGCM3, Braconnot et al., 2012) and removing the outlier values. This yields $f_{\mathrm{ocn},\mathrm{LGM}}=\mathrm{0.18}\pm \mathrm{0.02}$ for the LGM and $f_{\mathrm{ocn},\mathrm{Hol}}=\mathrm{0.5}\pm \mathrm{0.2}$ for the mid-Holocene. These values show therefore considerable spread. We chose f_ocn=0.22 as it falls within this range and is closer to LGM values, which is appropriate since the ocean forcing has greater influence in the temperature range we investigate. Equation (9) can be rewritten as a function of ΔT_JJA:

$$\begin{array}{}\text{(10)}& \mathrm{\Delta }{T}_{\mathrm{ocn}}=\mathrm{0.22}\mathrm{\Delta }{T}_{\mathrm{ann}}=\mathrm{0.33}\mathrm{\Delta }{T}_{\mathrm{JJA}}.\end{array}$$

Therefore, considering the values from Eq. (7), basal melting in ice shelves (B_gl≥0) is activated at $\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{10.1}$ K.

Here, in contrast to Tabone et al. (2018) but following subsequent studies (Tabone et al., 2019a, b, 2024), B_gl is limited to positive values, preventing refreezing at the base of ice shelves. While local variations in basal melting and refreezing can occur in reality, given our simplification of applying spatially homogeneous melting, we also prevent refreezing at the base of ice shelves. This is consistent with recent studies in the Antarctic ocean showing that basal melting can still be active under glacial conditions (Obase et al., 2017; Kusahara et al., 2015).

2.4 Experimental design

To study the stability of the GrIS, warming (retreat) and cooling (regrowth) simulations (hereafter branches) were performed, starting at different initial states, respectively: an LGM-like state and a virtually ice-free state (Fig. 1). To obtain these initial states, we performed two spin-up experiments (one for each branch), initialized with present-day topography and ice thickness (Morlighem et al., 2017, 2022), and ran them for 60 kyr with a constant forcing. The ice-sheet model has ten vertical layers in sigma-coordinates and a horizontal resolution of 16 km, the same as in other recent studies examining the stability of the GrIS (Bochow et al., 2023; Höning et al., 2023). While a 16 km resolution may misrepresent certain processes, it enables long-term simulations that otherwise would be not possible due to their high computational cost (see Appendix Sect. A). At this horizontal resolution, the interaction with the ocean is not resolved in the narrowest fjords. For this reason, as in Tabone et al. (2024), the bedrock elevation is reduced by 1–1.5σ, where σ represents its standard deviation, in the areas in contact with the ocean and higher uncertainty in the value of the bedrock elevation.

The LGM-like state (Fig. 1a) was obtained by imposing a regional summer air temperature anomaly of $\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{12}$ K (Buizert et al., 2018) and an ocean temperature anomaly of $\mathrm{\Delta }{T}_{\mathrm{ocn}}=-\mathrm{4}$ K, according to Eq. (9). It is important to note that this LGM-like state is not meant to be a reconstruction of the LGM, but an idealized configuration designed to analyze the effect of temperature on ice-sheet stability in isolation from other major paleoclimatic forcings that are relevant to achieve a realistic LGM state, such as changes in sea level, insolation, or atmospheric CO₂. Nevertheless, even though the boundary conditions do not fully represent those of the LGM, we obtain a configuration that closely resembles simulations with full LGM conditions (see Appendix B). We obtain an ice sheet that is considerably larger than the present-day one. It has a simulated extension of 3.07 million km², a total volume of 5.15 million km³ and a volume above flotation of 4.6 million km³, corresponding to an anomaly of 3.66 m of sea level equivalent (SLE). The latter value is within the range reported in literature, ranging between 2.6 and 4.7 m SLE (Buizert et al., 2018; Lecavalier et al., 2014; Tabone et al., 2018; Bradley et al., 2018). Low atmospheric and ocean temperatures allow for the expansion of the GrIS up to the continental shelf break, with floating ice shelves appearing in the southeast and in Baffin Bay, fed by multiple ice streams. Large agreement is found between the simulated and the reconstructed margins for the LGM (Lecavalier et al., 2014). An exception is the northeastern area, where our simulation reaches the continental shelf break while the reconstruction by Lecavalier et al. (2014) does not. Whether the GrIS reached the shelf break in that region or not is still under debate (Evans et al., 2009; Arndt et al., 2017; Leger et al., 2024) but a more recent study based on high-resolution hydro-acoustic data suggests that was the case (Arndt et al., 2017), thus matching our simulations. What is almost certain is that the GrIS expanded beyond the limits suggested by Lecavalier et al. (2014), whose reconstruction for the LGM appears to be a lower limit as compared to the former studies.

The virtually ice-free state was obtained through the simulation of a warmer state (Fig. 1b), with a regional summer air temperature anomaly of +4 K and an ocean temperature anomaly of ca. +1.3 K. Under these conditions, the ice sheet is drastically reduced both in volume and extent (Fig. 1b). Ice persists only in the high-mountain areas: at the southern tip of Greenland and in the easternmost region, around the Watkins range. The mass balance of the remnant ice caps is almost zero in their interior (all incoming accumulation is melted away) and negative at their margins, especially in areas in contact with the ocean (not shown).

Finally, we performed a third spin-up for present-day conditions, for which we find a close agreement between our simulation and the observations (see Fig. A1 in the Appendix).

Figure 1 Final state of the spin-up simulation for (a) the LGM-like state (initial state of the warming branch), obtained through a regional summer air temperature anomaly of −12 K, and (b) a virtually ice-free state of +4 K (initial state of the cooling branch). The black contour lines indicate the ice-sheet surface elevation every 500 m starting from 0 and with a thicker line every 1000 m. In (a) the blue line indicates the reconstructed ice-sheet margin of Lecavalier et al. (2014); the green line is the maximum extent of grounded ice from the full glacial extent (18–16 kyr BP) in the northeast region by Leger et al. (2024). Note that the lightly shaded area in the northwest (over Ellesmere Island) indicates the part of the simulated ice sheet that is not taken into account for the volume calculation in order to be consistent with the usual GrIS domain.

To study the wider stability of the GrIS, we perform a series of simulations, outlined in Fig. 2. First, quasi-equilibrium simulations were performed with different rates of forcing (between $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\text{yr}}^{-\mathrm{1}}$) in a range of regional summer air temperature anomalies from −12 to +4 K relative to the average value between 1958–2001, which is used for the boundary conditions of the model as in Robinson et al. (2012) and Bochow et al. (2023). For the warming experiments, a positive constant temperature anomaly rate is applied to the LGM-like state, until a temperature anomaly of +4 K, relative to the average value, is reached. In parallel, for the cooling experiments a negative constant temperature anomaly rate is applied, starting from the virtually ice-free state until a temperature anomaly of −12 K is reached. Therefore, to achieve the full temperature anomaly range of 16 K from the initial state, the durations of the simulations vary between 3200 years and 1.6 million years, depending on the forcing rate.

In addition, equilibrium states are derived through two types of simulations: (1) applying an instantaneous summer air-temperature anomaly to the initial state during 60 kyr across different temperature levels ranging from −12 to +4 K, in 1 K increments; and (2) similar to Garbe et al. (2020), branching off from the $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{K}{\text{yr}}^{-\mathrm{1}}$ quasi-equilibrium simulation, then maintaining a constant temperature at different forcing levels (ranging from −12 to +4 K for the atmospheric temperature, with 1 K increments and 0.2 K increments near the tipping points) for 80, 120 or 400 kyr, depending on the time required for the simulations to stabilize.

By decreasing the forcing rates in the quasi-equilibrium simulation, we bring the transient simulations closer to the equilibrium states. Ideally, this would be achieved with an infinitely slow forcing rate, which is however impossible in practical terms. Nonetheless, the combination of the transient quasi-equilibrium simulations with the equilibrium states allows us to perform an exhaustive analysis that characterizes the stability phase space over the chosen range of study – where the tipping point values are taken from the branching-off experiments, as these are the most reliable due to their gradual convergence toward equilibrium.

Figure 2 Diagram of the forcing in the different simulations performed to obtain the stability diagram for the warming branch. The quasi-equilibrium simulations consist of ramping-up at different rates of forcing. In contrast, each equilibrium state is the final state (represented by a triangle) of one simulation with (1) instantaneous forcing and (2) branching off from quasi-equilibrium simulations with a rate of forcing of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{K}{\text{yr}}^{-\mathrm{1}}$. For the cooling branch simulations, the forcing is the same but with negative trends. The initial state of the three experiments is obtained from a steady-state spin-up simulation of the LGM-like (Fig. 1a) for the warming branch and a steady-state spin-up simulation under +4 K warming for the cooling branch (Fig. 1b).

Finally, in order to differentiate the atmospheric and the oceanic effects, we perform two more quasi-equilibrium experiments with a forcing rate of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}y{r}^{-\mathrm{1}}$ under the following conditions:

• Only OCN: ΔT_JJA remains constant at −12 K during the warming branch and +4 K during the cooling branch, corresponding to initial state values.

• Only ATM: ΔT_ocn remains constant at −4 K during the warming branch and +1.3 K during the cooling branch, corresponding to initial state values.

The initial states of both experiments are the same as in the original stability diagram (Fig. 1a for the warming branch and Fig. 1b for the cooling branch). We selected this forcing rate, as it allows us to maintain simulations close to equilibrium while reducing computational costs compared to slower forcing rates (these experiments require 533 kyr versus 1.6 Myr for the quasi-equilibrium simulation with a forcing rate of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$).

Figure 3 (a) Stability diagram of the GrIS: volume above flotation versus regional summer air temperature anomaly. The solid lines indicate the quasi-equilibrium simulations at different rates of forcing. The triangles indicate the final states of the instantaneous-forcing simulations (gray) and of the branching-off-from-quasi-equilibrium simulations (black); the direction of the triangles indicates the sense of the forcing (to the right in the warming or retreat branch and to the left in the cooling or regrowth branch). The light red shaded area above the warming branch represents the projection of the basin of attraction for the warming branch, i.e., the points for which the mass balance is negative under the given forcing and the ice sheet is expected to retreat until reaching the equilibrium state. The light blue shading indicates the analogous region for the cooling branch, identifying the points where the ice sheet is expected to grow. (b)–(f) GrIS configuration for the equilibrium states marked in (a). The blue thick lines in (b) and (c) indicate the ice-sheet margin of Lecavalier et al. (2014) and the green line is the maximum extent of grounded ice from the full glacial extent (18–16 kyr BP) in the northeast region by Leger et al. (2024). The red thick lines in (d)–(f) indicate the present-day margin of Morlighem et al. (2017). The black contour lines indicate the ice-sheet surface elevation every 500 m starting from 0 and with a thicker line every 1000 m.

3 Results

3.1 Bifurcation diagram

Through the different equilibrium and quasi-equilibrium experiments of the GrIS, we obtain two different branches (warming and cooling, or retreat and regrowth) that constitute the bifurcation diagram (Fig. 3a). As the forcing rates decrease (and therefore the forcing periods increase) the transient quasi-equilibrium response converges to the equilibrium states. This is a consequence of the large inertia of the ice sheet: only long forcing periods allow for the ice-sheet response to the forcing to fully develop and to approach the equilibrium states. In the slowest simulation ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$), which takes 1.6 million years, the volumes coincide with the equilibrium values for the entire interval, except in the vicinity of the tipping points, where the ice sheet takes even longer to stabilize.

We identify two bifurcation points in the warming branch. As noted earlier, their values are best determined from the branching-off equilibrium simulation, which provides the most reliable representation of the equilibrium states. The first one, between −10 and −9 K, begins with the loss of most of the GrIS ice shelves, followed by the retreat of the northeast margin (the state of the ice sheet before and after the transition is shown in Fig. 3b and c and is analyzed in greater depth in Sect. 3.3). Once a temperature anomaly of −9 K is reached, a gradual decrease in volume begins, showing an almost linear trend with a slope close to zero relative to the forcing. This behavior persists until a temperature anomaly of 1.2 K is reached (Fig. 3d), at which point the positive melt-albedo and melt-elevation (atmospheric) feedbacks are triggered (Sect. 3.2). This second tipping point is identified in this study between +1.2 and +1.8 K.

Figure 4 shows a zoomed-in view around this tipping point. Near the bifurcation point – particularly at +1.4 and +1.6 K – in both equilibrium branches of the branching-off experiments, the ice sheet oscillates between two states (shown in more detail in Fig. C1). These oscillations show periods ranging from 70 to 120 kyr. At +1.8 K, the GrIS reaches a stable virtually ice-free state with a volume of 1.58 m of SLE (Fig. 3e). In this transition, we move from an ice sheet that almost completely covers the surface of Greenland to a virtually ice-free state, with some residual ice only at the highest altitudes, mainly in the Watkins Range and the southern tip.

Figure 4 Enlarged section of the stability diagram shown in Fig. 3, focusing on the temperature range from −1 to 2.5 K. Solid lines indicate the slowest quasi-equilibrium simulation ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$). Red (blue) triangles represent the branching-off equilibrium states of the warming (cooling) branch. In cases with regional summer temperature anomalies of 1.2, 1.4, 1.6, and 1.8 K, the simulations show stationary oscillations and stabilize in a dynamic state. In these cases, unlike in Fig. 3, the triangles represent mean values, and the whiskers indicate the amplitude of the oscillation.

In the regrowth (cooling) branch, the GrIS starts at its virtually ice-free state (Fig. 1b), with $\mathrm{\Delta }{T}_{\mathrm{JJA}}=+\mathrm{4}$ K, and it is cooled. The ice sheet remains nearly ice-free until 1.6 K, where oscillations resume, as in the warming branch. The ice sheet remains in an intermediate state (Fig. 3f) until 0 K, with no ice in the north due to the high temperatures and the low precipitation in the area.

In this cooling branch, quasi-equilibrium simulations are farther away from equilibrium than in the warming branch for the same rates, especially during the critical transition (Fig. 3a, quasi-equilibrium simulations with rates of forcing between 10⁻⁵–${\mathrm{10}}^{-\mathrm{3}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$ and the temperature range 0–+1.8 K). Starting from the ice-free state, the recovery of the ice sheet requires lower temperatures and more time. This is due to the fact that melting or calving can occur almost instantaneously with high rates, while the processes of accumulation and regrowth are slower and can require thousands of years. Between 0 and −9 K the ice sheet does not show large variations, but just a linear regrowth with small sensitivity to temperature changes. Finally, from −9 K onwards, the sensitivity is increased again: the ice shelves are recovered, and the northeastern margin expands until it reaches the continental shelf break. By −12 K, it attains a volume above flotation of 4.45 million km³, almost converging to the initial state of the warming branch.

Hysteresis persists throughout almost the entire temperature range (between −10 and +1.4 K), with the largest volume difference for temperature anomalies between 0 and +1.4 K. It is important to highlight how, for the branching-off equilibrium states, the hysteresis is narrower and the regrowth in the cooling branch allows intermediate equilibrium states like that shown in Fig. 3f. At $\mathrm{\Delta }{\mathrm{T}}_{\mathrm{JJA}}=+\mathrm{1.2}$ K multistability becomes evident, with two different ice-sheet configurations (Fig. 3d and f) and respective volumes of 2.99 and 1.82 million km³. These results confirm the multistability induced by atmospheric feedbacks related to elevation and albedo identified in previous studies (Robinson et al., 2012; Höning et al., 2023), as well as the irreversible nature of this transition on long timescales: if the present-day GrIS collapses, even if temperatures return to present-day values, there would still be a 1.4 m SLE difference.

Within the range from −9 to 0 K, the differences between branches especially affect the ice thickness in the margins. For example, for a summer air temperature anomaly of −5 K (Fig. 5a–c), the equilibrium ice sheet on the warming branch extends well beyond the margin of the equilibrium state in the cooling branch for the same temperature and has a 12 % larger volume. This nonlinearity originates from different mechanisms. To explain these, we analyze the ice-sheet evolution at selected transects (1 and 2; Fig. 5d, e). These are located in areas with the greatest differences between the grounding lines, following the ice-flow streamlines. Along these transects, the Disko Island (Qeqertarsuaq in Greenlandic) and a bathymetric peak in the northeast act as pinning points for the warming branch equilibrium states (analogous to those found for the Antarctic ice sheet in  Favier et al., 2016). These slow down the outflow and allow the ice sheet to have a higher volume.

On the other hand, in the cooling branch the ice sheet shows a more irregular margin than in the warming branch. As it grows from the interior, it encounters the coastline (which is highly irregular) as a boundary, where ocean-induced melting occurs, preventing further growth. In contrast, in the warming branch, since it originates from an ice sheet in equilibrium and of greater extent (the LGM-like state), the continuity of the ice allows basal melting caused by the ocean to melt the margin more uniformly, resulting in a more compact and regular ice sheet shape. This leads to a larger perimeter in the cooling branch: the cooling branch state has a margin of ca. 40 000 km, whereas the warming branch state's margin is ca. 27 000 km long. As a result, even though some areas in the cooling branch have not yet reached the coastline, the total margin in contact with the ocean is greater (ca. 10 000 km vs. ca. 7000 km in the warming branch). This generates a higher basal melting and calving, making the regrowth difficult and preventing the ice sheet from reaching the size of the warming branch at the same temperature.

Figure 5 Differences between the two equilibrium states at $\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{5}$ K (from the branching-off experiments). (a) and (b) show the surface velocities and the grounding line for the warming (red line) and the cooling (blue line) branch, respectively. The black contour lines indicate the ice-sheet surface elevation every 500 m starting from 0 and with a thicker line every 1000 m. (c) Ice thickness difference between the two branches (cooling branch minus warming branch). (d) and (e) show the bedrock elevation (brown), the ice thickness (gray) and the sea level (blue line) for both branches at transects 1 and 2 indicated in (a)–(c).

3.2 Forcing mechanisms: the ocean and the atmosphere

In order to analyze the mechanisms leading to ice loss and recovery along the stability diagram, we show the surface, basal and lateral mass fluxes in Fig. 6. This analysis can only be performed with the quasi-equilibrium simulations, so we selected the simulation with the lowest forcing rate ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$), as it is the closest to equilibrium. As we saw in the previous section, for each quasi-equilibrium simulation, the system is farthest from equilibrium near the tipping points (shaded gray). This is reflected in the total mass balance (Fig. 6), which remains stable over the entire bifurcation diagram except in the proximity of the tipping points, where it becomes negative for the warming branch and positive for the cooling branch.

For the warming branch, the SMB increases between −12 and −4 K as higher temperatures allow for more precipitation. However, from −4 K onwards, ablation increases and the SMB gradually decreases. Finally, at ca. +1.5 K SMB drops abruptly and becomes negative, showing that atmospheric feedbacks have been triggered. The basal mass balance (BMB) is negligible between −12 and −10.1 K due to the fact that low ocean temperatures do not allow for sub-shelf melting, and to the small value of the melting under grounded ice. When temperatures rise enough to produce sub-shelf melting, an ice-shelf reduction, a margin retreat, and a decrease in the ice-ocean contact zones is induced (see Appendix D). This eventually causes a reduction in calving. BMB continues the increase with ocean warming until around −5 K. After that point, sub-shelf melting begins to decrease as the contact with the ocean strongly diminishes, becoming almost negligible by +2 K. Calving is highest in the temperature ranges where ice shelves are present (−12 to −10.1 K). It decreases, as the BMB does, when the contact with the ocean decreases.

Therefore, before the first transition (at ca. −9.5 K for this quasi-equilibrium simulation), SMB is increasing, so only the sub-shelf basal melting – which intensifies with rising ocean temperatures – can drive the response, directly triggering ice-mass loss and accelerating ice dynamics. In contrast, during the second transition (at ca. +2 K), the influence of the ocean diminishes, as BMB decreases in amplitude and therefore cannot account for the observed mass loss. Instead, SMB decreases sharply, implying a positive feedback mechanism where surface forcing becomes the dominant driver.

In the cooling branch, surface processes are responsible for the ice-sheet growth, with accumulation increasing and ablation decreasing as temperature decreases. Basal melting and calving also increase as both the ice sheet and its contact zones with the ocean increase. Calving, in particular, strongly increases at around −9 K, when ice shelves recover and temperatures are too low to allow for basal melt. The transitions in this regrowth branch are more gradual than in the retreat branch, since the ice growth processes are slower.

Figure 6 Mass balance contributions of the quasi-equilibrium simulation of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$. (a) warming; (b) cooling branch. The gray shaded bars indicate the regions in which, despite the slow forcing rates, the system deviates from equilibrium due to the proximity to the tipping points. The black dashed line represents the total volume, with its axis located on the right-hand side.

Figure 7 (a) Ice volume above flotation for three quasi-equilibrium simulations with a forcing rate of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$: (1) OCN-only simulation, with fixed $\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{12}$ K (corresponding to the LGM-like state) during the warming branch and $\mathrm{\Delta }{T}_{\mathrm{JJA}}=+\mathrm{4}$ K (corresponding to the virtually ice-free state) during the cooling branch; (2) ATM-only simulation, maintaining ocean temperature anomaly constant at −4 K (corresponding to the LGM-like state) during the warming branch and at +1.3 K (corresponding to the virtually ice-free state) during the cooling branch; and (3) simulation with both atmospheric and oceanic forcing, same simulation presented in Fig.  3 with the same rate of forcing. (b)–(d) The three snapshots of the OCN-only simulation marked with blue dots in (a). (e)–(g) The three snapshots of the ATM-only simulation marked with red dots in (a). The black contour lines in (b)–(g) indicate the ice-sheet surface elevation every 500 m starting from 0 and with a thicker line every 1000 m.

To further investigate this issue, in Fig. 7 we show the results of the simulations with atmosphere and ocean separately. Regarding the warming branch, the OCN-only simulation follows the default run almost exactly until $\mathrm{\Delta }{T}_{\mathrm{ocn}}=-\mathrm{1}$ K, when the margin retreats to the coast and reaches a constant state regardless of the increase in ocean temperature. There is no second tipping point in this simulation, which confirms that the latter was due to the atmospheric feedbacks. Conversely, the ATM-only simulation shows minimal changes until $\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{3}$ K ($\mathrm{\Delta }{T}_{\mathrm{ocn}}=-\mathrm{1}$ K). At around $\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{1}$ K the ice sheet loses the northeast sector, then retreats to the present-day margin and subsequently reaches a virtually ice-free state.

Regarding the cooling branch, in the OCN-only experiment the ice sheet remains in its virtually ice-free state as expected, given that although oceanic temperatures decrease, it is impossible for the ice sheet to grow without a decrease in atmospheric temperatures. In contrast, the ATM-only experiment follows the control run almost exactly until $\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{7}$ K, when in the control run the ice sheet recovers the marine sectors and in the ATM-only experiment the warm ocean does not allow the grounding line to advance.

These results confirm the findings shown in Fig. 6: ocean warming drives the abrupt transition at the first bifurcation point, while atmospheric warming triggers positive surface feedbacks, specifically melt-elevation and melt-albedo feedbacks (Robinson et al., 2012; Bochow et al., 2023; Höning et al., 2023; Petrini et al., 2025; Pattyn et al., 2018; Noël et al., 2021; Boers et al., 2025). In the absence of ocean forcing, an abrupt loss in the northeastern region still occurs, but substantially higher atmospheric temperatures are needed to initiate such a margin retreat.

3.3 Regional study

We next analyze the regional changes in ice-sheet volume by dividing the GrIS into three main regions: north, west, and southeast (Fig. 8). This division is based on grouping basins with similar topographical, atmospheric, and oceanic characteristics. The northern region is formed by the basins at high latitudes around the Arctic Ocean, including also the North East Greenland Ice Stream (NEGIS) area (basins 1 and 2 in Zwally et al., 2012); the western region comprises basins discharging ice into Baffin Bay (basins 6, 7 and 8 in Zwally et al., 2012); and the southeastern region includes basins with mountain chains, where the ice persists for higher temperatures (basins 3, 4 and 5 in Zwally et al., 2012).

Figure 8 Warming branch of the stability diagram, separated into three regions (defined in inset) adapted from Zwally et al. (2012)'s basins. (a) Total volume (grounded and floating) of each region and (b), (c) and (d) the corresponding mass fluxes for each region separately. Solid lines represent the quasi-equilibrium simulation with a rate of forcing of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$, the triangles in (a) show the equilibrium states from the branching-off experiments.

In Fig. 8, we can see that while the first tipping point has a regional character, affecting only the northern zone (the only one with $\mathrm{d}M/\mathrm{d}t>\mathrm{10}{\mathrm{Gt}\mathrm{yr}}^{-\mathrm{1}}$), the second one is a global tipping point that influences the entire ice sheet. Additionally, when comparing the mass fluxes in the three regions in the quasi-equilibrium simulation with the lowest forcing rate as above (Fig. 8b–d), we find that in the northern area there is very little accumulation throughout the phase space. In contrast, the southeast area has a higher elevation and the SMB is significantly higher, with positive values even after exceeding the second tipping point.

Figure 9 MISI in the northeast GrIS. (a) and (b) show the surface velocities before ($\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{9.4}$ K) and after ($\mathrm{\Delta }{T}_{\mathrm{JJA}}=-\mathrm{9.3}$ K) crossing the tipping point in the quasi-equilibrium warming branch simulation (rate of forcing of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$). Black and red lines indicate the grounding line in both states, and the black contours represent the ice-sheet surface elevation every 500 m and with a thicker line every 1000 m. (c) shows the bedrock elevation (brown), the ice thickness (gray) and the sea level (blue line) for the transect 1 indicated in (a) and (b) with a white line.

At about −10 K, the ocean warming triggers sub-shelf melting, leading to a gradual volume decline (see Appendix D). Volume losses are concentrated in the western and southeastern regions, where the ice at the margin is floating and exhibits a higher sensitivity to oceanic forcing at lower temperatures. While these zones experience nearly linear losses up to approximately +2 K, in the northern region – where the ice is grounded and its thickness is higher – the volume remains nearly constant until the temperature anomaly reaches ca. −9.5 K (in this quasi-equilibrium simulation), when the abrupt ice loss occurs. The increased basal melting generates an imbalance that causes grounding-line retreat. Before the tipping point (Fig. 9), the ice sheet remains grounded on a flat bedrock, but when retreat begins, the ice-sheet margin retreats across a mostly retrograde slope. There, the grounding line becomes prone to marine ice-sheet instability (MISI) as described by Schoof (2007) and studied in different areas of Greenland (Khan et al., 2020). This triggers enhanced ice flux toward the margin, progressive thinning, transition from grounded to floating conditions, and eventually slightly increased calving, as shown by the spike at −9.5 K in Fig. 8c (left shaded area). The major retreat is not compensated by precipitation, as accumulation is minimal at these latitudes. This is evident in Fig. 8c (north), where the SMB is lower than in the southeast and west (Fig. 8b and d), preventing the ice-sheet recovery in this region.

At the second tipping point, the quasi-equilibrium simulation shows that losses also occur slightly earlier in the northern and western regions, where the decline of the SMB is more notable than in the southeast (Fig. 8b to d). However, when the ice sheet reaches equilibrium for $\mathrm{\Delta }{T}_{\mathrm{JJA}}=+\mathrm{1.8}$ K, the feedbacks that have triggered the ice losses result in the collapse of the entire ice sheet. Therefore, once ice loss begins in these regions, if the forcing persists in time, the complete loss of the GrIS becomes inevitable. Finally, as we saw in Fig. 1b, at 4 K the ice only remains in the southeast zone.

4 Discussion

We have mapped the stability diagram of the GrIS for the full range of temperatures it would experience, from the Last Glacial Maximum to future warming. We show the existence of hysteresis for almost the whole temperature range. This indicates that, once the GrIS retreats, cooling back to the original temperature is not sufficient to recover the original ice-sheet configuration, in agreement with previous studies (Robinson et al., 2012; Höning et al., 2023).

In our diagram, we identified two bifurcation points: one for colder climates and another one for warmer climates. These tipping points divide the bifurcation diagram into three ranges, each with a distinct almost linear trend in ice-volume changes with respect to regional summer air temperatures. The transition from one regime to another depends on the rate and duration of the applied forcing once the threshold has been crossed. If the warming persists over time, the ice sheet will tend to the equilibrium state corresponding to that temperature anomaly. However, a rapid increase in temperature, if not sustained over time, does not lead to a total loss of ice. This allows for overshooting of the thresholds without losing the complete ice sheet as shown by Bochow et al. (2023). For example, when $\mathrm{\Delta }{T}_{\mathrm{JJA}}=+\mathrm{4}$ K is reached in the simulation with the fastest rate of forcing ($\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$, pink line in Fig. 3a), the threshold is crossed by several degrees but the ice sheet still has a (transient) volume of 4.21 million km³, even greater than the present-day value.

To our knowledge, the existence of a MISI-driven bifurcation between −10 and −9 K has not been described in the literature to date. Reconstructions of the last deglaciation investigate the transition from the LGM to the present (Lecavalier et al., 2014; Tabone et al., 2018; Leger et al., 2024). These are very different experiments from ours. Nevertheless, they provide useful information for comparative analysis and further support our findings. Lecavalier et al. (2014) show that the GrIS margin experienced a rapid general retreat at ∼16 to 10 kyr BP. It began gradually in the south, while in the north it was abrupt and occurred slightly later, with higher temperatures. These changes are also abrupt in the simulations of Tabone et al. (2018), additionally highlighting the importance of oceanic forcing in the transitions between glacial and interglacial periods. Winsor et al. (2015) further emphasize the critical role of oceanic forcing as the main driver of retreat in the southwest margin, where marine-terminating glaciers are prevalent. Although these are transient simulations that deviate from equilibrium, they also show an abrupt transition largely influenced by ocean forcing above a certain temperature. In addition, as in our diagram, the southern region shows the greatest vulnerability, retreating at lower temperatures than the northern one. Moreover, although the retreat in the northern margin begins somewhat later, it occurs more rapidly. In addition, Larsen et al. (2018) show that the northeastern margin has experienced large fluctuations in the last 45 kyr, showing a large sensitivity to changes in ocean temperatures. In our study, the retreat in this area is driven by the ocean forcing, which accelerates the ice flux and triggers the MISI mechanism along the regions where the bedrock has a retrograde slope. In studies of the last deglaciation, instability is not necessary to explain the rapid retreat of the northeastern margin, as the temperature increase after the Younger Dryas was rapid and very large. However, through our equilibrium and quasi-equilibrium experiments, we find that nonlinear behavior is indeed triggered by oceanic forcing.

The exact value of this bifurcation point is heavily influenced by the representation of the ocean-ice interaction. On the one hand, the activation of basal melting is strongly determined by the parameters κ and B_ref of Eq. (7) (see Appendix E). On the other hand, it also depends on the ocean temperature forcing ΔT_ocn, which in turn depends on the scaling factor f_ocn (Eq. 10) that carries substantial uncertainty, as discussed in Sect. 2.3. However, B_ref corresponds to observable physical data (Wilson et al., 2017), and f_ocn falls within the range of values from PMIP3 LGM and mid-Holocene reconstructions. Furthermore, this parameter combination produces a configuration capable of accurately representing the GrIS during both the LGM and the present day, lending confidence to our results.

The second bifurcation point identified here is located between +1.2 and +1.8 K of regional summer air temperature anomaly. Between these two states, the GrIS equilibrium states oscillate. Similar behavior on long time scales was previously shown by Zeitz et al. (2022), and attributed to the interaction between the positive elevation feedback and the negative isostasy feedback. We do not analyze this effect and its consequences in detail, as it is beyond the scope of this work, but it appears to follow the mechanisms thoroughly characterized by Zeitz et al. (2022). In the vicinity of the bifurcation point, the ice sheet is unstable, and in all cases, once the +1.2 K threshold is exceeded, it undergoes a substantial retreat.

The hysteresis width found here, in the range of 0–2 K, is somewhat narrower than in previous studies. Although this is not reflected in the quasi-equilibrium simulations, equilibrium states (black triangles in Fig. 3) display an ice regrowth that starts at higher temperatures than in previous studies. Between 0 and +1.4 K of regional summer air temperature anomaly, in the cooling branch, there are intermediate states for which the GrIS shows a partial recovery (Fig. 3f). Ice-sheet regrowth in other studies was found at +0.3 K (Robinson et al., 2012), or between −0.5 and +0.5 K depending on the level of previous warming (Höning et al., 2023). We have not investigated in depth the reasons for these differences, but it is conceivable that they are due to structural differences between models. For instance, Robinson et al. (2012) employed the SICOPOLIS ice-sheet model, which relies on a different approximation of ice dynamics (SIA vs. DIVA), a different sliding law parametrization, and a different treatment of GIA, the elastic lithosphere and relaxed asthenosphere (ELRA) vs. ELVA. Additional differences include the representation of basal hydrology and the inherent numerical schemes. Moreover, our bedrock topography currently includes much more detail than the older studies, which can also influence the ice-sheet stability.

Despite these differences, Robinson et al. (2012) also identified an intermediate branch (in addition to warming and cooling branches conceptually similar to ours) initialized with intermediate initial state conditions. This constitutes a third branch of equilibrium states that, while not simulated in our study, corresponds to intermediate ice-sheet configurations across the temperature range +0.4 to +1.6 K of regional summer temperature anomaly. Although this represents a different experimental setup, it is noteworthy that the intermediate states in both studies display a similar spatial pattern (covering southern and central Greenland) and comparable volume in relative terms. This demonstrates that the combination of albedo reduction following considerable ice retreat and low precipitation in northern Greenland hampers the ice recovery in this region.

Observations indicate that the western and southeastern regions have experienced the most significant ice-mass loss in recent decades, under ongoing global warming (Otosaka et al., 2023; Khan et al., 2020; Mouginot et al., 2019). In particular, Jakobshavn Isbrae has been the main contributor to the GrIS-related sea-level rise to date. The outlet glaciers in the northern margin are susceptible to greatly increasing their ice discharge, as their current velocities are slow due to their buttressing ice shelves (Mouginot et al., 2019). In recent decades, they have experienced a significant increase in their losses (Khan et al., 2022), and if they lose their ice shelves, their contribution to sea-level rise will be similar to that of the western and southeastern regions. Our simulations agree with these observations, illustrating that these areas are the first to show changes with warming (except for the southeast, where losses occur slightly later in our simulations). However, other studies (Höning et al., 2023) suggest a possible equilibrium state in which only the southern sector of the GrIS is lost. They identified two tipping points: one between +0.6 and +0.9 K (of global warming relative to pre-industrial), where only the southern GrIS melts, and another one between +1.6 and +2 K, where the ice sheet almost completely disappears, with ice remaining only in the northeast. In contrast, in our stability diagram, intermediate states during this transition are only found at +1.4 and +1.6 K (of summer regional warming), and they oscillate between 1.51 and 1.97 million km³ and between 0.60 and 1.60 million km³. These oscillations alternate between an intermediate ice-sheet configuratio – in which, due to low precipitation in the north and rising temperatures, ice disappears from the northern region – and a near ice-free state. Therefore, in our case, we do not consider this transition to involve two distinct tipping points.

To relate these values to global mean temperature anomalies with respect to the pre-industrial period, we follow Bochow et al. (2023) and apply the scaling relationship $\mathrm{\Delta }{\mathrm{GMT}}_{\mathrm{PI}}=f\cdot \mathrm{\Delta }{T}_{\mathrm{JJA}}+\mathrm{0.5}\mathrm{K}$, where $f=\mathrm{1}/\mathrm{1.19}$ K⁻¹ represents the mean value derived from the CMIP6 SSP5-8.5 experiments (Bochow et al., 2023). We use the SSP5-8.5-derived factor rather than that from historical simulations because our target regional temperature anomalies are above present-day values, making this relationship more appropriate for future warming scenarios. Note that we apply this conversion only to the second bifurcation point, as it would not be appropriate at glacial temperatures. With this approach, we obtain a tipping point between +1.5 and +2 K, a value within the range reported in some previous studies (Höning et al., 2023; Bochow et al., 2023; Robinson et al., 2012) and below that of others (Noël et al., 2021; Gregory and Huybrechts, 2006). This temperature anomaly will likely be exceeded during this century unless drastic measures are taken (IPCC, 2021), and even the goal of the Paris Agreement to limit global warming to +2 K (+1.5 K if possible) does not guarantee safety (McKay et al., 2022).

However, equilibrium states should not be confused with future projections. As we have seen, the transient convergence of the GrIS to its equilibrium states depends on the rate and time of forcing. Nevertheless, these states act as attractors, providing an approximation of the ice sheet's tendency under a given forcing scenario. Thus, if anthropogenic climate change continues and this temperature threshold is exceeded for a sustained period of time, the GrIS could be subject to the largest bifurcation and reach a virtually ice-free state over the next millennia.

5 Conclusions

In this study, we perform a stability analysis of the GrIS, incorporating both quasi-equilibrium simulations and equilibrium states. Covering a range of ice-sheet states, from LGM-like to near ice-free conditions, we identify two tipping points: one between −10 and −9 K and the other one between +1.2 and +1.8 K of regional summer air-temperature anomaly. For the first tipping point, the ocean plays a dominant role in driving ice loss, particularly through basal melting in combination with MISI in some locations. In later stages, as temperatures rise and the GrIS crosses the first tipping point, the diminishing contact between the ice sheet and the ocean shifts the primary drivers of mass loss to atmospheric processes. However, throughout this entire phase between −9 and +1.2 K, the GrIS shows almost no sensitivity to temperature increases and exhibits a linear behavior with a slope close to zero. Once temperature anomalies exceed +1.2 K, atmospheric feedbacks are triggered and the ice sheet undergoes a drastic reduction, eventually stabilizing at a nearly ice-free state with only residual ice persisting in the highest-altitude regions. Consistent with previous work, this transition suggests that exceeding this temperature threshold under ongoing anthropogenic climate change would result in a virtually ice-free GrIS, with severe implications for global sea-level rise.

There is substantial uncertainty in the value of the first bifurcation point due to its dependence on ocean forcing and the basal melting parameterization, both of which are highly sensitive to parameter choice. However, despite this uncertainty in the exact threshold value, the retrograde character of the bedrock creates an unstable region for the grounding line. This instability therefore has a regional character and affects only the northeast sector of the ice sheet. In contrast, the second one, although beginning at slightly lower temperatures in the west and north in the quasi-equilibrium simulations, has a continental-scale reach. Once this second threshold is exceeded, the GrIS eventually loses practically all of its volume. Conversely, the regrowth of the ice sheet is much slower and requires colder temperatures, illustrating the multistability in the processes of ice loss and recovery. Hysteresis persists across almost the entire temperature range. Its amplitude is greatest between 0 and 1.4 K, due to climate feedbacks.

Throughout this study, we have mapped the stability of the GrIS for temperatures corresponding to the last glacial cycle and those expected over the coming centuries. In this way, two distinct mechanisms affecting its stability – and giving rise to two different bifurcation points – have been identified, further underscoring the importance of considering ice-ocean interactions in colder climates. The bifurcation point that we found for future anthropogenic warming is within the lower range of previous estimates and even if the Paris Agreement is met, this limit will be surpassed in the coming decades.

Appendix A: Present-day performance

In the following, we present the comparison of the Yelmo-REMBO simulation for present day conditions to observations (Fig. A1). As the initial states of the bifurcation diagram (Fig. 1), this third spin-up is initialized with present-day topography and ice thickness (Morlighem et al., 2022), but it is forced with a constant climate corresponding to the average conditions between 1958 and 2001, running for 60 kyr.

Moreover, Fig. A2 shows the surface mass balance simulated by Yelmo-REMBO compared to that simulated by MAR at 10 km resolution (Gallée and Schayes, 1994; Grailet et al., 2025), version v3.14.3, forced at its lateral boundaries by ERA5. These simulations cover the period 1958–2001 and represent the latest results from a state-of-the-art model for surface mass balance over the GrIS (Xavier Fettweis, personal communication, 2025).

For the present-day conditions, we obtain a volume of 3.16 million km³, which is 5 % higher than present-day observations (Morlighem et al., 2017). The main discrepancies in the ice-surface elevation (Fig. A1) are in the southwestern margin, where the ice margin is too advanced, and along the northeastern and eastern coasts, where ice thickness is overestimated. In addition, this excess ice, which starts from the margins, spreads slightly inland. The overestimation of ice in the north is a common problem in simulating the GrIS, as noted in several studies (Stone et al., 2010; Robinson et al., 2012; Born and Nisancioglu, 2012; Tabone et al., 2018; Höning et al., 2023). Additionally, REMBO shows an excess in accumulation in the northeast compared to other regional models (Robinson et al., 2010), as reflected in Fig. A2 by higher SMB values in that region relative to the MAR results. Even though a scaling correction has been applied to the precipitation, small changes in atmospheric forcing can lead to significant changes in ice-sheet evolution (Niu et al., 2019). Moreover, the grid size of 16 km results in misrepresentations of ice-ocean interactions in fjords, where the ocean penetrates the continent through narrower channels than the grid resolution allows. Therefore, the present-day simulation shows ice extending slightly beyond its present day margin. The same applies in the east, where the topography shows abrupt variations. Regarding surface ice velocities, our model performs well against the observations. The only exceptions are the Northeast Greenland Ice Stream (NEGIS), which is generally challenging for ice sheet models to represent, and the innermost regions of glaciers, where our model reduces velocities too early when moving inland. Even so, this is a very good overall representation for simulations in which the basal friction has not been optimized.

The aim of this study is not to achieve a perfect representation of the present day ice sheet, but to investigate its stability and thresholds over a broad range of temperatures. The experimental set-up allows both a representation of the present-day and the LGM-like state of the ice sheet which is adequate for this task. The choice of a 16 km grid size means that some processes may not be represented. However, it allows for quasi-equilibrium simulations extending up to a million of model years, which would not be computationally possible with a higher resolution grid.

Finally, we present the annual precipitation from this simulation to compare the original field with the corrected one obtained through the bias correction procedure described in Sect. 2.2. Figure A3 shows both precipitation fields at the downscaled resolution of 16 km used to calculate the surface mass balance applied in Yelmo. In annual terms, the corrected precipitation field shows higher precipitation in the southern tip and lower precipitation in the north and in the southwestern margin.

Figure A1 Present-day performance. First row: ice-surface elevation, observations from Morlighem et al. (2022) at 16 km resolution. Second row: the ice-surface velocities, observations from from Joughin et al. (2018). The black contour lines indicate the ice-sheet surface elevation every 500 m starting from 0 and with a thicker line every 1000 m.

Figure A2 Surface mass balance in (a) Yelmo-REMBO; and (b) MAR version 3.14.3 forced with ERA5, mean for the period 1958–2001. Values are expressed in meters of water equivalent per year.

Figure A3 Annual precipitation from REMBO in the present-day spin-up simulation in millimeters of water equivalent per year. Panel (a) shows the original field of precipitation and panel (b) shows the corrected field. The black contour lines indicate the ice-sheet surface elevation every 500 m starting from 0, with a thicker line every 1000 m.

Appendix B: Sensitivity of the LGM-like state to sea level and insolation

To explore the GrIS stability with respect to temperature changes in the glacial-interglacial range of temperatures, we use a simulation with LGM-like conditions as the initial state for the warming branch. This simulation corresponds to the final state of a 60 kyr spin-up simulation forced with the LGM atmospheric temperature (Buizert et al., 2018). However, we maintain insolation and sea level at present-day (PD) values to isolate the effect of temperature changes on GrIS stability. The validity of this simulation as an LGM representation may be questioned, as it does not incorporate all LGM boundary conditions. Here, we present additional spin-up simulations with different LGM boundary conditions to examine their effects on the resulting initial state.

According to the relative sea level (RSL) reconstruction of Waelbroeck et al. (2002), RSL during the LGM was 120 ± 13 m below present-day values. During the subsequent 5 kyr, regional summer air temperatures increased by only approximately 2 K (Buizert et al., 2018), while sea level rose to 80 m below present-day values. Given this substantial change in boundary conditions with minimal temperature variation, we performed six experiment permutations under identical temperature forcing (−12 K):

• PD insolation with PD sea level (the simulation presented as the LGM-like state in the original manuscript)

• PD insolation with sea level at −120 m

• PD insolation with sea level at −80 m

• LGM insolation with PD sea level

• LGM insolation with sea level at −120 m

• LGM insolation with sea level at −80 m

Figure B1 shows the equilibrium cold states for each case under −12 K regional summer air temperature forcing. The primary difference among simulations is the extent of floating ice along the western and southeastern margins. The simulation with LGM insolation and lowest sea level (Fig. B1e) most closely represents full LGM boundary conditions. Under these conditions, most of the previously floating ice becomes grounded, yielding the highest volume above flotation of 5.06 million km³, equivalent to +4.81 m SLE compared to present day, and consistent with previous estimates (Buizert et al., 2018; Lecavalier et al., 2014; Tabone et al., 2018; Bradley et al., 2018).

Importantly, ice-sheet geometry remains broadly similar across all simulations and lies within reconstructed paleoclimatic constraints. The northeastern margin, where the first tipping point occurs, exhibits nearly identical configuration in all cases.

Since the simulation in Fig. B1e represents a more complete LGM reconstruction with full representative boundary conditions, and given the similar configuration of the simulation in Fig. B1a, we can characterize the latter as an “LGM-like state”. While not identical to a full LGM reconstruction, it shows similar ice-sheet geometry and stability characteristics, making it a suitable initial condition for the warming branch of the hysteresis experiments.

Figure B1 Final state of spin-up simulations under −12 K regional summer air temperature anomaly forcing. The upper row shows simulations with present-day insolation; the lower row shows simulations with 20 kyr BP insolation. The first column represents present-day global sea level conditions, while the second and third columns show simulations with global sea levels 80 and 120 m below present-day values, respectively. Volume above flotation (in millions km³) is displayed at the bottom of each panel. The black contour lines indicate the ice-sheet surface elevation every 500 m starting from 0 and with a thicker line every 1000 m. The blue lines indicate the ice-sheet margin of Lecavalier et al. (2014) and the green line is the maximum extent of grounded ice from the maximum full glacial extent (18–16 kyr BP) in the northeast region by Leger et al. (2024).

Appendix C: Branching-off equilibrium states

Figure C1 Results of the branching-off simulations. Panels (a) and (b) represent the warming branch, while (c) and (d) show the cooling branch. The first row displays the forcing: the summer regional air-temperature anomaly, and the second row shows the volume above flotation. Starting from the quasi-equilibrium simulation with a forcing rate of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$, each simulation is extended for an additional 80 kyr after reaching the corresponding temperature level – except for the +1.2 and +1.8 K simulations, which are extended for 120 kyr, and the 1.4 and 1.6 K simulations, which are extended for 400 kyr.

Appendix D: Surface and basal mass balance in the first bifurcation point

In this appendix, we examine in detail the ice mass changes at the surface and base along the first bifurcation point. Figure D1 shows the grounded and floating area in the warming branch of the quasi-equilibrium simulation at $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$ (black line in Fig. 3). After a very small loss of grounded ice, there is a significant loss of floating ice, followed by an abrupt retreat of grounded ice and the grounding line. Figure D2 shows the evolution of ice-sheet surface and basal mass balance before the first tipping point (Fig. D2a–e and f–i, respectively) and immediately after (Fig. D2e and j). Basal melting occurs at the grounding line, at the base of ice shelves, and in some regions beneath the ice sheet where a water layer exists. As temperatures increase, basal melting intensifies both at the grounding line around the ice-sheet margin and below the ice shelves. This causes ice shelf retreat until they remain only in reduced form in the southeast. Eventually, basal melting becomes sufficiently intense to trigger abrupt retreat of the northeast sector.

Figure D1 Ice area (grounded and floating) in the warming branch of the quasi-equilibrium simulation of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$ in the temperature range of the first bifurcation point.

Figure D2 (a)–(e) Surface mass balance and (f)–(j) basal mass balance at different temperatures from the warming branch of the quasi-equilibrium simulation of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$ before the first tipping point. The black contour lines indicate the ice-sheet surface elevation every 500 m starting from 0 and with a thicker line every 1000 m. Note that the higher values in basal melting in (f)–(j) highlight the grounding line.

Appendix E: Sensitivity to the basal melting parameters

The first bifurcation point, being triggered by ocean warming, is strongly influenced by the parameters of the basal melting at the grounding line, κ and B_ref, following the equation:

$$\begin{array}{}\text{(E1)}& B_{\mathrm{gl}}=\mathit{\kappa }\mathrm{\Delta }{T}_{\mathrm{ocn}}+B_{\mathrm{ref}}\end{array}$$

As explained in the main text, ΔT_ocn is the ocean temperature anomaly relative to the present-day value, and it is related to the regional summer air temperature anomaly by $\mathrm{\Delta }{T}_{\mathrm{ocn}}=\mathrm{0.22}\mathrm{\Delta }{T}_{\mathrm{ann}}=\mathrm{0.33}\mathrm{\Delta }{T}_{\mathrm{JJA}}$.

In Fig. E1, it can be seen that as κ increases, higher temperatures are required to activate basal melting, allowing both the ice shelves and the northeastern margin to persist at higher temperatures. The opposite occurs with B_ref: as it increases, the tipping-point value decreases. In these two ensembles, parameter values were chosen to allow the GrIS to grow at LGM conditions while remaining consistent with observations (Wilson et al., 2017). This results in a range for the first bifurcation point between −11 and −8 K, making it highly sensitive to the parameterization of ocean interactions. However, despite the different temperatures at which this bifurcation point occurs, the process is the same across all simulations: the loss of ice shelves followed by an abrupt retreat in the northeast, accelerated by the MISI. In agreement with our mechanistic attribution of the bifurcation points, the second one is not affected by changes in the ocean parameterization.

Figure E1 Warming branch quasi-equilibrium simulations (rate of forcing of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$) for different values of κ and B_ref. (a) The total ice volume (grounded and floating) and (b) the basal mass balance.