Role of elevation feedbacks and ice sheet–climate interactions on future Greenland ice sheet melt

Abstract

The Greenland ice sheet (GrIS) stores freshwater equal to more than 7 m of potential sea level rise (SLR) and strongly interacts with the Arctic, North Atlantic and global climate. Over the last few decades, the ice sheet has been losing mass at a rate that is projected to increase. Interactions between the GrIS and the climate have the potential to amplify or reduce GrIS mass balance responses to ongoing and projected warming. Here, we investigate the impact of ice sheet–climate interactions on the climate and mass balance of the GrIS using the Community Ice Sheet Model version 2 (CISM2) coupled with the Community Earth System Model version 2 (CESM2). To this end, we compare two idealized multi-century simulations with a non-evolving and evolving ice sheet topography in which we apply an annual 1 $\mathrm{\%}$ increase in ${\mathrm{CO}}_{\mathrm{2}}$ concentrations, starting from pre-industrial (PI) until stabilization at $\mathrm{4}\times$PI ${\mathrm{CO}}_{\mathrm{2}}$ concentrations ($\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$). By comparing the one- and two-way coupled simulations, we find significant changes in atmospheric blocking, precipitation and cloud formation over Greenland as the GrIS topography evolves, acting as negative feedbacks on mass loss. We also attribute part of the overestimation of mass loss in the one-way coupled simulation to an overestimation of melt in the ablation area caused by the use of a uniform temperature lapse rate to reflect the elevation differences between the atmospheric and ice sheet grids. Furthermore, we investigate ice sheet–climate interactions in a simulation branched in year 350 from our two-way coupled simulation in which we annually reduce atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ by 5 $\mathrm{\%}$ until PI concentrations are reached. During the 350-$\mathrm{year}$ $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ forcing period, the ice sheet loses a total mass of 1.1 $\mathrm{m}$ sea level equivalent, and part of its margins retreat landward. When the PI ${\mathrm{CO}}_{\mathrm{2}}$ concentration is restored, melt decreases rapidly, leading to a small positive surface mass balance. Combined with the strongly reduced ice discharge resulting from the widespread retreat of the ice sheet margin, this halts GrIS mass loss despite a remaining global warming of 2 $\mathrm{K}$. The GrIS, Arctic and North Atlantic strongly interact, causing a complex transitional phase towards a colder climate during the century following the ${\mathrm{CO}}_{\mathrm{2}}$ reduction. Elevated atmospheric temperatures, larger ocean heat transport and deteriorated state of the snowpack, compared to the initial pre-industrial state, result in limited regrowth of the ice sheet under reintroduced PI ${\mathrm{CO}}_{\mathrm{2}}$ conditions.

1 Introduction

Over the past century, the rate of global mean sea level rise (SLR) has exceeded any previous period and has been accelerating since the late 1960s . Melt of the Greenland ice sheet (GrIS), storing more than 7 m of potential SLR , has caused $\sim \mathrm{18}$ $\mathrm{\%}$ of contemporary SLR , and its contribution is expected to increase . The Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change IPCC; attributes the largest uncertainties in future sea level rise to the Greenland and Antarctic ice sheets. Therefore, understanding the physical drivers for GrIS mass loss and adequate numerical modeling of their effects on the GrIS mass balance are crucial to provide reliable projections of global and regional SLR.

Contemporary climate change, caused by anthropogenic greenhouse gas emissions , influences surface mass processes over the GrIS in connection with rising Northern Hemisphere summer temperatures . Moreover, changes in the GrIS have the potential to influence local and global climate because of their interaction with key components of the Earth system through various processes and feedbacks . Important feedback mechanisms include the positive melt–elevation feedback , the positive melt–albedo feedback and the negative melt–discharge feedback . The strong coupling between the GrIS and the atmosphere potentially influences atmospheric circulation and precipitation patterns . The GrIS heavily interacts with the ocean, as the North Atlantic Meridional Overturning Circulation (NAMOC) influences GrIS temperatures, while freshwater influxes as a result of GrIS mass loss can influence the strength of the NAMOC . Accounting for these interactions between the Earth system and the GrIS in climate models is critical to obtain reliable projections of SLR and GrIS evolution . It has been shown that one-way coupled simulations, in which changes in the state of the GrIS are not communicated to the Earth system model (ESM), overestimate multi-centennial GrIS mass loss as several feedbacks are not represented in these simulations. identify more precipitation and less melt in a two-way coupled simulation, in which the GrIS is an interactive component within the Earth system model and attribute this to a change in atmospheric circulation patterns over Greenland. In contrast, attribute the reduced mass loss to an increase in cloud fraction, causing an increase in reflected shortwave radiation. Next to that, they find a land-inward migration of precipitation patterns caused by topographic changes.

Considering that by 2022, global temperatures have increased by 1.09 °C (0.95–1.20) compared to the 1850–1900 baseline , reaching the 1.5 °C goal of the Paris Agreement is becoming less likely over time . The rapidly increasing global temperatures call for the investigation of “overshoot” scenarios in which this temperature threshold is surpassed. Such an overshoot could have important implications for the evolution of the GrIS-induced SLR as global mean sea level (GMSL) could rise substantially under the higher temperatures during such an overshoot period . Due to the nonlinear and tipping nature of the Earth system, there might be non-reversible effects associated with an overshoot scenario. Consequently, investigating the response of the GrIS to these kinds of scenarios is of great importance for assessing the implications of following an overshoot pathway. Previous work shows that temperature and volume thresholds for irreversible GrIS mass loss might exist. These model-based studies show that if temperature overshoots are limited and the ice sheet does not lose a critical ice volume, mass loss can be halted or might even be reversible. However, a thorough model-based assessment of the role played by ice sheet–climate feedbacks in the reversibility of enhanced deglaciation rates is currently lacking. As ice sheet–climate interactions could potentially accelerate or slow down the changes in GMSL caused by the response of the GrIS to ${\mathrm{CO}}_{\mathrm{2}}$ reduction, it is important to account for them to obtain a reliable projection of SLR in an overshoot scenario.

Considering these ice sheet–climate feedbacks and interactions can be done by coupling an ice sheet model with an Earth system model or a regional climate model. Over the past few years, there has been a lot of development regarding these kinds of models e.g.,.

In this paper, we use a coupled ice sheet–Earth system model setup to assess the impact of ice sheet-climate interactions on the mass balance of the GrIS. We evaluate the impact of these interactions under an idealized extreme warming scenario, in which we annually increase atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ by 1 $\mathrm{\%}$ from pre-industrial (PI) concentrations until 4 times the PI ${\mathrm{CO}}_{\mathrm{2}}$. We begin by comparing two simulations with this idealized forcing, of which one allows for changes in the GrIS topography to be communicated to the Earth system model (two-way coupling) and one does not (one-way coupling). We then branch a third simulation from the two-way coupled simulation in year 350 (i.e., after 2 centuries of sustained $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ conditions), where we bring back PI ${\mathrm{CO}}_{\mathrm{2}}$ concentrations with an annual 5 $\mathrm{\%}$ decrease (covering 27 $\mathrm{years}$), to assess the impact of ice sheet–climate interactions on the response of GrIS deglaciation rates to ${\mathrm{CO}}_{\mathrm{2}}$ reduction.

In Sect. 2, we describe the coupled model setup, the simulations conducted and our analysis methods. Section 3 assesses the influence of the coupling on the simulated GrIS mass loss, and Sect. 4 gives an overview of the most important elevation feedbacks that influence the GrIS mass balance. In Sect. 5, we investigate other ice sheet–climate interactions and address GrIS mass balance response resulting from the reintroduction of PI ${\mathrm{CO}}_{\mathrm{2}}$ conditions. We discuss our results and draw conclusions in Sects. 6 and 7, respectively.

2 Method

2.1 Model description

We use the Community Earth System Model version 2 CESM2; and the Community Ice Sheet Model version 2 CISM2;, which are coupled to account for an evolving ice sheet. CESM2 is a state-of-the-art community-developed ESM, developed by the US National Science Foundation (NSF) National Center for Atmospheric Research (NCAR). CESM2 can produce a realistic surface mass balance SMB; over the GrIS using a surface energy balance (SEB) scheme. CESM2 consists of different component models for land and land biochemistry, atmosphere, river runoff, surface waves, ocean and marine biochemistry, sea ice, and land ice, which are coupled to exchange states and fluxes via a hub-and-spoke architecture .

