Long‐term impact of the proglacial lake Jökulsárlón on the flow velocity and stability of Breiðamerkurjökull glacier, Iceland

Proglacial lakes are becoming ubiquitous at the termini of many glaciers worldwide due to continued climate warming and glacier retreat, and such lakes have important consequences for the dynamics and future stability of these glaciers. In light of this, we quantified decadal changes in glacier velocity since 1991 using satellite remote sensing for Breiðamerkurjökull, a large lake‐terminating glacier in Iceland. We investigated its frontal retreat, lake area change and ice surface elevation change, combined with bed topography data, to understand its recent rapid retreat and future stability. We observed highly spatially variable velocity change from 1991 to 2015, with a substantial increase in peak velocity observed at the terminus of the lake‐terminating eastern arm from ~1.00 ± 0.36 m day−1 in 1991 to 3.50 ± 0.25 m day−1 in 2015, with mean velocities remaining elevated from 2008 onwards. This is in stark comparison to the predominately land‐terminating arms, which saw no discernible change in their velocity over the same period. We also observed a substantial increase in the area of the main proglacial lake (Jökulsárlón) since 1982 of ~20 km2, equating to an annual growth rate of 0.55 km2 year−1. Over the same period, the eastern arm retreated by ~3.50 km, which is significantly greater than the other arms. Such discrepancies between the different arms are due to the growth and, importantly, depth increase of Jökulsárlón, as the eastern arm has retreated into its ~300 m‐deep reverse‐sloping subglacial trough. We suggest that this growth in lake area, forced initially by rising air temperatures, combined with the increase in lake depth, triggered an increase in flow acceleration, leading to further rapid retreat and the initiation of a positive feedback mechanism. These findings may have important implications for how increased melt and calving forced by climate change will affect the future stability of large soft‐bedded, reverse‐sloped, subaqueous‐terminating glaciers elsewhere. © 2020 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd


