Modelled dynamic retreat of Kangerlussuaq Glacier, East Greenland, strongly influenced by the consecutive absence of an ice mélange in Kangerlussuaq Fjord

Abstract

Mass loss at the Greenland Ice Sheet is influenced by atmospheric processes controlling its surface mass balance, and by submarine melt and calving where glaciers terminate in fjords. There, an ice mélange – a composite matrix of calved ice bergs and sea ice – may provide a buttressing force on a glacier terminus and control terminus dynamics. Kangerlussuaq Glacier is a major outlet of the Greenland Ice Sheet, for which recent major retreat events in 2004/2005 and 2016–2018 coincided with the absence of an ice mélange in Kangerlussuaq Fjord. To better understand the response of Kangerlussuaq Glacier to climatic and oceanic drivers, a 2D flowline model is employed. Results indicate that an ice mélange buttressing force exerts a major control on calving frequency and rapid retreat. When an ice mélange forms in Kangerlussuaq Fjord, it provides stabilising forces and conditions favourable for winter terminus re-advance. When it fails to form during consecutive years, model results indicate that Kangerlussuaq Glacier is primed to retreat into the large overdeepenings in Kangerlussuaq Fjord, and to terminus positions more than 30 km farther inland, implying that excessive mass loss from Kangerlussuaq Glacier by the year 2065 cannot be excluded.

1. Introduction

The Greenland Ice Sheet (GrIS) has contributed 0.70 ± 0.15 mm a⁻¹ to global sea level rise between 2007 and 2017 (The IMBIE Team, 2020). This is not only a tenfold increase compared to the ice sheet's contribution between 1992 and 1997 (The IMBIE Team, 2020), and twice the contribution of the Antarctic Ice Sheet over the same period (Chen and others, 2017; The IMBIE Team, 2018), but also a lower bound for the rates to be observed in decades to come (Oppenheimer and others, 2019; Briner and others, 2020; Goelzer and others, 2020).

Mass loss from the GrIS is influenced by atmospheric processes controlling its surface mass balance (SMB) and interactions at the ice-ocean interface of the GrIS's outlet glaciers. There, submarine melt (SMM) and calving trigger changes to ice dynamics such as acceleration and dynamic thinning resulting in increased discharge (Enderlin and others, 2014; Fürst and others, 2015; Mouginot and others, 2019; Choi and others, 2021).

Spatio-temporal variations in the dominance of one mass loss mechanism over the other reveal that increasing discharge rates accounted for 66 ± 8% of the GrIS's net mass loss between 1972 and 2018 (Mouginot and others, 2019), with the near synchronous acceleration and retreat of the GrIS's four largest outlet glaciers – Jakobshavn Isbræ, Helheim Glacier, Petermann Glacier and Kangerlussuaq Glacier – accounting for 42% of the total GrIS discharge between 2000 and 2012 (Enderlin and others, 2014). In relation to rising air temperatures over the GrIS, SMB processes accounted for 84% of total mass loss since 2009 (Enderlin and others, 2014), and the behaviour of outlet glaciers has since become asynchronous: Helheim Glacier and Petermann Glacier have continued to retreat steadily since 2010 (Kehrl and others, 2017; Rückamp and others, 2019), Jakobshavn Isbræ began to readvance in 2016 (Khazendar and others, 2019), and Kangerlussuaq Glacier underwent a further episode of acceleration and retreat in the same year (Brough and others, 2019).

Observations at Kangerlussuaq Glacier (KG) show interannual variations in both the presence and strength of an ice mélange (Kehrl and others, 2017). This composite matrix of calved ice bergs and sea ice is suspected to impact calving frequency and terminus dynamics by exerting a buttressing force on the terminus of outlet glaciers (Amundson and others, 2010; Foga and others, 2014; Moon and others, 2015; Burton and others, 2018; Xie and others, 2019). To better understand the impacts of these variations is one of the motivations for simulating the future behaviour of KG using a 2D ice dynamical flowline model. In doing so, we subject KG to various ice mélange buttressing (IMB), SMB and SMM scenarios to determine the current drivers behind the recent rapid retreat events observed at the glacier and to assess which factors hold the greatest influence over glacier dynamics in a warming climate.

2. Regional setting

Kangerlussuaq Glacier (KG) is located on the GrIS's east coast, stretching c. 325 km from the GrIS's central parts, with an ice thickness up to c. 3 km, towards Kangerlussuaq Fjord (KF), where it terminates with a calving front c. 1 km tall (Morlighem and others, 2017) (Figs 1a,b). KF is c. 10 km wide, 200–900 m deep, and extends c. 75 km towards the Irminger Sea (Murray and others, 2010). It continues as an up to c. 300 m deep trough across the continental shelf, acting to funnel Atlantic Water and associated oceanic heat towards KG (Christoffersen and others, 2011; Inall and others, 2014; Bevan and others, 2019). Between 1972 and 2018, KG delivered 158 ± 51 Gt of ice from the GrIS central-east drainage basin to the ocean, surpassed only by Steenstrup-Dietrichson Glacier (northwest Greenland, 219 ± 11 Gt) and Jakobshavn Isbræ (central-west Greenland, 327 ± 40 Gt) (Mouginot and others, 2019). At KG's terminus, velocities regularly exceed 10,000 m a^-1 (Fig. 1c) (Joughin and others, 2008b; Murray and others, 2010; Brough and others, 2019). KG's retreat from its Little Ice Age (LIA) maximum extent, marked by a submarine moraine rising to 200 m below sea level c. 19 km downstream of the current terminus position (Kjeldsen and others, 2015; Batchelor and others, 2019) (Figs 1d, e), occurred shortly after its discovery in 1930 (Watkins, 1932; Koch, 1933), and was likely triggered by increasing ocean temperatures, but subsequently controlled by bathymetry (Vermassen and others, 2020).