The atmosphere model is the Community Atmosphere Model version 6 CAM6;, which runs on a nominal 1° (1.25° in longitude and 0.9° in latitude) grid and has 32 vertical levels. The Community Land Model version 5 CLM5; shares the horizontal grid of the atmosphere model. In this land model, every grid cell is divided into one or multiple fractions of land units, which can be glacier, lake, wetland, urban, vegetation and crop surface. Calculations are carried out separately over the different land units. The model has a fixed number of vertical layers for the soil, whereas there is a variable number of layers for snow and firn, with a maximum snow depth of 10 $\mathrm{m}$ water equivalent (w.e.). The model allows for the compaction of snow into firn. Accumulation of snow over 10 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ in a grid cell is transferred as positive SMB to the ice sheet model. If snow and firn are melted away, further melt of ice is transferred as negative SMB (ice ablation) to the ice sheet model. Changes within the 10 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ snow layer are not communicated to the ice sheet model and are only considered in the land model. The snow albedo is calculated with the SNow, ICe and Aerosol Radiation Model SNICAR;. The surface runoff from melt and rain is routed to the ocean using the Model for Scale Adaptive River Transport MOSART;. Modeling the ocean is done using the Parallel Ocean Program version 2 POP2;. POP2 has a resolution of a nominal 1$\mathrm{°}$, with a uniform resolution equal to 1.125° in the zonal direction. The sea ice model is the Los Alamos National Laboratory sea ice model version 5 CICE5;. This component shares the horizontal grid with the ocean model. The model consists of an ice thermodynamics model, an ice dynamics model, and a transport model that computes both horizontal transport and transport in thickness space.

The GrIS is modeled using the Community Ice Sheet Model version 2 CISM2;. CISM2 is a parallel 3D thermodynamics model, which solves the momentum balance and computes the ice sheet thickness and temperature . This ice sheet model has a 4 $\mathrm{km}$ rectangular grid, with 11 vertical sigma levels. The model uses approximations of the Stokes equations for incompressible viscous flow and a pseudo-plastic sliding law to parameterize basal sliding. Floating ice at marine margins is immediately discharged to the ocean using a flotation criterion .

2.2 Coupling description

By applying two-way coupling between CISM2 and CESM2, interactions and feedbacks between the ice sheet and the climate are accounted for in the projected evolution of GrIS mass loss. The coupling of the ice sheet model is bi-directional with the land and atmosphere models and uni-directional with the ocean model (Fig. ) and has an annual frequency, except for the communication of the GrIS topography to the atmosphere model, which is done every 5 or 10 model years. This means that the ocean model will receive updated freshwater fluxes from ice discharge and that the land and atmosphere models will receive an updated topography from the ice sheet model. Therefore, topography-related feedbacks will affect the state of the GrIS and the climate. In the land model, the SMB is computed by subtracting the runoff and sublimation from the precipitation over the ice sheet. Precipitation is the sum of snow and rain. Snow directly contributes positively to the SMB, while rain can have either a net SMB contribution of zero if it runs off or a positive contribution if it refreezes. Melted snow and ice can also be divided into a fraction that contributes to the runoff and a fraction that refreezes. Therefore, additionally, we can express the SMB as follows:

1$$\begin{array}{rl}\text{SMB}& =\text{snow}+\text{refreezing}-\text{melt}\\ & -\text{sublimation}.\end{array}$$

Figure 1

Schematic of four elements of coupling between ice sheets and other Earth system components, courtesy of Michele Petrini . The land model (CLM5) receives the location of the ice sheet margin from the ice sheet model (CISM2). The land and the atmosphere (CAM6) models receive the dynamic GrIS topography from the ice sheet model. The ice sheet model provides the freshwater influx from ice discharge to the ocean model (POP), while the land model provides the freshwater fluxes from surface runoff. The ice sheet model itself receives the SMB from the land model to compute changes in the GrIS topography and extent.

[Figure omitted. See PDF]

The model computes the energy available for melt using the surface energy balance (SEB): 2$$\begin{array}{rl}\text{melt}& =\text{net shortwave}+\text{net longwave}\\ & +\text{latent heat}+\text{sensible heat}\\ & +\text{ground heat}.\end{array}$$

The SMB and SEB are computed for multiple elevation classes , which allows for resolving the large heterogeneity around the steep ice sheet margins and for the remapping of the SMB on the ice sheet model grid, which is done using trilinear interpolation and renormalization . For every elevation class, the grid cell temperature is downscaled using a uniform lapse rate of $-$6 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$. Precipitation falls as snow if the elevation-corrected near-surface temperature is below $-$2 °C and as rain if this temperature is above 0 °C. Between $-$2 and 0 °C, precipitation falls as a mix of rain and snow . Based on the downscaled temperature, precipitation and surface energy fluxes, the SMB is computed for every elevation class. Then, the SMB for a grid cell in the land model is the area-weighted average of these elevation classes. Finally, in the ice sheet model, the total GrIS mass balance (MB) is computed as the sum of the SMB downscaled from the elevation classes, basal mass balance and ice discharge. As the ice sheet loses or gains mass, the updated topography and ice sheet margins are communicated to the land model. Here, the land units and topographic height are updated according to the updated ice sheet margins and topography, redistributing the weights of the different elevation classes . The changing ice sheet topography influences not only the computations in the land model but also those in the atmosphere component as the topography influences atmospheric circulation. Therefore, the ice sheet topography computed in the ice sheet model is communicated to the atmosphere model, where it is remapped to the atmospheric grid in the form of surface geopotential height . Regarding the uni-directional coupling with the ocean, the ice sheet model computes the annual ice discharge and basal melting, which is communicated to the ocean model and then supplied to the ocean at a constant rate throughout the following year. Coupling from the ocean to the ice sheet model is not yet implemented. Therefore, ocean-forced melting of marine-terminating glaciers is not accounted for , which could lead to biases in the computed ice discharge. Runoff as a result of melt and rain, computed in the land model, is distributed to the ocean using the river transport model and is therefore seasonally variable.

In this study, we compare with a one-way coupled simulation to quantify elevation feedbacks. The one-way coupled simulation is a simulation in which changes in ice sheet elevation and extent are not communicated to the other model components. The GrIS topography is constant in the land and the atmosphere model, and the freshwater fluxes from the GrIS into the ocean model are fixed. The topography and freshwater fluxes in the one-way configuration are the same as the initial state of the two-way coupled simulation after the spin-up procedure . The resulting climate and SMB are downscaled from the fixed land and atmosphere model topography using elevation classes , with a temperature lapse rate of $-$6 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$, and are interpolated onto the ice sheet model grid. The coupler communicates the downscaled SMB to the ice sheet model, where a new ice sheet topography is computed. However, this is not communicated to the other model components when using one-way coupling. This implies that in the one-way setup, the downscaling procedure is used to account for both spatial heterogeneity within a grid cell and elevation differences between the ice sheet model and land and atmosphere models caused by ice sheet mass loss. Since the surface albedo is computed in the land model and the land model allows for the albedo to change, the surface albedo is updated in both the one-way and two-way coupled simulation. Therefore, when snow and ice melt or snow accumulates, the albedo is updated.

2.3 Simulation design

We compare a one-way and a two-way coupled simulation, which are both forced with a $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ scenario (Fig. a). In this scenario, we apply a yearly 1 $\mathrm{\%}$ increase in ${\mathrm{CO}}_{\mathrm{2}}$ concentration until 4 times the PI ${\mathrm{CO}}_{\mathrm{2}}$ concentration (1140 $\mathrm{ppm}$) is reached in year 140. After year 140, we keep the ${\mathrm{CO}}_{\mathrm{2}}$ concentration constant at $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ until year 500. The $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ scenario is an extreme warming scenario and, after the year 140, has a similar radiative forcing to that of the SSP5-85 scenario at the end of the 21st century . The one-way coupled simulation is an extension of the simulation done by . Next, we analyze another simulation branched from our two-way coupled simulation, in which we decrease the ${\mathrm{CO}}_{\mathrm{2}}$ concentration by an annual 5 $\mathrm{\%}$ between years 350 and 377 until PI ${\mathrm{CO}}_{\mathrm{2}}$ conditions are reached, after which the ${\mathrm{CO}}_{\mathrm{2}}$ concentration remains constant until year 925 (Fig. a). The climate and ice sheet exchange information every year for the first 500 $\mathrm{years}$ of the simulation. After year 500, the climate does not change rapidly anymore. Therefore, to save computational resources, the coupling is done every 5 $\mathrm{years}$ after the year 500. Next to the ${\mathrm{CO}}_{\mathrm{2}}$ forcing experiments, we use a 300-$\mathrm{year}$ two-way coupled PI control simulation with a constant PI ${\mathrm{CO}}_{\mathrm{2}}$ concentration for comparison.

Figure 2

Comparison of forcing, climate and ice sheet evolution of one-way (blue) and two-way (red) coupled simulations for (a) ${\mathrm{CO}}_{\mathrm{2}}$ forcing [$\mathrm{ppm}$], (b) global warming [$\mathrm{K}$], (c) NAMOC index [$\mathrm{Sv}$], (d) GrIS total mass balance (MB) [$\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$] and its components ice discharge (ID) and surface mass balance (SMB) (the basal mass balance is not displayed as its contribution is limited), (e) GrIS-induced cumulative sea level rise [$\mathrm{mm}$], (f) GrIS area [${\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{km}}^{\mathrm{2}}$], (g) GrIS equilibrium line altitude [$\mathrm{m}$] and (h) GrIS ablation area compared to total GrIS area [$\mathrm{\%}$]. The mass balance (components) in (d) can be directly converted to the GrIS-induced SLR rate [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$]. The dashed grey line indicates the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up period and the start of the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ period. If the differences between the one- and two-way coupled simulations become significant throughout the simulation, the first year of significant difference is marked with a black cross. In (e) the year of significance for the annual SLR is given. In (g) and (h) the year of significant difference before the convergence between the two simulations is given as well.

