Propagating speedups during quiescence escalate to the 2020–2021 surge of Sít’ Kusá, southeast Alaska

Abstract

We use satellite image processing techniques to measure surface elevation and velocity changes on a temperate surging glacier, Sít’ Kusá, throughout its entire 2013–2021 surge cycle. We present detailed records of its dynamic changes during quiescence (2013–2019) and its surge progression (2020–2021). Throughout quiescence, we observe order-of-magnitude speedups that propagate down-glacier seasonally from the glacier's upper northern tributary, above a steep icefall, into the reservoir zone for the surging portion of the glacier. The speedups initiate in fall and gradually accelerate through winter until they peak in late spring, ~1 − 2 months after the onset of melt. Propagation distance of the speedups controls the distribution of mass accumulation in the reservoir zone prior to the surge. Furthermore, the intensity and propagation distance of each year's speedup is correlated with the positive degree day sum from the preceding melt season, suggesting that winter melt storage drives the seasonal speedups. We demonstrate that the speedups are kinematically similar to the 2020–2021 surge, differing mainly in that the surge propagates past the dynamic balance line at the lower limit of the reservoir zone, likely triggered by the exceedance of a tipping point in mass accumulation and basal enthalpy in the reservoir zone.

Introduction

Ice dynamics contribute substantially to accelerated global ice mass loss (King and others, 2020; Choi and others, 2021; Diener and others, 2021), yet the basal processes responsible for fast ice flow remain a significant source of uncertainty in ice-sheet models (Pattyn and Morlighem, 2020). One potent example of transient instabilities in glacier flow are glacier surges (Meier and Post, 1969; Herreid and Truffer, 2016). These cyclical, order-of-magnitude speedups have previously been characterized in a framework with two distinct phases: (1) a quiescent phase comprised of slow flow and localized glacier thickening and (2) a surge phase during which 10–100 times accelerated flow over a period of typically a year or more produces glacier advance and mass redistribution over ~10s of kilometers (Meier and Post, 1969; Clarke and others, 1986; Cogley and others, 2011; Jiskoot, 2011; Truffer and others, 2021). Only around 1% of the world's glaciers are confirmed to be surge-type and are concentrated in but not limited to the following regions: the Pamirs, Arctic Canada, Svalbard, Alaska-Yukon and the Karakoram (Sevestre and Benn, 2015). Although direct observations of surges are relatively rare, a growing number of observations of surging glaciers have revealed that they often undergo speedups on timescales ranging from hours to months (e.g. minisurges, wavy surges, pulses) outside of the defined surge phase (Dolgushin and Osipova, 1978; Kamb and Engelhardt, 1987; Abe and Furuya, 2015; Herreid and Truffer, 2016; Bhambri and others, 2017; Van Wychen and others, 2020; Beaud and others, 2022). These smaller-scale speedups are not unique to surging glaciers (Anderson and others, 2004; Vijay and others, 2019; Benn and others, 2023; Nanni and others, 2023), indicating that surges may be one end-member of a spectrum of ice flow variability.

Ice flow variability can manifest in several different manners. Speedups typically originate locally but propagate to other parts of the glacier (Thøgersen and others, 2019). Some speedups propagate up-glacier (Anderson and others, 2004; Bartholomew and others, 2010; Sevestre and others, 2018; Nanni and others, 2023) while others spread up- and down-glacier from a nucleus (Dolgushin and Osipova, 1978; Kamb and Engelhardt, 1987; Pritchard and others, 2005; Nanni and others, 2023). Upward-propagating speedups have been ascribed to up-glacier migration of the area affected by surface melting and crevassing (Kessler and Anderson, 2004; Sevestre and others, 2018; Nanni and others, 2023). Meanwhile, down-glacier propagating speedups have been associated with down-glacier propagating pulses of high water pressure and changes in subglacial drainage efficiency (Kamb and Engelhardt, 1987; Nanni and others, 2023) or migration of a thermal transition at the bed (Murray and others, 2000; Frappé and Clarke, 2007).

The mechanisms that drive these variations in speed are still poorly understood, in part due to the difficulty of monitoring spatiotemporal changes at the glacier bed. How do the mechanisms differ for speedups occurring on different temporal and spatial scales? Satellite observations are one approach to investigate these speedups. Remote sensing can be used to detect changes in supraglacial lakes and thus infer subglacial changes (Jiskoot and others, 2001). While satellite observations do not allow for direct observation of changes in subglacial conditions, they provide the ability to broadly map glacier behavior and provide insight to controls on the spatial evolution of flow instabilities. For example, satellite image feature-tracking techniques were used to map changes in flow velocity on Thwaites Glacier that were linked to shifts in a pinning point of the glacier's eastern ice shelf (Wild and others, 2022). Advances in the quality and quantity of satellite-derived glacier surface velocity data have increased the number of observations of speedup behavior (Paul and others, 2022; Benn and others, 2023) and may be harnessed to better understand the processes that drive their spatiotemporal evolution.

Here, we observe unstable flow on Sít’ Kusá (Turner Glacier), southeast Alaska, throughout its 2013–2021 surge cycle. The launch of the Landsat 8, Sentinel-1 and -2, and WorldView-3 satellites since the end of the 2013 surge provided the opportunity to construct extremely detailed records of surface velocity and elevation change throughout the surge cycle using pixel-tracking and stereophotogrammetry techniques. Our work builds upon the prior surge history compiled by Nolan and others (2021), which indicated that the surges of Sít’ Kusá are hydrologically driven. We observe order-of-magnitude seasonal accelerations on Sít’ Kusá that propagate down-glacier, even during the 2020–2021 surge. We investigate the spatial and temporal characteristics of the seasonal speedups, link them to Sít’ Kusá's geometry and hydrology, and discuss their similarities with the surge.