Introduction
Continued and more intensive global climate warming, particularly over the last decade, is driving patterns of glacier recession and dynamics, and consequently it is now widely established that almost all glaciers worldwide are undergoing widespread retreat (Björnsson et al., 2013;Zemp et al., 2015;Glasser et al., 2016). This has important consequences for their meltwater contribution to global sea level rise (s.l.r.) (Huss and Hock, 2015;Cazenave and WCRP Global Sea Level Budget Group, 2018;Rossini et al., 2018;Zemp et al., 2019), as well as for regional hydrology due to the strong control glacier meltwater has on modulating down-glacier streamflow, which in turn affects freshwater availability, hydropower operations and sediment transport (Immerzeel et al., 2014;Huss and Hock, 2018;Shannon et al., 2019).
There are several different mechanisms which it is believed are causing the widely observed dynamic changes occurring at land-, lake-and marine-terminating outlet glaciers. For example, some land-terminating glaciers globally are slowing down due to a reduction in ice thickness and surface slope, which has caused a corresponding decrease in the driving stress, and consequently ice flow over recent decades (Heid and Kääb, 2012;Mernild et al., 2013;Dehecq et al., 2019). In contrast, increased surface melt may lead to increased seasonal velocities, via basal lubrication and till deformation (e.g. Boulton et al., 2001;Zwally et al., 2002;Schoof, 2010;Bartholomew et al., 2011;Hart et al., 2011;Andrews et al., 2014;Bougamont et al., 2014;Hart et al., 2019a). This process can vary considerably over the melt season, however, depending on several different factors (Bartholomew et al., 2012;Tedstone et al., 2013;Hart et al., 2019b). These include: (i) the nature of the basal drainage system in place (e.g. inefficient cavities vs. efficient channels) (Röthilsberger, 1972;Schoof, 2010); (ii) both the amount and the rate at which meltwater reaches the bed (Meierbachtol et al., 2013;; and (iii) the area covered by one of these drainage systems at the expense of the other, and the interconnectivity between them Tedstone et al., 2015). This means that the seasonal response of each glacier to surface meltwater varies between individual glaciers and between glaciated regions (e.g. Andrews et al., 2014;Tedstone et al., 2015). Equally important, however, is the observed increase in the number of glaciers worldwide which are responding rapidly to sub-aqueously driven dynamic changes (e.g. Hart et al., 2011;Carrivick and Tweed, 2013;Sakakibara et al., 2013;King et al., 2018).
Calving fluxes, and the behaviour of calving glaciers in general, are strongly controlled by ice velocity (Meier and Post, 1987;Howat et al., 2005;Benn et al., 2007b;Pritchard et al., 2009;Sugiyama et al., 2011;Sakakibara et al., 2013). This is because the majority of calving events are a consequence of crevasse propagation in response to longitudinal stress caused by fast ice flow (Meier and Post, 1987;Naruse and Skvarca, 2000;Van der Veen, 2002;Benn et al., 2007b). As a glacier retreats into deeper water (e.g. down a reverse slope), the water depth and therefore the calving flux increase as a result (Schomacker, 2010). This is because glaciers flow faster when entering deeper water, due to the inverse relationship between effective pressure and basal drag, leading to rapid ice flow, thinning and retreat (Howat et al., 2007;Benn et al., 2007a), and such a relationship has been observed for a number of glaciers worldwide (e.g. Brown et al., 1982;Pelto and Warren, 1991;Naruse and Skvarca, 2000;Warren and Kirkbride, 2003;Nick et al., 2007).
Many of Iceland's glaciers have been retreating rapidly sincẽ 1990 in response to increasing mean air temperatures, with many more showing a particularly heightened rate of retreat over the last decade (Sigurđsson et al., 2007;Bradwell et al., 2013). Such a response can be attributed to a mean annual temperature increase of 1°C since 2000, which is three to four times higher than the Northern Hemisphere average over the same period (Jones et al., 2012), as well as a shift in atmospheric and oceanic circulation patterns around Iceland (Björnsson et al., 2013;Foresta et al., 2016). This has resulted in~9.5±1.5 Gt a À1 of mass loss between the mid-1990s and 2010, and a loss of 5.8±0.7 Gt a À1 between 2010-11 and 2014-15, equating to a s.l.r. contribution of~0.03 and 0.016 mm a À1 , respectively (Björnsson et al., 2013;Foresta et al., 2016). This retreat has been accompanied by a sharp increase in the formation and growth of proglacial lakes at the termini of many of the southerly flowing outlet glaciers of the country's main ice caps (Hannesdóttir et al., 2015;Dell et al., 2019). The effect these lakes are having on these glaciers is above and beyond that of climate alone (Carrivick and Tweed, 2013;Staines et al., 2015), because their development causes accelerated terminus retreat by influencing ice dynamics and velocity through calving (Haeberli et al., 2016;Nie et al., 2017;King et al., 2018). Furthermore, the size and number of such lakes is projected to increase in the future due to continued climate warming and glacier retreat, with the resulting effect on glacier dynamics in the region likely to be significant (Björnsson et al., 2001;Flowers et al., 2005;Schomacker, 2010).
One glacier which has seen such proglacial lake formation and development at its margin is Breiðamerkurjökull, a large soft-bedded temperate glacier in south-east Iceland. During its Little Ice Age (LIA) advance, the glacier excavated a 200-300m-deep proglacial trough under one of its larger, more dynamic ice lobes (Björnsson, 1996;Björnsson et al., 2001), leading to a marked overdeepening under the glacier. Since this maximum, however, the glacier has retreated rapidly (Hannesdóttir et al., 2015;Guðmundsson et al., 2017), but as a result of this overdeepening this dynamic lobe of the glacier is now retreating down a reverse slope below sea level (Björnsson, 1996;Guðmundsson et al., 2017).
Although there have been several studies undertaken on Breiðamerkurjökull in recent years (e.g. Schomacker, 2010;Voytenko et al., 2015;Guðmundsson and Björnsson, 2016;Storrar et al., 2017), these studies have tended to focus on small-scale velocity changes occurring over short time scales (e.g. Voytenko et al., 2015), or on one aspect of glacial dynamic change such as the influence of glacier topography on ice dynamics (e.g. Storrar et al., 2017).
Therefore, the aim of this study is to investigate recent velocity changes at Breiðamerkurjökull, south Iceland using satellite remote sensing to understand how the glacier is responding to recent change and evaluate the different key glaciological and morphological forcing factors. We quantify decadal changes in glacier velocity since 1991 in relation to frontal retreat, lake area change, bedrock topography and ice surface elevation change. In doing so, we have examined glacier velocity at a spatial and temporal resolution not yet attempted on Breiðamerkurjökull, enabling us to fully quantify the longterm, glacier-wide changes in dynamics that have occurred over recent decades. The conclusions from this site may be used to predict the response of other soft-bedded glaciers with subglacial reverse slopes and aquatic margins (e.g. in Iceland and Alaska) to future warming.
The glacier had a combined area of ∼⃒ 906km 2 in 2010, with an elevation range from~10to 1700ma.s.l. (Guðmundsson et al., 2017). The glacier is composed of three main lobes, or 'arms' (after Guðmundsson et al., 2017), separated by two large medial moraines. The larger, eastern arm (Norðlingalaegðarjökull) drains the large ice dome of Breiðabunga deep within the ice cap (Evans and Twigg, 2002). The central arm (Esjufjallajökull) flows out from the two large nunataks Máfabyggðir (1440ma.s.l.) and Esjufjöll (1770ma.s. l.), with the latter also giving the name to the large medial moraine (Esjufjallarönd) that separates the eastern and central arms (Guðmundsson and Björnsson, 2016). The western lobe (Máfabyggðajökull) collects ice from several glaciers on the eastern flanks of Öraefajökull before draining south of Máfabyggðir, divided from the central arm by the medial moraine Mávabyggðarönd (Guðmundsson and Björnsson, 2016;Guðmundsson et al., 2017). The most western lobe is composed of two arms which we have defined as western A and B (Figure 1). Two lakes of differing sizes are also present at the glacier terminusthe large~24km 2 , 300m-deep Jökulsárlón (Voytenko et al., 2015) adjacent to the eastern arm, and the smaller 5.8km 2 Breiðárlón adjacent to western A. Most of glacier bed sits at, or just above, sea level, however 2648 N. R. BAURLEY ET AL.
a small, shallow trough runs back from the margin of western A, while a large, deep trough is found under the eastern arm, which extends back~20km from Jökulsárlón into the glacier interior (Björnsson, 1996). This region has also seen a 1.5°C rise in the mean annual air temperature between 1980 and 2017, based on data collected by the Icelandic Met Office from their station at Fagurhólmsmýri (63°52′ N, 16°38′ W), located 5km from Breiðamerkurjökull at an elevation of 16ma.s.l.  Table 2. For these analyses, data were downloaded from their respective archives in image precision (IMP) mode for ERS-1/2, interferometric wide-swath (IW) mode for Sentinel-1 and strip-map (SM) mode for TerraSAR-X. Where possible, data acquired in September were downloaded to allow direct comparisons to be made between years, however, this was not always possible due to the availability of data and so in these instances the nearest possible dates to September were selected. Acquisition dates for each , the centrelines of each arm (deep red), the start and end point of each centreline (e.g. W to W′), and the location and area (as of 2019) of Jökulsárlón and Breiðárlón (light green). Glacier outlines obtained from the GLIMS database. Centrelines derived manually using the geometry of the outlines. Nunatak outlines taken from Guðmundsson et al. (2017). (d) Bedrock topography of Breiðamerkurjökull, which extends down to À300m below sea level at its deepest (darkest green). Interpolated from contours provided in Björnsson (1996). Lake areas and glacier outlines as before. Background image is a Sentinel-2 acquisition (5, 3, 2 RGB) from 06/07/2019. [Colour figure can be viewed at wileyonlinelibrary.com] 2649 LONG-TERM FLOW VELOCITY AND STABILITY OF BREIÐAMERKURJÖKULL, ICELAND SAR image used in the velocity analyses, as well as the time between each pair for each year, are given in Table 3. Overall, the dataset spans 1991-2018.