[Figure omitted. See PDF]

We consider 20 data points in time for centered moving averages and periodic means to obtain a climatological mean state and variability. This means that before the year 500, 20 $\mathrm{years}$ are considered, whereas after the year 500, 100 $\mathrm{years}$ are considered. As the climate evolves less rapidly after the year 500, a 100-$\mathrm{year}$ mean will be able to represent the mean state of the climate while being consistent in terms of variability compared to the period before the year 500.

2.4 Emergence and recovery

We use the emergence concept applied by to assess whether the differences between the one- and two-way simulations are significant. We define the first year of significant difference as the first year when the 20-$\mathrm{year}$ centered moving average of the differences between the one- and two-way simulation emerged from the natural variability of the differences. We define the natural variability of the differences as all absolute values smaller than 1 standard deviation of the differences (1${\mathit{\sigma }}_{\text{differences}}$). To obtain ${\mathit{\sigma }}_{\text{differences}}$, we use the standard deviations computed from the two-way coupled PI control simulation (${\mathit{\sigma }}_{\text{control}}$). We assume that this control simulation can represent the mean state of a similar one-way control simulation as the ice sheet is nearly in equilibrium, with a limited mean SLR rate of 0.03 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$, and its variance only represents natural variability. Then, using error propagation, we can express the 1${\mathit{\sigma }}_{\text{differences}}$ interval in terms of ${\mathit{\sigma }}_{\text{control}}$ as $[-{\mathit{\sigma }}_{\text{control}}\sqrt{\mathrm{2}}$, ${\mathit{\sigma }}_{\text{control}}\sqrt{\mathrm{2}}]$. When the 20-$\mathrm{year}$ centered moving average of the difference has migrated permanently outside the 1$\mathit{\sigma }$ interval, we consider the signal emerged.

To assess whether climate and ice sheet variables can recover to their PI state after the reintroduction of PI ${\mathrm{CO}}_{\mathrm{2}}$ levels, we apply a modified approach of the aforementioned emergence concept. A variable is considered recovered if the 20-$\mathrm{year}$ centered moving average returns within 1 standard deviation (1$\mathit{\sigma }$) of the PI mean and stays within the 1$\mathit{\sigma }$ interval (obtained from the PI control simulation) for the remainder of the simulation. In case a variable has not recovered by the end of the simulation, we compute an extrapolated year of recovery from the remaining linear trend over the last 100 $\mathrm{years}$ (825–925) and the remaining anomaly of the mean state of the variable in the last 100 $\mathrm{years}$ of the simulation with respect to the closest boundary of the 1$\mathit{\sigma }$ confidence interval. We compute the relative amount of recovery by comparing the remaining anomaly in the final 100 $\mathrm{years}$ of the simulation to the maximum response during the $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ forcing period.

3 Simulated mass loss in the one- and two-way coupled simulation

To assess the impact of simulating an interactive ice sheet, we first consider differences in the global and regional climate to see whether this affects or is affected by GrIS mass loss. As a response to the ${\mathrm{CO}}_{\mathrm{2}}$ forcing (Fig. a), global temperatures rise, with a global warming response of nearly 10 $\mathrm{K}$ by year 500, which is not affected by the coupling (Fig. b). Under the $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ forcing, the North Atlantic Meridional Overturning Circulation (NAMOC; the maximum strength of the overturning stream function north of 28° N and below 500 $\mathrm{m}$ depth) collapses in around 150 $\mathrm{years}$ in both simulations (Fig. c). A near-zero GrIS integrated mass balance in the first 120 $\mathrm{years}$ (Fig. d) results in limited mass loss in this period (Fig. e). During these years, the NAMOC evolution in the one- and two-way coupled simulations is similar, showing that this NAMOC collapse is not strongly related to increased freshwater fluxes from the GrIS.

After year 120, there is an accelerating decrease in SMB, which is the strongest in the one-way coupled simulation (Figs. d and g–j), leading to the separation of the SLR curves at the year 282 (Fig. e). In contrast, the contribution from ice discharge decreases in both simulations due to the melt–discharge feedback as the ice sheet starts retreating (Fig. b–e). Because of the small ice discharge contribution in both simulations, the total mass balance differences emerge in the same year as the SMB. After emergence, the one-way coupled simulation consistently projects $\mathrm{17.0}\pm \mathrm{0.4}$ $\mathrm{\%}$ more mass loss ($p<\mathrm{0.01}$). GrIS area loss acceleration is slightly delayed compared to the mass loss, and the differences become significant in year 300 (Fig. f). Ice mass is mainly lost at the margins, although less strongly in the two-way coupled simulation (Fig. b–e). The increased melting of the ice sheet is reflected in the rise in the equilibrium line altitude (ELA; Fig. g) and the expansion of the ablation area (the area with negative SMB; Fig. h). As melting is stronger in the one-way case, the ablation area expands faster and the ELA heightens faster. However, the ELA and percentage of the ice sheet that experiences ablation in the two simulations converge in years 320 and 400, respectively. As the GrIS loses mass, the equilibrium line moves more inland, but its altitude does not increase strongly anymore as the ice sheet's elevation decreases simultaneously. As the ice sheet extent becomes smaller, parts of the ablation area are lost, slowing down the increase in the total percentage of the ablation area. The period of this slow increase is reached later in the two-way simulation but results in a period (years 400–464) during which ablation and accumulation areas are similar in both simulations.

4 Climate feedback response to a dynamic GrIS topography

The differences in the mass balance evolutions of the two simulations indicate that the local climate and GrIS surface processes are affected by the coupling. We first consider the positive melt–elevation feedback, which enhances mass loss, and assess its representation in the one-way coupled configuration as we compute the climate using a fixed topography in the one-way simulation. Then, we analyze feedbacks related to atmospheric processes by considering atmospheric circulation, clouds and precipitation.

4.1 The melt-elevation feedback

The melt response to the ${\mathrm{CO}}_{\mathrm{2}}$ forcing causes a lowering of the topography (Fig. a–e). This triggers the melt–elevation feedback (Fig. ), as lower elevations coincide with higher temperatures, and consequently enhances the melt–albedo feedback, as bare ice and melting snow have a lower albedo. In the two-way simulation, we account for this by allowing the topography to evolve in the land model, while in the one-way simulation, the topography is constant. The effect of incorporating a dynamic topography on the melt becomes significant in year 215, and accounting for this leads to 66 $\mathrm{\%}$ more melt in the land model in the two-way simulation compared to the one-way simulation by year 500. Since both simulations include an interactive calculation of the albedo, the melt differences include both the melt–elevation and the resulting enhancement of the melt–albedo feedback. A quantification of the strength of this feedback has been done by . Next to the albedo, the melt is affected by atmospheric-related feedbacks as well.

Figure 3

Change in melt [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] as modeled by the land model (CLM5) for years 480–500 with respect to the control simulation over the ice sheet extent of the two-way simulation for the (a) one-way and (b) two-way coupled simulation. (c) difference between (a) and (b), showing the effect of topographic change on melt evolution. The dashed black lines depict contour lines for every 500 $\mathrm{m}$ elevation of (a) the one-way and (b–c) the two-way simulation. The light blue line represents the ELA of (a) the one-way and (b–c) the two-way simulation. All changes are significant according to the emergence criterion (Sect. ).

[Figure omitted. See PDF]

The 2 $\mathrm{m}$ air temperature ($T_{\mathrm{2}\mathrm{m}}$) increase as a response to elevation changes causes changes in the GrIS surface energy fluxes. The downscaling of the $T_{\mathrm{2}\mathrm{m}}$ and the downwelling longwave radiation (LW_in) has a large influence on the computed SEB . To assess the effect of the elevation change on these variables, we compute lapse rates resulting from elevation changes (Fig.  – monthly means, Fig.  – seasonal maps) and compare them with the applied lapse rates for elevation change in the one-way coupled simulation (Table ). These lapse rates are computed by dividing the difference in the mean state of a variable between the one- and two-way simulation by the difference in elevation between the one- and two-way simulation for the years 480–500. Thresholds of 250 $\mathrm{m}$ elevation difference and an ice sheet fraction of 90 $\mathrm{\%}$ of the grid cell are taken to exclude grid cells in which most of the changes are not caused by elevation changes. It should be noted that these lapse rates describe the change in a variable caused by GrIS surface changes rather than the rate of change over a changing pressure level as described by the free-atmospheric lapse rate.

Figure 4

Monthly mean lapse rates (dark red) for (a) 2 $\mathrm{m}$ air temperature [$\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$] and (b) incoming longwave radiation [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{km}}^{-\mathrm{1}}$] for years 480–500. The light blue line shows the monthly mean surface temperature [$\mathrm{K}$] of the locations considered for computing the lapse rates and is the same in both (a) and (b). The locations considered for the mean lapse rates and surface temperatures are shown in Fig. .

[Figure omitted. See PDF]

Table 1

