Journal topic
The Cryosphere, 13, 723–734, 2019
https://doi.org/10.5194/tc-13-723-2019
The Cryosphere, 13, 723–734, 2019
https://doi.org/10.5194/tc-13-723-2019

Research article 01 Mar 2019

Research article | 01 Mar 2019

Modeling the response of northwest Greenland to enhanced ocean thermal forcing and subglacial discharge

Modeling the response of northwest Greenland to enhanced ocean thermal forcing and subglacial discharge
Mathieu Morlighem1,*, Michael Wood1, Hélène Seroussi2, Youngmin Choi1, and Eric Rignot1,2 Mathieu Morlighem et al.
• 1Department of Earth System Science, University of California, Irvine, 3218 Croul Hall, Irvine, CA 92697-3100, USA
• 2Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109-8099, USA
• * Invited contribution by Mathieu Morlighem, recipient of the EGU Arne Richter Award for Outstanding Early Career Scientists 2018.

Correspondence: Mathieu Morlighem (mathieu.morlighem@uci.edu)

Abstract

Calving-front dynamics is an important control on Greenland's ice mass balance. Ice front retreat of marine-terminating glaciers may, for example, lead to a loss in resistive stress, which ultimately results in glacier acceleration and thinning. Over the past decade, it has been suggested that such retreats may be triggered by warm and salty Atlantic Water, which is typically found at a depth below 200–300 m. An increase in subglacial water discharge at glacier ice fronts due to enhanced surface runoff may also be responsible for an intensification of undercutting and calving. An increase in ocean thermal forcing or subglacial discharge therefore has the potential to destabilize marine-terminating glaciers along the coast of Greenland. It remains unclear which glaciers are currently stable but may retreat in the future and how far inland and how fast they will retreat. Here, we quantify the sensitivity and vulnerability of marine-terminating glaciers along the northwest coast of Greenland (from 72.5 to 76 N) to ocean forcing and subglacial discharge using the Ice Sheet System Model (ISSM). We rely on a parameterization of undercutting based on ocean thermal forcing and subglacial discharge and use ocean temperature and salinity from high-resolution ECCO2 (Estimating the Circulation and Climate of the Ocean, Phase II) simulations at the fjord mouth to constrain the ocean thermal forcing. The ice flow model includes a calving law based on a tensile von Mises criterion. We find that some glaciers, such as Dietrichson Gletscher or Alison Glacier, are sensitive to small increases in ocean thermal forcing, while others, such as Illullip Sermia or Cornell Gletscher, are remarkably stable, even in a +3C ocean warming scenario. Under the most intense experiment, we find that Hayes Gletscher retreats by more than 50 km inland by 2100 into a deep trough, and its velocity increases by a factor of 3 over only 23 years. The model confirms that ice–ocean interactions can trigger extensive and rapid glacier retreat, but the bed controls the rate and magnitude of the retreat. Under current oceanic and atmospheric conditions, we find that this sector of the Greenland ice sheet alone will contribute more than 1 cm to sea level rise and up to 3 cm by 2100 under the most extreme scenario.

1 Introduction

Over the past 2 decades, many glaciers along the northwest coast of Greenland have been retreating and accelerating, sometimes dramatically (e.g., Moon et al.2012; Wood et al.2018). It has been suggested that the retreat of these glaciers is initiated by the presence of warm and salty subsurface Atlantic Water (AW) in the fjords . This water is typically found 200 to 300 m below the surface (e.g., Rignot et al.2016a; Holland et al.2008). Surface runoff has also been increasing over the past decades , which enhances subglacial water discharge at the base of calving fronts. This freshwater flux enhances the circulation of the ocean in the fjord , which in turn further increases the melting rate and therefore the rate of undercutting at the calving face of marine-terminating glaciers. While we expect both surface runoff and ocean heat content to continue to increase over the next century, it remains unclear how they are going to affect ice dynamics and ice discharge into the ocean.

While geographically close, individual outlet glaciers along the coast respond differently to frontal forcing. It has been proposed (e.g., Wood et al.2018; Catania et al.2018) that this heterogeneity in glacier behavior may be due to differences in bed topography and fjord bathymetry, which may prevent AW from interacting with calving fronts due to the presence of sills in the fjord. It has also been suggested that many glaciers are currently resting on pronounced ridges, or in regions of lateral constrictions, which stabilizes the glaciers' calving fronts and prevents warm water from dislodging them from their current position . The idea that ice front dynamics is, to a large extent, controlled by subglacial topography was first investigated in Alaska and was more recently extended to Greenland . It is not certain to which degree the glaciers of the northwest coast remain sensitive to enhanced thermal forcing from the ocean: some glaciers are on the verge of a fast and extensive retreat, others may continue retreating at the same rate, and some may remain stable. Numerical modeling can help us assess the sensitivity of these individual glaciers to ocean temperature along the coast and their potential for fast retreat and mass loss, which affects sea level rise.

While many model-based studies have focused on the response of the Greenland ice sheet to climate change, they either (i) did not include moving calving fronts or (ii) were based on flow-line models (e.g., Nick et al.2013) that do not capture changes in lateral drag well (since lateral drag is parameterized) or the complex three-dimensional shape of the bed that affects the retreat rate , and (iii) did not consider undercutting. Here, we want to overcome these limitations by using a plan-view model with a moving calving front. The calving-front position is allowed to move and is a function of ice velocity, calving rate, and rate of undercutting. While much progress has been made in terms of capturing ice flow through improved datasets and through the development of new stress balance solvers not based on the Shallow Ice Approximation, calving and undercutting remain areas of active research. We use two existing parameterizations of ocean undercutting and calving . While these parameterizations are approximations and do not include all the physics involved in ice–ocean interactions, they have been tested with reasonable success on several glaciers of Greenland . The objective of this study is not to make projections, as we are not forcing the model with given representative concentration pathway (RCP) scenarios, but to assess the sensitivity of northwestern Greenland using existing parameterizations for iceberg calving and undercutting.

