Velocity response of Petermann Glacier, northwest Greenland, to past and future calving events

Dynamic ice discharge from outlet glaciers across the Greenland Ice Sheet has increased since the beginning of the 21st century. Calving from floating ice tongues that buttress these outlets can accelerate ice flow and discharge of grounded ice. However, little is known about the dynamic impact of ice tongue loss in Greenland compared to ice shelf collapse in Antarctica. The rapidly flowing (∼1000 m a⁻¹) Petermann Glacier in northwest Greenland has one of the ice sheet's last remaining ice tongues, but it lost ∼50 %–60 % (∼40 km in length) of this tongue via two large calving events in 2010 and 2012. The glacier showed a limited velocity response to these calving events, but it is unclear how sensitive it is to future ice tongue loss. Here, we use an ice flow model (Úa) to assess the instantaneous velocity response of Petermann Glacier to past and future calving events. Our results confirm that the glacier was dynamically insensitive to large calving events in 2010 and 2012 (<10 % annual acceleration). We then simulate the future loss of similarly sized sections to the 2012 calving event (∼8 km long) of the ice tongue back to the grounding line. We conclude that thin, soft sections of the ice tongue >12 km away from the grounding line provide little frontal buttressing, and removing them is unlikely to significantly increase ice velocity or discharge. However, once calving removes ice within 12 km of the grounding line, loss of these thicker and stiffer sections of ice tongue could perturb stresses at the grounding line enough to substantially increase inland flow speeds (∼900 m a⁻¹), grounded ice discharge, and Petermann Glacier's contribution to global sea level rise.

1 Introduction

Dynamic ice discharge from marine-terminating outlet glaciers is an important component of recent mass loss from the Greenland Ice Sheet (GrIS) (Enderlin et al., 2014; van den Broeke et al., 2016). Since the 1990s, outlet glaciers in Greenland have been thinning (Krabill et al., 2000; Pritchard et al., 2009), retreating (Carr et al., 2017; Jensen et al., 2016; Moon and Joughin, 2008), and accelerating (Joughin et al., 2010; Moon et al., 2012) in response to climate–ocean forcing. Marine-terminating glaciers are influenced by ocean warming (Holland et al., 2008; Mouginot et al., 2015; Straneo and Heimbach, 2013), increased surface air temperatures (Moon and Joughin, 2008), and reduced sea ice concentration in the fjords (Amundson et al., 2010; Reeh et al., 2001; Shroyer et al., 2017). However, glacier response to ocean–climate forcing is highly variable among regions and among individual glaciers due to differences in glacier topography and fjord geometry (Bunce et al., 2018; Carr et al., 2013; Porter et al., 2014). Moreover, changes at the terminus of these glaciers (i.e., calving or thinning) can reduce basal and lateral resistance which alters the force balance at the terminus and causes inland ice flow to accelerate. Indeed, 21st century retreat at two large outlet glaciers in southeast Greenland (Helheim and Kangerlussuaq) was followed by acceleration and ice surface thinning (Howat et al., 2005, 2007; Nick et al., 2009).

Floating ice shelves or tongues that extend out from outlet glacier grounding lines can also control a glacier's response to calving events (Schoof et al., 2017). Floating ice adjacent to the glacier grounding line can buttress inland ice, depending on the amount of shear and lateral resistance provided along the ice shelf margins (Haseloff and Sergienko, 2018; Pegler, 2016; Pegler et al., 2013). Consequently, thinning and retreat of ice shelves can reduce backstress, which can perturb the stresses at the grounding line and propagate increases in driving stress inland causing accelerated ice flow. Understanding how glaciers may respond to ice shelf loss is therefore important for estimating future flow speeds and, ultimately, their increased contributions to grounded ice discharge and global sea level rise. Considerable work has focused on the role of buttressing ice shelves on grounded ice dynamics in Antarctica (Goldberg et al., 2009; Gudmundsson et al., 2012; Reese et al., 2018b; Schoof, 2007), but less work has been done on floating ice tongues in Greenland, where large calving events have recently taken place (Box and Decker, 2011; Hill et al., 2017; Rignot et al., 2001).

Figure 1 Study location, Petermann Glacier, northwest Greenland. Observed ice speeds from the MEaSUREs program (winter 2009/2010, Joughin et al., 2010) across the Petermann Glacier catchment, which corresponds with our model domain. White lines show 300 m ice surface contours across the catchment. The thick black line is the glacier centerline and the thick red line is the glacier grounding line. Inset shows the location of newly prescribed terminus positions for each diagnostic perturbation experiment (A–H). Light turquoise shows floating ice, and white is grounded ice.

One of the last remaining ice tongues in Greenland is at Petermann Glacier, northwest Greenland. Petermann Glacier is fast flowing (∼1000 m a⁻¹, Fig. 1) and drains approximately 4 % of the GrIS by area (Hill et al., 2017; Rignot and Kanagaratnam, 2006). Mass loss is predominantly via high melt rates (10–50 m a⁻¹) beneath the ice tongue (Rignot and Steffen, 2008; Wilson et al., 2017) and also occurs via large episodic calving events (Johannessen et al., 2013). Formerly the glacier terminated in a 70 km floating ice tongue, but two well-documented large calving events in 2010 and 2012 removed ∼40 km of the tongue (Falkner et al., 2011; Johannessen et al., 2013; Münchow et al., 2014; Nick et al., 2012). Contrary to the behavior of glaciers terminating in floating ice elsewhere, large calving events at Petermann Glacier were noted to be followed by minimal glacier acceleration (<100 m a⁻¹) (Ahlstrøm et al., 2013; Münchow et al., 2014, 2016; Nick et al., 2012), suggesting that calving from the seaward parts of the ice tongue appear to have limited impact on flow upstream of the grounding line.

