Accelerated Greenland Ice Sheet Mass Loss Under High Greenhouse Gas Forcing as Simulated by the Coupled CESM2.1‐CISM2.1

Abstract

The Greenland ice sheet (GrIS) is now losing mass at a rate of 0.7 mm of sea level rise (SLR) per year. Here we explore future GrIS evolution and interactions with global and regional climate under high greenhouse gas forcing with the Community Earth System Model version 2.1 (CESM2.1), which includes an interactive ice sheet component (the Community Ice Sheet Model v2.1 [CISM2.1]) and an advanced energy balance‐based calculation of surface melt. We run an idealized 350‐year scenario in which atmospheric CO₂ concentration increases by 1% annually until reaching four times pre‐industrial values at year 140, after which it is held fixed. The global mean temperature increases by 5.2 and 8.5 K by years 131–150 and 331–350, respectively. The projected GrIS contribution to global mean SLR is 107 mm by year 150 and 1,140 mm by year 350. The rate of SLR increases from 2 mm yr⁻¹ at year 150 to almost 7 mm yr⁻¹ by year 350. The accelerated mass loss is caused by rapidly increasing surface melt as the ablation area expands, with associated albedo feedback and increased sensible and latent heat fluxes. This acceleration occurs for a global warming of approximately 4.2 K with respect to pre‐industrial and is in part explained by the quasi‐parabolic shape of the ice sheet, which favors rapid expansion of the ablation area as it approaches the interior “plateau.”

Introduction

The Greenland ice sheet (GrIS) is the largest freshwater reservoir in the Northern Hemisphere, storing around 7.4‐m potential global mean sea level rise (SLR) (Bamber et al., 2013; Morlighem et al., 2017). Between 2007 and 2017, the GrIS has contributed to global mean SLR at a rate of 0.7 mm yr⁻¹ (Shepherd et al., 2019), as a result of increased surface melt, runoff, and ice discharge to the ocean (van den Broeke et al., 2016). The GrIS contribution to global mean SLR is expected to increase, with large uncertainty ranges (Bamber et al., 2019; Church et al., 2013; Oppenheimer et al., 2019).

The contemporary Greenland and Antarctic ice sheets are integral components of the Earth system. They are sensitive to climate change, and they influence climate through changes in topography, albedo, and freshwater fluxes to the ocean (Fyke et al., 2018). Among the major ice sheet/atmosphere feedbacks and interactions are the elevation feedback on melt (Edwards et al., 2014; Oerlemans, 1981), the ice/albedo feedback (Box et al., 2012), the coupling between the surface mass balance (SMB) and ice discharge (Goelzer et al., 2013; Lipscomb et al., 2013), and effects of orographic change on atmospheric circulation (Ridley et al., 2005).

Earth system models (ESMs) coupled to ice sheet models are important tools to better understand how these interactions affect GrIS mass loss. There has been much recent progress in improving these tools (Goelzer et al., 2017; Hanna et al., 2020; Pollard, 2010; Rybak et al., 2018; Vizcaino, 2014). Vizcaíno et al. (2013) and Alexander et al. (2019) highlighted the importance of computing SMB based on the surface energy balance; melt parameterizations based on temperature cannot simulate climate feedbacks and interactions in a physically realistic way. The realism of the SMB calculation in an ESM depends on the snow physics scheme (van Kampenhout et al., 2017), the albedo calculation (Helsen et al., 2017), and model resolution (Gregory & Huybrechts, 2006; Lofverstrom & Liakka, 2018; van Kampenhout et al., 2019), among other factors. The computational demand of coupling large‐scale climate processes with local‐scale ice sheet processes, combined with the long response time of ice sheets, is an additional challenge, as is the conservation of mass and energy in coupled climate/ice sheet models (Vizcaino, 2014).

Recognizing these challenges, the Community Ice Sheet Model version 2.1 (CISM2.1) has recently been included as an interactive component in the Community Earth System Model 2.1 (CESM2.1). The model includes bidirectional coupling of the ice sheet with the land and atmosphere through an energy‐based calculation of surface melt, downscaled through elevation classes to the ice sheet model grid (Lipscomb et al., 2013; Sellevold et al., 2019), along with dynamic ice sheet topography and glacier cover (Muntjewerf et al., 2020). In addition, the GrIS provides time‐evolving freshwater fluxes to the ocean component, although ocean forcing of marine‐terminating glaciers is not included in this study.

This paper presents the results of a multicentury (350‐year) simulation of GrIS evolution in an idealized scenario. Atmospheric CO₂ increases by 1% annually to four times the preindustrial value and then is held at this high value. Our goal is to investigate multicentury GrIS mass loss with a focus on ice sheet/climate interaction and feedbacks. We assess the timing and magnitude of the GrIS response in relation to global and regional climate change, and the relative importance of the processes that govern GrIS behavior. We also investigate nonlinearity in the response to the forcing and apparent accelerations and tipping points in the GrIS contribution to SLR. In a companion study, CESM2.1‐CISM2.1 was used to explore the GrIS response under historical and SSP5‐8.5 forcing from 1850 to 2100 (Muntjewerf et al., 2020).

Section 2 describes the coupled model CESM2.1‐CISM2.1 and the experimental setup. The main results are analyzed in section 3, which includes subsections on (1) overall changes in ice sheet mass and climate,(2) global, Arctic, and North Atlantic climate change and evolution of ocean circulation in relation to the timing of Greenland freshwater forcing, (3) expansion of GrIS ablation areas and changes in SMB components, (4) changes in the surface energy budget and summer Greenland climate, and (5) changes in ice sheet dynamics. Section 4 contains a discussion and conclusions.

Method: Model Description and Experimental Setup

Model Description

Model simulations were carried out with CESM2.1 (Danabasoglu et al., 2020), a fully coupled, global ESM with prognostic components of the atmosphere, ocean, land, sea ice, and land ice. CESM2 is one of the models contributing to the Coupled Model Intercomparison Project phase 6 (CMIP6; Eyring et al., 2016) and the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6; Nowicki et al., 2016). The atmosphere is simulated with the Community Atmosphere Model version 6, using the finite volume dynamical core (CAM6‐FV; Lin & Rood, 1997), on a nominal 1° horizontal grid with 32 vertical levels. The ocean component is the Parallel Ocean Program version 2 (POP2; Smith et al., 2010), which runs on a nominal 1° displaced‐pole grid with 60 vertical levels. Sea ice is represented by the CICE5 sea‐ice model (Hunke et al., 2015), which runs on the same horizontal grid as POP2. Land processes are simulated by the Community Land Model version 5 (CLM5; Lawrence et al., 2019), which runs on the same horizontal grid as CAM6. The CLM5 snowpack includes up to 10 vertical layers with a maximum total depth of 10‐m water equivalent. CLM5 also includes the Model for Scale Adaptive River Transport (MOSART) to handle land surface runoff based on topographic gradients.

