Investigating the transition in ice dynamics from a land- to lake-terminating glacier using a simple glacier dynamics model

Glaciers interact with adjacent proglacial lakes through a range of thermomechanical processes. These interactions occur in addition to climate-driven ablation and are capable of amplifying or modifying climatic effects through various feedback mechanisms. In particular, the connection between lake water-level and the subglacial hydrological system can reduce basal friction, which in turn leads to increased glacier flow and dynamic thinning. This creates a positive feedback loop in which decreased effective pressure, also driven potentially by negative surface mass balance, enhances flow velocity, in line with similar processes observed at marine-terminating glaciers.

Our aim in this study is to introduce a simple model that can be used to understand the critical controls on the dynamic behaviour of glaciers as they transition from land- to lake-terminating systems. Here, we investigate the behaviour of Skaftafellsjökull in Iceland, which has undergone such a transition over the past twenty-five years. More specifically, we use the Shallow Shelf Approximation (SSA) in Elmer/Ice to model ice velocities, incorporating a water pressure-dependent friction law to model basal sliding with a simple parameterization of basal water pressure.

The model successfully reproduces the observed velocity patterns, capturing the shift from downstream deceleration near the glacier terminus in 2010 to acceleration in 2018, reflecting the growing influence of the proglacial lake. We find a threshold behaviour between basal water pressure and ice velocity, whereby small increases in water pressure beyond a critical value led to strong acceleration, consistent with previous empirical observations. Furthermore, our results imply that surface thinning exerts a stronger control on the near terminus acceleration than the observed terminus retreat. Our results suggest that the modelling framework developed provides a valuable tool for simulating these complex interactions in a computationally efficient manner.

1 Introduction

Lake-terminating glaciers interact with adjacent proglacial lakes through a series of thermomechanical processes (Carrivick and Tweed, 2013; Truffer and Motyka, 2016; Tsutaki et al., 2011). These processes operate in addition to climate-driven surface ablation, but these systems may also interact via multiple feedback mechanisms. Calving events can temporarily destabilize the ice margin, partially decoupling glacier dynamics from climatic forcing (Kirkbride, 1993). Furthermore, since there is a hydraulic connection between the proglacial lake and the subglacial drainage system, basal friction is reduced, facilitating higher glacier flow speeds, through a reduction in basal resistance (Carrivick et al., 2020). The resulting acceleration promotes dynamic thinning (Trüssel et al., 2013; Tsutaki et al., 2019), initiating positive feedback in which thinning decreases effective pressure which in turn enhances basal motion and further increases flow velocity (Sugiyama et al., 2011). Increased basal water pressure not only promotes sliding, but can also result in basal crevasse development and propagation, favouring calving (cf. Van Der Veen, 1998). In addition, thermally-induced melting at the lake-ice interface can cause notches to develop at the water line, and this thermal undercutting can favour calving (Röhl, 2006).

Many examples of contrasting dynamics between land- and lake- terminating glaciers, with associated impacts on mass balance, can be found across the glaciated regions on Earth. For example, more negative mass balance at lake- compared with land-terminating glaciers is observed in the Himalaya (King et al., 2019) where centre-line flow velocities of lake-terminating glaciers are on average more than double those of land-terminating glaciers in the central and eastern Himalaya (Pronk et al., 2021). In Greenland, annual surface velocity at the terminus of the lake-terminating glacier Isortuarsuup Sermia increased by a factor of 2.5 during the period 2013–2021 while the near-terminus surface elevation thinned at a rate of −2.3 ± 1.1 m yr⁻¹, whereas velocity change at the nearby land-terminating Kangaasarsuup Sermia was minimal, and surface thinning rates were approximately half at comparable elevations (−1.2 ± 0.3 m yr⁻¹) (Holt et al., 2024).

Furthermore, in Iceland, pronounced differences in surface velocity, front retreat and surface elevation were found between the lake-terminating eastern arm and land-terminating central arm of Breiðamerkurjökull, a glacier draining Vatnajökull Ice Cap. Baurley et al. (2020) found the eastern arm's surface velocities increased from 1.0 ± 0.36 m d⁻¹ in 1991 to 3.5 ± 0.25 m d⁻¹ in 2015, while the central arm showed no discernible change in velocity over the same period. In addition, the eastern arm retreated ∼ 3.5 km between 1982–2018, whereas the central arm retreated ∼ 2 km over the same period (Baurley et al., 2020).

Icelandic glaciers have been experiencing mass loss since the end of the Little Ice Age (∼ 1890). However, in recent decades, rates of mass loss have accelerated, with nearly half of the total mass loss since 1890 having occurred since 1994 (Aðalgeirsdóttir et al., 2020). While glacier mass loss has generally increased in recent decades, there is substantial interannual variability (Foresta et al., 2016). For instance, due to regional cooling in the North Atlantic, the average rate of mass loss since 2011 has been about 50 % lower than between 1995–2010 (Noël et al., 2022a), yet 2018–2019 marked one of the most extreme years of loss on record (−15 ± 1.6 Gt yr⁻¹) (Aðalgeirsdóttir et al., 2020). Recent studies reveal that processes beyond surface melting, such as geothermal heating, volcanic activity and frontal ablation (via glacier calving), have played a significant role in shaping these recent patterns of glacier mass change (Jóhannesson et al., 2020). These non-surface driven processes are estimated to account for roughly 20 % of total glacier mass loss since 1994, and at certain ice caps, particularly in southern Vatnajökull, they may contribute nearly 40 % (Aðalgeirsdóttir et al., 2020).

Among these processes, frontal ablation at lake-terminating glaciers can be especially important, causing glaciers to lose mass rapidly and retreat independently of direct climatic influences (Benn et al., 2007b; Carrivick and Tweed, 2013). Although calving played only a minor role in glacier change during the early 20th century, its importance has risen since the mid-1990s (Aðalgeirsdóttir et al., 2020). This is due primarily to outlet glaciers retreating into overdeepened bedrock troughs, thereby promoting the formation and growth of proglacial lakes (Guðmundsson et al., 2019).

