Kinematics of the exceptionally-short surge cycles of Sít’ Kusá (Turner Glacier), Alaska, from 1983 to 2013

Abstract

Glacier surges are periodic episodes of mass redistribution characterized by dramatic increases in ice flow velocity and, sometimes, terminus advance. We use optical satellite imagery to document five previously unexamined surge events of Sít’ Kusá (Turner Glacier) in the St. Elias Mountains of Alaska from 1983 to 2013. Surge events had an average recurrence interval of ~5 years, making it the shortest known regular recurrence interval in the world. Surge events appear to initiate in the winter, with speeds reaching up to ~25 m d⁻¹. The surges propagate down-glacier over ~2 years, resulting in maximum thinning of ~100 m in the reservoir zone and comparable thickening at the terminus. Collectively, the rapid recurrence interval, winter initiation and down-glacier propagation suggest Sít’ Kusá's surges are driven by periodic changes in subglacial hydrology and glacier sliding. Elevation change observations from the northern tributary show a kinematic disconnect above and below an icefall located 23 km from the terminus. We suggest the kinematic disconnect inhibits drawdown from the accumulation zone above the icefall, which leads to a steady flux of ice into the reservoir zone, and contributes to the glacier's exceptionally short recurrence interval.

1. Introduction

Glacier surges are periodic instabilities that result in increased surface velocities, mass redistribution and terminus advance (Meier and Post, 1969). Although <1% of glaciers worldwide are classified as surge-type (Jiskoot and others, 1998; Sevestre and Benn, 2015), Alaska and western Canada are home to 113 confirmed surge-type glaciers (Sevestre and Benn, 2015). There are also notably high concentration of surge-type glacier in Svalbard, East Greenland, the Pamir and Karakoram Mountains (Sevestre and Benn, 2015). Surge events have the potential to strongly influence glacier mass balance, as they transport mass from the reservoir zone, across the dynamic balance line, to the lower elevation receiving zone (Meier and Post, 1969; Dolgoushin and Osipova, 1975; Raymond, 1987), where ablation rates are higher and mass loss is accelerated (Aðalgeirsdóttir and others, 2005). Surge events have been documented at both marine-terminating (e.g. Murray and others, 2003; Sevestre and others, 2018) and land-terminating glaciers (e.g. Kamb and others, 1985), after the collapse of ice-shelves (e.g. De Angelis and Skvarca, 2003), and at outlet glaciers of Arctic ice caps (e.g. Dunse and others, 2015; Willis and others, 2018). Glacier surges are one example of a spectrum of fast-flow events (Clarke, 1987a; Herreid and Truffer, 2016) including ice streaming (e.g. Blankenship and others, 1987), the tidewater glacier cycle (e.g. Meier and Post, 1987), ice-avalanching (e.g. Gilbert and others, 2015) and glacier collapse (e.g. Kääb and others, 2018; Jacquemart and others, 2020). While each of these instabilities have unique characteristics and physics that control them, they all are marked by their stark departure from steady-state dynamics and are inextricably linked to basal processes. Glacier surges present a prime natural laboratory to study ice instabilities given the high frequency at which they occur (cf., ice streaming, tidewater glacier cycle) and the wealth of existing research into the mechanisms that control these instabilities (e.g. Kamb, 1987; Fowler and others, 2001; Jay-Allemand and others, 2011).

There are presently two main hypotheses to explain why glacier surges occur. One hypothesis is the polythermal switch; when cold ice frozen to the bed rapidly transitions to warm ice detached from the bed, triggering acceleration (Fowler and others, 2001). The polythermal switch mechanism is commonly proposed for surge-type glaciers in Svalbard (Murray and others, 2003), but has also been documented at polythermal glaciers in Yukon, Canada (Clarke, 1976) and smaller surge-type glaciers in East Greenland (Jiskoot and Juhlin, 2009). According to the polythermal switch model of glacier surging, for glaciers with little to no sliding at the terminus, surging will initiate up-glacier from the terminus and the surge front will propagate down-glacier. This down-glacier propagation is supported by observations at land terminating surge-type glaciers Bakaninbreen (Murray and others, 1998) and Usherbreen (Hagen, 1987) in Svalbard. In contrast, for marine-terminating thermally regulated surge-glaciers, the ‘activation wave’ at the onset of a surge is faster than ice flow, and there is no observed down-glacier propagation of the surge front (Fowler and others, 2001; Murray and others, 2003). For surge-type polythermal glaciers in Svalbard, which represent the most extensively observed polythermal surge-type glaciers, the average repeat interval is estimated to be ~50–100 years based on the few glaciers with repeat surge observations.

The second hypothesis, the hydrologic switch, explains surge motion through the transition from an efficient drainage system with low water pressure, to an inefficient system with high basal water pressure and enhanced basal sliding (Kamb, 1987). The hydrologic switch mechanism for surging was initially proposed for a temperate hard bedded glacier (Kamb, 1987) but subsequent research suggests the plastic deformation of subglacial till can cause surge initiation (e.g. Truffer and others, 2000; Minchew and Meyer, 2020). For hydraulically regulated surges, the surge front propagates down-glacier (Meier and Post, 1969; Dolgoushin and Osipova, 1975). The hydrologic switch mechanism is supported by observations of subglacial water pressure and proglacial discharge at Variegated Glacier (Kamb and others, 1985; Kamb, 1987) as well as speed and thickness change observations for Variegated Glacier (Eisen and others, 2005), Lowell Glacier (Bevington and Copland, 2014), Bering Glacier (Fatland and Lingle, 2002), and Donjek Glacier (Kochtitzky and others, 2019) in Yukon/Alaska. For land-terminating Alaskan surge-type glaciers, where this is the dominant mechanism, the average surge recurrence interval is ~15 years (Harrison and Post, 2003; Sevestre and Benn, 2015). Repeat surge cycle observations of both speed and elevation change are only available for the four aforementioned glaciers, limiting our understanding of temperate glacier surge kinematics.