Study area

Sít’ Kusá is located on the southern side of the St. Elias mountain range in southeast Alaska, in a maritime climate characterized by ~2–6 m a⁻¹ of annual snowfall (Marcus and Ragle, 1970; Newman and others, 2020) and sea level air temperatures that range from annual extremes of −12 to 26 °C from 2014–2021 (calculated using data from the Haenke Island weather station shown in Fig. 1). Surge-type glaciers are common in the St. Elias range (Meier and Post, 1969; Clarke and others, 1986), while the greater Alaska-Yukon region contains 322 confirmed or probable surge-type glaciers, the fourth largest cluster of surging glaciers in the world (Sevestre and Benn, 2015). High accumulation and substantial summer melt promote rapid mass turnover and change in subglacial drainage, likely resulting in the abundance of surging glaciers and short surge recurrence intervals in this region (Eisen and others, 2001; Bevington and Copland, 2014; Abe and others, 2016; Kochtitzky and others, 2019; Nolan and others, 2021).

Figure 1. Map of Sít’ Kusá with 250 m surface elevation contours from the 2012 IfSAR-Alaska digital elevation model (DEM) overlain on a reduced-contrast Sentinel-2 near infrared (NIR) image from 2020-09-14. Glacier centerlines used for extracting surface elevation and velocity profiles are displayed with 2 km markings. Spatial zones defined for the 2020–2021 surge and colored points labeled by kilometer-distance from the terminus are used for the surface speed time series in Figure 3.

The Tlingit name Sít’ Kusá means ‘narrow glacier’ (Thornton, 1997), an apt name considering that the glacier is 30 km long but less than 2 km wide along its main trunk. Sít’ is the Tlingit word that means ‘glacier’ while Kusá means ‘narrow’. Sít’ Kusá is located between Sít’ Tlein (Malaspina Glacier) and Sít’ Tlein (Hubbard Glacier), both named ‘large glacier’. Sít’ Kusá can advance into Disenchantment Bay during surges but typically terminates on a pro-glacial sediment shoal when retracted in between surges (Goff and others, 2012; Nolan and others, 2021). The main trunk of the glacier (0–15 km from the terminus) rises up to 900 m elevation above sea level (a.s.l.). Two large tributaries that flow into the main trunk comprise the glacier's accumulation area, which reaches up to 2750 m a.s.l. The northern tributary contains a steep icefall ~20.5–22.5 km from the terminus, where elevation rises by ~500 m over a 2 km distance (Fig. 1).

The glacier's surge history from 1983 to 2013, compiled by Nolan and others (2021), indicates that the surge recurrence interval is relatively consistent at ~8 years. The quiescent periods typically last ~6 years. During the ~2-year surges, mass is transferred down-glacier from the reservoir zone to the receiving zone and speeds increase from ~1 to 25 m d⁻¹ (Nolan and others, 2021). The dynamic balance line separating the reservoir and receiving zone was estimated to be ~9 km from the 2010 terminus by Nolan and others (2021), based on temporally sparse surface elevation change patterns across four DEMs spanning the two previous surge cycles. During surges, Sít’ Kusá's terminus advances and increases the terminal area by 2–5 km². We study the most recent surge cycle, which began in 2013 (Nolan and others, 2021).

Data and methods

We analyzed Sít’ Kusá's kinematic evolution throughout the 2013–2021 surge cycle using detailed satellite observations of surface elevation, surface velocity and terminus area produced using the methods described in the following sections. Air temperature data from a nearby weather station (location in Fig. 1), which were used as a proxy for melt, are described below as well.

Surface elevations

We measured Sít’ Kusá's mass movement throughout the surge cycle using remotely sensed elevation data from the sources listed in (Table 1). We assumed no change in bed elevation throughout this time period and thus use surface elevation change as a proxy for ice thickness change. The majority of the elevation data are ~2-m resolution DEMs produced from Maxar WorldView-1, −2 and −3 stereo image pairs. Many of these DEMs do not cover the full extent of the glacier, as described in Table 1. Two partial DEMs in January 2018 were mosaicked, with average elevations used in overlapping regions. The 2013–2019 DEMs provided by the Polar Geospatial Center were produced using the SETSM algorithm and have random vertical uncertainties of ±2 m (Noh and Howat, 2017). The 2021 DEM was processed using the NASA Ames Stereo Pipeline (ASP), which has random vertical uncertainties of ±3 m (Shean and others, 2016). The Interferometric Synthetic Aperture Radar – Alaska (IfSAR-Alaska) DEM from 2012-08-26 has a vertical accuracy within 3 m for slopes of 0–10 degrees (90% confidence). The 2012 DEM, although outside of the 2013–2021 surge cycle, was the most spatially continuous DEM available for Sít’ Kusá and was used solely as an input for georeferencing the velocity maps. The centerline elevation dataset from June 2020 was collected by IceBridge UAF Lidar scanner L1B Geolocated Surface Elevation Triplets, Version 1 (ILAKS1B) which has a vertical accuracy of ±0.3 m (Larsen, 2010). To compare with the other DEMs, we generated a 2 m resolution DEM of the centerline elevations using the SAGA Rasterize function in QGIS, which averages the lidar point elevations in each 2 m pixel. We did not coregister the DEMs, because they are primarily standard WorldView products which are precisely coregistered via SETSM (Noh and Howat, 2017). We estimate that horizontal offsets of several pixels over a gradually sloped glacier like Sít’ Kusá would result in several meters of vertical offset at most, which is much smaller than magnitude of elevation changes measured throughout the surge (10s of meters). The elevation change maps indicate no systematic biases between the DEMs (Fig. S1).