Fig. 1. (a) Map of Greenland adapted from Moon and others (2021). Blue shaded region: drainage basin of Kangerlussuaq Glacier (KG), from Mouginot and Rignot (2019). Orange line: flowline of KG. (b) Ice surface and bedrock topography data from the Bedmachine v3 dataset (Morlighem and others, 2017) along the flowline, cf. (a), with horizontal mesh resolution highlighted at the ice divide and terminus. Sea level is represented by the blue line. (c) KG in KF. Colour scale represents average surface ice velocities in 2015 from the MEaSUREs project (Joughin, 2020). (d) The terminus of KG taken from Landsat 8 on 23/08/2020. Coloured lines: terminus positions of KG. Ice fronts from 2014 to 2020 were sourced for this study, whilst ice fronts predating this are from Khan and others (2014). (e) Zoom into (b), with terminus positions and colour scale from (d).

For several ensuing decades, KG featured a rather stable terminus position before the glacier gradually thinned, retreated and accelerated during the 1990s (Thomas and others, 2000; Khan and others, 2014). A rapid retreat of 5 km during 2004/2005 saw terminus velocities increase from c. 7500 to c. 13 000 m a^-1 (Luckman and others, 2006; Howat and others, 2007; Joughin and others, 2008a). This retreat was attributed to markedly increased ocean temperatures (Cowton and others, 2016; Millan and others, 2018), coinciding with extensive retreat and acceleration of many glaciers along Greenland's east coast (Seale and others, 2011). At KG, deep warm water masses in the fjord, combined with elevated surface water temperatures and strong katabatic winds, contributed to enhanced calving, and, notably, hindered the formation of a terminus-stabilising and winter-advance facilitating ice mélange in 2004 (Joughin and others, 2008b; Christoffersen and others, 2011; Kehrl and others, 2017; Bevan and others, 2019; Brough and others, 2019).

The annual formation of a winter ice mélange returned in 2005 and consistently prevented winter calving as KG followed a seasonal pattern of winter advance and summer retreat, entailing fluctuations in its terminus position of 2–5 km, while steadily readvancing into KF until c. 2015 (Kehrl and others, 2017; Bjørk and others, 2018). During 2016–2018, KG retreated a further 5 km to a position atop a reverse sloped section of bedrock, priming KG for yet another retreat (Weertman, 1974; Joughin and Alley, 2011). The retreat coincided with the lack of a protective ice mélange which failed to form over two successive winters, resulting in terminus velocities increasing from c. 7000 to 10 000 m a^-1 (Brough and others, 2019) and a rapid acceleration in ice discharge (Hansen and others, 2021). When an ice mélange materialised again, in the winter 2018/19, KG began to readvance into KF (Bevan and others, 2019; Brough and others, 2019).

3. Methods