Recent work has focused on the development of unifying theories of surging able to explain the variety of dynamics across regions and irrespective of mechanism. From their statistical study on climatic and geometric controls on the distribution of surge-type glacier, Sevestre and Benn (2015) were the first to propose an enthalpy formulation for glacier surges. The enthalpy of a glacial system is defined by the internal energy, which is a function of liquid water content and temperature within the ice (Aschwanden and others, 2012; Sevestre and Benn, 2015). Sevestre and Benn (2015) suggest that climate regimes that inhibit a steady-state enthalpy balance are prone to surging. Benn and others (2019) expanded on the enthalpy-based framework of Sevestre and Benn (2015) by coupling equations for glacier thickness and enthalpy in a lumped parameter model. The enthalpy balance model produces the self-sustained oscillation in thickness and velocity typical of glacier surges (Benn and others, 2019). The abundance of non-surging glaciers in climate regimes favorable to surging suggests, however, that local geology and geometry (i.e. slope, substrate, etc.) exert a strong control on surging (e.g. Post, 1969; Clarke, 1991; Jiskoot and others, 2000; Crompton and others, 2018).

Here we make use of the exceptionally short surge recurrence interval of Sít’ Kusá in southeast Alaska, which has surged five times between 1983 and 2013, to investigate the kinematics of multiple surge events. Post (1969) noted the possible surge-type nature of Sít’ Kusá based on surface morphology, which was later corroborated by McNabb and Hock (2014) based on its multiyear cycles of advance and retreat. Our analysis of Sít’ Kusá is the first to look at the kinematics of the individual surge events in detail and represents the densest record of surge events in the satellite-era for a single glacier to date. Using terminus position change and velocity data for all surge cycles, and elevation data from 2001 to 2013, we show that surge events initiate in the winter and propagate down-glacier, providing support for the hydrologic switch model of surge initiation.

2. Study site