Several ice tongues have been lost from neighboring glaciers in northern Greenland since the early 2000s (C. H. Ostenfeld, Zacharaie Isstrøm, Hagen Bræ; Hill et al., 2017; Mouginot et al., 2015; Rignot et al., 2001), and with Arctic air and ocean temperatures predicted to increase in a warming climate (Gregory et al., 2004), the question remains, at what point will Petermann Glacier lose its ice tongue and how might its complete removal impact ice dynamics? Petermann's tongue has been retreating from the end of the fjord (∼90 km from the present grounding line) since the beginning of the Holocene (Jakobsson et al., 2018) and currently resides at its most retreated position in recent history (Falkner et al., 2011; Hill et al., 2018; Jakobsson et al., 2018). More recently (2016), another large rift formed across the ice tongue (Münchow et al., 2016), suggesting another large calving event is imminent. As Petermann Glacier is fast flowing and drains a large area of the GrIS (∼4 %), it has the potential to contribute to increased ice discharge and ultimately sea level rise, once it becomes grounded. Here, we attempt to answer the question at what point do large calving events from the Petermann ice tongue cause substantial acceleration (i.e., >100 m a⁻¹ that propagates inland of the grounding line) and increased ice discharge? To do this we use the community finite-element ice flow model Úa (Gudmundsson et al., 2012) to

• i.

  infer the stress conditions beneath the glacier catchment and along the ice tongue walls;

• ii.

  test whether the recent small changes in velocity following calving events in 2010 and 2012 can be replicated;

• iii.

  assess the future response (acceleration and ice discharge) of Petermann Glacier to further calving events, and eventual entire loss of the remainder of the ice tongue (Fig. 1).

First, we initialize the model using observational datasets of surface and basal topography. We then invert observed ice velocities prior to the 2010 calving event to determine the initial basal conditions (slipperiness and rheology of the ice). Finally, we perform a series of diagnostic perturbation experiments in which we remove sections of the Petermann Glacier ice tongue and assess the instantaneous glacier acceleration and increase in grounding line ice flux.

2 Methodology

2.1 Data input

To initialize our model, we used several observational datasets of glacier geometry and ice velocity. Ice surface topography of Petermann Glacier was taken from the Greenland Ice Sheet Mapping Project (GIMP) digital elevation model (DEM) (Howat et al., 2014). We also used this to calculate the surface drainage catchment, using flow-routing hydrological analysis in MATLAB TopoToolbox (Schwanghart and Kuhn, 2010). The defined catchment is ∼85 000 km², extends approximately 550 km inland of the grounding line, and encompasses the tributary glaciers flowing into the east side of the Petermann Glacier ice tongue (Fig. 1). Ice thickness and basal topography were taken from the Operation IceBridge BedMachine v3 dataset (Morlighem et al., 2017). These data were generated from radar ice thicknesses, ice motion, and the mass conservation method to resolve the basal topography and ice thickness of the GrIS (Morlighem et al., 2017). This version also includes high resolution bathymetry of Petermann Glacier fjord seaward of the ice tongue (Jakobsson et al., 2018; Morlighem et al., 2017). All three topographic datasets have a resolution of 150 m and a nominal date of 2007, which precedes the large calving event in 2010.

Firstly, to initialize our model via inversion (see Sect. 2.3), we required annual ice velocities prior to our first experiment, which is the calving event in 2010 (Table 1). These were taken from winter 2009/2010 from the Greenland MEaSUREs dataset (Table 1; Joughin et al., 2010). This dataset has a resolution of 500 m, and an average error of 8 m a⁻¹ across the entire Petermann catchment, which increases to 18 m a⁻¹ along the floating ice tongue (Table 1). To validate our modeled velocity changes in response to calving events in both 2010 and 2012, we required observed velocities from the years preceding and succeeding these events. Velocities from winter 2009/2010 (used for inversion) also acted as our baseline velocities that we compared with observed velocities after each calving event. However, Greenland-wide velocities for the winter following the 2010 calving event (2010/2011) were not readily available. Instead, a series of datasets exist that cover select regions of the ice sheet, derived from feature or intensity tracking of optical Landsat imagery (Howat, 2017; Rosenau et al., 2015) or synthetic aperture radar (SAR) imagery (Joughin et al., 2010). We do not use optical Landsat 7 ETM+ derived velocities because their coverage is restricted to midsummer, which may reflect seasonal speedups rather than the interannual impact of large calving events on ice velocity. Additionally, Landsat-derived velocities may have errors associated with cloud cover and/or the scanning line correction image banding from May 2003 onwards. Instead, we used a combination of SAR-derived velocities (TSX and PALSAR), which are not limited to the summer months and benefit from higher resolution, and more frequent repeat pass imagery (Table 1). High-resolution TSX imagery (100 m) was acquired from the MEaSUREs program (Joughin et al., 2010), but is limited to 11–45 km inland of the grounding line at Petermann Glacier. To supplement this we also used PALSAR-derived velocities (Table 1; Nagler et al., 2016), which provided additional coverage along the western half of the floating ice tongue. Average errors across the catchment are 4 and 16 m a⁻¹ for TSX and PALSAR, respectively (Table 1). Greenland-wide winter 2012/2013 (after the 2012 calving event) velocities were also acquired from the MEaSUREs program (Table 1; Joughin et al., 2010).