The finite-element model Elmer/Ice (Gagliardini and others, 2013) is used to model KG. Elmer/Ice is an open source software (https://elmerice.elmerfem.org/) designed to simulate ice-sheet dynamics, glacier flow and calving at marine terminating margins by solving the Full Stokes equations (e.g. Greve and Blatter, 2009), here combined with the crevasse depth criterion (Nye, 1957; Benn and others, 2007; Nick and others, 2020; Otero and others, 2010; James and others, 2014; Todd and Christoffersen, 2014) to model calving. Full Stokes ice dynamics including calving can be modelled in 3D (Todd and others, 2018), however, a 2D flowline approach is chosen here because the considerable reduction in computational time and cost compared to 3D simulations allows for more flexibility and a larger number of runs in exploring KG's response to a forcing that so far has received little attention, viz. IMB.

3.1 Model setup

A 2D flowline model is set up for KG, stretching c. 350 km through KG's entire drainage basin (sourced from Mouginot and Rignot (2019)) and embedded in the region where the highest ice velocities were observed from the MEaSUREs project (Joughin, 2020) (Fig. 1c). Ice thickness and bed topography along the flowline were extracted from the Bedmachine v3 dataset (Morlighem and others, 2017) and used to render a closed 2D domain (Fig. 1b) within which a mesh was produced using GMSH (Geuzaine and Remacle, 2009). The mesh for KG is composed of 21 008 triangular elements; ranging from c. 2 km (length along the flowline) in the upstream regions to c. 20 m when approaching the terminus where high spatial resolution is needed to accurately capture calving dynamics. Because KG is more than 200 km wide in the upper parts of its drainage basin and narrows to less than 10 km at its front, our setup accounts for thickening of the glacier where the flowpath is narrowing – this is one manifestation of the lateral effects otherwise neglected in 2D flowline models. This 3D effect of convergent flow has previously been accounted for by modifying the surface mass balance (Cook and others, 2014) but is here approximated by adding a glacier-width and downstream velocity-related source term that also scales with element size to the mass conservation equation, as suggested by Todd and Christoffersen (2014) and Passalacqua and others (2016). This modification renders the 2D flowline model to what could be referred to as 2.5D model.

Boundary conditions applied include a stress-free ice surface and subaerial part of the calving front, and Weertman-type sliding where KG is grounded using a spatially uniform friction coefficient of 0.025 and the exponent of basal sliding velocity taken as 1/3 (Weertman, 1974; Gagliardini and others, 2013). Whether KG is grounded or floating is determined by using a height-above-buoyancy criteria (Nick and others, 2020; Favier and others, 2012) viz. the comparison of water pressure and ice overburden pressure. At floating points (such as the underside of a glacier tongue lifted by buoyant forces) and at the submerged calving front, hydrostatic pressure is applied as boundary condition. At the calving front, hydrostatic pressure can be replaced by a buttressing force to mimic an ice mélange at KG's terminus if it is determined to be the greater force. At the upstream margin of the modelling domain, ice velocities are set to zero. KG's thermal state is given from a steady-state simulation of the glacier. A modelled temperature profile of Jakobshavn Isbræ  by Funk and others (1994) provides the initial conditions for the steady-state simulation, forced with a geothermal heat flux of 40 mW m^-2 at the ice base, and ice surface temperatures that are derived from the 2000–2019 average surface ice temperatures along KG's flowline (Hall, 2019). This thermal state is then fixed and maintained throughout our prognostic runs.

3.2 Mesh evolution and calving

During prognostic runs, spatio-temporal changes in modelling domain – and hence, mesh – are a reflection of ice flow and the boundary conditions imposed. Generally mesh evolution is handled by standard algorithms in Elmer/Ice discussed in Gagliardini and others (2013), yet any attempts to model marine-terminating glaciers must solve the mechanisms behind calving, a process in which a glacier expels large quantities of ice at its terminus aided by the propagation of crevasses in response to stresses. Calving together with submarine melt is referred to as frontal ablation, and relates to ice velocity and changes in frontal position as follows (Luckman and others, 2015):

(1)$${{\rm d} l\over {\rm d} t} = v_{{\rm ice}} - \dot{a} ,\; $$

where $\dot {a}$ is the frontal ablation rate, v _ice the terminus ice speed, and dl/dt the change in terminus position over time. In our setup, melt rates will be parameterised, while calving rates are determined based on investigations of crevasse depths: The depth of surface crevasses, C _s, is described following Nye (1957):

(2)$$\textit{C}_{{\rm s}} = {{2{\sigma}_{xx}\over \rho_{\rm w} g}} ,\; $$

where σ _xx is the longitudinal Cauchy stress and ρ _w is the density of water with g representing gravitational acceleration. Similarly, the height of basal crevasses, C _b, is defined by Nick and others (2020):

(3)$$\textit{C}_{{\rm b}} = {2{\sigma}_{xx} + p_{\rm w}\over \rho_{\rm w} g} ,\; $$

where p _w is water pressure. The crevasse depth criterion, first proposed by Benn and others (2007) and numerically investigated by Nick and others (2020), aids numerical solutions of the calving problem by stating that a calving event occurs when surface and basal crevasses merge, viz. penetrate through the entire glacier thickness and thus define a new calving front. The crevasse depth criterion was implemented in Elmer/Ice by Todd and others (2018) in 3D, and by Todd and Christoffersen (2014) in a 2D flowline model, to replicate seasonal calving of Store Glacier. In the 2D setting, they found the crevasse depth criterion to underestimate crevasse penetration, likely because a 2D setting does not account for the full stress regime in the ice. This can be compensated for by multiplying longitudinal stresses in the crevasse depth criterion with a tuning factor. For KG, a calving tuning factor of 1.13 was applied in all simulations. To solve Eqn (1), changes of the terminus position in response to submarine melt rates are hence evaluated, after which the crevasses depth criteria is the assessed. If fulfilled then re-meshing takes place, defining a new calving front position at the furthest inland point where surface and basal crevasses meet. Such re-meshing maintains the spatial resolution of the domain.

3.3 Model forcings

Two different SMB profiles are used, SMB_ref and SMB_RCP8.5. Both are based on 5 km spatial resolution SMB maps of Greenland, derived from the regional climate model MAR (Fettweis and others, 2017). For SMB_ref, MAR is forced with ERA-Interim reanalysis data (Dee and others, 2011) to create a single SMB map that represents the averaged annual SMB over the years 1979–2014, and from which SMB values are extracted along KG's flowline (Fig. 2a). SMB_RCP8.5 is comprised of 50 individual SMB profiles each along KG's flowline, representing the years 2015–2065, an extracted from SMB maps obtained from MAR, forced with output data from MIROC5 (Watanabe and others, 2010) for RCP pathway 8.5 (Fig. 2a).

Fig. 2. (a) Average SMB along the KG flowline between the years 1979–2014 from MAR (Fettweis and others, 2017) used in the SMB_ref (blue) alongside annual RCP 8.5 profiles from 2015 (light orange) to 2065 (dark orange). (b) Annual mélange backstress pressure (blue) and submarine melt rate (red) for profiles IMB_ref and SMM_ref respectively. (c) Average velocities along the flow line during the final year of the spin up (green) compared with average annual velocities along the flow line in 2015 from the MEaSUREs dataset (dotted purple). (d) Original 2D domain (stippled red) created from the Bedmachine v3 dataset (Morlighem and others, 2017), and modelled minimum summer (orange) and maximum winter (blue) extents as in the final year of the spin-up forced with SMB_ref, SMM_ref, and IMB_ref.

Further, four different SMM forcings are used, SMM_ref, SMM_×2, SMM_×3 and SMM_×4. While underwater instruments measuring water temperature can be placed near calving tidewater glacier fronts to infer melt rates (Holmes and others, 2019), fjord waters at the terminus of KG being continually littered with calved ice bergs make such a task difficult. Instead, instruments have been placed away from the terminus in KF and are used in conjunction with fjord circulation and plume models to provide SMM estimates (Cowton and others, 2015; Carroll and others, 2016; Fraser and others, 2018). As SMM are dependent on water temperatures and subglacial discharge (Motyka and others, 2003; Xu and others, 2013), we vary SMM rates in SMM_ref between a winter and summer state. Winter (December to June) melt rates, in association with cool waters and minimal subglacial discharge, are set to 0.2 m d^-1 as found by Fraser and others (2018) modelling oceanic circulation in KF. Summer melt rates, associated with greater heat transport in KF and the presence of melt water plumes from increased subglacial discharge, are set to peak 4.0 m d^-1 following a melt water plume modelling study by Carroll and others (2016). As melt water plumes are not spatial consistent across a glacier's calving front, and thus difficult to represent in 2D, we apply a sinusoidal transition from winter to summer melt that allows extreme SMM from a melt plume to peak in late August (Fig. 2b). To account for submarine melt rate variations with water depth (Carroll and others, 2016), SMM_ref is interpolated linearly from values of 0 m d^-1 at the waterline on the calving front to the defined melt rates (winter, peak summer, interpolated) at the depth of the grounding line. For SMM_×2, SMM_×3 and SMM_×4, the winter melt rate remains unchanged, however, the peak summer melt rate is doubled, tripled and quadrupled from the originally defined summer melt rate in SMM_ref, to 8, 12 and 16 m  d^-1 respectively. This allows to determine the impact of a likely doubling in SMM by 2065 at KG (Slater and others, 2020) and the response of the glacier to further extreme melt rates. For simulations spanning 50 years, submarine melt rates in SMM_×2, SMM_×3 and SMM_×4 are increased linearly (from 4 m d^-1) during simulation years 1-40, and then maintained at their maximum values during simulation years 41–50.

Also, four different IMB forcings are used, IMB_ref, IMB_×0.5, IMB_×0.25 and IMB_×0. While the behaviour and presence of KG's ice mélange is well documented (Joughin and others, 2008a; Christoffersen and others, 2011; Kehrl and others, 2017), its buttressing force on KG's front is yet to be quantified. The only in-situ quantification of an ice mélange's buttressing force is from Store Glacier, West Greenland (Walter and others, 2012), where IMB is reported to range from 30 to 60 kPa. With an IMB of 45 kPa, Todd and Christoffersen (2014) were able to reproduce the seasonal dynamics observed at Store Glacier in a 2D flow line model. Similarly, Vieli and Nick (2011) found the value of 40 kPa suitable for modelling ice mélange at Jakobshavn Isbræ. Hence, a backstress force of 45 kPa, σ _fb, is chosen for IMB_ref. This is then used to produce a per metre of mélange backstress, σ _M, of 198 kPa following Todd and Christoffersen (2014):

(4)$$\sigma_{\rm M} = \sigma_{{\rm fb}}{H_{{\rm T.Avg}}\over H_{\rm M}} ,\; $$

where H _T.Avg is the average height of the terminus, 597 m (Morlighem and others, 2017), and H _M is average mélange thickness. While no information regarding KF's ice mélange thickness is available, a LiDAR survey found the mélange thickness at the nearby Helheim Glacier to be 135 m (Cook and others, 2014). Here, we thus employ H _M = 135 m, and assume that σ _M acts on KG's submerged terminus (0–135 m below sea level), following a seasonal cycle observed by Kehrl and others (2017) during which the mélange becomes increasingly rigid during December, maintains its strength and presence from January till June, before disintegrating (Fig. 2b).

For IMB_×0.5, IMB_×0.25 and IMB_×0, σ _M is decreased from 198 to 100, 50 and 0 kPa, respectively, in simulations running 50 years into the future. This reduction is chosen to explore KG's sensitivity to a reduced ice mélange buttressing force, where, however, the relative importance of possible future changes in σ _bf, H _T.Avg and H _M cannot be resolved by simply reducing σ _M, see (4). The decrease in IMB forcing is applied linearly over years 1–40, and then maintained at its final value during years 41–50. Mimicking observations at KG of an ice mélange failing to materialise, a further three IMB forcings scenarios are used within which σ _fb = 0 kPa (‘skipped’), and where hydrostatic pressure is instead applied during years 11, 21 and 22, 31 to 33, and 41 to 44.

3.4 Spin-up and prognostic runs

To ensure that changes in ice dynamics occur in response to model forcings and not initial adjustments to the prescribed basal topography, ice thickness, ice temperature and boundary conditions, a spin-up simulation is conducted. During the spin-up, KG is forced by SMB_ref, SMM_ref and IMB_ref, for a duration of 50 years at a time step of 2 days, with the crevasse depth criterion activated. During the final 10 years of the spin-up, KG settles on a seasonal advance and retreat cycle of c. 3 km, between a winter maximum in May and a summer minimum in November (Fig. 2d), replicating the position and scale of fluctuations observed at the terminus in 2015 (Kehrl and others, 2017; Bevan and others, 2019). Furthermore, spun-up annually averaged ice surface velocities of 8140 m a^-1 at KGs terminus also match observations of 8516 m a^-1 in 2015, retrieved from the MEaSUREs dataset (Joughin, 2020), well (Fig. 2c). This spun-up state of KG serves as initial condition for all prognostic simulations that aim to asses the possible behaviour of KG in response to a warming climate (Table 1). All prognostic runs start from the glaciers state at the end of the final spin-up year and are again run for 50 years at a 2-day timestep with the final state thus reflecting a modelled state of KG in 2065.

Table 1. A list of prognostic simulations and their associated SMB, SMM and IMB forcings

4. Results

Results from 13 prognostic runs (cf. Table 1) departing from KG's spun-up state are shown in Figure 3. It is seen that altering the annual SMB in run SMB8.5 has very little effect on KG, with the calving front retreating only c. 0.5 km during the final 20 years (first row in Fig. 3).

Fig. 3. Evolution of KGs geometry and terminus position during prognostic simulations (Table 1). Glacier configurations are shown by colour coded lines (blue: 0 model years; red: 50 model years) and are plotted along bed topography (Morlighem and others, 2017)) cf. Figures 1 and 2. The vertical and horizontal scales on each panel are consistent with the exception of IMBx0s, where the horizontal scale is extended, and COMB1, COMB2 and COMB3, where both horizontal and vertical scales are extended to accommodate for the larger retreat of KG.