Notable examples include Svínafellsjökull, Breiðamerkurjökull, Hoffellsjökull and Skaftafellsjökull, all southerly-flowing outlets of Vatnajökull with bedrock troughs reaching between 150 and 200 m below sea level (Guðmundsson et al., 2019; Magnússon et al., 2012). These deep troughs facilitate rapid proglacial lake expansion, enhancing calving activity and thus accelerating terminus retreat and mass loss (Dell et al., 2019; Hannesdóttir et al., 2015). Such trends are expected to persist into the future, not only in Iceland but globally given trends in increasing proglacial lake extent (Zhang et al., 2024b), with major implications for glacier dynamics and regional hydrology (Schomacker, 2010). As a result, there is a pressing need to improve our ability to model the current and future dynamics of glaciers terminating in these evolving proglacial lake systems.

There have been several recent efforts to model lake-terminating glaciers and to compare them with land-terminating glacier systems. Sutherland et al. (2020) used the vertically integrated ice flow model, BISICLES, to analyse the effects of a potential proglacial lake on the retreat of the Pukaki Glacier, New Zealand, since the Last Glacial Maximum. The proglacial lake produced a maximum effect of a four-fold increase in grounding line recession and an eight-fold increase in ice velocities, relative to simulations of a land‐terminating glacier forced by the same climate. Elsewhere, a depth integrated synthetic flow line model, that implicitly employed the Shallow Shelf Approximation (SSA), and assumed that subglacial water pressure follows a piezometric surface equivalent to 60 % of the ice thickness away from the ice front, suggested that changes in force balance at the glacier termini are likely to play a key role in accelerating near-terminus glacier flow, with variations in basal friction being of only minimal importance (Pronk et al., 2021).

The modelling approaches used for marine-terminating glaciers also have considerable relevance and potential for investigating the dynamics of lake-terminating glaciers given the numerous similar characteristics of both systems (Benn et al., 2007b). However, certain mechanisms that concern marine-terminating glaciers are not relevant on lake-terminating glaciers; for example, submarine melting produced by buoyant freshwater plumes from subglacial discharge (Jenkins, 2011), are absent at lake-terminating glaciers due to the lack in density contrast between the subglacial discharge and freshwater proglacial lake.

In this study, our aim is to investigate the changes in dynamics that occur during the transition from a land- to lake-terminating glacier and the mechanisms responsible for such changes, using a simple and computationally efficient model. Initial experiments were undertaken to determine the relationship between basal water pressure, effective pressure and ice velocities. The work then sought to determine the contrasting roles of surface thinning and ice front retreat on the force balance and thus ice motion. Finally, the work aimed to determine the overall rate of frontal (lake-terminus) ablation and the potential impact of lake-influenced ice-dynamics in amplifying glacier mass loss.

2 Data

2.1 Study site

Skaftafellsjökull is a small valley outlet glacier located on the western side of Öræfajökull Ice Cap, which itself is situated on the southern margin of Vatnajökull Ice Cap, SE Iceland (Fig. 1a). Skaftafellsjökull, with an area of 84.1 km² in 2010 (Baurley et al., 2025), descends from 1500 m elevation on the north-western margin of Öræfajökull to ∼ 100 m above sea level (Fig. 1b) before terminating in an unnamed proglacial lake which started forming around 2001. The lake had expanded to ∼ 1.3 km² by 2018, with an average depth of 23 m and surface elevation at 87 m a.s.l., and developed initially as multiple small lakes that grew rapidly and merged into a single lake in 2011 (Guðmundsson et al., 2019). The lake has an outlet river, which drains from the centre of the topographic dam and keeps the water level stable (Fig. 3a). The glacier has a substantial high-altitude accumulation area which drains to a narrow ∼ 2 km wide tongue, via an ice fall flowing between several nunataks (Fig. 1b). The lower glacier tongue flows across an overdeepened bed extending up to 200 m b.s.l. (Fig. 2).

Figure 1 (a) Location of Skaftafellsjökull (in red) and Vatnajökull Ice Cap with glacier outlines from the Randolph Glacier Inventory (RGI). Inset shows the location of Vatnajökull Ice Cap in Iceland. (b) 2020 Surface elevation map from ArcticDEM. Background images from Google Maps (Imagery © 2026 NASA, Map data © 2026 Google).

Two key features make Skaftafellsjökull an ideal candidate for developing a model seeking to understand the critical processes operating when a glacier transitions from land- to lake-terminating. Firstly, the glacier is in the early stages of such a transition with ongoing retreat since the initial appearance of the proglacial lake in 2000, and secondly, there is a substantial volume of high resolution spatial and temporal data available, including bed geometry, surface velocity, surface elevation and terminus position during the transition period.

Figure 2 Bed topography of Skaftafellsjökull derived from RES surveys (Magnússon et al., 2012, 2014) showing a substantial overdeepening extending to below sea-level (dark blue zone) in the central part of the glacier tongue.

2.2 Observational Data

We present glaciological observations from recent decades at Skaftafellsjökull which was used in our modelling study and gives an overview of recent changes. Glacier bed elevation, surface elevation, and front position provide ice model geometry in our experiments while ice velocities are used to tune and validate our model parameters. Surface mass balance is not used directly in our model, as we do not evolve the model in time, but it is used in order to characterise the glacier's frontal ablation (Sects. 3.5, 4.4).

2.2.1 Glacier terminus positions

Ice-front positions were derived by manually delineating the terminus of the glacier for each year between 1985 and 2024 using the Google Earth Engine Digitisation Tool (GEEDiT: Lea, 2018). The delineations were performed on ASTER, Landsat 5–8, and Sentinel-2 satellite imagery. To ensure consistency across the time series, termini were preferentially digitized from images acquired near the end of the melt season (i.e., August–September), capturing the glacier's minimum annual extent. Width-averaged terminus position change was subsequently quantified using the rectilinear box method (Moon and Joughin, 2008).

Skaftafellsjökull's terminus was approximately stable between 1985 and 1999 with a mean retreat rate of 3.5 m yr⁻¹. Glacier retreat then accelerated, especially across the eastern margin, at a mean rate of 50 m yr⁻¹ from 2000 to 2024. Over the whole period from 1985 to 2024, the front retreated on average more than 1.2 km of which ∼ 1.1 km was after 2000, following retreat of the terminus into the proglacial lake (Fig. 3).