(Joughin et al., 2010)(Joughin et al., 2010)(Nagler et al., 2016)(Joughin et al., 2010)

Table 1 Velocity data sources for Petermann Glacier.

To determine observed velocity change, we differenced velocity fields after each calving event (2010/2011 and 2012/2013) from initial baseline velocities (2009/2010) (Fig. S1 in the Supplement). For 2009/2010 to 2010/2011 speedup along the ice tongue averaged 29 m a⁻¹, but further inland of the grounding line, noisy and unphysical velocity differences (Fig. S1) indicate that there was no coherent velocity change, and it is unlikely velocity changes propagated far inland. For clarity we present centerline velocity profiles in Fig. 2, which show that increases in speed were limited to the lower portions of the ice tongue. Observed velocity estimates presented here after the calving event in 2010 are within the range of previous studies, which showed a 30–125 m a⁻¹ speed increase along the ice tongue (Johannessen et al., 2013; Münchow et al., 2016; Nick et al., 2012) and limited change further inland (Nick et al., 2012). After the calving event in 2012, velocity increases averaged 79 m a⁻¹ along the ice tongue and propagated further towards the grounding line (Figs. 2 and S2).

2.2 Model initialization

To model the response of Petermann Glacier to ice tongue loss we used the finite-element model Úa (Gudmundsson et al., 2012). Úa solves equations of ice dynamics using the shallow ice-stream approximation (SSA) (Macayeal, 1989; Morland, 1987), a Weertman sliding law, and Glen's flow law. The momentum equation of the vertically integrated SSA can be written in the form

$$\begin{array}{}\text{(1)}& {\mathrm{\nabla }}_{xy}\cdot \left(h\mathbf{T}\right)-{\mathit{\tau }}_{\mathrm{bh}}={\mathit{\rho }}_{\mathrm{i}}gh{\mathrm{\nabla }}_{xy}s+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}{\mathrm{\nabla }}_{xy}{\mathit{\rho }}_{\mathrm{i}},\end{array}$$

where

$$\begin{array}{}\text{(2)}& {\mathrm{\nabla }}_{xy}=({\partial }_x,{\partial }_y{)}^{\mathrm{T}},\end{array}$$

and T is the resistive stress tensor defined as

$$\begin{array}{}\text{(3)}& \mathbf{T}=\left(\begin{array}{cc}\mathrm{2}{\mathit{\tau }}_{xx}+{\mathit{\tau }}_{yy}& {\mathit{\tau }}_{xy}\\ {\mathit{\tau }}_{xy}& {\mathit{\tau }}_{xx}+\mathrm{2}{\mathit{\tau }}_{yy}\end{array}\right).\end{array}$$

In the above equation s is the surface topography, h is the ice thickness, p_i is vertically averaged ice density, g is the gravitational acceleration, and τ_bh is the horizontal part of the bed-tangential basal traction τ_b.

Úa has been previously used to understand glacier behavior following ice shelf loss in Antarctica (De Rydt et al., 2015) and ice tongue collapse at the Northeast Greenland Ice Stream (Rathmann et al., 2017). It has also been used to assess the impact of buttressing ice shelves around Antarctica (Reese et al., 2018a). Previous modeling studies at Petermann Glacier have been conducted using a one-dimensional flow line approach (Nick et al., 2012). Here we aim to expand on earlier work at Petermann to assess if Úa, a two-horizontal-dimensional vertically integrated approach, can also replicate the observed velocity response to calving events in 2010 and 2012 and then be used to estimate the impact of future calving events on ice flow and discharge. This model is advantageous over a flow line approach, as it allows us to account for stresses in both horizontal dimensions, which can better assess the impact of ice shelf changes on the force balance at the glacier grounding line (Gudmundsson, 2013).

To set up the model, we used the surface velocity catchment (Fig. 1) as the outer computational boundary and imposed Dirichlet (essential) boundary conditions by fixing velocities to zero inland of the ice divide. Nunataks and rock outcrops along the east side of the ice tongue were digitized using Landsat 8 imagery and treated as holes within the mesh. Ice velocities along the nunataks' boundaries were set to zero, i.e., a no-slip boundary condition. Initially, it was less clear what type of boundary condition to impose along the margins of the ice shelf where it is in contact with the side walls. We began by imposing a free-slip boundary condition along the ice tongue margins but conducted further runs with no-slip boundary conditions along the side walls to determine which boundary conditions were most appropriate (see Sect. 2.4).

The initial calving front boundary was the location of the terminus in 2009, digitized from Landsat 7 ETM+ imagery, and in all cases a Neumann (natural) boundary condition was imposed along the terminus. Using this computational domain, and the finite element mesh generator Gmsh (Geuzaine and Remacle, 2009), we generated a high-resolution mesh, with 58 000 linear (three-node) elements, and ∼30 000 nodes (Fig. S3). The unstructured mesh capabilities of Úa allowed us to refine the mesh based on the observed velocity field. Where ice speeds are fastest (>500 m a⁻¹), primarily along the ice tongue, element sizes are 0.75 km, whereas element sizes inland have a maximum size of 2.7 km. Overall the mean element size is 1.52 km, with a median of 1.4 km. We also increased the mesh resolution of the slower-flowing (<500 m a⁻¹) tributary glaciers to the east of the Petermann Glacier ice tongue to 0.75 km. Topographic datasets (surface, bed, and ice thickness) and ice velocities observed before the calving (winter 2009/2010) were mapped onto this mesh using linear interpolation.