The southern coast of the St. Elias Mountains is characterized by large annual snowfall totals of more than 3 m water equivalent (Marcus and Ragle, 1970). These high snowfall totals contribute to fast-flowing glaciers (e.g. Burgess and others, 2013), many of which reach the ocean (e.g. McNabb and others, 2015). Sít’ Kusá (60° 02′N, 139° 39′W) is one of only a few known surge-type tidewater glacier in Alaska (Post, 1969; McNabb and Hock, 2014). The indigenous Tlingit name Sít’ Kusá means Narrow Glacier (Thornton, 2010). Sít’ Kusá lies within Wrangell-St Elias National Park in the St. Elias Mountains of Alaska, USA. The glacier is ~30 km long, ~2 km wide and covers an area of 177 km² (RGI Consortium, 2017). Terminating on the western bank of Disenchantment Bay, it has a maximum elevation of 4144 m a.s.l. and a median elevation of 1297 m a.s.l. (RGI Consortium, 2017). The climate is sub-polar maritime, as supported by mean annual temperature of 4.7°C and mean annual precipitation of 385 cm from 1980 and 2015 at the coastal town of Yakutat ~50 km to the southeast (http://climate.gi.alaska.edu/acis_data).

The surface of Sít’ Kusá is characterized by extensive crevassing and abundant debris cover. Heavy rock fall activity and erosion, characteristics of the St. Elias Mountains, are likely responsible for the abundant debris cover. The glacier's main trunk (lower 17 km) is fed by two sources: the northern (32.5 km long) and southern (24 km long) tributaries (Fig. 1). The northern tributary is marked by an ice-fall located 23 km from the terminus. Sít’ Kusá flows over a hanging valley, 2 km from its terminus, to sea-level, creating the unique terminus lobes to the north and south of the valley walls (Figs 1, 2). A partially-subaerial moraine is often evident at the terminus, limiting ocean access to the glacier. At present (2020), calving is observed only along the southern half of the terminus (Fig. 2a). Preliminary observations from a September 2020 field campaign suggest the glacier was surging at that time (Bartholomaus, 2020). Long-term progradation of Sít’ Kusá's morainal bank into Disenchantment Bay suggests that the sediment production and transport rates beneath the glacier are high throughout both the quiescent and active phases of the surge events (Goff and others, 2012).

Fig. 1. Sít’ Kusá (60° 02′N, 139° 39′W), St. Elias Mountains, Alaska. Sít’ Kusá is located in Disenchantment Bay, at the toe of Hubbard Glacier to the east. Flowlines shown here are used for velocity and elevation profiles. Black circles mark the distance from the maximum terminus position (ξ) in kilometers, plotted at 1 km intervals. Background imagine is from Sentinel-2 on 31 August 2018 projected on a UTM 7N grid.

Fig. 2. Area and velocity observations of Sít’ Kusá from 1984 to 2017. (a) Worldview 2 image (Imagery copyright 2016 DigitalGlobe, Inc.) from 10 May 2016 overlain by maximum glacier extent during the five (A₁– A₅) surge events (Table 2). Black circles mark the distance from the maximum terminus position (ξ) in kilometers, plotted at 1 km intervals. (b) Area change calculated relative to the first observation (4 June 1980). Uncertainty due to mixed pixels is indicated by gray shading. Vertical lines correspond to the date of the maximum terminus area mapped in (a). (c) Surface velocity. Sampling points are denoted by the flow following coordinate system ξ. Points are connected with straight lines only to make the visualization of the velocity variations easier and should not be interpreted as indicative of trends over time.

3. Data and methods

Our analysis used time series of glacier terminus position, surface velocity and surface elevation extracted from optical satellite imagery. We used the Landsat archive to map terminus position from 1980 to 2017 and surface velocity from 1984 to 2013. These observations were paired with sparser digital elevation models (DEM) constructed from a variety of satellite and airborne platforms to quantify surface elevation changes through time. Details on each of these datasets and the methodologies used for their analysis are presented below.

3.1. Terminus delineation

Four Landsat 2 Multi Spectral Scanner (MSS), one Landsat 3 MSS, seven Landsat 4 MSS, one Landsat 4 Thematic Mapper (TM), 163 Landsat 5 TM, 178 Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and 78 Landsat 8 Operational Land Imager (OLI) SWIR-1, NIR, Red false color composites were used to map the terminus position of Sít’ Kusá from 24 June 1980 to 4 June 2017. The orthorectified L1 T scenes were downloaded from the United States Geological Survey's Earth Explorer (https://earthexplorer.usgs.gov). Landsat 7 products containing scan-line errors were used if the scan line gaps were small enough that linear interpolation across the gaps did not distort the terminus shape with respect to the terminus geometry closest in date. Because of Sít’ Kusá's lobate terminus, we manually delineated the terminus position between two fixed tie points located on the lateral margins of the glacier (Fig. 2a) and then closed the terminus polygons to analyze changes in terminus area over time.

For the Landsat images used in our analysis, we assumed the mixed pixel effect (i.e. error in delineations due to pixels that overlap both the glacier terminus and the surrounding terrain/embayment) to be the largest source of uncertainty in terminus delineations. To account for this effect, we assumed that the mapped positions could be in error by up one pixel. We summed the number of pixel intersections along each delineation and multiplied by the pixel resolution of the corresponding image to produce an estimate of uncertainty (Silverio and Jaquet, 2005; Rivera and others, 2007). Thus, uncertainty in the terminus position is a function of the area of each pixel (3600 m² for Landsat 2–3 MSS and 900 m² for Landsat 4–5, 7 and 8; Fig. 2).

3.2. Velocity mapping

Surface flow velocities were mapped from 23 April 1984 to 7 October 2013 using 31 pairs of orthorectified level one Landsat 5 TM, 7 ETM+ and 8 OLI scenes. We used Band 4 (30 m) scenes for Landsat 5 and Band 8 (15 m) scenes for Landsat's 7 and 8 to produce our velocities. We did not use Landsat 7 products containing scan-line errors in our velocity analysis. Extensive cloud cover limited the imagery from which velocities could be extracted, requiring scenes from different path/row combinations to be used in our analysis.

Manual examination of nine scene pairs revealed registration errors that needed to be addressed prior to velocity extraction. To fix these, we matched the poorly coregistered scenes to the ground control scenes used by the USGS for the corresponding path/row (LE07_L1TP_061018_20010719_20160929_01_T1 for 61/18, LE07_L1TP_062018_20010608_20160929_01_T1 for 62/18). To find potential matches, we used normalized cross-correlation on a grid with 400 pixel spacing, matching a 101-pixel kernel within a 401-pixel search window. To avoid erroneous coregistration matches on the moving glacier surfaces, clouds or shadow, we masked glaciers using the Randolph Glacier Inventory 6.0 outlines (RGI Consortium, 2017), and cloud/shadow using the Landsat Quality Assessment (BQA) band provided with each scene. Potential static point matches were filtered based on the strength of the correlation and their fit to a 2-D affine transformation between the images estimated using random sample consensus (Fischler and Bolles, 1981), implemented using the scikit-image python package (Van der Walt and others, 2014). Using the successfully-matched control points and the TanDEM-X 90 m Global DEM (Rizzoli and others, 2017), we then transformed each scene using a first-order rational function model (RFM) transform for L1TP scenes, and a third-order RFM for L1GT/L1GS scenes (e.g. Tao and Hu, 2001). If necessary, the L1GS scenes were then orthorectified based on the methods described in Gao and others (2009).

After all scene-pairs were adequately coregistered, we used NASA's Ames Stereo Pipeline to perform normalized cross-correlation of repeat satellite imagery following the approach of Shean and others (2016). Normalized cross-correlation was computed in the spatial domain using a Gaussian pyramid approach, where correlation is computed on sub-sampled images and the disparities from sub-sampled images are used to seed finer-resolution disparity maps (Shean and others, 2016). As this Gaussian pyramid approach automatically determines the search window size, for each pair we chose the size of the correlation window (i.e. kernel) between 9 and 35 pixels based on the surface conditions and observed flow velocities that produced the most spatially-extensive velocity maps (Table 1). Images with limited surface features and/or low contrast (e.g. snow on surface of glacier or limited illumination) produced the best results with a kernel > 21 pixels. Images with more distinct surface features (e.g. heavy crevassing, debris cover) produced the best results when correlated with a kernel size < 21 pixels. Sensitivity tests of kernel size indicate that variations in kernel size have limited effect on quiescent velocity maps, as any difference in velocity between kernel choices were well within our uncertainty bounds (Fig. S1). For the active phases, larger kernel sizes resulted in slower velocities because more of the glacier's static margins were included in the kernels (Fig. S2). Since spatial smoothing decreases with kernel size, smaller kernel sizes resulted in faster but noisier active-phase velocities (Shean and others, 2016). Therefore, we manually selected the smallest kernel size that returned the most accurate velocities (i.e. sufficient spatial averaging to reduce random uncertainties but not enough to bias estimates; Fig. S2; Table 1).

Table 1. Landsat 5 TM (LT05), 7 ETM+ (LE07) and 8 OLI (LC08) scenes and the associated correlation window (i.e. kernel) size used to produce the velocity maps

Scene pair dates (DD/MM/YYYY) and paths (all scenes were from row 18) are also listed. Bias and uncertainty metrics for the magnitude components of the velocity are listed.

We estimated the error of our velocity products by measuring the apparent displacement (in the x and y directions) over bare ground, as determined by the union of a glacier mask (RGI Consortium, 2017), a bare ground mask (Hansen and others, 2013) and the BQA bands for the Landsat scenes. The effects of displacements in the x and y directions on our velocities were then corrected by subtracting the median displacement over the identified bare ground surfaces from the displacement measurements used to compute velocities. We used the median apparent displacement (over bare ground surfaces) as our bias metric and one median absolute deviation (MAD) as our uncertainty metric (Table 1), as they are robust against outliers caused by clouds, shadows and/or rising snow-lines. Nonetheless, the average bias (plus or minus average MAD) prior to displacement corrections was − 0.04 ± 0.60, − 0.06 ± 0.53 and 0.74 ± 0.55 m d⁻¹ in the x, y and magnitude components, respectively (Fig. 3; Table 1). These uncertainties apply to regions with abundant surface features (e.g. chaotic crevasses), including the ablation zone below ξ _N = 21 km and the accumulation zone above ξ _N = 25 km. Uncertainties are likely larger across the icefall, however, where the ‘standing-wave’ of crevasses causes the correlation algorithm to fail; a phenomena known as ‘surface-locking.’ Across the icefall ‘surface-locking’ results in erroneous speed estimates of ~0 m d⁻¹.

Fig. 3. Median offset (m d⁻¹) over bare ground in the x and y directions prior to displacement corrections for all velocity pairs. Marker symbol denotes the path/row combination of the scene-pair and color denotes the kernel size used for correlation.

After we corrected for median offsets over stable ground, we filtered our velocity products using a signal to noise ratio > 3.5 threshold (Paul and others, 2017). Then, velocities (x and y) were filtered for orientation following the approach used by Burgess and others (2012). First, we calculated the angle θ between each velocity vector v _i,j and the median vector within a neighborhood of width w = 7 pixels, $\tilde {v}^{w}_{i, j}$:

(1)$$\theta = \cos^{-1}\left({v_{i, j} \cdot \tilde{v}^{w}_{i, j}\over \left\vert v_{i, j}\right\vert \left\vert \tilde{v}^{w}_{i, j}\right\vert } \right).$$

We iteratively removed vectors whose θ was above a threshold of 24, 18 and 12 degrees. The remaining velocity vectors were filtered for magnitude so that $\left \vert v_{i, j} - \tilde {v}^{w}_{i, j}\right \vert$ $> 30\percnt$ of the mean velocity were removed. Finally, these filtered velocity products (in the x and y directions) were used to create the velocity maps, which were smoothed with a Gaussian kernel of size 7 × 7 pixels to minimize the influence of noise in our interpretation. Velocity measurements were manually extracted from individual pixels along the northern, central and southern flowlines (Fig. 1). We defined a flow following coordinate system ξ, originating (ξ = 0 km) at Sít’ Kusá's maximum terminus position. At the confluence of the two tributaries (ξ = 16 km), we split ξ into two separate coordinates, ξ _N and ξ _S, following the center flowlines of the north and south tributaries, respectively. In this coordinate system, the heads of the northern and southern tributaries are at ξ _N = 33 km and ξ _S = 24 km, respectively. We report the average and standard deviation in observed velocities for all three flowlines from the terminus to ξ = 16 km. The reported velocities were restricted to ξ ≤ 16 km so that the averages were unbiased by differences in flowline length and data coverage. On average, velocities were mapped for $86\percnt$ of the three flowlines over the lower 16 km, with poorer data coverage farther inland at higher elevations. Additionally, as described below, the upper portion of the northern tributary (ξ ≥ 23 km) was not affected by the surges and the inclusion of velocities from this portion of the glacier would bias inter-surge velocity analysis (Figs 6, 8).

3.3. Digital elevation models

Surface elevation changes were documented using four Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), one Satellite Pour l'Observation de la Terre 5 (SPOT5), one Interferometric Synthetic Aperture Radar Alaska (IFSAR-Alaska) and two Digital Globe WorldView-1 DEMs. We used MMASTER to create ASTER DEMs with ~10 m vertical uncertainty and 30 m spatial resolution from April 2001, May 2003, March 2006 and July 2012 stereo imagery (Girod and others, 2017). We used two ArcticDEM strips, derived from WorldView stereo pairs acquired in May and December 2013 (Porter and others, 2018). These 2 m-resolution DEMs have an estimated vertical accuracy of ~3 m (Noh and Howat, 2015; Shean and others, 2016). We downloaded a 5 m-resolution, 3 m vertical uncertainty IFSAR-Alaska DEM acquired in late summer 2012 from Earth Explorer (https://earthexplorer.usgs.gov/). Finally, we downloaded a September 2007 SPOT5 DEM with a vertical uncertainty of 6 m from the SPIRIT project (https://theia-landsat.cnes.fr; Korona and others, 2009). See Table S2 for full granule names and acquisition dates. All of the DEMs were resampled to 30 m spatial resolution using bilinear interpolation, and coregistered to the 2013 WorldView-1 DEM using the methods of Nuth and Kääb (2011). Average elevation uncertainties along the flowlines were calculated as:

(2)$$\varepsilon_{\rm h} = {1.96 \sigma_{\Delta \rm{z}}\over \sqrt{L/C}},\; $$

where σ_Δz is the standard deviation of the elevation change measurements, L is the length over which they were measured, and C is the autocorrelation length (500 m; Howat and others, 2008).

4. Results

4.1. Surge chronology overview

To explore the general characteristics of the surge events (A₁–A₅), we used terminus area changes to identify the initiation and termination of surge events (Fig. 2), since the terminus change record is more temporally dense than the velocity and elevation datasets (Table 1; Figs 2, 4). Because surges are defined by changes in flow, not necessarily terminus position, we would ideally use our velocity record to constrain the timing of surge events. For properly tidewater surge-type glaciers in Svalbard terminus retreat does not necessarily mean the end of the active phase, just a increase in frontal ablation (Mansell and others, 2012). However, cloud cover in coastal Alaska prevents the construction of a velocity record with sufficiently dense temporal sampling to confidently identify surge initiation and termination from changes in flow speed.

Fig. 4. Timing of surge initiation and termination for five surge events (A₁–A₅). The shaded area is the temporal uncertainty of when the changes occurred. The time when the terminus is advancing is denoted by the black brackets.

The rate of terminus position change is calculated as

(3)$${{\rm d}L\over {\rm d}t} = \bar{v}_{\rm s} - \bar{v}_{\rm a},\; $$

where L is the glacier length, $\bar {v}_{\rm s}$ is the width-averaged speed at the terminus, and $\bar {v}_{\rm a}$ is the width averaged rate of frontal ablation (sum of calving and submarine melting). We observe that terminus advance (A₁–A₅) consistently begins between February and June (Table 2; Fig. 4). Assuming there is no dramatic periodic decrease in $\bar {v}_{\rm a}$ that initiates the terminus advance every ~5 years, then increases in L (A₁–A₅) must be driven by an increase in advection toward the terminus ($\bar {v}_{\rm s}$). Thus, we hypothesize that the periodic advance of the terminus through the winter and into spring indicates changes in ice flow toward the terminus. Based on the terminus change time series, paired with velocity and elevation time series when available, we identify five surge events (A₁–A₅) since 1980 (Table 2; Figs 2, 4). The surge events have an average recurrence interval of ~5 years and an active phase of ~2 years (Table 2). Observations from each surge event are described in detail below.

TABLE 2. Transition between active (A) and quiescent (Q) phases of Sít’ Kusá from 1984 to 2013

The duration (in years) of the surge phase (as inferred by terminus change) and associated terminus change (km²) are listed with corresponding uncertainties. For the active phase, speeds reported are the maximum observed speeds with location of corresponding speeds. For the quiescent phase, we report mean velocities over the trunk ($\xi = 0 - 16\, {\rm km}$) of the glacier. Transition between phases are separated by a solid line. Date format is DD/MM/YYYY.

4.2. 1983–1986 Surge event (A ₁)

