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

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


Ice sheets lose mass through runoff of surface meltwater and ice discharge into the surrounding ocean, and increase 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 ( 4 ]). A strong 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 Bay and Kangerdlugssuaq Glacier on the southeast coast ( 7 ). Between 1

990 and 2008, surface melt nearly doubled in magnitude, most notably in southwest Greenland ( 8 ), resulting in additional widespread thinning at lower elevations. Thinning the ice surface and exhibiting it to higher air temperatures, leading to enhanced melting surface mass (elevation feedback), establishing a second positive feedback for mass loss. The existence of such positive feedbacks can lead to strongly nonlinear responses to the ice sheet to environmental forcings ( 9] )

As recognized in previous reports by the Intergovernmental Panel on Climate Change, the ability to track outlet glacier behavior is necessary to accurately project ice sheet evolution, yet previous studies 11 ) had limited success due to a lack of accurate ice thickness ( 12 ), a leading order constraint on ice flow. Therefore, 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 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) emission scenarios ( 15 ): a pathway with reduced greenhouse gas emissions that aligns with the goals of the 2015 Paris Climate Agreement (RCP 2.6), an intermediate pathway ( RCP 4.5), and a high-emission 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 (.55.5 km) regional climate model HIRHAM5 ( 16 ). Because most GCM projections are only available until 2300, we extend the temperature anomalies at 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 the submarine (ie, melt the ice and the ice shelves from the water in contact) with parameterizations informed by observations ( 17 18 , ), and theoretical considerations ( 21 ), and we assume that ocean temperatures rise to the same rate as atmospheric temperatures. To estimate the impact of parametric uncertainties for each of the three emission scenarios, we perform a rigorous 500-member ensemble or 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. 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 lines (thick lines) of RCP 2.6, 4.5, and 8.5 temperature anomalies with respect to calculations 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 the 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 ˙

. () Contribution of ice discharge to

M ˙

given as

D 19 M ˙ = D ˙ / ( D ˙ R ˙; [19659028]) M ˙


D 19


R ˙

are ice discharge rate and surface runoff rate, respectively. ( E] ) Ratio of ice discharge rate and 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 will very likely experience significant margin recession in the west and north (Fig. 2B). For the control simulation, ice flow 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

A ) Observed 2008 ice extent ( 53 ) 2 Observed 2008 state and simulations of the Greenland Ice Sheet at year 3000. ( B to D ) Likelihood (percentiles) or ice cover as percentage of the ensemble simulations with nonzero ice thickness. Likelihoods are less than the 16th percentile are masked. ( E] ) Multiyear composite of observed surface speeds ( 61 ). ( F to H] ) Surface speeds from the control simulation. Basic 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 asl) [(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) that 1.8-km horizontal grid resolution is given as the 16th and 84th percentile range. In addition, the GMSL contribution from the control simulation (CTRL) is shown at 900-m horizontal grid resolution. 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 glacial 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

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 melting) 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 emission 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. In this scenario, mass loss is sufficiently large for the 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 market increase in the variance of the relative contribution of ice discharge to mass loss at the beginning of the 23rd century.

A A A A A A A A A A. ( D to F ) Partition 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. (A) to (I) are plotted as 11-year running means.

In a warm 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). During RCP 8.5, 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). Towards 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 exhibits the ice to higher air temperatures. To assess the role of this feedback in driving increases in surface melt, we assume 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 because of the surface mass balance elevation feedback. This committed sea level rise ( 24 ) demonstrates that stabilization of mass change will be more difficult than time passes.

Additional feedbacks on play

Besides the surface mass balance – elevation feedback and the outlet glacier – retrieval 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 basic stickiness and decreases mass loss at reducing flow rates. This feedback complies 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 indirectly 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 do not consider, but are potentially relevant, are the effect of enhanced surface melt affecting basal motion (). ), and ice cliff failure ()]and ice cliff failure () [), and ice cliff failure. 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); ( 29 ).

Outlet glacier retreat