Figure 3 (a) Front positions of Skaftafellsjökull between 1985 and 2024 (white lines) derived from ASTER, Landsat 5–8, and Sentinel-2 satellite imagery with positions in 2000 (orange), 2010 (pink) and 2020 (red) highlighted. Background image from Google Satellite (Imagery © 2026 Airbus, Maxar technologies, CNES/Airbus, Map data © 2026). (b) Width-averaged front retreat from 1985.

2.2.2 Glacier surface elevation

The surface DEM for 2010 was obtained from the surface DEM of Öræfajökull Ice Cap derived from airborne LiDAR surveys (Jóhannesson et al., 2013). The 2018 and 2020 ice surface elevations were derived from ArcticDEM strips, using high-resolution (0.32 to 0.5 m) Maxar satellite imagery with a spatial resolution of 2 m (Porter et al., 2022). These DEMs were co-registered to the stable, ice-free terrain of the ArcticDEM Mosaic (Porter et al., 2022). The DEMs were corrected for horizontal and vertical biases using the method proposed by Nuth and Kääb (2011), implemented by the Python package xDEM (xDEM contributors, 2024).

Between 2010 and 2018, the glacier showed a slight thickening in the accumulation area (mean value of 1.1 ± 0.45 m yr⁻¹) and a thinning across the tongue (mean value of −1.9 ± 1.5 m yr⁻¹), which is more rapid near the glacier terminus (more than −7 m yr⁻¹) (Fig. 4c).

Figure 4 (a) 2010 surface elevation DEM from airborne LiDAR surveys. (b) 2018 surface elevation DEM from ArticDEM. (c) Rate of surface elevation change (2010–2018) showing thinning near the glacier front.

2.2.3 Bed elevation

The bed elevation at Skaftafellsjökull is obtained from the bedrock DEM of Öræfajökull Ice Cap, derived from radio echo sounding (RES) profiles at higher elevations on the ice cap (> 1250 m) and RES point measurements in the lower tongue (below ∼ 900 m elevation), and carried out between 1991 and 2012 (Magnússon et al., 2012, 2014).

The bed data derived from the RES point measurements have a higher error than areas surveyed with RES profiles because first, migration of this type of dataset is not possible due to its discrete nature, and second, the data discontinuity makes misinterpretation of the bedrock reflection more likely than for the continuous RES measurements (Magnússon et al., 2012). For these reasons, the error in determining the bed at the tongue of the glacier could be substantial (and is estimated by Magnússon et al., 2012 as ± 20 m).

The most important feature of the glacier's basal geometry, in terms of potential proglacial lake development, is a substantial overdeepening present in the central part of the tongue. This depression reaches up to 200 m b.s.l., with a maximum depth 5.2 km upglacier from the 2000 terminus, and extends from immediately downglacier of the ice fall to beyond the current glacier terminus, where the proglacial lake is now present. The elevation of the barrier impounding the lake to the south is 87 m a.s.l. (Figs. 2 and 3).

2.2.4 Surface velocities

Surface velocities were generated using the offset tracking algorithm within the European Space Agency (ESA) Sentinel Application Platform (SNAP) by Baurley et al. (2025). These velocity data were generated from imagery acquired in late summer (August–September); hence they only represent late-summer velocities, not annual averages. This study is focused on the role of lake-ice interactions on near-terminus ice-dynamics, so only surface velocities across the glacier tongue downstream of the icefall are used to compare observed and modelled velocities.

Velocities decrease along the 5 km glacier tongue from the base of the ice fall to the terminus in 2010, reaching ∼ 18 m yr⁻¹ near the front (Figs. 5a and 7a). However, in 2018, while velocities decrease initially downglacier from the ice fall, they remained higher than in 2010 and in the terminus region, velocities actually started to increase ∼ 1 km upglacier from the terminus, reaching ∼ 78 m yr⁻¹ (Figs. 5b, 7b and 8), near the glacier front (Baurley et al., 2025).

Figure 5 Maps of surface velocity derived from Sentinel (Baurley et al., 2025). (a) 2010 surface velocity. Black line is the central flowline used in Fig. 8. (b) 2018 surface velocity. (c) Change in velocity between 2010 and 2018.

The difference in velocity between 2010 and 2018 suggests a “spotty” patterning (Fig. 5c) in the lower glacier tongue; this pattern might be realistic but given its structure is perhaps more likely to be a processing artefact. As a result, we use averaged swaths of velocity, at a given distance upglacier, in our analyses rather than using specific at-a-point observations. Importantly, both the differenced velocity fields (Fig. 5c) and averaged swaths confirm the pronounced increase in velocity, and thus velocity difference between 2010 and 2018, as one approaches the glacier terminus.

2.2.5 Surface mass balance

Surface mass balance (SMB) is derived from the daily output of RACMO2.3, forced by ERA reanalyses (1958–2019) statistically downscaled to a 500 m ice mask derived from the Randolph Glacier Inventory (RGI Consortium, 2017) version 6 and the 100 m spatial resolution ArcticDEM (Porter et al., 2018), re-sampled onto a 500 m grid by Noël et al. (2022a).

In the 2010–2019 period, estimated SMB for Skaftafellsjökull was positive every year except 2010 when it was slightly negative (Fig. 6), concurrent with the massive eruption of Eyjafjallajökull that covered surrounding ice caps with dark volcanic tephra, reducing surface albedo (Gudmundsson et al., 2012). The 2011–2019 period represents an anomalous period of surface mass balance where the glaciers across Iceland experienced a slowdown in their rates of mass loss that coincided with a regional cooling across the North Atlantic Ocean (Noël et al., 2022a). Furthermore, while Icelandic glacier mass loss continued but at a slower rate during this period, Skaftafellsjökull's positive SMB reflected its unusual hypsometry with gains in the large high altitude accumulation area dominating over mass loss in the narrow ablation zone (Fig. 4c).

Figure 6 SMB of Skaftafellsjökull during the period 2010–2019 derived from RACMO.

3 Methods

3.1 The ice-flow model