While surge events presumably occurred prior to 1983, we are not confident in the data quality and temporal density prior to 1980. We first observe terminus advance between September 1981 and September 1982 (Table 2, Fig. 2 and 4), with the terminus advancing until late summer 1984. Our first speed observation, in fall 1984, was 0.46 ± 0.37 m d⁻¹ over the trunk of the glacier (Fig. 5). In summer 1985 the mean speed was 2.58 ± 2.05 m d⁻¹, with the largest increase at the confluence ($\xi = 10 - 15\, {\rm km}$) of the tributaries (Fig. 5). Terminus advance began again in spring 1986 accompanied by a 11.4-fold increase in mean speed by fall 1986 (Table 2; Figs 4, 5). Terminus advance ceased in spring/summer 1987, ending a 5.6 ± 1.9 year surge event that caused the glacier to advance 2.0 ± 1.0 km² (Figs 2, 4; Table 2). Velocity maps from August 1987 show quiescent speeds over the whole glacier, except for slightly elevated speeds for ξ < 5 km (Fig. 5).

Fig. 5. 1984–1986 surge speed observations. Speeds extracted along the northern (a), central (b) and southern (c) flowlines with dates for all profiles shown in subplot d. Active phase speeds are denoted in bold on the legend in subplot (b). The vertical grey line at ξ = 9.25 km is the approximate location of the dynamic balance line. The vertical grey line at ξ _N = 23 km in (a) is the approximate location of the icefall located on the northern tributary. Date of image pairs plotted on the timeline, with the color of the boxes corresponding to the lines in subplots a–c.

4.3. 1991–1993 Surge event A ₂

The glacier appeared to be in quiescence in spring 1990, with a mean speed of 0.61 ± 0.81 m d⁻¹ over its trunk (Figs 6a–c; Table 2). Terminus advance began in May 1991, after retreating 3.4 ± 0.6 km² since 1986 (Fig. 4; Table 2). Terminus advance was accompanied by a 7.8-fold increase in mean speed in June 1991, as compared to spring 1990 (Figs 6a–c). Elevated surface speeds continued through 1992 and reached a maximum of 26.15 ± 0.82 m d⁻¹ at the terminus in April 1993 (Table 2; Figs 6a–b). Terminus retreat commenced again in 1993 and we observed quiescent speeds through summer 1994 (Table 2; Figs 2, 4, 6a–c). The surge event lasted 2.1 ± 0.2 years, and resulted in a 3.2 ± 0.6 km² increase in glacier area (Table 2; Fig. 2).

Fig. 6. 1991–1993 surge speed observations. Speeds extracted along the northern (a), central (b) and southern (c) flowlines with dates for all profiles shown in subplot b. Quiescent (active) speeds are colored in green (purple) and become darker in time. Active phase speeds are denoted in bold on the legend in subplot (b). The vertical grey line at ξ = 9.25 km is the approximate location of the dynamic balance line. The vertical grey line at ξ _N = 23 km in (a) is the approximate location of the icefall located on the northern tributary. Date of image pairs plotted on the timeline (d), with the color of the boxes corresponding to the lines in subplots a a–c.

4.4. 1999–2002 Surge event (A ₃)

The glacier was in quiescence through 1998, with a mean speed of 0.61 ± 0.61 m d⁻¹ over its trunk in June 1998. The terminus began to advance in April 1999, ending 3.8 ± 0.6 km² of terminus retreat since the termination of the last surge (Table 2; Figs 4, 2). By summer (14 July–6 August) 1999, the glacier speed increased 15-fold, reaching a peak of 12.52 ± 0.36 m d⁻¹ at ξ = 7.9 km (Fig. 7.) Elevated speeds continued through July 2000, reaching a peak of 17.74 ± 1.03 m d⁻¹ at the terminus (Fig. 7; Table 2). In June 2001, we observed terminus retreat, and by fall of 2001, we recorded quiescent-like speeds of 0.39 ± 0.36 m d⁻¹ over the glacier trunk (Figs 2, 4, 7). This portion of the surge lasted 2.1 ± 0.1 years, and resulted in a 4.1 ± 0.6 km² increase in glacier area (Figs 2, 4; Table 2).

Fig. 7. 1999–2002 surge speed and elevation change observations. Speeds and elevations extracted along the northern (a, d), central (b, e) and southern (c, f) flowlines with dates for profiles shown in subplots b and e. Quiescent (active) speeds are colored in green (purple) and become darker in time. DEMs are plotted with the color of the velocity map closest to the date range of the differenced DEMs. Date pairs exhibiting active phase speeds are denoted in bold on the legend in subplot (b). The vertical grey line at ξ = 9.25 km (ξ _N = 23 km) is the approximate location of the dynamic balance line (northern tributary icefall). Date of speed scene pairs plotted on the timeline (g), with the color of the (solid) boxes corresponding to the lines in subplots a–c. The thinner (dotted) boxes are the time span of DEM pairs.

Terminus retreat continued until early spring (3 February–16 March) 2002 (Figs 2 and 4), resulting in a total retreat of 1.1 ± 0.7 km² (Table 2). There was a ~14-fold increase in mean speed by early April 2002, as compared to quiescent speeds from summer 2000 (Fig. 7). Active-phase speeds continued through mid-April 2002, reaching 20 times the quiescent speeds, with down-glacier propagation of the surface speed maximum (Fig. 7). Terminus advance ceased during the summer of 2002, for a total advance of 1.2 ± 0.7 km² since March 2002 (Table 2; Figs 2, 4), and the glacier slowed to quiescent speeds by spring 2003 (Fig. 7).

From April 2001 to 2003, surface elevations decreased by an average of 6.9 ± 1.4 m from $\xi = 0 - 5\, {\rm km}$ across all three flowlines (Figs 7d–e, pink shading). Over the same time period, there was an average increase in elevation of 8.8 ± 0.4 m for all three flowlines above ξ = 5 km (Figs 7d–f, blue shading). Surface elevation increase along the northern flowline was confined to below the icefall at ξ = 23 km (Fig. 7d).

4.5. 2006–2008 Surge event (A ₄)

Terminus advance in winter 2006 ended a 5.6 ± 0.1-year quiescent period during which the terminus retreated 2.2 ± 0.8 km² (Table 2; Figs 2, 4). By May 2006, the glacier had sped up 22.6-fold over its trunk (Fig. 8), reaching a peak of 27.6 ± 1.2 m d⁻¹ at ξ = 10.9 km (Table 2; Figs 8a–c). Active-phase speeds continued through summer 2007, with a surface speed maximum (of 18.49 ± 0.32 m d⁻¹) at ξ = 2.7 km (Fig. 8). By summer 2008, velocities decreased ~65% relative to summer 2007, and the terminus began to retreat (Table 2; Figs 2, 4, 8a–c). The surge event lasted 2.4 ± 0.2 years and was associated with a 3.0 ± 0.8 km² increase in area (Fig. 2; Table 2).