We focus here on the northwest coast of Greenland between 72.5 and 76 N: from Upernavik Isstrøm to Sverdrup Gletscher (Fig. 1). This is one of the regions of Greenland where the bed is remarkably well constrained by ice thickness measurements from NASA's Operation IceBridge mission and where NASA's Oceans Melting Greenland mission has been collecting multibeam bathymetry data in most fjords.

We first describe the numerical model and then run the model to 2100 under different scenarios of increase in ocean thermal forcing and subglacial discharge. We then discuss the implications of these experiments and the model limitations, as well as make recommendations for future model studies.

Figure 1Ocean bathymetry (m, blue color scale) and ice velocity (m a−1Joughin et al.2010) of northwestern Greenland. The white line shows the 2007 ice sheet extent and white crosses indicate the locations of CTD data from NASA's Oceans Melting Greenland campaign that are used to calibrate the thermal forcing.

2 Method and data

2.1 Ice flow model setup

We use the Ice Sheet System Model (ISSM, Larour et al.2012) and initialize the model with conditions similar to 2007, which is the nominal year of the surface digital elevation map used here (Greenland Ice Mapping Project Digital Elevation Model, Howat et al.2014). The ice surface elevation and bed topography are from BedMachine v3 , and we use satellite-derived surface velocities from to invert for basal friction, following . We use the shelfy-stream approximation (SSA, MacAyeal1989) for the ice stress balance. While not accurate in slow moving regions, this model is an excellent approximation for the fast outlet glaciers (i.e., >200 m yr−1) that we are focusing on here, where sliding velocities are significantly larger than deformational velocities (e.g., Rignot and Mouginot2012). We assume a depth-averaged viscosity equivalent to a temperature of −8C, which is consistent with , and we use a linear viscous basal friction law following :

$\begin{array}{}\text{(1)}& {\mathbit{\tau }}_{\mathrm{b}}=-{C}^{\mathrm{2}}N{\mathbit{v}}_{\mathrm{b}},\end{array}$

where τb is the basal friction, vb is the ice basal velocity, C is a friction coefficient that is inverted for using surface velocities, and N is the effective pressure. For simplicity, we assume that N is equal to the ice pressure above hydrostatic equilibrium, as if the subglacial hydrological system was forming a sheet connected to the ocean. The model mesh is comprised of 380 000 elements, and its resolution varies between 100 m near the coast and 1 km inland. The model time step is 1 week.

In order to capture the dynamic motion of the calving front, we rely on the level-set method , where the velocity at which the calving front moves is defined as follows:

$\begin{array}{}\text{(2)}& {\mathbit{v}}_{\mathrm{front}}=\mathbit{v}-\left(c+\stackrel{\mathrm{˙}}{M}\right)\mathbit{n},\end{array}$

where v is the ice horizontal velocity vector, c is the calving rate, $\stackrel{\mathrm{˙}}{M}$ is the rate of undercutting at the calving face, and n is a unit normal vector that points outward from the ice domain. Much research is currently being dedicated to derive parameterizations for c and $\stackrel{\mathrm{˙}}{M}$; here we chose to use two recent parameterizations, which are described below.

2.2 Undercutting parameterization

We rely on the undercutting parameterization from , where the rate of undercutting (in m day−1) at the calving face is assumed to follow

$\begin{array}{}\text{(3)}& \stackrel{\mathrm{˙}}{M}=\left(A\phantom{\rule{0.125em}{0ex}}h\phantom{\rule{0.125em}{0ex}}{q}_{\mathrm{sg}}^{\mathit{\alpha }}+B\right){\stackrel{\mathrm{̃}}{T}}^{\mathit{\beta }},\end{array}$

where h is the water depth at the calving front (in m), $A=\mathrm{3}×{\mathrm{10}}^{-\mathrm{4}}$ mα dayα−1Cβ, α=0.39, B=0.15 m day−1Cβ, and β=1.18. $\stackrel{\mathrm{̃}}{T}$ is the ocean thermal forcing (in C), defined as the difference in temperature between the potential temperature of the ocean and the depth-dependent freezing point of sea water:

$\begin{array}{}\text{(4)}& \stackrel{\mathrm{̃}}{T}=T-{T}_{\mathrm{F}},\end{array}$

where T is the ocean temperature at a given depth, and TF is the temperature of the local freezing point, which is assumed to be a linear function of salinity and pressure, following Eq. (1) of . qsg is the subglacial discharge at the glacier terminus (in m day−1). Both $\stackrel{\mathrm{̃}}{T}$ and qsg are monthly averaged. The coefficients α and β are close to the ones expected from the plume theory but were determined from a high-resolution ocean modeling study. The introduction of B is necessary to account for the presence of melt in the case where there is no subglacial discharge. The dependence on h was determined from model experiments with different depths and seems to reflect an acceleration of the melt plume when it rises from greater depths .

To estimate the subglacial discharge of melt water, qsg, we use the results from the downscaled 1 km regional atmospheric climate model (RACMO) runoff field with the subglacial melt rates from and assume for simplicity that the discharge is uniformly distributed across the calving face. showed that the assumption of uniformly distributed melt only generates a 15 % difference in melt compared to a distributed source of qsg.

The ocean thermal forcing, $\stackrel{\mathrm{̃}}{T}$, is derived from the Estimating the Circulation and Climate of the Ocean, Phase II (ECCO2, 2007–2011) and Phase IV (2007–2015), following the procedure described in . To account for the presence of sills in the fjord, $\stackrel{\mathrm{̃}}{T}$ is depth averaged between sea level and the deepest point for which there is a direct horizontal connection to the fjord mouth. The calculated effective depth assumes a perfectly stratified ocean and decreases as we get closer to the calving front where ocean currents are potentially blocked by the bathymetry. Figure 2 illustrates the effective depth for the case of Sverdrup Gletscher. Note that we define the effective depth over the entire model domain, even under currently ice-covered regions. If the modeled ice front retreats past a high bump, it will be accounted for in the calculation of the thermal forcing and the rate of undercutting will be reduced (see Fig. 3b and c). This undercutting parameterization facilitates the definition of the rate of undercutting everywhere in the model domain, and its magnitude depends on the ice front location. The ice sheet model is forced by the surface mass balance of RACMO 2.3 averaged between 1961 and 1990: the increase in runoff (due to the anomaly applied) is assumed to not affect the surface mass balance but only the rate of undercutting through the parameterization provided by Eq. (3).

