The Greenland Ice Sheet holds 7.2 m of sea level equivalent and in recent decades, rising temperatures have led to accelerated mass loss. Current ice margin recession is led by the retreat of outlet glaciers, large rivers of ice ending in narrow fjords that drain the interior. We pair an outlet glacier–resolving ice sheet model with a comprehensive uncertainty quantification to estimate Greenland’s contribution to sea level over the next millennium. We find that Greenland could contribute 5 to 33 cm to sea level by 2100, with discharge from outlet glaciers contributing 8 to 45% of total mass loss. Our analysis shows that uncertainties in projecting mass loss are dominated by uncertainties in climate scenarios and surface processes, whereas uncertainties in calving and frontal melt play a minor role. We project that Greenland will very likely become ice free within a millennium without substantial reductions in greenhouse gas emissions.

( A ) Ensemble minimum and maximum (thin lines) and mean (thick lines) of RCP 2.6, 4.5, and 8.5 temperature anomalies with respect to 2006–2015 derived from four GCM simulations that extend until 2300. Beyond 2300, the linear 2200–2300 trend was extrapolated to 2500, after which the 2500 value was kept constant (see Materials and Methods). The area between ensemble minimum and maximum is shaded. ( B ) Cumulative contribution to global mean sea level (GMSL) since 2008 (ΔGMSL). ( C ) Rates of GMSL in millimeter sea level equivalent (SLE) per year, M ˙ . ( D ) Contribution of ice discharge to M ˙ given as D ˙ M ˙ = D ˙ / ( D ˙ + R ˙; ) M ˙ , where D ˙ and R ˙ are ice discharge rate and surface runoff rate, respectively. ( E ) Ratio of ice discharge rate and the total of ice discharge rate and surface runoff rate, D ˙ % = D ˙ / ( D ˙ + R ˙ ) . (B to E) Uncertainties are shaded between 16th and 84th percentile of the 500 ensemble members, the solid line is the median, and the thin dashed line is the control simulation. Some simulations under RCP 8.5 lose all ice, thus the 84th percentile of the cumulative contribution tapers out (B) and the rates decline (C). Rates in (C) to (E) are 11-year running means. The ensemble mean is used for the control simulation.

Here, we present a new assessment of the response of the Greenland Ice Sheet to future warming using the Parallel Ice Sheet Model [PISM; ( 14 )], which is capable of reproducing the complex flow patterns evident in Greenland’s outlet glaciers ( 12 ). PISM simulates the evolution of ice geometry and ice flow. We use three extended Representative Concentration Pathways (RCPs) emissions scenarios ( 15 ): a pathway with reduced greenhouse gas emissions that aligns with meeting the goals of the 2015 Paris Climate Agreement (RCP 2.6), an intermediate pathway (RCP 4.5), and a high-emissions pathway (RCP 8.5). For these pathways, we derive air temperature anomalies from global climate model (GCM) realizations ( Fig. 1A ), which we use to adjust present-day temperatures and precipitation based on climatologies from the high-resolution (∼5.5 km) regional climate model HIRHAM5 ( 16 ). Because most GCM projections are only available until 2300, we extend the temperature anomalies by linearly extrapolating the 2200–2300 trend to the year 2500 and keeping the 2500 monthly values constant afterward (see Materials and Methods). At the ice-ocean interface, we prescribe submarine melt (i.e., melt at the ice front and below ice shelves from the water in contact) with parameterizations informed by observations ( 17 , 18 , 19 ), numerical modeling ( 20 ), and theoretical considerations ( 21 ), and we assume that ocean temperatures rise at the same rate as atmospheric temperatures. To estimate the impact of parametric uncertainties for each of the three emissions scenarios, we perform a rigorous 500-member ensemble of simulations in which we vary 11 key model parameters governing atmospheric forcing, surface processes, submarine melt, calving, and ice dynamics. In addition, we identify a control simulation for each scenario, which we use for inter-scenario comparisons.

As recognized in previous reports of the Intergovernmental Panel on Climate Change, the ability to track outlet glacier behavior is necessary to accurately project ice sheet evolution, yet previous studies ( 10 , 11 ) had limited success due to a lack of accurate ice thickness ( 12 ), a leading order constraint on ice flow. Consequently, major progress was not possible until 2014, when a new high-resolution ice thickness map ( 13 ) became available that now allows capturing outlet glacier flow with high fidelity ( 12 ).