Observations indicate that Greenland's 200+ major exhibit a market asynchronicity in retreat behavior ( 30 ). Even adjacent glaciers that experience similar environmental conditions may differ 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 glacier, 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 in lower emission scenarios (not shown). Upernavik Isstrom 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 to 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 Great Glacier and Upernavik Isstrom South undergo several episodes of fast retreat. During these episodes, maximum speeds are comparable between individual episodes, but maximum fluxes (product or 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] ) Large Glacier (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 or 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 coastal glaciers on the west coast, present-day observed surface speeds ( 61 ), and flow lines of Upernavik Isstrøm South and Store Glacier (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 of 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 into the ice sheet interior ( 13 ), and these glaciers are responsible for much of the asymmetric retreat (Fig. 2, F and G). At 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 essentially 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. Ash walls 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 (19459006).

To characterize the importance of capturing outlet glacier flow, we performed in simulation where flow was calculated because of vertical shearing at a coarse horizontal grid resolution of 18 km, ignoring longitudinal stresses relevant to outlet glacier flow. For RCPs 2.6 and 4.5, this approach underestimates mass loss by a meter sea level equivalent to year 3000 compared to the control simulation (Table 1), with a projected mass gain or 25 cm sea level equivalent for RCP 2.6. While under RCP 8.5, Greenland will be free of charge if 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.

Uncertainty quantification

Our large ensemble of simulations allow 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 variables produced by all parameters 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 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 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.

Surface processes such as melt and refreezing are the result of a complex surface energy balance. Here, we make the first-order assumption that 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 melting and refreezing, little to uncertainties are contributed to 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 mean 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, basic motion is generally high near marine termini, leading to 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 basic 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 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 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. Due to these feedbacks, delaying mitigation of greenhouse gas emissions is likely to increase sea level rise, even with future significant reductions or 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 m, with the control simulations, and the ensemble simulations using resolutions of 900 and 1800 m, respectively. At basal boundary, geothermal flux varies in space ( 40 ), and pseudoplastic power law refers to 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 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 overlay by an elastic plate lithosphere ( 41 ).

Atmospheric forcing

The spatial and seasonal distributions or 2-m air temperature [ T 0

( x y t] )] and precipitation [] P 0 ( x y t )] to force PISM are provided by the regional climate model HIRHAM5 ( 16 ) at a resolution of 0.05 (.55.55 km), which was forced at the lateral boundaries using the European Center 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 P 0 and T 0 are periodic in time with a period of 1 year.

To perform future simulations, we adjust 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 Δ T air ( 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 table monthly from the time series to create anomalies Δ T air (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 by

T ( x y t ) = T 0 ( x y [19659009] t ) Δ T Γ ( x y t ) + t t t 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 Γ or 6 K km −1 . Precipitation is increased as a function of air temperature increase

P ( x ] y t ) = P 0 ( x y t ) x exp ( c [19659009]⋅Δ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 (ie, 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 (1846), numerical modeling (20), and theoretical considerations (21). We parameterize the present-day submarine melt rate


using a linear function of latitude L between 71°N and 80°N


where the values of




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−1whereas 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”:


m year−1“MID”:


m year−1and “HIGH”:


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


(1)with 0.5 ≲ α ≲ 1 and 1 ≲ β ≲ 1.6 (18). From Eq. 1, we derive an expression for submarine melt anomalies


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



We use three scenarios for


: 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−1which should be viewed as an upper bound (fig. S4D).

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 ub

tb=τ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-Coulombcriterion”(50)]

τc=(tanϕ) N

(4)where ϕ is the till friction angle and the effective pressure N is given in (5152)



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 N0Cc 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




(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



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 fipositive degree-day snow melt factor fsand refreezing proportion ψ. Literature values for the ice melt factor fi range from 8 to 40 mm K−1 day−1 (5758); 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 fswe 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


and the temporal


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

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 (www.numpy.org), CDO (https://code.mpimet.mpg.de/projects/cdo/), matplotlib (https://matplotlib.org/), QGIS3 (https://qgis.org/), and TimeManager (https://github.com/anitagraser/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 https://github.com/pism/pism/tree/v1.0-millennium-study. The model simulations are archived at https://arcticdata.io (https://doi.org/10.18739/A2Z60C21V) and HIRHAM5 RCM output at http://polarportal.dk/en/groenlands-indlandsis/nbsp/links/. All datasets used in this study are publicly available except the uplift rates, which may be requested from S.A.K.

Source link