Terminus thinning drives recent acceleration of a Greenlandic lake-terminating outlet glacier

Abstract

Ice-contact proglacial lakes affect ice dynamics and the transition of glacier termini from land- to lake-terminating has been shown to cause ice flow acceleration. In recent decades, the number and size of Greenlandic ice-marginal lakes has increased, highlighting the need to further understand these lake-terminating ice-margins as their influence on ice sheet mass balance increases. Here, time series of satellite-derived observations of ice velocity, surface elevation, and terminus position were generated at a lake-terminating outlet glacier, Isortuarsuup Sermia, and the nearby land-terminating Kangaasarsuup Sermia in south-west Greenland. At Isortuarsuup Sermia, annual surface velocity at the terminus increased by a factor of 2.5 to 214 ± 4 m yr⁻¹ (2013–2021), with the magnitude of this acceleration declining with distance up-glacier. Meanwhile, near-terminus surface elevation changed at a rate of −2.3 ± 1.1 m yr⁻¹ (2012–2021). Conversely, velocity change at Kangaasarsuup Sermia was minimal, while surface elevation change was approximately half at comparable elevations (−1.2 ± 0.3 m yr⁻¹). We attribute these dynamic differences to thinning at Isortuarsuup Sermia and subsequent retreat from a stabilising sublacustrine moraine, and emphasise the potential of proglacial lakes to enhance future rates of mass loss from the Greenland Ice Sheet.

1. Introduction

Rates of mass loss from the Greenland Ice Sheet (GrIS) increased six-fold between the 1980s and 2018 (Mouginot and others, 2019), raising sea levels by 10.8 ± 0.9 mm (1992–2018) (The IMBIE Team, 2020). The GrIS is projected to continue losing mass, and estimates of GrIS sea level rise (SLR) contributions vary with emissions scenario. By 2100 an additional 70 ± 40 mm of SLR is projected under RCP2.6, increasing to between 80 and 270  mm under RCP8.5 (Fox-Kemper and others, 2021). Refining these SLR estimates requires greater understanding of the controls on ice sheet mass loss.