The GrIS is simulated using the CISM2.1 (Lipscomb et al., 2019), using a 4‐km rectangular grid with 11 terrain‐following vertical levels. The velocity solver uses a depth‐integrated higher‐order approximation (Goldberg, 2011) of the Stokes equations for ice flow. A pseudo‐plastic sliding law and simple basal hydrology model as described by Aschwanden et al. (2016) is used to parameterize basal sliding. In this scheme, the yield stress is determined by the till friction angle and the effective pressure. The former depends on the bedrock elevation via a fixed piecewise linear relationship. The bedrock evolves according to an Elastic plate Lithosphere plus Relaxing Asthenosphere (ELRA) model. Calving in this study is parameterized with a flotation criterion, whereby floating ice is discharged instantaneously to the ocean.

Coupling Description

In the default CESM2 configuration, ice sheets do not evolve, but the simulations described here with CESM2.1‐CISM2.1 have a dynamic GrIS that is interactively coupled to other Earth system components. The model features an SMB calculation with a surface energy balance calculation of melt. The SMB (the sum of snowfall and refreezing minus sublimation and melt) is computed in CLM5 in multiple elevation classes (ECs). In this study (and by default), there are up to 10 ECs with fixed elevation ranges per glaciated grid cell (Lipscomb et al., 2013; Sellevold et al., 2019). The EC downscaling uses a fixed near‐surface temperature lapse rate of 6 K km⁻¹, and specific humidity is downscaled assuming constant relative humidity with elevation. Rainfall is converted to snowfall when the near‐surface temperature is lower than ‐2°C and snowfall to rainfall when the near‐surface temperature is higher than 0°C. In the case of snowfall‐to‐rainfall conversion, the resulting rainfall is routed directly to the ocean. The EC scheme features interactive coupling to CAM6 and explicit modeling of albedo, refreezing, and snow and firn compaction (van Kampenhout et al., 2017, 2020). The SMB is then downscaled by the coupler to the higher‐resolution CISM grid using a trilinear remapping scheme (bilinear in the horizontal and linear in the vertical), with corrections of accumulation and ablation to conserve global water mass.

The freshwater flux received by POP2 from the GrIS is the sum of surface runoff from CLM5 and basal melt and ice discharge from CISM2.1. These fluxes are supplied as salinity anomalies. Surface runoff is routed to the ocean via MOSART based on topographic gradients. These gradients are static; that is, they are not updated during the simulations to account for changes in the ice sheet topography. In the ocean, this flux together with basal melt from CISM2.1 is distributed by an estuary box model over the upper 30 m of the grid cell (Sun et al., 2017). Ice discharge as calculated by CISM2.1 is delivered to the nearest ocean grid cell and spread horizontally in the surface layer with a Gaussian distribution over a distance of 300 km. This ice is melted instantaneously, with the corresponding latent heat removed from the ocean.

In CLM5, land surface types (known as land units) change dynamically from glaciated to vegetated as the ice sheet retreats and vice versa when the ice sheet advances. The ice sheet surface topography from CISM2.1 is used to recompute the fractional glacier coverage in each CLM5 grid cell, modifying the albedo and soil and vegetation characteristics. The evolving ice sheet topography is coupled offline, at regular intervals, to the atmosphere model, which enables orographic circulation feedbacks. For this study, mean surface elevation and subgrid‐scale surface roughness fields of CAM6 are updated every 10 model years based on CISM2.1 topography.

Experimental Setup

We analyze two simulations: a 300‐year control simulation of the pre‐industrial era (year 1850 CE Danabasoglu, 2019a) and a 350‐year transient simulation with an idealized CO₂ scenario (Danabasoglu, 2019b). The atmospheric CO₂ concentration increases by 1% per year until reaching 1,140 ppmv, four times the pre‐industrial value (hereafter 4×CO₂) in year 140 (Figure 1a). The 4×CO₂ level is maintained for the remaining 210 years. The purpose of this idealized scenario is to distinguish direct responses of climate and ice sheet from lagged responses and feedbacks, while eliminating the effects of aerosol and land use changes. These simulations are part of ISMIP6. Further details on the forcing scenarios are provided by Eyring et al. (2016), and details on the experimental setup are provided by Nowicki et al. (2016).

1 Figure. Evolution of (a) CO2 (ppm), (b) top‐of‐atmosphere (TOA) radiation imbalance (W m−2), (c) global, GrIS, and Arctic region (north of 60°N) annual near‐surface temperature anomaly with respect to pre‐industrial mean, (d) cumulative GrIS contribution to global mean sea level rise (mm), (e) mass balance (MB, black) and components (SMB, orange; ice discharge [ID], green; basal mass balance [BMB], blue) in mm yr−1 (left) and Gt yr−1 (right), and (f) GrIS total area (left vertical axis, million km2) and percentage with respect to initial area (right vertical axis). Ice discharge and BMB are defined as negative in (e) for visualization purposes.

Both simulations start from the spun‐up pre‐industrial Earth system/ice sheet state described by Lofverstrom et al. (2020). A near‐equilibrium state was obtained by alternating between a fully coupled configuration with all model components active and a more efficient partially coupled configuration with a data atmosphere component. After almost 10,000 years of ice sheet spin‐up, the residual drift in the near‐equilibrated GrIS volume is 0.03 mm yr⁻¹, with GrIS volume and area overestimated by 12% and 15%, respectively, compared to present‐day observations. The simulated pre‐industrial ice sheet velocities and SMB compare reasonably well with present‐day observations and regional modeling reconstructions. van Kampenhout et al. (2020) evaluated the SMB and Greenland climate in CESM2 simulations with GrIS topography fixed at observed present‐day values, and Muntjewerf et al. (2020) evaluated the simulated present‐day ice discharge.

Results

The analysis focuses on two main periods, broadly representative of one‐century and multicentury change in the GrIS. We will refer to years 131–150 (around the time the model reaches 4×CO₂) as “CO₂ stabilization” and years 331–350 as “end‐of‐simulation” (indicated by blue shading in Figure 1). The CO₂ stabilization period has similar radiative forcing to the CMIP RCP8.5 and SSP5‐85 scenarios at the end of the 21st century (as in Muntjewerf et al., 2020; Vizcaino et al., 2015). The end‐of‐simulation period illustrates the multicentury GrIS response to elevated greenhouse gas forcing comparable to year 2300 under high scenario forcing (e.g., Aschwanden et al., 2019; Vizcaino et al., 2015).

GrIS Contribution to SLR