2 $\mathrm{m}$ air temperature [$\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$] and downwelling longwave radiation [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{km}}^{-\mathrm{1}}$] lapse rates for annual, summer (June, July, August) and winter (December, January, February) means, compared to the applied lapse rates in the one-way coupled simulation, computed using the years 480–500. The number within brackets represents 1 standard deviation in the spatial domain.

Mean lapse rate	$T_{\mathrm{2}\mathrm{m}}$	LWin
Annual	$-$7.3 (2.4) $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$	$-$40.0 (19.4) $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{km}}^{-\mathrm{1}}$
Summer (JJA)	$-$2.4 (1.1) $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$	$-$22.1 (13.9) $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{km}}^{-\mathrm{1}}$
Winter (DJF)	$-$12.0 (4.1) $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$	$-$55.8 (26.4) $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{km}}^{-\mathrm{1}}$
Applied in one-way	$-$6 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$	0 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{km}}^{-\mathrm{1}}$

The one-way lapse rate of $-$6 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$ corresponds relatively well to the annual mean temperature lapse rate (Table ). However, applying these pre-defined lapse rates leads to an overestimation of the $T_{\mathrm{2}\mathrm{m}}$ in summer and an underestimation in winter. The overestimation in summer is partly compensated by not applying an LW_in lapse rate, resulting in an underestimation of LW_in in the one-way coupled simulation. In summer, the surface in the ablation area, which is the area in which most of the elevation change happens, is at the melting point (Fig. ). Therefore, part of the available energy will be used for melting instead of heating the atmosphere, leading to a limited $T_{\mathrm{2}\mathrm{m}}$ increase in the ablation area. The LW_in lapse rates show a similar pattern as these are largely influenced by atmospheric temperatures.

To assess the effect of the difference in applied lapse rates, we consider a point that transitions from the accumulation to the ablation area (66.44° N, 45° E, shown in Fig. ). Here, the one- and two-way simulations show a similar response as long as the elevation does not change (Fig. ). As soon as the mean summer surface temperature reaches the melting point (Fig. c, dashed orange line) and elevation starts to lower (Fig. a), $T_{\mathrm{2}\mathrm{m}}$ and melt responses start to diverge between the simulations. Since in the one-way simulation the temperature and surface fluxes correspond to a fixed elevation, melt does not increase as strongly as in the two-way case. However, adding the temperature change due to the lowering of the elevation by considering the $-$6 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$ lapse rate (dashed blue line in Fig. c) results in a large overestimation of the one-way $T_{\mathrm{2}\mathrm{m}}$, potentially causing an overestimation of the sensible and latent heat flux contribution to the melt energy. Although this is partially compensated by an underestimation of the incoming longwave flux, this could lead to an overestimation of the melt–elevation and melt–albedo feedback when computing mass loss in the ice sheet model, which could explain part of the overestimated mass loss.

Figure 5

Evolution of (a) elevation [$\mathrm{m}$], (b) annual melt [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] and (c) mean summer (June, July, August) 2 $\mathrm{m}$ air and surface temperature [$\mathrm{K}$] for the one-way (blue, 1 $\mathrm{W}$) and two-way (red, 2 $\mathrm{W}$) coupled simulations for a point that transitions to the ablation area (66.44° N, 45° E, shown in Fig. ). The dashed grey line indicates the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up period and the start of the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ period. The dashed orange line indicates the timing of reaching surface melt conditions throughout the whole summer (273 $\mathrm{K}$). In (c) the dashed blue line shows the evolution of one-way $T_{\mathrm{2}\mathrm{m}}$ when applying the $-$6 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$ lapse rate to the two-way elevation change.

[Figure omitted. See PDF]

4.2 Elevation feedbacks related to changes in precipitation, atmospheric circulation and clouds

Besides affecting temperature, changes in the GrIS topography can change the local climate in several other ways. First of all, the effect of the presence of the cold and high-elevation GrIS on atmospheric circulation might change as its surface lowers and warms, which might have large implications for the regional climate. Atmospheric blocking events, associated with persistent high-pressure conditions over Greenland, can be related to higher summer temperatures and reduced cloud cover over Greenland . Therefore, increases in the amount of summer atmospheric blocking have been linked to increased melt . We compute the Greenland blocking index (GBI) using the method proposed by (GB2) for both the summer (JJA) and winter (DJF). We subtract the area-weighted mean of the 500 $\mathrm{hPa}$ geopotential height over the Arctic region (60–80° N) from the Greenland region (60–80° N, 20–80° W) and normalize using the variability of the control simulation. We find a large reduction in the Greenland blocking index resulting from a changing GrIS surface (Fig. a and b), especially in summer. This aligns with the positive relationship between orography and blocking events . As the large changes in GrIS topography might have large influences on the atmospheric circulation around Greenland as a whole, our computed changes in the blocking index may not only consist of changes in blocking but are influenced by other changes in atmospheric circulation as well. Additionally, we find that the changing GrIS geometry has no significant effect on the phase of the North Atlantic Oscillation (NAO; the principle component corresponding to the leading empirical orthogonal function of the seasonal sea level pressure in the North Atlantic region (20–80° N, 90° W–40° E) normalized using the variability of the control simulation; ) (Fig. d and e). However, the summer NAO index shows a trend towards a positive phase as a response to the warming in both simulations ($p<\mathrm{0.001}$).

Figure 6

Evolution of (a, b) Greenland blocking index (GBI) [–] and (d, e) North Atlantic Oscillation (NAO) index [–] for (a, d) June–July–August (JJA) mean and (b, e) December–January–February (DJF) mean for the one-way (blue) and two-way (red) coupled simulations. The dashed grey line indicates the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up period and the start of the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ period. If the differences between the one- and two-way coupled simulations become significant throughout the simulation, the first year of significant difference is marked with a black cross. (c) shows the difference in SMB [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] caused by the coupling that can be explained by the difference in summer blocking index, computed using linear regression for the area with significant Pearson correlation ($p$ = 0.01) for the years 480–500, for grid points that have not deglaciated in both simulations. The light blue line represents the ELA in the two-way simulation in years 480–500.

[Figure omitted. See PDF]

To relate the decrease in blocking index to the projected mass loss, we regress the GBI differences between the one- and two-way coupled simulations onto the SMB differences (from CISM2; Fig. c). A positive regression coefficient indicates a linear relationship between the differences in the blocking index and the differences in SMB reduction, and its magnitude represents the slope of this linear relationship. The relationship between the GBI and SMB differences is the strongest in the ablation area. In the accumulation area, the surface melt is limited, which is likely the reason why we do not find a strong relationship with changes in GBI. As the topographies of the two simulations evolve differently, so does the SMB pattern, resulting in small areas in which the correlation signal is negative, indicating a slightly larger SMB reduction in the two-way coupled simulation. This might not be caused by the differences in blocking index but rather by differences in the topography. Considering areas with significantly strong correlation (positive or negative, as tested with a $t$ test using $p$ = 0.01), we can link 49 $\mathrm{\%}$ of the difference in SMB reduction between the two simulations at the final simulation years to the decrease in GBI (Fig. c). As the SMB differences are mainly dominated by melt differences in this period, we hypothesize that the decrease in blocking as a result of GrIS surface elevation changes (Fig. a–e) acts as a negative feedback on surface melting.

We find another negative feedback, this time related to clouds and water vapor. Under increasing temperatures, the atmosphere can contain more moisture (as described by the Clausius–Clapeyron relationship), which can lead to an increase in clouds and enhanced precipitation . The amount of precipitable water in the atmospheric column (Fig. d) increases in both simulations, although more strongly in the two-way coupled simulation, as the atmospheric column warms more and becomes larger as a result of elevation change. The increase in precipitable water translates to a stronger reflection and absorption of incoming solar radiation in the atmosphere (Fig. a). The differences in the atmospheric effect on incoming shortwave radiation and precipitable water both become significant in year 262, indicating they are strongly related. There is no increase in cloud cover (Fig. b) in either of the simulations or even a slight decrease in medium-level clouds. However, cloud cover itself does not strongly control the incoming shortwave radiation. Instead, the cloud optical thickness, influenced by the presence of liquid water and ice in the atmosphere, is an important determining factor for the incoming shortwave radiation . As the cloud water path (Fig. c) increases throughout the simulation, meaning clouds become thicker, more shortwave radiation is reflected. This effect is stronger than the effect caused by the slight decrease in cloud cover, resulting in a stronger cloud effect on incoming surface solar radiation (Fig. a) in the two-way coupled simulation. The differences in cloud water path and cloud effect on downwelling shortwave radiation both become significant in year 401, showing that the cloud thickness is the determining factor for incoming shortwave radiation rather than cloud cover.

Figure 7

Time series over the first-year ice sheet extent of one-way (blue, 1 $\mathrm{W}$) and two-way (red, 2 $\mathrm{W}$) simulations for (a) effect on downwelling shortwave radiation [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$] of the atmospheric column and clouds. The effect of the atmospheric column is computed by comparing the incoming solar radiation at the top of the atmosphere with the received radiation at the surface for cloud-free conditions. The effect of clouds is computed by comparing the incoming solar radiation at the top of the atmosphere with the received radiation at the surface for all-sky conditions and correcting for the effect of the atmospheric column. (b) cloud fraction [–] for low-, medium- and high-level clouds, (c) cloud water path [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$] and (d) precipitable water in the atmospheric column [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$]. The dashed grey line indicates the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up period and the start of the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ period. If the differences between the one- and two-way coupled simulations become significant throughout the simulation, the first year of significant difference is marked with a black cross.