Figure 2 Observed and modeled velocities along the Petermann Glacier centerline before and after calving events in 2010 (a) and 2012 (b). Observed initial speed (MEaSUREs InSAR winter 2009/2010) is shown with a solid black line. Observed post-calving speeds after the 2010 calving event (PALSAR and TerraSAR-X winter 2010/2011) and the 2012 calving event (MEaSUREs InSAR winter 2012/2013) are shown with solid purple lines. Modeled initial speeds are shown in dashed black, and speeds after each respective modeled calving event are shown in dashed purple. Note the small change in both observed and modeled velocities, particularly inland of the grounding line.

2.3 Model inversion

Before modeling changes in the flow speed of Petermann Glacier due to perturbations in the calving front position, we must first estimate the prior stress conditions. We used Úa to invert the known velocity field (winter 2009/2010) before the calving event to simultaneously estimate parameters of basal slipperiness (C) and ice rate factor (A) across the catchment. To estimate these parameters, Úa uses a standard methodology whereby a cost function involving a misfit term and a regularization term is minimized. The gradients of the cost function with respect to A and C are determined in a computationally efficient way using the adjoint method. Here we used Tikhonov regularization involving both amplitude and spatial gradients of A and C. Values of regularization parameters were varied by orders of magnitude between 1 and 10 000 and then within range. We also experimented with different sliding law exponent values of m ($\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5},\mathrm{7},\mathrm{9}$) and found the results of our diagnostic experiments to be insensitive to the value of m (Fig. S3). We set the stress exponent in Glen's flow law to n=3.

To begin with, we inverted the model using a fixed zero velocity condition along the outer catchment boundary only, and we allowed the Petermann Glacier ice tongue to have free-slip boundary conditions (Sect. 2.2). In the following section we discuss two additional inversion experiments in which we varied the boundary conditions along the ice tongue. The model was inverted until the misfit converged, which was after 120 iterations. Resultant model velocities (U_mod) are in good agreement with observations (U_obs) as shown in Fig. 3. The mean percentage difference between observed and modeled velocities is 26 %, which equates to an absolute difference of 11 m a⁻¹ (Table 2). Absolute mean velocity difference increases to 28 m a⁻¹ in areas flowing faster than 300 m a⁻¹ and to 66 m a⁻¹ along the floating ice tongue, which is only 7 % of the average ice tongue speed (967 m a⁻¹). Flow velocities are not well resolved along the far northeastern tributary glacier, which is due to small ice thicknesses and poorly resolved bed topography from interpolation (errors of ∼150 m).

Figure 3 Observed (U_obs) winter velocities 2009/10 (a) in the proximity of the grounding line. (b) Modeled (U_mod) velocities. Panel (c) shows a normalized bivariate histogram of the velocity residuals, which are the difference between modeled and observed velocities in the vicinity of the grounding line (magenta). $\mathrm{\Delta }u=u_{\mathrm{mod}}-u_{\mathrm{obs}}$ and $\mathrm{\Delta }v=v_{\mathrm{mod}}-v_{\mathrm{obs}}$. u and v are x and y components of the velocity vectors, respectively.

By inverting the known velocity field (winter 2009/10), we can infer the basal conditions beneath Petermann Glacier. To our knowledge, a catchment-scale assessment of the basal slipperiness and ice stiffness has not been previously documented for this region. Some studies have examined the basal thermal state (Chu et al., 2018; MacGregor et al., 2016) or provided Greenland-wide slipperiness estimates (Lee et al., 2015). However, aside from these, little is known on the stress conditions of Petermann Glacier. Here, we provide a new record of the basal conditions beneath Petermann Glacier, which are important for understanding dynamic glacier behavior. For our initial inversion, the distribution of basal slipperiness (C), ice rate factor (A), and the misfit between observed and modeled velocities are shown in the first line of Fig. 4. Basal slipperiness was on average 2 orders of magnitude greater within 10 km of the grounding line ($C\approx \mathrm{7.2}\times {\mathrm{10}}^{-\mathrm{2}}$ m a⁻¹ kPa⁻³) than the rest of the grounded glacier catchment (mean $C\approx \mathrm{1.47}\times {\mathrm{10}}^{-\mathrm{4}}$ m a⁻¹ kPa⁻³, Fig. 4). Grounded ice across the Petermann catchment is on average stiffer ($A\approx \mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{8}}$ a⁻¹ kPa⁻³) than along the ice tongue ($A\approx \mathrm{7.4}\times {\mathrm{10}}^{-\mathrm{8}}$ a⁻¹ kPa⁻³). However, the misfit between observed and modeled velocities is highest along the ice tongue (Fig. 4g), which suggests that the ice rheology parameter (A) may not reflect the true stress conditions along the ice tongue. The distribution of basal slipperiness and the ice rheology parameter (A) for our initial inversion are discussed in more detail in the following section.

Table 2 Misfit between observed and modeled velocity for each boundary condition scenario.

Figure 4 Inversion experiments using three sets of boundary conditions along the floating ice tongue. Panels (a)–(c) show logarithmically calculated basal slipperiness (C) for boundary condition Scenarios 1 to 3, respectively, in which orange represents highly slippery areas. Glen's flow law rate factor A (d–f), for which light blue represents soft ice, and brown is stiff ice. The final column (g–i) is the absolute difference between observed (U_obs) and modeled velocities (U_mod) after inversion using each set of boundary conditions.

2.4 Boundary conditions