Figure 2(a) Effective depth (m) of the fjord of Sverdrup Gletscher. The effective depth decreases as we go from the fjord mouth (x=80 km) to the glacier terminus (x=0 km). (b) Thermal forcing at the fjord's mouth (C) for Sverdrup Gletscher from ECCO2.

Figure 3(a) Bed topography (m), (b) effective depth (m), (c) calculated mean rate of undercutting from 2007 to 2017 (m day−1), and (d) calibrated σmax (kPa).

2.3 Calving parameterization

We assume that the calving rate follows the parameterization proposed by , for which the calving rate is proportional to the tensile von Mises stress:

$\begin{array}{}\text{(5)}& c=‖\mathbit{v}‖\frac{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}{{\mathit{\sigma }}_{max}},\end{array}$

where $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ is the tensile von Mises stress, as defined in , and σmax is a threshold that needs to be calibrated for each basin. This calving law is obviously a simplification that may not capture all modes of calving as it only relies on tensile stresses. It is also assumed here that c and $\stackrel{\mathrm{˙}}{M}$ are independent, which is a simplification but has shown some promising results in real-world applications .

To calibrate the calving threshold, we run the model for 10 years (from 2007 to 2017), using the thermal forcing from ECCO2, and adjust σmax in order to match the extent of Landsat-derived ice front retreat: we try to match the observed retreat rather than the retreat rate from 2007 to 2017 along a central flow line for each glacier. This calving threshold is uniform by basin and held constant through time in all simulations. Another possible approach would be to calibrate σmax during a period of ice front stability. One of the problems with this alternative approach is that stable glaciers generally have their termini on distinct basal features, such as ridges or ledges. The numerical model is also stable for a wide range of σmax under these conditions, as shown in and . The threshold σmax is easier to calibrate for retreating glaciers, as it directly constrains the rate of retreat.

2.4 Experiments

After this calibration phase, we run the model forward, from 2017 to 2100, under different scenarios of ocean forcings and different scenarios of increase in subglacial discharge. analyzed the results of 19 climate models to quantify ocean warming around the coast of Greenland over the coming centuries. They found that western Greenland's subsurface ocean temperature reaches between 0.5 and 4 C, with a mean of 1.5 C by 2100. Coupled Model Intercomparison Project Phase 5 (CMIP5) results suggest similar rates of warming by the end of the century under RCP8.5 (Donald Slater, personal communication, 2018). A 2 C increase is also in line with the global atmospheric temperature rise target of the Paris Agreement. Even though there will be a lag in the response of the ocean to atmospheric warming, we do expect that polar amplification could increase ocean temperature further at high latitudes. We therefore consider here a range in $\stackrel{\mathrm{̃}}{T}$ increase from 0 to +3C.

In terms of subglacial discharge, observations over the past decade have shown that surface melting has increased over the entire Greenland ice sheet . showed that meltwater runoff could be multiplied by a factor of 10 by the end of the century. We therefore multiply the subglacial discharge by a factor of up to 10, starting in 2017.

Overall, we perform 40 experiments here: we increase the ocean thermal forcing, $\stackrel{\mathrm{̃}}{T}$, instantly by increments of 1 up to 3 C and multiply the ocean subglacial discharge by up to a factor of 10. We then run the model forward from 2017 to 2100. The rate of undercutting (Eq. 3) is therefore modified as follows:

$\begin{array}{}\text{(6)}& \stackrel{\mathrm{˙}}{M}=\left(A\phantom{\rule{0.125em}{0ex}}h\phantom{\rule{0.125em}{0ex}}{\left({q}_{\mathrm{sg}}×{q}_{\mathrm{a}}\right)}^{\mathit{\alpha }}+B\right){\left(\stackrel{\mathrm{̃}}{T}+{\stackrel{\mathrm{̃}}{T}}_{\mathrm{a}}\right)}^{\mathit{\beta }},\end{array}$

where the subglacial discharge anomaly factor qa varies from 1 to 10, and the thermal forcing anomaly, ${\stackrel{\mathrm{̃}}{T}}_{\mathrm{a}}$, varies from 0 to 3 C. From 2007 to 2016, we rely on the thermal forcing ($\stackrel{\mathrm{̃}}{T}$) and subglacial discharge (qsg) from ECCO2 and RACMO. For 2017 to 2100, as we do not run a coupled model, we repeat the thermal forcing and subglacial discharge of the year 2016 until the end of the century, with the anomalies described above. While a gradual increase in ocean thermal forcing and subglacial discharge would be more realistic, we want to perform a sensitivity analysis in order to determine the glaciers that are more at risk.

Additionally, we perform a control experiment where the ice front is kept fixed. This control experiment is designed to quantify the impact of including moving boundaries in future simulations.

3 Results

Figure 3d shows the chosen value of the stress threshold over the model domain. For the southern half, we find a stress threshold within 20 % of 1 MPa, which is consistent with what was found in other studies . Over the northern side of the domain, however, the stress threshold has to be decreased to ∼650 kPa in order to match the pattern of retreat. This would suggest that the ice is less resistant to tensile stress, but this is more likely an artifact that is due to our underestimation of the rate of undercutting in this region. noted that the north–south temperature gradient in the ocean model was poorly represented in this region and that the resulting thermal forcing was too cold. The model therefore requires a decrease in the stress threshold, thereby increasing the calving rate, c, in order to capture the correct amount of ice retreat over the past 10 years. We could have kept σmax constant, equal to 1 MPa, and optimized the ocean thermal forcing instead, but the spatial and temporal variability in $\stackrel{\mathrm{̃}}{T}$ makes its calibration difficult. Optimizing a single scalar parameter per glacier is more practical.