[Figure omitted. See PDF]

Finally, a lowering topography leads to an overall increase in precipitation, especially in the southeast (Fig. a). Higher temperatures in the two-way simulation result in more evaporation from the surrounding oceans. An increase in the transport of this moisture leads to increased precipitation over the GrIS (Fig. d). In some locations in the western and northern margins, precipitation decreases, although not significantly in all locations, as orographic precipitation moves more land inward, following the margin. In the western margins, we see a corresponding strong increase in integrated vapor transport (IVT; specific humidity times the wind velocity, integrated from 1000 to 300 $\mathrm{hPa}$ and normalized by the gravity; ). The difference in precipitation is mainly caused by increased rainfall (Fig. c) rather than increased snowfall (Fig. b). Snowfall moves more towards the interior as temperatures in the margins become too high for precipitation to fall as snow. As a result, snowfall decreases in the ablation area and increases in the accumulation area (Fig. ). The increase in snowfall in the interior contributes positively to the SMB, acting as a negative feedback. However, the elevation effect on temperature results in a larger fraction of rainfall, which can enhance the melt–albedo feedback. Therefore, the changes in rainfall can be considered a positive feedback.

Figure 8

Annual mean differences (two-way minus one-way) resulting from the coupling for (a) total precipitation [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] (the sum of snowfall and rainfall) (b) snowfall [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$], (c) rainfall [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] and (d) mean integrated vapor transport (IVT) [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$] for years 480–500 over the first-year ice sheet extent. The dashed black lines depict contour lines for every 500 $\mathrm{m}$ elevation of the two-way simulation. The light blue line represents the ELA of the two-way simulation. Hatches denote areas in which the differences between the one- and two-way coupled simulations have not become significant before the year 500.

[Figure omitted. See PDF]

5 Ice sheet–climate interactions and mass balance response under reintroduction of PI ${\mathrm{CO}}_{\mathrm{2}}$ conditions

The $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ forcing triggers large responses in global, Arctic and GrIS temperatures (Figs. b and b). By year 350, the GrIS has a strongly negative mass balance in the two-way coupled simulation (Figs. d and d), causing an SLR rate of $\mathrm{6.6}\pm \mathrm{1.0}$ $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$, with a cumulative GrIS SLR contribution of 1.13 $\mathrm{m}$ (Figs. e and e). The discharge contribution to the mass balance is limited ($-\mathrm{78}\pm \mathrm{8.8}$ $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$) as the ice sheet has strongly retreated. Besides, the NAMOC has collapsed (Figs. c and c), reducing the amount of northward heat transport. We branch another simulation from the two-way coupled simulation, starting in year 350, and apply an annual 5 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ reduction until PI ${\mathrm{CO}}_{\mathrm{2}}$ conditions are reached in year 377 (Fig. a) to assess ice sheet–climate interactions and the resulting mass balance as a response to ${\mathrm{CO}}_{\mathrm{2}}$ reduction. We first consider global, Arctic and North Atlantic climate change and evaluate the GrIS mass loss response to ${\mathrm{CO}}_{\mathrm{2}}$ reduction. We zoom in to the first century after the ${\mathrm{CO}}_{\mathrm{2}}$ reduction, which is characterized by a complex transition phase towards a colder climate, and look at the NAMOC evolution and its impact on regional climate. Finally, we assess GrIS surface processes by considering the response of the snowpack.

Figure 9

Evolution of ${\mathrm{CO}}_{\mathrm{2}}$ reduction to pre-industrial ${\mathrm{CO}}_{\mathrm{2}}$ (green, reduc) and the full $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ (red, $\mathrm{4}\times$) simulations for (a) ${\mathrm{CO}}_{\mathrm{2}}$ forcing [$\mathrm{ppm}$], (b) global, evolving GrIS and Arctic near-surface air temperature ($T_{\mathrm{2}\mathrm{m}}$) anomalies with respect to PI [$\mathrm{K}$], (c) North Atlantic Meridional Overturning Circulation (NAMOC) index [$\mathrm{Sv}$], (d) GrIS total mass balance (MB) [$\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$] and its components ice discharge (ID) and surface mass balance (SMB) (the basal mass balance is not displayed as its contribution is limited) and (e) GrIS cumulative sea level rise contribution [$\mathrm{mm}$]. The mass balance (components) in (d) can be directly converted to the GrIS-induced SLR rate [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$]. The dashed grey lines indicate the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up, the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$, the ramp-down (RD) to PI ${\mathrm{CO}}_{\mathrm{2}}$ and the continuous PI ${\mathrm{CO}}_{\mathrm{2}}$ periods. If a variable has recovered throughout the simulation, the first year of recovery is marked with a black cross.

[Figure omitted. See PDF]

The response to ${\mathrm{CO}}_{\mathrm{2}}$ reduction can be divided into two periods. During the first period, spanning the 27 $\mathrm{years}$ during the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-down and the following 85 $\mathrm{years}$, the GrIS experiences a complex transitional phase because of strong interactions with the NAMOC. During this phase, the NAMOC is still weak as its recovery is delayed compared to the atmospheric response, giving rise to relatively cold temperatures in the Arctic and over the GrIS (Fig. b). Under these largely reduced air temperatures resulting from the ${\mathrm{CO}}_{\mathrm{2}}$ decrease and weak state of the NAMOC, the amount of GrIS melt decreases strongly, leading to a small positive SMB and MB (Fig. d). During the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-down and first years thereafter, there is an additional SLR of 1.18 $\mathrm{m}$. This is followed by a short period (years 413–461) of positive mass balance, leading to a 0.13 $\mathrm{m}$ sea level drop, after which GrIS SLR contribution is halted (Fig. d and e). The second period, spanning from the end of the transitional phase (year 462) to the end of the simulation, is characterized by initial temperature increases in the Arctic and over the GrIS as a response to an overshooting rebound of the NAMOC followed by a slow continuous decrease in atmospheric temperature and NAMOC strength. The mass balance becomes smaller than during the transitional phase, nearing zero, and recovers to its PI state ($\mathrm{0.04}\pm \mathrm{0.23}$ $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$) by the year 452 (Table ). However, from year 735 onward, the mass balance becomes slightly positive again, leading to an annual sea level drop of $\mathrm{0.06}\pm \mathrm{0.19}$ $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ at the end of the simulation. This indicates that there might be potential for ice sheet regrowth, although at this rate this would take over 10 000 years. At the end of the simulation, the SMB does not recover (Table ) as global and GrIS temperatures remain elevated compared to the initial state. However, the coinciding retreated ice sheet margins result in a smaller contribution of ice discharge to the mass balance compared to the first 2 centuries (Fig. d). Together with a small positive SMB, the small negative contribution of the ice discharge allows for a small positive total mass balance. We compare this state to the state of the GrIS under the same global $T_{\mathrm{2}\mathrm{m}}$ anomaly (2 $\mathrm{K}$) during the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up (year 70), where we see an ice sheet that is largely marine-terminating and is out of balance as a result of this temperature anomaly (Table , Fig. ). Although the ice sheet has a stronger positive SMB in year 70, the large ice discharge leads to net mass loss, while for the retreated ice sheet at the end of the ${\mathrm{CO}}_{\mathrm{2}}$ reduction simulation, the total mass balance is slightly positive despite the smaller SMB.

Table 2

Comparison of the state of the GrIS under an annual global $T_{\mathrm{2}\mathrm{m}}$ anomaly of 2.0 $\mathrm{K}$ during the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up period in year 70 (average over years 60–80) and at the end of the ${\mathrm{CO}}_{\mathrm{2}}$ reduction simulation (average over years 825–925). The number within brackets represents 1 standard deviation. A topographic representation of these states can be found in Fig. .

	Year 70	End of simulation
	(60–80 mean)	(825–925 mean)
Global $T_{\mathrm{2}\mathrm{m}}$ anomaly	2.0 (0.22) $\mathrm{K}$	2.0 (0.08) $\mathrm{K}$
GrIS $T_{\mathrm{2}\mathrm{m}}$ anomaly	1.9 (0.68) $\mathrm{K}$	2.9 (1.0) $\mathrm{K}$
Total mass balance	$-$105 (76) $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$	21 (70) $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$
SMB	438 (81) $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$	161 (70) $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$
Integrated SMB	224 (40) $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$	103 (43) $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$
Ice discharge	$-$520 (7.0) $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$	$-$124 (4.8) $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$