The GrIS is currently losing mass via both surface and dynamic processes (The IMBIE Team, 2020), with dynamics being, in part, dependent on surface processes. For example, long-term negative surface mass balance induced thinning (negative surface elevation change) may lead to acceleration at lake or marine-terminating margins if thinning causes greater reductions in resistive stresses than driving stresses (Pfeffer, 2007). Furthermore, at lake- and marine-terminating outlets, acceleration and surface lowering may be self-sustaining if accompanied by retreat into deeper water and further dynamic thinning (Weertman, 1974; Meier and Post, 1987; O'Neel and others, 2005; Pfeffer, 2007).

At present, however, the influence of proglacial lakes on ice mass loss and SLR are either absent or poorly represented in ice sheet models (Carrivick and others, 2020). Since proglacial lakes modify ice dynamics, by altering terminus profile, subglacial hydrology and local force balance (e.g. Warren and Kirkbride, 2003; Sugiyama and others, 2011; Baurley and others, 2020), there is a need to determine the extent to which proglacial lakes will impact ice sheet mass loss over the coming century and beyond (Carrivick and others, 2022).

There are many subglacial bedrock overdeepenings beneath the GrIS (Patton and others, 2016; Morlighem and others, 2017) which fill with meltwater runoff during margin recession, forming ice-marginal proglacial lakes (Costa and Schuster, 1988; Carrivick and Tweed, 2013). In recent decades, the area and number of ice-marginal lakes has increased in south-west Greenland, as well as globally (Carrivick and Quincey, 2014; Shugar and others, 2020; How and others, 2021; Rick and others, 2022). Such changes in ice margin configuration are significant because outlet glaciers that terminate in lakes typically have greater rates of mass loss and terminus retreat than their land-terminating counterparts (Kirkbride, 1993; Warren and Kirkbride, 2003; Schomacker, 2010; Tsutaki and others, 2011; Carr and others, 2017; King and others, 2018; Mallalieu and others, 2021). This change in dynamics reflects differences in the boundary conditions at (a) the bed, and (b) the terminus (Pronk and others, 2021). The presence of a lake leads to a reduction in the effective pressure at the terminus, and up-glacier, enabling greater rates of basal sliding (Bindschadler, 1983; Benn and others, 2007; Sugiyama and others, 2011). Effective pressure is the difference between ice-overburden pressure and basal water pressure, hence thinner ice at the terminus and or deeper lake water will promote greater ice velocities (Kirkbride and Warren, 1997; Tsutaki and others, 2013). Lake depth at the terminus sets the base level and thus minimum basal water pressure for the near-terminus subglacial hydraulic system, and the influence of the pressure head due to the lake declines with distance from the terminus at a rate dependent on the bed topography (Meier and Post, 1987; Benn and others, 2007). Water at the terminus also initiates a suite of complementary processes: sub-aqueous melt, thermal-notch erosion, and calving (both sub-aerial and sub-aqueous) (Röhl, 2006; Sugiyama and others, 2011, 2019; Mallalieu and others, 2020). Water-depth impacts terminus buoyancy with calving rates increasing as water depth increases (Benn and others, 2007; Boyce and others, 2007; Dykes and others, 2011). These processes may contribute to terminus retreat, surface steepening, and further increases in velocity and longitudinal strain rates (e.g. Warren and Kirkbride, 2003; Tsutaki and others, 2013; King and others, 2018). Furthermore, modelling results suggest that glacier response to the development of proglacial lakes may be partially decoupled from short-term changes in climate and contribute to rapid and sustained retreat (Sutherland and others, 2020).

The effect of proglacial lakes, and their development, on ice dynamics has been observed in many glaciated regions. For example, at Breiðamerkurjökull, Iceland (which has both lake- and land-terminating distributaries), between 1991–2015 there was no change in ice velocity adjacent to the land-terminating margins, whereas velocity increased by a factor of three to 3.5  m day⁻¹ proximate to the terminus of the lake-terminating arm (Baurley and others, 2020). This change was linked to increases in surface air temperature initiating terminus retreat into a 200–300 m deep subglacial trough, triggering a positive feedback mechanism (Pfeffer, 2007; Nick and others, 2009). Retreat into deeper water enabled ice flow acceleration, dynamic thinning (in addition to thinning from changes in surface mass balance), greater rates of calving, and further retreat into deeper water.

The contrasting patterns in ice dynamics between lake- and land-terminating glaciers have also been observed in the Himalaya (e.g. King and others, 2019; Tsutaki and others, 2019; Pronk and others, 2021). For example, greater rates of mass loss and terminus retreat have been observed at lake-terminating outlets (King and others, 2019), and their centreline velocities are typically double those measured at land-terminating glaciers (18.83 vs. 8.24 m yr⁻¹, for the period 2017–2019) (Pronk and others, 2021). One notable difference between lake- and land-terminating glaciers is the centreline velocity profile, whereby lake-, like marine-terminating outlet glaciers, accelerate toward the terminus (Tsutaki and others, 2019; Pronk and others, 2021). This extensional flow contributes to surface thinning, whereas the compressive flow regime at land-terminating glaciers can lead to thickening, which may offset surface mass balance induced surface lowering (Tsutaki and others, 2013, 2019). Additionally, glacier geometry influences both terminus stability and velocity, with valley constrictions or submerged sills and their associated impacts on lateral and back-stress causing reductions in velocity and enhanced terminus stability (van der Veen and Whillans, 1989; O'Neel and others, 2005; Benn and others, 2007). These observations are supported by sensitivity modelling which suggests that thicker ice, a wider terminus and steeper surface slopes lead to elevated near-terminus velocities at lake-terminating glaciers (Pronk and others, 2021).

Lake-terminating glacier dynamics can also be affected by surface meltwater input to the glacier bed (Sugiyama and others, 2011). For example, at Glacier Perito Moreno, a lake-terminating glacier in Patagonia, Argentina, hourly variations in measured basal water pressures 4 km from the terminus correspond closely with changes in surface temperature and ice velocity measured in-situ using dGPS (Sugiyama and others, 2011). These observations suggest surface meltwater can reach the bed rapidly, and that ice velocity is sensitive to small changes in basal water pressures. This sensitivity is evidenced by the observed difference in relative change: over a ten day period ice velocity varied by 37% about its mean (1.43 m day⁻¹), whereas basal water pressures only varied by 5%. This corresponded to an increase in velocity of 0.053 m day⁻¹ per 1 $^\circ$C (Sugiyama and others, 2011).

In addition to the above geometric and climatic controls, the characteristics of individual proglacial lakes will influence the behaviour of lake-terminating glaciers (e.g. Sugiyama and others, 2016; Mallalieu and others, 2020; Watson and others, 2020; Dye and others, 2021). For example, a reduction in calving frequency and volume corresponds to the timing of lake-ice freeze up (Mallalieu and others, 2020). Additionally, sub-aqueous melt and associated thermal-notch erosion are both functions of the thermal structure of the lake (e.g. Haresign and Warren, 2005; Röhl, 2006; Minowa and others, 2017). Observations in Patagonia (Sugiyama and others, 2016, 2019) revealed a layer of cold turbid water derived from subglacial discharge underlying warmer surface waters. This stratification can prohibit the upwelling of meltwater that is seen in glacial fjords, and allows for the formation of ice terraces below the waterline (Kirkbride and Warren, 1997; Sugiyama and others, 2019). The thermal state of a proglacial lake is strongly coupled to climate, and is dependent on incident shortwave radiation, surface air temperatures, winds, precipitation and runoff (e.g. Schomacker, 2010; Richards and others, 2012).

The observed impacts of proglacial lakes on ice dynamics (Kirkbride, 1993; Warren and Kirkbride, 2003; Sugiyama and others, 2011; Baurley and others, 2020; Pronk and others, 2021) are furthermore supported by modelling studies (e.g. Sutherland and others, 2020). Collectively, these works suggest the presence of proglacial lakes leads to greater rates of terminus retreat and mass loss than those at land-terminating glaciers, thereby contributing to accelerated rates of deglaciation.

Given the clear potential of lakes to perturb ice dynamics (e.g. Kirkbride, 1993), and the projected prevalence of ice-marginal lakes in Greenland (Carrivick and others, 2022), it is important to evaluate how, against a backdrop of a warming climate, lakes are impacting ice motion and terminus positions. In south-west Greenland mean annual changes in margin position have varied since the 1980s, with notable differences between the lake- and land-terminating sectors. Average annual rates of margin recession increased by an order of magnitude from 1.1  m yr⁻¹ (1987–1992) to 11.5  m yr⁻¹ (2010–2015) along lacustrine margins, whereas the magnitude of changes at terrestrial margins was more modest: from advance of 1.2  m yr⁻¹ to recession of 2.8  m yr⁻¹ (Mallalieu and others, 2021). Furthermore, observations from across Greenland indicate ice-marginal lakes enhance the flow of adjacent ice by ~25% (Carrivick and others, 2022).

This study aims to investigate recent (2013–2021) changes in ice velocity, surface elevation, and terminus retreat at two proximate but contrasting outlet glaciers in south-west Greenland; the lake-terminating Isortuarsuup Sermia, and the nearby land-terminating Kangaasarsuup Sermia. We additionally consider the varying processes that impact dynamics at this lake-terminating system, and the potential significance of these at the ice-sheet scale.

2. Data and methodology

2.1 Study site

We generated ice velocities, rates of surface elevation change, and terminus position from two outlet glaciers in south-west Greenland: Isortuarsuup Sermia (IS) (63$^\circ 50'$ N, 49$^\circ 59'$ W) and Kangaasarsuup Sermia (KS) ($64^\circ 07'$ N, $49^\circ 54'$ W) (Fig. 1). Isortuarsuup Sermia is a lake-terminating glacier that drains into one of the largest proglacial lakes in south-west Greenland, the ~60 km² Isortuarsuup Tasia, and KS is a nearby land-terminating glacier. These outlet glaciers were selected due to their close proximity (~30 km) and their similar morphological characteristics including terminus elevation (~315 m at KS and ~500 m at IS), surface slope, and valley width. Based on meltwater routing via subglacial hydraulic potential (Shreve, 1972), KS drains 660  km², whereas IS drains 122 km² (Mankoff, 2020) (Fig. 1a).

Figure 1. (a) The lake-terminating Isortuarsuup Sermia (IS) and land-terminating Kangaasarsuup Sermia (KS) with their respective centrelines (solid lines) and runoff catchments from Mankoff and others (2020) (dashed lines). Ice surface contours at 250 m intervals generated from BedMachine v5 (Morlighem and others, 2022); background image: Sentinel-2 acquisition from 19 September 2022 (ESA Copernicus, 2022); Location of study area, south-west Greenland, shown by pink box in inset (upper right). Dashed red box at terminus of IS denotes extent of (b); (b) terminus region of IS illustrating the trim-line (white arrows), terminal moraine (black arrow), grounded icebergs (orange arrow). Sentinel-2 acquisition from 18 September 2019 (ESA Copernicus, 2022).

There is clear evidence of a distinct terminal moraine at IS (black arrow in Fig. 1b), which likely formed during a period of prolonged stability, as evidenced by the pronounced trim-line (Fig. 1b) inferred to be of Little Ice Age origin (during the 19^th century) by Weidick and others (2012). Furthermore, icebergs are often grounded in the lake approximately 400 m from the glacier terminus suggesting a sublacustrine extension of the visible terminal moraine (Fig. 1b).

2.2 Ice velocity

Ice velocities were obtained from the NASA MEaSUREs ITS_LIVE version 2 data-cubes (Gardner and others, 2018, 2022). This data set is derived from optical (Landsat 8 and Sentinel-2A/B) and radar acquisitions (Sentinel-1A/B). Velocities are determined using the autonomous repeat image feature tracking algorithm applied to pairs of overlapping images from a given sensor, separated by time (date_dt) (Gardner and others, 2018; Lei and others, 2021). Velocity fields obtained from image pairs with small time separations (date_dt ≤ 30 days) reveal short-term changes in velocity while long time separations (date_dt ≥ 300 days) provide better estimates of annual averages. Calculated velocities are posted onto a uniform 120 m grid, with a spatially variable effective resolution of 240–1920 m (Lei and others, 2021). The resulting data-cube has dimensions easting (x), northing (y), and time (t), where t is the mid-date between the two satellite acquisitions used to generate the velocity field. Each grid cell contains the velocity component in the x (v _x) and y (v _y) directions, and image-pair time dependent error estimates ($\sigma _{v_t}$) are also supplied. To minimise the effects of point sampling, and to allow for the spatially variable effective resolution, which is a function of the search window size used when feature tracking (Lei and others, 2022), the v _x and v _y fields were first spatially averaged using a 3 × 3 window, and subsequently sampled every 250  m along the glacier centrelines, and the resultant velocity calculated (Eqn (1)):

(1)$$v_t = \sqrt{\bar{v_{x_t}}^2 + \bar{v_{y_t}}^2}$$

Error-weighted average velocities ($\bar {v}$) were determined (Eq. (2)):

(2)$$\overline{v} = {\sum_{t}{{\bar{v_t}}/ \sigma_{v_t}^2}\over \sum_{t}{1/\sigma_{v_t}^2}}$$

With the uncertainty in $\bar {v}$ calculated as (Eq. (3)):

(3)$$\sigma_{v} = \sqrt{{1\over \sum_{t}{1/\sigma_{v_t}^2}}}$$

The error-weighted average annual velocities were calculated using all velocity fields derived from image-pairs separated by between 300 and 430 days, with assignment to a given year on the basis of the mid-date of the image-pair. Annual average velocity profiles were constructed for 2013–2021, and the velocity trends computed from these averages for each point along the centreline using linear regression. To further assess differences in velocity regime, seasonal averages were constructed from velocity fields where date_dt ≤ 30 days, and rates of acceleration along the centreline determined for each season. Seasons were defined as: winter (December, January, February), spring (March, April, May), summer (June, July, August), and autumn (September, October, November). Due to the limited availability of velocity fields with small image separations prior to 2016, the seasonal analysis presented here is confined to the period 2016–2021. For both the annual and seasonal trend analysis, the null hypothesis that there is no trend (i.e. a regression slope coefficient of zero), was evaluated using a two-tailed Wald test, as implemented in the scipy python package. Following convention, when the returned p-value was ≤ 0.05, the trends are taken to be significant.

2.3 Surface elevation change

Rates of surface elevation change were derived from ArcticDEM 4.1 (Porter and others, 2022). ArcticDEM is a collection of high resolution (2 m) time dependent (2007–2021) digital elevation models (DEMs). These are constructed from stereo auto-correlation methods (Noh and Howat, 2015), applied to sub-metre resolution optical satellite imagery from the Maxar constellation. The absolute accuracy of individual ArcticDEM strips is approximately 4 m both horizontally and vertically (Porter and others, 2022). The volume of data processed in the construction of ArcticDEM allows the accidental inclusion of errors in the data set (Błaszczyk and others, 2019). Despite this, once DEMs have been co-registered, accuracy has been shown to be, in many cases, better than 4 m and with precision of approximately 1  m, making ArcticDEM suitable for measuring changes ≥ 1 m (Błaszczyk and others, 2019).

To co-register ArcticDEMs, a 5 km buffer was constructed for each glacier centreline, and all ArcticDEM strips that intersected this region were selected. To avoid detecting changes in elevation as a function of seasonal snow cover, the list of DEMs was filtered to include only those constructed from images acquired in June, July, August and September. After this filtering, there were 20 DEMs spanning the period September 2012–June 2021 over the study region at IS, and 51 (June 2011– September 2021) at KS. The supplied bit-mask was applied to each DEM to preserve only those pixels marked as good data; cloud, water and edge pixels were masked. All DEMs were visually inspected and the DEM with the best coverage of each glacier was selected to be the reference DEM, to which all others were co-registered. At IS the reference DEM was from 15/06/2016; at KS 04/08/2014. A cloud-free Sentinel-2 scene (18/08/2022 at IS, and 24/07/2022 at KS) was used to identify regions of stable terrain (SI. 1). During the co-registration process, differences in elevation over these regions of stable terrain are minimised, enabling changes in ice surface elevation to be observed accurately. Each DEM was co-registered to the reference using the method proposed by Nuth and Kääb (2011), as implemented in the python package XDEM (Xdem Contributors, 2021). Once co-registered each DEM was down-sampled from 2 m to 20 m using bi-linear interpolation, to reduce both file size and the effects of point sampling. To provide a measure of width-averaged rates of thinning, the co-registered DEMs were sampled every 50 m along a series of parallel lines spaced every 250 m. Rates of change were computed from these elevation samples using linear regression. As per the velocity trend analysis, significant trends were evaluated using a two-tailed Wald test. For each co-registered DEM a quadratic surface (Eq. (4)) was fitted to elevation values (z) within a 21 × 21 (420 × 420 m) window using least squares regression, and the surface slope (S) calculated from the coefficients (Eq. (5)) (e.g. Hurst and others, 2012). These slope surfaces were also sampled every 100 m along offset parallel lines, and the width-averaged (median) change in slope between the earliest and latest DEM calculated.

(4)$$z = Ax^2 + Bx + 2 + Cxy + Dx + Ey + F$$

(5)$$S = \sqrt{d^2 + e^2}$$

Uncertainty in the estimated rates of surface elevation change was minimised through the co-registration process. Prior to co-registration, the median difference in elevations over stable terrain ranged between −2.2 and 5.0 m at IS and −8.6 and 2.0 m at KS. After co-registration, the median differences in elevation over stable terrain were much closer to zero (−0.02 − 0.02 m, and −0.03 − 0.01 m, at IS and KS, respectively). The normalised median absolute deviation (NMAD), which measures a sample's dispersion, of elevation differences was similarly reduced in the co-registration process (SI. 2) At IS, NMADs were reduced from between 0.41–1.71 m to 0.24–0.69 m, with similar improvements at KS (0.34–2.41 m, to 0.23–1.97 m). We are therefore confident in our ability to detect changes in elevation ≥ 1 m, which is equivalent to annual rates of change of ~0.1 m yr⁻¹ over the period which ArcticDEM is available.

2.4 Terminus positions

Terminus positions were manually digitised from optical imagery captured by Landsat 8 and Sentinel-2 using the Google Earth Engine Digitization Tool (Lea, 2018) from 2013 to 2022. Relative changes in terminus position were determined using the rectilinear box method (Moon and Joughin, 2008) which allows for uneven retreat across the terminus. Uncertainty in relative terminus positions arise from image co-registration errors and manual digitisation errors (e.g. Carr and others, 2014). Image co-registration errors are a function of poor spatial alignment between satellite image acquisitions. The Landsat scenes used in this study had a median image registration accuracy of 4.6 m, and the Sentinel-2 technical reference indicates geolocation uncertainties are ≤ 11 m (S2 MSI ESL Team, 2022). To quantify the digitisation precision, termini were re-digitised five times (Paul and others, 2013), and the standard error was determined to be 16.7 m, which yields an uncertainty in rate of terminus position change of ±3.7 m yr⁻¹ over the study period, and is a lower bound on the measurable rate of terminus position change. At IS, digitising errors arose due to difficulties in discriminating between the terminus and either recently calved icebergs or lake ice cover. Digitising errors at KS were principally due to snow cover at the beginning of the melt season, debris cover at the end of the melt season, and shadows cast by the ridge to the south. As such, for KS, the results presented here are average relative terminus positions over each summer.

2.5 Runoff

The runoff data set used (Mankoff and others, 2020) comprises liquid water discharge estimates at hydrological outlets derived from two regional climate models: MAR (Modele Atmospherique Regional) and RACMO (Regional Atmospheric Climate Model) (Noël and others, 2016; Fettweis and others, 2017). The subglacial stream network is determined from ice surface elevations (Porter and others, 2018) and ice thickness estimates from BedMachine (Morlighem and others, 2017) using a model for the subglacial pressure head (Shreve, 1972). This data set accounts for surface melt, rainfall, meltwater retention and refreezing. Supraglacial flow is discounted and meltwater is assumed to generate and reach the bed within the same model grid cell (Mankoff and others, 2020). At the study site, the basin output closest to each glacier termini was selected and cumulative runoff calculated for the years 2011–2021. Linear regression was used to evaluate trends in cumulative annual runoff, for both MAR and RACMO, and tested for significance using a two-tailed Wald test. Similarly, to assess the relationship between runoff and ice velocity, average annual ice velocity was regressed against cumulative annual runoff (2013–2021) (derived from the mean average of MAR and RACMO).

This data set is not supplied with uncertainty estimates, however it contains three principal sources of uncertainty. (1) Temporal uncertainty which is a function of how the routing model handles the time lag between meltwater generation and runoff within each grid cell; (2) basin uncertainty, which arises from the data used in computing the subglacial stream network and catchments; and (3) uncertainties in the regional climate models from which the liquid water discharge estimates are derived (Mankoff and others, 2020). Here, temporal uncertainties are mitigated by computing cumulative annual totals.

3. Results

3.1 Ice velocity

Ice velocities show clear and contrasting patterns in behaviour at the two outlet glaciers between 2013 and 2021 (Fig. 2). Annual average ice velocity within 500 m of the terminus more than doubled from ~80 to 220 m yr⁻¹ at IS, while there was minimal change at KS (from ~25 to 20 m yr⁻¹) (Fig. 2). The magnitude of the acceleration at IS decreases from 15.0 ± 2.4 m yr⁻² at the terminus to 1.4 ± 0.5 m yr⁻² 15 km up-glacier (Fig. 2e); furthermore, the increase in velocity extends across the full width of the terminus (Fig. 3a). Notably, between 2018 and 2019, there was a substantial (~30%) increase in near terminus surface velocity at IS from 147 to 193  m yr⁻¹ (Fig. 2a & c). By contrast, barring a region of low magnitude (1.0–1.9 m yr⁻²) acceleration 2.5–8.5 km from the terminus, there is no discernible trend in ice velocity at KS (Figs. 2f and 4a). The differences in velocity magnitude and velocity trend between these two glaciers declines with distance from the terminus, such that at a distance of ~15 km from their respective termini, differences in the average annual velocity change over the study period at IS (from 127 to 141 m yr⁻¹) and KS (from 117 to 125  m yr⁻¹) are minimal (Figs. 2 and 3).

Figure 2. Annual average ice velocity (2013–2021) profiles along centrelines shown in Figure 1 at (a) Isortuarsuup Sermia, and (b) Kangaasarsuup Sermia. (c) and (d) show percentage change relative to 2013. (e) and (f) show linear trends where regression slope coefficients are significant at p $\leq 0.05$; error bars denote 95% confidence interval.

Figure 3. Isortuarsuup Sermia: (a) change in average annual velocity between 2013–2021; (b) rate of surface elevation change (September 2012– June 2021) from ArcticDEM (negative denotes thinning); manually digitised ice margin shown in black in (a) and (b); (c) terminus positions 2014–2021. Red box in (a) denotes extent of (c). White arrows in (c) indicates the Little Ice Age trim-line and the black arrow points to the associated terminal moraine with icebergs grounded on its sublacustrine extension indicated by the orange arrow.

Figure 4. Kangaasarsuup Sermia: (a) change in average annual velocity between 2013–2021; (b) rate of surface elevation change (June 2011– September 2021) from ArcticDEM (negative denotes thinning); manually digitised ice margin shown in black in (a) and (b); (c) terminus positions 2014–2021. Red box in (a) denotes extent of (c).

Seasonal cycles were observed at both glaciers, however, there were key differences between IS and KS with regards to seasonal velocity trends over the five year period 2016–2021 (Figs. 5 and 6). At IS, along the entire lower 10 km of the glacier, winter ice surface velocity increased, with the magnitude of acceleration declining from 20.2 (±12) m yr⁻² at the terminus to 4.4 (±4.3)  m yr⁻² at 10 km (Fig. 6). There were also significant positive trends in autumn and summer within the 2 km closest to terminus, although these were of slightly lower magnitude than those in winter. Further up-glacier (6–9  km), summer ice surface velocity also accelerated (~10 m yr⁻²). Conversely, at KS, where velocities decline toward the terminus (Fig. 2b), there have been minimal changes in seasonal velocity with just a few locations along the centreline exhibiting statistically significant trends. Specifically, in autumn at ~4 km from the terminus, velocities were decreasing at approximately 6.2 m yr⁻² (Fig. 6). The absence of clear changes in seasonal velocities at KS is consistent with the minimal change in annual velocity (Figs. 2d and 4a).

Figure 5. Time series of ice surface velocity at (a) 1  km and (b) 10 km from the terminus at the lake-terminating Isortuarsuup Sermia (blue) and the land-terminating Kangaasarsuup Sermia (orange, dashed). These velocities are computed from image-pairs separated by ≤ 30 days. Lines show the rolling 28 day median.

Figure 6. Seasonal velocity trends (2016–2021) along glacier centrelines at Isortuarsuup Sermia (blue) and Kangaasarsuup Sermia (orange) for winter (DJF, circles), summer (JJA, triangles) and autumn (SON, crosses). Seasonal trends derived from velocity fields where date_dt ≤ 30 days. Trends are linear fits, and only significant (p ≤ 0.05) trends are shown; error bars denote 95% confidence interval. Points at the same distance from the terminus have been offset from one another to aid readability.

3.2 Surface elevation change

There is a clear thinning signal at both IS and KS since 2011 (Figs. 3 and 4) with the rate and magnitude of thinning increasing towards both termini (Fig. 7). At ~15 km from the terminus, KS was thinning at a rate of 0.8 ± 0.3 m yr⁻¹, and IS at 0.3 ± 0.2 m yr⁻¹. These increased to their greatest width-averaged rates of thinning of 3.1 ± 0.7 m yr⁻¹ at the terminus of KS, and 2.1 ± 0.6 m yr⁻¹ 900 m from the terminus at IS, where the elevation is 550 m. However, at the equivalent altitude at KS, the thinning rate is ~1.1  m yr⁻¹ and accounting for differences in elevation, and the associated changes in lapse rate and thus surface melt processes, rates of thinning at IS are typically between 0.33–0.65 m yr⁻¹ (interquartile range; median difference: 0.5 m yr⁻¹) greater than at KS (Fig. 7c). The net differences in elevation change (up to −21 m and −13 m, at KS and IS, respectively) are substantially more than the differences measured over stable terrain after co-registration, giving us confidence in these observations.

Figure 7. Rates of surface elevation change at (a) Isortuarsuup Sermia (September 2012–June 2021) and (b) Kangaasarsuup Sermia (June 2011– September 2021) from ArcticDEM. Rates determined from linear regression. Only significant trends (p $\leq 0.05$) are shown. Error bars denote 95% confidence interval. ArcticDEM was sampled every 100 m along 9(7) parallel lines spaced every 250 m across the glacier at IS(KS). Coloured lines show 500 m rolling width-averaged median. Purple line (right hand axis) denotes number of ArcticDEMs with valid elevation measurements at each point along centreline. Shading represents standard deviation of number of ArcticDEMs at each point along centreline to account for the parallel offsets. (c) Rate of surface elevation change from ArcticDEM, against surface elevation at Isortuarsuup Sermia (blue) and Kangaasarsuup Sermia (orange); coloured lines show rolling median over 50 m bins. In (a), (b) and (c) negative values indicate surface thinning. (d) Net change in surface slope (positive indicates surface steepening) between first and last DEM; line represents median change across parallel lines, shading denotes median absolute deviation.

Due to greater rates of thinning at lower elevations, changes in surface gradient are generally positive at both IS and KS. Whilst these increases are of low magnitude (Fig. 7d), they do indicate some surface steepening. At IS, the median change in gradient along the centreline was $0.03\pm 0.2^\circ$, and at KS was $0.08\pm 0.3^\circ$. Immediately proximate to the terminus at IS, there is evidence of a lessening of the surface gradient toward the terminus, whereas at KS surface slopes increased steadily over the lower ~3 km by approximately 0.5$^\circ$ (Fig. 7d). This change in gradient at IS corresponds with a clear decrease in rate of surface elevation change over the 1 km closest to the terminus (Fig. 7a)

3.3 Terminus position

Terminus retreat differed between the two outlets over the course of the study period (Fig. 8). Between August 2014 and September 2021 KS retreated 210 ± 46 m, at an average rate of 30 ± 4 m yr⁻¹, which is greater than the uncertainty in our method. By contrast, IS retreated more slowly at a rate of 9 ± 4 m yr⁻¹.

Figure 8. Relative terminus position at Isortuarsuup Sermia (blue) and KS (orange) July– September 2014–2022. Blue circles denote relative position measured using box method; dashed blue lines representative of winter advance. Orange line illustrates average annual relative terminus position, shading denotes 95% confidence interval.

A seasonal cycle is observed at IS with a median winter advance of 30 m (median absolute deviation (MAD): 18 m), and median summer retreat of 19 m (MAD: 12 m). Furthermore, the terminus at IS showed distinct across-glacier spatial variability with the most pronounced localised retreat (~200  m) occurring in 2019 when a more advanced section of the northern terminus retreated from the sublacustrine terminal moraine (orange arrow in Fig. 3c). Following retreat from this moraine, the glacier did not re-advance back on to it during the remainder of our observation period. By contrast, retreat at KS appears progressive and sustained.

3.4 Runoff

Modelled runoff showed liquid water discharge (2011–2021) at KS was greater than that at IS by a factor of ~3 (Fig. 9), which is consistent with its larger catchment. At IS, there is good agreement between the two climate models whereas at KS, cumulative annual runoff is 5–25% greater in RACMO. Annual peaks in average daily runoff were typically between 60 and 90 m³ s⁻¹ at IS, and 230–320 m³ s⁻¹ at KS. Linear regression of cumulative annual runoff, from both RACMO and MAR, against time showed no significant trend at either IS or KS (SI. 3). Over the study period there was no consistent change in the timing of runoff onset or cessation. Runoff typically started between the end of April and mid-May at KS, and a week later at IS, as expected given its higher elevation, and had generally ceased by early- and mid-October at IS and KS, respectively.

Figure 9. Cumulative (a & b) and average daily runoff rate (5-day rolling average) (c & d) at Isortuarsuup Sermia (a & c) and Kangaasarsuup Sermia (b & d). Colours denote regional climate model with MAR in turquoise, and RACMO in orange.

4. Discussion

The near-terminus increase in ice surface velocity at IS (Fig. 2a) is similar to the accelerations seen at many other lake- and marine-terminating outlet glaciers in recent years (e.g. Joughin and others, 2018; Baurley and others, 2020). Our findings suggest that the presence of water at the terminus of IS, and the associated effects on near-terminus force balance, has enabled the observed changes in ice dynamics. This is reflected in the shape of the surface velocity profile, which for any given year shows velocity increasing towards the terminus from ~10  km up-glacier (Fig. 2a), and the significant increase in velocity over time along the entire 15 km centreline. We consider below evidence for potential drivers of the dynamic changes observed at IS in contrast to the behaviour at the neighbouring land-terminating KS.

In agreement with previous work (e.g. Tedstone and others, 2015), these data show no statistically significant relationship between cumulative runoff and average annual velocity at either glacier (SI. 3), and no significant melt trend was observed, either with respect to melt volume or timing. This result does not preclude a relationship between runoff and ice velocity over shorter timescales than are frequently measured from remotely sensed observations of velocity and gridded runoff estimates derived from regional climate models. Indeed, clear seasonal velocity cycles are evident at both glaciers, reflecting the commonly observed coupling between seasonal runoff, the hydraulic evolution of the subglacial drainage system and ice-dynamics (Davison and others, 2019). Nevertheless, the absence of a relationship between annual runoff and ice velocity at the lake-terminating IS is important, and suggests that the observed acceleration is not directly attributable to enhanced sliding due to increased meltwater input to the bed.

The observed acceleration is therefore likely the result of a dynamic feedback (Benn and others, 2007) driven by sustained negative surface mass balance induced thinning (The IMBIE Team, 2020), which is seen at both IS and KS (Fig. 7). Thinning will have brought the terminus closer to flotation and enhanced rates of basal sliding (Pfeffer, 2007; Tsutaki and others, 2019). This suggestion is supported by both the glacier-wide acceleration (Fig. 3a) and by the acceleration in all seasons (Fig. 6) at IS. Furthermore, due to minimal changes in lake water level, basal water pressures proximate to the terminus will be held approximately constant while ice-overburden pressure will decrease due to surface thinning (Figs. 3 and 7), leading to a long-term decrease in effective pressure. While a sustained decrease in effective pressure promotes the observed increase in ice motion at IS (Fig. 2a), a pronounced acceleration at the terminus occurs between 2018 and 2019 (Fig. 2b). This results from ongoing thinning at the glacier terminus promoting flotation and the subsequent retreat of the northern part of the terminus away from the Little Ice Age sublacustrine moraine (Fig. 3c). The timing of this flotation and retreat is further evidenced by the substantial change in rate of surface elevation change post-2018 (SI. 4) in conjunction with the pronounced increase in velocity, presumably in response to the associated removal of buttressing. The change in surface elevation suggests a clear hinge (above the grounding line) ~1 km up-glacier from the terminus, akin to what has been observed at Helheim Glacier (James and others, 2014). Furthermore, acceleration following retreat from a sublacustrine moraine replicates the behaviour observed at the lake-terminating Yakutat Glacier, Alaska (Trüssel and others, 2013).

The acceleration and associated extensional flow has led to enhanced rates of surface lowering (2.1 ± 0.6 m yr⁻¹) across the whole width of the terminus region (Figs. 3b and 7c). These spatially variable rates of thinning change the glacier surface slope, which have generally steepened up-glacier (Fig. 7d), likely leading to an increase in driving stresses, which may promote further acceleration and thinning (Howat and others, 2005). This is consistent with modelling work that demonstrated near-terminus velocities at lake-terminating glaciers increased with surface slope (Pronk and others, 2021). The absence of bed topography estimates near the terminus inhibit efforts to quantify stress-coupling lengths (e.g. Enderlin and others, 2016); however, it is suggested here that a high degree of longitudinal coupling allows these thinning-driven near-terminus accelerations to propagate ~15  km up-glacier (Fig. 2e).

The observed rates of change in terminus position at these two outlets (Fig. 8) are consistent with those previously documented (Warren, 1991). Between 1943 and 1983 KS retreated at an average rate of 38 m yr⁻¹, whilst the terminus at IS was shown to be stable (1949–1985). Additionally, there has only been ~400 m of retreat at IS since the Little Ice Age during which time KS has retreated approximately 3.4 km (Weidick and others, 2012). This observation of greater retreat at the land-terminating KS, as opposed to the lake-terminating IS, differs from the regional average pattern along the entire south-west margin where recent rates of margin recession were typically greater along lacustrine sections (mean annual rates of margin change: −11.5 m yr⁻¹, 2010–2015) than land-terminating margins (−2.8 m yr⁻¹) (Mallalieu and others, 2021). Nevertheless, this anomalous result is not unexpected given that the range of annual rates of margin recession at lacustrine (n=374) and terrestrial margins (n=3325) are approximately equivalent (Mallalieu and others, 2021, Figure 4b), and we are only presenting results from two such points. Furthermore, sustained long term stability at IS is evidenced by the minimal retreat and ongoing proximity of the terminus to the Little Ice Age maximum (Weidick and others, 2012) (leftmost white arrow in Fig. 3c) and the large sublacustrine moraine. This suggests that the topographic configuration at IS has enabled this stable terminus position in a manner similar to the long-term stability observed at other lake- (Trüssel and others, 2013) and marine- (Catania and others, 218) terminating glaciers. For example, the bed topography at Store Glacier (e.g. Catania and others, 218; Box and Decker, 2011) has promoted stability at its terminus, contrary to the regional trend, and in spite of being in sustained negative balance. Additionally, with respect to terminus position at IS, the recent dramatic increase in ice surface velocity (Figs. 2a and 3a), may in part offset any retreat from frontal ablation.

At Breiðamerkurjökull, Iceland, between 1982 and 2018 there was terminus retreat and a corresponding increase in proglacial lake area at Jökulsárlón and Breiðárlón (Baurley and others, 2020). However, the net retreat at these two proximate lake-terminating margins differed by a factor of three. Furthermore, and in line with our observations, cumulative retreat at Breiðárlón over this period was less than at an adjacent land-terminating section. These authors suggest that the relative stability of the terminus at Breiðárlón is attributable to the shallow subglacial trough (60 m vs. 300 m at Jökulsárlón). At present, there are no available lake bathymetry estimates at Isortuarsuup Tasia, and ice thickness estimates from BedMachine v5 (Morlighem and others, 2022) are not available proximate to either terminus. Additionally, the long-term stability of the terminus at IS (Weidick and others, 2012) suggests that bed topography, lateral support from the valley sides and the development of the large terminal moraine, have all exerted a strong control on terminus position (Warren, 1991).

The differences in velocity between IS and KS support recent work demonstrating that ice-contact lakes amplify near-terminus velocities (Baurley and others, 2020; Pronk and others, 2021). In Greenland, this velocity uplift has been estimated to be ~25% (Carrivick and others, 2022). However, this value of 25% contrasts ice-motion adjacent to all styles of ice-marginal lake, including ice-dammed lakes along valley sides tangential to the main ice-flow, with that of all land margins, regardless of whether the margin is terminal or not. However, it is likely that this velocity difference is greater when considering only terminal ice-contact proglacial lakes (i.e. those where the main flow unit within an outlet glacier is flowing directly into the lake). While ice- and moraine-dammed lakes are susceptible to periodic draining and catastrophic outburst floods (Costa and Schuster, 1988; Carrivick and Tweed, 2013), bedrock-dammed lakes are inherently stable and maintain a greater level of hydraulic connectivity to the up-glacier system (Sugiyama and others, 2011; Carrivick and Tweed, 2013). This is illustrated by the ratio of ice surface velocities at IS and KS (2 km from their respective termini) differing by a factor ~1.5 in 2013 and ~3.4 in 2021. Furthermore, this ratio increases toward the terminus, with velocities at IS an order of magnitude greater than those at KS at the end of the study period. Additionally, bedrock overdeepenings are typically sited in regions where ice flow is laterally constrained by topography, and their association with outlet glacier confluences means that they are often within large ice catchments (Patton and others, 2016). Consequently, bedrock-dammed proglacial lakes occupying glacially eroded valley bottom overdeepenings, are likely to be of greater importance in controlling ice sheet mass balance than ice-dammed marginal lakes, and future investigations should prioritise these.

In summary, the near terminus thinning and acceleration observed at IS highlight the potential importance of proglacial lakes, as the combined effects on ice dynamics reach inland and can lead to greater rates of mass loss. The behaviour replicates the expected positive feedback effects associated with sustained terminus thinning at calving glaciers (Benn and others, 2007) and we argue that as both ice-marginal melt-rate (The IMBIE Team, 2020) and lake number (Carrivick and Quincey, 2014; Shugar and others, 2020) increase, so too will the significance of ice-marginal lake processes for GrIS mass loss. The importance of this behaviour can currently be seen in Alaska (Trüssel and others, 2013; Larsen and others, 2015), Novaya Zemlya (Carr and others, 2017), and the Patagonian ice fields where ice mass loss is strongly controlled by fast-flowing lake-terminating outlets (Sakakibara and Sugiyama, 2014).

5. Conclusion

As the margin of the Greenland Ice Sheet recedes, ice-marginal lakes are expected to increase in both number and area in the coming decades, with an attendant increase in their influence on the wider ice sheet (Carrivick and others, 2022). Our findings suggest that the distinct dynamic differences between the land- and lake-terminating outlets in this study are largely attributable to the presence of the lake. More specifically, we argue that the doubling of near-terminus ice velocity at Isortuarsuup Sermia is likely driven by ongoing negative surface mass balance and glacier thinning. This has reduced ice-overburden pressures near the terminus, where the lake maintains high basal water pressures year-round, and facilitated ice acceleration in all seasons. Furthermore, ongoing thinning and subsequent flotation off a sublacustrine moraine has instigated retreat, thereby promoting enhanced acceleration across the terminus region through the removal of buttressing. In contrast, reductions in ice thickness at the land-terminating KS have not led to flow acceleration due to the profound differences in terminus processes and subglacial hydrological setting. The acceleration and attendant extensional flow at Isortuarsuup Sermia has also led to enhanced rates of thinning near-terminus of between 0.33 and 0.65 m yr⁻¹.

Our observations show the effect of recent mass balance change on the ice-dynamics of a lake-terminating glacier reaches ~15 km up-glacier, highlighting the ability of proglacial lakes to perturb inland ice. This supports earlier observations (Kirkbride, 1993; Warren and Kirkbride, 2003; Sakakibara and Sugiyama, 2014; Tsutaki and others, 2019; Mallalieu and others, 2021) and modelling work (Sutherland and others, 2020) that stress the importance of proglacial lakes on glacier and ice sheet mass loss. We suggest future work should discriminate between dam type and lake setting (as per Rick and others, 2022) when evaluating ice-marginal lake impacts on ice dynamics, as we contend that the relative importance of proglacial bedrock-dammed lakes on ice sheet mass loss is likely greater than ice-dammed marginal lakes. Additionally, there is a need to establish whether the recent pattern of behaviour seen at Isortuarsuup Sermia is typical for other Greenlandic lake-terminating outlets. Accurately quantifying the effect of ice-marginal lakes on these glaciers demands greater knowledge of ice-marginal lake characteristics, including bathymetry. This work is timely, as climate warming is seeing the ice margin retreat towards the many glacially eroded overdeepenings beneath the Greenland Ice Sheet. An increased incidence of lake-terminating glaciers would likely enhance the dynamic mass loss from Greenland due to accelerated glacier flow, in line with expected positive feedbacks associated with melt induced thinning of these glacier termini (Benn and others, 2007), and as witnessed already across numerous glaciated regions including Alaska (Trüssel and others, 2013), Iceland (Baurley and others, 2020), and Patagonia (Sugiyama and others, 2011, 2019).