Similarly, altering SMM profiles has little impact on KG's behaviour for low-melt runs SMMx2 and SMMx3, although some thinning is observed in SMMx3 compared to SMMx2. For the high-melt run SMMx4 (melt rate 16 m d^-1), KG retreats by c. 2 km, combined with a thinning of more than 100 m. Retreat is initiated at around year 20, when the melt rate reaches 10 m  d^-1 (see Model forcings; second row in Fig. 3).

Further, KG's terminus position remains essentially unchanged in scenario IMBx0.5 where reducing the IMB by 50% to 100 kPa per metre of mélange continues to prevent winter calving and preserves the seasonal cycle of advance and retreat (third row in Fig. 3; Fig. 4). During run IMBx0.25, KG retreats c. 1.5 km and thins by c. 150 m. These changes occur after c. 30 model years, where the IMB reduces below 90 kPa per metre of mélange and winter calving begins to occur, thus weakening the seasonal winter advance. By reducing the IMB by 100% (run IMBx0), KG thins by more than 200 m, retreats by c. 5 km, and settles on a small pinning point at c. 324 km along the flowline. Already before the complete ice mélange-free state is attained in model year 41 (see Model forcings), KG loses its seasonal winter advance and summer retreat cycle (by model year 34) and retreats to a terminus position at c. 327 km along the flowline. There, it remains stable for some years, before the final retreat onto the bedrock pinning point takes place.