In an attempt to improve the misfit between observed and modeled velocities, and accurately replicate the lateral resistive stresses along the ice tongue margins, we conducted runs with both no-slip and free-slip boundary conditions along the side walls. We then tested which produced the best fit to observed velocities. Alongside our additional inversion (Scenario 1), in which no velocity condition was imposed along the ice tongue margins, we inverted the model using two further sets of boundary conditions. These are Scenario 2 with fixed velocities along the western margin of the ice tongue to zero (no-slip) and the east margin left as free-slip and Scenario 3 with velocities fixed to zero (no-slip) along both margins of the floating ice tongue. We then based our assessment of these boundary condition scenarios on three criteria: (i) the misfit between observed and modeled velocities, (ii) observations of the confinement and attachment of the ice tongue to the fjord walls in satellite imagery (i.e., heavy rifting or not), (iii) the ability of each set of boundary conditions to replicate the observed velocity response following the 2010 calving event. The first two criteria are discussed in this section, and the third criterion is discussed alongside our model experiments in the following section (Sect. 2.5). As before, we perform each inversion for 120 iterations until the misfit has converged, and we use the same values of m, n, and regularization parameters for each scenario. The slipperiness (C), ice rate factor (A), and misfit distributions ($\left|U_{\mathrm{obs}}\right|-\left|U_{\mathrm{mod}}\right|$) are shown for each boundary condition scenario in Fig. 4. Mean misfits between observed and modeled velocities for each scenario are in Table 2.

Slipperiness values showed a similar spatial distribution across all boundary condition scenarios, i.e., increasing towards the grounding line and decreasing further inland (Fig. 4a–c), and do not vary substantially within 10 km of the grounding line (range of (3.8–$\mathrm{7.2})\times {\mathrm{10}}^{-\mathrm{2}}$ m a⁻¹ kPa⁻³) or across the entire grounded catchment (range of (1.47–$\mathrm{2.18})\times {\mathrm{10}}^{-\mathrm{4}}$ m a⁻¹ kPa⁻³). In contrast, spatial variations in ice stiffness (A) were more obvious between each scenario, which corresponds to differences in the misfit distributions (Fig. 4). Average A values along the ice tongue varied by 3 orders of magnitude between Scenarios 1 ($A\approx \mathrm{7.4}\times {\mathrm{10}}^{-\mathrm{8}}$ a⁻¹ kPa⁻³) and 3 ($A\approx \mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{5}}$ a⁻¹ kPa⁻³). Scenarios 1 and 2 show stiff ice across the entire ice tongue, which does not reflect the stiff ice tongue center and weaker margins we would expect from lateral reductions in longitudinal strain rates and ice velocity associated with shearing along the fjord walls (Raymond, 1996). In both cases, velocities do not reproduce observations along the lateral margins and lower portion of the ice tongue well (criterion i, Fig. 4g–h).

In contrast, areas of softer ice ($A\approx \mathrm{4.6}\times {\mathrm{10}}^{-\mathrm{5}}$ a⁻¹ kPa⁻³) exist along the lateral margins of the lower portion of the ice tongue (Fig. 4f) when we used a no-slip boundary condition along both ice tongue margins during inversion (Scenario 3). In accordance with previous studies (Nick et al., 2012) and our own observations of satellite imagery, this replicates the apparent weak attachment of floating ice to the fjord walls in the lower and eastern parts of the tongue prior to the 2010 calving event (criterion ii). In Scenario 3, overall mean percentage difference between U_obs and U_mod was also improved by 6 % and the absolute difference was reduced by 15 m a⁻¹ along the ice tongue (Table 2). We find that imposing a no-slip boundary condition along both side walls of the ice tongue (Scenario 3) allowed the inversion procedure to automatically resolve the weak margins of the tongue. Based on criteria (i) and (ii), Scenario (3) therefore provides the most realistic distribution of ice softness along the ice tongue (criterion ii) and the best model fit to observed velocities (criterion i). These experiments have shown the importance of considering boundary conditions, particularly along floating ice margins in this study, for accurately determining lateral resistive stresses and replicating the observed velocity field.

2.5 Model experiments

Following model initialization and inversion, we performed a series of diagnostic experiments (Fig. 1) that perturb the calving front position to replicate previous large calving events and potential future loss from the Petermann Glacier ice tongue. We then examine the instantaneous velocity change with respect to our initial modeled velocities. As the focus of this paper is the impact of large calving events on glacier velocity, we do not incorporate ice loss via surface and/or basal melting. For each perturbation experiment, we removed all elements from the mesh downstream of the new calving front position and mapped all topographic datasets onto the new mesh. We then performed a forward-diagnostic model run, which solves the shallow ice stream equations (see Eq. 1) independent of time. In each case, the model is restarted from the previous experiment setup. During all experiments, grounding line position, boundary conditions, and ice thickness remained fixed. We also use our initially inverted parameters of basal slipperiness and ice rheology as inputs to all experiments. We acknowledge that in reality there is likely to be a period of relaxation and geometric adjustment after each individual calving event. Thus, by performing these sequential calving experiments in a diagnostic time-independent mode, we are modeling the immediate velocity change in response to a perturbation in buttressing at the terminus. We therefore do not allow any relaxation, adjustment, or deviation from our initial stress conditions (slipperiness or rheology) in between subsequent events. It is therefore possible that these estimates of ice flow and discharge may be higher than the transient glacier response over longer timescales.