Table 1. Summary information for glacier surface elevation datasets used for analysis of mass movement throughout the 2013–2021 surge cycle

We extracted surface elevation profiles at 500 m increments along manually delineated centerlines for the trunk and two major tributaries. At each centerline point, we averaged elevations within a 250 m radius circle to reduce uncertainties associated with the datasets and potential biases introduced by crevassing. From these elevation profiles, we used the slope of the line of best fit within each of the spatial zones shown in Figure 1 to calculate mean surface slope for that zone.

Surface velocities

We used data from the automated pixel-tracking and speckle-tracking algorithm, autoRIFT (Gardner and others, 2018; Lei and others, 2021), applied to satellite image pairs from Landsat 8, Sentinel-1 and Sentinel-2 to measure surface velocities from June 2013 (after the Landsat 8 launch) through mid-2022. The autoRIFT algorithm applies a normalized image cross-correlation with iteratively larger chip sizes until pixel displacements are successfully calculated (Lei and others, 2021). For each image pair, autoRIFT uses the pixel displacements and the date separation of the two images to calculate the time-averaged surface velocity.

The global ice velocities provided by the NASA MEaSUREs ITS_LIVE project are produced using autoRIFT with standard parameters for mapping ice motion at the regional scale (Gardner and others, 2022). While the ITS_LIVE data successfully characterize Sít’ Kusá's quiescent motion, there are gaps during the surge when the glacier's main trunk (0–15 km from the terminus) accelerates (Fig. S2). To fill the gaps in the ITS_LIVE speed time series during the surge, we applied autoRIFT with custom parameters to Sentinel-2A, Sentinel-2B and Landsat 8 images that were 90% cloud-free over the glacier area. To optimize feature tracking during the surge, we ran autoRIFT using custom pixel-tracking parameters that include minimum chip size, maximum search range and a reference velocity map from the 2006 surge from Nolan and others (2021). The date separations of 5–60 days between images were chosen based on successful results from Nolan and others (2021) and adjusted based on the quality of the resulting velocity maps. We kept all ITS_LIVE velocities that were within the same date separation range (5–60 days) from 2013 to 2021.

We produced georeferenced velocity maps from all pairs of 10 m resolution NIR images from the Sentinel-2 satellites and 15 m resolution panchromatic images from Landsat 8. The 2012-08-26 IfSAR DEM was resampled to the target chip size and used for georeferencing. For each image pair, we ran autoRIFT with minimum chip sizes ranging from 100 to 300 m. The maximum search range was automatically calculated using the reference velocity map. For each date separation, we then chose the resulting map that contains the greatest fraction of data within the glacier outline and the least error over stable surfaces (Table S1). Over autumn and winter months, when clouds and darkness hinder feature-tracking with optical images, we used the OnDemand autoRIFT service available through the Alaska Satellite Facility (ASF) HyP3 pipeline (Hogenson and others, 2020) to produce velocity maps from Sentinel-1A and -1B image pairs with temporal separation of exactly 12 days. We constructed speed time series using the ITS_LIVE, ASF and custom autoRIFT velocity datasets throughout the entire surge cycle.

Velocity uncertainties

For the standard velocity datasets, ITS_LIVE and ASF autoRIFT, we used the reported errors. At each geographic coordinate, the ITS_LIVE velocity dataset reports speed errors associated with each data point in the time series. For each ASF autoRIFT velocity map, we used the per-pixel speed uncertainties provided with the product. For the custom autoRIFT velocity maps, we calculated the error over stable surfaces. We defined the stable surface error for each velocity map as the root-mean-square error (RMSE) in surface speed over non-ice pixels, which should have no substantial motion over the spatial ranges and temporal separations used for velocity-mapping. The stable surface mask comprised all areas outside of a 200 m buffer around all glacier outlines available through the Global Land Ice Measurements from Space (GLIMS) Randolph Glacier Inventory (RGI, 2017). This mask was manually edited to remove additional area around the termini of Sít’ Kusá and Sít’ Tlein (Hubbard Glacier) in order to account for possible terminus advance since the date that the GLIMS outlines were created. The mean stable surface error for all custom velocity maps was <1 m d⁻¹ and never exceeded ±4 m d⁻¹ for any individual map (Table S1).

Speedup extent

Surges propagate over several km in distance, comprising large portions of a glacier or even an entire glacier (Raymond, 1987; Benn and others, 2009). On Sít’ Kusá, the steep icefall in the northern tributary is thought to be an important kinematic barrier (Nolan and others, 2021), as confirmed by the observations presented below. Thus, we use the icefall as the reference point for measuring the down-glacier propagation distance of all observed speedups, including the 2020–2021 surge. We define the speedup propagation distance in the velocity maps as the centerline distance from the icefall to where the surface speed remains below 1 m d⁻¹. Following the principle of mass conservation, the leading edge of a propagating speed perturbation is also associated with a perturbation in mass (e.g. surge bulge) (Jiskoot, 2011). Therefore, we also inferred speedup propagation distance from surface elevation changes and compared them with the maximum speedup propagation distance measured from the velocity maps. The down-glacier propagation distance of all speedups measured using the two datasets agree within 500 m.

Terminus area