Ice models simulate the dynamic evolution of ice across multiple scales, solving ice flow equations through different approximations. One of the simplest approaches is the shallow shelf approximation (SSA; Weis et al., 1999). The SSA is suited to floating ice shelves and low-drag ice streams dominated by longitudinal (or membrane) stresses, where the driving stress is low and the ice column has an almost uniform velocity profile, i.e., plug flow (Robinson et al., 2022). The SSA model assumes a vertically constant velocity profile, such that it only accounts for the longitudinal components of the stress while neglecting the vertical components; it is then depth-integrated and thus two-dimensional (Jouvet, 2014).

As elsewhere, we infer that the glacier motion in the near-terminus region of lake-terminating Skaftafellsjökull is dominated by rapid flow due to low basal resistance, such that the deformation velocity is very small compared to sliding velocity. As a result, the approximation of generating velocities using the SSA, which only takes account of sliding, may be sufficiently accurate in the most important part of the glacier for this study, i.e., the tongue region. In this paper, all experiments are conducted using the SSA implementations of Elmer/Ice (Gagliardini et al., 2013), an open-source finite-element ice dynamic model. An unstructured finite-element mesh with triangular elements typically of 150 m edge length, but denser near the front (as small as 50 m edge length), is created using Gmsh (Geuzaine and Remacle, 2009) for the footprint of the glacier.

Basal water pressure plays a key role in controlling ice motion (Iken and Bindschadler, 1986) including at lake-terminating glaciers (Sugiyama et al., 2011). Here, basal sliding is parameterized using a regularized Coulomb-type friction law (Gagliardini et al., 2007; Schoof, 2005), which is a non-linear water pressure-dependent friction law, implemented in Elmer/Ice. The dependence from the water pressure is given through the effective pressure at the ice-bed interface N, defined as the difference between ice-overburden and basal-water pressure.

Based on Gagliardini et al. (2007), the friction law can be written as

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{b}}=N{\left[{\displaystyle \frac{u_{\mathrm{b}}}{u_{\mathrm{b}}+N^n{A}_{\mathrm{s}}}}\right]}^{\frac{\mathrm{1}}{n}}\end{array}$$

where τ_b is the basal drag, u_b is the basal velocity, n is the flow-law exponent (n=3) and A_s is the sliding coefficient.

The effective pressure N can be expressed as

$$\begin{array}{}\text{(2)}& N={\mathit{\rho }}_{\mathrm{i}}gH-{\mathit{\rho }}_{\mathrm{w}}g{P}_{\mathrm{wp}}\left(z_{\mathrm{l}}-z_{\mathrm{b}}\right)\end{array}$$

where ρ_i is the ice density, g is the acceleration due to gravity, H is the ice thickness, ρ_w is the water density, z_l is the lake surface elevation, z_b is the glacier bed elevation, both above sea level, and P_wp is a scaling parameter to take in to account the increasing water pressure above its “minimum level”. This minimum level of water pressure is the pressure which corresponds with the water height equal to the elevation of the lake surface, if we assume a “perfect” hydrological connection between the lake and the upglacier subglacial hydrological system. This idealised assumption is based on the consideration that the hydrostatic pressure exerted by the lake water against the calving face prevents the subglacial water level from dropping below the lake's surface level. This assumption is commonly used when modelling both marine (e.g., Bondzio et al., 2016; Cook et al., 2014) and lake-terminating glaciers (e.g., Carrivick et al., 2020; Sugiyama et al., 2011). Furthermore, the limited evidence available from observations of subglacial water pressures, derived from borehole measurements at a lake-terminating glacier, suggest that the hydraulic head was in fact higher than the surface level of the proglacial water body (Sugiyama et al., 2011).

P_wp, the water pressure parameter, is a free parameter of the model that was added to parameterize the influence of water pressure on sliding velocity through the friction law. Such a parameterization is a simplified way of representing the basal hydrology. Another way of incorporating water pressure in models would be to consider water pressure as a fraction of ice overburden pressure for bedrock below sea level, as used by Martin et al. (2011) in an Antarctic Ice Sheet model, or in a model of lake terminating glaciers in the Himalaya, in which subglacial water pressure is assumed to follow a piezometric surface equivalent to a fraction of the ice thickness away from the calving front (Pronk et al., 2021).

We additionally note that, in using P_wp, we are not suggesting that the complexity of subglacial hydrology can be represented by a single scalar parameter. However, a more physical model, such as one that represents evolving subglacial cavities and channels (Gagliardini and Werder, 2018), would be computationally expensive and involve a large number of unknown parameters, which would be very difficult to constrain due to the lack of direct observations. The use of a single scalar parameter which represents broad-scale effects on basal pressure, and which can be estimated through simple tuning, allows more extensive investigation of the effect of lakes on near-terminus velocities.

Near the ice front, the vertical force due to the static weight of the ice exceeds the horizontal static force due to the water (Cuffey and Paterson, 2010). Therefore, a depth-integrated Neumann boundary condition was applied, following Cuffey and Paterson (2010), to model the margin stress as

$$\begin{array}{}\text{(3)}& \mathit{\varphi }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{i}}g{H}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{w}}g{\left(z_{\mathrm{l}}-z_{\mathrm{b}}\right)}^{\mathrm{2}}\end{array}$$

The free parameters, water pressure parameter P_wp, sliding parameter A_s and viscosity parameter A were tuned running simulations with both the 2010 and 2018 geometries. These times correspond with a period which overlaps with the land- to lake-terminating transition, but also years when both surface velocity and surface elevation datasets were available. For each year's geometry, 133 000 simulations (190 A_s values × 100 A values × 7 P_wp values) were run (see Table 1a), with the aim to minimize the mismatch between the observed and modelled surface velocities in both 2010 and 2018 and also, to minimize the difference between the observed and modelled acceleration, in a region of the tongue that is 1 km wide (centred) and ∼ 5 km long extending upglacier from the terminus; this is the so called “tuning” region. The key aim of this tuning process is to obtain a combination of parameters where the model is able to reproduce the observed behaviour of glacier surface velocities, during the transition from land- to lake-terminating conditions.