We started by removing sections of the ice tongue that calved in 2010 and 2012 (Fig. 1), for which the new terminus positions were digitized from Landsat ETM+ imagery from 31 August and 21 July, respectively. We then assumed that the next iceberg to calve from the tongue will follow the path of the rift that formed in 2016. Then, we estimated the glacier response to future calving events, following two assumptions: (i) Petermann Glacier will continue to calve episodically, via rift propagation, back to the grounding line and (ii) future icebergs will be sized similarly to previous calving events (∼8 km long, Fig. 1). Each segment along the ice tongue acted as the new prescribed terminus position, which has a natural boundary condition. In reality, the size and nature of future calving events may vary (e.g., may be a series of small icebergs), but we conduct these experiments to assess the impact of future events similar in magnitude to previous calving. After each diagnostic experiment we calculated the vertically and horizontally integrated flux across the grounding line in gigatonnes per year with respect to pre-calving flux. We then convert to sea level equivalent (mm) by dividing by the volume of ice needed to raise global sea levels by 1 mm (361.8 Gt).

For the first diagnostic experiment, we used all three sets of boundary conditions (Scenarios 1–3) proposed in Sect. 2.4 to fulfill our third criterion (iii) of which scenario best replicates the small increase in velocity observed after the 2010 calving event (Fig. 5). Basal slipperiness (C), ice rate factor (A), and boundary conditions for each scenario were input into this first experiment (2010 calving) and the instantaneous increase in speed is presented in Fig. 5. We found that the differences in modeled velocity changes due to calving, using different sets of boundary conditions, were relatively small. Hence, our results are insensitive to the type of boundary condition applied (see also Sect. 2.4). This insensitivity to the type of side-wall boundary conditions can be understood to be related to our inverse methodology (see Sect. 2.3), in which A is inferred from measured velocities. In all cases the inversion was able to converge and provide a good fit to observations. Despite this, applying no resistance along the ice tongue margins (Scenario 1, free slip) produces no change in velocity along the tongue, which does not reflect observations (Figs. 2 and S1). However, when applying no-slip side-wall boundary conditions, modeled speed increases along the ice tongue (Scenarios 2 and 3) produced a better fit to observed changes, i.e., acceleration at the terminus that did not propagate far inland, and average 45 and 37 m a⁻¹, respectively. While both Scenarios 2 and 3 appear to adequately reproduce the observed velocity response after the 2010 calving event, we discount Scenario 2 due to the high misfit between modeled and observed velocities (Table 2) and unrealistic ice stiffness along the tongue (Fig. 4). Thus, boundary conditions and parameters of basal slipperiness (C) and ice rate factor (A) calculated in Scenario 3 are input into our subsequent diagnostic experiments after 2010.

Figure 5 Diagnostic perturbation experiments for the 2010 calving event, using three scenarios of boundary conditions applied along the floating ice tongue (a–c). Graduated white to red shows speed increase between initial modeled velocities (U_initial) and modeled velocities after the 2010 calving (U_calving). Bottom plots show observed velocities (U_obs) after calving (winter 2010/2011) and modeled velocities (U_mod) after the 2010 calving along the glacier centerline.

3 Results

3.1 Response to 2010 and 2012 calving

We have shown that the two-horizontal-dimensional model Úa can reproduce the flow of Petermann Glacier before the large calving event in 2010. Following this we removed sections of the ice tongue to replicate large calving events in 2010 and 2012, and we compare the model results with observed changes in flow speeds (Fig. 6a and b).

Figure 6 Diagnostic perturbation experiments at Petermann Glacier. The top panel shows a cross-sectional centerline profile of Petermann Glacier from the BedMachine v3 dataset (Morlighem et al., 2017). Letters A to H represent the points along the glacier centerline at which sections of the terminus were removed for each experiment. A is the 2010 calving event, B is the 2012 calving event, and C is the location of a large rift that formed in 2016. D to G are successive 8 km splices and H is the current grounding line location. Panels (a)–(h) show the modeled instantaneous increase in speed after each experiment with respect to initial pre-calving (before 2010) speeds.

The iceberg that calved away from Petermann Glacier in 2010 was ∼214 km² and on average 83 m thick (Fig. 6). Inverted ice rate factor A reveals that this section of the tongue is softer ($A\approx \mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{5}}$ a⁻¹ kPa⁻³) by 3 orders of magnitude than the rest of the ice tongue (Fig. 7b) or grounded glacier catchment ($A\approx \mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{8}}$ a⁻¹ kPa⁻³). Model results show that removing this section of the tongue was followed by a slight instantaneous increase in speed, ranging from ∼65 m a⁻¹ (6 % increase in tongue speed) at the terminus to <20 m a⁻¹ between 60 and 80 km along the centerline (Fig. 6a). These modeled velocity changes are in very good agreement with observed velocities presented here (Fig. 2) and documented in previous studies (Münchow et al., 2016; Nick et al., 2012). After the 2010 calving event, increases in speed across the entire ice tongue averaged 29 and 37 m a⁻¹ in observed and modeled velocities, respectively (Fig. 2). Modeled perturbations in flow speeds did not propagate far inland and averaged only +6 m a⁻¹ inland of the grounding line, which is below the average misfit between observed and modeled velocities (9.4 m a⁻¹, Table 2) and therefore indistinguishable from errors. Prior to the calving event, our modeled grounding line flux of 10.12 Gt a⁻¹ was within the range of previous estimates by Rignot and Steffen (2008) (12±1 Gt a⁻¹) and Wilson et al. (2017) (10.8±0.52 Gt a⁻¹). Limited changes in speed following the 2010 calving event were accompanied by little change in modeled grounding line flux (+0.14 Gt a⁻¹) and a negligible increase in sea level rise contribution (Fig. 7b).