During the $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ forcing period, the North Atlantic warms less rapidly due to the weakening of the NAMOC, resulting in a warming hole southeast of Greenland (Fig. a). During the transitional ${\mathrm{CO}}_{\mathrm{2}}$ reduction phase, the North Atlantic cools rapidly as this warming hole persists under a weak NAMOC (Fig. c). However, the delayed NAMOC overshoot causes the North Atlantic temperatures to subsequently rise, causing the warming hole to disappear (Fig. e). Besides, the timing of changes in GrIS, North Atlantic and Arctic processes strongly coincides with the timing of NAMOC index changes (Fig. ), highlighting the complexity of the interactions between the ocean, atmosphere, Arctic sea ice and GrIS.

Figure 10

(a, c, e) Mean anomalies of near-surface temperature ($T_{\mathrm{2}\mathrm{m}}$) [$\mathrm{K}$] for the ${\mathrm{CO}}_{\mathrm{2}}$ reduction simulation compared to the control simulation for (a) years 280–300 during the constant $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ period, (c) years 380–400 just after the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-down in the transition period and (e) years 480–500 after the transition period. (b, d) Evolution of reduction to pre-industrial ${\mathrm{CO}}_{\mathrm{2}}$ (green) and the full $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ (red) simulations for (b) sea ice extent [${\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{km}}^{\mathrm{2}}$] and (d) sea ice thickness [$\mathrm{m}$] of the evolving sea ice extent in March and September. (f) Evolution of the mean mixed layer depth (MLD) of January, February and March (JFM) in the Labrador Sea (blue), Irminger Sea (red), Norwegian Sea (green) and Nansen Basin (pink) for the reduction to pre-industrial ${\mathrm{CO}}_{\mathrm{2}}$ simulation. In (b), (d), (f) the dashed grey lines indicate the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up, the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$, the ramp-down (RD) to PI ${\mathrm{CO}}_{\mathrm{2}}$ and the continuous PI ${\mathrm{CO}}_{\mathrm{2}}$ periods. If a variable has recovered throughout the simulation, the first year of recovery is marked with a black cross.

[Figure omitted. See PDF]

The Arctic interacts strongly with the NAMOC and the atmosphere. Due to a large drop in Arctic $T_{\mathrm{2}\mathrm{m}}$ as a response to reverting ${\mathrm{CO}}_{\mathrm{2}}$ conditions, the March Arctic sea ice extent (the area north of 60° N where sea ice concentration is greater than 15 $\mathrm{\%}$) recovers quickly (Fig. b), amplified by the thin ice feedback , as an area with thin sea ice expands faster. In contrast, the September sea ice does not recover since the newly formed sea ice is relatively thin (Fig. d), and ocean and atmosphere summer temperatures remain elevated. Besides, the sea ice melt is amplified by the sea ice–albedo feedback , leading to the nearly complete melting of the recently formed sea ice in summer. The regrowth of sea ice in winter influences the NAMOC strength as brine rejection leads to more saline high-latitude waters, enhancing overturning. The overshooting response of the NAMOC leads to a subsequent decrease in Arctic sea ice extent and thickness as a result of increased northward heat transport.

The mixed layer depth (MLD; Fig. f), which is an indicator of the amount of overturning and thus of the strength of the NAMOC, decreases in the Labrador Sea, Irminger Sea and Norwegian Sea during the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up period. In contrast, the MLD in the Nansen Basin starts to increase slightly. This indicates a movement of the area with strong overturning towards higher latitudes, which is possibly caused by the decrease in Arctic sea ice extent (Fig. b). The MLDs in the Labrador Sea, Irminger Sea and Norwegian Sea all experience an overshoot after the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-down, similarly to the NAMOC index. Only in the Labrador Sea, where we see the large remaining temperature anomalies as well (Fig. c) – indicating elevated heat transport, is there no complete recovery from the overshoot.

Looking at GrIS surface processes, we find that melt increases due to the ${\mathrm{CO}}_{\mathrm{2}}$ forcing are partly countered by enhanced refreezing (Fig. a). At the end of the ramp-up period, in year 130, the peak in the amount of refreezing is reached despite the ongoing increase in the amount of available water for refreezing. This implies that around this year, the pore space and/or available energy for melt starts to decrease. This peak refreezing behavior is typical of high-emission scenarios . The refreezing capacity (amount of refreezing divided by the amount of available water) peaks earlier as the snowpack is not able to refreeze a similar fraction of the larger amount of available water since it does not become thicker concurrently (Fig. a). After the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-down, the amount of refreezing slightly peaks again in year 478 as the snowpack partially recovers. However, the poorer state of the snow compared to its state before the $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ forcing was applied, characterized by higher snow temperatures and a thinner snowpack (Fig. ), prevents the snowpack from returning to refreezing rates similar to year 130 despite the similar amount of available water. Besides, the poorer state of the snowpack, combined with a larger amount of available water than in the initial state, results in a refreezing capacity that does not recover (Table ).

Figure 11

Evolution of the snow properties (a) refreezing capacity (black) [–], defined as the refreezing divided by the amount of available water for refreezing (meltwater and rain) and its components meltwater (dark green), rainfall (green) and refreezing (light green) [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] over the evolving ice sheet and (b) summer (JJA) albedo anomalies [–] with respect to the control simulation as a function of time and distance to the ice sheet margin for the reintroduction of PI ${\mathrm{CO}}_{\mathrm{2}}$ simulation and the location of the equilibrium line (black). In (b), the dark green line shows the relative recovery (in %) for each distance class. The dashed grey lines indicate the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up, the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$, the ramp-down (RD) to PI ${\mathrm{CO}}_{\mathrm{2}}$ and the continuous PI ${\mathrm{CO}}_{\mathrm{2}}$ periods. None of the variables in (a) have recovered by the end of the simulation.

[Figure omitted. See PDF]

Next to the refreezing, the ice sheet albedo does not completely recover either (Fig. b). In the cooler period around years 400–450, the albedo nearly recovers to its initial state, even around the margins, as melt strongly reduces. However, subsequent temperature increases due to the NAMOC recovery lead to more melt and rain and a corresponding decrease in albedo, aligning with the expansion of the ablation area, resulting in an ELA that does not recover. The lower surface albedo causes the net shortwave flux to remain larger, which is the largest contributor to the fact that the SMB and SEB (Fig. ) do not recover completely.

6 Discussion

The results presented in this study highlight the importance of accounting for interactions between the GrIS and the climate for SLR projections. We find that uncoupled simulations of ice sheet and climate evolution can lead to an overestimation of the SLR contribution. For CESM2, we find an overestimation of the GrIS SLR contribution of 17 $\mathrm{\%}$ under a $\mathrm{4}\times$PI ${\mathrm{CO}}_{\mathrm{2}}$ forcing when using the one-way coupled configuration. Similarly, found an overestimation of mass loss in an uncoupled AOGCM-ISM simulation, although these only become apparent after a total mass loss of around 2.5 $\mathrm{m}$ sea level equivalent (SLE) as opposed to 0.7 $\mathrm{m}$ SLE in our study. Our results are close to those of , who found a relative overestimation of mass loss of 13 $\mathrm{\%}$ in an uncoupled abrupt $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ simulation using the Glimmer ice sheet model coupled with FAMOUS-ice. attributed part of the smaller mass loss in the coupled simulation to a larger precipitation increase as precipitation in the southwest moves land inward, which is similar to what we find. However, the larger amount of rainfall does not contribute positively to the SMB, and, although the snowfall in the accumulation area is larger in the two-way coupled simulation, its contribution is limited (28 $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$ by year 500, Fig. ), compared to other feedbacks. Besides, related the smaller mass to a negative cloud feedback on the downwelling shortwave radiation due to a larger cloud fraction in the coupled simulation. Although we do not find significant differences in cloud fraction, we find significantly thicker clouds in the two-way coupled simulation. indicated that cloud thickness, rather than cloud cover, is the most important cloud-related control on incoming shortwave radiation as cloud cover is not a measure of transmissivity. This aligns with our results but is in contrast to .

This study is the first to look into the effect of GrIS topographic changes on atmospheric blocking. In our two-way coupled simulation, the blocking index is significantly reduced when the GrIS topography lowers. Besides, we find a strong relationship between the atmospheric blocking index and SMB evolution, linking nearly half of the difference in SMB reduction to the decrease in the blocking index. Therefore, we hypothesize that the topographic control on Greenland blocking occurrence can act as a negative feedback on surface melt. However, further investigation into the causes of these changes in blocking and their relationship with melt is necessary to make a definite attribution, especially since the robustness of the blocking index towards large topographic changes has not been evaluated. Observations and climate projections show that blocking strongly affects GrIS melting . However, and showed that climate models are not able to capture the present-day increase in GBI and consistently project a future decrease in blocking. Blocking projections should therefore be treated with caution. Nevertheless, we relate the significant decrease in GBI in our two-way simulation to large topographic changes, which are not the main driver of present-day GBI changes.