In order to derive optimal parameters, a number of metrics were defined and applied to the experimental ensemble, quantifying fit to various aspects of the data. The types of metrics are listed below according to the corresponding observations; note that for each type, the Mean Absolute Error (MAE) between model and observations, as well as the Mean Squared Error (MSE), is applied, unless otherwise indicated. Results and optimising parameter pairs are shown in Table 1b.

• 2010 and 2018: the observed velocities in each year.

• Acceleration: the change in velocity between 2010 and 2018.

• Mean MAE: the algebraic mean of the 2010, 2018, and Acceleration MAEs.

• Weighted Mean: equivalent to $\frac{\mathrm{2}}{\mathrm{3}}{\mathrm{MAE}}_{\text{Acceleration}}+\frac{\mathrm{1}}{\mathrm{6}}{\mathrm{MAE}}_{\mathrm{2010}}+\frac{\mathrm{1}}{\mathrm{6}}{\mathrm{MAE}}_{\mathrm{2018}}$.

• Dist Zero: Root mean square of the 2010, 2018 and Acceleration MAEs.

• Mean MSE: the algebraic mean of the 2010, 2018, and Acceleration MSEs.

Table 1 (a) Minimum and maximum values of the free parameters used in the tuning process, number of values taken for each parameter (N) and best-fit set of parameters selected after the process. For each parameter, values are evenly spaced between the maximum and minimum. (b) Set of free parameters found to minimize the various metrics. Bold value of each column corresponds to the minimized metric.

Results show that similar sets of parameters consistently minimized the metrics tested. There is little variation in each parameter relative to the range considered, especially for the metrics that take multiple observations into account (first three rows of Table 1b). Therefore, the set of parameters A_s=101, P_wp=1.5 and A=20 was chosen as the “best-fit” parameters, and used subsequently for the rest of the simulations discussed in the results below.

3.2 2030 projection

In a final experimental projection, the geometry, surface elevation and front position of the glacier was extrapolated to 2030 using the data and trends observed between 2010 and 2018 or 2020, in order to predict how the glacier's dynamics might evolve over the coming decade (Fig. 9a). The surface DEM for 2030 is estimated assuming a linear extrapolation of elevation change from the 2010 and 2020 DEMs. To estimate the 2030 front position, the 2010–2020 width averaged retreat rate (Fig. 3b) is used to derive the front position for 2030. The lateral boundaries were kept fixed in the projection, and the depth of the lake was estimated from the basal geometry (Fig. 2). The model with the “best-fit” set of parameters was thus applied to this extrapolated geometry to estimate glacier velocities in 2030.

3.3 Basal water pressure

Sensitivity analyses were undertaken to evaluate the behaviour of the modelled ice velocities to different values of basal water pressure, and thus different values of effective pressure. These experiments were run for both 2010 and 2018 geometries, with the sliding and the viscosity parameters fixed, while varying the water pressure parameter P_wp, to account for different basal water pressures. We used 7 values of P_wp, ranging from 1.1 to 1.7 for each date, and in order to compare this parameter with a representative velocity scale, averages of the values of velocity and water pressure from all of the nodes within 1.5 km of the terminus were used.

3.4 Influence of front retreat and thinning on surface velocities

It is recognised that frontal ablation at tidewater glaciers can lead to acceleration due to loss of resistive stresses (i.e. backstress) in the glacier (Howat et al., 2005; Thomas, 2004). As the terminus of Skaftafellsjökull retreated from 2010 to 2020, and ongoing retreat is expected, we must consider the possibility that any changes in the observed velocity could be attributed to loss of backstress resulting from calving front retreat. In addition, sustained terminus thinning, driven by both ongoing negative surface mass balance and ice acceleration (and thus dynamic thinning) will reduce near-terminus effective pressure thereby facilitating acceleration via basal sliding (Benn et al., 2007a; Holt et al., 2024; Pfeffer, 2007).

Experiments were therefore undertaken to test the competing influences of front terminus retreat and ice-thinning on surface velocities. These experiments control for either terminus retreat or ice thinning by combining front position change and thickness change from different epochs using the 2010 and 2020 observed geometries and 2030 extrapolated geometry. Five simulations were run: (1) with 2010 geometry (S10F10); (2) with 2010 surface elevation and 2020 front position (S10F20); (3) with 2020 geometry (S20F20); (4) with 2020 surface elevation and 2030 front position (S20F30); and (5) with 2030 geometry (S30F30) (see Table 2a). These velocity calculations allow for a number of comparisons: for example, S20F20 vs S10F20 examines the effect of the 2010–2020 thinning only. Four such comparisons are defined in Table 2b where the comparisons labelled EXP1 isolate velocity changes driven by surface thinning, while those labelled EXP2 isolate velocity changes driven by front retreat.

Table 2 (a) Simulations conducted with different geometries (surface elevation and front position) to facilitate comparative analyses between the first and second experiments. SXXFYY represents a configuration with geometry from 20XX and front position corresponding to 20YY. (b) Outline of the first, EXP1, and second, EXP2, experiments for each time interval.

3.5 Frontal ablation

Quantifying frontal ablation is of critical importance for water-terminating glaciers, either marine- or lake-terminating, in order to determine the respective proportions of mass lost by either surface melt processes or by frontal ablation, the latter of which comprises both calving and subaqueous frontal melt. For example, across Greenland, recent estimates suggest that mass loss between 1992 and 2018 has come equally from both surface melt and frontal ablation at tidewater glaciers (Shepherd et al., 2020) while at Patagonian icefields, frontal ablation from 2000 to 2019 over a set of 30 lake-terminating and 15 ocean-terminating glaciers represented ∼ 34 ± 6 % of the total ablation, while for the subset of lake-terminating glaciers, frontal ablation accounted for ∼ 25 % of total ablation (Minowa et al., 2021).