We measured the change in Sít’ Kusá's terminus area by manually tracing the entire terminus lobe using satellite images, starting from a reference line defined by two tie points on the lateral margins of the glacier (Fig. S3b) as done by Nolan and others (2021). All Collection 2 Level-1 Landsat 8 Operational Land Imager (OLI) panchromatic (B8) images from 2013 to 2022 that covered the full terminus without cloud obstruction were used to trace terminus areas. There are 10–20 terminus area delineations per year. Seventy-one of the terminus traces in this study were coincident with traces made by Nolan and others (2021) from 2013 to 2016 (Fig. S3c). We used the RMSE in terminus area calculated from overlapping traces (i.e. the inter-analyst uncertainty) to estimate the error in the terminus area measurements.

Air temperatures

We used air temperatures above freezing as a proxy for surface melt on the glacier, focusing on the timing of meltwater input and relative differences in total annual melt. We calculated the average air temperature from two independent temperature sensors on a weather station maintained by the US Army Corps of Engineers Cold Regions Research and Engineering Laboratory (CRREL), located on Haenke Island near the terminus of Sít’ Tlein (Hubbard Glacier) (Fig. 1), situated approximately at sea level. The data were available from May 2014 through September 2022 at sub-hourly resolution. In order to estimate when temperatures rise above freezing at different locations on the glacier, we corrected the air temperatures by elevation at those locations using the mean atmospheric lapse rate of −7.1 $^\circ$C km⁻¹ calculated by Bartholomaus and others (2015) in the adjacent Icy Bay. We count a positive degree day (PDD) when the elevation-corrected daily average temperature is above 0. There were no gaps in the daily average temperature record.

Results

On a surging glacier, the reservoir and receiving zones are characterized by the following patterns in thickness change: the reservoir zone is the area on the glacier that thickens during quiescence and thins during the surge while the receiving zone is the region that thins during quiescence and thickens during the surge (Dolgushin and Osipova, 1978). Thus, for Sít’ Kusá's 2020 surge, the reservoir and receiving zones in Figure 1 are defined using the surface elevation time series (Fig. 2). The dynamic balance line separates these two zones, approximately 8 km from the glacier terminus. The steep icefall in the northern tributary at ~22 km, identified as a likely kinematic barrier by Nolan and others (2021), separates the upper northern tributary from the reservoir and receiving zones (Fig. 1).

Figure 2. Surface elevation profiles along Sít’ Kusá's northern centerline below the icefall (a,c) and southern centerline (b) during quiescence (a,b) and the surge (c). Panel d shows surface slope in degrees within the spatial zones outlined in Figure 1, calculated along the northern glacier centerline. Panel e shows the elevation deviation from the 2013-12-07 surface along the northern centerline, including the icefall (yellow box) and area up-glacier of it. Dotted vertical lines in panels a, b, c and e indicate where the tributaries converge. Black arrows mark the location of the dynamic balance line, at ~8 km.

Evolution during quiescence, 2013–2019

The 2012 surge described in Nolan and others (2021) terminated by the end of the 2013 summer, confirmed by the velocity data we present here (Video S1, Fig. 3). Thus, the DEM from December 2013 captures the post-surge geometry, acquired during the beginning of a quiescent period. The following DEMs from 2014 to 2019 resolve thickening up-glacier in the reservoir zone and thinning in the receiving zone down-glacier (Figs 2a, b; Fig. S1), producing a 1.25x increase in the glacier surface slope during the quiescent period. The centerline elevation profiles indicate that the quiescent thickening was concentrated in the northern tributary 16–18 km from the terminus and down-glacier of the confluence area of the two tributaries (Fig. 2c). Between 16 and 18 km, the glacier thickened by up to 120 m (Figs 2a, e). Meanwhile, the terminal lobe thinned by ~60 m and shrunk by ~3 km². Above the icefall, the maximum observed elevation change was approximately +20 m throughout the entire surge cycle (Fig. 2e) and the surface slope remained virtually constant in this region (Fig. 2d).

Figure 3. Glacier surface speed from 2013 to 2022 at points along the centerlines, labeled by km-distance from the terminus and grouped by spatial zones (Fig. 1). Surface speeds are plotted as rectangular patches with widths corresponding to the time period covered by the velocity maps and heights corresponding to the stable surface errors. Vertical grid lines mark 1 January of each year. The gray box bounds the 2020–2021 surge and the shaded portion corresponds to the period of surge front propagation to the terminus. Arrows in panels b and c emphasize multiyear patterns in speedup. Inset in panel b is an idealized diagram of the down-glacier propagation of seasonal speedups in the reservoir zone. Panel d shows the terminus area (blue, secondary y-axis) along with receiving zone speeds.

The surface velocity time series resolved seasonal increases in speed that appear to propagate down-glacier from above the icefall in the northern tributary throughout the quiescent period, which we refer to as seasonal speedups (Video S1). Surface velocities in the upper northern tributary (>22 km, orange in Fig. 1) remained steady throughout the surge cycle, ranging from 1 to 4 m d⁻¹ at 28 km and 5 to 8 m d⁻¹ at 24 km, closer to the icefall. The receiving zone (<6 km, green in Fig. 1) velocities remained steadily low (<3 m d⁻¹) throughout quiescence. Meanwhile, in the reservoir zone (6–20 km, purple in Fig. 1), we observe speedups that propagated down into the reservoir zone, reaching peak speed (Fig. 3) and extent (Fig. 4) in late spring. The timing of the seasonal speedups is relatively consistent from year to year. Each year, the speedup initiated around September/October, continued to increase gradually through winter and spring until it reached peak speed around May/June, then quickly returned to slow speeds of ~1 m d⁻¹ by July (Fig. 3b). The peak speed attained during seasonal speedups (as measured at 18 km) varied from year to year, ranging from 2 to 10x the slowest annual speeds (Figs 3b–d). In some years, speedups with similar peak timing were resolved in the southern tributary as well (Fig. 3c), although they were weaker than those in the northern tributary and do not appear to reach the main trunk (Video S1). From 2013 to 2017, the speedups propagated further down-glacier until the 2017 speedup is arrested at the dynamic balance line at ~7–8 km (Fig. 4a). Interestingly, the peak speed and propagation distance dropped in 2018, and then increased monotonically again until the surge (Figs 3b, 4a).