Ice sheets lose mass through runoff of surface meltwater and ice discharge into the surrounding ocean, and increases in both over the past two decades have resulted in accelerated mass loss from the Greenland Ice Sheet ( 1 , 2 ). Between 1995 and 1998, subsurface ocean temperatures rose by about 1.5°C along the west coast of Greenland as a result of increasing subsurface water temperatures in the subpolar gyre ( 3 ). These warm waters led to a disintegration of buttressing floating ice tongues ( 3 ), which triggered a positive feedback between retreat, thinning, and outlet glacier acceleration (outlet glacier–acceleration feedback) ( 4 ). A stark example is Jakobshavn Isbræ, Greenland’s largest and fastest flowing outlet glacier. Following the loss of its floating tongue between 2000 and 2003, the glacier’s flow speeds doubled as the ice thinned ( 5 , 6 ). Since then, rapid increases in ice discharge have been observed in many outlet glaciers around Greenland, including Sverdrup Gletscher and Ummiamaku Isbræ on the northwest coast and Køge Bugt and Kangerdlugssuaq Gletscher on the southeast coast ( 7 ). Between 1990 and 2008, surface melt nearly doubled in magnitude, most notably in southwest Greenland ( 8 ), resulting in additional widespread thinning at lower elevations. Thinning lowers the ice surface and exposes it to higher air temperatures, thereby leading to enhanced melt (surface mass balance–elevation feedback), establishing a second positive feedback for mass loss. The existence of such positive feedbacks can lead to strongly nonlinear responses of the ice sheet to environmental forcings ( 9 ).

RESULTS AND DISCUSSION

Greenland in a thousand years In a thousand years, the Greenland Ice Sheet will look significantly different than today (Fig. 2). Depending on the emission scenario, the Greenland Ice Sheet will have lost 8 to 25% (RCP 2.6), 26 to 57% (RCP 4.5), or 72 to 100% (RCP 8.5) of its present-day mass, contributing 0.59 to 1.88 m, 1.86 to 4.17 m, or 5.23 to 7.28 m to global mean sea level, respectively, where ranges refer to the 16th and 84th percentiles (Fig. 1B and Table 1). We illustrate the range of possible ice sheet trajectories by computing the probability, in our 500-run ensemble, that a given location is ice covered after 1000 years (Fig. 2, B to E). The ice sheet following the RCP 2.6 emissions path will very likely experience significant margin recession in the west and north (Fig. 2B). For the control simulation, ice flow in this ice sheet shows patterns similar to present day in the southeast and southwest (Fig. 2F). RCP 4.5 results in large retreats around all of Greenland, with the exception of high-elevation areas in the east of the southern and central domes (Fig. 2C). Fast flow in this reduced ice sheet configuration is still topographically concentrated in channels that reach below sea level. A wide swath of outlet glaciers feed fjords in west-central Greenland, and the contemporary large outlet glaciers of northwestern and northeastern Greenland merge into a few ice stream–like features (Fig. 2G). For RCP 8.5, 67% of the ensemble members lose >90% of the initial volume, and 58% of ensemble members lose >99% of initial ice volume, including the control simulation (Fig. 2D). Evidently, by continuing on the RCP 8.5 path, it is very likely that Greenland will become ice free within a millennium. Fig. 2 Observed 2008 state and simulations of the Greenland Ice Sheet at year 3000. (A) Observed 2008 ice extent (53). (B to D) Likelihood (percentiles) of ice cover as percentage of the ensemble simulations with nonzero ice thickness. Likelihoods less than the 16th percentile are masked. (E) Multiyear composite of observed surface speeds (61). (F to H) Surface speeds from the control simulation. Basin names shown in (A) in clockwise order are southwest (SW), central-west (CW), northwest (NW), north (NO), northeast (NE), and southeast (SE). RCP 2.6 (B and F), RCP 4.5 (C and G), and RCP 8.5 (D and H). Topography in meters above sea level (m a.s.l.) [(A) to (H)]. Table 1 Contribution to GMSL in centimeters relative to 2008 for the years 2100, 2200, 2300, and 3000. Uncertainties for the ensemble analysis (ENS) at 1.8-km horizontal grid resolution are given as the 16th and 84th percentile range. In addition, the GMSL contribution from the control simulation (CTRL) at 900-m horizontal grid resolution is shown. To study the sensitivity of mass loss to grid resolution, we run additional simulations at 18 km (G18000), 9 km (G9000), 4.5 km (G4500), 3.6 km (G3600), 1.8 km (G1800), and 600 m (G600). NGIA is a simulation without glacio-isostatic adjustment, and NTLR is a simulation without a temperature lapse rate; both simulations were performed at 900 m. We also performed a simulation at 18 km that used the shallow-ice approximation (SIA18000). G600–18000, SIA18000, NGIA, and NTRL use the same parameters as CTRL. View this table:

Partitioning mass loss At the beginning of the 21st century, mass loss was roughly equally split between surface meltwater runoff and ice discharge (sum of mechanical calving and frontal melt) into the surrounding ocean, albeit with high year-to-year variations (22). Our simulated mass loss for the same period is consistent with this observation (Fig. 1). During the 21st and 22nd centuries, ice discharge remains a major contributor to mass loss regardless of the emissions scenario (Figs. 1, D and E, and 3), contributing 2 to 39 cm to sea level by 2200, corresponding to 6 to 45% of the total mass loss. Over time, however, the relative importance of ice discharge diminishes, except for RCP 8.5. Under this scenario, mass loss is sufficiently large for the ice margin to retreat into interior areas below sea level, resulting in large calving fronts and increased ice discharge. The exact timing of this increase varies across RCP 8.5 ensemble members, resulting in a marked increase in the variance of the relative contribution of ice discharge to mass loss at the beginning of the 23rd century. Fig. 3 Evolution of ice sheet area and mass balance components for the control simulation for each RCP scenario. (A to C) Ice area evolution. (D to F) Partitioning of ice sheet wide mass balance rates into snow accumulation, runoff, and ice discharge into the ocean shown in Gt year−1 (D to F) and kg m−2 year−1 (left axis) and m year−1 ice equivalent (right axis) (G to I). We distinguish between runoff due to climate warming and runoff due to surface elevation lowering. (A) to (I) are plotted as 11-year running means. In a warming climate, precipitation (and thus snow accumulation) is expected to increase because of the higher moisture holding ability of warmer air. Here, we increase precipitation by 5 to 7% for each degree of warming (23). We find that in our simulations, total snow accumulation over the ice sheet remains relatively steady over time for RCPs 2.6 and 4.5 (Fig. 3, D and E), since decreasing ice sheet area is offset by an increase in snow accumulation per unit area due to increasing precipitation (Fig. 3, J and K). Under RCP 8.5, however, the decrease in snow accumulation due to accumulation area reduction (Fig. 3C) outpaces the increase in snow accumulation due to warming, and thus, total snow accumulation decreases after around year 2200 (Fig. 3F). Toward the end of the millennium, snow accumulation per unit area increases as the Greenland Ice Sheet is reduced to a few high-elevation areas in the southeast (Fig. 3I). For RCP 8.5, ice sheet–wide surface meltwater runoff rates decrease after passing a maximum around year 2500 despite continuously increasing runoff rates per unit area. At the time of the maximum, runoff exceeds snow accumulation by a factor of about 17. Surface melt is amplified by the positive surface mass balance–elevation feedback as surface lowering exposes the ice to higher air temperatures. To assess the role of this feedback in driving increases in surface melt, we perform a simulation assuming no changes in surface elevation for the calculation of the air temperatures (i.e., setting the temperature lapse rate to zero). We find that in all scenarios, the increased melt rates per unit area caused by higher air temperatures due to climate change exceeds the impact of higher temperatures due to surface lowering (Fig. 3, G to I). In our simulations, temperatures are kept constant beyond the year 2500, yet Greenland continues to lose mass in part because of the surface mass balance–elevation feedback. This committed sea level rise (24) demonstrates that stabilization of mass change will be more difficult as time passes.