Following McNabb et al. (2015), the total rate (or volume) of frontal ablation FA can be defined as the sum of the integral over the front surface of the ice velocity at the terminus u_t and the change of terminus position, FR (negative if the front is retreating). Given that the front surface is constantly changing, the integral of u_t, i.e., flux through the front, is estimated as a finite sum over a cross section (flux gate) of the glacier, located some distance upstream of the front (Fig. 12b). The flux gate is divided into vertical sections with flux through each section calculated as section surface velocity u_si, multiplied by a factor f to take account of decreasing velocity with depth, multiplied by section area S_i, and by cosine of the angle α_i between the velocity vector and the vector normal to the flux gate surface. Finally, the flux through the flux gate (F_fg) is the sum of the flux through each section

$$\begin{array}{}\text{(4)}& F_{\mathrm{fg}}=\sum _i{\mathrm{fu}}_{\mathrm{si}}{S}_i\cos {\mathit{\alpha }}_i\end{array}$$

As the flux gate is located some distance upstream of the front, mass changes due to surface mass balance ($\dot{b}$) after the ice passes through the fluxgate must be taken into account. The change of terminus position, FR, is calculated following the “box method” (McNabb and Hock, 2014; Moon and Joughin, 2008). Therefore, frontal ablation is defined as

$$\begin{array}{}\text{(5)}& \mathrm{FA}=F_{\mathrm{fg}}+\mathrm{FR}-\dot{b}\end{array}$$

4 Results and discussion

4.1 Surface velocity

Figure 7 Detail of surface velocity maps of the glacier tongue. Velocity contours are plotted at 10 m intervals. (a) Observed (left) and modelled (right) velocities for 2010 and (b) observed (left) and modelled (right) velocities for 2018.

Using the “best-fit” set of parameters found in the tuning process, the model reproduced surface velocities well for both 2010 and 2018. The observed and modelled 2010 velocities decrease along the 5 km glacier tongue from the ice fall to the terminus, reaching minimum values of ∼ 20 (observed) and ∼ 30 m yr⁻¹ (modelled) at the terminus (Figs. 7a and 8). By 2010, the glacier already terminated in the lake (Fig. 3), which is likely why the terminus velocity is nonzero, but it was shallow (∼ 12 m maximum deep) and discontinuous across the ice front (Guðmundsson et al., 2019). As such, it seems clear that the lake was influencing glacier dynamics, but not significantly at this time. By 2018, the observed and modelled velocities also decrease initially downglacier from the ice fall, but ca. 1.5 km from the terminus, both observed and modelled velocities begin to increase slightly towards the front, reaching ∼ 78 (observed) and ∼ 77 m yr⁻¹ (modelled). The observed increase in near-terminus velocity in 2018, in conjunction with the pronounced terminus acceleration between 2010 and 2018, is consistent with the behaviour typical of a glacier terminating in a deepening proglacial lake (e.g., Baurley et al.,2020; Wu et al., 2025) (Figs. 7b and 8). By 2018, the proglacial lake extended across the full glacier terminus width, covering 1.3 km² and reaching ∼ 50 m deep at its deepest point (Guðmundsson et al., 2019).

As noted in Sect. 3.2, we applied the model to the 2030 extrapolated geometry in order to investigate how the glacier's dynamics might evolve over the coming decade. The results predict a sustained terminus acceleration of ∼ 12 m yr⁻²; an approximately three-fold higher rate than was observed between 2010–2018, resulting in near terminus velocities of up to ∼ 265 m yr⁻¹ by 2030 (Fig. 9b). The increase in velocities shown by the model near the front in 2030 highlights the behaviour that this glacier may exhibit in the near future with potentially devastating consequences for glacier mass balance.

Figure 8 Comparison between observed (solid lines) and modelled (dotted lines) velocities for 2010 (red), 2018 (blue) and the 2010–2018 acceleration (orange) in a flowline along the centre of the glacier tongue (Fig. 5a). Velocities at each node of the flowline, ∼ 200 m apart, were calculated as an average of the velocities of nine points located in a 400 × 400 m grid centred on the node in order to avoid local effects.

Figure 9 (a) The extrapolated glacier geometry (surface elevation and front position) in 2030 (b) The modelled surface velocities of the glacier tongue based on the extrapolated 2030 geometry with contour lines plotted at 20 m yr⁻¹ intervals. Velocity scale is different from that of Fig. 7 to reflect the higher velocities.

4.2 Basal water pressure

Our tuning experiments were undertaken to find an optimal set of model parameters, but they also provide a means of analysing sensitivities to the basal pressure parameter, or alternatively to basal pressure itself. Here, we examine velocity behaviour in different simulations as pressure varies.

In 2010, water pressure values ranged from 0.51 to 0.8 MPa, with P_wp ranging from 1.1 to 1.7, and velocities varied linearly with water pressure between 26 and 51 m yr⁻¹ (Fig. 10). However, for 2018, water pressure values ranged from 0.68 to 1.06 MPa, with P_wp again ranging from 1.1 to 1.7, and velocities varying from 40 to 182 m yr⁻¹.

The results demonstrate that for values of water pressure less than 0.9 MPa, velocity increases linearly with water pressure with a slope of ∼ 7 m yr⁻¹ per 0.1 MPa. However, at values greater than 0.9 MPa, surface velocity scales exponentially with water pressure whereby small increments in water pressure increase velocity dramatically (Fig. 10).

While it is widely recognized that subglacial water pressure and ice velocity exhibit a nonlinear relationship (e.g., Cuffey and Paterson, 2010), a universal formulation of this interaction has yet to be established. However, despite all the limitations inherent in the simplicity of the model used, the results of this study may give us an idea of what this relationship might look like. Furthermore, these findings support previous work suggesting that there is a critical threshold in the relationship between water pressure/effective pressure, and ice velocity beyond which there is a shift from an approximately linear to a highly nonlinear relationship (Iken and Bindschadler, 1986; Sugiyama et al., 2011).

Figure 10 Mean modelled surface velocity within 1.5 km from the glacier front as a function of basal water pressure for both 2010 and 2018 and constrained by different values of the water pressure parameter P_wp.

4.3 Influence of front retreat and thinning on surface velocities