Figure 1b shows the evolution of the cumulative top‐of‐atmosphere (TOA) radiation imbalance. The TOA imbalance increases until CO₂ stabilization and decreases thereafter. In response, the global annual average near‐surface temperature increases at an approximately constant rate in the first 140 model years. By CO₂ stabilization, the warming is 5.2 K (σ = 0.3 K; Figure 1c), and the long‐term trend is about 0.5 K per decade. In the next two centuries, the temperature increases by an additional 3.3 K as the oceans continue to warm. By the end of simulation, the global warming trend has decreased to around 0.1 K per decade. However, a weak warming trend is expected to continue long after CO₂ stabilization, as it can take several millennia for the Earth system to reach a new equilibrium state (Rugenstein et al., 2020). Temperatures in the Arctic (defined as north of 60°N) follow a broadly similar trajectory. Polar amplification (the ratio of Arctic and global temperature increase) is 1.6, with much of the signal coming from Arctic sea ice loss. GrIS amplification (the ratio of GrIS and global temperature increase) is about 1.2; it is smaller than the Arctic amplification because the GrIS is a terrestrial region with a perennial ice/snow cover that holds the summer surface temperature near the melting point.

The simulated pre‐industrial GrIS is close to equilibrium, with a global mean SLR contribution of 0.03 mm yr⁻¹ and a relatively large standard deviation of 0.23 mm yr⁻¹ over 300 simulation years (Table 1). The trajectory of mass loss in the 1% simulation can be separated into three distinct periods (Tables 1 and S1 and Figure 1e). In the first period (years 1–119), the GrIS responds slowly and the mass loss increases by 2.4 Gt yr⁻². The modern observed mass loss of 250 Gt yr⁻² (∼0.7 mm yr⁻¹, Shepherd et al., 2019) is not reached until early in the second century. Then, from year 120 at a global warming of 4.2 K, the mass loss increases by 11.3 Gt yr⁻² until year 225. The average SLR contribution in years 131–150 is 764 Gt yr⁻¹ (+2.2 mm yr⁻¹). Finally, during years 226–350, the mass loss increases at a slower rate of 4.6 Gt yr⁻². The average mass loss by the end of the simulation is 2,350 Gt yr⁻¹ (+6.6 mm  yr⁻¹). The cumulative contribution to SLR is 107 mm by year 150 and 1,140 mm by year 350 (Figure 1d). The major GrIS mass budget components, SMB and dynamic ice discharge, are discussed in sections 3.3 and 3.5. The basal melt rate is not discussed further, since its contribution to the mass budget is small (between 14 and 24 Gt yr⁻¹, Table 1).

1 TableAnnual Rate of GrIS Mass Loss (mm yr⁻¹), Cumulative GrIS Mass Loss (mm), Mass Balance Components (Gt yr⁻¹), and GrIS Area (10⁶ km²)

	Pre‐industrial	Years 131–150	Years 331–350
Annual mass loss	0.03 (0.23)	2.16 (0.47)	6.58 (1.04)
Cumulative mass loss	11	107	1,140
MB	–13 (84)	–764 (160)	–2,350 (358)
SMB	585 (85)	–367 (166)	–2,259 (357)
ID	574 (5)	378 (26)	77 (8)
BMB	–24 (0)	–19 (4)	–14 (0)
GrIS area	1.966	1.918	1.598

Note. Mass Balance (MB) = Surface Mass Balance (MB) ‐ Ice Discharge (ID) + Basal Melt Balance (BMB). Three time periods are shown, corresponding to pre‐industrial (years 1–300) and 1% simulations (years 131–150 and 331–350). Values are time average (with standard deviation between parentheses if applicable).

Global, Arctic, and North Atlantic Climate Change

Global and Regional Climate Change

The annual mean near‐surface temperature increases globally (Figure S1), with the greatest warming (>18 K by end of simulation) in the Arctic basin, the Canadian archipelago, and Antarctica. The North Atlantic warms the least, in connection with changes in the ocean circulation and associated meridional heat transport (see section 3.2.2). The Arctic becomes seasonally ice free by the end of the first century (Figure S2) and almost completely ice free by year 270, with March sea ice extent declining to less than 2 × 10⁶ km².

Zonal averages of June–August (JJA) and December–February (DJF) near‐surface temperatures are shown in Figure 2. Before year 140, the high Arctic (>80°N) warms less than lower Northern Hemisphere latitudes. This is possibly connected with the decreasing sea ice cover. From years 140 to 350, the high Arctic warms more than other Northern Hemisphere latitudes due to the lack of sea ice and a generally reduced (seasonal) snow cover. By end of simulation, the interior of the GrIS is the only region in the Northern Hemisphere where near‐surface temperatures remain below freezing throughout the summer months (Figure 2b).

2 Figure. Zonal mean (top) and maps (bottom) of summer (JJA; left) and winter (DJF; right) near‐surface temperature (°C). The lines show the pre‐industrial control (CTRL; black), 4×CO2 years 131–150 (blue), and 4×CO2 years 331–350 (red). The maps show the seasonal averages at end of simulation (years 331–350).

The zonally averaged near‐surface temperature in the Northern Hemisphere winter (Figure 2c) highlights the polar amplification for years 131–150 and 331–350. The meridional temperature gradient reverses at ∼70°N in both periods, though more pronounced in the latter period. This reversal reflects the sea ice thinning and retreat by 131–150, and sea ice‐free conditions by 331–350, given the Arctic land‐ocean distribution. By end of simulation, the surface temperature is below freezing in most Arctic land areas in boreal winter; however, the Arctic Ocean remains temperate through the winter months and thus stays largely ice free throughout the year (Figures 2d and S2). The GrIS is the coldest region in the Northern Hemisphere; this climatic signature of the GrIS by 331–350 is illustrated in Figure 2d.

Changes in Ocean Circulation and Relationship with Greenland Meltwater

The North Atlantic Meridional Overturning Circulation (NAMOC) weakens significantly during the period of CO₂ increase (Figure 3a). The NAMOC index—defined as the maximum of the overturning stream function north of 28°N and below 500‐m depth—decreases at a rate of about 0.12 Sv yr⁻¹ until year 140 and by 0.06 Sv yr⁻¹ from 140 to 170, reaching values below 6 Sv. The thickness of the upper, poleward‐moving branch of the overturning cell decreases by 1 km by CO₂ stabilization (Figure 3b). By end of simulation, the location of the maximum overturning has migrated equatorward and towards the surface. Significant NAMOC reduction and shoaling of the overturning cell in elevated CO₂ scenarios are common in CMIP6 ocean models, which generally simulate greater NAMOC decline than CMIP5 models (Weijer et al., 2020).

3 Figure. (a) Evolution of mean January‐February‐March mixed layer depth at the local maximum within four regions: Labrador Sea (red), Irminger Sea (purple), Iceland Basin (green), and Barent Sea (blue), and the NAMOC Index (Sv, black). Dashed lines represent pre‐industrial. (b) Mean North Atlantic Meridional stream function (Sv) for (upper panel ) pre‐industrial, (middle panel ) years 131–150, and (bottom panel) end of simulation (years 331–350).