In the next experiment, we removed a 96 km² section of the ice tongue to replicate a subsequent large calving event in July 2012. Importantly this calving event removed a thicker section of the ice tongue that averaged 111 m (Fig. 7a). Similar to the modeled dynamic response after 2010, ice flow speeds increased along the ice tongue after the 2012 calving event and did not propagate far inland of the grounding line (Fig. 6b). These modeled velocity changes are consistent with observed velocities in 2012/2013 (Figs. 2 and S2). Both modeled and observed speeds increased within the range of 3 %–5 % at the terminus and showed limited change inland of the grounding line (Fig. 2). The 2012 calving event was <50 % of the 2010 iceberg area and almost 4 times softer ($A\approx \mathrm{9.7}\times {\mathrm{10}}^{-\mathrm{5}}$ a⁻¹ kPa⁻³), suggesting it should provide less resistive stress. However, speed increases were 46 % greater along the ice tongue than in 2010 (averaging 59 m a⁻¹) and propagated further towards the grounding line. Despite some acceleration, the 2012 calving event had a limited impact on grounding line flux, increasing it by only 0.35 Gt a⁻¹ (3.4 %) compared to initial ice flux and increasing sea level rise contribution to 0.029 mm a⁻¹ (Fig. 7b).

Figure 7 Modeled experiment parameters for Petermann Glacier, for each diagnostic experiment (A–H), also shown in Fig. 6. (a) The black line shows iceberg area lost between each experiment (km²) and magenta is the average ice thickness (m) of each section of ice removed. (b) Green is the average ice rate factor (log¹⁰A) across the section of ice removed in each diagnostic experiment. Orange and purple lines represent the cumulative grounding line flux with respect to an initial ice flux of 10.12 Gt a⁻¹ and sea level equivalent contribution after the removal of each section of ice. (c) Average increases in ice speed across the entire ice tongue (red) and average ice speed within 10 km inland of the grounding line (blue) after each diagnostic experiment.

3.2 Response to future calving events

Calving events in 2010 and 2012 had a limited dynamic impact on the ice flow of Petermann Glacier and were followed by <10 % acceleration along the ice tongue and <2 % at the grounding line (Fig. 2). These modeled findings are consistent with observed velocity change (Fig. 2) and previous modeling of the 2010 calving event (Nick et al., 2012). After accurately replicating the observed velocity response to large calving events in 2010 and 2012, we were confident in the model's ability to estimate the instantaneous velocity response of Petermann Glacier to a change in stress conditions at the grounding line associated with removing large sections of the ice tongue.

We conducted six further experiments to analyze the glacier dynamic response (instantaneous change in flow speeds) and grounding line flux to large calving events. Each of the new calving front positions was approximately 8 km apart along the tongue, and ice loss area averages 125 km (Fig. 7a). First, we assume that the next calving event from Petermann Glacier will fracture along the path of a large rift that formed in 2016, removing a ∼154 km² section of the ice tongue (Experiment C, Fig. 6). This segment is also on average 35 m thicker and sturdier ($A\approx \mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{6}}$ a⁻¹ kPa⁻³) than the downstream section of the tongue that collapsed in 2010 and 2012 (Fig. 7b). In this case, average increases in speed along the ice tongue were greater (94 m a⁻¹) and propagated further (∼30 km from the terminus) towards the grounding line than after previous calving events (Fig. 6). Acceleration 10 km inland of the grounding line more than doubled to 24 m a⁻¹ compared to acceleration after the 2012 calving event (10 m a⁻¹). Despite this, acceleration inland of the grounding line remained 75 % smaller than increases along the ice tongue (Fig. 7c) and did not propagate far into the glacier catchment (Fig. 6c). After this experiment grounding line flux increased to 11 Gt a⁻¹ (+0.87 Gt a⁻¹) and sea level equivalent rose to 0.03 mm a⁻¹ (Fig. 7b).

Over the subsequent diagnostic calving experiments (D–H), there was a linear increase in average speed change across the ice tongue, as well as a near doubling of average speed increases immediately inland of the grounding line (within 10 km) between each experiment (Fig. 7). After removing the 164 m thick section D from the tongue, average ice tongue speeds increased by 21 % compared to initial velocities (∼954 m a⁻¹), but increases inland of the grounding line (within 10 km) remained small in comparison ($\sim +\mathrm{57}$ m a⁻¹). During the following three experiments (E–G), instantaneous average velocity increases across the tongue were more substantial than after previous calving events ranging from 304 m a⁻¹ after removing segment E to increasing by +1000 m a⁻¹ (>100 % of initial flow speeds) across the small remaining section of the ice tongue after experiment G. Throughout these experiments, higher-magnitude increases in speed propagated further into the catchment (∼10–15 km inland of the grounding line) than after previous calving events (Fig. 6). Simultaneous to increases along the ice tongue, average speed increases inland of the grounding line (10 km) went from 103 m a⁻¹ (experiment E) to 453 m a⁻¹ (experiment G, Fig. 7). Once the last remaining section of the ice tongue was removed (54 km², H), speed increases were double those of experiment G, reaching +900 m a⁻¹ immediately inland of the grounding line (10 km). Removing the entire ice tongue, and consequently detaching it from any tributary glaciers, also led to a ∼530 m a⁻¹ speed up at the terminus of Porsild Glacier (Fig. 6h).

