Greenland Ice Sheet Contribution to 21st Century Sea Level Rise as Simulated by the Coupled CESM2.1‐CISM2.1

The Greenland Ice Sheet (GrIS) mass balance is examined with an Earth system/ice sheet model that interactively couples the GrIS to the broader Earth system. The simulation runs from 1850 to 2100, with historical and SSP5‐8.5 forcing. By the mid‐21st century, the cumulative GrIS contribution to global mean sea level rise (SLR) is 23 mm. During the second half of the 21st century, the surface mass balance becomes negative in all drainage basins, with an additional SLR contribution of 86 mm. The annual mean GrIS mass loss in the last two decades is 2.7‐mm sea level equivalent (SLE) year−1. The increased SLR contribution from the surface mass balance (3.1 mm SLE year−1) is partly offset by reduced ice discharge from thinning and retreat of outlet glaciers. The southern GrIS drainage basins contribute 73% of the mass loss in mid‐century but 55% by 2100, as surface runoff increases in the northern basins.


Introduction
During the past two decades, the Greenland Ice Sheet (GrIS) has lost mass at an increasing rate (Shepherd et al., 2019) and has become a major contributor to global mean sea level rise (SLR) (Chen et al., 2017). The polar ice sheets were identified in the IPCC Fifth Assessment Report as large sources of uncertainty in 21st century SLR projections (Church et al., 2013). A recent expert assessment estimates that by 2100, the GrIS will contribute between 20-and 990-mm sea level equivalent (SLE), with a median of 230-mm SLE (Bamber et al., 2019). The uncertainty in SLR stems partly from insufficient understanding of the complex interactions and feedbacks among ice sheets, surface mass balance (SMB), and climate (Fyke et al., 2018), highlighting the importance of studies with coupled Earth system/ice sheet models (Goelzer et al., 2017;Vizcaino, 2014).
Many studies on future ice sheet evolution have been done using regional climate models with static ice sheet topography (Lenaerts et al., 2015;Rae et al., 2012;Tedesco & Fettweis, 2012) or dynamic ice sheet models with prescribed climate forcing (e.g., Aschwanden et al., 2019;Graversen et al., 2011;Ruckamp et al., 2019). By accounting only for climate-SMB interactions or ice sheet-SMB interactions, these studies miss interactions between ice sheet dynamics and climate. Further, some studies with partially coupled ice sheet-climate models have used positive-degree-day schemes, which parameterize melt based on the number of days the temperature is above the melting point (e.g., Huybrechts et al., 2011;Mikolajewicz et al., 2007;Ridley et al., 2005;Vizcaino et al., 2008). In new climate regimes, however, empirical relations between melt and temperature may no longer hold, and an explicit surface energy balance (SEB) is necessary (Bauer & Ganopolski, 2017).
Physics-based models that explicitly resolve the SEB and snowpack processes are needed to simulate feedbacks, such as the ice-albedo feedback (Bougamont et al., 2007). A main challenge is the fine horizontal grid resolution needed to realistically simulate SMB gradients (Van den Broeke et al., 2008). There has been recent progress in coupling an SEB-based SMB calculation with ice dynamics using regional climate models of both intermediate complexity (Born et al., 2019;Krapp et al., 2017) and higher complexity (e.g., regional climate model MAR Le clec'h et al., 2019). However, the use of global models with SEB-based SMB to simulate Earth system interactions with ice sheets has been limited (Vizcaino et al., 2008. In this study, the Community Earth System Model version 2 (CESM2 Danabasoglu et al., 2020) with an interactive GrIS model (CISM2.1 Lipscomb et al., 2019) is used to simulate the period 1850-2100 under historical and SSP5-8.5 forcing. The model includes an SEB-derived SMB calculation with multiple elevation classes (Fyke et al., 2011;Lipscomb et al., 2013;Sellevold et al., 2019) and, thus, is able to capture ice sheet-SMB-climate feedbacks. The simulations follow the protocols for coupled ice sheet-climate model experiments in the framework of the Ice Sheet Modeling Intercomparison Project for CMIP6 (ISMIP6 Nowicki et al., 2016). This paper presents the CESM2.1-CISM2.1 projection of 21st century global climate change and GrIS response, as well as the projected GrIS contribution to global mean sea level. Results are compared to CESM2.1 simulations without an interactive ice sheet.

The Community Earth System Model
The CESM2 (Danabasoglu et al., 2020) is a comprehensive, fully coupled Earth system model that is contributing simulations of past, present, and future climates to the Coupled Model Intercomparison Project phase 6 (CMIP6 Eyring et al., 2016). CESM2 includes component models of the atmosphere (CAM6), land (CLM5 Lawrence et al., 2019), ocean (POP2 Smith et al., 2010;Danabasoglu et al., 2012), sea ice (CICE5 Hunke et al., 2017), river transport ( MOSART Li et al., 2015), and land ice (CISM2.1 Lipscomb et al., 2019). The simulations described here were run with nominal 1-degree horizontal resolution (∼110 km at the equator) in the atmosphere, land, ocean, and sea ice components. The ice sheet model was run on a 4-km limited-area grid centered on Greenland.

Interactive Earth System/Ice Sheet Coupling
CESM2.1-CISM2.1 supports an evolving Greenland Ice Sheet that is interactively coupled to other Earth system components. The SMB is computed in CLM5 as the difference between annual snow accumulation and surface ablation, derived from the surface energy balance and snow pack hydrology. The SMB is calculated for multiple elevation classes in glaciated grid cells to account for subgrid-scale variations in surface climate and SMB components (Fyke et al., 2011;Lipscomb et al., 2013;Sellevold et al., 2019). The SMB is remapped from the CLM5 grid to the higher resolution CISM2 grid, using a trilinear interpolation scheme that conserves both the total ablated and the total accumulated mass (Leguy et al., 2018).
Freshwater fluxes from the GrIS to the ocean are the sum of surface runoff from CLM5 and basal meltwater and ice discharge from CISM2. Surface runoff is routed to the ocean via MOSART, accounting for the surface slope, and basal meltwater reaches the ocean by nearest neighbor routing. In the ocean, meltwater is distributed over the upper 30 m (Sun et al., 2017). Solid ice discharge (i.e., calving) is also routed to the nearest ocean neighbor but is spread diffusively in the ocean surface layer over a region extending up to 300 km from the coast and is melted instantaneously using energy from the ocean surface. Ocean thermal forcing to calving fronts is not simulated. Greenland's floating ice shelves, such as Petermann, are not simulated, and thus, the basal melt term refers only to grounded ice. This term is small compared to mass loss from surface runoff and calving and is disregarded in the remaining discussion.
Dynamic land units in CLM5 enable the transition from glaciated to non-glaciated land cover, consistent with the evolving ice sheet margin in CISM2. The ice sheet surface topography in CISM2 is used to recompute the glacier coverage in CLM5 at runtime, modifying the albedo, soil, and vegetation characteristics. Surface elevation and topographic roughness fields in CAM6 are updated every 10 years to incorporate ice sheet geometry changes in atmospheric calculations. The coupling is further described in the supporting information, Text S1.

Experimental Setup
We analyze two simulations: a historical run from 1850 to 2014 and its continuation to 2100 following the SSP5-8.5 scenario (Nowicki et al., 2016;O'Neill et al., 2016). The historical forcing is a reconstruction based on observations of greenhouse gas concentrations, stratospheric aerosol data (volcanoes), land use change, and solar insolation. The preindustrial (1850 CE) CO 2 concentration is 287 ppmv (parts per million by volume), increasing to 397 ppmv in 2014. Further details on the forcing protocol can be found in Eyring et al. (2016). The SSP5-8.5 scenario starts in 2015 at the end of the historical period and ends in 2100 when the atmospheric CO 2 concentration is 1,142 ppmv, and the radiative forcing is 8.5 W m −2 relative to preindustrial.
The initial condition of the historical simulation is the spun-up model state. A description of the spin-up procedure is provided in Text S2. At the beginning of our simulation, the GrIS is in near equilibrium with the simulated preindustrial climate of CESM2.1, with residual drift of about 0.03-mm SLE year −1 . This initial GrIS overestimates the present-day observed area and volume by about 15% and 12%, respectively.

Basin-Scale Analysis
For regional-scale analysis, we use six major Greenland drainage basins defined by Rignot and Mouginot (2012). The basins are directly remapped onto the CISM grid since the modeled drainage divide locations roughly agree with observations (see velocity fields of Figure S2a, showing observations from Joughin et al., 2010, 2015, and Figure S2b). In regions where the ice sheet extent is overestimated, drainage basins are extended to the ice sheet margin. By extending each drainage basin from the margin into the ocean based on the ice flow direction, we also define six major ice-ocean sectors that receive freshwater from each basin.

Results
The analysis focuses on three climatological periods: the contemporary period (average over 1995-2014) from the historical simulation and the mid-century (2031-2050) and end of century (2081-2100) from the SSP5-8.5 simulation.

Evolution of Global Climate and GrIS Mass Budget
The atmospheric CO 2 concentration increases from 287 ppmv in 1850 to 397 ppmv in 2014, 566 ppmv in 2050, and 1142 ppmv in 2100 ( Figure 1, Table S1). The simulated end-of-century global mean near surface temperatures increase is 5.4 K relative to preindustrial. With respect to the contemporary period (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), the simulated global temperature increases by 1.4 K at mid-century and by 4.6 K at end of century. The Arctic amplification (AA) factor, defined as the ratio of temperature change north of 60 • N to the global mean, is 2.0 by mid-century and 1.8 by end of century. The decrease in AA between mid-century and end of century is significant (p < 0.01) and is likely the result of reduced sea ice cover, since AA depends on the magnitude of sea ice change (e.g., Dai et al., 2019). At mid-century, sea ice concentration is still declining throughout the year, while at end of century, the Arctic is ice free in late summer and early fall and thus experiences no further loss.
The simulated North Atlantic Meridional Overturning Circulation (NAMOC; defined as the maximum of the overturning stream function north of 28 • N and below 500 m depth in the North Atlantic basin) is relatively stable during the historical period, with a mean of 24 Sv (blue line in Figure 1b). The only exception is an anomalously stronger overturning circulation in the 1960s, when the NAMOC increases by about 3 Sv. RAPID-AMOC measurements at 26 • N in the Atlantic basin give a mean overturning strength of 17 ± 3.3 Sv for 2000(Frajka-Williams et al., 2019, in good agreement with the simulated value (18.7 ± 2.4 Sv) at this latitude over the same period. During the 21st century, the simulated overturning cell progressively weakens, collapsing to 8.6 Sv by end of century (Table S3).
The simulated climate change drives a positive GrIS contribution to global mean SLR (Figure 1c). The rate of mass loss increases from the preindustrial near equilibrium (0.03-mm SLE year −1 ) to 0.08-mm SLE year −1 during the contemporary period, 0.55-mm SLE year −1 by mid-century, and 2.68-mm SLE year −1 by end of century (Table S1). The associated SLR is 23-mm SLE by 2050 and 109-mm SLE by 2100. Estimates of 2003-2013 Greenland mass loss from the Gravity Recovery and Climate Experiment satellite mission range from 0.77-to 0.79-mm SLE year −1 (Ran et al., 2018;Schrama et al., 2014;Velicogna et al., 2014). The simulated SLR in the contemporary period is underestimated, likely for several reasons. First, a substantial part of the observed mass loss is attributed to ocean forcing (Wood et al., 2018), which is missing in the model. Second, anomalous atmospheric circulation, in particular blocking over Greenland, has contributed to the observed mass loss but is not simulated correspondingly in global climate models (Delhasse et al., 2018;Hanna et al., 2018). Other model biases (e.g., initial ice sheet topography) can explain the remaining differences. The temperature change at the time of mass loss acceleration is approximately 2.7 K with respect to preindustrial. By 2100, the simulated GrIS area and volume decrease by 3% and 1.2%, respectively, compared to the preindustrial ice sheet.
The SMB is the main contributor to GrIS mass loss (Figure 1d, Table S1). By mid-century, the GrIS-integrated SMB is still positive (350 Gt year −1 ) and 214 Gt year −1 lower than in the contemporary period. Based on a linear, long-term trend, the integrated SMB becomes negative by 2077. The expansion of ablation areas (i.e., areas with average SMB < 0) accelerates with similar timing. By mid-century, the SMB is strongly reduced in southern Greenland ( Figure S1), and by end of century, the ablation areas extend far inland around the entire ice sheet, including the northern periphery. At the northern margins, the ablation zone advances later, and progresses faster inland, than at the southern margins. By end of century the northern surface mass loss also has intensified. In the ice sheet interior, the SMB moderately increases due to greater snowfall.
As the SMB becomes more negative, the ice sheet thins (Figures S1b-S1c). Most of the thinning occurs in the south and below the 2000-m elevation contour, while the ice sheet interior moderately thickens (Figures S1e-S1f). Surface velocities increase in the intermediate areas between the high interior and the ice sheet margins (Figures S1h-S1i) where the surface is steeper because of SMB-induced margin thinning. Conversely, surface velocities decrease at the thinner margins. As a result, total ice discharge is reduced by 8% (45 Gt year −1 ) at mid-century and by 33% (189 Gt year −1 ) at end of century compared to the contemporary period (Table S1 and Figure S2). The reduced discharge partly compensates for the mass loss due to more negative SMB (214 and 1,129 Gt year −1 at mid-century and end of century, respectively; see Table S1).

Sea Level Rise Contribution by Drainage Basin
This section presents mid-century and end-of-century mass balance changes for individual drainage basins, with the contemporary mass budget as a reference ( Figure S3). Text S3 gives a brief comparison of simulated and observed outlet glacier ice discharge (Enderlin et al., 2014).
By mid-century, the mean total GrIS mass loss is 169 Gt year −1 (Figure 2), as the result of an SMB decrease partly offset by reduced ice discharge. In all six drainage basins the SMB decreases but remains positive ( Figure S3), with the largest reductions in the SW (66 Gt year −1 ) and SE (80 Gt year −1 ), there is a modest mid-century decrease in the NE basin (32 Gt year −1 ), with smaller reductions in the CW, NW, and NO basins. Taken together, the SMB in the northern basins (NO, NW, NE) decreases by 55 Gt year −1 or about 26% of the six basin total.
At the end of the century (right panel in Figure 2), the mean mass budget of the GrIS has decreased by 935 Gt year −1 compared to the contemporary budget, and the total surface mass balance is reduced by 1,129 Gt year −1 . The mass loss from ice discharge is reduced by 189 Gt year −1 , partly (17%) canceling the surface mass loss. The largest end of century decreases in SMB (290 Gt year −1 ) and ice discharge (81 Gt year −1 ) are in the SE basin, which make the second largest SLR contribution (215 Gt year −1 ; right panel in Figure S3). The SW basin, with the next largest SMB decrease (272 Gt year −1 ) and a smaller reduction in discharge (27 Gt year −1 ), is the largest net contributor to SLR (235 Gt year −1 ). The three northern basins also have large SMB decreases: 203, 145, and 141 Gt year −1 , respectively, for the NE, NO, and NW basins. Together, the northern basins contribute 43% of the total end-of-century SMB decrease and 45% of the total SLR.

Freshwater Fluxes
This section presents the mid-century and end-of-century freshwater fluxes to the ocean from individual drainage basins (Figure 3). The change in total freshwater fluxes is moderate at mid-century, with runoff increasing by less than 200 Gt year −1 (from 427 to 619 Gt year −1 ), and a small decrease in the solid flux (i.e., ice discharge), from 481 to 430 Gt year −1 (Figure 3). By this time, the NAMOC has already decreased by 6 Sv (Figure 1), suggesting that increased Greenland freshwater fluxes play only a minor role in weakening the NAMOC.
By end of century, runoff has more than tripled compared to the contemporary period, increasing to 1,445 Gt year −1 (Figure 3). The reduced ice discharge (260 Gt year −1 ) results in a total freshwater flux that is not quite double the contemporary flux (1,705 Gt year −1 compared to 908). Per basin, the SW and SE regions contribute the most to the total runoff during all periods. However, their relative contribution decreases from 63% in the contemporary period to 48% at end of century. Thanks to a large increase in runoff, the relative freshwater contribution from the northern basins (NW, NE, and NO) increases from 29% in the contemporary period to 44% at end of century.
Throughout the simulation, the SE region has the largest GrIS ice discharge in absolute terms (Figure 3). Relative to the contemporary discharge, all basins have similar reductions of around 10% at mid-century, but by end of century the northern basins have the greatest reductions (62%).

Comparison to CESM2 Simulations Without an Interactive Greenland Ice Sheet
Finally, we compare the SMB and NAMOC responses to the historical ensembles and a suite of scenario simulations conducted using CESM2.1 without an interactive GrIS but with the same SEB-based SMB calculation. We consider 11 ensemble members from the historical period and the following scenario simulations between 2015 and 2100: SSP1-2.6 (two members), SSP2-4.5 (three members), SSP3-7.0 (two members), and SSP5-8.5 (two members), with scenario details provided by O'Neill et al. (2016).
At end of century under SSP5-8.5 forcing, the SMB is almost 400 Gt year −1 lower in CESM2.1 compared to CESM2.1-CISM2.1 (−906 vs. −511 Gt year −1 ). The smaller SMB reduction in CESM2.1-CISM2.1 is likely the result of dynamic removal of high-melt areas on the margins. These areas are allowed to remain in the nonevolving CESM2.1 simulation, favoring an overestimate of SMB decline. The CESM2.1-CISM2.1 10.1029/2019GL086836 response to SSP5-8.5 does, however, exceeds the CESM2.1 response to less extreme scenarios, including SSP3-7.0. Despite the differences in SMB magnitude, the timing of SMB evolution in the two models is similar for SSP5-8.5. NAMOC evolution also is similar in CESM2.1-CISM2.1 and CESM2.1 simulations (Figure 4 and Table S3). The peak NAMOC strength in the second half of the 20th century, mentioned in section 3.1, is present in both simulations. Generally, the NAMOC in CESM2 is sensitive to external forcing in the 21st century, with large reductions in all scenarios including SSP1-2.6. In the SSP5-8.5 projections, there is no discernible differences in NAMOC weakening between simulations with and without dynamical GrIS.

Discussion and Conclusions
This study presents the first simulations performed with CESM2.1-CISM2.1 including an interactive Greenland Ice Sheet. The projected GrIS contribution to SLR by 2100 in the SSP5-8.5 scenario is 109 mm, in general agreement with pre-AR5 multimodel results (Bindschadler et al., 2013) and the AR5 assessment (Church et al., 2013). The latter gives a likely range of 70 to 210 mm. Our projection also lies within the range of post-AR5 estimates from Fürst et al. (2015), Calov et al. (2018), and Golledge et al. (2019), which are of 102 ± 32, 46-130, and 112 mm, respectively. Vizcaino et al. (2015) gave a lower estimate (58 mm), based on a coupled Earth system/ice sheet model of coarse resolution (3.75 • ) with an energy-balance-based melt calculation. Aschwanden et al. (2019) estimated a higher range of 140-330 mm, using an ice sheet model forced with spatially uniform warming. In our projection, the negative SMB trend increases sharply at mid-century, so that most of the SLR contribution is after 2050.
This study projects an increased contribution of the northern basins to total GrIS mass loss after mid-century, in agreement with previous studies.  found a 46% expansion of ablation area in northern Greenland between the early 1990s and 2017, almost twice as much as in the south, because runoff in the north and northeast increased relative to the total runoff due to longer exposure of bare ice. The high response in the north is also explained by the relatively wide ablation zones from low surface slopes ( Van den Broeke et al., 2016). In 21st century RCP4.5 and RCP8.5 projections with the regional climate model MAR with static GrIS topography, Tedesco and Fettweis (2012) also found a high sensitivity of SMB to higher temperatures in the north. Similarly, Lenaerts et al. (2015), in a 21st century RCP4.5 projection with the regional climate model RACMO, projected a strongly nonlinear response of runoff in northern and northeastern Greenland, in connection with low meltwater buffering capacity of the snow pack.