Figure 4 shows ice front positions that were manually digitized from Level 1 Landsat imagery, together with modeled ice front positions between 2007 and 2017 for four glaciers along the coast. The first two columns of Table 1 list the observed and modeled retreat for the same time period along a central flow line for the chosen value of the stress threshold. By manually tuning the stress threshold (σmax) for each basin, we are able to match the retreat of the past 10 years for all 17 glaciers for which a change has been documented, except for Ussing Bræer N (Table 1), for which we model a retreat of almost 3 km instead of an advance of 300 m. This inconsistency may be due to errors in the bed topography near the front. We note, however, that under all scenarios this glacier remains remarkably stable at its 3 km retreated position, which coincides with a large bump in the bed topography. Overall, we find that with a unique scalar parameter constant in time for each glacier, the modeled ice front retreat is in very good agreement with observations, which is consistent with . The retreat rate of Dietrichson Gletscher is well captured (Fig. 4a vs. b). While the model overestimates the retreat on the southern side of the fjord, there is nonetheless an overall good agreement between the modeled and observed retreat between 2007 and 2017. The front of Illullip Sermia is remarkably stable in both observations and the model (Fig. 4c and d), as it is currently located on a pronounced sill in the bed topography. The modeled ice front of Upernavik Isstrøm retreats more in the southern half of the fjord than the northern half compared to the observations, but the increase in ice retreat over the past 2 years is captured (Fig. 4e and f). The complex pattern of ice front retreat of Kakivfaat Sermiat is also reproduced with a slight difference in timing (Fig. 4g and h). The 2017 modeled front position is also more retreated than what has been observed, but we find the same strong control of the bed topography in the pattern of retreat.

Table 1Observed and modeled ice front retreat (in km along a centerline) between 2007 and 2017 under current forcing (first two columns) and modeled retreat between 2007 and 2030 and between 2007 and 2100, under different scenarios of ocean forcing with today's qsg for individual glaciers along the northwest coast. A more complete table is provided in the Supplement.

Figure 4Observed (a, c, e, g) and modeled (b, d, f, h) ice front position for Dietrichson Gletscher (a, b), Illulip Sermia (c, d), Upernavik Isstrøm C (e, f), and Kakivfaat Sermiat (g, h) under current conditions ($\stackrel{\mathrm{̃}}{T}$ +0C, qsg×1). Yellow to red lines are annual ice front position from 2007 to 2017, and light blue to pink are the model projections for 2017 to 2100.

If we now look at projections, Tables 1 and S1 in the Supplement list the modeled retreated distance compared to the 2007 position for all 40 experiments along a central flow line, and Fig. 5 shows velocity profiles for the different experiments in 2030. Under today's oceanic conditions ($\stackrel{\mathrm{̃}}{T}$ +0C and qsg×1), Sverdrup Gletscher is predicted to continue to retreat for another 5 km (i.e., 8 km upstream of its 2007 position) by 2030 and yet another 5 km by 2100. Under the strongest scenario (i.e., $\stackrel{\mathrm{̃}}{T}+\mathrm{3}$C and qsg×10), Sverdrup Gletscher retreats by 23 km compared to 2007 by 2030 and remains there until the end of the century. We find that Sverdrup Gletscher has three distinct stable positions: ∼8, 13, and 23 km upstream of the 2007 terminus are the ice front positions that we find for a majority of simulations, and they coincide with clear features in the bed topography. Further south, Dietrichson Gletscher will retreat another 1–3 km under the current thermal forcing and may retreat by up to 55 km by 2100 compared to 2007 if $\stackrel{\mathrm{̃}}{T}$ increases by 1 C or more, or if the subglacial discharge increases by a factor of 8 or more. Again, we find clear common retreated positions, 5, 8, 30, 38, and 55 km upstream of the 2007 position, which coincide with topographic features in the bed. Steenstrup Gletscher remains somewhat stable without further ocean warming but retreats by more than 30 km upstream, where the bed rises above sea level if the ocean temperature warms by 1 C or more, or if the subglacial discharge is doubled. Kjer Gletscher exhibits almost the same behavior for all scenarios: it will continue to retreat another ∼40 km upstream over the coming 2 decades in a region of prograde bed slope and remain stable there. Hayes Gletscher N slightly readvanced over the past 10 years but the model suggests that it will retreat by up to 70 km upstream to where the bed is higher than sea level. Hayes Gletscher would retreat 13 km by 2030, in a marked overdeepening of the bed, and continues to retreat another 17 km to reach a position 30 km upstream of its 2007 position by the end of the simulation. If the thermal forcing increases by 2 or 3 C, the glacier retreats 20 km further inland. The different branches of the unnamed glacier south of Hayes Gletscher also retreat; the northern branch retreats 45 km by 2100 in all scenarios to reach a position where the bed rises above sea level. The middle branch (M) retreats by about 40 km by the end of the century in all cases except if the thermal forcing increases by +3C, in which case its ice front retreats by 64 km by 2100. The southern branch shows a more binary behavior: it retreats another 3–7 km, depending on the warming scenario but for enhanced thermal forcing simulations it may retreat 43 km upstream or even 65 km upstream in the case of a +3C warming in $\stackrel{\mathrm{̃}}{T}$. Alison Gletscher has been retreating by 2.5 km over the past 10 years, and the model projects that by 2030, in all cases, it will retreat another 7–8 km upstream due to the lack of features in the bed topography that may stop the retreat. By 2100, the glacier may retreat another 5 km if the thermal forcing increases by +2C or more.

Figure 5Modeled ice velocities (solid lines) and ice front positions (dashed vertical lines) in 2030 for all 40 scenarios. The black dashed line is the current ice velocity (m) and the x axis shows the distance to the current calving-front position.

Illullip Sermia also has a binary behavior. For the strongest forcing, it retreats by 17–18 km, but in the more conservative scenarios it stays at its current position that coincides with a large bump in the bed topography. Cornell Gletscher is one of the most stable glaciers of the model: under all scenarios, it retreats another kilometer upstream of its 2017 position and remains stable there, except in the case of +3C increase in $\stackrel{\mathrm{̃}}{T}$, for which it could retreat by another ∼10 km.