Figure 4. (a) Surface speed evolution along the northern centerline of Sít’ Kusá. Black indicates no data. (b–i) Surface speed maps from 2014 to 2021 showing the peak glacier speeds each year with 2 km markings and the same color scale as in (a). The quiescent speed maps (b–g) show the maximum propagation distance of the seasonal speedups, also marked by the short white vertical bars in (a). Vertical dotted white line denotes the point of convergence between the two tributaries while the dashed white line denotes the dynamic balance line. The dates of the speed maps vary and are listed at the top of each subplot.

Progression of the 2020–2021 surge

The area 16–18 km from the terminus along the northern centerline, the region of greatest thickness buildup during quiescence, showed the first signs of drawdown after the surge initiated in 2020 (Fig. 2e). Thus, the surge bulge nucleated there, just below the icefall in the northern tributary. The 2019–2020 speedup passed through the 16–18 km area and continued to accelerate through the rest of the glacier beginning in March 2020, marking the initiation of the surge. The high velocities rapidly propagated down-glacier from the reservoir zone, such that the down-glacier limit of the propagating high velocity anomaly (i.e. surge front) reached the terminus in the fall of 2020 (Fig. 4).

When the surge front neared the terminus in late spring 2020, we observed a 20-fold increase in speed at the terminus compared to quiescent velocities (Fig. 3d). Field observations indicated that the surge front had not yet reached the terminus by August of 2020 (Fig. S4), although the terminus was advancing as fast as 7 m d⁻¹. The velocity maps between July and October 2020 contain data gaps near the terminus that prevent the precise identification of the date of surge front propagation to the terminus. The surge velocities entered the terminal lobe by mid-September although there are gaps in the velocity map near the terminus boundary (Fig. S5j). The glacier covered its sediment shoal and reached the bay by mid-October 2020 according to the satellite image record, indicating that the surge front had reached the terminus before then. Thus, the surge front likely reached the terminus between September and October of 2020. The propagation of the surge front by ~18 km down-glacier over the 6- to 7-month period yields a mean down-glacier surge propagation velocity of ~30 km a⁻¹ (~82 m d⁻¹).

While the surge propagated, the down-glacier ice (e.g. at 9 km) began to move faster than the up-glacier ice (e.g. 18 km), reversing the longitudinal stress gradient present throughout the localized quiescent speedups (Fig. 3). The reservoir zone reached peak speeds in late spring of the surge (in 2020) that were 2–4x larger than peak seasonal speeds during quiescence (Fig. 3b). The speed in the southern tributary remained slow until the confluence area (e.g. 15 km) accelerated, indicating that surge speeds propagated upwards through that tributary, via down-glacier extension.

The entire glacier besides the upper northern tributary slowed over the 2020–2021 winter, but maintained speeds above the typical quiescent winter speeds everywhere below the icefall (Fig. 3). Following a 2-month gap in the velocity record, in February 2021 speeds were higher than in 2020 and continued to accelerate, reaching peak speeds several meters per day higher than in 2020. In the reservoir zone, the peak speeds ranged from 10 to ~25 m d⁻¹ (Fig. 3b) while the terminus area had peak speeds of 27 m d⁻¹ in 2021 (Fig. 3d). By the end of summer 2021, the glacier terminus increased in area by ~4 km² (Fig. 3d) and in thickness by ~90 m (Fig. 2b). Surface elevation and slope in the reservoir zone decreased to approximately the values from December 2013, after the end of the previous surge (Fig. 2). Meanwhile, the speeds slowed to less than 2 m d⁻¹ throughout the majority of the glacier, indicating the termination of the surge (Fig. 4). Afterwards, the glacier appears to resume quiescent patterns of motion in the reservoir zone, with indication of a 2022 speedup at 18 km (Fig. 3b).

Discussion

Speedup origin

The speedups that propagated down-glacier through the reservoir zone originated from the fast-moving ice above the icefall in the northern tributary (Video S1, Fig. 4). Steep icefalls longitudinally pull ice down-glacier (Jiskoot and others, 2017), and at Sít’ Kusá, this produced speeds that were consistently 1–2 m d⁻¹ faster near the icefall at 24 km compared to 28 km (Fig. 3a). The late summer snowline, an approximation of the equilibrium line altitude (ELA), observed in the satellite imagery from September 2020 extends to the base of the icefall, ~1200 m a.s.l. Thus, subglacial hydrology is not likely to play a major role in the subtle seasonal variation in speed (amplitude of 1–2 m d⁻¹) for the area above the icefall, several hundred meters above the ELA. Once the speedups propagated past the icefall, the amplitude of seasonal speed variation in the reservoir zone increases substantially, to 1–10 m d⁻¹. The seasonal down-glacier propagation of high speeds from above the icefall provided pulses of high mass flux to the reservoir zone, primarily thickening the glacier along the northern centerline at an average rate of 20 m a⁻¹ and up to a maximum of 50 m a⁻¹ between 2014 and 2020. The 2020 surge initiated at the base of the icefall as a typical winter speedup but propagated first through the entire reservoir zone, breaching the dynamic balance line (8 km from the terminus) at which previous quiescent speedups had arrested (Fig. 4). Acceleration continued as the surge propagated through the receiving zone to the glacier terminus. Drawdown occurred throughout the entire reservoir zone but did not extend above the icefall (Fig. 2). Thus, the tributary above the icefall remained kinematically separate from the surging portion of the glacier below, as suggested for previous surges by Nolan and others (2021).