Figure 3a shows the simulated evolution of the 20‐year mean January–March mixed layer depth (MLD) in the deep convection regions in the Labrador Sea, Irminger Sea, Iceland Basin, and Barents Sea. The mixed layer in all regions becomes much shallower in the first 70 to 100 years, due to enhanced buoyancy from surface warming and additional freshwater fluxes from sea ice reductions, enhanced GrIS runoff, and enhanced Arctic precipitation. The Labrador Sea is the first region where the mixed layer shallows to a minimum depth indicative of negligible deep convection. Next, deep convection ceases in the Irminger Sea, Iceland Basin, and Barents Sea. By year 150, convection has stopped in all four regions.

Next, we compare NAMOC evolution with the evolution of the freshwater flux from the GrIS (Figure 4). In this way, we tentatively explore a causal relationship between the simulated NAMOC decline and the accelerated melt over the GrIS, in the absence of a conclusive paired one‐way coupled simulation that isolates the role of the bidirectional coupling (as in e.g., Mikolajewicz et al., 2007). As such, we consider the freshwater fluxes from both GrIS surface runoff and ice discharge. The ice discharge (solid freshwater flux) decreases during the simulation, as shown in Figure 1e and examined in detail in section 3.5. From year 110, the liquid freshwater flux (from surface runoff and basal melt) accelerates, with a similar timing of the net mass loss acceleration discussed in section 3.1. At this time, the magnitude of the runoff is less than 0.3 × 10⁵ m³ s⁻¹, but the average NAMOC strength has already declined to below 10 Sv. The NAMOC decline begins before year 50, at a time where the GrIS melt signal has not yet emerged from background variability. By end of simulation, runoff reaches values more than 0.8 × 10⁵ m³ s⁻¹, and the NAMOC index has remained near 5 Sv for almost two centuries. The spatial map of runoff indicates an increasing contribution from the northern GrIS basins (Figure S3), as also noted in Muntjewerf et al. (2020) for SSP5‐8.5 scenario forcing. The southern GrIS basins, however, remain as primary contributors to the total runoff throughout the simulation.

4 Figure. Comparison of the evolution of solid (solid line) and liquid (dashed line) freshwater fluxes from the GrIS (m3 s−1, black lines) and NAMOC index (Sv, blue line). Solid fluxes correspond to ice discharge; liquid fluxes are the sum of runoff and basal melt as computed by the land model and the ice sheet model, respectively.

A similar relationship between NAMOC decline and GrIS freshwater fluxes is found under SSP5‐8.5 forcing. Muntjewerf et al. (2020) compared NAMOC index evolution in two CESM2.1 simulations, with and without an interactive GrIS, and found similar indexes for the two simulations.

SMB Evolution and Ablation Area Expansion

The pre‐industrial GrIS SMB is 585 Gt yr⁻¹ (Table 1), which is higher than present‐day SMB (Fettweis et al., 2017; Noël et al., 2015, 2016), primarily due to a larger ice sheet and overestimated snowfall (Lofverstrom et al., 2020; van Kampenhout et al., 2020). In the 1% simulation, the evolution of surface mass loss can be roughly divided into three periods (Figure 5c, black line, and Table S1), similar to the total mass loss in section 3.1. SMB decreases by 3.5 Gt yr⁻² until year 119, by 13.9 Gt yr⁻² for years 120–226, and by 5.4 Gt yr⁻² during years 226–350. The anthropogenic SMB signal emerges from background variability (i.e., the 20‐year running mean is 2 standard deviations below the pre‐industrial control mean Fyke et al., 2014) by year 84, when the global mean temperature anomaly is 2.5 K. The SMB becomes negative by year 96, at a warming of 2.9 K, close to the 2.7 K threshold found by Gregory et al. (2004).

5 Figure. Change in annual SMB and JJA surface energy balance over the GrIS. (a) Ablation area (%; as percentage of total ice sheet area), (b) cumulative ice sheet area (x axis) below a certain elevation (m; y axis), (c) annual GrIS‐integrated SMB components (Gt yr−1), (d) JJA area‐weighted averages of surface energy balance components (W m−2). The thick lines in (a), (c), and (d) show the 20‐year running means.

Figure 5a shows the evolution of the ablation area, defined as the area with average SMB < 0. The ablation area is expressed as a percentage of the total GrIS area, which decreases with time during the simulation. Starting from a value of 5.5% (1.1 × 10⁵ km²) in the pre‐industrial simulation, the ablation area grows rapidly, with three distinct trends whose timing is different from that of the overall SMB trends. In the first 98 years, the ablation area expands by around 0.1% yr⁻¹. The anthropogenic‐forced signal emerges from background variability in year 46, when the global mean temperature anomaly is 1.1 K. In a CESM2.1‐only simulation (without an interactive ice sheet) under the same scenario forcing (Sellevold & Vizcaino, 2020), this ablation area signal emerges before the SMB signal due to lower variability. During years 99–192, the rate of expansion triples to 0.3% yr⁻¹; by years 131–150, the ablation region makes up 24.2% (4.8 × 10⁵ km²) of the total ice sheet area. Between years 193 to 350, the trend is again 0.1% yr⁻¹, and by end of simulation, the ablation area is 60.1% (10.1 × 10⁵ km²) of the total GrIS area.

In order to relate this ablation area expansion to the ice sheet hypsometry, Figure 5b shows the cumulative area below a certain elevation. Because of the quasi‐parabolic shape of the GrIS, this area increases rapidly with elevation. Mapping ice sheet elevation to the mean equilibrium line altitude (ELA, where SMB = 0), and the area below this mean ELA to the ablation area, it is clear that a given increase in ELA represents more ablation area expansion after a certain ELA threshold is reached. This in part explains the rapid expansion of the ablation area after year 99, when the mean ELA is approximately 800 m. The change in ELA from 1,000 to 2,000 m with the pre‐industrial topography (black line) expands the ablation area with 33% (this percentage is similar for the slightly different hypsometries at years 131–150 and 331–350).

Next, we consider ELA evolution for different regions of the ice sheet. We examine three 400 km long cross sections following the topographic gradient in the southwest (67°N) and in the north (Petermann and NEGIS catchments) (Figure 6).

6 Figure. Hovmöller plots of annual mean SMB (m yr−1; colors) and annual mean topography (m; contours) for three transects on the GrIS. The transects are shown in the maps to the right, where the black line corresponds to the x axis on the left side. Note the different scaling of the positive and negative colors.