Ussing Bræer N is the glacier for which we do not capture the advance, but under all scenarios the model projects that it will remain stable 3 km upstream of its current position, where the bed is very shallow. Ussing Bræer has been stable over the past 10 years, and the model suggests that it may retreat by 9 to 15 km if the ocean thermal forcing increases by 2 to 3 C, but the glacier does not retreat even when the subglacial discharge is multiplied by 10 in the case of no additional increase in $\stackrel{\mathrm{̃}}{T}$. Qeqertarsuup Sermia is also one of the stable glaciers of this region: the model marginally retreats and under the strongest forcing (+3C) retreats by about 10 km. Kakivfaat Sermiat, on the other hand, has retreated more than 4 km since 2007. The model suggests that, in all cases, it will retreat another 15 km, where a pronounced feature in the bed topography keeps the ice front stable (Fig. 5). Our simulations suggest that the glacier may reach this position by 2030 and remain stable there. Upernavik Isstrøm N retreats by 4 or 11 km depending on the forcing, by 2100. Upernavik Isstrøm C continues to retreat about 3–6 km upstream of its 2007 position, except in the case of a +3C ocean warming under which it would retreat by 23 km. Finally, Upernavik Isstrøm S would remain stable if the current conditions of qsg and $\stackrel{\mathrm{̃}}{T}$ are maintained but may retreat between 17 and 29 km if the subglacial discharge is multiplied by a factor of 6 or if the thermal forcing increases.

Figure 6 shows the contribution to sea level rise of the entire domain for the 40 different scenarios. In all cases, even under current conditions, our simulations suggest that this region will continue to lose mass. The mass loss is significantly higher than in the control experiment, in which we kept the ice front fixed. We also notice that the spread in mass loss due to temperature change (with a fixed qsg) is significantly larger than the spread in mass loss due to an increase in subglacial discharge (with fixed $\stackrel{\mathrm{̃}}{T}$). Note that we rely here on a 1960–1991 average surface mass balance, and the projections of ice loss do not account for the increase in surface melt. Our simulations are therefore conservative and should not be used as actual projections.

Figure 6Contribution to sea level rise (mm) for all 40 scenarios. The black dashed line is the modeled contribution to sea level with a fixed calving front. All simulations rely on a constant surface mass balance.

4 Discussion

Our simulations suggest that all glaciers of the northwest coast, except for four (Illullip Sermia, Ussing Bræer, Qeqertarsuup Sermia and Upernavik Isstrøm S), will continue to retreat several kilometers inland under today's thermal forcing and subglacial discharge. Under these conditions, we do not find any glacier which advances.

In all scenarios, we find that the rate and extent of ice front retreat is strongly dependent on the bed topography: ice fronts are stable on topographic bumps and prograde bed slopes and unstable on retrograde bed slope, which is consistent with previous studies . This is, for example, illustrated in Fig. 4h, where the ice front jumps from basal bump to basal bump and retreats rapidly in overdeepenings. We find this behavior common to all glaciers in the model domain. There is, however, no “intuitive” way to predict where the glaciers will stabilize without running a model. In most cases, the fjords are not symmetrical or ridges do not go all the way across the fjord walls, which makes it difficult to determine whether the ice front will stabilize or not.

We find that some glaciers, such as Alison Gletscher or Upernavik Isstrøm S, are more sensitive to small increases in ocean thermal forcing, while others, such as Cornell Gletscher or Qeqertarsuup Sermia, are very difficult to destabilize, even under a +3C increase in ocean thermal forcing. On the other hand, we find that Hayes Gletscher retreats more than 30 km inland into a deep trough once it goes past a ridge, and its velocity increases by a factor of 3 over only 23 years, before restabilizing, under all warming scenarios.

We show here that calving dynamics is an important control on the ice sheet mass balance that should not be ignored. It has been driving the recent dynamic thinning of several Greenland outlet glaciers , and our model study shows that it may continue to control the mass balance of Greenland. Figure 6 shows, for example, that in all cases the system loses a significant amount of mass, and this mass loss is not captured by the model that keeps a fixed calving front. Models that keep the ice boundary fixed will consistently underestimate ice sheet mass loss as they do not capture the effect of ocean warming. These conservative projections should therefore be treated with caution and efforts should be made to include moving boundaries in continental-scale simulations of the Greenland ice sheet in order to account for ice–ocean interactions, despite the complexity and high grid resolution needed to resolve moving boundaries (∼1 km, ) of such simulations. It is also important to note that the future evolution of Greenland is strongly influenced by the ocean (through the ocean thermal forcing). It is important not only to force predictive ice sheet models with projections of surface mass balance but also to include projections of ocean thermal forcing at the fjord mouth. There may also be some positive or negative feedbacks between changes in surface mass balance and calving. More surface melt, for example, could enhance calving through hydrofracture, while at the same time reducing the ice thickness at the calving front, hence reducing the stress. Ideally, the community should move towards ice–ocean–climate coupled models to fully understand the processes that control the stability of the ice sheet .

Another interesting aspect of this analysis is that glaciers are more sensitive to an increase of 1–2 C in ocean thermal forcing than in a 5- to 10-fold increase in subglacial discharge. This is actually a result of the parameterization of undercutting used here (Eq. 3), which is itself more sensitive to $\stackrel{\mathrm{̃}}{T}$ than qsg: the parameterization is sublinearly dependent on qsg and above linear in $\stackrel{\mathrm{̃}}{T}$. The effect of surface runoff is also limited to summer months, while the ocean thermal forcing affects the glacier year-round. That being said, we do not account for other effects that surface runoff may have on ice dynamics, such as enhanced damage due to hydrofracture, which may lead to a decrease in the stress threshold σmax. Glaciers might therefore be more prone to retreat as qsg increases than what is captured by the current model.