The steep icefall appears to behave as a one-way kinematic barrier. It allows speedup down through the icefall that can propagate and distribute mass down-glacier, but appears to dampen up-glacier propagation of speedup and dynamic thinning, shielding the upper northern tributary from velocity and mass fluctuations associated with the surge (Felikson and others, 2017). Icefalls that kinematically shield upper regions of glaciers from surges have been documented on Medvezhiy Glacier and Variegated Glacier (Raymond, 1987; Raymond and Harrison, 1988). At Medvezhiy Glacier, the surges also initiate from the base of the icefall that feeds its reservoir zone (Raymond, 1987). The arrest of speedup at icefalls has been documented on several non-surging glaciers as well, in South Central Alaska (Armstrong and others, 2017). At Sít’ Kusá, the icefall is an important geometrical feature that regulates the rapid filling of the reservoir zone and promotes the glacier's short surge recurrence interval.

Speedup propagation

The seasonal down-glacier propagating speedups on Sít’ Kusá are most similar to those documented on Medvezhiy Glacier, described as ‘wavy surges’ that advanced the thickness anomaly down-glacier during quiescence and appeared to differ from its surges only in intensity (Dolgushin and others, 1974; Dolgushin and Osipova, 1978). Reservoir zones on surging glaciers are thought to form due to an increase in friction that slows ice velocity relative to the balance velocity, promoting thickening behind the dynamic balance line (Benn and others, 2019; Terleth and others, 2021). The increase in friction can be associated with spatial differences in glacier geometry and bed material (Fowler and others, 2001; Flowers and others, 2011; Abe and others, 2016; Crompton and others, 2018; Lovell and others, 2018), changes to subglacial water storage and drainage efficiency (Kamb and others, 1985; Kamb and Engelhardt, 1987) and basal thermal regime (for polythermal surging glaciers) (Murray and Porter, 2001). At Sít’ Kusá, there are insufficient observations of the glacier bed to determine the cause of locally high basal resistance at the dynamic balance line, but the surface elevation bump observed at that location (black arrow in Fig. 2) supports the presence of a bedrock ridge, change in bed lithology or till with anomalous shear strength.

While locally elevated basal resistance at the 7–8 km feature likely bounds the maximum propagation distance of the seasonal speedups, it does not explain the substantial drop in propagation distance in 2018. The velocity maps indicate that the 2018 speedup propagated only 6 km past the icefall (14 km from the terminus) and 8 km up-glacier of the dynamic balance line. As a result, mass piled up in the upper portion of the reservoir zone just below the icefall, 14–18 km from the terminus (Fig. 5a). If increased speedup propagation distance was controlled only by the gradual increase in thickness as the reservoir zone fills, the 2018 propagation distance would at least match the propagation distance in 2017, which arrested at the dynamic balance line. The less extensive speedup in 2018 indicates that another factor modulates speedup propagation distance.

Speedup seasonality and hydrologic control

Speedups that last over seasonal timescales and propagate on both surging and non-surging glaciers, e.g. in the Pamir Mountains (Nanni and others, 2023) including at Medvezhiy Glacier (Dolgushin and Osipova, 1978), have been attributed to seasonal changes in subglacial drainage. The comparison of timing of seasonal speedup and surface runoff can yield insight to the evolution of subglacial drainage systems on glaciers (Vijay and others, 2019; Solgaard and others, 2022). At Sít’ Kusá, we used elevation-corrected daily temperatures above freezing (i.e. PDD) recorded at the Haenke Island weather station as a proxy for surface melt delivered to the subglacial hydrologic network. We focus on interannual differences in cumulative melt and the timing of melt and discuss hydrological control in terms of relative changes in efficiency of subglacial drainage and its impact on effective pressure. We used the surface speed time series at 18 km from the terminus (see location in Fig. 1), where the speedups first propagate through, to evaluate the seasonal evolution in glacier speed with respect to annual cumulative PDD (Fig. 6). The surface speeds at 18 km were averaged weekly and gaps were filled using a simple linear interpolation in order to execute cross-correlation analysis of the variables.

Figure 5. (a) Surface elevation deviation along the northern centerline of Sít' Kusá since 2013-12-07, the end of the previous surge. Yellow box outlines the icefall area. Triangles mark the extent of thickening in the reservoir zone, an estimate of the extent of down-glacier propagating speedups that originate above the icefall. Down-glacier propagation distance is measured from the icefall (22 km from the terminus). (b) Down-glacier propagation distance of the speedups from velocity maps temporally closest to the elevation datasets. The solid black line is the 1-1 line indicating perfect agreement in estimated propagation distance from the elevation profiles and velocity maps. Black dashed lines in both panels mark the dynamic balance line.