In the southwest, the initial ELA is at ∼1,400 m. The ELA moves upwards by 400 m around year 120 and by another 300 m around year 150, in a quasi step‐wise fashion. This 700‐m increase brings the ablation area 100 km farther inland. The next 400‐m increment, raising the ELA to above 2,500 m, takes longer (90 years) and brings the ablation area another 100 km inland. During the entire simulation, the topography lowers and steepens as a result of the more negative SMB. For example, at the start of the simulation, the distance from the 2,000‐m contour to the 2,500‐m contour is about 100 km, while at end of simulation, this distance is only 80 km. For the contour line at 1,500 m (near x = 100 km at the start of the simulation), the SMB decreases from near zero at the start to ‐4 mm yr⁻¹ at year 350. If we follow the SMB evolution at x = 100 km, the SMB declines to ‐8 mm yr⁻¹ by end of simulation. This shows similar contributions of the direct climate change and the elevation feedback to increased melting during the simulation.

In the northern sections, the topographic gradient is less than in the southwest. The ELA increase starts earlier (around year 70 in NEGIS and year 100 in Petermann) and is larger than in the southwest, as the initial ELA is much lower (∼500 m in NEGIS and ∼800 m in Petermann) and the final ELA is only slightly lower. The ablation area migrates inland by an additional 60–80 km compared to the southwest.

Figure 5c and Table 2 show the time evolution of the SMB components. The total precipitation rate increases over the course of the simulation, but the signal emerges relatively late (year 202, when the global mean temperature increase is 6.8 K). This is due to the combination of global warming and reduced NAMOC signals (Figures 1c and 3) in the Greenland region, with the latter reducing the precipitation in the southern parts of the ice sheet (through its influence on the NAO Berdahl et al., 2018; Peings & Magnusdottir, 2014) and partly compensating the moderate precipitation increases elsewhere (Sellevold & Vizcaino, 2020). Snowfall, unlike total precipitation, decreases during the simulation, but does not emerge from background variability. This decrease results from a greater fraction of precipitation falling as rain (from 9% pre‐industrial to 39% by end of simulation), as a result of warming. The slightly reduced snowfall together with increased rain impacts refreezing, as more rain can be refrozen while the refreezing capacity remains nearly constant. Sellevold and Vizcaino (2020) give a more detailed analysis with spatial maps for a 150‐year CESM2.1‐only simulation under the same scenario forcing but with prescribed GrIS topography.

2 TableAnnual Ice Sheet‐Integrated Surface Mass Balance and Components Mean (Standard Deviation) and Anomalies of the Mean With Respect to Pre‐Industrial (Gt yr⁻¹)

		Years 131–150	Years 331–350
Component	Pre‐industrial	Absolute	Anomaly	Absolute	Anomaly
SMB (4 km)	585 (85)	–367 (166)	–952	–2,259 (357)	–2,844
SMB (1°)	544 (103)	–521 (217)	–1,065	–2,589 (442)	–3,133
Precipitation	846 (83)	986 (97)	140	1,122 (97)	276
Snowfall	780 (80)	750 (74)	–30*	683 (71)	–97
Rain	72 (12)	235 (38)	163	439 (59)	367
Refreezing	223 (54)	693 (73)	470	534 (43)	311
Melt	415 (92)	1,914 (251)	1,499	3,804 (443)	3,389
Sublimation	45 (4)	50 (6)	5	3 (11)	–42
Rain (%)	8 (1)	24 (3)	16	39 (4)	31
Refreezing (%)	46 (4)	32 (3)	–14	13 (1)	–33

Note. SMB (1°) values are calculated as the sum of components in CLM5. SMB (4 km) values are in CISM, after downscaling and remapping. SMB (1°) = Snowfall + Refreezing – Melt – Sublimation. Rain (%) = Rain * 100 / (Snowfall + Rain). Refreezing (%) = Refreezing * 100 / (Rain + Melt). Differences with the downscaled SMB used by CISM2.1 (Table 1) are due to mass definition across components, for mass conservation purposes (see, e.g., Vizcaino et al., 2013). *All changes in the mean are significant (p < 0.05) except snowfall by 131–150.

Melt increases from the start of the simulation and accelerates after the first century. By CO₂ stabilization, the total GrIS melt is five times greater than the pre‐industrial value (Table 2). Melt continues to increase until year 280, reaching nine times the pre‐industrial value by end of simulation. Refreezing also increases from the start of the simulation. This is mostly due to increased available liquid water from surface melt and rainfall, with melt representing the largest contribution (90% by end of simulation). The refreezing capacity—defined as the fraction of refreezing to available melt water—decreases from 46% pre‐industrial (in agreement with estimates from RACMO, Noël et al., 2018) to 32% by years 131–150, in broad agreement with RCP4.5 projections (van Angelen et al., 2013). After stabilization, the refreezing no longer increases, despite further increases in available water, because the capacity of snow to store meltwater is saturated. From year 200 to end of simulation, refreezing rates decrease. By end of simulation, the refreezing capacity is only about 13%. The maximum refreezing has values close to but below the total snowfall rate (93% for years 131–150 and 79% for years 331–350), confirming the validity of parameterizations that estimate potential refreezing as a fraction of total snowfall (Aschwanden et al., 2019).

On the Energy Sources of Melt Acceleration

The surface energy balance components (Figure 5d and Table 3) are necessary to explain the melt acceleration after year 120 (Figure 5c) and the subsequent acceleration in SLR contribution (Figure 1). In the first century, the primary source of additional melt energy is the increase in net longwave radiation from increasing atmospheric CO₂ concentration. The net shortwave radiation at the surface does not increase, because reduced incoming radiation from enhanced cloudiness (Sellevold & Vizcaino, 2020) balances the reduced reflected shortwave radiation as the surface albedo decreases. By CO₂ stabilization, longwave radiation is still the largest source of melt (40% of the total melt). By end of simulation, albedo decreases make solar radiation the largest source (39%), followed by the turbulent fluxes (34%).

3 TableSummer GrIS‐Averaged Albedo (–), Near‐Surface Temperature and Skin Temperature (°C), Incoming Shortwave Radiation at the Surface, Incoming Longwave Radiation at the Surface, and Surface Energy Balance Components (W m⁻²), Mean (Standard Deviation)

	Pre‐industrial	Years 131–150	Years 331–350
Albedo	0.78 (0.01)	0.72 (0.01)	0.62 (0.01)
T2m	–7.1 (0.8)	–1.5 (0.5)	0.6 (0.3)
Tskin	–7.6 (0.8)	–2.3 (0.4)	–0.8 (0.2)
SWin	289.6 (3.7)	264.4 (5.2)	252.6 (6.2)
LWin	231.3 (3.7)	266.6 (3.5)	279.7 (3.4)
Melt energy	8.2 (2.0)	38.2 (5.0)	83.1 (9.1)
SWnet	62.5 (2.3)	71.3 (3.4)	91.4 (4.4)
LWnet	–49.8 (2.0)	–37.7 (2.7)	–31.4 (2.8)
SHF	5.0 (1.0)	9.6 (1.9)	20.8 (2.9)
LHF	–7.8 (0.4)	–6.3 (1.0)	2.1 (2.1)
GHF	–1.7 (0.3)	1.2 (0.5)	0.2 (0.4)