Glacier-wide velocities
All data were analysed and processed using the Sentinel Application Platform (SNAP) (ESA, 2016), the freely available ESA-developed software toolbox for analysing ESA products (//step.esa.int/main/download/). All sets of data were processed using the well-established offset-tracking methodology (Lu and Veci, 2016). Precise orbital vectors were first applied to each individual SAR acquisition to provide accurate satellite position and velocity information for each product (Serco Italia SPA, 2018), before sub-setting each image to the area of Breiðamerkurjökull to streamline processing. Digital elevation model (DEM)-assisted co-registration was then undertaken on pairs of images for each period using the LiDAR DEM of Iceland (Landmaelingar Íslands, 2016). The DEM has a resolution of 10×10m and was updated in 2016 to include elevation data for all glaciers in Iceland (Landmaelingar Íslands, 2016). The LiDAR DEM was preferred to ArcticDEM, as the former was deemed to be of better quality and, therefore, would be capable of ensuring more accurate image co-registration (Heid and Kääb, 2012). This process was performed on each image pair before offset tracking was undertaken on the resulting stacked images. All images for each sensor type were processed using the exact same parameters to ensure spatial coverage did not vary between images, allowing direct comparisons to be made between sensor types (Table 4).
Offset tracking estimates the movement of glacier surfaces between master and slave images in both the slant range and azimuth direction through cross-correlation on selected ground control points (GCPs) Nagler et al., 2015;Fahnestock et al., 2016;Serco Italia SPA, 2018). The glacier velocity is computed based on the offsets estimated by the cross-correlation algorithm before the final velocity map is generated through interpolation of the velocities computed on the GCP grid (Lu and Veci, 2016). Displacement vectors larger than 5m day À1 or with a correlation coefficient less than 0.1 were removed. The displacement vectors were then averaged over a 20×20 pixel grid to create the velocity raster. The holes left by the anomalous values were then filled by replacing each missing point with a new offset computed by locally weighted average (Serco Italia SPA, 2018). The final velocity maps were subsequently terrain-corrected using the LiDAR DEM of Iceland to convert between slant-range and ground-range coordinates. To allow for a more robust comparison of the velocity data between years to be made, we present mean velocity measurements taken at 1km intervals along the glacier centrelines (shown in Figure 1). Furthermore, although we used images from approximately the same period each year, we wanted to investigate the effect of seasonality on our data, and test whether the velocity change we observed for the eastern arm reflected actual change, or simply seasonal variation. To do so, we utilized several TerraSAR-X images covering 24 July-11 November 2008 and processed them following the same method as was used for the velocity analyses.

Lake area and terminus position changes
To assess lake area change for both Jökulsárlón and Breiðárlón, the individual lake areas were manually digitized in ArcGIS 10.7 for 11 time steps between 1982 and 2018 using a  1982,1988,1989,1991,1992,1994,1998  combination of black and white orthorectified aerial photographs and Landsat and Sentinel-2 imagery. Information about each of the aerial photos used in this study is given in Table 5. Manual digitization was preferred to automated normalized difference water index (NDWI) image classification because turbid water can often have varying NDWI values, resulting in lakes being classified as glacier ice and thus hindering our ability to accurately delineate lake area through time (Racoviteanu et al., 2009;Gjermundsen et al., 2011). The individual photos were first orthorectified and then georeferenced to a base image with known coordinates. This was undertaken in ArcGIS by marking on stable points within each individual aerial photo and then finding the same points in a Landsat OLI acquisition from 2018. A Landsat OLI acquisition was utilized as the base image due to its low geolocation errors (Dell et al., 2019). All aerial photos were georeferenced with a root mean square error of less than 0.5 pixels. Lake areas were then digitized at a scale of 1:6000 for the aerial photographs and 1:10000 for the Landsat and Sentinel-2 imagery. These scales allowed the accurate mapping of the lake areas and prevented pixelated images hindering reliable interpretation (Lovell, 2016;Dell et al., 2019). Channels that exited the lake (such as the one flowing from Breiðárlón into the neighbouring Fjallsárlón) were ignored during the digitization at the point where the channel began to form.
Terminus position changes for each of the main arms of Breiðamerkurjökull were assessed using the same imagery, time steps and digitizing scale used to quantify lake area change. At each time step, the position of the terminus for each of the main arms was manually digitized. To calculate the positional change through time for each arm, the rectilinear box method was implemented (Moon and Joughin, 2008;Howat and Eddy, 2011), which quantifies the change in area between different terminus positions using a fixed-width rectilinear box drawn over the glacier trunk. This was then converted to a 1-D figure by dividing the change in area by the change in terminus width (Howat and Eddy, 2011;Lea et al., 2014). This method was chosen for this study due to its ability to account for asymmetrical changes at the calving front (Lea et al., 2014;Larsen et al., 2016;Dell et al., 2019), while other methods, such as automatic ice classification of Landsat images, can be influenced by turbid water, and it cannot be used for aerial photos, which would hinder its use in this anal-ysis (Gjermundsen et al., 2011). When assessing the change occurring at the calving front of the eastern arm and western A, the width of the box used encompassed the maximum delineated width of the lake-terminating portion of the front (1992 for Jökulsárlón, 1989 for Breiðárlón) rather than the whole terminus of both arms. This ensured that the captured rate of positional change actively related to what was occurring at the calving front.