Figure 6. (a) Weekly averaged surface speeds at 18 km from the terminus, annual cumulative PDD from air temperature corrected by annual elevations at that point and annual PDD gradient from 2014 to 2021. Annual grid lines mark the 1st day of the year. Vertical dashed lines mark the annual timing of the onset of temperatures above freezing, i.e. melt onset (black), and deceleration to the slowest fall speeds (purple). Purple horizontal bars along the top of the panel mark the duration of the speedup. (b) Temporal cross-correlogram of the speed and annual cumulative PDD. The square marks the maximum correlation, which corresponds to a time lag of –6 months.

Surface speeds at the 18 km point were slowest in mid-summer (late June to July) when there is ~1 m d⁻¹ of surface motion each year, except for the first year of the surge, which is discussed separately below. Speeds increased from their minimum values as each melt season waned, and gradually increased through fall and winter without the addition of any melt (Fig. 6a). The acceleration continued for more than 6 months. The dense crevasse fields (Fig. S4), especially in the reservoir zone below the icefall, provide a means for surface meltwater to enter internal cavities. The gradual closure of the drainage system over winter may then increase water pressure and drive speedup (Nanni and others, 2023). The reduction of meltwater input to these cavities over winter allows for creep closure and the subglacial drainage system becomes more disconnected and less efficient, as is typical of temperate glaciers (Iken and Truffer, 1997; Fountain and Walder, 1998; Bartholomaus and others, 2011). Gradual creep closure of increasingly isolated water pockets throughout winter sustains a continuous decline in bed effective pressure and in turn, ice acceleration. The increase in frictional heating associated with ice acceleration generates additional basal water, thus increasing basal water pressure and storage that likely further contributes to ice acceleration (Benn and others, 2019, 2023). Another possible mechanism proposed for winter speedup is the englacial storage of water that continues to travel to the bed throughout winter via gravity-driven flow, as proposed by Lingle and Fatland (2003). Whether the water is stored in englacial or subglacial openings such as basal crevasses (Zhan, 2019) still needs to be investigated. Here, we generalize the water storage as ‘internal’.

The speedups nearly reach their maxima prior to the melt onset that occurs between mid-April and the end of May each year (Fig. 6a). After the melt onset, surface speeds either continue to increase briefly or plateau until the June/July deceleration (Fig. 6a). The lack of a secondary acceleration upon the onset of surface melt suggests that the new meltwater does not suddenly overwhelm the subglacial drainage system present at the end of winter. In the case of continued gradual surface speed increase, the additional meltwater input likely produces gradual increase in basal water pressure until the the subglacial drainage system becomes efficient enough to flush all water out in June/July. In contrast, plateaued surface speeds suggest additional meltwater merely sustains the level of basal water pressurization that was established at the end of winter. In either case, the substantial annual snowpack (up to 6 m or more) present on Sít’ Kusá in April/May may buffer the response of the subglacial drainage network to seasonal melt onset because melt transport through the snowpack may slow entry of surface meltwater into the glacier hydrological system, allowing for a more gradual response.

The sudden deceleration just 1- to 2-months after the melt onset (Fig. 6a) implies that the subglacial drainage system adjusts rapidly to accommodate the influx of surface meltwater, evolving into a more efficient drainage system that is able to quickly flush out water and relieve the buildup in water pressure throughout the rest of the melt season (Fountain and Walder, 1998). In other words, there is a remarkably rapid evolution to a well-connected drainage network with low water pressure. The timing of these seasonal speed changes aligns with those observed on outlet glaciers in Greenland categorized as having type III velocity patterns in Moon and others (2014) and Vijay and others (2019) and cluster 2 velocity patterns in Solgaard and others (2022), more prevalent in regions with greater amounts of surface melt such as southern Greenland.

The cumulative annual melt (PDD) appears to be proportional to the peak speed of the following winter's speedup for all quiescent years and the first year of the surge (Fig. 6a). Cross-correlation between weekly speeds and PDD reveals annual peaks in correlation that are offset by –6 months (Fig. 6b), indicating that each winter's speedup is most correlated with the PDD from the previous melt season. Furthermore, the speedups that attain faster peak speeds also last longer and propagate further (Figs 6a, S6). Thus, more spring and summer melt produces a more intense speedup the following winter and spring. These observations suggest that the amount of summer melt modulates the amount or pressure of water that is stored internally over winter that in turn is responsible for the majority of the seasonal speedup. This is opposite to the inverse relationship between summer melt and winter speed that was found across 160 glacier systems in Alaska by Burgess and others (2013), and suggests that speedup mechanism may not be ubiquitous across Alaska.

Comparison of speedup seasonality to the surge seasonality

The initiation of the surge in 2020 follows a similar seasonal evolution to the quiescent speedups, continuing from a winter speedup that peaks in late spring. However, instead of the rapid deceleration to ~1 m d⁻¹ in June/July during the quiescent years, the 2020 speeds decelerate more gradually, reaching an annual low in mid-October that is far higher (~4 m d⁻¹) than that of a quiescent fall (<1 m d⁻¹). Furthermore, the speed during the surge winter remains faster than in quiescent winters (Fig. 6a). Analysis of glaciohydraulic tremor from the seismometer network deployed around Sít’ Kusá indicates that tremor levels associated with subglacial hydraulic activity remain higher through the surge winter (2020–2021) compared to the post-surge winter (2021–2022) (Terleth and others, in review). The drainage system remains more active throughout winter, possibly due to the increase in the amount of water retained after the melt season that did not fully drain in 2020. After the end of the melt season in mid-October 2020, a new winter acceleration builds on this fall minimum due to gradual closure of isolated cavities. The surge termination occurs over the course of the 2021 summer (Terleth and others, in review), far more gradually and later than the spring quiescent slowdowns.