Besides, we hypothesize that a large part of the increased mass loss in the one-way simulation arises from the application of a uniform $T_{\mathrm{2}\mathrm{m}}$ lapse rate. Compared to the lapse rates computed from our one- and two-way simulations, the applied uniform lapse rate results in an overestimation of summer $T_{\mathrm{2}\mathrm{m}}$ in the ablation area, and, therefore, melt might be overestimated. Although the use of uniform lapse rates in offline simulations is common e.g.,, it has been pointed out that lapse rates are not temporally and spatially uniform and that the melt–elevation feedback is sensitive to the choice of the lapse rate . found that applying a seasonally and spatially varying lapse rate for downscaling SMB from one-way coupled CESM simulations of the MIS-11c Greenland ice sheet gives the most accurate representation of the melt–elevation feedback, showing that it is likely that offline corrections in ice sheet models can be improved by accounting for the temporal and spatial variability in the lapse rate in the ablation area by considering whether the surface temperature has reached melting point. It should be noted that our computed $T_{\mathrm{2}\mathrm{m}}$ lapse rates are affected by not only the surface reaching melting point but also the other negative feedbacks influencing temperature, such as increasing cloud thickness and reduced atmospheric blocking.

To investigate other interactions between ice sheet and climate, and their effect on the mass balance, we apply a 5 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ reduction to PI ${\mathrm{CO}}_{\mathrm{2}}$ and find that the GrIS SLR contribution can be halted. However, the chosen ${\mathrm{CO}}_{\mathrm{2}}$ forcing, the timing, and the rate of the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up and ramp-down likely have a large influence on the mass balance evolution and influence whether mass loss can be halted or reversed. In our simulation, a small positive mass balance during the transitional phase just after ${\mathrm{CO}}_{\mathrm{2}}$ reduction, as well as at the end of the simulation, allows for a small regrowth of the ice sheet. However, reversing the 1.1 $\mathrm{m}$ (equal to 14 $\mathrm{\%}$ of the initial ice sheet volume) SLE of mass lost during the forcing period would likely take tens of thousands of years. Our simulation is designed to investigate the response of the GrIS to ${\mathrm{CO}}_{\mathrm{2}}$ reduction and the climate interactions that play a role therein rather than finding thresholds for irreversible mass loss and potential equilibrium states in contrast to previous work. Despite the long timescales used in previous work, the importance of certain feedbacks and interactions has been pointed out before. Using a coupled AOGCM-ISM configuration, investigated the potential regrowth of the GrIS from 11 different initial states when reverting to PI climate. They showed that for states with an ice sheet area reduction of 10 $\mathrm{\%}$ and 20 $\mathrm{\%}$, which are on the same order as this study, climate–ice sheet interactions, resulting from the temperature dependence on elevation, make the difference in whether or not the ice sheet could regrow. Although in the study by the ice sheet is coupled with the ocean model, interactions with the NAMOC are likely not captured due to the experimental setup with asynchronous coupling. Using the FAMOUS-ice–GLIMMER coupled configuration, performed ${\mathrm{CO}}_{\mathrm{2}}$ removal experiments starting from multiple GrIS (steady) states and highlighted the importance of the melt–albedo feedback on reversibility, which our results agree with. Their mass loss trajectories do not show a complex transitional phase as the model setup does not allow for interaction with the ocean, meaning there is no influence from an overshooting NAMOC recovery. used two ice sheet models to explore the GrIS mass balance response to a linear temperature ramp-up and subsequent ramp-down. Their experiment with a 100-$\mathrm{year}$ ramp-up to 7 $\mathrm{K}$ warming followed by a 100-$\mathrm{year}$ ramp-down to a remaining 2 $\mathrm{K}$ warming and an additional simulated 100 $\mathrm{kyr}$ at this 2 $\mathrm{K}$ warming level is closest to our experiment. Depending on the model, they find a notably larger mass loss of 2–6 $\mathrm{m}$ SLE, although over a larger timescale. However, at the end of our simulation, the ice sheet has a small but positive mass balance, and therefore we do not expect the ice sheet to lose more mass if we were to extend the simulation. The study by only includes limited interactions between the ice sheet and the atmosphere (melt–elevation and melt–albedo feedback and precipitation changes) and no interactions with the ocean. These differences show the large impact of complex climate–ice sheet interactions on projected mass loss.

The large impact of interactions with the NAMOC in the transitional phase towards a colder climate is apparent in this study and is not discussed in other studies on GrIS mass loss behavior in ${\mathrm{CO}}_{\mathrm{2}}$ reduction scenarios. However, the response of the (N)AMOC to ${\mathrm{CO}}_{\mathrm{2}}$ reduction has been studied under similar 1 $\mathrm{\%}$ ramp-up and ramp-down $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ scenarios . There is agreement on the existence of overshoot behavior of the AMOC, reaching AMOC strengths above those in the PI climate, and this response is delayed compared to the atmospheric response to ${\mathrm{CO}}_{\mathrm{2}}$ reduction. This behavior is mainly attributed to the enhanced salt–advection feedback. Our study only focuses on interactions with the Arctic and North Atlantic. Hence, changes in salinity in the subtropics resulting from an intensification of the hydrological cycle are not taken into consideration in our analysis. Therefore, our assessment of NAMOC changes related to changes in the Arctic and North Atlantic describes interactions between and dependencies on the different processes in these areas rather than attributing causes and effects of NAMOC changes.

7 Conclusions

In this study, we compare the one- and two-way coupled configuration of CESM2-CISM2 to investigate GrIS elevation-related feedbacks and assess to what extent one-way coupled simulations can capture the elevation feedbacks. We find that in a $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ scenario, the sum of topography-related feedback responses results in 66 $\mathrm{\%}$ more melt by year 500 when considering a dynamic GrIS using two-way coupling. However, the offline topography correction in CISM2 in the one-way configuration overestimates these feedbacks, causing an additional mass loss of 17 $\mathrm{\%}$ in the one-way simulation. This overestimation partly arises from the overestimation of the positive melt–elevation feedback as well as from not including several negative atmospheric feedback mechanisms. In the one-way coupled simulation, a uniform temperature lapse rate of $-$6 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$ is used to correct for elevation changes. However, the two-way simulation indicates that the temperature lapse rate has a large seasonal variability, as the surface in the ablation area is at melting point in summer. Overestimated temperature lapse rates in summer cause an overestimation of the melt–elevation feedback in the one-way simulation. Taking this seasonal dependency into account when doing offline corrections for elevation could help resolve part of the discrepancy between the one- and two-way coupled simulations. For other models that heavily rely on temperature-related parameterizations in particular, it is of great importance to carefully consider the applied lapse rate as a uniform lapse rate might cause a large overestimation of melt. Besides enhancing melt, the changes in GrIS topography give rise to increased precipitation, reduced atmospheric blocking and decreased solar radiation reaching the surface arising from increased atmospheric water vapor content and cloud thickness. These atmospheric responses act as negative feedbacks on melt, resulting in smaller mass loss in the two-way coupled simulation. These findings stress the importance of considering ice sheet–climate interactions for SLR projections using two-way coupled simulations or a more realistic offline topography correction (e.g., spatially and temporally varying lapse rate).

A rapid reduction in atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ from $\mathrm{4}\times$PI to $\mathrm{1}\times$PI conditions leads to strong interactions between the global, GrIS and Arctic climate. Despite a remaining global temperature anomaly of 2 $\mathrm{K}$, compared to the initial PI climate, GrIS mass loss halts, following from the small contribution of discharge as the ice sheet has partially retreated under a total mass loss of 1.1 $\mathrm{m}$ SLE. The collapsed NAMOC induces a warming hole south of Greenland during the forcing period and the first years after the ${\mathrm{CO}}_{\mathrm{2}}$ ramp-down but eventually recovers and overshoots under the reintroduced PI ${\mathrm{CO}}_{\mathrm{2}}$, resulting in a cooling hole in this area. This gives rise to a complex transitional phase in which temperatures first decrease rapidly but subsequently increase. Arctic sea ice extent in winter recovers well, but the ice remains thin and largely melts away in summer. Consequently, Arctic temperatures remain high. Although the GrIS mass balance recovers, its SMB does not, resulting from elevated temperatures, elevated ocean heat transport, and snowpack that does not recover compared to the initial PI state. The refreezing capacity of the snowpack reduces as a result of a smaller snow depth and higher snow temperatures. Besides, the albedo of the snowpack remains lower, leading to a stronger melt–albedo feedback.

This idealized ${\mathrm{CO}}_{\mathrm{2}}$ ramp-down scenario shows the importance of ice sheet–climate interactions on the response of the GrIS mass balance and climate to ${\mathrm{CO}}_{\mathrm{2}}$ reduction. Overshoot scenarios are becoming more important in policy-making, while studies regarding GrIS mass loss reversibility focus on multi-millennial behavior and do not always account for interactions with the climate e.g., in. Although this study is a first step towards understanding GrIS mass balance evolution under ${\mathrm{CO}}_{\mathrm{2}}$ reduction and the ice sheet–climate interactions involved, considering shorter-term simulations and less rapid ramp-down scenarios is crucial for SLR projections. In short-term simulations, SLR might not be halted when decreasing ${\mathrm{CO}}_{\mathrm{2}}$ concentrations, as discharge contributes more to the mass balance, and the interactions with the NAMOC could be smaller, as the magnitude of the NAMOC response to ${\mathrm{CO}}_{\mathrm{2}}$ reduction depends on how large the ${\mathrm{CO}}_{\mathrm{2}}$ forcing is and for how long it is applied .