Fig. 8. 2006–2008 surge speed and elevation change observations. Speeds and elevations extracted along the northern (a, d), central (b, e) and southern (c, f) flowlines with dates for profiles shown in subplots b and e. Quiescent (active) speeds are colored in green (purple) and become darker in time. Differenced DEMs (d–f) are plotted with the color of the velocity map closet to the date range of the differenced DEMs, with thickening shaded in blue and thinning in red. Date pairs exhibiting active phase speeds are denoted in bold on the legend in subplot (b). The vertical grey line at ξ = 9.25 km (ξ _N = 23 km) is the approximate location of the dynamic balance line (northern tributary icefall). Date of speed scene pairs plotted on the timeline (g), with the color of the (solid) boxes corresponding to the lines in subplots a–c. The thinner (dotted) boxes are the time span of DEM pairs.

From 20 May 2003 to 3 March 2006, a time period including both the quiescent phase and surge initiation, the glacier surface elevation increased by an average of 28.7 ± 1.1 m above ξ = 9.3 km (Figs 8d–f, blue shading). For all three flowlines, the surface elevation decreased over $\xi = 0 - 10\, {\rm km}$ by an average of 28.9 ± 1.5 m (Figs 8d–e, pink shading). The along-flow patterns of surface elevation change were reversed from March 2006 to September 2007, during the active phase of the surge. Surface elevations increased over $\xi = 0 - 7\, {\rm km}$ by an average of 54.5 ± 3.8 m, with a maximum increase of 141 ± 12.0 m at the terminus (Figs 8d–f, dark blue line). Above the dynamic balance line (ξ = 9.25 km), the surface lowered by 56.9 ± 2.0 m on average, with a maximum thinning of 123.6 ± 12.0 m at ξ _N = 18.1 km (Fig. 8d). Over both time periods (May 2003 to March 2006 and March 2006 to September 2007), little change occurred above ξ _N = 23 km, the approximate location of the icefall (Fig. 8d).

4.6. 2011–2013 Surge event (A₅)

The glacier was in quiescence for 3.7 ± 0.2 years, through spring of 2011, during which the glacier retreated 2.6 ± 0.8 km² (Table 2; Figs 2, 9a–c). While quiescent speeds were uniform below ξ < 10 km, we see an increase in speed at the confluence of the two tributaries ($\xi = 10 - 15\, {\rm km}$) through the spring of 2011 (Figs 9a–c). Terminus advance began in spring 2012.

Fig. 9. 2011–2013 surge speed and elevation change observations. Speeds and elevations extracted along the northern (a, d), central (b, e) and southern (c, f) flowlines with dates for profiles shown in subplot b and e. Quiescent (active) speeds are colored in green (purple) and become darker in time. Differenced DEMs (d–f) are plotted with the color of the velocity map closet to the date range of the differenced DEMs, with thickening shaded in blue and thinning in red. Date pairs exhibiting active phase speeds are denoted in bold on the legend in subplot (b). The vertical grey line at ξ = 9.25 km (ξ _N = 23 km) is the approximate location of the dynamic balance line (northern tributary icefall). Date of speed scene pairs plotted on the timeline (g), with the color of the (solid) boxes corresponding to the lines in subplots a-c. The thinner (dotted) boxes are the time span of DEM pairs.

We first observed active-phase speeds in June 2013, a ninefold increase from early-spring 2011, with a surface speed maximum of 19.04 ± 0.27 m d⁻¹ at ξ = 2 km (Figs 9a–c). Speeds decelerated, propagating downglacier, through the summer of 2013 (Figs 9a–c). Our final observation from fall 2013 shows quiescent speeds over the entire glacier, except ξ < 5 km where mean speeds were 1.68 ± 0.24 m d⁻¹. Terminus retreat began in early summer 2013, signaling the end of a 1.2 ± 0.2 year surge event, where the glacier advanced 2.6 ± 0.8 km².

Elevation observations from September 2007 to July 2012 (spanning the quiescent period) show a mean surface lowering of 59.93 ± 0.4 m from $\xi = 0 - 9\, {\rm km}$ and an average thickening of 31.6 ± 1.1 m above ξ = 9 km (Figs 9d–f). From July 2012 to late August 2012, the glacier surface lowered by 3.0 ± 0.4 m above ξ = 10 km and thickened by 32.9 ± 0.7 m below ξ = 10 km on average (Figs 9d–f). The May 2013 DEM only covers the lowest 9 km of the glacier, but shows an average thickening of 45.9 ± 4.3 m as compared to August 2012 (Figs 9d–f). From May to December 2013, surface elevation decreased by an average of 26.6 ± 1.2 m over ξ = 0 − 9 km (Figs 9d–f).

4.7. Inter-surge comparison

Sít’ Kusá surged five times between 1980 and 2017, making it one of the most active surge-type glacier currently known in the world (Table 2; Fig. 2; Sevestre and Benn, 2015). We find that, on average, surge events have an active phase of ~2 years and a recurrence interval of ~5 years (Table 2; Fig. 2). We observe as much as 20-fold increases in speed (e.g. A ₂, A ₄) during the active-phase. Preceding surge-initiation, we observe increased surface speeds at the confluence ($\xi = 10 - 15\, {\rm km}$) of the two tributaries (e.g. Figs 5, 6, 8, 9). Terminus advance consistently begins between February and May (Table 2; Fig. 4). During the active-phase of surge events, we observe surface-speed maximums propagate down-glacier (Figs 5–9). Elevation change observations suggest the boundary between the reservoir zone and the receiving zone (i.e. dynamic balance line) is at ξ ≈ 9 km (Figs 8 and 9). For the northern tributary, the inland extent of the reservoir zone coincides with the base of the icefall located at ξ _N = 23 km (Figs 7, 8). Surface elevation profiles from the 2006–2008 surge event (A₄) suggest the reservoir zone extends over the entire southern tributary above ξ = 9 km (Fig. 8f). Based on the extent of surface elevation gain during the 2006–2008 (A ₄) and 2011–2013 (A ₅) surge events, we infer that the receiving zone extends from ξ ≈ 9 km to the terminus.