As outlined in Sect. 3.4, we undertook a suite of experiments in order to isolate the respective effects of ice front retreat versus surface thinning on the ice dynamic response of the near terminus region (Table 2). Our results demonstrate that during the first period, 2010–2020, the observed thinning exerts a far greater influence on glacier acceleration than terminus retreat, resulting in near-terminus acceleration of up to ∼ 8 m yr⁻². In contrast, the front retreat, and the associated impact on buttressing and near-terminus stress-regime, produced only localized and minimal near-terminus acceleration of up to ∼ 2 m yr⁻²; as such, the influence of frontal retreat on ice motion was only 25 % of that driven by surface thinning for the period 2010–2020 (Fig. 11a). During the second period, the impact associated with front retreat becomes more important, with near-terminus acceleration of up to ∼ 6 m yr⁻², while the impact of surface thinning causes acceleration of up to ∼ 14 m yr⁻². The influence of front retreat therefore increases during this later period, but surface thinning clearly remains the dominant driver of acceleration during the period 2020–2030 (Fig. 11b). These findings suggest that changes in the stress-regime associated with thinning, rather than front retreat, play the dominant role in controlling the observed velocity evolution at Skaftafellsjökull during the transition from a land- to a lake-terminating system.

One potential issue with these experiments is the fact that the timescales on which these two processes operate can be different. Whilst thinning is a more gradual sustained process that operates most intensely during the summer months, the retreat of the glacier front is an episodic process often characterised by brief periods of intense iceberg calving. Nevertheless, we suggest these differences do not significantly affect the experiment as prescribed, because the modelled time scales under consideration, span a decade, and are therefore integrating change over a much longer timeframe than is typical of either of these processes.

Figure 11 (a) Surface acceleration of the glacier tongue for EXP1 (left) and EXP2 (right) for the 2010–2020 period. (b) Detail of surface acceleration map of the glacier tongue for EXP1 (left) and EXP2 (right) for the 2020–2030 period.

4.4 Frontal ablation

In order to estimate changes in frontal ablation at Skaftafellsjökull through time, ice discharge through a flux gate located ∼ 1.5 km upglacier from the 2010 terminus (Fig. 12b) was calculated for each of 2010, 2018, 2020 and 2030 (Fig. 12a). As the model assumes ice flow is primarily driven by sliding, we prescribe f=1 , and α_i=1 assuming the flux gate is perpendicular to the flow. The results from the modelled estimates of ice flux show an approximately linear increase between 2010 and 2030, with a factor three increase in flux (Fig. 12a), driven by the near-terminus acceleration, and despite glacier wide surface thinning (Fig. 4c) which results in a reduction in the cross-sectional area at the flux gate.

Figure 12 (a) Estimates of ice flux through a flux gate, located ∼ 1.5 km upglacier from the 2010 terminus, between 2010 and 2030. The best fit regression line (orange) reports a statistically significant increase (R² = 0.91) in flux with time. (b) Location of the flux gate and the front positions for the different years.

The annual mean surface mass balance (or surface ablation $\dot{b})$ below the fluxgate for the period 2010–2018 was estimated using RACMO data, taking the cells from the RACMO grid located below the fluxgate for each year from 2010–2018, giving a mean value of mass loss of −11.8 × 10⁶ m³ w.e. yr⁻¹. The mean annual frontal ablation rate between 2010 and 2018 is subsequently calculated following Eq. (5). Here, the annual mean ice flux through the flux gate (−18.3 × 10⁶ m³ w.e. yr⁻¹ (F_fg), negative because it is a mass loss) between 2010 and 2018 is added to the annual mean ice-volume lost due to frontal retreat over the same period (−37.2 × 10⁶ m³ w.e. yr⁻¹); the annual mean surface mass balance below the fluxgate (−11.8 × 10⁶ m³ w.e. yr⁻¹) is then subtracted giving an annual value of frontal ablation of −43.7 × 10⁶ m³ w.e. yr⁻¹ (Table 3). All ice volumes are converted to water equivalent assuming a density of 900 kg m⁻³ (following Osmanoglu et al., 2014).

Annual mean surface ablation for the whole glacier for the period 2010–2018 was retrieved from RACMO data giving a value of −74.5 × 10⁶ m³ w.e. yr⁻¹. Total mean annual ablation for the whole glacier between 2010–2018 is therefore −118.2 × 10⁶ m³ w.e. yr⁻¹ of which frontal ablation, driven by terminus acceleration and retreat, represents ∼ 37 % (Table 3). Frontal ablation thus emerges as a critical driver of mass loss during the study period emphasizing the growing importance of the evolving lake-induced dynamic processes, relative to surface mass balance, in governing the glacier's mass budget.

Table 3 Estimates of the different annual mean flux and ablation terms for the period 2010–2018 including, in bold, annual mean frontal ablation and total ablation.

4.5 Limitations of our modelling protocol

As with all modelling studies, there are potential limitations with the work that merit consideration. One possible limitation of the model concerns the fact that the free parameters have been treated as constant across the domain, even though they are known not to be. This could mean that spatial variations in velocity are not adequately represented in the model. This may explain, for example, why the more complex spatial variations in the observed velocities shown in Fig. 7 are not reflected in the modelled velocities, which are smoothed and more uniform. We would contend however that this limitation does not significantly affect the key findings of the model which are focused more on the broad velocity trends in the near terminus region and along the glacier centreline (e.g. Fig. 8). In addition, by applying our tuned parameters outside of the 2010–2018 range, we are implicitly assuming the parameters are constant in time, or rather that any change is unimportant. While there are insufficient observations to constrain the temporal behaviour of these parameters, it must be kept in mind that unresolved changes in ice stiffness and bed or hydrological characteristics could bias our results.

Another limitation concerns the timing of the derived ice velocities, and their uncertain degree of representation of annual means. The available ice velocities were from August and September, that is, from late summer. These late-summer velocities can be considered a reasonable approximation of annual velocities, because the high velocities typical of early summer have already passed, but they have not yet dropped to the low velocities typical of autumn and winter (Davison et al., 2019; Yang et al., 2025). It is clear however that the seasonal variability of velocities could limit the applicability of the velocity data used.

