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

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.

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.

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

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

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.

Perspectives

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 [T₀(x, y, t)] and precipitation [P₀(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 P₀ and T₀ are periodic in time with a period of 1 year.

To perform future simulations, we adjust T₀(x, y, t) by scalar temperature anomalies ΔT_air(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 2006–2015 monthly means 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

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⁻¹. Precipitation is increased as a function of air temperature increase

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⁻³ 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 ${\dot{m}}_x^o$ using a linear function of latitude L between 71°N and 80°N

where the values of ${\dot{m}}_{71}$ and ${\dot{m}}_{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⁻¹ in 1985 with an increase of 57 ± 12 m year⁻¹ 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⁻¹ (18). Inferred contemporary subshelf melt rates at Nioghalvfjerdsfjorden (79 North Glacier; 79°0′N, 20°0′W) are 13.3 ± 4 m year⁻¹, whereas at Zachariæ Isstrøm (78°0′N, 30°0′W), melt rates increased from 14.6 ± 4.1 m year⁻¹ (1999–2010) to 25 ± 12 m year⁻¹ (2010–2014) (19). To capture the observed spatial and temporal variability and associated uncertainties, we consider three scenarios: “LOW”: $({\dot{m}}_{71},{\dot{m}}_{80})=(\mathrm{300,10})$ m year⁻¹, “MID”: $({\dot{m}}_{71},{\dot{m}}_{80})=(\mathrm{400,20})$ m year⁻¹, and “HIGH”: $({\dot{m}}_{71},{\dot{m}}_{80})=(\mathrm{500,30})$ m year⁻¹.

The parameterization of the temporal variability is inspired by Xu et al. (18), who calculate submarine melt $\dot{m}$ as a function of subglacial discharge Q_sgl and thermal forcing T_h

with 0.5 ≲ α ≲ 1 and 1 ≲ β ≲ 1.6 (18). From Eq. 1, we derive an expression for submarine melt anomalies ${\dot{m}}_t^o(\mathrm{\Delta }{T}_{\text{air}})$ as a function air temperature anomalies ΔT_air(t). First, we assume subglacial discharge Q_sgl to equal surface meltwater runoff Q_r. Second, we express Q_r as function of air temperature anomalies ΔT_air(t). We fit a linear function Q_r = aT_a + b (a = 0.50, b = 0.8, r² = 0.93) (fig. S3C) through mean summer (June-July-August) derived air temperatures T_a and surface meltwater runoff Q_r 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 T_h in terms of ΔT_air. While the ocean warms slower than the atmosphere, we ignore this lag and assume that ocean thermal forcing T_h increases by 1 K for every 1-K increase in air temperature. We thus scale submarine melt by a function of temperature anomalies ΔT_air

We use three scenarios for ${\dot{m}}_t^o$: 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 ΔT_air = 8 K, submarine melt rates increase by a factor of ≈6 to 24 m year⁻¹ to 3000 to 12,000 m year⁻¹, which should be viewed as an upper bound (fig. S4D).

Calving

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 h_min 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 h₇₄ and h₇₆, respectively. To capture the observed spatial and temporal variability and associated uncertainties, we consider three scenarios; LOW: h_min = (h₇₄, h₇₆) = (400,50) m, MID: h_min = (h₇₄, h₇₆) = (500,100) m, and HIGH: h_min = (h₇₄, h₇₆) = (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, η

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 t_b and the sliding velocity u_b

where τ_c is the yield stress, q is the pseudoplasticity exponent, and u₀ = 100 m year⁻¹ is a threshold speed. We assume that yield stress τ_c is proportional to effective pressure N [“Mohr-Coulomb criterion” (50)]

where ϕ is the till friction angle and the effective pressure N is given in (51, 52)

Here, δ = 0.02 is a lower limit of the effective pressure, expressed as a fraction of overburden pressure, e₀ is the void ratio at a reference effective pressure N₀, C_c is the coefficient of compressibility of the sediment, W is the effective thickness of water, and W^max is the maximum amount of basal water. We use a nonconserving hydrology model that connects W to the basal melt rate ${\dot{B}}_b$ (51)

where ρ_w is the density of water and C_d = 1 mm year⁻¹ is a fixed drainage rate. The basal stickiness β is a measure of the glacier bed’s resistance to basal motion and is defined as

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⁻¹ (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 f_i, positive degree-day snow melt factor f_s, and refreezing proportion ψ. Literature values for the ice melt factor f_i range from 8 to 40 mm K⁻¹ day⁻¹ (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⁻¹ day⁻¹ significantly overestimate melt (not shown). We thus use a truncated normal distributions with mean μ = 8 mm K⁻¹ day⁻¹ and SD σ = 4 mm K⁻¹ day⁻¹. For the snow melt factor f_s, we also use a truncated normal distribution with a mean μ = 4.1 mm K⁻¹ day⁻¹ and an SD of 1.5 mm K⁻¹ day⁻¹ (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 ${\dot{m}}_x$ and the temporal ${\dot{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 h_min. 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).