5. Discussion

5.1. Probable surge mechanism

Our sparse sampling of velocities is dictated by extensive cloud-cover in southeast Alaska. Therefore, we note that our maximum speeds only represent what we are able to observe, and do not coincide with the actual maximum speed of the glacier during each surge event. Nonetheless, elevated speeds over the entire glacier (excluding the area above the icefall) during the active phase support our use of terminus position change as an indicator of surge advance, since enhanced flow toward the terminus will result in terminus advance in the absence of an equivalent change in calving and/or submarine melting (Eqn (3)). While McNabb and others (2015) document variations in frontal ablation for Sít’ Kusá, their estimates are based solely on surface velocities and changes in terminus position. Because McNabb and others (2015) did not have coincident elevation change observations, we cannot partition their rates of frontal ablation between oceanic forcing and glacier dynamics. We observe terminus advance consistently beginning between February and May (Table 2; Fig. 4), suggesting surge events initiate in late winter at Sít’ Kusá. Late winter surge initiation has been observed for a number of glaciers, including Variegated Glacier (Kamb and others, 1985; Eisen and others, 2005), Medvezhiy Glacier (Dolgoushin and Osipova, 1975), West Fork Glacier (Harrison and others, 1994), Bering Glacier (Roush and others, 2003) and Sortebræ (Pritchard and others, 2005), and is commonly attributed to pressurization of an inefficient subglacial hydrologic network and enhanced basal sliding (Kamb, 1987). For hydrologically-regulated surges, the surge front propagates down-glacier (Meier and Post, 1969; Dolgoushin and Osipova, 1975), which we also observe at Sít’ Kusá (Figs 5–9). Thus, we interpret the apparent winter initiation in conjunction with down-glacier propagation to suggest the hydrologic switch as the mechanism responsible for Sít’ Kusá's surges. Both observational (e.g. Clarke and others, 1984; Truffer and others, 2000; Woodward and others, 2003) and theoretical studies (e.g. Clarke, 1987b; Minchew and Meyer, 2020) suggest the deformation of till to be a critical component in the initiation and propagation of surge events, while till has also been proposed as the mechanism controlling the fast flow of icecaps in Iceland (e.g. Boulton and Hindmarsh, 1987; Kjær and others, 2006; Minchew and others, 2016) and ice streams in Antarctica (e.g. Blankenship and others, 1987; Tulaczyk and others, 2000). Given the high sediment production rates at Sít’ Kusá (Goff and others, 2012), it likely overlies a soft bed of deformable sediments, which could play an important role in surge initiation and propagation. We currently lack sufficient observations to make conclusions about the till, highlighting the importance of future work here to understand the subglacial environment.

5.2. Controls on surge recurrence intervals

While the kinematics of Sít’ Kusá's surge cycle are similar to its land-based Alaskan counterparts, its recurrence interval is exceptionally short. The typical surge recurrence interval in Alaska for temperate, hydrologically-regulated surge events is about 15 years (Sevestre and Benn, 2015). Variegated Glacier's surge recurrence interval is 13–18 years (Eisen and others, 2001), for Lowell Glacier between 11 and 18 years (Bevington and Copland, 2014), and 9–12 years for Donjek Glacier (Kochtitzky and others, 2019). Environmental conditions (thermal regime and precipitation) have been suggested as an explanation for regional differences in the length of the surge cycles between Svalbard (50–100 year quiescent period and 3–10 year active phases) and Alaska/Yukon (10–15 year quiescent period and 1–2 year active phases; Dowdeswell and others, 1991; Murray and others, 2003). However, environmental and geometric factors that influence reservoir evacuation also likely influence surge recurrence intervals. For example, the 1995 surge of Variegated terminated early with respect to previous surges (Eisen and others, 2005) and the next surge occurred only 9 years later in 2003/2004, a recurrence interval much shorter than the previously observed 13–18 years (Harrison and others, 2008). This anomalously short surge interval has been attributed to the early termination of the 1995 surge, which resulted in only partial evacuation of the reservoir zone and reduced the mass accumulation required to reach the critical stress threshold for surge initiation (Eisen and others, 2005; Harrison and others, 2008). While environmental conditions are clearly important, it remains difficult to isolate their impact on the surge cycle.

Given Sít’ Kusá's location in southeast Alaska, an area characterized by high accumulation rates (Marcus and Ragle, 1970), it is important to understand the influence of climate and mass balance on Sít’ Kusá's surge kinematics. While glacier surges are known to be caused by internal dynamics, climate is an important control on the distribution of surge-type glaciers (Sevestre and Benn, 2015), surge recurrence intervals (Eisen and others, 2001; Striberger and others, 2011), the vigor of the active phase (Frappé and Clarke, 2007) and regional changes in the number of surge-type glaciers (Dowdeswell and others, 1991; Copland and others, 2011). The link between climate and surge recurrence is supported by geologic evidence of glaciers that were formerly surge-type under more suitable climatic conditions (Hoinkes, 1969) and the correlation of surge recurrence intervals with centennial to millennial climatic variability (Striberger and others, 2011). Accumulation data, from one ablation stake, for Variegated Glacier shows a threshold of 43.5 ± 1.2 m must be reached for surge initiation to occur, a quantity that remains constant despite variations in the length of the recurrence interval (Eisen and others, 2001). Sít’ Kusá is located ~ 20 km away from Variegated Glacier, just across Disenchantment Bay. While Sít’ Kusá is a much larger glacier (177.0 versus 36.3 km²; RGI Consortium, 2017) and extends to higher elevations (Fig. 10) where we know little about the mass-balance profile, it experiences similar regional climatic conditions to Variegated. Therefore, if mass balance were the sole factor dictating a glacier's ability to reach the critical basal shear stress needed for surge-initiation, we would expect Sít’ Kusá and Variegated to surge at similar frequencies. The fact that they do not demonstrate that external forcings alone are insufficient to cause a glacier to surge or explain the heterogeneity in surge kinematics within a geographic region (i.e. Sít’ Kusá versus Variegated).

Fig. 10. Hypsometry for Sít’ Kusá and Variegated Glacier. Curves show the normalized distribution of ice area with elevation, such that the area under each curve is equal. Elevation data comes from IFSAR-Alaska DEM from August 2012. Long dashed lines indicate the mean elevations and short dashed lines indicate plus or minus one standard deviation.

