Introduction
Ice dynamics contribute substantially to accelerated global ice mass loss (King and others, Reference King2020; Choi and others, Reference Choi, Morlighem, Rignot and Wood2021; Diener and others, Reference Diener2021), yet the basal processes responsible for fast ice flow remain a significant source of uncertainty in ice-sheet models (Pattyn and Morlighem, Reference Pattyn and Morlighem2020). One potent example of transient instabilities in glacier flow are glacier surges (Meier and Post, Reference Meier and Post1969; Herreid and Truffer, Reference Herreid and Truffer2016). 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, Reference Meier and Post1969; Clarke and others, Reference Clarke, Schmok and Ommanney1986; Cogley and others, Reference Cogley2011; Jiskoot, Reference Jiskoot2011; Truffer and others, Reference Truffer2021). 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, Reference Sevestre and Benn2015). 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, Reference Dolgushin and Osipova1978; Kamb and Engelhardt, Reference Kamb and Engelhardt1987; Abe and Furuya, Reference Abe and Furuya2015; Herreid and Truffer, Reference Herreid and Truffer2016; Bhambri and others, Reference Bhambri, Hewitt, Kawishwar and Pratap2017; Van Wychen and others, Reference Van Wychen2020; Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022). These smaller-scale speedups are not unique to surging glaciers (Anderson and others, Reference Anderson2004; Vijay and others, Reference Vijay2019; Benn and others, Reference Benn, Hewitt and Luckman2023; Nanni and others, Reference Nanni2023), 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, Reference Thøgersen, Gilbert, Schuler and Malthe-Sørenssen2019). Some speedups propagate up-glacier (Anderson and others, Reference Anderson2004; Bartholomew and others, Reference Bartholomew2010; Sevestre and others, Reference Sevestre2018; Nanni and others, Reference Nanni2023) while others spread up- and down-glacier from a nucleus (Dolgushin and Osipova, Reference Dolgushin and Osipova1978; Kamb and Engelhardt, Reference Kamb and Engelhardt1987; Pritchard and others, Reference Pritchard, Murray, Luckman, Strozzi and Barr2005; Nanni and others, Reference Nanni2023). Upward-propagating speedups have been ascribed to up-glacier migration of the area affected by surface melting and crevassing (Kessler and Anderson, Reference Kessler and Anderson2004; Sevestre and others, Reference Sevestre2018; Nanni and others, Reference Nanni2023). 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, Reference Kamb and Engelhardt1987; Nanni and others, Reference Nanni2023) or migration of a thermal transition at the bed (Murray and others, Reference Murray2000; Frappé and Clarke, Reference Frappé and Clarke2007).
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, Reference Jiskoot, Pedersen and Murray2001). 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, Reference Wild2022). 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, Reference Paul2022; Benn and others, Reference Benn, Hewitt and Luckman2023) 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021), 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−1 of annual snowfall (Marcus and Ragle, Reference Marcus and Ragle1970; Newman and others, Reference Newman, Clark, Wood and Arnold2020) 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, Reference Meier and Post1969; Clarke and others, Reference Clarke, Schmok and Ommanney1986), 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, Reference Sevestre and Benn2015). 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, Reference Eisen, Harrison and Raymond2001; Bevington and Copland, Reference Bevington and Copland2014; Abe and others, Reference Abe, Furuya and Sakakibara2016; Kochtitzky and others, Reference Kochtitzky2019; Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021).
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, Reference Thornton1997), 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, Reference Goff, Lawson, Willems, Davis and Gulick2012; Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021). 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021), 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−1 (Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021). The dynamic balance line separating the reservoir and receiving zone was estimated to be ~9 km from the 2010 terminus by Nolan and others (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021), 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 km2. We study the most recent surge cycle, which began in 2013 (Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021).
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, Reference Noh and Howat2017). The 2021 DEM was processed using the NASA Ames Stereo Pipeline (ASP), which has random vertical uncertainties of ±3 m (Shean and others, Reference Shean2016). 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, Reference Larsen2010). 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, Reference Noh and Howat2017). 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, Reference Gardner2018; Lei and others, Reference Lei, Gardner and Agram2021), 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, Reference Lei, Gardner and Agram2021). 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, Reference Gardner, Fahnestock and Scambos2022). 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021). The date separations of 5–60 days between images were chosen based on successful results from Nolan and others (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021) 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, Reference Hogenson2020) 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, Reference Consortium2017). 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−1 and never exceeded ±4 m d−1 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, Reference Raymond1987; Benn and others, Reference Benn, Kristensen and Gulley2009). On Sít’ Kusá, the steep icefall in the northern tributary is thought to be an important kinematic barrier (Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021), 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−1. 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, Reference Jiskoot2011). 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021). 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021) 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−1 calculated by Bartholomaus and others (Reference Bartholomaus2015) 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, Reference Dolgushin and Osipova1978). 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021), 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021) 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 km2. 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−1 at 28 km and 5 to 8 m d−1 at 24 km, closer to the icefall. The receiving zone (<6 km, green in Fig. 1) velocities remained steadily low (<3 m d−1) 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−1 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−1. 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−1 (~82 m d−1).
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−1 (Fig. 3b) while the terminus area had peak speeds of 27 m d−1 in 2021 (Fig. 3d). By the end of summer 2021, the glacier terminus increased in area by ~4 km2 (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−1 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, Reference Jiskoot, Fox and Van Wychen2017), and at Sít’ Kusá, this produced speeds that were consistently 1–2 m d−1 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−1) 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−1. 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−1 and up to a maximum of 50 m a−1 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 (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021).
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, Reference Felikson2017). Icefalls that kinematically shield upper regions of glaciers from surges have been documented on Medvezhiy Glacier and Variegated Glacier (Raymond, Reference Raymond1987; Raymond and Harrison, Reference Raymond and Harrison1988). At Medvezhiy Glacier, the surges also initiate from the base of the icefall that feeds its reservoir zone (Raymond, Reference Raymond1987). The arrest of speedup at icefalls has been documented on several non-surging glaciers as well, in South Central Alaska (Armstrong and others, Reference Armstrong, Anderson and Fahnestock2017). 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, Reference Dolgushin, Osipova and Stulov1974; Dolgushin and Osipova, Reference Dolgushin and Osipova1978). 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, Reference Benn, Fowler, Hewitt and Sevestre2019; Terleth and others, Reference Terleth, Van Pelt, Pohjola and Pettersson2021). The increase in friction can be associated with spatial differences in glacier geometry and bed material (Fowler and others, Reference Fowler, Murray and Ng2001; Flowers and others, Reference Flowers, Roux, Pimentel and Schoof2011; Abe and others, Reference Abe, Furuya and Sakakibara2016; Crompton and others, Reference Crompton, Flowers and Stead2018; Lovell and others, Reference Lovell, Carr and Stokes2018), changes to subglacial water storage and drainage efficiency (Kamb and others, Reference Kamb1985; Kamb and Engelhardt, Reference Kamb and Engelhardt1987) and basal thermal regime (for polythermal surging glaciers) (Murray and Porter, Reference Murray and Porter2001). 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, Reference Nanni2023) including at Medvezhiy Glacier (Dolgushin and Osipova, Reference Dolgushin and Osipova1978), 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, Reference Vijay2019; Solgaard and others, Reference Solgaard, Rapp, Noël and Hvidberg2022). 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−1 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, Reference Nanni2023). 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, Reference Iken and Truffer1997; Fountain and Walder, Reference Fountain and Walder1998; Bartholomaus and others, Reference Bartholomaus, Anderson and Anderson2011). 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, Reference Benn, Fowler, Hewitt and Sevestre2019, Reference Benn, Hewitt and Luckman2023). 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 (Reference Lingle and Fatland2003). Whether the water is stored in englacial or subglacial openings such as basal crevasses (Zhan, Reference Zhan2019) 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, Reference Fountain and Walder1998). 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 (Reference Moon2014) and Vijay and others (Reference Vijay2019) and cluster 2 velocity patterns in Solgaard and others (Reference Solgaard, Rapp, Noël and Hvidberg2022), 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 (Reference Burgess, Larsen and Forster2013), 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−1 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−1) than that of a quiescent fall (<1 m d−1). 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, Reference Terlethin 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, Reference Terlethin 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, Reference Benn, Fowler, Hewitt and Sevestre2019, Reference Benn, Hewitt and Luckman2023). 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, Reference Eisen, Harrison and Raymond2001, Reference Eisen2005). 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 (Reference Eisen, Harrison and Raymond2001) 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 (Reference Eisen, Harrison and Raymond2001) 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, Reference Kamb and Engelhardt1987) as well as the annual down-glacier propagating ‘wavy surges’ at Medvezhiy Glacier (Dolgushin and Osipova, Reference Dolgushin and Osipova1978). 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−1.
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á.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2023.99.
Data
The codes used to produce georeferenced velocity maps and process the glacier velocity and elevation datasets are available through a public GitHub repository: https://doi.org/10.5281/zenodo.8361508. All velocity maps, terminus traces and DEMs produced are available through the Cryosphere Geophysics and Remote Sensing (CryoGARS) data repository: https://scholarworks.boisestate.edu/cryogars_data/3/.
Acknowledgements
This research is funded by NASA FINESST Award 80NSSC21K1640 and NSF Award ANS1954006. WorldView DEMs were provided by the Polar Geospatial Center under NSF-OPP awards 1032681, 144691 and 1542736 and created from Maxar imagery. David Finnegan from the US Army Corps of Engineers Cold Regions Research and Engineering Lab (CRREL) shared the Haenke Island weather station data. The constructive feedback from Associate Chief Editor Hester Jiskoot, Scientific Editor Rakesh Bhambri, Doug Benn and a second, anonymous reviewer improved the manuscript. We are extremely grateful to all the people that have contributed to fieldwork for the project, particularly the local businesses in Yakutat, Alaska, such as Coastal Air Service (formerly Yakutat Coastal Airlines), Yakutat Hardware, Glacier Bear Lodge, Monti Bay Lodge, Slacktide Café and the Yakutat Charter Boat Company.
Open access