Note. Melt energy = net shortwave radiation SW_net + net longwave radiation LW_net + sensible heat flux SHF + latent heat flux LHF + ground heat flux GHF. All changes in the mean are significant (p < 0.05).

A threshold—or tipping point—in the melt energy is reached around year 120. The net solar and turbulent heat fluxes substantially increase, while the net longwave continues a smoother increase, as the global atmosphere warms (Figure 1c). The strong increases in solar and turbulent fluxes result from the combination of two processes. First, the ice‐albedo feedback is triggered and amplifies the melt increase as the ablation area expands (Figure 5a), exposing more bare ice. Bare ice has an albedo of ∼0.4, much lower than fresh snow and wet snow, which have albedos of 0.85–0.90 and 0.65–0.75, respectively. Second, the global mean temperature increase exceeds a certain threshold (4.2 K) that is regionally translated into summer GrIS mean temperatures close to the melting point (Table 3). Large parts of the ice sheet surface are at the melting point, while near‐surface air temperatures can be above the melting point. This results in a stronger surface temperature inversion and enhanced turbulent fluxes.

In the following, we examine the spatial changes in the near‐surface Greenland summer (JJA) climate, including the nonglaciated regions, with a focus on albedo and turbulent heat fluxes. The model accounts for the land cover change involved in the transition from glacier to bare land or vegetation as the margins retreat. The GrIS loses 3% of its pre‐industrial area by years 131–150 and 19% by years 331–350 (Figure 1f). The ablation area expansion was shown in Figure 5a.

In the pre‐industrial summer when the ablation regions are narrow, most of the ice sheet area is covered with snow, and the island of Greenland has a mean albedo of 0.71 (Figure 7a). More bare ice is exposed as the ablation regions widen. By 131–150 and 331–350, the overall Greenland albedo decreases to 0.64 and 0.50, respectively, from GrIS retreat and the low albedo of the expanding tundra and bare land (Figures 7b and 7c). The GrIS summer albedo decreases from 0.78 to 0.72 and 0.62 over the same periods (Table 3).

7 Figure. JJA Greenland climate for (left row) pre‐industrial (CTRL) (1–300), (middle row) years 131–150, and (right row) 331–350, with (a–c) albedo (–), (d–f) latent heat flux (W m−2), and (g–i) sensible heat flux (W m−2). Positive heat flux indicates heat transfer from atmosphere to land. Pink contour encloses areas with >70% glacier fraction within the grid cell.

The latent heat flux represents energy transfer due to the phase change of water—defined as positive when directed to the surface. The sign depends on the humidity gradient in the surface layer and the temperature‐dependent saturation point. The pre‐industrial summer latent heat flux (Figure 7d) is of similar magnitude to that obtained by Ettema et al. (2010) with the RACMO regional model during 1958–2008: evaporation of about –40 W m⁻² over the west tundra and sublimation of about –10 W m⁻² over the ablation zones. The latent heat flux over tundra in boral summer (JJA) is negative and becomes more negative as the simulation progresses, indicating an increased evaporation as the near‐surface air warms (Figures 7d to 7f). The latent heat flux over the ice sheet interior is negative as well, indicating sublimation in heated but nonmelting regions. The ablation areas show pronounced sublimation in the pre‐industrial era. However, as near‐surface air warms compared to the surface, and the temperature inversion in the surface layer develops and strengthens, there is deposition or condensation. This implies that moist air cools as it flows over the cold surface and saturates, such that excess water vapor directly condenses (hoarfrost) or deposits on the ice. Integrated over the entire GrIS, the deposition becomes the dominant process as the melt increases and the ablation areas expand: in the time series of GrIS summer latent heat (Figure 5d, red line), the flux sign changes from negative (sublimation) to positive (deposition). Deposition increases the melt flux as it provides extra energy to the surface, but little extra mass.

The sensible heat flux represents energy transfer associated with warming or cooling of the surface—defined as positive when directed to the surface. The sensible heat flux over the tundra is negative: the tundra warms the near‐surface air (Figures 7g to 7i). As ice‐covered surfaces transition to tundra over the course of the simulation, the region of negative fluxes expands. The ice sheet interior sensible heat flux is negative as well, indicating warming of the air, though less than over the tundra. At the ice sheet margins, where the air is warmer than the ice, the sensible heat flux is positive. The margin surface heat flux becomes more positive during the simulation and maps well with the expansion of the ablation areas (Figures 6a to 6c). This process dominates the long‐term increase in GrIS summer sensible heat flux (Figure 5d, green line).

Change in Ice Flow and Discharge

Map of Mass Loss and Velocity Changes

Figures 8 and S4 show the spatial distribution of mass change and its components (SMB and ice discharge change), as well as changes in surface velocity. The ablation zone expands throughout the simulation, raising the ELA to around 2,000 m by years 131–150 and 2,500 m by years 331–350 (Figures 8b and 8c). The expansion of the ablation zone results in extensive ice sheet thinning (Figures 8e and 8f), which in turn leads to increased ice flow from the interior towards the margins (Figures 8h and 8i). This ice flow acceleration, which further contributes to ice thinning in the expanded ablation zone, is due to the formation of steeper slopes as the ELA rises. In fact, the ice sheet thickens somewhat in the interior, as a result of local increases in snowfall from an enhanced hydrological cycle (spatial maps of changes in snowfall and precipitation for a fixed present‐day GrIS topography are shown in Sellevold & Vizcaino, 2020). As a result, the slope angle increases substantially where the ablation and accumulation zones meet. This increases the driving stress, giving higher velocities in the transition area between the high interior and the thinning margins. While faster flow partly reduces mass loss at the margins as ice advection from the interior increases, it favors inland migration of the equilibrium line and thinning upstream (e.g., Vizcaino et al., 2015). The gross pattern of SMB, velocity, and thickness change by year 140 is similar to the changes by 2081–2100 under the SSP5‐8.5 scenario (Muntjewerf et al., 2020), when the atmospheric CO₂ concentration is similar. The SSP5‐8.5 simulation, however, reaches a more negative SMB (–565 vs. –367 Gt yr⁻¹) due to a stronger increase in CO₂ forcing late in the simulation (after 2070). As the ice sheet margins thin and retreat, outlet glaciers decelerate, resulting in almost 200 Gt yr⁻¹ lower discharge by years 131–150 (Table 1 and Figure S4). By end of simulation, ice sheet retreat results in a generally terrestrial margin in all basins. The northwestern outlet glacier termini become terrestrial, and discharge exceeds 5 Gt yr⁻¹ only for the Jakobshavn, Petermann, Helheim, and Kangerlussuaq glaciers and the NEGIS.

