1. Introduction
The Greenland Ice Sheet (GrIS) has contributed 0.70 ± 0.15 mm a−1 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, Reference Chen2017; The IMBIE Team, 2018), but also a lower bound for the rates to be observed in decades to come (Oppenheimer and others, Reference Oppenheimer2019; Briner and others, Reference Briner2020; Goelzer and others, Reference Goelzer2020).
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, Reference Enderlin2014; Fürst and others, Reference Fürst, Goelzer and Huybrechts2015; Mouginot and others, Reference Mouginot2019; Choi and others, Reference Choi, Morlighem, Rignot and Wood2021).
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, Reference Mouginot2019), 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, Reference Enderlin2014). In relation to rising air temperatures over the GrIS, SMB processes accounted for 84% of total mass loss since 2009 (Enderlin and others, Reference Enderlin2014), 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, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017; Rückamp and others, Reference Rückamp, Neckel, Berger, Humbert and Helm2019), Jakobshavn Isbræ began to readvance in 2016 (Khazendar and others, Reference Khazendar2019), and Kangerlussuaq Glacier underwent a further episode of acceleration and retreat in the same year (Brough and others, Reference Brough, Carr, Ross and Lea2019).
Observations at Kangerlussuaq Glacier (KG) show interannual variations in both the presence and strength of an ice mélange (Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). 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, Reference Amundson2010; Foga and others, Reference Foga, Stearns and van der Veen2014; Moon and others, Reference Moon, Joughin and Smith2015; Burton and others, Reference Burton, Amundson, Cassotto, Kuo and Dennin2018; Xie and others, Reference Xie, Dixon, Holland, Voytenko and Vankova2019). 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, Reference Morlighem2017) (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, Reference Murray2010). 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, Reference Christoffersen2011; Inall and others, Reference Inall2014; Bevan and others, Reference Bevan, Luckman, Benn, Cowton and Todd2019). 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, Reference Mouginot2019). At KG's terminus, velocities regularly exceed 10,000 m a-1 (Fig. 1c) (Joughin and others, Reference Joughin2008b; Murray and others, Reference Murray2010; Brough and others, Reference Brough, Carr, Ross and Lea2019). 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, Reference Kjeldsen2015; Batchelor and others, Reference Batchelor, Dowdeswell, Rignot and Millan2019) (Figs 1d, e), occurred shortly after its discovery in 1930 (Watkins, Reference Watkins1932; Koch, Reference Koch1933), and was likely triggered by increasing ocean temperatures, but subsequently controlled by bathymetry (Vermassen and others, Reference Vermassen2020).
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, Reference Thomas2000; Khan and others, Reference Khan2014). 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, Reference Luckman, Murray, de Lange and Hanna2006; Howat and others, Reference Howat, Joughin and Scambos2007; Joughin and others, Reference Joughin2008a). This retreat was attributed to markedly increased ocean temperatures (Cowton and others, Reference Cowton2016; Millan and others, Reference Millan2018), coinciding with extensive retreat and acceleration of many glaciers along Greenland's east coast (Seale and others, Reference Seale, Christoffersen, Mugford and O'Leary2011). 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, Reference Joughin2008b; Christoffersen and others, Reference Christoffersen2011; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017; Bevan and others, Reference Bevan, Luckman, Benn, Cowton and Todd2019; Brough and others, Reference Brough, Carr, Ross and Lea2019).
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, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017; Bjørk and others, Reference Bjørk2018). 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, Reference Weertman1974; Joughin and Alley, Reference Joughin and Alley2011). 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, Reference Brough, Carr, Ross and Lea2019) and a rapid acceleration in ice discharge (Hansen and others, Reference Hansen2021). When an ice mélange materialised again, in the winter 2018/19, KG began to readvance into KF (Bevan and others, Reference Bevan, Luckman, Benn, Cowton and Todd2019; Brough and others, Reference Brough, Carr, Ross and Lea2019).
3. Methods
The finite-element model Elmer/Ice (Gagliardini and others, Reference Gagliardini2013) 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, Reference Greve and Blatter2009), here combined with the crevasse depth criterion (Nye, Reference Nye1957; Benn and others, Reference Benn, Warren and Mottram2007; Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2020; Otero and others, Reference Otero, Navarro, Martin, Cuadrado and Corcuera2010; James and others, Reference James, Murray, Selmes, Scharrer and O'Leary2014; Todd and Christoffersen, Reference Todd and Christoffersen2014) to model calving. Full Stokes ice dynamics including calving can be modelled in 3D (Todd and others, Reference Todd2018), 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 (Reference Mouginot and Rignot2019)) and embedded in the region where the highest ice velocities were observed from the MEaSUREs project (Joughin, Reference Joughin2020) (Fig. 1c). Ice thickness and bed topography along the flowline were extracted from the Bedmachine v3 dataset (Morlighem and others, Reference Morlighem2017) and used to render a closed 2D domain (Fig. 1b) within which a mesh was produced using GMSH (Geuzaine and Remacle, Reference Geuzaine and Remacle2009). 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, Reference Cook2014) 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 (Reference Todd and Christoffersen2014) and Passalacqua and others (Reference Passalacqua2016). 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, Reference Weertman1974; Gagliardini and others, Reference Gagliardini2013). Whether KG is grounded or floating is determined by using a height-above-buoyancy criteria (Nick and others, Reference Nick, Van Der Veen, Vieli and Benn2020; Favier and others, Reference Favier, Gagliardini, Durand and Zwinger2012) 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 (Reference Funk, Echelmeyer and Iken1994) 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, Reference Hall2019). 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 (Reference Gagliardini2013), 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, Reference Luckman2015):
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 (Reference Nye1957):
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 (Reference Nick, Van Der Veen, Vieli and Benn2020):
where p w is water pressure. The crevasse depth criterion, first proposed by Benn and others (Reference Benn, Warren and Mottram2007) and numerically investigated by Nick and others (Reference Nick, Van Der Veen, Vieli and Benn2020), 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 (Reference Todd2018) in 3D, and by Todd and Christoffersen (Reference Todd and Christoffersen2014) 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, SMBref and SMBRCP8.5. Both are based on 5 km spatial resolution SMB maps of Greenland, derived from the regional climate model MAR (Fettweis and others, Reference Fettweis2017). For SMBref, MAR is forced with ERA-Interim reanalysis data (Dee and others, Reference Dee2011) 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). SMBRCP8.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, Reference Watanabe2010) for RCP pathway 8.5 (Fig. 2a).
Further, four different SMM forcings are used, SMMref, 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, Reference Holmes, Kirchner, Kuttenkeuler, Krützfeldt and Noormets2019), 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, Reference Cowton, Slater, Sole, Goldberg and Nienow2015; Carroll and others, Reference Carroll2016; Fraser and others, Reference Fraser, Inall, Magaldi, Haine and Jones2018). As SMM are dependent on water temperatures and subglacial discharge (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003; Xu and others, Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013), we vary SMM rates in SMMref 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 (Reference Fraser, Inall, Magaldi, Haine and Jones2018) 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 (Reference Carroll2016). 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, Reference Carroll2016), SMMref 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 SMMref, 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, Reference Slater2020) 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, IMBref, 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, Reference Joughin2008a; Christoffersen and others, Reference Christoffersen2011; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017), 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, Reference Walter2012), where IMB is reported to range from 30 to 60 kPa. With an IMB of 45 kPa, Todd and Christoffersen (Reference Todd and Christoffersen2014) were able to reproduce the seasonal dynamics observed at Store Glacier in a 2D flow line model. Similarly, Vieli and Nick (Reference Vieli and Nick2011) 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 IMBref. This is then used to produce a per metre of mélange backstress, σ M, of 198 kPa following Todd and Christoffersen (Reference Todd and Christoffersen2014):
where H T.Avg is the average height of the terminus, 597 m (Morlighem and others, Reference Morlighem2017), 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, Reference Cook2014). 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 (Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017) 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 SMBref, SMMref and IMBref, 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, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017; Bevan and others, Reference Bevan, Luckman, Benn, Cowton and Todd2019). 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, Reference Joughin2020), 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.
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).
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.
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.
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.
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.
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 (Reference Nick2013) and Beckmann and others (Reference Beckmann2019) 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, Reference Fettweis2013; Noél and others, Reference Noél, van Kampenhout, Lenaerts, van de Berg and van den Broeke2021). 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, Reference Joughin, Tulaczyk, Fahnestock and Kwok1996; Zwally and others, Reference Zwally2002; Chudley and others, Reference Chudley2019; Yang and others, Reference Yang2019; Lampkin and others, Reference Lampkin, Koenig, Joseph and Box2020). Also, basal glacier mass balance is suspected to interact with SMB, and to augment overall mass loss (Karlsson and others, Reference Karlsson2021). 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, Reference Beaird, Straneo and Jenkins2018), 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, Reference Motyka, Hunter, Echelmeyer and Connor2003; O'Leary and Christoffersen, Reference O'Leary and Christoffersen2013; Ma and Bassis, Reference Ma and Bassis2019). By 2065, SMM rates at KG are expected to have doubled to 8 m d-1 (Slater and others, Reference Slater2020), 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, Reference Slater2020)) 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, Reference Carroll2016), 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, Reference Todd and Christoffersen2014; Cook and others, Reference Cook2014; Krug and others, Reference Krug, Durand, Gagliardini and Weiss2015), 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, Reference Slater2019) used to parameterise ice-ocean dynamics in the ISMIP6 project (Nowicki and others, Reference Nowicki2016), 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, Reference Joughin2008a, Reference Joughinb; Amundson and others, Reference Amundson2010; Moon and others, Reference Moon, Joughin, Smith and Howat2012; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017; Amundson and Burton, Reference Amundson and Burton2018). 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 (Reference Krug, Durand, Gagliardini and Weiss2015), where a buttressing force between 130 and 90 kPa was found sufficient to impede calving. Similarly, Amundson and others (Reference Amundson2010) 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 (Reference Walter2012)'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, Reference Joughin2008a; Christoffersen and others, Reference Christoffersen2011) and winters 2016/17 and 2017/18 (Bevan and others, Reference Bevan, Luckman, Benn, Cowton and Todd2019; Brough and others, Reference Brough, Carr, Ross and Lea2019). 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, Reference Bevan, Luckman, Benn, Cowton and Todd2019). 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, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017) and Jakobshavn Isbræ (Amundson and others, Reference Amundson2010; Joughin and others, Reference Joughin, Shean, Smith and Floricioiu2020), 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, Reference Howat, Box, Ahn, Herrington and McFadden2010; Carr and others, Reference Carr, Vieli and Stokes2013; Moon and others, Reference Moon, Joughin and Smith2015). Sea surface temperatures (SST) are suggested to have a large influence on the rigidity and buttressing capacity of ice mélanges (Carr and others, Reference Carr, Vieli and Stokes2013; Bevan and others, Reference Bevan, Luckman, Benn, Cowton and Todd2019; Joughin and others, Reference Joughin, Shean, Smith and Floricioiu2020), and have been correlated to glacier retreat in southeastern Greenland (Jensen and others, Reference Jensen, Box and Hvidberg2016), 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, Reference Straneo2013; Inall and others, Reference Inall2014). Other local factors such as katabatic winds (Christoffersen and others, Reference Christoffersen2011) and fjord geometry (Robel, Reference Robel2017) 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, Reference Robel2017; Cassotto and others, Reference Cassotto, Burton, Amundson, Fahnestock and Truffer2021) 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, Reference Nowicki2016) 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 (Reference Morlighem2016). 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, Reference Benn, Warren and Mottram2007; Favier and others, Reference Favier, Gagliardini, Durand and Zwinger2012; Todd and Christoffersen, Reference Todd and Christoffersen2014). 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 (Reference Todd and Christoffersen2014) and Rignot and others (Reference Rignot2016) 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, Reference Vermassen2020). 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, Reference Weertman1974; Schoof, Reference Schoof2007; Joughin and Alley, Reference Joughin and Alley2011; Schoof, Reference Schoof2012).
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, Reference Jamieson2012), highlighting shortcomings of a 2D marine ice-sheet instability framework when contrasted with 3D geological evidence. Indeed, Gudmundsson (Reference Gudmundsson2013) 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 (Reference Nick2013) and Beckmann and others (Reference Beckmann2019) 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 (Reference Beckmann2019), Nick and others (Reference Nick2013) 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 (Reference Nick2013) and Beckmann and others (Reference Beckmann2019) (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 (Reference Nick2013) 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, Reference Morlighem2017) used in this study (Fig. 1b) and by Beckmann and others (Reference Beckmann2019), 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, Reference Slater2019).
-
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, Reference Amundson2010; Krug and others, Reference Krug, Durand, Gagliardini and Weiss2015).
-
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.
Acknowledgements
The work was partially supported by the Swedish Research Council FORMAS under grant 2017-00665 awarded to N.K. The numerical simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre at Linköping University, Sweden. We thanks two anonymous reviewers, whose constructive feedback helped to improve the manuscript, and provided inspiration for future work.