Among other limitations in this study, no numerical ocean model is included: the thermal forcing is prescribed and dictates the rate of undercutting. Similarly, the calving law does not capture all the modes of calving and requires more validation. This study indeed relies on two parameterizations that drive the response of the model to ocean forcings. It is therefore critical to further validate these parameterizations or develop new ones that include more physics and better capture the transfer of heat from the fjord mouth to the calving face as well as iceberg calving. We also assumed that the subglacial discharge was distributed uniformly across the calving front, but observations show that the majority of discharge is routed to one or more large channel outlets (e.g., Fried et al.2015). Frontal undercutting is therefore not distributed uniformly either, even though numerical experiments suggest that the uncertainty in melt is on the order of 15 % . We have also shown how our results were strongly influenced by the bed topography. While the bed is pretty well constrained in this region , it is not free of error, and we have shown again here how important features in the bed topography are for calving front stability.

More importantly, this study paves the way for a Greenland-wide projection that includes realistic parameterizations of moving boundaries, which will provide more reliable estimates than current models that do not include calving. This work also suggests that development of more accurate parameterization of undercutting and calving should be developed as they control the response of the model and its stability in future scenarios. While this work is a first step in this direction, more validation should be performed on these parameterizations, and future parameterizations of undercutting and calving will make models more reliable.

5 Conclusions

In this study, we modeled the response of the northwest coast of Greenland to enhanced oceanic forcing and subglacial discharge and found that this sector will continue to lose mass over the coming decades, regardless of the scenario adopted. The model confirms that ice–ocean interactions have the potential to trigger extensive glacier retreat over a short amount of time (i.e., decades), but the bed topography controls the magnitude and rate of retreat. Overall, the model showed greater sensitivity to enhanced thermal forcing compared to subglacial discharge but did not account for other effects that runoff may have on ice flow. While more work on validating this parameterization of undercutting and the calving law employed here is needed, we showed that accounting for ice front dynamics can lead to significantly more ice loss than with a fixed calving front. Under the current oceanic and atmospheric conditions, this sector alone will contribute more than 1 cm to sea level rise by the end of this century and up to 3 cm in the worst-case scenario.

Code and data availability
Code and data availability.

The data used in this study are freely available from the National Snow and Ice Data Center or upon request to the authors. ISSM is open source and freely available at https://issm.jpl.nasa.gov/ (last access: 15 November 2018, ).

Supplement
Supplement.

Author contributions
Author contributions.

MM set up the model, designed the experiments, ran the simulations, and wrote the manuscript. MW provided the data required to compute the rate of undercutting; HS and YC assisted in conducting the numerical experiments. All authors participated in the writing of the manuscript.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

This work was performed at the University of California Irvine under a contract with the National Aeronautics and Space Administration, Cryospheric Sciences Program (no. NNX15AD55G), and the National Science Foundation's ARCSS program (no. 1504230). 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. This work would not have been possible without data from the NASA Ocean Melting Greenland EVS-3 mission, and NASA Operation IceBridge. We thank Andy Aschwanden, an anonymous reviewer, and the editor, Andreas Vieli, for their helpful and insightful comments.

Edited by: Andreas Vieli
Reviewed by: Andy Aschwanden and one anonymous referee

References

Aschwanden, A., Fahnestock, M. A., and Truffer, M.: Complex Greenland outlet glacier flow captured, Nat. Commun., 7, 1–8, https://doi.org/10.1038/ncomms10524, 2016. a

Bassis, J. N.: Diverse calving patterns linked to glacier geometry, Nat. Geosci., 6, 833–836, https://doi.org/10.1038/ngeo1887, 2013. a

Bindschadler, R. A., Nowicki, S., Abe-Ouchi, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Saito, F.and Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W. L.: Ice-Sheet Model Sensitivities to Environmental Forcing and Their Use in Projecting Future Sea-Level (The SeaRISE Project), J. Glaciol., 59, 195–224, https://doi.org/10.3189/2013JoG12J125, 2013. a, b

Bondzio, J., Morlighem, M., Seroussi, H., Kleiner, T., Ruckamp, M., Mouginot, J., Moon, T., Larour, E., and Humbert, A.: The mechanisms behind Jakobshavn Isbræ's acceleration and mass loss: A 3-D thermomechanical model study, Geophys. Res. Lett., 44, 6252–6260, https://doi.org/10.1002/2017GL073309, 2017. a, b

Bondzio, J. H., Seroussi, H., Morlighem, M., Kleiner, T., Rückamp, M., Humbert, A., and Larour, E. Y.: Modelling calving front dynamics using a level-set method: application to Jakobshavn Isbræ, West Greenland, The Cryosphere, 10, 497–510, https://doi.org/10.5194/tc-10-497-2016, 2016. a, b

Bondzio, J. H., Morlighem, M., Seroussi, H., Wood, M., and Mouginot, J.: Control of ocean temperature on Jakobshavn Isbræ's present and future mass loss, Geophys. Res. Lett., 45, 12912–12921, https://doi.org/10.1029/2018GL079827, 2018. a, b

Budd, W. F., Keage, P. L., and Blundy, N. A.: Empirical studies of ice sliding, J. Glaciol., 23, 157–170, 1979. a

Carr, J. R., Vieli, A., Stokes, C. R., Jamieson, S. S. R., Palmer, S. J., Christoffersen, P., Dowdeswell, J. A., Nick, F. M., Blankenship, D. D., and Young, D. A.: Basal topographic controls on rapid retreat of Humboldt Glacier, Northern Greenland, J. Glaciol., 61, 137–150, https://doi.org/10.3189/2015JoG14J128, 2015. a, b

Catania, G. A., Stearns, L. A., Sutherland, D. A., Fried, M. J., Bartholomaus, T. C., Morlighem, M., Shroyer, E., and Nash, J.: Geometric Controls on Tidewater Glacier Retreat in Central Western Greenland, J. Geophys. Res.-Earth, 123, 2024–2038, https://doi.org/10.1029/2017JF004499, 2018. a, b, c