Additional feedbacks at play Besides the surface mass balance–elevation feedback and the outlet glacier–thinning feedback, additional negative and positive feedbacks are at play. A negative feedback that can reduce mass loss is glacio-isostatic adjustment that results from unloading of the lithosphere due to ice loss. The uplift of bedrock is caused by the viscoelastic response of the underlying mantle and, as a consequence, will reduce the surface area that is exposed to surface melt and the ice area in direct contact with ocean water. However, because of the high mantle viscosity, this process is relatively slow, causing millimeters to centimeters per year rates of uplift. For each RCP scenario, we perform a simulation without glacio-isostatic adjustment and find that this effect is negligible on the centennial time scale, and after a millennium, mass loss is reduced by only about 2%, which is much less than the variance of the ensemble simulations on the 16th or 84th percentile (Table 1). Another negative feedback is the coastward advection of deep cold ice, which increases ice viscosity and basal stickiness and decreases mass loss by reducing flow rates. This feedback competes with the positive feedback of inland migration of decreased basal stickiness due to outlet glacier acceleration and thinning (fig. S1). Acceleration of outlet glaciers not only leads to enhanced ice discharge but also contributes indirectly to increased surface runoff via the surface mass balance–elevation feedback: Thinning induced by outlet glacier acceleration lowers the ice surface in the vicinity of the glacier terminus, resulting in enhanced melt. While our model takes these three feedbacks into account, their impact is, however, difficult to quantify. Positive feedbacks that we have not considered, but are potentially relevant, are the effect of enhanced surface melt affecting basal motion (25), warming of the ice by refreezing of meltwater [“cryo-hydrologic warming” (26)], and ice cliff failure (27). These three feedbacks, if taken into account, could further increase our mass loss estimates. Last, geomorphological processes such as bedrock erosion, sediment transport, and deposition could affect ice sheet evolution on centennial and longer time scales (28); however, whether these processes would amplify or slow down mass loss is not clear (29).

Outlet glacier retreat Observations indicate that Greenland’s 200+ major outlet glaciers have displayed a marked asynchronicity in retreat behavior (30). Even adjacent glaciers that experience similar environmental conditions may behave differently because retreat is strongly controlled by glacier geometry (31). Our simulations produce similar behavior, as illustrated by two examples. In west-central Greenland, Store Gletscher’s current extent is strongly controlled by geometry because its terminus resides on a local bedrock high, and high frontal melt rates are required to thin the glacier to flotation and cause retreat (32). In our simulations, once Store Gletscher loses contact with the bedrock high, it retreats within a decade or less over a distance of about 25 km until it becomes land terminating (Fig. 4). This behavior is observed in all three RCP scenarios but occurs later under lower emissions scenarios (not shown). Upernavik Isstrøm South, in the same region, exhibits a more gradual retreat in the simulations due to a smoother bed topography. Both glaciers accelerate during retreat and then slow down at stable front positions on bedrock highs while thinning continues. Once thinned to flotation at this local high point, the glaciers retreat again. In our simulations, both Store Gletscher and Upernavik Isstrøm South undergo several episodes of fast retreat. During these episodes, maximum speeds are comparable between individual episodes, but maximum fluxes (product of ice thickness and vertically averaged velocity) decline over time because of thinning, a behavior that was also reported by Nick et al. (33). Fig. 4 Retreat of two outlet glaciers in a similar climatic setting between 2015 and 2315 for the RCP 4.5 scenario. (A) Upernavik Isstrøm South (UIS) shows a gradual retreat of about 50 km over the next 200 years. (B) Store Gletscher (SG) is currently in a very stable position on a bedrock high. It takes almost a hundred years before substantial retreat happens. However, once the glacier loses contact with the bedrock high, retreat of 25 km occurs in less than a decade. The glacier retreats quickly until it is out of the water. Every line represents a year (A and B). (C) Location of the two outlet glaciers on the west coast, present-day observed surface speeds (61), and flow lines of Upernavik Isstrøm South and Store Gletscher (white dashed lines). Small inset shows area where the two glaciers are located. Over time, ice sheet–wide ice discharge decreases because of outlet glacier thinning and outlet glaciers becoming land terminating (Fig. 3). Because the flow of outlet glaciers is strongly controlled by geometry (34) and most of the submarine channels beneath large outlet glaciers extend only to about 100 km inland, their potential for sustained fast retreat (and large ice discharge) is limited. Jakobshavn Isbræ, Humboldt Gletscher, and Petermann Gletscher, however, have channels that extend far into the ice sheet interior (13), and these glaciers are responsible for much of the asymmetric retreat (Fig. 2, F and G). By the year 2300 (RCP 8.5) or 2500 (RCP 4.5), almost all outlet glaciers in northwest Greenland have become land terminating, and ice discharge there is greatly reduced. In contrast, in southeast Greenland, outlet glaciers retreat substantially only for the RCP 8.5 scenario, and ice discharge remains an important contributor to mass loss even for the RCP 8.5 scenario until the 23rd century. However, future ice discharge in the southeast may be underestimated in our simulations because of inaccurate subglacial topography. Aschwanden et al. (12) already reported poor agreement between observed and modeled surface speeds for several glaciers in southeast Greenland and attributed the mismatch to poorly constrained ice thickness. This hypothesis was recently corroborated by inversion of gravimetric data that revealed glacial fjords several hundreds of meters deeper than previously assumed (35). To characterize the importance of capturing outlet glacier flow, we performed a simulation where flow was calculated because of vertical shearing at a coarse horizontal grid resolution of 18 km, ignoring longitudinal stresses relevant for outlet glacier flow. For RCPs 2.6 and 4.5, this approach underestimates mass loss by a meter sea level equivalent at year 3000 compared to the control simulation (Table 1), with a projected mass gain of 25 cm sea level equivalent for RCP 2.6. While under RCP 8.5, Greenland will become ice free whether outlet glaciers are resolved or not, larger mass loss occurs in the early centuries in the control simulation, with the nonresolving simulation underestimating mass loss by >10% by the year 2300. Large mass loss in early centuries may make a recovery harder, even if climate warming were to reverse. Consequently, it is crucial to resolve outlet glaciers to reduce uncertainties in mass loss projections.