8 Figure. Spatial change over the GrIS as simulated in the ice sheet model (CISM2.1) for pre‐industrial (left column) and differences with respect to the former by model years 131–150 (middle column) and 331–350 (right column). (a) Surface mass balance (kg m−2 yr−1) with accumulation zones (SMB > 0) and ablation zones (SMB < 0), (b) ice sheet thickness (m), and (c) surface velocity (m yr−1).

Changes in Major Outlet Glaciers

To further analyze changes at the GrIS margins, we examine flowline sections of seven major outlet glaciers draining different basins: Nioghalvfjerdsfjord glacier and Zachariae glacier in the northeast (NE) basin, Petermann glacier and Humboldt glacier in the north (NO) basin, Kangerlussuaq glacier and Helheim glacier in the southeast (SE) basin, and Jakobshavn glacier in the central‐west (CW) basin. Figure 9 shows the location of these flowlines, as well as the timing of margin retreat and the comparison between surface velocity maps at pre‐industrial and end of simulation.

9 Figure. A flowline section of seven selected glaciers. Evolution of (bottom panels) ice thickness, (central panels) ice velocity, and (top panels) SMB. Each line corresponds to a simulation year, and years 0, 140, and 351 are highlighted in blue, yellow, and red, respectively (see time label bar). For clarity, different scales for ice velocity are used for different glaciers. The transect locations of the seven selected glaciers are illustrated by the light blue lines in the two maps at the bottom right. Left map: initial ice velocity (m yr−1); blue line indicates the initial ice margin position. Right map: ice velocity in year 350 (m yr−1); blue lines as in the left panel, light to dark red shading between the initial and final margin indicate the ice margin position throughout the simulation (see time label bar).

In the southeast basin, the margins of Kangerlussuaq and Hellheim glaciers do not retreat during the simulation (Figure 9 and Table S2). At year 140, the SMB is still positive over a large portion of the glaciers. At this time, only the lower part has a negative SMB, with a rapid downstream decline and values as low as –1.5 m yr⁻¹ at the glacier terminus. A nearly identical, relatively steep downstream SMB gradient is simulated until year 350, when values range from –0.5 m yr⁻¹ inland to –2.5 m yr⁻¹ at the terminus. At the beginning of the simulation, ice velocity is more than 2 km yr⁻¹ in regions with steep bedrock topography. At this time, in Helheim glacier, the ice velocity increases smoothly from 1 to 4 km yr⁻¹ in the lower part of the glacier. Similar velocities are simulated for Kangerlussuaq glacier, although the downstream increase is not as smooth as for Helheim, but rather has individual peaks between 2 and 4 km yr⁻¹. For both glaciers, the ice velocity smoothly declines with time, and at end of simulation, the peaks are below 2 km yr⁻¹.

In the north basin, Petermann glacier begins to retreat relatively late, after year 246. In the last 100 years, the glacier margin migrates inland by 37 km, without losing contact with the ocean. The lower part of the glacier already has a negative SMB at year 140. At this time, however, the SMB is slightly positive in the upper part of the glacier (0–40 km along the transect), with a gentle downstream decline to –1 m yr⁻¹ at the terminus. At end of simulation, a negative SMB of around –1 m yr⁻¹ is simulated everywhere along the glacier transect. The terminus velocity peaks around 1 km yr⁻¹ at pre‐industrial and declines over time, with episodic increases after margin retreats. In the same basin, Humboldt glacier starts retreating much earlier than Petermann, with the first retreat episode around year 184. By end of simulation, Humboldt has become land terminating, with a margin retreat of around 60 km. The SMB is already negative at year 140 over the entire glacier length, with a smooth gradient ranging from negative values near zero to –1 m yr⁻¹ at the terminus. At end of simulation, the negative SMB is around –1 m yr⁻¹ over the whole glacier length. At the beginning of the simulation, Humboldt glacier's maximum ice velocity is relatively low (700 m yr⁻¹) compared to other Greenland major drainage systems. The pattern of overall velocity decrease, with only episodic speed‐ups after retreats, is similar to Petermann.

In the northeast basin, the Nioghalvfjerdsfjord glacier starts to retreat around year 159. At end of simulation, the glacier margin has retreated by 46 km, with the largest part of the retreat occurring between years 270 and 350. A similar retreat of around 50 km is simulated for the Zachariae glacier at year 350, although the initial retreat occurs later, around year 180. Both glaciers remain in contact with the ocean at end of simulation, as the fjord extends upstream by several tens of kilometers. The SMB along the flowline is negative at year 140 for both glaciers, varying smoothly from negative values close to zero to values around –0.5 m yr⁻¹ downstream. In the last two decades, the model simulates a negative SMB close to –1 m yr⁻¹ everywhere along the transect. For both glaciers, the ice velocity near the margin peaks around 1 km yr⁻¹ at pre‐industrial and declines over time, with episodic increases in terminus velocity after margin retreat.

In the central‐west basin, Jakobshavn glacier starts to retreat relatively late, at year 271. At end of simulation, the glacier margin has retreated by only 20 km. At year 140, the SMB is negative everywhere along the glacier length, with a relatively sharp downstream gradient from slightly negative values inland to values below –1 m yr⁻¹ at the terminus. By end of simulation, the SMB is negative over the whole glacier length, ranging between –1 and –2 m yr⁻¹. At pre‐industrial, simulated ice velocities are larger than 1 km yr⁻¹ in the lower part of the glacier, with a sharp peak of 4 km yr⁻¹ near the terminus. Ice velocities exceed 1 km yr⁻¹ downstream of a local high in the bedrock, after which the bedrock topography is more steep. Similar to other glaciers, the overall velocity decreases with time, with episodic speed‐ups after margin retreats. At end of simulation, peaks in ice velocity are lower than 2 km yr⁻¹.

In summary (Table S2), the sensitivity of these major outlet glaciers to the simulated climate change is heterogeneous, with the relatively slower, drier basin‐draining northern glaciers retreating the most, and the relatively faster, wetter basin‐draining southeastern glaciers not retreating by the end of the simulation. Humboldt retreats the most (60 km) and is the only one of the seven glaciers to become fully terrestrial. Petermann is the latest northern glacier to start retreating. Jakobshavn Isbrae, in the central‐west, retreats less than the northern glaciers and starts retreating later.

Summary and Discussion