Choi, Y., Morlighem, M., Rignot, E., Mouginot, J., and Wood, M.: Modeling the response of Nioghalvfjerdsfjorden and Zachariae Isstrøm glaciers, Greenland, to ocean forcing over the next century, Geophys. Res. Lett., 44, 11071–11079, https://doi.org/10.1002/2017GL075174,2017. a, b, c

Choi, Y., Morlighem, M., Wood, M., and Bondzio, J. H.: Comparison of four calving laws to model Greenland outlet glaciers, The Cryosphere, 12, 3735–3746, https://doi.org/10.5194/tc-12-3735-2018, 2018. a, b, c

Felikson, D., Bartholomaus, T. C., Catania, G. A., Korsgaard, N. J., Kjaer, K. H., Morlighem, M., Noel, B., van den Broeke, M., Stearns, L. A., Shroyer, E. L., Sutherland, D. A., and Nash, J. D.: Inland thinning on the Greenland Ice Sheet controlled by outlet glacier geometry, Nat. Geosci., 10, 366–371, https://doi.org/10.1038/ngeo2934, 2017. a

Fettweis, X., Franco, B., Tedesco, M., van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., and Gallée, H.: Estimating the Greenland ice sheet surface mass balance contribution to future sea level rise using the regional atmospheric climate model MAR, The Cryosphere, 7, 469–489, https://doi.org/10.5194/tc-7-469-2013, 2013a. a

Fettweis, X., Hanna, E., Lang, C., Belleflamme, A., Erpicum, M., and Gallée, H.: Brief communication “Important role of the mid-tropospheric atmospheric circulation in the recent surface melt increase over the Greenland ice sheet”, The Cryosphere, 7, 241–248, https://doi.org/10.5194/tc-7-241-2013, 2013b. a, b

Fried, M. J., Catania, G. A., Bartholomaus, T. C., Duncan, D., Davis, M., Stearns, L. A., Nash, J., Shroyer, E., and Sutherland, D.: Distributed subglacial discharge drives significant submarine melt at a Greenland tidewater glacier, Geophys. Res. Lett., 42, 9328–9336, https://doi.org/10.1002/2015GL065806, 2015. a

Gillet-Chaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model, The Cryosphere, 6, 1561–1576, https://doi.org/10.5194/tc-6-1561-2012, 2012. a, b

Holland, D. M., Thomas, R. H., De Young, B., Ribergaard, M. H., and Lyberth, B.: Acceleration of Jakobshavn Isbrae triggered by warm subsurface ocean waters, Nat. Geosci., 1, 659–664, https://doi.org/10.1038/ngeo316, 2008. a, b

Howat, I. M., Negrete, A., and Smith, B. E.: The Greenland Ice Mapping Project (GIMP) land classification and surface elevation data sets, The Cryosphere, 8, 1509–1518, https://doi.org/10.5194/tc-8-1509-2014, 2014. a

Jenkins, A.: Convection-Driven Melting near the Grounding Lines of Ice Shelves and Tidewater Glaciers, J. Phys. Oceanogr., 41, 2279–2294, https://doi.org/10.1175/JPO-D-11-03.1, 2011. a

Jenkins, A., Dutrieux, P., Jacobs, S., McPhail, S., Perrett, J., Webb, A., and White, D.: Observations beneath Pine Island Glacier in West Antarctica and implications for its retreat, Nat. Geosci., 3, 468–472, 2010. a

Joughin, I., Smith, B. E., Howat, I. M., Scambos, T., and Moon, T.: Greenland flow variability from ice-sheet-wide velocity mapping, J. Glaciol., 56, 416–430, 2010. a, b

Khan, S. A., Kjaer, K. H., Bevis, M., Bamber, J. L., Wahr, J., Kjeldsen, K. K., Bjork, A. A., Korsgaard, N. J., Stearns, L. A., van den Broeke, M. R., Liu, L., Larsen, N. K., and Muresan, I. S.: Sustained mass loss of the northeast Greenland ice sheet triggered by regional warming, Nat. Clim. Change, 4, 292–299, https://doi.org/10.1038/NCLIMATE2161, 2014. a

Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, 1–20, https://doi.org/10.1029/2011JF002140, 2012. a, b

Lüthi, M. P., Vieli, A., Moreau, L., Joughin, I. R., Reisser, M., Small, D., and Stober, M.: A century of geometry and velocity evolution at Eqip Sermia, West Greenland, J. Glaciol., 62, 640–654, https://doi.org/10.1017/jog.2016.38, 2016. a

MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment: Theory and application to Ice Stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, 1989. a

Meier, M. F. and Post, A.: Fast tidewater glaciers, J. Geophys. Res., 92, 9051–9058, 1987. a

Mercer, J. H.: The response of fjord glaciers to changes in the firn limit, J. Glaciol., 3, 850–858, 1961. a

Moon, T., Joughin, I., Smith, B., and Howat, I.: 21st-Century Evolution of Greenland Outlet Glacier Velocities, Science, 336, 576–578, https://doi.org/10.1126/science.1219985, 2012. a

Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Ben Dhia, H., and Aubry, D.: Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica, Geophys. Res. Lett., 37, 1–6, https://doi.org/10.1029/2010GL043853, 2010. a

Morlighem, M., Bondzio, J., Seroussi, H., Rignot, E., Larour, E., Humbert, A., and Rebuffi, S.-A.: Modeling of Store Gletscher's calving dynamics, West Greenland, in response to ocean thermal forcing, Geophys. Res. Lett., 43, 2659–2666, https://doi.org/10.1002/2016GL067695, 2016. a, b, c, d, e, f, g

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., O'Cofaigh, C., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multi-beam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11051–11061, https://doi.org/10.1002/2017GL074954, 2017. a, b, c

Nick, F. M., Vieli, A., Howat, I. M., and Joughin, I.: Large-scale changes in Greenland outlet glacier dynamics triggered at the terminus, Nat. Geosci., 2, 110–114, https://doi.org/10.1038/NGEO394, 2009. a