Alongside linear increases in speed after large calving events, we also note positive trends in the thickness of each calved iceberg, grounding line discharge, and sea level equivalent (Fig. 7). Ice thickness along the Petermann Glacier tongue increases from ∼50 m towards the terminus to ∼500 m at the grounding line (Fig. 6; Münchow et al., 2014). In our experiments, the ice thickness of each segment increased by an average of 221 m (Fig. 7a). At the same time, grounding line discharge increased by an average of +2.17 Gt a⁻¹ after each experiment and once the entire ice tongue was removed, cumulative grounding line flux reached 25 Gt a⁻¹ and the glacier contribution to sea level rise increased to approximately 0.07 mm yr⁻¹ (Fig. 7b). As well as increases in ice thickness along the ice tongue, there is also a general increase in the stiffness of the ice back towards the grounding line. The ice is generally soft in the lower ∼40 km of the tongue (Fig. 4f) before ice rate factor (A) values decrease by 1–2 orders of magnitude during the last five experiments (E–H: Fig. 7b). Importantly, grounded ice immediately inland of the grounding line (within 10 km) is stiffer than the entire ice tongue.

4 Discussion and conclusions

Here, we expand on previous work and provide new insight into the velocity response of Petermann Glacier to past and future large calving events and eventual ice tongue collapse. In contrast to the removal of buttressing ice shelves elsewhere in Greenland (e.g., Jakobshavn Isbræ, Joughin et al., 2008; Zacharaie Isstrøm, Mouginot et al., 2015) and from the Antarctic Peninsula (e.g., Larsen B, De Rydt et al., 2015; Scambos et al., 2004), we show that Petermann Glacier was dynamically insensitive to the removal of ∼310 km² of the ice tongue via calving events in 2010 and 2012 (Fig. 6). After both calving events there was a limited increase in speed (<10 % of initial flow speeds, Fig. 2) that remained below the ∼22 %–25 % seasonal variability in flow speeds observed between 2006 and 2017 (Lemos et al., 2018; Nick et al., 2012). This insensitivity of ice velocities to large calving events can be explained by weak resistance provided by the lower portion of the ice tongue along its lateral margins (Fig. 4f; Nick et al., 2012). From this, we can conclude that the section of the ice tongue that calved away in 2010 and 2012 provided little frontal buttressing on grounded ice.

Given that several floating ice tongues have been lost from neighboring glaciers in northern Greenland (Hill et al., 2018), and the rapid nature of Petermann Glacier's Holocene retreat from the fjord mouth (Jakobsson et al., 2018), it is possible that the ice tongue will continue to calve episodically and in the not-too-distant future collapse entirely. We set out to determine at what point future calving events at Petermann Glacier (similar in magnitude to past calving) will cause substantial acceleration and increased ice discharge. The key conclusion of this work is that future calving events (C–E) from the lower portions of the ice tongue (>12 km from the grounding line) appear to be passive. We attribute the small modeled velocity response (<100 m a⁻¹ increase at the grounding line) to calving events from this lower portion of the tongue to be due to ice that is thinner (<200 m) and an order of magnitude softer, which provides limited buttressing on grounded ice. Indeed, if the next calving event takes the path of the 2016 rift formation, it is unlikely to substantially accelerate ice flow (Fig. 6c) or increase the glacier contribution to grounded ice discharge and sea level rise (Fig. 7). However, we find that removing sections of the ice tongue within 12 km of the grounding line (F–H) has a larger impact on ice flow speeds, increasing them by an average of 900 m a⁻¹ (96 %) after the entire ice tongue is removed (Fig. 6h). Alongside this, cumulative ice flux across the grounding line increases from 11 Gt a⁻¹ (Experiment C) to 25 Gt a⁻¹ (H, Fig. 7b) and cumulative sea level rise could reach 0.07 mm a⁻¹ (for event H). Importantly, within 12 km of the grounding line, the thickness and stiffness of the ice tongue increase dramatically (Fig. 7). As such, this part of the ice tongue provides greater lateral resistance along the fjord walls and is therefore more effective at buttressing grounded ice. Removing these sections of ice is thus likely to alter the resistive stresses at the grounding line enough to cause a greater increase in flow speeds that propagate further inland.

Overall, our findings show that Petermann Glacier has not responded dynamically to previous calving events in 2010 and 2012 and is unlikely to accelerate substantially after imminent future calving events (Fig. 6c). However, future large episodic calving events closer to the grounding line have the potential to perturb the stresses acting on grounded ice and substantially increase flow speeds and ice discharge (Fig. 7). Despite substantial increases in speed forecast after the ice tongue is removed, the question remains as to whether acceleration will be short-lived and the glacier will restabilize at the current grounding line position or retreat inland. Similar to when the glacier was buttressed by an ice shelf at the end of the fjord (Jakobsson et al., 2018), it may be that the current ice tongue has allowed grounding line stability, and its collapse will similarly lead to unstable grounding line retreat. Indeed, the eastern portion of the current grounding line lies within a deep bedrock canyon (Bamber et al., 2013; Morlighem et al., 2017), which may allow for marine ice sheet instability. However, we cannot discount the possibility that an ice tongue may regrow in the future (Nick et al., 2012). Here, we have estimated the instantaneous response of Petermann Glacier's ice flow to the immediate removal of sections of the ice tongue. Importantly, future work incorporating the transient evolution of the glacier geometry and grounding line position in between these large calving events is needed to assess the long-term response of Petermann Glacier to future ice tongue loss.