Fig. 4. Change in KG terminus position during prognostic runs IMBx0.5 (orange), IMBx0.5s (red) and COMB1 (blue). Green shaded time period: years when IMB forcing is skipped, affecting runs IMBx0.5s and COMB1 but not IMBx0.5 (cf. Table 1).

In IMBx0.5s (fourth row in Fig. 3) a decreasing IMB is coupled with scenarios of successive skipped mélange years mimicking the years when an ice mélange fails to form. During the first year of a skipped mélange scenario, KG retreats c. 5 km from its previous winter maximum to c. 327 km along the flowline while subsequent ice mélange-free years do not induce further retreat (Fig. 4). Instead, KG maintains its position at c. 327 km along the flowline and begins to readvance once the IMB is reapplied. The same behaviour is observed in IMBx0.25s, albeit KG thins by a further c. 20 m. For run IMBx0s, KG behaviour is similar to that in IMBx0.5s and IMBx0.25s until model year 34, when the IMB is no longer strong enough to produce a winter advance and KG stabilises at 327 km along the flowline, cf. IMBx0.5s. In year 41 (thus, 2 years earlier than in IMBx0), KG retreats suddenly to the bedrock pinning point at 324 km along the flowline, from where it retreats even further in year 48 into the overdeepening basin before the simulation ends with the terminus at c. 309 km along the flowline, representing a c. 25 km retreat.

During COMB1, COMB2 and COMB3 (fifth row in Fig. 3), where alterations of SMB, SMM and IMB-forcings are combined (cf. Table 1), KG undergoes extensive retreat. In COMB1, KG follows the same pattern of behaviour as IMBx0.5s until year 33 when the glacier is unable to maintain its position at c. 327 km along the flowline during a skipped mélange year and retreats to a bedrock pinning point at 324 km (Fig. 4). By year 41, further retreat is initiated and continued along the reverse-sloping bedrock before the terminus reaches a very pronounced pinning point elevating 230 m above the surrounding seafloor, at 298 km along the flowline, and settles there until the end of the simulation. KG's behaviour in COMB2 and COMB3 initially follows the same pattern as COMB1, although KG begins its retreat from the pinning point at 324 km earlier, at years 37 and 33 of the simulation respectively. Unlike COMB1, KG is unable to maintain its terminus position at the pronounced pinning point at 298 km along the flowline, and instead retreats a further c. 10 km until stabilising on a steep forward slope at c. 288 km along the flowline.

Figure 5 shows the temporal evolution of average monthly calving extent, and associated changes in the terminus position, for runs SMMx4, IMBx0, IMBx0s and COMB3. Initially, all runs follow a similar annual calving pattern; the IMB prevents winter calving, large post-mélange calving events in excess of 200 m in extent occur in June and regular calving events averaging c. 150 m in extent occur for the remainder of the year. Such a pattern is observed throughout SMMx4 although regular calving events increase in size to a minimum average of c. 175 m during the final 10 years. During IMBx0, large calving events averaging c. 275 m begin to occur during winter when the IMB reduces below 90 kPa in year 21, with winter calving becoming increasingly regular as the ice mélange strength decreases and fails to provide a winter advance. A similar pattern is found in IMBx0s, yet calving rates increase in magnitude to average c. 375 m in extent during year 48 as KG retreats from the pinning point at 324 km along the flowline into a reverse sloped bedrock. The same increase in calving magnitude is found in COMB3 as KG retreats into a reverse bedrock during year 33, with calving extent decreasing in magnitude when KG reaches a steep forward slope at c. 288 km along the flowline.