Nick, F. M., Vieli, A., Andersen, M. L., Joughin, I., Payne, A., Edwards, T. L., Pattyn, F., and van de Wal, R. S. W.: Future sea-level rise from Greenland's main outlet glaciers in a warming climate, Nature, 497, 235–238, https://doi.org/10.1038/nature12068, 2013. a, b

Noël, B., van de Berg, W. J., Machguth, H., Lhermitte, S., Howat, I., Fettweis, X., and van den Broeke, M. R.: A daily, 1 km resolution data set of downscaled Greenland ice sheet surface mass balance (1958–2015), The Cryosphere, 10, 2361–2377, https://doi.org/10.5194/tc-10-2361-2016, 2016. a

Nowicki, S. and Seroussi, H.: Projections of future sea level contributions from the Greenland and Antarctic Ice Sheets: Challenges beyond dynamical ice sheet modeling, Oceanography, 31, 107–117, https://doi.org/10.5670/oceanog.2018.216, 2018. a

Osher, S. and Sethian, J. A.: Fronts Propagating with Curvature-Dependent Speed – Algorithms Based on Hamilton-Jacobi Formulations, J. Comput. Phys., 79, 12–49, 1988. a

Petrovic, J. J.: Review Mechanical properties of ice and snow, J. Mater. Sci., 38, 1–6, https://doi.org/10.1023/A:1021134128038, 2003. a

Rignot, E. and Mouginot, J.: Ice flow in Greenland for the International Polar Year 2008-2009, Geophys. Res. Lett., 39, L11501, https://doi.org/10.1029/2012GL051634, 2012. a

Rignot, E., Fenty, I., Menemenlis, D., and Xu, Y.: Spreading of warm ocean waters around Greenland as a possible cause for glacier acceleration, Ann. Glaciol., 53, 257–266, https://doi.org/10.3189/2012AoG60A136, 2012. a

Rignot, E., Fenty, I., Xu, Y., Cai, C., Velicogna, I., Ó Cofaigh, C., Dowdeswell, J. A., Weinrebe, W., Catania, G., and Duncan, D.: Bathymetry data reveal glaciers vulnerable to ice-ocean interaction in Uummannaq and Vaigat glacial fjords, west Greenland, Geophys. Res. Lett., 43, 2667–2674, 2016a. a

Rignot, E., Xu, Y., Menemenlis, D., Mouginot, J., Scheuchl, B., Li, X., Morlighem, M., Seroussi, H., den Broeke, M. v., Fenty, I., Cai, C., An, L., and de Fleurian, B.: Modeling of ocean-induced ice melt rates of five West Greenland glaciers over the past two decades, Geophys. Res. Lett., 43, 6374–6382, https://doi.org/10.1002/2016GL068784, 2016b. a, b, c, d, e

Seroussi, H., Morlighem, M., Rignot, E., Khazendar, A., Larour, E., and Mouginot, J.: Dependence of century-scale projections of the Greenland ice sheet on its thermal regime, J. Glaciol., 59, 1024–1034, https://doi.org/10.3189/2013JoG13J054, 2013. a, b, c

Straneo, F. and Heimbach, P.: North Atlantic warming and the retreat of Greenland's outlet glaciers, Nature, 504, 36–43, https://doi.org/10.1038/nature12854, 2013. a

Straneo, F., Hamilton, G. S., Sutherland, D. A., Stearns, L. A., Davidson, F., Hammill, M. O., Stenson, G. B., and Rosing-Asvid, A.: Rapid circulation of warm subtropical waters in a major glacial fjord in East Greenland, Nat. Geosci., 3, 182–186, https://doi.org/10.1038/NGEO764, 2010. a

Tedesco, M., Fettweis, X., Mote, T., Wahr, J., Alexander, P., Box, J. E., and Wouters, B.: Evidence and analysis of 2012 Greenland records from spaceborne observations, a regional climate model and reanalysis data, The Cryosphere, 7, 615–630, https://doi.org/10.5194/tc-7-615-2013, 2013. a, b

van den Broeke, M., Bamber, J., Ettema, J., Rignot, E., Schrama, E., van de Berg, W. J., van Meijgaard, E., Velicogna, I., and Wouters, B.: Partitioning Recent Greenland Mass Loss, Science, 326, 984–986, https://doi.org/10.1126/science.1178176, 2009. a, b

Warren, C. R.: Terminal environment, topographic control and fluctuations of West Greenland glaciers, Boreas, 20, 1–15, https://doi.org/10.1111/j.1502-3885.1991.tb00453.x, 1991. a, b

Warren, C. R. and Glasser, N. F.: Contrasting response of south Greenland glaciers to recent climatic change, Arctic Alpine Res., 24, 124–132, 1992.  a

Wood, M., Rignot, E., Fenty, I., Menemenlis, D., Millan, R., Morlighem, M., Mouginot, J., and Seroussi, H.: Ocean-induced melt triggers glacier retreat in Northwest Greenland, Geophys. Res. Lett., 45, 8334–8342, https://doi.org/10.1029/2018GL078024, 2018. a, b, c, d, e

Xu, Y., Rignot, E., Menemenlis, D., and Koppes, M.: Numerical experiments on subaqueous melting of Greenland tidewater glaciers in response to ocean warming and enhanced subglacial discharge, Ann. Glaciol., 53, 229–234, https://doi.org/10.3189/2012AoG60A139, 2012. a, b

Xu, Y., Rignot, E., Fenty, I., Menemenlis, D., and Flexas, M. M.: Subaqueous melting of Store Glacier, west Greenland from three-dimensional, high-resolution numerical modeling and ocean observations, Geophys. Res. Lett., 40, 4648–4653, https://doi.org/10.1002/grl.50825, 2013. a, b

Yin, J., Overpeck, J. T., Griffies, S. M., Hu, A., Russell, J. L., and Stouffer, R. J.: Different magnitudes of projected subsurface ocean warming around Greenland and Antarctica, Nat. Geosci., 4, 524–528, https://doi.org/10.1038/NGEO1189, 2011. a