The main novelty of this study comes from the use of an advanced coupled Earth system/ice sheet model, CESM2.1‐CISM2.1, to investigate future, multicentury, climate and GrIS evolution. The model has relatively high resolution (∼1°) in the ocean and atmospheric components and an advanced SMB calculation with energy fluxes computed directly in the climate component. This study therefore bridges the gap between multicentury and multimillennia ice sheet model‐only SLR projections (e.g., Aschwanden et al., 2019) and century‐scale SMB projections from global and regional climate models.

This study investigates a multicentury (350‐year) idealized 4×CO₂ scenario. By year 150—one decade after CO₂ stabilization at 4× the pre‐industrial concentration—the projected GrIS contribution to global mean SLR is 107 mm SLR, with a global mean near‐surface temperature anomaly of 5.3 K relative to pre‐industrial. The polar amplification factor is 1.6 for the Arctic (north of 60°N), but only 1.2 for the GrIS. The NAMOC strongly declines, starting before a substantial increase in GrIS melting and runoff. After another 210 years of stable, elevated atmospheric CO₂, the total projected GrIS contribution to SLR is approximately 1,140 mm, and the global mean temperature anomaly is around 8.5 K.

We compare our projected GrIS SLR contribution by year 150 with previous estimates under the CMIP RCP8.5 and SSP5‐85 scenarios, in which the CO₂ concentration is close to 4×CO₂ by year 2100. We focus specifically on CESM studies, including coupled simulations and ice sheet‐only projections. Using CESM1.0, Vizcaíno et al. (2014) projected 55 mm by 2100 from the SMB contribution alone, under a global mean warming of 3.7 K. When forcing a stand‐alone ice sheet model (CISM1.0; a comparatively low‐complexity predecessor to CISM2.1) with the aforementioned CESM1.0 SMB field, the total mass loss—including both SMB and ice discharge—was 76 mm (Lipscomb et al., 2013). The coarse‐resolution Earth system/ice sheet model study by Vizcaino et al. (2015) projected 2100 SLR and warming of 67 mm and 4.3 K, respectively. Moreover, in an ice sheet model‐only study with PISM, Aschwanden et al. (2019) estimated a GrIS SLR contribution of 140–330 mm by year 2100, with a global mean temperature change around 5 K. Finally, Muntjewerf et al. (2020) projected 109 mm GrIS SLR under the SSP5‐85 scenario with CESM2.1‐CISM2.1 (same model as in this study).

The lack of direct ocean forcing of the GrIS is a limitation of this study. Fürst et al. (2015) included ocean temperature forcing at calving fronts along with surface forcing and estimated 102 mm SLR under the CMIP5 RCP8.5 scenario. Mass loss due to enhanced ice dynamics from ocean forcing is estimated to be an order of magnitude smaller: Price et al. (2011) found 6 mm of committed dynamic mass loss from ocean forcing by 2100. Nick et al. (2013) projected a dynamic contribution between 11 and 18 mm SLR by 2100 from the four largest outlet glaciers. On the other hand, RCP8.5 simulations without ocean forcing give reduced ice discharge at 2100 relative to present day (Ruckamp et al., 2019; Vizcaino et al., 2015), illustrating the important role of the ocean. Although ocean interactions are second order compared to SMB for overall GrIS mass loss (on the whole ice sheet scale), the above studies suggest that ocean forcing enhances solid ice discharge, which otherwise decreases. The question remains what the net effect is; answering this question requires the development of models with GrIS ice‐ocean interactions. To better understand ocean terminus and ocean shelf processes, the ISMIP6 stand‐alone ice sheet model GrIS experiments are provided with scenario ocean boundary forcing (Slater et al., 2019) to accommodate parameterizations of marine terminus retreat and submarine melt.

To date, relatively few studies have investigated GrIS behavior on multicentury time scales with forcing from ESMs. This requires either an expensive simulation with a coupled Earth system/ice sheet model or an ice sheet‐only simulation with appropriate SMB downscaling to account for the growing divergence between the static GrIS geometry in the climate model and the evolving geometry in the ice sheet model. An example of the former is the relatively coarse‐resolution, coupled model in Vizcaino et al. (2015) that projected a GrIS SLR contribution of 536 mm by year 2300 under a simulated temperature anomaly of 9.4 K. In contrast, the RCP8.5 ice sheet model simulations by Aschwanden et al. (2019) estimated a GrIS SLR contribution of 940–3740 mm under a 10‐K warming in year 2300. Our projection lies at the lower end of this range.

In our CESM2.1‐CISM2.1 projection, the rate of SLR increases from 2 mm yr⁻¹ in year 150 to almost 7 mm yr⁻¹ by year 350. The accelerated mass loss is caused by a rapidly declining SMB. Part of this decline is compensated by reduced ice discharge, as the GrIS margins thin and retreat and many outlet glaciers become land terminating. The final GrIS area reduction is around 20%. The SMB decreases slowly at the start of the simulation, when increased longwave radiation is the primary source of melt as the atmosphere CO₂ concentration is increasing. Surface ablation and overall mass loss accelerate after about 120 years of increasing CO₂, corresponding to a global warming of 4.2 K, as a result of increasingly high surface melt. During this time, the ablation areas expand enough to trigger a large albedo feedback as the net shortwave radiation at the surface increases. In addition, increased turbulent heat fluxes (sensible and latent) in the summer seasons further amplify GrIS melt, and near‐surface temperature inversions expand and strengthen. With the increase in meltwater, the GrIS refreezing capacity sharply decreases, enhancing runoff. By the end of the 350‐year simulation, the ablation region covers around 60% of the total GrIS area.

In agreement with previous work, our results show how sustained high CO₂ results in large, accelerating mass loss. We have shown how coupled climate and ice sheet modeling via advanced SMB calculations provides novel insights into both short‐ and long‐term coupling of global and regional climate change and ice sheet melt change. Adequate modeling of global climate as well as snow and ice processes is key to capturing the interactions and feedbacks that will determine future changes.

Acknowledgments

L. M., M. P., M. V., and C. E. D. acknowledge funding from the European Research Council (grant ERC‐StG‐678145‐CoupledIceClim). R. S. acknowledges funding from the Dutch Research Council (NWO) (grant ALWOP.2015.096). The CESM project is supported primarily by the National Science Foundation (NSF). This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under Cooperative Agreement No. 1852977. Computing and data storage resources, including the Cheyenne supercomputer (doi:10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF. We thank the ISMIP6 steering committee, the ISMIP6 model selection group, and ISMIP6 data set preparation group for their continuous engagement in defining ISMIP6. This is ISMIP6 contribution no. 18.

Data Availability Statement

CESM2 is an open source model, available at https://www.cesm.ucar.edu/.

The World Climate Research Program (WGCM) Infrastructure Panel is the official CMIP document home: https://www.wcrp-climate.org/wgcm-cmip. The CMIP6 and ISMIP6 simulations are freely available and accessible via the Earth System Grid Federation (ESGF) data portals https://esgf.llnl.gov/nodes.html.