Appendix A One- and two-way coupled simulations

Figure A1

Ice thickness [$\mathrm{m}$] and surface mass balance (SMB) [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] as downscaled to the CISM2 grid (4 $\mathrm{km}$) for (a, f) the mean of the control simulation, (b, c, g, h) anomalies of the two-way coupled simulation with respect to the control simulation and (d, e, i, j) anomalies of the two-way with respect to the one-way coupled simulation for grid points that have not deglaciated in both simulations, for (b, d, g, i) the years 240–260 and (c, e, h, j) the years 480–500. The dashed black lines depict contour lines for every 500 $\mathrm{m}$ of elevation of the control (a, f) and two-way coupled (b–e, g–j) simulations. Hatches denote areas in which the differences between the one- and two-way coupled simulations have not become significant before the last year of the corresponding period. Because of a different evolution of ice sheet topography, the two simulations show a different SMB pattern, leading to a few locations in which the one-way simulation shows less SMB reduction.

[Figure omitted. See PDF]

Figure A2

Lapse rate of (a–c) near-surface (2 $\mathrm{m}$) air temperature [$\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$] and (d–f) downwelling shortwave radiation [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{km}}^{-\mathrm{1}}$] for (a, d) June–July–August (JJA) mean, (b, e) December–January–February (DJF) mean and (c, f) annual mean. The black lines depict contour lines for every 500 $\mathrm{m}$ elevation of the two-way simulation. The light blue line represents the ELA of the two-way simulation. The black cross shows the location of the point shown in Fig. .

[Figure omitted. See PDF]

Figure A3

Evolution of snowfall [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] in the one-way (blue) and two-way (red) coupled simulation for the end-of-simulation ablation and accumulation area. The dashed grey line indicates the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up period and the start of the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$ period. If the differences between the one- and two-way coupled simulations become significant throughout the simulation, the first year of significant difference is marked with a black cross.

[Figure omitted. See PDF]

Appendix B CO₂ reduction

Table B1

Year of recovery. In case the variable has not recovered by the end of the simulation (year 925), the remaining anomaly of the mean state of the variable in the last 100 $\mathrm{years}$ of the simulation with respect to the closest boundary of the 1$\mathit{\sigma }$ confidence interval is computed, as well as the remaining linear trend over the last 100 $\mathrm{years}$. A positive anomaly means that the end-of-simulation state of the variable is larger than in the control simulation. If the trend and remaining anomaly are in the opposite direction, an extrapolated recovery year is given. The relative amount of recovery [$\mathrm{\%}$] is given for variables that have not recovered by the end of the simulation.

	Year of	Remaining	Remaining	Extrapolated year	Percentage
	recovery	anomaly	trend	of recovery	recovered
Mean global $T_{\mathrm{2}\mathrm{m}}$	–	1.94 $\mathrm{K}$	$-$0.00022 $\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$	9567	76 $\mathrm{\%}$
Mean GrIS $T_{\mathrm{2}\mathrm{m}}$	–	2.30 $\mathrm{K}$	$-$0.0099 $\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$	1157	68 $\mathrm{\%}$
Mean Arctic $T_{\mathrm{2}\mathrm{m}}$	–	3.71 $\mathrm{K}$	0.0017 $\mathrm{K}{\mathrm{yr}}^{-\mathrm{1}}$	–	70 $\mathrm{\%}$
Total mass balance	452	–	–	–	–
Surface mass balance	–	$-$339 $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$	$-$0.23 $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{2}}$	–	89 $\mathrm{\%}$
Calving flux	–	$-$444 $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$	$-$0.088 $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{2}}$	–	16 $\mathrm{\%}$
March sea ice extent	–	$-\mathrm{1.04}\times {\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{km}}^{\mathrm{2}}$	$\mathrm{0.00076}\times {\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{km}}^{\mathrm{2}}$ ${\mathrm{yr}}^{-\mathrm{1}}$	2291	90 $\mathrm{\%}$
September sea ice extent	–	$-\mathrm{5.39}\times {\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{km}}^{\mathrm{2}}$	$-\mathrm{0.00092}\times {\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{km}}^{\mathrm{2}}$ ${\mathrm{yr}}^{-\mathrm{1}}$	–	14 $\mathrm{\%}$
March sea ice thickness	–	$-$0.64 $\mathrm{m}$	$-$0.00047 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$	–	55 $\mathrm{\%}$
September sea ice thickness	–	$-$0.74 $\mathrm{m}$	$-$0.000099 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$	–	9 $\mathrm{\%}$
NAMOC index	–	4.84 $\mathrm{Sv}$	$-$0.0062 $\mathrm{Sv}{\mathrm{yr}}^{-\mathrm{1}}$	1707	–
Labrador Sea MLD	–	115.5 $\mathrm{m}$	0.80 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$	–	–
Irminger Sea MLD	700	–	–	–	–
Norwegian Sea MLD	459	–	–	–	–
Nansen Basin MLD	915	–	–	–	–
Refreezing capacity	–	$-$0.08	0.000024 ${\mathrm{yr}}^{-\mathrm{1}}$	4352	62 $\mathrm{\%}$

Figure B1

Ice thickness [$\mathrm{m}$] downscaled to the CISM2 grid (4 $\mathrm{km}$) for (a) the mean of the control simulation, (b) the anomalies of the ${\mathrm{CO}}_{\mathrm{2}}$ reduction simulation with respect to the control simulation in year 70 (60–80 mean) and (c) the same as in (b) for the end of the simulation (year 825–925).

[Figure omitted. See PDF]

Figure B2

Evolution of (a) GrIS annual mean SMB [$\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$] and its components and (b) GrIS JJA (June, July, August) mean SEB [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$] and its components for the simulation in which we reintroduce PI ${\mathrm{CO}}_{\mathrm{2}}$. The dashed grey lines indicate the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up, the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$, the ramp-down (RD) to PI ${\mathrm{CO}}_{\mathrm{2}}$ and the continuous PI ${\mathrm{CO}}_{\mathrm{2}}$ periods.

[Figure omitted. See PDF]

Figure B3

The 20-$\mathrm{year}$ moving average of the derivative of processes interacting with the NAMOC (black). These are the winter Arctic $T_{\mathrm{2}\mathrm{m}}$; the winter Arctic surface heat flux; the GrIS freshwater flux; the March and September sea ice extent; and the winter MLDs in the Labrador Sea, Irminger Sea, Norwegian Sea and Nansen Basin. The derivatives are computed using the central difference method and thereafter scaled to their maximum response. The dashed grey lines indicate the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up, the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$, the ramp-down (RD) to PI ${\mathrm{CO}}_{\mathrm{2}}$ and the continuous PI ${\mathrm{CO}}_{\mathrm{2}}$ periods.

[Figure omitted. See PDF]

Figure B4

Anomalies of (a) annual mean snow depth [$\mathrm{m}\mathrm{w}.\mathrm{e}.$] and (b) annual mean snow temperature [$\mathrm{K}$] over the GrIS with respect to the control mean as a function of time and distance to the ice sheet margin for the reintroduction of PI ${\mathrm{CO}}_{\mathrm{2}}$ simulation. The black line denotes the mean distance of the equilibrium line to the margin, and the dark green line shows the relative recovery (in %) for each distance class. The dashed grey lines indicate the end of the 1 $\mathrm{\%}$ ${\mathrm{CO}}_{\mathrm{2}}$ ramp-up, the continuous $\mathrm{4}\times$${\mathrm{CO}}_{\mathrm{2}}$, the ramp-down (RD) to PI ${\mathrm{CO}}_{\mathrm{2}}$ and the continuous PI ${\mathrm{CO}}_{\mathrm{2}}$ periods.

[Figure omitted. See PDF]

Code availability

CESM2-CISM2 is a community-developed open-source model available at https://www.cesm.ucar.edu/models/cesm2 and https://github.com/ESCOMP/CESM/tree/release-cesm2.1.5 .

Data availability

Simulation results are stored on the Cheyenne supercomputer and available on request.

Author contributions

TF designed and performed the analysis and wrote the manuscript; MV and BW supervised this research; MV, RS, MP and KTC designed and performed the simulations; and all authors read and provided inputs to the manuscript.

Competing interests

At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

Leo van Kampenhout is thanked for providing the code to prescribe ice sheet freshwater fluxes for the one-way simulation. Sotiria Georgiou is thanked for implementing prescribed freshwater fluxes in the one-way simulation and providing tools to analyze ocean model output. This material is based upon work supported by the National Center for Atmospheric Research (NCAR), which is a major facility sponsored by the National Science Foundation (NSF) under cooperative agreement no. 1852977. Computing and data storage resources, including the Cheyenne supercomputer (10.5065/D6RX99HX, ), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR.

Financial support

Miren Vizcaino, Michele Petrini and Raymond Sellevold received funding from the European Research Council EU H2020 program (grant no. RC-StG40 678145-CoupledIceClim). The CESM project is supported primarily by the NSF.

Review statement

This paper was edited by Ruth Mottram and reviewed by Jorge Bernales and two anonymous referees.