Breaking News
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

## INTRODUCTION

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 ˙$

where

$D 19$

and

$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

## RESULTS AND DISCUSSION

### Greenland in a thousand years

In a thousand years, the Greenland Ice Sheet will look significantly different than today (Fig. 2). Depending on the emission scenario, the Greenland Ice Sheet will have lost 8 to 25% (RCP 2.6), 26 to 57% (RCP 4.5), or 72 to 100% (RCP 8.5) of its present-day mass, contributing 0.59 to 1.88 m, 1.86 to 4.17 m, or 5.23 to 7.28 m to global mean sea level, respectively, where ranges refer to the 16th and 84th percentiles (Fig. 1B and Table 1). We illustrate the range of possible ice sheet trajectories by computing the probability, in our 500-run ensemble, that a given location is ice covered after 1000 years (Fig. 2, B to E). The ice sheet following the RCP 2.6 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.

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

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

### Perspectives

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.

## MATERIALS AND METHODS

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