Research ArticleCLIMATOLOGY

Contribution of the Greenland Ice Sheet to sea level over the next millennium

See allHide authors and affiliations

Science Advances  19 Jun 2019:
Vol. 5, no. 6, eaav9396
DOI: 10.1126/sciadv.aav9396


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.


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).

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).

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.

Fig. 1 Time series of air temperature anomalies, cumulative contribution to GMSL since 2008, and rate of GMSL rise due to mass changes of the Greenland Ice Sheet.

(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.


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.


Our simulations suggest that by following the RCP 8.5 scenario, the Greenland Ice Sheet will disappear within a millennium and that the contribution of discharge from outlet glaciers will remain a key component of mass loss over the next centuries. Better quantification of feedbacks that interlink and amplify mass loss processes will help in understanding the extent to which these feedbacks make changes in the Greenland Ice Sheet nonreversible. Because of these feedbacks, delaying mitigation of greenhouse gas emissions is likely to increase sea level rise, even with future significant reductions of emissions.


Ice sheet model

Simulations are performed with the open-source PISM (14), which is thermomechanically coupled and polythermal (38) and uses a hybrid stress balance (39), making PISM well suited for Greenland (12). Horizontal grid resolution ranges from 450 m to 18 km, with the control simulations, and the ensemble simulations using resolutions of 900 and 1800 m, respectively. At the basal boundary, geothermal flux varies in space (40), and a pseudoplastic power law relates bed-parallel shear stress and sliding; details are given in (12). Compared to (12), subglacial topography was updated to version 3; submarine melting at vertical ice fronts was implemented; and a new stress-based calving law suitable for outlet glaciers was implemented (32). To track the bedrock response to a changing ice load, PISM uses the model of the viscous half-space overlain by an elastic plate lithosphere (41).

Atmospheric forcing

The spatial and seasonal distributions of 2-m air temperature [T0(x, y, t)] and precipitation [P0(x, y, t)] to force PISM are provided by the regional climate model HIRHAM5 (16) at a resolution of 0.05 (∼5.55 km), which was forced at the lateral boundaries using the European Centre for Medium-Range Weather Forecasts ERA-Interim reanalysis product (42). We use fields of daily mean temperature and total precipitation averaged over 2000–2015 (43). Here, both P0 and T0 are periodic in time with a period of 1 year.

To perform future simulations, we adjust T0(x, y, t) by scalar temperature anomalies ΔTair(t) derived from a GCM from the fifth phase of the Coupled Model Intercomparison Project (CMIP5). As most CMIP5 simulations are only available until 2100, we select the four GCMs that extend to year 2300: GISS-E2-H, GISS-E2-R, IPSL-CM5A-LR, and MPI-ESM-LR. To calculate ΔTair(t), we first extract monthly near-surface temperature averaged over the domain containing Greenland (12°W to 73°W and 59°N to 84°N). We then extend the resulting temperature series beyond 2300 by linearly extrapolating the 2200–2300 trend to the year 2500, after which we keep the 2500 monthly values constant. Our choice of linear trend extrapolation is supported by GCM simulations with Goddard Institute for Space Studies (GISS) ModelE2 (44). Last, we subtract the 2006–2015 monthly means from the time series to create anomalies ΔTair (fig. S2). The resulting temperature anomalies fall within the ±1σ of the CMIP5 ensemble, which consists of 27 (RCP 2.6), 38 (RCP 4.5), and 40 (RCP 8.5) model realizations. Future air temperatures T(x, y, t) are then computed byT(x,y,t)=T0(x,y,t)+ΔTΓ(x,y,t)+ΔTair(t)where ΔTΓ(x, y, t) is the temperature adjustment due to the differences between the time-evolving modeled ice surface elevation and the fixed topography used by HIRHAM5 via a standard atmospheric lapse rate Γ of 6 K km−1. Precipitation is increased as a function of air temperature increaseP(x,y,t)=P0(x,y,t)×exp(cΔTair(t))where the value of c corresponds to the change of 5 to 7% per Kelvin (23).

Climatic mass balance

The ice flow model requires climatic mass balance (i.e., the balance of accumulation, melt, and refreezing) as a surface boundary condition for mass conservation and temperature for conservation of energy. Accumulation is computed from precipitation and 2-m air temperature; precipitation falls as snow at temperatures below 0°C and as rain at temperatures above 2°C, with a linear transition in between. Surface melt is computed by a temperature index model (45), which uses 2-m air temperature as input. Melt factors for snow are used if snow or firn is exposed at the surface, and melt factors for ice are used otherwise.

We use a simple firn model that only allows removal of firn. At the beginning of the simulation, the firn layer thickness is set to the average depth of the 750 kg m−3 isopycnal calculated by the HIRHAM5 for the period 2000–2015 after an offline spin-up of the snow and firn pack for a 100-year period using the downscaled 1979 ERA-Interim atmospheric forcing (fig. S3A) (43). Refreezing and retention of liquid water in the firn model is based on the control run parametrization given in (43), where the density and cold content of the snow layers determine the refreezing capacity. Where firn is present, the firn layer thickness is reduced by the amount of melted firn.

Submarine melt

We define submarine melt as the sum of melt occurring at front of marine terminating outlet glaciers and melt below ice shelves. Frontal melt is primarily controlled by subglacial discharge and thermal forcing, as supported by observations (18, 46), numerical modeling (20), and theoretical considerations (21). We parameterize the present-day submarine melt rate ṁxo using a linear function of latitude L between 71°N and 80°Nṁxo(L)={ṁ71,L<71°ṁ71+19(L71)(ṁ80ṁ71),71°L80°,ṁ80,L>80°where the values of ṁ71 and ṁ80 were chosen using the following observations. At Jakobshavn Isbræ (69°10′N, 49°50′W), Motyka et al. (17) estimated an area-averaged subshelf melt rate of 228 ± 49 m year−1 in 1985 with an increase of 57 ± 12 m year−1 afterward; however, melt rates close to the grounding line are considerably higher than area averages. August 2010 frontal melt rates at Store Glacier (70°22′N, 50°8′W) were 1100 ± 365 m day−1 (18). Inferred contemporary subshelf melt rates at Nioghalvfjerdsfjorden (79 North Glacier; 79°0′N, 20°0′W) are 13.3 ± 4 m year−1, whereas at Zachariæ Isstrøm (78°0′N, 30°0′W), melt rates increased from 14.6 ± 4.1 m year−1 (1999–2010) to 25 ± 12 m year−1 (2010–2014) (19). To capture the observed spatial and temporal variability and associated uncertainties, we consider three scenarios: “LOW”: (ṁ71,ṁ80)=(300,10) m year−1, “MID”: (ṁ71,ṁ80)=(400,20) m year−1, and “HIGH”: (ṁ71,ṁ80)=(500,30) m year−1.

The parameterization of the temporal variability is inspired by Xu et al. (18), who calculate submarine melt ṁ as a function of subglacial discharge Qsgl and thermal forcing Thṁ=(AQsglα+B)Thβ(1)with 0.5 ≲ α ≲ 1 and 1 ≲ β ≲ 1.6 (18). From Eq. 1, we derive an expression for submarine melt anomalies ṁto(ΔTair) as a function air temperature anomalies ΔTair(t). First, we assume subglacial discharge Qsgl to equal surface meltwater runoff Qr. Second, we express Qr as function of air temperature anomalies ΔTair(t). We fit a linear function Qr = aTa + b (a = 0.50, b = 0.8, r2 = 0.93) (fig. S3C) through mean summer (June-July-August) derived air temperatures Ta and surface meltwater runoff Qr extracted from simulations with the regional climate model MARv3.5.2 2006–2100 forced with CanESM2, MIROC5, and NorESM2, each with the RCP 4.5 and 8.5 scenarios (47). Last, we express the thermal forcing Th in terms of ΔTair. While the ocean warms slower than the atmosphere, we ignore this lag and assume that ocean thermal forcing Th increases by 1 K for every 1-K increase in air temperature. We thus scale submarine melt by a function of temperature anomalies ΔTairṁ(x,y,t)=ṁxo(L(x,y))ṁto(ΔTair(t)),ṁto(ΔTair)=1+C(aΔTair)αΔTairβ

We use three scenarios for m˙to: LOW uses α = 0.5 and β = 1.0, MID uses α = 0.54 and β = 1.17, and HIGH uses α = 0.85 and β = 1.61. We tuned C = 0.15 so that for a warming of ΔTair = 8 K, submarine melt rates increase by a factor of ≈6 to 24 m year−1 to 3000 to 12,000 m year−1, which should be viewed as an upper bound (fig. S4D).


Calving, using a von Mises yield criterion, occurs when the maximum tensile stress exceeds σmax (32). In addition, calving occurs below a minimum floating ice thickness, motivated by the following observations. In 2008, only six perennial ice shelves remained [Petermann Gletscher, Steensby Gletscher, Ryder Gletscher, Nioghalvfjerdsfjorden (79 North), Zachariæ Isstrøm, and Storstrømmen], all north of 76°N. Minimum shelf thickness near the calving front of these shelves ranges from 60 to 100 m (13). In contrast, Jakobshavn Isbræ’s (69°10′N, 49°50′W) floating tongue had a minimum thickness of ∼600 m at the calving front (17). Guided by these observations, we prescribe the minimum floating ice thickness hmin as a linear function of latitude between 74°N and 76°N. South of 74°N and north of 76°N, values are kept constant at h74 and h76, respectively. To capture the observed spatial and temporal variability and associated uncertainties, we consider three scenarios; LOW: hmin = (h74, h76) = (400,50) m, MID: hmin = (h74, h76) = (500,100) m, and HIGH: hmin = (h74, h76) = (600,150) m.

Ice dynamics

Aschwanden et al. (12) performed a thorough calibration of parameters related to ice dynamics, such that simulated surface speeds agree well with observed winter 2008 speeds (48), and we adopt these parameters. The two key parameters here are the Shallow Ice Approximation enhancement factor and the exponent of the sliding law. The enhancement factor E appears in the effective viscosity of glacier ice, η2η=(E A)1(τe2+ε2)1n2n(2)where τe is the effective stress, A is the enthalpy-dependent rate factor, and n is the exponent of the power law. The small constant ε (units of stress) regularizes the flow law at low effective stress, avoiding the problem of infinite viscosity at zero deviatoric stress. The pseudoplastic power law (49) relates bed-parallel shear stress tb and the sliding velocity ubtb=τc ubub(1q)u0q(3)where τc is the yield stress, q is the pseudoplasticity exponent, and u0 = 100 m year−1 is a threshold speed. We assume that yield stress τc is proportional to effective pressure N [“Mohr-Coulomb criterion” (50)]τc=(tanϕ) N(4)where ϕ is the till friction angle and the effective pressure N is given in (51, 52)N=δPo10(e0/Cc)(1(W/Wmax))(5)

Here, δ = 0.02 is a lower limit of the effective pressure, expressed as a fraction of overburden pressure, e0 is the void ratio at a reference effective pressure N0, Cc is the coefficient of compressibility of the sediment, W is the effective thickness of water, and Wmax is the maximum amount of basal water. We use a nonconserving hydrology model that connects W to the basal melt rate B˙b (51)Wt=B˙bρwCd(6)where ρw is the density of water and Cd = 1 mm year−1 is a fixed drainage rate. The basal stickiness β is a measure of the glacier bed’s resistance to basal motion and is defined asβ=τcub(1q)u0q(7)

Ice sheet initial conditions

Ice sheet initial conditions are provided by the “calibrated” experiments in (12), which includes a 125-ka paleoclimate simulation followed by calibration of ice dynamical parameters to minimized differences between observed (48) and simulated surface velocities. Using this as a starting point, we then ran a 100-year long relaxation simulation to account for differences/updates in model physics, but we kept the ice surface close to observations using a flux correction (12). The result was an initial state that is both close to the observed geometry (53) and surface speeds (48) of 2008. On the time scale considered here, predictions of ice sheet behavior are similar to climate forecasts where initial conditions have little impact on the long-term evolution (54). The bedrock deformation model is initialized with present-day subglacial topography (13) and viscous uplift rates (fig. S3B) (55).

Ensemble analysis

Many physical processes governing ice sheet evolution are not precisely constrained. To assess the sensitivity of volume change estimates to these uncertainties, we performed a large ensemble of simulations, with 11 key parameters drawn from a priori probability distributions of plausible values (table S2). Instead of probing such a large parameter space, for each RCP scenario, we produced 500 parameter combinations using Latin hypercube sampling (56), which ensures that the ensemble is optimally representative of the joint distribution of parameters, while gaining efficiency relative to factorial or Monte Carlo designs by eliminating potential redundancies. We vary parameters drawn from five broadly defined groups: (i) Climate: We select the four GCM projections with equal probability and model the percentage increase in precipitation per 1 K increase in air temperature, ω, as a uniform distribution bounded by 5 and 7% K−1 (23). (ii) Surface processes: Parameters governing the interaction between the climate and the ice sheet surface mass balance include a positive degree-day ice melt factor fi, positive degree-day snow melt factor fs, and refreezing proportion ψ. Literature values for the ice melt factor fi range from 8 to 40 mm K−1 day−1 (57, 58); however, a comparison with the mean 2000–2015 surface mass balance simulated with HIRHAM5’s energy balance model reveals that values ≫8 mm K−1 day−1 significantly overestimate melt (not shown). We thus use a truncated normal distributions with mean μ = 8 mm K−1 day−1 and SD σ = 4 mm K−1 day−1. For the snow melt factor fs, we also use a truncated normal distribution with a mean μ = 4.1 mm K−1 day−1 and an SD of 1.5 mm K−1 day−1 (59). The amount of annual snow fall that is allowed to refreeze ψ is also modeled as a truncated normal, but with a mean μ = 50% and SD σ = 20% to cover reported literature values (60). (iii) Ocean: We specify three suites of parameters corresponding to low, moderate, and high regimes for both the spatial m˙x and the temporal m˙t variability in submarine melt. These parameter suites are sampled with equal probability. We specify three suites of parameters corresponding to low, moderate, and high regimes of minimum shelf thickness hmin. These three scenarios are selected with equal probability. The maximum tensile stress of a marine terminus σmax is modeled as a normal distribution symmetrically truncated above 1.4 MPa and below 0.7 MPa (32). (iv) Ice dynamics: Aschwanden et al. (12) performed a calibration of ice dynamical parameters and found values of q = 0.6 for the exponent of the sliding law and E = 1.25 for the shallow ice enhancement factor to result in good agreement with observed flow speeds. Here, we use a truncated normal distribution for q with μ = 0.6 and σ = 0.2, whereas E is modeled with a γ distribution centered around 1.25.

To attribute uncertainties in mass loss to parametric uncertainties, we compute main-effect Sobol indices for each input variable (36). Main-effect Sobol indices can be interpreted as the fraction of output variance that is explained by the variance in a given input parameter, neglecting second-order effects due to parameter interactions. Thus, the Sobol indices produced for a given output (in our case total mass change) always sum to less than unity. Typically, a large fraction of unity signifies that interactions between parameter uncertainty are of second-order importance.

Additional simulations

We also identified an optimal (control) simulation (CTRL) that best reproduces the 2000–2015 mean surface mass balance calculated by HIRHAM5. To characterize model behavior, we performed additional simulations at horizontal grid resolutions ranging from 450 m (G450) to 18 km (G18000) and found that outlet glacier behavior is well captured at horizontal grid resolutions less than 2 km (fig. S4), in agreement with (12). The simulation at 450 m was only performed for RCP 4.5 from 2008 until 2200 because of large file sizes. We also performed two additional simulations at a resolution of 900 m, one without glacio-isostatic adjustment (“NGIA”) and one without a lapse rate (“NTLR”). We used the surface meltwater runoff from the NTLR simulation to approximate the meltwater runoff due to anomaly air temperature warming. To characterize the importance of resolving outlet glacier flow, we performed a simulation that calculates flow due to horizontal shearing only at a horizontal grid resolution of 18 km (SIA18000).


Supplementary material for this article is available at

Fig. S1. Coastward migration of basal cold ice competing with the inland migration of outlet glacier acceleration and thinning.

Fig. S2. Time series of temperature anomalies for the four GCMs that extend until 2300.

Fig. S3. Initial forcing and boundary conditions.

Fig. S4. Ice discharge as a function of horizontal grid resolution.

Table S1. Partitioning of mass fluxes for the control simulation.

Table S2. Parameters and their distributions used in the ensemble analysis, source for distributions, and values for the control simulation CTRL.

Movie S1. Evolution of the Greenland Ice Sheet over the next millennium.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank T. Moon for comments on an earlier version of the manuscript, X. Fettweis for providing MAR simulations, and P. Langen for assistance with the HIRHAM snowpack model. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and by the University of Alaska’s Research Computing Systems (RCS). Open-source software was used at all stages, in particular NumPy (, CDO (, matplotlib (, QGIS3 (, and TimeManager ( Funding: A.A. was supported by NASA grant NNX16AQ40G and NSF grant PLR-1603799; A.A., M.A.F., and D.J.B. were supported by NASA grant NNX17AG65G; R.H. was supported by NSF grant PLR11603815 and NASA grant NSSC17K0566; R.M. was supported by the European Research Council under the European Community Seventh Framework Programme (FP7/ 2007-2013)/ERC grant agreement 610055 as part of the ice2ice project; and S.A.K. was funded, in part, by Carlsbergfondet and the Danish Council for Independent Research. Author contributions: A.A., M.A.F., and M.T. jointly developed the experimental design with contributions from R.H. D.J.B. designed and performed the uncertainty analysis, and A.A. performed and analyzed the simulations. C.K. implemented all project-related code changes. R.M. processed firn depth and climate forcing from HIRHAM5, and S.A.K. provided uplift rates. All authors jointly wrote the manuscript and contributed with interpretation. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. The PISM code is available at The model simulations are archived at ( and HIRHAM5 RCM output at All datasets used in this study are publicly available except the uplift rates, which may be requested from S.A.K.

Stay Connected to Science Advances

Navigate This Article