Finally, the use of the SSA implies that vertical stresses are being neglected in favour of horizontal ones, which means that ice deformation is not being taken into account and the only movement considered is sliding. This could lead to significant errors in a glacier with a large area of slow flow, dominated by deformation, as might be expected in the accumulation zone of Skaftafellsjökull. However, this cannot be confirmed by the observed velocities, which are not well resolved in this area (Fig. 5), likely because of their slow flow and lack of surface features. Nevertheless, in the fast-flowing areas of the glacier and particularly in the lower glacier below the icefall, is the region of most importance to this study, the glacier's movement is indeed well modelled by the SSA velocities (Fig. 8).

4.6 Broader implications

The study supports the growing body of evidence demonstrating that lake-terminating glaciers are governed by a complex suite of mechanical feedbacks that distinguish them from land-terminating systems (e.g., Benn et al., 2007b). The proposed model highlights how the hydraulic connection, understood as a connection in the sense of a pressure influence rather than water exchange, between the proglacial lake and the subglacial drainage system significantly reduces basal friction in the near-terminus region when compared to a land-terminating glacier. For regions with glaciers in potential transition from land to lake-terminating, such as Iceland (Aðalgeirsdóttir et al., 2020) and Greenland (Patton et al., 2016; Shugar et al., 2020), this implies that as proglacial lakes expand due to ongoing retreat in a warming climate, we can expect an acceleration in ice flow and enhanced mass loss as a result. The reduction in back-stress and the increase in basal water pressure, often approaching overburden pressure, facilitate hydro-fracturing and calving, creating a self-sustaining cycle of dynamic thinning and terminus recession (Sugiyama et al., 2011). However, as this model represents basal water pressure in a very simple way, with a constant parameter value across the entire glacier, further research is needed, using models that represent basal hydrology in a more complex manner, to investigate subglacial hydrology in lake-terminating glaciers and determine the extent to which our simplification is problematic.

Specifically for the Icelandic context, such as the Vatnajökull Ice Cap margin, the findings suggest that the depth of proglacial lakes are primary drivers of glacier stability. As these lakes deepen as glaciers retreat in to overdeepened basins, the increased buoyancy will likely trigger more frequent and higher-magnitude calving events (Baurley et al., 2020; Benn et al., 2007b). This study provides an important framework for predicting the future evolution of Icelandic ice caps, where the rapid formation of new lakes is currently transforming relatively stable land-terminating margins into highly sensitive lake-terminating ones, potentially leading to faster-than-predicted mass loss (Schomacker, 2010).

Both the observations and model demonstrate that frontal ablation plays an important role in total ablation at Skaftafellsjökull, highlighting the influence of the lake on glacier mass balance. Given the increasing number and area of proglacial lakes worldwide due to terminus retreat (Zhang et al., 2024a), these processes have the potential to become increasingly important to global glacier mass balance and thus sea-level rise.

The implications for Greenland may be particularly significant given the extent of the ice-sheet's land-terminating margin and since its peripheral glaciers are retreating and ice marginal lake area is increasing (Zhang et al., 2024a). Furthermore, bedrock over-deepenings are widespread beneath the ice-sheet margins, especially where outlet systems follow deeply incised troughs (Patton et al., 2016), hence the transition to lake-terminating conditions in these systems could serve as a major catalyst for future accelerated discharge of inland ice. Our study demonstrates that lake-terminating glaciers experience distinct boundary conditions that impact glacier dynamics and associated mass loss, highlighting the need to consider the influence of proglacial lakes in future projections of Greenland's contribution to global sea level.

Despite the valuable insights provided by this modelling framework, several inherent limitations highlight critical avenues for future research. First, the parameterization of subglacial hydrology via a spatially uniform water pressure parameter, P_wp, simplifies a highly localized and dynamic system. Future iterations should incorporate fully coupled, distributed-channelized subglacial drainage models to capture spatial variations in effective pressure and to compare these with simpler models, such as the one presented in this paper. Second, the reliance on late-summer velocity datasets introduces a degree of uncertainty as to whether those velocities are representative of the annual velocity. Integrating annual velocity observations will therefore be essential to refine model tuning. Finally, while the depth-integrated Shallow Shelf Approximation (SSA) successfully captures rapid terminus sliding, upgrading to higher-order or Full-Stokes mechanics will better resolve three-dimensional stress fields and vertical shear variations in the transition zones up-glacier.

5 Conclusions

This study presents a computationally efficient modelling framework to analyse the transition of an Icelandic outlet glacier, Skaftafellsjökull, from a land- to a lake-terminating glacier. Through the tuning of three parameters, the model successfully reproduces the observed velocity patterns and captures the pronounced near-terminus acceleration between 2010 and 2018, reflecting the growing influence of the proglacial lake. The results indicate a threshold behaviour between effective pressure and ice velocity, whereby small decreases in effective pressure beyond a critical value result in strong accelerations, consistent with previous empirical observations (Iken and Bindschadler, 1986; Sugiyama et al., 2011). Furthermore, frontal ablation at the lake-terminating margin, accounted for roughly approximately 37 % of the total ablation between 2010 and 2018, highlighting the importance of the proglacial lake in accelerating mass loss via ice-dynamic processes.

Sensitivity experiments demonstrate that surface thinning exerts a much stronger control (3–4× greater) on near terminus velocity than terminus retreat. However, the influence of front retreat on terminus acceleration increases with time, and is expected to be greater in the future, as the glacier terminus retreats into deeper water, in a similar process to that observed across large parts of West Antarctica where the bedrock of marine-terminating glaciers deepens inland resulting in grounding line retreat and ice sheet instability (e.g., Rignot et al., 2014). Conversely, considering behaviour over longer timeframes, many lake-terminating mountain glaciers, as well as those in Greenland, will eventually, as they continue to retreat, become disconnected from their proglacial lakes and once again terminate on land, with consequences for ice dynamics.

Overall, our results provide further evidence highlighting the sensitivity of glacier dynamics to evolving proglacial lake conditions and the associated impact on near-terminus effective pressures. The modelling framework developed here provides a valuable tool for simulating these complex interactions in a computationally efficient manner. By capturing threshold behaviour in the effective pressure–velocity relationship, it contributes to an improved understanding of how lake-terminating glaciers will respond to ongoing climatic warming and the associated expansion of proglacial lakes. These insights could be essential for refining projections of future ice mass loss and its contribution to global sea-level rise.