Uncertainty quantification Our large ensemble of simulations allows us to attribute the fraction of mass loss uncertainty due to poorly constrained model parameters using Sobol indices (36). The sum of Sobol indices must be less than unity because the combination of variance produced by all parameters simultaneously must be less than the total variance. Across all RCP scenarios, we find that 26 to 53% of mass loss uncertainty in the 21st century is caused by underconstrained ice dynamics parameters, particularly uncertainty in basal motion (Table 2). This percentage declines to 5 to 38% in the 22nd century and steadily decreases to 2 to 33% by the year 2300 (Table 2). Over time, uncertainty in ice dynamics explains a progressively smaller fraction of mass loss variance, while the uncertainty in climate forcing (temperature projections in particular) explains an increasingly large fraction of the total ensemble variance, reaching 7 to 45% by the year 2300. We note that beyond year 2300, the calculated variance is likely to underestimate the true variance because temperature projections are not available and we instead produce temperature projections via extrapolation (see Materials and Methods). Table 2 Sobol indices computed from large ensemble of simulations. Values represent the percentage of variance in mass loss attributable to the variance in a given parameter. Large values imply that uncertainty in that parameter is responsible for a commensurately large uncertainty in mass loss. Small values imply that uncertainty in a given parameter has relatively little effect on uncertainty in total mass loss. Numbers for the variance in air temperature for year 3000 are in parentheses because they do not reflect the GCM intermodel variability but the choice of extrapolation. View this table: Surface processes such as melt and refreezing are the result of a complex surface energy balance. Here, we make the first-order assumption that melt is proportional to the sum of days with temperatures above freezing. While more sophisticated approaches exist, they are not computationally tractable on the long time scales and the number of ensembles considered here; they also suffer from underconstrained parameters (37). We find that uncertainties in surface processes are the dominant source of uncertainties across RCPs until year 2300, ranging from 14 to 50%. Between melt and refreezing, refreezing contributes little to uncertainties in surface processes. The mass loss variance explained by uncertainties in submarine melt and calving parameters is <5% for all scenarios and at all times. We emphasize that this does not necessarily imply that these mechanisms are unimportant for ice sheet evolution. Rather, variability in parameter values over their plausible ranges produces relatively little corresponding variability in simulated mass loss. However, such variability can have large regional impacts, and many of the frontal mass loss parameters are empirical and based on recent observations. Furthermore, basal motion is generally high near marine termini, leading to a tight coupling between ocean and basal processes. Thus, it may be difficult to separate the variance in mass loss explained by basal motion from the variance explained by ocean forcing. Our findings suggest that better constraints on basal motion modeling are critical for reducing uncertainties in prediction of mass loss over the next two centuries, and a reduction in uncertainty of temperature projections and for surface mass balance will lead to better predictions of Greenland mass loss over multicentennial and longer time scales.