Fig. 5. The average calving extent for each month of the prognostic simulations (vertical stipples, left axis) and changes in terminus position (solid line, right axis) during prognostic runs SMMx4 (green), IMBx0 (orange), IMBx0s (red) and COMB3 (blue). Note the change in right x-axis scale between the two graph panels.

Figure 6 highlights similarities in the timing and specifics of KG's retreat behaviour, exemplified by the retreat from the pinning point at c. 324 km along the flowline during IMBx0s, COMB1, COMB2 and COMB3. Retreat coincides with the absence of a winter advance, either because the IMB is no longer strong enough to produce an advance (COMB2 and IMBx0s) or due to a skipped mélange year (COMB1, COMB3). Furthermore, the timing of retreat occurs in July, when melt rates are ramping from subdued winter levels towards the summer maximum (cf. Fig. 2b). The pattern and scale of retreat from the pinning point is again similar. All runs are found to retreat c. 20 km from the pinning point in the first year, before briefly stabilising then undertaking a further c. 7 km retreat to a large pinning point at 298 km along the flowline.

Fig. 6. Change in KG terminus positions during runs IMBx0s (red), COMB1 (light blue), COMB2 (medium blue) and COMB3 (dark blue), 4 years before and 4 years after KG retreats from a stabilising pinning point into a reverse sloped section of bedrock topography that deepens c. 250 m (cf. Fig. 1e).

Figure 7 summarises the final results of all prognostic runs (cf. Table 1, Fig. 3), showing the total retreat of KGs front position, both quantitatively (in km only) and in relation to bedrock configuration, viz. position along the flowline, after 50 model years, representing year 2065.

Fig. 7. (a) A comparison of KG's terminus retreat for the prognostic runs (mimicking 2015–2065) described in Table 1. (b) KG's final (2065) configuration for runs COMB3 (blue), IMBx0s (red), IMBx0 (orange), SMMx4 (green) and SMB8.5 (purple) along side the initial spin-up state (dashed-black).

5. Discussion

A first indication of the importance of oceanic drivers for the dynamic behaviour of KG is obvious from Figure 3: KG is virtually unaffected when its SMB is derived from RCP 8.5 (run SMB8.5, top row), instead of the continued use of an average SMB from the years 1979–2014 (results not shown). This is in agreement with Nick and others (2013) and Beckmann and others (2019) who, through 2D flowline investigations, found a SMB component forced by a RCP 8.5 scenario to have little impact of KG's retreat over the coming century. A possible explanation for the quasi steady state of KG under this type of forcing is that any possible changes to SMB in southeastern Greenland are likely restricted to low elevations affecting the fastest flowing ice (Fig. 2a), while changes in SMB are forecast reach far inland, and to high elevations, in West Greenland (Fettweis and others, 2013; Noél and others, 2021). However, through a coupling with hydrological processes not considered in our model set-up, SMB might influence glacier dynamics. For instance, the drainage of supraglacial meltwater lakes, and associated routing of meltwater from the glacier surface to the glacier base, has been suggested to have possible implications for seasonal and inter-annual changes in ice flow velocities (Joughin and others, 1996; Zwally and others, 2002; Chudley and others, 2019; Yang and others, 2019; Lampkin and others, 2020). Also, basal glacier mass balance is suspected to interact with SMB, and to augment overall mass loss (Karlsson and others, 2021). Accounting for SMB and basal mass balance induced hydrological changes at KG and assessing their interplay with other drivers of calving, such as enhanced plume formation leading to glacially modified water masses and hence modified SMM (Beaird and others, 2018), remains to be explored.