5.3. Possible explanations for Sít’ Kusá's exceptional recurrence interval

Eisen and others (2005) hypothesize that surge events initiate when the reservoir zone fills sufficiently and reaches a critical basal shear stress. This assumes the driving stress equates to the basal shear (τ_b) stress:

(4)$$\tau_{\rm b} \approx F \rho g H \sin \alpha ,\; $$

where F is a shape factor, ρ is the ice density, g is the acceleration due to gravity, H is ice-thickness and α is the surface slope. This assumption has been shown to be invalid during surge-initiation (Schoof, 2005), when the driving stress and basal shear stress do not equate, due to rapid changes in effective pressure. Nonetheless, we use this crude approximation for a conceptual discussion of the quiescent period leading to surge-initiation. Changes in ice-thickness (H) with time are dictated by:

(5)$${\partial H\over \partial t} = \nabla \cdot Q + \dot b,\; $$

both the flux divergence ($\nabla \cdot Q$) and the mass balance ($\dot b$). Assuming the mass balance between Sít’ Kusá and Variegated are similar enough to be insufficient to explain the dramatically different surge-recurrence intervals, we suggest that geometric factors that control the flux-divergence ($\nabla \cdot Q$) play a role in Sít’ Kusá's exceptionally short recurrence interval.

Velocity and surface elevation observations, from the 2006–2008 surge in particular, show surge speed and mass redistribution confined to $\xi _N = 0 - 22\, {\rm km}$ by the icefall at ξ _N = 23 km (Figs 7–9). The observed disconnect in kinematics above and below the icefall is in line with the observations of large flow speed changes and mass redistribution below icefalls at a number of glaciers in Alaska (e.g. McNabb and others, 2012; Armstrong and others, 2017; Durkin and others, 2017; Enderlin and others, 2018) and previous observations of glacier surges confined by an icefall (e.g. Echelmeyer and others, 1987; Pritchard and others, 2003). This disconnect in kinematics likely inhibits draw down of ice volume from the accumulation zone during a surge event, as supported by the minimal observed elevation changes above ξ _N = 23 km during the 2006–2008 surge (Fig. 8). Although our velocity observations are temporally sparse, the data suggest that advection from the northern tributary's accumulation zone is either relatively steady or slightly elevated after surges (Figs 7–9), enabling rapid accumulation of ice in the reservoir zone between surge events. Due to the likelihood of surface locking across the icefall, which may bias our velocity time series in this location, we cannot confidently interpret temporal variations (or lack thereof) in speed across the icefall ($\xi _N\approx 22 - 24\, {\rm km}$). However, Sít’ Kusá's extensive crevassing enables successful offset tracking in the isolated accumulation zone above the icefall (Figs 5–9). It is unlikely that surge events initiate at the icefall, but the icefall does appear, in part, responsible for the rapid filling of the reservoir zone leading to surge initiation. The elevated surface speeds on the main trunk of the glacier (Figs 5, 6, 8, 9), which propagate down-glacier (Figs 8, 9), suggest surge events initiate somewhere along the confluence of the two tributaries ($\xi = 10 - 15\, {\rm km}$).

Flowers and others (2011) note a similar kinematic disconnect during surge events above and below a pronounced bedrock ridge of a small surge-type glacier in Yukon, Canada. While the kinematic disconnect does not alter the recurrence interval for the glacier in that study, prognostic numerical modeling simulations under a warming climate suggest the bedrock ridge contributes to the glacier's ability to form a reservoir zone even when climatic conditions are sufficiently negative to inhibit surging (Flowers and others, 2011). Therefore, glacier geometry, especially bedrock topography, not only affects surge characteristics, it may also influence (locally) if and how a glacier will surge under a warming climate.

In summary, our analysis of several surge cycles of Sít’ Kusá is consistent with a model of surging where surge propagation is facilitated by an inefficient subglacial drainage system (Kamb, 1987; Eisen and others, 2005; Harrison and others, 2008). We do not have sufficient data to fully assess the role of glacier geometry on Sít’ Kusá's surge characteristics, but we interpret the muted surface velocity and surface elevation changes inland of the northern tributary's icefall (at ξ _N = 23 km) as an indicator that the underlying topography strongly influences the characteristics of the glacier's surges. Although it is well known that an interplay between climate and geometry influence surge behavior, we recommend future investigation of Sít’ Kusá's exceptionally short recurrence interval with more detailed observational data to gain unique insights into controls on glacier surging.

6. Conclusion

From our analysis of optical satellite archives, we document five surge events of Sít’ Kusá between 1980 and 2013. Surge events lasted ~2 years on average with an average ~5-year recurrence interval between surges, representing the shortest known regular recurrence interval in the world (Table 2). We observe terminus advance to consistently begin between February and May, corresponding to hydrologic winter, for the 1991–1993, 1999–2002, 2006–2008 and 2011–2013 (A₂–A₅) surge events (Fig. 4). Velocity observations from the A₁–A₅ surge events (Figs 5–9) show down-glacier propagation of surface speed maximums over the course of each surge event. Surface elevation observations from the 2006–2008 (A₄) and 2011–2013 (A₅) surge events reveal the dynamic balance line to be 9 km from the terminus with the reservoir zone extending over the entire southern tributary while mass redistributions are confined below the icefall along the northern flowline at ξ _N = 23 km (Figs 8–9). Based on the winter initiation, down-glacier propagation, and geographic location, we believe Sít’ Kusá's surge events are hydrologically regulated. The kinematic discontent above and below the icefall (ξ _N = 23 km) may explain the rapid recurrence interval of Sít’ Kusá. In order to understand what drives the recurrence interval, and therefore surge initiation, more detailed investigation to the fluxes coming from the isolated reservoir zone is needed. Additionally, characterization of basal properties, glacier dynamics, and surface accumulation and melt during quiescent and active phases will facilitate an improved understanding of hydrologic controls on surging at Sít’ Kusá. Given similarities between Sít’ Kusá's surge characteristics and other Alaskan surge-type glaciers, we recommend that future studies of glacier surging leverage Sít’ Kusá's exceptionally short recurrence interval to advance not only the understanding of surge initiation and propagation at this particular glacier, but all glaciers that undergo dynamic instabilities attributed to rapid changes in subglacial hydrology.