Ice surface elevation change
Changes in the ice surface elevation of Breiðamerkurjökull were examined using the freely available ArcticDEM dataset distributed by the Polar Geospatial Centre (Porter et al., 2018) (https://www.pgc.umn.edu/data/arcticdem). This dataset provides digital surface models (DSMs) at a spatial resolution of 2m for areas north of 60°from 2011 in most regions (Morin et al., 2016). Data are typically downloaded as 17×110km strips (Barr et al., 2018). In the case of Breiðamerkurjökull, the data coverage is particularly good, with several DSMs available covering varying extents of the glacier for different years. We chose two DSMs for this analysis, one from 05/10/2012 and the other from 20/10/2017 as both DSMs covered the entire glacier area, while also affording us the opportunity to investigate surface elevation change associated with glacier recession. The DSMs were linearly co-registered following the method of Nuth and Kääb (2011).

Uncertainty analyses
Velocity uncertainty To determine the uncertainty of our velocity measurements, we measured displacements over terrain that we regarded as stable (Robson et al., 2018). We carried out this analysis over two different areas ( Figure 2) before calculating the combined stochastic error of both. The first was restricted to an area~1km away from the land-terminating portion of Breiðamerkurjökull, on a portion of presumably stable ground situated between both proglacial lakes that was free from shadow and had little or no surface slope. The second was placed on sloping, but stable, ground to the right of Breiðamerkurjökull to assess the accuracy of the co-registration process. Generally, we find that mean stochastic uncertainties for each satellite sensor are somewhat similar. Calculated values are greatest for the ERS-1/2 data, with a mean error of ±0.19m, while for TerraSAR and Sentinel-1 the error is ±0.14 and ±0.11m, respectively, indicating that the levels of uncertainty we estimate are not greater than the change in velocity over the duration of our study.

Digitization uncertainty
To quantify the uncertainty when digitizing the lake areas and glacier frontal positions, the areas of Jökulsárlón and Breiðárlón in 1982, 2014 and 2018 were repeatedly digitized 10 times at the same scale as used in the initial analyses. Key statistical values, including the coefficient of variation (CV) and standard uncertainty, or standard error (SE) were then calculated to allow comparisons to be made between the three different sets of aerial imagery (Table 6). Importantly, the calculated uncertainties for all image types are not greater than the change observed in our analyses.
Surface elevation change uncertainty The uncertainty of the elevation changes was calculated by combining the standard error with the mean vertical bias over the terrain that was assumed to be stable (i.e. on non-glacier terrain) (Nuth et al., 2012;Robson et al., 2018). The SE takes into account standard deviation over stable terrain (SD STABLE ) and the number of independent pixels included in the DEM differencing (n): n is determined by the original number of pixels in the DEM differencing analysis (N tot ), the pixel size (PS) and spatial autocorrelation (d), which we assumed to be 600m: DEM differencing uncertainty was then calculated by combining the standard error with the mean vertical bias on stable ground, following the method of Robson et al. (2018). This gave an elevation error on stable ground of ±0.45m over the period 2012-2017, equating to 0.09myear À1 , indicating that our elevation assessment is highly accurate.

Glacier-wide velocities
Surface velocity 1991-2018 We observe spatially variable velocity change for Breiðamerkurjökull over the period 1991-2018 ( Figure 3). The calculated velocity observations (Figure 3) show that the highest velocity values for the study period are consistently found for the eastern trunk, particularly from 2008 onwards.
Maximum values increased from~0.65 ±0.15m day À1 in 2001 to~3.50 ±0.25mday À1 by 2015, before dropping tõ 1.72 ±0.11m day À1 in 2018. Such a trend is also observed in the full velocity dataset (see Figure S1 in the online Supporting Information), despite maximum velocities first decreasing slightly between 1991 and 1996, from~1.00 ± 0.36 to~0.55 ±0.11m day À1 , before then increasing again until 2015.
The mean velocities for the years 2001-2018, calculated using the glacier flowlines, indicate that the greatest change in mean velocities over the period is observed for the eastern arm (Figure 4; full dataset shown in Figures S2-S5 of the online Supporting Information). Between 2001 and 2005, mean velocities over the whole trunk of the eastern arm remained stable, apart from a slight peak in velocities near the calving front of 0.58 ±0.15and 0.44 ±0.12mday À1 observed for 2001 and 2005, respectively. From 2008 onwards, however, the mean velocity remained elevated and increased year on year, particularly in the near-terminus region. Peak near-terminus velocities in 2008 were 1.58 ±0.09mday À1 , jumping to 1.63 ±0.11 mday À1 in 2012, increasing slightly to 1.64 ±0.15mday À1 in 2015, before peaking at 1.68 ±0.11m day À1 in 2018. No such dynamic changes were observed at the other three arms over the same period.
Mean velocities for the central arm remained stable between 2001 and 2018, with very little yearly variation observed. Indeed, apart from a peak in near-terminus velocities of 0.49 ±0.15mday À1 in 2001, the observed general trend from 2001 onwards is one of decreasing velocities down-glacier right up to the margin. Indeed, a similar scenario is also found for both  western A and western B, with the former having the lowest overall velocities for any arm of Breiðamerkurjökull, particularly from 2008 onwards.

Seasonality analysis
The results of our seasonality analysis ( Figure 5) indicate that, despite there being some variation (which we infer to be seasonal variation), overall the velocity signal is stable and does not show any large variations from the mean, which is particularly true for the terminus region (>22km along transect). Alongside this we have also plotted the velocity series covering our study period to investigate how these differ from the data for 2008. This clearly shows that in 2012 (the next closest data series to 2008), the velocity values are not significantly different to 2008, apart from in the near-terminus region where velocities are greater. Conversely, the data from 1991 to 2005 sit well outside the data margins from 2008, indicating that this data is significantly different from the 2008 data. Furthermore, it also reinforces the inference that the most significant changes are occurring in the near-terminus region (i.e. with increasing velocities from 2008). This indicates that the velocity change observed in 2008, and in the other years of our analysis, likely reflects actual long-term change, and not natural seasonal variability.

Lake area and frontal position change
The areas of both Jökulsárlón and Breiðárlón have increased over the research period, although at vastly different rates (Figure 6). The area of Jökulsárlón increased from 7.58to 27.20km 2 (Figure 7), an overall increase of nearly 270% between 1982 and 2018. This equates to an overall lake growth since 1982 of~0.55km 2 year À1 . In comparison, the area of Breiðárlón increased only very slightly, from 5.39km 2 in 1982 to 5.84 km 2 in 2018 (Figure 7), equating to an increase of 8.36% (the uncertainty values for this analysis were <1%, see Table 6).
Over the study period, the land-terminating arms of the glacier seemed to be retreating at a steady rate between 1982 and 2018, with these arms all displaying similar rates of retreat from 1982 onwards (Figure 8). However, while the retreat of the lake-terminating portion of the eastern arm between 1982 and 1992 seemed to be occurring at a similar rate to that of the rest of Breiðamerkurjökull, after this point rapid changes started to occur. Overall, this portion of the eastern arm has retreated by~3.50km, or 92.69myear À1 (Figure 7). Particularly rapid retreat occurred between 2002 and 2014, at a rate of 470.25myear À1 . Similar retreat rates are not evident for any other arms of Breiðamerkurjökull, land-terminating or otherwise. Interestingly, for western A, which flows into Breiðárlón, we observe a retreat of~1.11km, equating to a rate of 31.09 myear À1 (Figure 7), the smallest overall change for the study period.

Ice surface elevation change
Between 2012 and 2017, significant thinning has occurred at the margin of the eastern arm where it terminates into Jökulsárlón, with 100-130±0.45m of negative surface elevation change observed (Figure 9). Furthermore, this thinning has occurred~2-3km back from the margin along the entire width of the calving front, with 60-80 ±0.45m of negative change observed, which is greater than the change observed at the other three lobes over the same period. Comparatively small elevation changes are observed for western A and B, with 20 ±0.45m of negative change observed for the former in particular, while the margin of the central arm thinned by~50 ± 0.45m. Meanwhile, slightly positive changes of 0-10 ±0.45m are observed over the higher reaches of the glacier.

Velocity change 1991-2018
We have presented the first multi-year and highly detailed velocity analyses of Breiðamerkurjökull covering almost three decades using high-resolution SAR offset-tracking techniques. Our data illustrate a clear increase in surface velocities for the eastern arm since 1991, with heightened velocity increases occurring from~2008 up to 2015, where values of up to 3.50 ±0.25m day À1 were recorded in the near-terminus region. This is in stark comparison to the other three arms of Breiðamerkurjökull, where we do not detect any change beyond the uncertainty of our methods. Such a contrast is likely due to the highly dynamic nature of the eastern arm and its close relationship with the growth of the Jökulsárlón lagoon, which has caused an increase in ice acceleration and retreat (Benn et al., 2007a;Schomacker, 2010). This relationship will be discussed in more detail below. Our data are comparable with the few other studies that have investigated the highly dynamic nature of the eastern arm over the same period. For example, Howat et al. (2008), using an automated surface-feature tracking method on repeat ASTER images acquired in the summer of 2005, found that in the near-terminus region the eastern arm had an annual flow of 550myear À1 . This equates to~1.5m day À1 , putting it directly in line with the values of 1-3.5m day À1 recorded in this study between . Similarly, Nagler et al. (2012 found surface velocities of <2m day À1 for the eastern arm using TerraSAR-X data from 2010. Interestingly, Voytenko et al. (2015), in their analyses of the calving front between 2011 and 2013 using terrestrial radar interferometry (TRI), found daily values of up to 5m day À1 , which is much faster than anything we have found in this study, particularly for our ERS and Sentinel data. However, the authors attributed such values to the high temporal resolution of the TRI, where the short averaging times for data acquisition were likely capturing short-lived dynamic phenomena that are smoothed in the longer time-averaged satellite data, even for high-resolution TerraSAR-X data (Voytenko et al., 2015).
It is worth discussing why we seemingly observe a decrease in velocities for the eastern arm between 2015 and 2018 (Figure 3). This decrease may be due to glaciological or geomorphological factors that have caused glacier slowdown, or it could indicate that the change in resolution between different satellite sensors is influencing the velocity estimates during this time, rather than the decrease being representative of wholesale glacier slowdown.
The observed slowdown of some lake-terminating glaciers elsewhere over recent years has been attributed to a combination of reduced driving stress and increased resistive stresses as the terminus of these glaciers narrow, in association with surface thinning (Pimentel et al., 2010;Adhikari and Marshall, 2012;King et al., 2018). Increased thinning, forced by sustained periods of ice-front retreat and glacier acceleration, causes the surface slope to flatten, reducing the driving stress up-glacier and, therefore, reducing overall ice velocities (Cuffey and Paterson, 2010;King et al., 2018). At Breiðamerkurjökull, an area of flat surface topography immediately behind the calving front of the eastern arm has been observed in several studies (e.g. Evans and Twigg, 2002;Guðmundsson and Björnsson, 2016;Storrar et al., 2017). Here, the surface slope is only thought to be~2°, whereas at the calving front the slope is~14° (Voytenko et al., 2015). This area of flat topography is thought to coincide with the deepest part of the subglacial trough (Evans and Twigg, 2002;Storrar et al., 2017), and its area has grown in recent years as the glacier continues to thin and retreat at an increased rate (Guðmundsson and Björnsson, 2016;Storrar et al., 2017). This reduction in surface slope may have caused a reduction in the driving stress since~2015, which, therefore, might explain why we see a decrease in velocities between 2015 and 2018. However, this is beyond the scope of this paper and as such will not be addressed in any further detail.
The change in sensor resolution could also explain the observed decrease in velocities, because errors can persist when using Sentinel-1 to map slower-moving areas of ice compared to when using TerraSAR-X imagery due to the impact of ionospheric noise and temporal decorrelation (Nagler et al., 2015). Equally, where sharp velocity gradients are present, such as at the terminus of Breiðamerkurjökull, the higher resolution of TerraSAR-X allows such changes in velocity to be tracked, whereas Sentinel-1 will tend to smooth over such gradients due to coarser image resolution (Joughin et al., 2018). Similar issues can also persist when using ERS data due to its much coarser 30m resolution, as well as the fact that ERS data can often give noisier velocity estimates due to weaker intra-pair correlation caused by surface changes over slower-moving areas (Joughin et al., 2018). To test this, we resampled our TerraSAR-X outputs to 10and 30m resolution before comparing these to the originals. We found that the only difference was a loss of detail from the higher-resolution outputs and more smoothing in the coarser-resolution outputs. No significant differences between the outputs were observed. This suggests that although the change in pixel resolution has had an impact on our velocity estimates, the impact overall is low, further highlighting the validity of our observations. We believe, therefore, that our data strongly indicate that substantial changes in the velocity of Breiðamerkurjökull have occurred between 1991 and 2018.
Lake growth and frontal retreat 1982-2018 We have calculated lake area change and cumulative retreat for Breiðamerkurjökull for a period covering nearly four decades, indicating substantial lake growth and frontal retreat, particularly since the turn of the century. Our growth rate for Jökulsárlón of 0.55km 2 year À1 since 1982 shows close agreement with the values of 0.5and 0.61km 2 year À1 put forward by Björnsson et al. (2001) and Schomacker (2010), respectively. The greater value found by the latter is due to the author investigating lake change over a 10-year period from 1999 to 2009, rather than from 1982 to 2018 as we have done in this study. However, if we average our data over a similar period (from 1998 to 2010, equating to 12years, not 10), we find an annual growth rate of 0.6km 2 year À1 , indicating good agreement with the findings of Schomacker (2010). Furthermore, our data are also coincident with the continued growth and development of other proglacial lakes in Iceland, particularly for the outlets of Vatnajökull, over the same period (Schomacker, 2010;Magnússon et al., 2012;Dell et al., 2019). Our results also share similarities with the results of those studies which have investigated long-term proglacial lake development and growth in the Himalaya (e.g. Nie et al., 2017;King et al., 2018). One such study by Zhang et al. (2019) examined the development of 23 proglacial lakes in the Poiqu River basin, central Himalaya, from 1964 to 2017. The authors found that lake area increased at a rate of 0.23 ±0.01km 2 between 1964 and~2010, with the most rapid increase in area observed for the~30years immediately following initial lake development. From 2010 onwards, growth rates decreased significantly to 0.03 ±0.03km 2 , which the authors attributed to the fact that many lakes had reached their maximum extent by this time, becoming bounded by an abrupt steepening of valley topography, which limited further growth (Zhang et al., 2019). As we have seen in this study, the area of Jökulsárlón has increased steadily since 1982, with particularly large increases in area occurring since~2006 (Figure 7). These heightened rates of lake growth are likely to be sustained over the coming decades while the glacier continues to retreat through its 22km-long, 300m-deep subglacial trough, and will only start to decrease and stabilize once the glacier retreats out of the trough and the lake becomes constrained by its FIGURE 7. Lake area change for Jökulsárlón and Breiðárlón since 1982 versus the cumulative retreat for all terminus regions of Breiðamerkurjökull over the same period. Note that the two lake-terminating portions of the glacier have the same marker style as the lake that they terminate into.
[Colour figure can be viewed at wileyonlinelibrary.com] 2657 LONG-TERM FLOW VELOCITY AND STABILITY OF BREIÐAMERKURJÖKULL, ICELAND topography (Flowers et al., 2005;Nick et al., 2007). We suggest, therefore, that Jökulsárlón has recently entered a stage in its development characterized by rapid increases in area, and that these increased rates will likely continue over the decades to come.
Our retreat data also shows strong agreement with other rapidly retreating outlet glaciers in Iceland (e.g. Bradwell et al., 2013;Hannesdóttir et al., 2015;Einarsson, 2017), particularly those that terminate in proglacial lakes (Schomacker, 2010;Dell et al., 2019). Our value of~3.50km of retreat since 1982 for the eastern arm is in the same order of magnitude as the 2km measured by Einarsson (2017) for the land-terminating portion of the eastern arm for the period 1982-2014. Interactions between lake growth, frontal retreat and ice velocity The clear increase in velocity observed at Breiðamerkurjökull over the study period is likely to have been caused by the close relationship that exists between lake growth, ice acceleration and terminus retreat (Storrar et al., 2017). The fact that the observed velocity increase is most dramatically associated with the lake growth at the margin of the eastern arm indicates the importance of Jökulsárlón in forcing the behaviour of this part of Breiðamerkurjökull. As the glacier has continued to retreat over the last few decades in response to warming air temperatures, it has receded into the large, 200-300m-deep subglacial trough that it excavated during its LIA advance (Figure 1) (Björnsson et al., 2001;Guðmundsson and Björnsson, 2016).
Consequently, the area, and importantly the depth, of Jökulsárlón has increased in response as meltwater infills the overdeepened basin (Schomacker, 2010). This increase in lake depth caused the glacier to flow faster as it entered deeper water due to the inverse relationship that exists between effective pressure and basal drag (Warren and Kirkbride, 2003;Howat et al., 2005;Benn et al., 2007a). This would have resulted in an increase in longitudinal strain rates and fracture propagation at the terminus, leading to a corresponding increase in calving activity and, therefore, the rate of retreat (Nick et al., 2007;Benn et al., 2007a;King et al., 2018). Further retreat down this reverse slope into the subglacial trough has caused the lake depth to increase further, causing more rapid flow acceleration, calving and retreat to occur, driving a positive feedback mechanism (Howat et al., 2007). This mechanism is exacerbated by the fact that ice flow from the interior of Breiðamerkurjökull cannot balance the substantial losses occurring at the terminus, leading to further rapid acceleration, retreat and lake growth (Björnsson et al., 2001;Nick et al., 2007).
Moreover, it has been stated by Storrar et al. (2017) that this arm of Breiðamerkurjökull is currently retreating over the deepest part of the subglacial trench, and as such, water depth is also likely to be at its deepest as a result (Figure 1). This means that, all other factors considered, ice velocity is likely to have reached its maximum over the last several years in response to this significant deepening of the lake. Indeed, we found daily velocity values of up to 3.50 ±0.25mday À1 at the calving front in 2015 where the trench, and therefore the water depth, is likely to be at its deepest (Björnsson et al., 2001;Voytenko et al., 2015). Therefore, while initial retreat at Breiðamerkurjökull was instigated by rising air temperatures following the LIA, once Jökulsárlón increased to a sufficient size and was able to start influencing frontal retreat and ice flow, we suggest that this then became the dominant mechanism in causing the rapid retreat and flow velocities observed at Breiðamerkurjökull since the turn of the 21st century. This supports previous work that argues that the rapid retreat and highly dynamic changes occurring at calving glaciers cannot be solely attributed to climatic forcing alone, and, in exceptional cases, occurs completely independent of climatic influences (e.g. Meier and Post, 1987;Rignot et al., 2003;Benn et al., 2007a;Howat et al., 2007;Sugiyama et al., 2011;Sakakibara and Sugiyama, 2014). We can contrast the behaviour of the eastern arm, with its rapid velocity and retreat rate, with the much slower and more stable western A arm, despite the presence of Breiðárlón at its margin. The area of Breiðárlón has remained approximately constant over our research period (Figures 6 and 7). We suggest that this stability is due to the topography of the bed. Beneath western A there is a much shallower trough (Figure 1), which sits above sea level for much of this distance (Guðmundsson et al., 2017), while the depth of the lake itself is only thought to be approximately 60m (Howarth and Price, 1967). Furthermore, the margin in this area is being pinned by large areas of bedrock (Schomacker, 2010) which are being exposed as the glacier retreats (Evans and Twigg, 2002). As a result, the lake has remained a similar depth, and the calving rate has remained relatively constant (Nick et al., 2007), resulting in a more stable ice margin.

Implications and wider importance
The findings of this study show how soft-bedded, reversesloped, aquatic terminating glaciers in both Iceland and elsewhere may respond to future climate warming. In Iceland, widespread glacier retreat caused by climate warming has been occurring since 1985 (Sigurđsson, 1998;Sigurđsson et al., 2007;Björnsson and Pálsson, 2008), and this rate of retreat has accelerated over recent years due to high summer melt forced by a 1°C increase in mean annual temperature since the turn of the century (Björnsson and Pálsson, 2008;Jones et al., 2012). The southern outlets of Vatnajökull, which have retreated rapidly since 1890, are particularly vulnerable to such warming conditions (Hannesdóttir et al., 2015). This is because many carved deep subglacial troughs in the underlying soft sediments during their LIA advance, and consequently now have beds that sit 100-300m at their deepest below the current elevation of the terminus (Björnsson et al., 2001;Björnsson and Pálsson, 2008;Magnússon et al., 2012). Now, as they rapidly retreat into these deep troughs, meltwater can infill the newly exposed topographic depressions, leading to the formation of proglacial lakes at their termini (Schomacker, 2010;Carrivick and Tweed, 2013). These lakes cause additional melting and mass loss through calving, facilitating further retreat and, therefore, further lake growth, initiating a positive feedback mechanism (Carrivick and Tweed, 2013;Hannesdóttir et al., 2015).
As we have shown, such a feedback mechanism has been occurring at Breiðamerkurjökull since the mid-2000s, and as a result there is the strong possibility that other southern outlets of Vatnajökull will undergo a similar pattern of retreat and mass loss in future (Schomacker, 2010;Dell et al., 2019). For example, particularly deep troughs exist under the neighbouring outlets Svínafellsjökull (~320m), Skaftafellsjökull (230m) and Fjallsjökull (210m), and all have been undergoing rapid retreat since the 2000s . Indeed, at Fjallsjökull, the adjacent glacier to Breiðamerkurjökull, such a mechanism is already thought to be underway Dell et al., 2019).
Similar patterns of flow acceleration and retreat linked to increases in water depth have been observed at aquatic-terminating glaciers in other glaciated regions. For example, Brown et al. (1982) and Pelto and Warren (1991) found a strong relationship between an increase in water depth and an increase in calving speed for several Alaskan tidewater glaciers. At Mendenhall Glacier, Alaska, Motyka et al. (2003) found that as the glacier thinned and retreated into deeper water, the terminus reached floatation and destabilized. Once this occurred, the terminus began to calve at an increased rate into its proglacial lake, causing the glacier to retreat further into deeper water and initiating a positive feedback mechanism (Motyka et al., 2003). Nick et al. (2009), in their study of Helheim Glacier in Greenland, successfully modelled the rapid acceleration (to 12kmyear À1 ) and thinning (>100m in 2 years) of the glacier as it retreated down its reverse-bed slope into deeper water. The retreat of the glacier over a bedrock high and back into deeper water caused ice speed and discharge to increase, leading to further thinning, retreat and the initiation of a positive feedback mechanism (Nick et al., 2009). Finally, Sakakibara et al. (2013) attributed the observed recession and acceleration at Glacier Upsala, Patagonia to a change in the longitudinal stress exerted by the bed in response to the glacier retreating over a bedrock rise and down a reversed slope into deeper water.
Although there are differences between the processes influencing aquatic-terminating glaciers in these regions and those occurring on lake-terminating glaciers in Iceland, there are many similarities in their overall dynamics which could be the focus of future investigations at Breiðamerkurjökull. For example, the velocity time series could be extended back to 1982 to coincide with the period covered by our lake area and frontal position change data. The temporal gap between our velocity measurements could also be decreased (depending on data availability), which would allow us to investigate how the velocity of Breiðamerkurjökull varies over annual time scales. Finally, we could also extend the period covered by our surface elevation change analyses to match that of our other data series. This would allow for a more detailed comparison to be made between glacier dynamics, bedrock topography, surface thinning and lake depth. As a result, a more thorough understanding of the mechanisms and interactions occurring at Breiðamerkurjökull may help us to predict the response of these larger, more dynamic aquatic-terminating glaciers to climate change (e.g. Joughin et al., 2012;Enderlin et al., 2014).

Conclusions
We have presented the first multi-year, highly detailed velocity analyses of Breiðamerkurjökull, Iceland, covering almost three decades, using high-resolution SAR offset-tracking techniques. Our data demonstrate that the greatest retreat and increase in surface velocities is found for the larger, more dynamic eastern arm, with~3.50km of recession (92.69myear À1 ) since 1982 and velocity increases occurring from~2008 up to 2015, with values of up to 3.50 ±0.25mday À1 recorded in the near-terminus region in 2015. This is in stark comparison to the other three arms, which see considerably less retreat (between~1 and 2km), and where we do not detect any change in velocity beyond the uncertainty of our methods. We also observe a substantial increase in the area of Jökulsárlón since 1982 of~20km 2 , equating to an annual growth of 0.55 km 2 year À1 . We suggest that the differences between the eastern and other arms are due to the growth and, importantly, depth increase of Jökulsárlón as the eastern arm has retreated into the deep reverse-sloping subglacial trough it carved during the LIA. This triggered an increase in ice flow at the terminus due to the inverse relationship that exists between effective pressure and basal drag, leading to an increase in terminus retreat. Retreat down this reverse slope caused the water depth to increase further, initiating a positive feedback mechanism. As our results indicate, such feedback mechanisms are likely to be important for the future response of the other outlet glaciers in Iceland. As warming continues, these glaciers will retreat into their own subglacial troughs while lake development at their termini increases in situ, causing them to undergo a similar dynamic response to what we observe for Breiðamerkurjökull in this study. Our findings are also in agreement with the results of previous studies in other glaciated regions, which have investigated the link between flow speed, water depth and retreat, and therefore may give an indication of how the larger, soft-bedded, reverse-sloped, aquatic-terminating glaciers in other glaciated regions may respond in future.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1. SNAP-obtained velocity outputs for all years covering our study period. All outputs shown at 20% transparency overlaid on a hillshade of the LiDAR DEM of Iceland from 2016. Longitudinal flowlines are also shown. Glacier outlines (white) and lake areas (light green) are from each respective year. Figure S2. Velocity profiles for the eastern arm of Breiðamerkurjökull for the period 1991-2018, calculated by taking the mean velocity at 1 km segments along the entire length of the central flowline. Associated uncertainty margins for each year are also shown. Figure S3. Velocity profiles for the central arm of Breiðamerkurjökull for the period 1991-2018, calculated by taking the mean velocity at 1 km segments along the entire length of the central flowline. Associated uncertainty margins for each year are also shown. Figure S4. Velocity profiles for the western 'A' arm of Breiðamerkurjökull for the period 1991-2018, calculated by taking the mean velocity at 1 km segments along the entire length of the central flowline. Associated uncertainty margins for each year are also shown. Figure S5. Velocity profiles for the western 'B' arm of Breiðamerkurjökull for the period 1991-2018, calculated by taking the mean velocity at 1 km segments along the entire length of the central flowline. Associated uncertainty margins for each year are also shown.