Observations have supported the notion that rising SMM promotes calving through the undercutting of glacier termini, leading to increased mass loss (Motyka and others, 2003; O'Leary and Christoffersen, 2013; Ma and Bassis, 2019). By 2065, SMM rates at KG are expected to have doubled to 8 m d^-1 (Slater and others, 2020), which, according to Figure 3 (second row), does not have a significant impact on KG's terminus position. Further, SMM rates of 12 m d^-1 (expected at KG by 2100 (Slater and others, 2020)) also fail to initiate a retreat and it is only when SMM reaches an extreme rate of 16 m d^-1 in scenario SMMx4, a 400% increase from current estimated peak summer melt rates (Carroll and others, 2016), that a noticeable retreat and an increase in calving extent is observed (Figs 3 and 5). Such insensitivity in KG's dynamic behaviour to SMM, when varied in isolation from other climate forcings, finds agreement with other 2D flowline modelling studies operating on Store Glacier, Helheim Glacier and a synthetic glacier (Todd and Christoffersen, 2014; Cook and others, 2014; Krug and others, 2015), all of which rest on a flat or forward facing slope and remained impervious to rises in SMM. To what extent such factors, or the 2D flowline nature of the model environment, contribute to the low sensitivity of raising only the SMM forcing and oppose observations remains to be explored. KG's modelled behaviour is also in contrast to a proposed linear relationship between SMM and terminus position (Slater and others, 2019) used to parameterise ice-ocean dynamics in the ISMIP6 project (Nowicki and others, 2016), hinting that factors other than SMM, or SMM in combination with other forcings, as indicated by the results from the COMB1-COMB3 scenarios, are key in controlling KG's future retreat and mass loss.

A reduced winter calving flux is recognised to limit discharge from the GrIS and has been linked to the presence of sea ice or an ice mélange at glacier termini, providing sufficient backstress to suppress calving and promoting conditions for frontal advance (Joughin and others, 2008a, b; Amundson and others, 2010; Moon and others, 2012; Kehrl and others, 2017; Amundson and Burton, 2018). This is supported by the modelled behaviour of KG (Fig. 3, third row), where an ice mélange backstress decreasing to 100 kPa per metre of mélange in run IMBx0.5 continually impedes winter calving, thereby allowing KG's terminus to re-advance, countering any summer retreat and maintaining the glacier's position in KF (Fig. 4, IMB1). Reducing the IMB force to below 90 kPa per metre of mélange in runs IMBx0.25 and IMBx0.0 allows for winter calving. Such a threshold aligns with work on a 2D synthetic glacier by Krug and others (2015), where a buttressing force between 130 and 90 kPa was found sufficient to impede calving. Similarly, Amundson and others (2010) estimated that a per meter of mélange force of 100 and 160 kPa, for a grounded terminus and a floating glacier tongue, respectively, would impede calving. As such, the reference IMB (based on Walter and others (2012)'s study of Store Glacier) and the derived IMB scenarios (cf. Table 1) provide useful pathways to explore the impact of IMB on glacier dynamics.

The absence of IMB initiates a period of extended calving and retreat (Fig. 3, fourth row), as observed at KG in the winter of 2004/05 (Joughin and others, 2008a; Christoffersen and others, 2011) and winters 2016/17 and 2017/18 (Bevan and others, 2019; Brough and others, 2019). During the initial year when IMB is skipped (Fig. 4 IMBx0.5s), KG retreats 5 km from its previous winter maximum, mimicking what occurred in the winter 2004/05. Subjecting KG to multiple years without a winter IMB does not induce further retreat (Fig. 4 IMBx0.5s), aligning with observations of KG from 2016 to 2018 when the ice mélange failed to form over two successive winters, yet KG did not retreat during the second year (Bevan and others, 2019). This suggests that in KG's current state, the number of years without a winter ice mélange is not proportional to the extent of retreat.

Moreover, it is suggested that IMB plays a vital role in controlling glacier dynamics at KG. A consistently re-occurring winter ice mélange allows KG to maintain its position in KF, with only extreme SMM rates found to induce a small retreat (Fig. 3, SMMx4), and while a failed mélange induces rapid retreat, its re-appearance allows KG's re-advance into KF towards its previous maximum (Fig. 4 IMBx0.5s). Such behaviour has also been observed at Helheim Glacier (Kehrl and others, 2017) and Jakobshavn Isbræ (Amundson and others, 2010; Joughin and others, 2020), and raises the question whether open water conditions and glacial retreat, and the presence of an ice mélange and glacial advance, are correlated (Howat and others, 2010; Carr and others, 2013; Moon and others, 2015). Sea surface temperatures (SST) are suggested to have a large influence on the rigidity and buttressing capacity of ice mélanges (Carr and others, 2013; Bevan and others, 2019; Joughin and others, 2020), and have been correlated to glacier retreat in southeastern Greenland (Jensen and others, 2016), implying that rising SST around the GrIS should be given equal attention as the increasing presence of warm (subsurface) Atlantic Waters suspected to play a crucial role in SMM (Straneo and others, 2013; Inall and others, 2014). Other local factors such as katabatic winds (Christoffersen and others, 2011) and fjord geometry (Robel, 2017) have also been proposed as contributing factors to ice mélange rigidity, stressing both the complexity of ice mélange formation and the need for ice mélange models (Robel, 2017; Cassotto and others, 2021) to be coupled with ice flow and fjord circulation processes to better describe the future response of GrIS outlet glaciers. Current parameterisation of glacier terminus positions employed in ISMIP6 (Nowicki and others, 2016) utilises SMM that varies with future meltwater discharge and ocean thermal conditions. Because these factors likely also affect mélange formation and rigidity, disentangling the interwoven processes at the ice-ocean interfaces along the outlet glacier termini remains challenging, but necessary to improve predictions of the future behaviour of Greenland outlet glaciers.

Combining a decreasing/skipped IMB with an increasing SMM at KG produces an exaggerated retreat of > 30 km (Fig. 3, bottom row), not observed in previous scenarios where individual changes in either IMB or SMM induced little retreat (with the exception of IMBx0s). This suggests that KG's sensitivity to SMM is dependent on the glacier's terminus position along the flowline. In its current state, on a flat bedrock, KG is impervious to increased SMM (Fig. 3, second row). Yet skipped ice mélange years draw the glacier into a section of bedrock which deepens by c. 100 m (Fig. 3, fourth row), increasing the surface area of the terminus to SMM and promoting calving on a terminus already subject to increased longitudinal stress. This is best illustrated when comparing scenarios IMBx0.5s and COMB1 (Fig. 4); where both scenarios follow the same pattern of behaviour until the skipped mélange years 30–33, where rising summer SMM in COMB1 (now peaking at 7 m d^-1) induces a further retreat of 3 km to a small pinning point rising up c. 50 m at 324 km along the flowline. While the returning IMB in year 33 provides KG with conditions to re-advance in scenario IMBx0.5s, the glacier is found to consistently return to the pinning point in scenario COMB1 and undergoes further retreat when the IMB is skipped again in year 40. This varying influence of SMM aligns with findings from a 2D plan view modelling study of Store Glacier, West Greenland by Morlighem and others (2016). Their findings highlighted how summer SMM rates of 3 m d^-1 fail to induce retreat as the glacier rests at the front of a large shallow stabilising sill. Yet when the terminus is artificially retreated to the edge of the sill, a top of large over deepening in bedrock topography, a SMM of 3 m d^-1 induces a large scale retreat.

All scenarios that reach the pinning point at c. 324 km find their retreat interrupted or halted (Fig. 3), as the rise in bedrock topography offers stability through suppressed longitudinal stress on the stoss side of the pinning point (Benn and others, 2007; Favier and others, 2012; Todd and Christoffersen, 2014). Such pinning points are known to allow glaciers to maintain a stable terminus position for decades; modelling studies of Store Glacier, West Greenland, by Todd and Christoffersen (2014) and Rignot and others (2016) find the glaciers’ peerless stability since the 1960s related to a substantial sill upstream of the terminus, and KG's delayed retreat from its LIA maximum extent is likely explained by the stability offered by a high rising moraine (Fig. 1e) (Vermassen and others, 2020). Yet, retreat from such pinning points into a reverse sloped section of bedrock leaves a glacier susceptible to runaway mass loss through a positive feedback of terminus retreat and increased ice thickness and ice flux at the grounding line (Weertman, 1974; Schoof, 2007; Joughin and Alley, 2011; Schoof, 2012).

Such runaway mass loss occurs in all COMB scenarios (and IMBx0s), as KG retreats into a reverse sloped section of KF which deepens by c. 250 m causing the average calving extent to double in size (Fig. 5) and producing a retreat of c. 25 km in 3 years (Fig. 6). The same pattern of behaviour for each scenario (Fig. 6), indicates that the nature of KF's bedrock topography is crucial in controlling the long-term magnitude and timing of retreat: with a grounding line on reverse bedrock, exaggerated mass loss and thinning produces a final configuration of KG that is of stark contrast to scenarios where the glacier maintains its position on flat bedrock topography (Fig. 7). However, mapped submarine landforms such grounding zone wedges indicate that grounding lines may be stable even on reverse slopes (Jamieson and others, 2012), highlighting shortcomings of a 2D marine ice-sheet instability framework when contrasted with 3D geological evidence. Indeed, Gudmundsson (2013) stresses that the stability of grounding lines is a 3D problem, where horizontal buttressing forces can aid in reducing the ice flux as ice thickness increases. While we have accounted for the 3D convergence of ice in our flowline model we have not considered the effects of lateral friction, which may provide grounding line stability and control the otherwise runaway retreat.

Despite our model limitations, the extent and pattern of exceptional retreat observed at KG is in line with similar 2D flowline investigations that have account for the effects of lateral friction and forcings such as SMM and sea ice, where the latter can be argued to be similar to IMB. Both Nick and others (2013) and Beckmann and others (2019) found the large overdeepening in KF critical in producing a similar retreat extent as seen in COMB2 and COMB3 (Fig. 3) when exposing KG to a RCP 8.5 scenario until the years 2100 and 2200, respectively. While KG's annual retreat is not resolved in Beckmann and others (2019), Nick and others (2013) details rapid uncontrolled extensive retreat over 2–3 years when KG retreats into the large overdeepening section of bedrock topography, supporting our findings here. A caveat concerns modelled retreat rates which we find to be higher (30 km over 50 years) than Nick and others (2013) and Beckmann and others (2019) (30 km over 100 years). This may be related to our model limitations, in particular the lack of KG's stabilisation by lateral drag. Also the chosen bedrock topography dataset may also impact retreat rates, with Nick and others (2013) employing a dataset which captures an additional submarine sill upstream of KG's terminus, a feature not found in the Bedmachine v3 dataset (Morlighem and others, 2017) used in this study (Fig. 1b) and by Beckmann and others (2019), which likely provides KG with further stability and delays retreat towards the large overdeepening section of the fjord.

6. Conclusion

Here we have presented results from investigation into ice dynamics at Kangerlussuaq Glacier, east Greenland, utilising a 2D flowline model within Elmer/Ice. The novelty of our model experiments lies in considering the effect of skipped ice mélange backstress, which is found to be vital in initiating KG's retreat from flat bedrock towards the overdeepening section of bedrock in KF. Using a variety of different forcing scenarios applied over a 50 year period to gain insight into the glacier's state in 2065 we find that:

• SMB forced by an RCP 8.5 scenario has little to no influence on ice dynamics of KG.

• In its current position on flat bedrock topography, and with the consistent buttressing force of a winter ice mélange, KG is invulnerable to the predicted SMM rates of 8–12 m d^-1 expected by the end of the century (Slater and others, 2019).

• The buttressing pressure of a winter ice mélange is crucial to maintaining KG's stable position in KF. The mélange prevents winter calving and provides KG with conditions to re-advance into KF and counter any summer retreat. Winter calving only begins to occur when then IMB decreases below 90 kPa per metre of mélange, supporting existing evidence for this threshold in literature (Amundson and others, 2010; Krug and others, 2015).

• Subjecting KG to a skipped ice mélange year creates an extended calving season and produces a retreat c. 5 km from the glacier's previous winter maximum. Yet successive years of no buttressing mélange fail to induce further retreat. Such behaviour aligns with observed at KG during the retreat events in 2004/05 and 2016–18.

• Coupling skipped mélange years with rising SMM produces an excessive retreat of > 30 km. A failed ice mélange primes KG to retreat towards a reverse sloped section of bedrock where the glacier becomes increasingly vulnerable to submarine ablation that would otherwise render little impact.

• While the presence and rigidity of an ice mélange at KG hold a key influence on ice dynamics at KG in the near future, the long-term timing and magnitude of retreat is bound by the extensive over-deepenings of c. 250 m in KF that give KG the capacity to undertake rapid and uncontrolled mass loss.