Mechanism for surge initiation

In the context of the mass and enthalpy imbalance framework for the buildup of surge instabilities, the multiannual buildup of imbalances in mass (via quiescent thickening) as well as basal ice temperature and water (i.e. enthalpy) trigger surges (Benn and others, 2019, 2023). Increase in basal enthalpy that reduces basal friction throughout quiescence can come from many sources, including increased water storage and changes in drainage efficiency. Gradual thickening of the glacier increases the driving stress and therefore basal shear stress prior to surges. Studies of cumulative mass balance of Variegated Glacier suggest that surges initiate once a critical basal shear stress is attained (Eisen and others, 2001, 2005). However, at Sít’ Kusá, the annual order-of-magnitude speedups can produce substantial annual fluctuations in basal water pressure that alters this critical basal shear stress. The non-monotonic change in seasonal speedup amplitude on Sít’ Kusá appears to be correlated to annual meltwater input rather than a gradual buildup in water storage capacity or driving stress. Ultimately, Sít’ Kusá's 2020 surge may have been triggered by the combination of above-average annual melt in 2019 and the localized mass accumulation modulated by the seasonal speedups in the reservoir zone.

Speedup amplitude and propagation distance correlated with annual cumulative PDD. The cumulative PDD in 2016 and 2019 (directly prior to the surge) were similar; however, the 2016–2017 speedup was arrested at the dynamic balance line whereas in 2020, the surge propagated past it. The meltwater input for the 2016–2017 speedup may have provided enough enthalpy gains to drive speedup propagation all the way to the glacier terminus, except that the accumulated mass was likely not enough to sufficiently alter basal shear stress after only 3 years of mass accumulation in the reservoir zone. There was an additional ~40 m of thickening from 2017 to 2019 before the initiation of the surge (Fig. 6). Eisen and others (2001) were able to calculate cumulative mass balance in Variegated Glacier's reservoir zone using an empirical relationship between annual balance and temperature and precipitation data, tuned to measured annual balances for a 10-year period. Using the parameters determined by Eisen and others (2001) at the nearby Variegated Glacier, just across Disenchantment Bay, and climatic data from the down-scaled Weather Research and Forecasting product (Text S1), we calculate thickness gain from just surface mass balance of ~35 m during Sít’ Kusá's quiescent period, less than one-third of the total 120 m of thickening in the reservoir zone. Thus, at Sít’ Kusá, the annual mass flux into the reservoir zone from the area above the icefall far exceeds the cumulative surface mass balance directly in its reservoir zone. Thus, the role that annual melt plays on speedup propagation from the upper northern tributary likely produces variability in the surge recurrence interval, which is regulated overall by accumulated mass balance in the reservoir zone. Given these observations, it is likely that both high basal enthalpy and exceedance of a critical driving stress that triggers surges at Sít’ Kusá.

Conclusions

Satellite observations of Sít’ Kusá's surface velocities and elevations throughout the 2013–2021 surge cycle reveal that its northern tributary is extremely dynamic even during the quiescent period. The observed variability in velocity and elevation is akin to the ‘minisurges’ documented at non-surging alpine glaciers (Kamb and Engelhardt, 1987) as well as the annual down-glacier propagating ‘wavy surges’ at Medvezhiy Glacier (Dolgushin and Osipova, 1978). Throughout quiescence, 2013–2019, we observed seasonal speedups that last ~8 months through fall, winter and spring and that propagate down-glacier from the upper northern tributary. We found that the amplitude and propagation distance of the quiescent speedups and the surge are correlated with the cumulative PDD from the preceding melt season each year, suggesting that meltwater storage through winter fuels the propagating speedups. The surge propagates all the way to the glacier terminus in 2020. In winter 2020/2021, the glacier maintains elevated winter velocities compared to quiescent winters, resulting in a second pulse of acceleration that lasts until the surge terminates throughout summer 2021. The quiescent speedups appear to differ from the 2020 surge only in their magnitude and duration, raising the prospect that surging glaciers differ from their non-surging counterparts only in the extent to which they are sensitive to meltwater. On non-surging glaciers, imbalances in basal enthalpy and mass are effectively dissipated. On surging glaciers, localized accumulation of basal enthalpy and mass over multiannual timescales eventually trigger a surge but may manifest as smaller speedups in quiescent years.

Seasonal variation in velocities during the quiescent speedups and the 2020–2021 surge was likely related to seasonal evolution of the subglacial drainage system, although the timing differs slightly between the two. The quiescent speedups decelerate rapidly partway through the melt season (June/July), just 1–2 months after the melt onset. The winter subglacial drainage system configuration and water retention allow for the rapid adjustment of drainage efficiency after melt onset to flush out the additional meltwater input. During the surge, velocities decrease gradually through the end of the melt season, indicative of a less efficient subglacial drainage system, where glacier speed and effective pressure both respond to melt input. Elevated velocities during the winter of the surge indicate that the first year of the surge did not fully drain the subglacial water, which likely allows the renewed acceleration in the second year of the surge, when the glacier terminus reaches its highest velocities of 30 m d⁻¹.

The icefall in the northern tributary and a high-friction feature at the dynamic balance line (7-8 km from the glacier terminus) bound the reservoir zone of Sít’ Kusá. The active upper northern tributary above the icefall modulates high rates of filling of the reservoir zone. The quiescent speedups that propagate from that upper area control the redistribution of thickness throughout the reservoir zone. Our observations suggest that the mass accumulation in the reservoir zone modulated by the hydrologically-driven seasonal speedups trigger surges at Sít’ Kusá.