Transient subglacial water routing efficiency modulates ice velocities prior to surge termination on Sít’ Kusá, Alaska

Abstract

Glacier surges are opportunities to study large amplitude changes in ice velocities and accompanying links to subglacial hydrology. Although the surge phase is generally explained as a disruption in the glacier's ability to drain water from the bed, the extent and duration of this disruption remain difficult to observe. Here we present a combination of in situ and remotely sensed observations of subglacial water discharge and evacuation during the latter half of an active surge and subsequent quiescent period. Our data reveal intermittently efficient subglacial drainage prior to surge termination, showing that glacier surges can persist in the presence of channel-like subglacial drainage and that successive changes in subglacial drainage efficiency can modulate active phase ice dynamics at timescales shorter than the surge cycle. Our observations favor an explanation of fast ice flow sustained through an out-of-equilibrium drainage system and a basal water surplus rather than binary switching between states in drainage efficiency.

1. Introduction

Glacier surges are drastic, human-timescale changes in glacier behavior. They are characterized by semi-periodic, multi-year oscillations in ice velocities despite consistent, seasonal changes in melt input (e.g. Meier and Post, 1969; Truffer and others, 2021). Slow ice flow during quiescent phases (~5 − >100 years) alternates with 5- to 100-fold ice velocity increases during short active phases (~1–30 years). Surge-type glaciers cluster geographically within an envelope of climatic conditions (Sevestre and Benn, 2015). Surge-type ice flow behavior is diverse but occurs on a continuous spectrum rather than within distinct categories, suggesting there is a unifying physical mechanism underlying glacier surging (e.g. Sevestre and Benn, 2015). From a theoretical perspective, recent years have seen considerable progress toward uncovering such a universal model of glacier surging (Terleth and others, 2021). Approaches toward a universal model have included process-based considerations of evolving friction at the glacier bed (Thøgersen and others, 2019; Minchew and Meyer, 2020). Process-based models are promising avenues forward and more widely applicable models are emerging (Beaud and others, 2022). However, they do not yet explicitly include the important influence of changing water fluxes to and from the glacier bed during the surge cycle, although they do emphasize the importance of basal water pressure (Thøgersen and others, 2019; Minchew and Meyer, 2020; Beaud and others, 2022). A more systems-based approach toward a unifying model of glacier surging is the enthalpy framework outlined in Benn and others (2019). While the enthalpy framework incorporates both polythermal and temperate glaciers, it simulates ice flow acceleration through increased basal water pressure with a simplified sliding law. As such, the importance of hydraulic forcing to the surge mechanism is universally acknowledged in recent theories of glacier surging. Additionally, observational studies repeatedly note the influence of water presence and pressure at the glacier bed in driving surge dynamics (e.g. Kamb and others, 1985; Murray and others, 2000; Kotlyakov and others, 2004; Benn and others, 2019).

Hydraulic forcing on ice velocities is not specific to glacier surging and largely depends on the bed's ability to evacuate water influxes from surface runoff (e.g. Iken and Bindschadler, 1986), which in turn depends on the subglacial drainage system's configuration (e.g. Kamb, 1987). A variety of possible subglacial drainage systems exist within the literature, each with specific characteristics. While there is a spectrum of geometries and behaviors, most proposed drainage configurations fit loosely within one of two broad categories. The first grouping of drainage systems is distributed and inefficient, including flow through a water film between the ice base and the substrate (Weertman, 1972), flow through porous substrates (e.g. Clarke, 1996; Flowers and Clarke, 2002; Kyrke-Smith and others, 2014), flow through poorly connected cavities (e.g. Lliboutry, 1968; Walder, 1986) or flow from the surface to unconnected cavities (e.g. Rada and Schoof, 2018; Nanni and others, 2021). Inefficient drainage systems tend to promote fast ice velocities due to their ability to sustain high basal water pressure. Increases in water pressure at the glacier base promote basal sliding through two mechanisms: changes in the ice-contact area with the bed surface by water-filled cavity growth (e.g. Iken, 1981; Anderson and others, 2004; Zoet and Iverson, 2015) and the dependence of subglacial till strength on effective pressure (e.g. Truffer and others, 2000; Tulaczyk and others, 2000; Iverson, 2010; Zoet and Iverson, 2020).

The second grouping of drainage systems includes configurations that transport water through localized and efficient channels forming at the glacier sole (Röthlisberger, 1972) or within the substrate (Nye, 1976). Channelized drainage systems adjust their morphology to changes in water influx from surface runoff and thus undergo only short lived (hours to days) increases in basal water pressure (e.g. Bartholomaus and others, 2008; Beaud and others, 2018). Efficient channelized systems tend to grow into dendritic patterns with limited spatial extent below glacier beds, as larger, lower pressure channels draw from smaller, higher pressure channels (e.g. Walder, 1986; Church and others, 2021; Nanni and others, 2021). Channel-like drainage systems are thought to evolve from distributed systems under sustained water supply (e.g. Hock and Hooke, 1993; Sundal and others, 2011). Once formed, they increase basal drainage efficiency and decrease basal water pressures, leading to a reduction in glacier velocities. This temporal evolution of basal drainage is a widely accepted mechanism for seasonal ice velocity changes, based on modeling studies (e.g. Schoof, 2010) and observational evidence (e.g. Tedstone and Arnold, 2012; Moon and others, 2014; Andrews and others, 2014). However, there are examples of channel-like systems in soft substrates that do not clearly transition to low pressure and high discharge regimes and that have the ability to restrict water flow over prolonged periods of time (Hock and Hooke, 1993; Walder and Fowler, 1994; Gulley and others, 2012; Hart and others, 2022).

The association between distributed, low efficiency drainage systems and high ice velocities hints at a mechanism explaining surging through persisting distributed and inefficient drainage. The most detailed description of such a hydrologically driven model of glacier surging was derived for the conditions of Variegated Glacier, Alaska (Kamb and others, 1985; Kamb, 1987). It suggests the active phase is sustained as long as there is a stable distributed drainage system of linked cavities that sustains high basal water pressures and does not adapt its morphology to changes in water supply. The surge terminates with the destabilization and collapse of these linked cavities in favor of a channelized drainage system, causing an abrupt release of the subglacial water volume (Kamb and others, 1985). The model's specificity to hard beds and the requirement of low water supply conditions prior to surge initiation somewhat limit its wider applicability (Harrison and Post, 2003). Explaining a wider range of surging behavior, such as surging under the presence of soft substrates (e.g. Hamilton and Dowdeswell, 1996; Truffer and others, 2000) or surge initiation during the melt season (Dunse and others, 2015; Sevestre and others, 2018), through changes in basal hydrology requires a reconsideration of the hydrologically driven surge model (Benn and others, 2022). In their observations of the 82–83 surge of Variegated Glacier, Kamb and others (1985) note large variations in borehole water level, in ice velocity and in terminus stream discharge prior to surge termination (their Figs 5, 9, 10). This variability suggests a complexity in the evolution of the drainage system during a glacier's surge phase and a resilience of high basal pressure to temporary episodes of water release that is not yet fully captured in the Kamb (1987) single hydrological switch model. New, modern in situ observations of the evolution of the subglacial drainage system are a critical avenue toward a truly universal and more detailed reconsideration of the hydrologically driven surge model (e.g. Truffer and others, 2021).

Here, we present observations of a well-instrumented surge on a temperate glacier in Alaska. We combine time-series of seismic observations, ice velocities and fjord water turbidity toward a partial record of subglacial drainage efficiency. Following a description of our data collection and the observational and model results, we devote the first part of our discussion to careful interpretation of each of the collected time-series signals (e.g. what time-series of seismic observations or of remotely sensed fjord color reveal about glacier behavior). In the second discussion section, following the attribution of observations to processes, we consider these processes in relation to one another and discuss the role of successive changes subglacial drainage in modulating surge dynamics. We place the significance of our findings in the context of previous work, and suggest potential implications for the surge mechanism. The complex variability in drainage efficiency prior to surge termination hints that conceptual models of drainage system evolution during the surge cycle may need expanding.

2. Study site

Our work centers around the 2020–2021 surge of Sít’ Kusá (briefly known as Turner Glacier), located on Tlingit land in the St. Elias range in Wrangell-St. Elias National Park, Alaska (Fig. 1). The ~30 km long and ~2 km wide Sít’ Kusá, which translates from Tlingit to ‘Narrow Glacier’, consists of multiple tributaries, with two main branches merging into a main trunk. This main trunk flows to sea level and terminates on a sediment shoal between surges and at tidewater in Disenchantment Bay during surge-driven advances. Surges initiate in the northernmost main tributary and propagate downglacier toward the terminus. Sít’ Kusá exhibits active phases of 1–2 years and quiescent phases of ~6 years, making it the most frequently surging glacier described in the literature (Nolan and others, 2021). The most recent surge initiated in March 2020, with velocities increasing from ~3 to ~25 m d⁻¹ in the lower northern tributary (Liu and others, 2024). An extensive array of instrumentation was installed on and around the glacier in late August of 2020 (Fig. 1). The surge front reached the glacier terminus in October 2020 (Liu and others, 2024), and the surge remained active until termination during the 2021 melt season, when surface velocities decreased to <5 m d⁻¹ and generally remained at ~1 m d⁻¹.

Figure 1. (a) Main components of deployed instrument network. Glacier extent and centerlines from RGI database. Background imagery is a Landsat-8 OLI scene acquired on 18 July 2021. Inset shows location in Alaska. Grid is in the coordinate reference system (CRS): UTM 7N, EPSG:32606. The datum is WGS 84. (b) Gantt chart showing temporal coverage of deployed network.

3. Data acquisition and analysis

3.1 Glaciohydraulic tremor

We aim to identify change in the subglacial hydrological system by monitoring seismic tremor, i.e. low amplitude seismic signals with consistent spectral content and durations of hours to months, around the glacier (Bartholomaus and others, 2015). In light of topographical constraints, twelve broadband seismometers were deployed as evenly spaced as possible at locations surrounding the main trunk of the glacier (Fig. 1). Sensors are named according to their placement on the East or West side of the glacier, and their distance in kilometers from the glacier terminus. All sensors were buried at depths of ~40 cm in glacier-proximal sediment. Eight stations provide high-quality, continuous records over nearly 24 months, while four stations suffered either wildlife damage (similar to that described in Tape and others, 2019) or other instrument malfunction (Fig. 1b).

Seismic stations included Nanometrics Trillium Compact Posthole and Nanometrics Meridian Compact Posthole seismometers, sampling at 250 Hz and with 20 and 120 s low-frequency corners, respectively. Throughout this study, we analyze instrument-corrected vertical-component data. Prior to deployment, all sensors were tested for uniformity and show inter-comparable seismic power (within ±0.2  dB) at frequencies above 0.1 Hz. Here we draw on data from three stations with continuous records that span the instrumented reach of the glacier. The record of these three stations is representative of the rest of our deployed network (Figs S1.1, S1.2), and our findings are reproducible with other stations.

We follow the methodology outlined in Bartholomaus and others (2015) and shared via Bartholomaus and Terleth (2023) to quantify the strength of seismic tremor that has previously been associated with glaciohydraulic sources. For each station, we compute the power spectral density (PSD) of 20 s windows with 50% overlap. We then compute the median power over 1 h long time-windows with 50% overlap. This yields a median valued PSD every 30 min (i.e. 48 PSDs in 24 h), as illustrated in the example spectrograms in Figures 2a, b. In Bartholomaus and others (2015), power within the 1.5–10 Hz frequency range is attributed to glaciohydraulic tremor. However, the 0.5–3 Hz frequency range is also influenced by calving events (O'Neel and Pfeffer, 2007; Bartholomaus and others, 2012). The high number of calving events in Disenchantment Bay noticeably impact the median spectra below 3 Hz (Figs 2c, d, Text S2, Fig. S2.1), so we focus on the frequency band of 3–10 Hz to isolate glaciohydraulic tremor and sum and standardize the PSD power within this band to obtain the glaciohydraulic tremor time-series shown in Figure 2c and in Figure S1.2. Similar considerations of frequency ranges >3 Hz have been used successfully to monitor glaciohydraulic tremor in Nanni and others (2020) and Lindner and others (2020).

Figure 2. Illustration of the analysis process that yields time-lags and velocities of a seismic tremor pulse. (a) Median spectrogram for SE7. (b) Median spectrogram for SW2. Dotted lines show frequency bounds between which we consider glaciohydraulic tremor. (c) Time-series of PSD amplitudes for SE7 and SW2 summed between 3 and 10 Hz in each 1 h time-window, with 30 min overlap. (d) Wavelet-based lag time estimation between signals recorded at SE7 and SW2, through time for different periods of oscillation. Positive (red) lags mean SE7 signal occurs before SW2 signal. Lags are plotted if coherence >0.7. Dashed lines show oscillation period corresponding to synoptic variability, 3–5 d. (e) Time-series of median lag between SE7 and SW2 for coherent signals with periods between 3 and 5 d.

3.2 Lags in the downglacier tremor signals

Beyond temporal changes, we are interested in observing spatial variability in recorded tremor. The spatial extent and density of the deployed seismic array is too wide to effectively conduct precise location tracking of tremor sources over time (e.g. Nanni and others, 2021; Labedz and others, 2022), but it does allow for inter-station comparisons of the PSDs in terms of temporal lag between glaciohydraulic tremor signals. We favor this approach over dominant noise source tracking (Vore and others, 2019) because for our purposes we are interested in the spatial propagation of tremor signals rather than the location of the highest amplitude tremor source. We use wavelet coherence analysis (Grinsted and others, 2004) to determine similarity between the time-series of station power spectral densities in time frequency space. We assume that the median spectral power between 3 and 10 Hz received at any given station is dominated by source(s) proximal to the station. This assumption is more valid for stations that are spatially distant: glaciohydraulic tremor signals have not been detectable at ranges >1 − 3 km in previous studies (Bartholomaus and others, 2015; Vore and others, 2019), meaning signals recorded at stations >2 − 6 km apart are likely independent. Wavelet coherence analysis produces two outputs in time and frequency space: (1) coherence values between zero and one, which reflect the similarity between the two signals, and (2) time-lag values which reflect the time-shift needed to obtain the highest similarity between the two signals (Fig. 2d). The time-lags are masked when the coherence between the signals at a given periodicity is below 0.7, to ensure the obtained lags are based on signals with a high degree of similarity. We focus on oscillations with a 3–5 d period as this captures the main variability within the glacio-hydraulic tremor signals (Fig. 2c). In order to obtain a time-series of time-lags, we integrate the time-lags over the oscillation period by taking the median value of the lags between the 3 and 5 d period bounds, drawn in dotted lines in Figure 2d. This yields a single lag time value for each period over which there are high coherence lags (Fig. 2e).

3.3 Surface runoff

To estimate variation in surface meltwater supply to the subglacial environment, we apply the Energy Balance Firn Model (van Pelt and Oerlemans, 2012) to Sít’ Kusá. The model solves the surface energy balance to compute surface temperature and melt values. The energy balance model is dynamically coupled to a physically based multi-layer snow and firn model that accounts for snow and firn pack densities, temperatures, water content and vertical liquid water transport (e.g. van Pelt and others, 2012, 2021). The model thus accounts for water retention in snow and firn, which can significantly impact the timing and volume of surface melt delivery to the englacial water system (Vallot and others, 2017; van Pelt and others, 2018; Alexander and others, 2020). Sít’ Kusá is estimated to receive as much as 7.5 m of annual precipitation (Simpson and others, 2005), much of it as snow, thus the incorporation of water storage in snow is an important component of the energy balance firn model. We expect there to be little delay in englacial water transfer during the considered time period as the glacier was already heavily crevassed in August 2020 (cf. Dunse and others, 2015; Gong and others, 2018), making the modeled surface runoff a relatively good estimate of the variability in water supply to the subglacial drainage system.

We force the energy balance firn model with meteorological data acquired at an automatic weather station located on Haenke Island (Fig. 1) that is maintained and operated by the Cold Regions Research and Engineering Laboratory (Finnegan, pers. comm.). Precipitation, cloud cover and relative humidity are required as model input data but are not recorded at the automatic weather station. For these variables we use data from the corresponding grid cell of the European Centre for Medium-Range Weather Forecasts reanalysis version 5 product (Fig. 1; Table S3.1). The energy balance firn model is distributed onto ArcticDEM, a high-resolution digital surface model of the Arctic (Porter and others, 2018), with a 32 m spatial resolution. We simulate the period between January 2017 and August 2022 and we assess model performance through its ability to reproduce surface elevation change derived from Worldview high-resolution satellite imagery acquired during the 2022 melt season (Text S3b, Fig. S3.1). We further include downglacier water routing by using the flow accumulation tool from the Matlab-based topotoolbox (Schwanghart and Scherler, 2014). We assume transfer of surface runoff to the glacier bed happens instantaneously through the severely crevassed glacier; we then use the subglacial hydropotential (Shreve, 1972) as input for the flow accumulation computation in order to estimate total surface runoff upstream of any given point on the modeled grid (Text S3c). We assume uniform water pressure at the overburden pressure and use bed topography derived by subtracting ice thickness modeled in Millan and others (2022) from ArcticDEM.

3.4 Subglacial water discharge at the terminus

Previous work has shown the feasibility of using remote-sensing imagery to obtain a qualitative understanding of frontal water release through time (Chu and others, 2009; McGrath and others, 2010; Tedstone and Arnold, 2012; Schild and others, 2017; Benn and others, 2019). The Sít’ Kusá terminus sits on a sediment shoal with water depths <40 m (Goff and others, 2012), the edge of which limits its advance into Disenchantment Bay during active phases. While we do not have a quantitative record of proglacial water discharge, we use sea surface characteristics of the area in front of the terminus of Sít’ Kusá as a proxy for relative changes in subglacial discharge. During quiescent phases, a calving embayment consistently forms on the southern half of Sít’ Kusá's calving front (Fig. 3a). The formation of such embayments has been attributed to subglacial discharge release in previous work (Sikonia and Post, 1980; Fried and others, 2018). At Sít’ Kusá, the location of this embayment coincides with the most likely discharge channel based on hydropotential mapping (Fig. 3a). Therefore, we manually delimit a ~6.5  km² area of sea surface in front of the calving embayment and use this region as our area of interest (AOI) to assess water discharge. We use the average pixel values within the AOI for band 7 of the Sentinel-3 Ocean and Land Color Imager (OLCI) Level-1b product as a proxy for sediment loading (ESA, 2022c). Band 7 records radiance within wavelengths of 615–625 nm (orange light in the visible spectrum), the spatial resolution of the pixels is 300 × 300 m, and at the latitude of Sít’ Kusá there is a temporal resolution below 1 d. To increase the robustness of the derived time-series of sediment loading, we also investigate the surface reflectance provided by the Sentinel-2 Multi-Spectral Instrument Level-2A band 4 centered at 665 nm (ESA, 2022b). The reflectance of water surfaces at these wavelengths scales closely with water turbidity (e.g. Schild and others, 2017; Hossain and others, 2021) and should thus reflect relative changes in turbidity in our area of interest. All imagery is filtered for clouds using the provided data flags. The two time-series are well-correlated (r = 0.76) between January 2019 and June 2022 (Fig. 3).

Figure 3. (a) False color showing surface reflectance in band 4 of Sentinel-2 MSI instrument. Image acquired on 27 July 2020, ~5 months after surge initiation and 13 months before surge termination. Sít’ Kusá RGI outline shown in purple and shore shown in gray. Later in the surge the terminus advances entirely into the bay. Orange polygon shows area of interest over which observations are averaged. Blue pixels show hypothesized subglacial flow pathway based on flow accumulation analysis. (b) Modeled surface runoff. (c) Radiance in Sentinel-3 OLCI band 7. (d) Surface reflectance in Sentinel-2 band 4. (e) Sentinel-1 Synthetic Aperture Radar Ground Range Detected Vertical-Vertical-polarized back-scatter. Time-series show individual data points and a 10-point moving average.

To further assess the validity of our record as a proxy for sub-aqueous frontal discharge, we also create time-series of backscatter in vertical–vertical polarization from the Sentinel-1 C-band synthetic aperture radar (ESA (2022a); Fig. 3e). The intensity of back-scatter from a water surface is expected to increase primarily with the presence of ice within the water (e.g. Ferdous and others, 2018; Benn and others, 2019). As such, we expect the backscatter intensity to increase with calving during the melt season but to decrease when large volumes of meltwater are released as freshwater upwelling would push icebergs out of our defined area of interest (Bartholomaus and others, 2013).

3.5 Surface velocities

We include data from five Trimble Net R9 Global Positioning System (GPS) receivers and Septentrio PolaNt-x MF antennas provided by EarthScope Consortium which were deployed and recovered on the ice surface at various times and locations (Fig. 1). We processed the data with the Canadian Spatial Reference System Precise Point Positioning algorithm. We averaged the positions at daily intervals and differenced consecutive positions to compute daily ice surface velocities. Surface velocities are not spatially uniform throughout the surge (Liu and others, 2024), such that these point velocities cannot be spatially extrapolated. Nevertheless, the records provide insight into velocity fluctuations at high temporal resolution. We supplement the in situ velocity measurements with satellite-image derived velocity estimates made with the open-source autoRIFT package (Lei and others, 2021). We include surface velocities in the upper trunk (15 km from the terminus) from pixel displacements derived from pairs of optical images (Sentinel-2 and Landsat-8) with date separations between 5 and 60 d and Synthetic Aperture Radar Imagery (Sentinel-1A and -1B) with date separations of 12 d (Liu and others, 2024).

4. Results

We present time-series of the modeled surface runoff and seismic tremor (Fig. 4a), of coherence and time-lags between tremor signals (Fig. 4b), of glacier surface velocity (Fig. 4c) and of water turbidity in front of the terminus (Fig. 4d). The time-series start in late summer 2020, ~6 months after the start of the surge, and continue to August 2022, ~12 months after surge termination. The sections below describe these time-series in further detail.

Figure 4. (a) Modeled estimate of surface runoff on Sít’ Kusá at a location near SE7 (UTM zone 7, 571931 E, 6658040 N) and median seismic power recorded at SE7 in 3–10 Hz frequency range. Five-day moving averages of plotted as thicker lines. (b) Time-lags between glaciohydraulic tremor signals between station pairs SE15–SE7 and SE7–SW2. Only values with coherence above 0.7 are plotted. (c) Surface velocities recorded at various on ice GPS receivers and through satellite image pairs. G12 and G9 overlap with G15 and G11 and are difficult to discern on the figure. (d) Radiance recorded in Sentinel 3 OLCI band 7. Shaded area in the time-series show the extent of the 2020–2021 active phase.

4.1 Simulated surface runoff

The timing of the onset and end of surface runoff is relatively consistent from year to year in our study period, around 15 April and 15 October, respectively (Fig. 4a). The highest runoff rates are generally reached around mid-June and persist until late August. The runoff during the 2021 melt season is not abnormally high (Fig. 3b), with several spikes of ~37  m³s⁻¹ occurring on 30 June and on 13–15 August which is the maximum runoff rate during the 2021 melt season.

4.2 Seismic tremor signal

The long-term tremor signal shows high amplitude during the 2020 surge winter and during the 2021 melt season. A strong decrease in amplitude during the 2021–2022 winter is followed by renewed higher tremor amplitude during the 2022 melt season. There is a sudden but relatively small decrease in tremor amplitude (~−150–~−155  dB) in the glaciohydraulic tremor window in early October 2020, at the end of the melt season (Fig. 4a, Fig. S1.2). However, tremor levels remain comparatively high throughout the winter of 2020–2021 with small variations of ±2 dB and a gradual increase coinciding with the onset of the 2021 melt season. The variability in glaciohydraulic tremor amplitude increases and remains high from June 2021 until the end of our record (the standard deviation in the SE7 tremor signal between 15 September 2020 and 1 June 2021 is 1.63 dB and between 1 June 2021 and 1 August 2022 the standard deviation is 7.30 dB). The 13 August runoff event coincides with a strong spike in tremor, followed by a net drop in tremor amplitude relative to pre-termination levels (~−150–~−155 dB). The tremor signal decreases from ~−150 to ~−167 dB between 15 August 2021 and 1 November 2021, closely mirroring the gradual decrease in surface runoff from its summer maximum to the end of the melt season (Fig. 4a, Fig. S4.2). During this period, brief spikes in surface runoff coincide with spikes in the tremor signal. Tremor remains low (~−167 dB) throughout the winter but the signal continues to show variability of ±5  dB that appears unforced by surface runoff. After the onset of melt in 2022, the tremor amplitude increases to ~−153 dB, about 3 dB below levels recorded during the surge. The rate of increase closely follows the rate of increase of surface runoff with a more gradual increase between 15 April and 15 May followed by a more rapid rise between 15 May and 5 June. Seasonal median noise levels recorded at SE7 (Fig. 5) also indicate that melt season tremor amplitude is relatively similar during the surge (May–August 2021) and post-surge (May–August 2022) but that winter tremor amplitudes are much higher during the surge (December 2020–March 2021) than after the surge (December 2021–March 2022).

Figure 5. Power spectral density functions (McNamara and Buland, 2004) with 50^th percentile values plotted from seismic noise recorded at SE7 showing distribution and median power for time-periods (a) in winter (1 December–1 March) and (b) during the following melt season (1 May–1 August). Purple and blue lines differentiate between the active phase (winter 2020–2021 and summer 2021) and the subsequent quiescent phase (winter 2021–2022 and summer 2022). Green shading indicates frequency range within which we consider glaciohydraulic tremor.

4.3 Lags between seismic tremor time-series

There is high coherence between seismic tremor signals within both the SE15–SE7 and SE7–SW2 station pairs during the surge with brief (~4 weeks) disruptions in coherence for the lower trunk during March and early April 2021 (Fig. 4b, Fig. S5.2). Coherence in the upper trunk frequently breaks down after surge termination. In the lower trunk, coherence is maintained during the 2021–2022 winter but largely absent during the 2022 melt season.

The time-lags between coherent tremor signals observed at SE7 and SW2 (Fig. 4b), vary between ~8 and 20 h from August 2020 to April 2021, with a faint increasing trend. The disruptions in coherence for the lower trunk occurring in March–April are accompanied by a brief period of negative time-lags (~−3 h). Meanwhile, time-lags in the upper trunk show a more clear but non-monotonic increasing trend from ~2 h during the 2020 melt season to ~20 h by April 2021. Time-lags in both the upper and lower trunk decrease between April 2021 and July 2021 and vary between −10 and 10 h during the second half of the 2021 melt season.

After surge termination, limited moments of high coherence in the upper trunk are marked by a wide range in time-lags with values ranging between −10 and 20 h. Lower trunk values vary between ~−1 and ~10 h for the remainder of the melt season and during early winter. Between January and April 2022, we observe time-lags in the lower trunk that range between −5 and 10 h. Despite these variations, the time-lags during the 2021–2022 winter are consistently lower than those during the 2020–2021 winter for both the upper and lower trunk.

4.4 Remote-sensing time-series

There is a consistent seasonal signal in the radiance values recorded with Sentinel-3 OLCI Band 7 (615–625 nm), with turbidity decreasing shortly after the end of the melt season and increasing around mid-February of each year (including during the 2020 and 2021 surge winters), ~2 months before the onset of the melt season (Fig. 4d). This seasonal signal is also present in the temporal evolution of the Sentinel-2 MSI surface reflectance although the latter record is limited in temporal resolution. The Sentinel-1 record shows a seasonal peak in back-scattering that precedes the peaks in the other two signals. Both the Sentinel-1 and Sentinel-2 records show a spike during the first half of 2021, when calving rates are maximal (Fig. S2.1). In the higher temporal resolution Sentinel-3 signal, yearly average radiance is only slightly higher in 2021 (1.25 W m⁻¹ sr⁻¹ μ m⁻¹) than in 2019 (0.91 W m⁻¹ sr⁻¹ μm⁻¹) and 2020 (1.04 W m⁻¹ sr⁻¹ μ m⁻¹). The total amplitude of these variations is ~0.34 W m⁻¹ sr⁻¹ μ m⁻¹ relative to an average seasonal standard deviation of 0.76 W m⁻¹ sr⁻¹ μ m⁻¹ (standard deviation of 0.14 W m⁻¹ sr⁻¹  μ m⁻¹ between seasons). Finally, we do not observe abnormal peaks or steps in turbidity during the 2021 melt season.

4.5 Surface velocities

Figure 4c shows the evolution of surface velocities from 2020 to 2022 recorded through GPS measurements and displacement observed in satellite imagery (Liu and others, 2024). Velocity measured at the G12 GPS station decreases from ~20 to ~10 m d⁻¹ in early October 2020, coinciding with the end of the melt season. However, ice velocities increase again within a few days and the satellite record indicates velocities reached ~20 m  d⁻¹ by mid-February. The velocities derived from satellite observations further indicate a gradual slowdown from nearly ~20 m d⁻¹ in June 2021 to ~10 m  d⁻¹ in mid-August 2021. Thirteen August 2021 is marked by a brief peak in ice velocities to ~15 m d⁻¹ followed by a sudden decrease to velocities generally below 5 m d⁻¹.

It is noteworthy here that there are high amplitude (±5 m d⁻¹) variations in surface velocity during the surge on timescales of 1 to ~10 d that are captured by the GPS record but missed in the satellite-derived velocity estimates. The amplitude of this variability is reduced, but remains present, at the end of the 2021 melt season during which surge termination occurs. We note two very short lived spikes to over 10 m d⁻¹ in late September 2021 in the G11 record. Velocities remain very low (~1 m d⁻¹) during the 2021–2022 fall to early winter and increase to ~3 m d⁻¹ from February 2022 until early July 2022. The velocities then decrease to ~1 m d⁻¹ by mid-July 2022.

5. Discussion: signal attribution

This work relies on a wide range of observation sources that complement each other toward a picture of ice flow and subglacial drainage behavior. Many of the time-series are proxies for the glaciological quantities we are interested in. We devote this initial part of the discussion to an assessment of the reliability of, and possible caveats to, our observations.

5.1 Relation between modeled surface melt and surface runoff

The energy balance firn model used to estimate surface runoff has several free parameters, notably including the rain to snow transition temperature, the elevation-dependent precipitation gradient and broadband albedo values (e.g. van Pelt and Oerlemans, 2012). Additionally, energy balance models are strongly dependent on the quality of meteorological input data. To assess the accuracy of our model estimates we compare the modeled surface elevation change on Sít’ Kusá at several dates with five digital elevation models derived from Worldview images acquired during the 2022 melt season. We find a correlation coefficient of 0.88 and a root mean squared error of 1.3 m (Text S3b). We tolerate this error extent for our purposes as we focus on the relative variation in surface runoff rather than absolute surface mass-balance estimates. The advantage of including a consideration of runoff buffering in the snow and firn pack outweighs the simplicity of a positive degree day approach.

5.2 Relation between surface runoff and subglacial water delivery to the glacier bed

We observe that peaks in modeled surface runoff consistently coincide with spikes in seismic tremor power (Fig. 4, Fig. S4.2). Several studies have noted how increases in tensile stresses and extensive crevassing during surges promotes penetration of surface melt to the bed (Dunse and others, 2015; Sevestre and others, 2018; Gong and others, 2018). Our data cover broadly the latter two-thirds of the surge: extensive crevassing was already present in August 2020 at the time of sensor installation, allowing a direct and widespread connection between the glacier surface and the bed. Supraglacial water routing would have been extremely limited by the crevassed nature of the surface during the surge (Fig. 6, Text S3c). Figures 4 and 8 provide perspective on the short-term response of ice velocity to changes in surface melt supply. An example here is the peak in surface melt driven by a rainstorm on 13 August 2021 that resulted in a brief threefold increase in ice velocity (Fig. 4). Such a direct response of velocity to water supply echoes observations made on non-surge-type alpine glaciers (e.g. Iken and Truffer, 1997; Bartholomaus and others, 2008), and shows that volumes of surface melt contributions to the subglacial water budget are sufficient to impact the system. Thus, while delivery of surface melt to the bed might still be lagged by several hours relative to estimated runoff time, we take modeled surface runoff as the best available proxy for water delivery to the glacier bed.

Figure 6. Photo of the glacier surface taken on 8 September 2020, ~10 km from the terminus along the main trunk of Sít’ Kusá, looking southeast toward the terminus. Inset shows location photo was taken. Photo by T.C. Bartholomaus.

5.3 Relation between seismic tremor and subglacial water flow

As with surface streams, the mechanism driving ground motion in subglacial conduits is thought to be a combination of the drag between turbulent water flow and conduit roughness and the rolling and saltation of sediment within the conduit (e.g. Tsai and others, 2012; Gimbert and others, 2016). As a result, the amplitude of the seismic tremor varies with both the water velocity and the sediment flux through the conduit, with the respective contributions of these processes difficult to disentangle (Gimbert and others, 2014). Understanding the drivers of seismic noise produced in streams remains an active field of study (e.g. Bakker and others, 2020) but previous work has shown that useful information on subglacial water flow can be derived from seismic tremor by neglecting contributions from sediment motion (e.g. Nanni and others, 2020, 2022) or by remaining agnostic regarding the exact source mechanism (e.g. Bartholomaus and others, 2015).

A possible alternate source of tremor is frictional stick-slip tremor at the ice-bed interface (Lipovsky and Dunham, 2017), echoing behavior observed along subduction zones (e.g. Shelly and others, 2006). Podolskiy and others (2021) find tidally modulated changes in seismic noise within the 3–14 Hz frequency range that are best explained as sourced from changes in basal sliding speed along a glacier bed/till interface. Stick-slip tremor, correlated with surface motion, has also been identified using geophones installed within 50 m of the glacier bed (Köpfli and others, 2022). There, the spectral content (chiefly >10 Hz) was tightly banded and varied with fluctuating basal water pressures. If the tremor signal observed on Sít’ Kusá was modulated by sliding rates, we would expect the correlation between tremor and sliding to be most clear in the absence of strong glaciohydraulic tremor. Contemporaneous recordings of GNSS-based surface velocity at G14 and seismic tremor at the adjacent SW14 at the close of the melt season (from 1 October 2020 to 15 November 2020) do not reveal a significant correlation between the two time-series (r = 0.17; Fig. S4.1). We also do not identify any shifts in frequency content (i.e. gliding) within the tremor we record that would be consistent with stick-slip tremor (Köpfli and others, 2022; Lipovsky and Dunham, 2017).

Instead, we record clear increases in tremor amplitude during the melt season (Figs 4a, 5) and find good agreement between variations in modeled surface runoff and seismic tremor (Fig. 4a; r = 0.72 for the 2021 melt season, 15 April to 15 October 2021; and r = 0.52 for the whole record), as reported at other glaciers (Bartholomaus and others, 2015; Vore and others, 2019). Even at the close of the melt season, with waning water influx, the correlation between tremor and melt is far stronger than that between tremor and ice flow velocity (r = 0.69 as compared with r = 0.17; Fig. S4.1). Furthermore, the tremor time-series at upglacier and downglacier stations reveal coherent variations in power that lag each other and propagate downglacier with celerities expected for water flow (see next section). As such, we infer that the amplitude of the tremor signal is driven primarily by the hydraulics of subglacial waterways beneath Sít’ Kusá and scales with water velocity through these waterways (Bartholomaus and others, 2015). We have not attempted to resolve which hydraulic process(es) may specifically control variations in tremor power at Sít’ Kusá – whether turbulent water flow or sediment transport – and expect that stick-slip tremor may be present within our time-series at some level. However, from the preponderance of evidence we interpret seismic tremor variations as a proxy for subglacial water flow.

5.4 Relation between seismic tremor time-lags and subglacial drainage system configuration

The time-lags between glaciohydraulic tremor recorded at different locations along the glacier require careful interpretation. Following Grinsted and others (2004), we compute time-lags only when there is a high level of coherence between the seismic tremor time-series. Despite their similarity, the source mechanisms and source locations of the tremor time-series must be independent: the time-lags between the time-series are on the order of hours, which is much longer than if the tremor time-series recorded at each station were sourced by a single process and the lags were caused by travel times of seismic waves. Glacio-hydraulic tremor generally attenuates to undetectable levels within ~3  km (Bartholomaus and others, 2015; Vore and others, 2019), leading us to expect that the recorded tremor was generally sourced in close proximity to the respective receivers. We focus on the SE15–SE7 and SE7–SW2 station pairs to maximize the distance between likely tremor source locations.

While Figure 2 illustrates a straight-forward time lag signal extracted from a portion of the seismic record, the full time-series show more complex behavior. Notably, the occurrence of negative lags suggests the tremor time-series recorded at the downglacier station can be ahead of the tremor recorded at the up-glacier station. A closer inspection of the tremor time-series during key time-periods provides some insight into the various situations producing these time-lags between SE15 and SE7 (Fig. 7), as described below. Supplementary material Text S5 provides a similar evaluation for SE7–SW2 along with the full time-frequency plots of the time-lags between both station pairs (Figs S5.3, S5.4).

Figure 7. (a) Three to five-day period phase lags between seismic tremor signals recorded at SE15 and SE7 shown on left hand axis in purple. Modeled surface runoff shown in gray on right-hand axis. (b)–(i) Detail views of four highlighted periods. Each colored frame of two panels shows the two tremor signals and the modeled melt signal in the upper panels (b, c, c, g) along with their associated time/frequency lag plots in the lower panels (d, e, f, i).

During March 2021, in the absence of surface runoff, tremor at SE15 shows a remarkable pattern including a gradual increase of ~1 dB lasting at least several days followed by a sudden ~4 dB spike in tremor that lasts ~48 h. This pattern is mimicked at SE7 ~22 h later and with lower amplitude (~2 dB peaks). If the seismic tremor is indeed predominantly hydraulic in origin, this similar tremor pattern strongly hints at a pulse in water velocities that travels downglacier. These time-lags are some of the longest in our record and we note a lack of coherence between signals before and after the highlighted period, suggesting a lack of connection in the drainage system during those surrounding time periods. The pulse-like pattern within the tremor time-series could be the signature of water being released out of overwinter storage (Liu and others, 2024), likely under high pressure, as this could generate the strong increases in tremor power (e.g. Gimbert and others, 2016; Nanni and others, 2021). Such events would likely have an expression on the ice surface, including a brief speedup and surface uplift (e.g. Iken and Bindschadler, 1986). Unfortunately we do not have surface GPS measurements during this time period. Echoes of similar behavior are present sporadically in our record (Text S5) but none are as clear as shown here.

During early February 2022, we find evidence of a likely similarly poorly drained subglacial environment. Short lived and strong pulses (>5 dB) in tremor recorded at SE15 do not materialize at SE7. Furthermore, the tremor pulses that do exist in the records of both seismic stations occur at SE7 5–10 h earlier than at SE15. The travel of tremor pulses in the up-glacier direction might be caused by a pressurized drainage system that ‘backs-up’ water, where the drainage system cannot accommodate the water influx and forces a water pressure pulse to migrate up-glacier (e.g. Barrett and Collins, 1997; Bartholomaus and others, 2008). Another example of this suspected behavior is shown in Figure S5.2. The behavior changes after 5–10 February 2022, a brief warm period during which there is up to 10 m s⁻¹ of surface runoff on the glacier, an infrequent event during the winter time. Both tremor signals respond to the water supply and the SE15 signal leads the SE7 signal afterwards, suggesting that perhaps the additional water led to an increase in drainage efficiency.

During June and July 2021, there is a transition from ~20 to ~5 h lags. Both tremor signals show daily fluctuations that follow the diurnal melt cycle and time-lags at a 1 d period are close to zero. This suggests a well-connected system where the tremor signals are a combined product of subglacial water flow downglacier and distributed water input from surface runoff. Interestingly, the long period (>5 d) increase in both tremor signals, forced by a long-term increase in surface runoff, occurs earlier at SE7 than at SE15. This difference in timing is likely caused by differences in snow cover thickness as a thinner snowpack would have a lower capacity to retain surface melt (e.g. van Pelt and others, 2016), leading meltwater reach the subglacial environment more quickly lower on the glacier.

Similar behavior is more obvious during June and July 2022, with the SE7 signal consistently leading the SE15 signal. Both tremor time-series show diurnal variability that follows the runoff signal but the amplitude of these oscillations is higher for SE7 than for SE15. Occasional peaks in SE15 tremor do not materialize in the SE7 record. We suggest that during June and July 2022 the time-lag signal is dominated by distributed surface water supply to the subglacial system: localized runoff quickly penetrates to the glacier bed and affects the local tremor signal. When we detect negative lags, the effect of the distributed water supply on the drainage system is strong enough to ‘drown out’ the tremor generated by downglacier water transport. This would be aided by poor along-flow connection within the drainage system but could be caused by very high surface runoff rates. There are frequent examples of this occurring on diurnal timescales (Fig. 7, Fig. S5.2).

The above examples, along with further examples and the full time lag time-series shown in Figures S5.3 and S5.4, shape the interpretation of the time lag time-series. Strongly positive lags (~5–25 h) hint at connected but slow moving subglacial drainage in the downglacier direction. Slightly positive lags point at a connected drainage system with some efficient component through which water moves downglacier quickly. Negative lags are often driven by distributed surface runoff supply but could also reflect up-glacier migration of pressure pulses in a connected but poorly draining system. Finally, negative lags that persist for longer periods of time during the melt season are most likely driven by surface melt reaching the subglacial environment in a distributed manner.

We note that the time-lags are less varied and more consistently positive during the surge. Following our interpretation of the time-lags outlined above, this implies that the drainage system is more frequently connected and downglacier motion of water or pressure pulses is more prevalent during the surge than after termination. Nevertheless, the disruption in coherence during the late winter of 2021 shows that there is a continued seasonal evolution, including a poorly connected phase, in the drainage system during the surge.

5.5 Relation between remote sensing time-series and fjord conditions

The seasonal fluctuation in the Sentinel-3 and Sentinel-2 signals (Fig. 3) provides confidence that the longer term changes in the recorded radiance and reflectance are driven by the release of meltwater (e.g. McGrath and others, 2010). The Sentinel-1 signal reflects the back-scattering at the surface, which is largely driven by the presence of small icebergs at the water surface (Benn and others, 2019). Our observations reflect this as the amount of back-scattering strongly increases in October of 2020. This coincides with the time at which Sít’ Kusá has covered its sediment shoal and is advancing into Disenchantment Bay, promoting more consistent calving (Liu and others, 2024). The Sentinel-2 signal shows a similar increase in late 2020 to levels similar to those recorded for the permanently iceberg-filled water in front of Sít’ Tlein (Hubbard Glacier), suggesting that the surface reflectance is also affected by the increase in icebergs in the bay. The Sentinel-3 signal seems less affected by iceberg presence, which we suspect might be due to a combination of higher spectral resolution and lower spatial resolution relative to Sentinel-2.

The Sentinel-1 back-scattering intensity is maximal early in the melt season and decreases during the summer as the turbidity derived from Sentinel-3 increases (Fig. 3), consistent with the suggestion that increased subglacial discharge pushes away icebergs from the area of interest. This pattern gives us confidence that the Sentinel-3 OLCI band 7 time-series capture changes in water turbidity in front of the Sít’ Kusá terminus and is a plausible proxy observation for the relative intensity of subglacial water release at the terminus. While our observations do not provide a definitive measure of sub-aqueous discharge volumes, they are mutually consistent and we have good confidence that they reflect the relative variability in sub-aqueous frontal discharge. The close proximity of the tidewater terminus of Sít’ Tlein likely affects our observations to some extent, as shown by seasonality in the inferred iceberg production before the Sít’ Kusá terminus reaches the ocean. A similar analysis for an area chosen to capture discharge from the Sít’ Tlein terminus shows commonalities with the data recorded for the Sít’ Kusá terminus but also notable differences (Fig. S6.1). We have chosen the area of interest as far removed from the Sít’ Tlein terminus as possible (Figs 1, 3a) and thus expect the signals to be dominated by changes at the Sít’ Kusá terminus.

6. Discussion: signal interpretation

In this section, we aim to disentangle the behavior of the subglacial drainage system and how it affects ice dynamics through the interpretation of our various time-series. Each subsection is titled with an assertion that we subsequently support by our observations.

6.1 An efficient component of the subglacial drainage system exists intermittently prior to surge termination

Our observations show that tremor-generating glaciohydraulic sources exist prior to surge termination at least at several locations along the main trunk of Sít’ Kusá (Figs S1.1, S5.5): we find tremor during the mid-surge summer that has similar power to that of the quiescent summer and the tremor power of the mid-surge winter is within several dB of summer levels (cf. quiescent winter tremor which is ~15 dB lower than mid-surge winter). Additionally, the tremor spectral pattern remains consistent between hydraulically active summers and the active phase winter (Fig. 5), suggesting similar source mechanisms. Previous work widely attributes such tremor to turbulent flow and sediment transport through a channelized subglacial drainage system (e.g. Bartholomaus and others, 2015; Gimbert and others, 2016; Nanni and others, 2020; Lindner and others, 2020). Furthermore, Nanni and others (2021) inferred that inefficient, broadly distributed, linked cavities can also produce detectable hydraulic tremor, albeit with power 20 dB less than that of mid-summer (Gimbert and others, 2021). However, the mid-surge 2020–2021 winter tremor we record remains within ~5 dB of peak power (Fig. 4a), which is a much smaller gap than the change in tremor power associated with the transition from linked cavities to channelized drainage on Argentiére Glacier (Nanni and others, 2020, 2021).

These tremor-generating locations have hydrologic connections to the surface (Fig. 4a), as the tremor power consistently shows spikes that coincide with spikes in modeled surface runoff (Fig. 4a). Additionally, these tremor-generating locations are connected in the along flow direction, as we can detect lagged coherence between the tremor signal observed at various seismic stations (Fig. 4b). We do not know the precise location of the tremor sources and thus cannot infer the exact travel distance of the water velocity pulse through the subglacial drainage system. Nevertheless, the source locations are almost certainly within ±3 km of the centerline distance between station pairs, meaning ~8.2±3 km for SE15–SE7 and 5±3 km for SE7–SW2. Despite the large uncertainties, these distances are within the same order of magnitude as those used for dye tracing experiments on Variegated Glacier in Kamb and others (1985) (8 and 10 km, their Fig. 11). Our time-lags are universally below 25 h; and the median positive time-lag for SE15–SE7 is 8.5 h (Fig. 4b) over the whole record and 9.4 h when computed for just the surge phase (prior to 15 August 2021). These values are relatively close to those that Kamb and others (1985) attribute to post-surge efficient drainage (main dye concentration peak after ~4 h) and shorter than those found for the surge phase linked cavity system on Variegated Glacier (first dye concentration peak after ~50 h).

The upper range of the time-lags we observe (20–25 h) and the detailed inspection of the lags in Figure 7 show that there are moments of inefficient, or disrupted drainage during the surge, notably from 1 March to 15 April 2021 and during July 2021 (Fig. 4b). However, considerable time periods of the surge with short time-lags and the absence of changes in the tremor frequency content lead us to suggest that there is at least intermittent efficient water transport occurring in an along flow component of the drainage system prior to surge termination.

This component is frequently present throughout the lower trunk of the glacier from September 2020 to April 2022 and in the upper trunk at least from September 2020 to October 2021. However, we cannot fully ascertain how prevalent the efficient components are spatially. It is possible and perhaps even likely that the tremor-generating components of the subglacial drainage system do not coincide spatially with the area(s) of the glacier base regulating ice velocities or surge propagation.

The absence of coherent tremor variations in the upper trunk between October 2021 and August 2022, and between April 2022 and August 2022 suggests that the subglacial drainage system is gradually disconnecting during the winter after surge termination. This is in line with the lower tremor amplitude during that period, as lower basal water volumes and slower water flow would allow for creep closure to gradually close off waterways and reduce the connectivity of the drainage system (e.g. Hart and others, 2022).

6.2 Subglacial water-flow continues during the surge winter

A comparison of the power spectral density probability for the winter (December–March) and summer seasons (May–August) (Fig. 5) shows that the median spectral power across the glaciohydraulic tremor frequency range is similar to but slightly lower during the post surge melt season compared to the surge melt season (Fig. 5b). Both spectra have a similar shape, suggesting that there is a stable and consistent seismic source process (Gimbert and others, 2014) and further hinting that the subglacial hydrological systems during the summer of 2021 and 2022 are comparable in behavior. Meanwhile, median spectral power is considerably lower during the post surge winter compared to the active phase winter (Fig. 5a). A second notable difference between the two winters in our record lies in the time-lags. The surge winter is marked by relatively consistent positive lags that reflect a connected drainage system where most seismic noise is produced by water or pressure pulses that move downglacier (Figs 4b, 7). The 2021–2022 post surge winter sees less consistent coherence and much more variable lags (Fig. 4b). This behavior points to a less connected system where distributed inputs from surface melt are relatively more important (Fig. 7).

Together, our interpretation of the time-lags and the higher noise levels suggest that the drainage system during the 2020–2021 winter is ‘higher volume’ than during the 2021–2022 winter. This winter volume is likely considerably lower than melt-season values, as there are disruptions in drainage during late winter 2021 (Fig. 4b, Text S5) and bay turbidity decreases during the fall of 2020, in a similar pattern to other years (Fig. 4d). Nevertheless, some continued base level of water availability during the winter would allow the drainage system to remain more connected and explain the higher tremor levels observed in 2020–2021. In the absence of water supply from surface runoff, water might be sourced from sub- or englacial reservoirs. Such overwinter storage likely occurs annually on Sít’ Kusá (Liu and others, 2024) and has been suggested to play a role in the surge mechanism on other Alaskan glaciers (Humphrey and Raymond, 1994; Lingle and Fatland, 2003; Barrett and others, 2008; Zhan, 2019; Hart and others, 2022). Additionally, water could be sourced through strain heating and basal melt (Benn and others, 2019). Our interpretation of a ‘high water-volume’ drainage system during the surge and ‘low water volume’ drainage system post-surge is also consistent with the enthalpy-based model of surging, where surge termination is driven by the draining of the glacier base (Benn and others, 2019, 2022).

6.3 Indistinguishable change in frontal discharge before and after the surge

The seasonal cycle contained in the Sentinel-3 record seems largely unchanged throughout Sít’ Kusá's surge cycle (Figs 3, 4e). The Sentinel-2 surface reflectance does see a spike starting in October 2020, but as discussed earlier its similarity to the Sentinel-1 signal suggests that it is driven by iceberg presence rather than water turbidity (Fig. 3). As such, our indirect observations do not point to a significant disruption in frontal discharge during the surge build up or during the active phase and while the seasonally averaged turbidity in the bay is slightly higher in 2021 relative to 2019, we are unable to observe a sudden single abnormally high-volume discharge event during the 2021 melt season (Fig. 4d).

Previous work frequently notes retention of water below the glacier during the active phase (Clarke and others, 1984; Lingle and Fatland, 2003), while termination coincides with release of large volumes of water from the subglacial environment (e.g. Kamb and others, 1985; Benn and others, 2019). We are unaware of descriptions of prior Alaskan surge terminations specifically noting an absence of abrupt water releases from the terminus.

While it is possible that our proxy record simply missed a spike due to cloudy conditions, the bay turbidity record is consistent with our earlier interpretations of the drainage system. Each year, turbidity in front of the terminus starts increasing before the spring onset of surface runoff (Figs 3, 4d). This early increase might be associated with the frontal release of water stored overwinter, linked to the annual early spring speedups described in Liu and others (2024). The apparent lack of disruption in frontal discharge during the surge is also consistent with our interpretation that efficient drainage occurs intermittently during the surge. The absence of a single major water release event could suggest a rather gradual termination of the surge on Sít’ Kusá, which we discuss further in the next section.

6.4 Short-term velocity fluctuations modulated by subglacial drainage overlay gradual surge termination

Here we focus on the various time-series related to hydrology for the 2021 melt season that marks surge termination (Fig. 8). Lags in the upper trunk decrease from 20 h in early June to alternating slightly negative and positive 5 h lags, which we have interpreted earlier to reflect an intermittently connected and efficient drainage system component with significant contributions from distributed surface melt. From mid-June to early-July, the tremor pulse velocities in the lower trunk gradually decrease from >5 to ~1 h. The lag time decrease coincides with a gradual decrease in ice surface velocities from >15 to ~10 m d⁻¹. These changes hint at a shift toward a more efficient subglacial drainage system in the lower trunk and resemble early summer slowdowns due to hydrological changes observed during quiescent years on Sít’ Kusá (Liu and others, 2024) and on other Alaskan surge-type glaciers during quiescence (Abe and Furuya, 2015). During this time, correlations between the modeled surface runoff and the bay turbidity record are generally around r = 0.5, implying some component of the surface runoff travels through the glacier quickly and affects the bay turbidity levels. From early July to early August 2021, there are disruptions in coherence in the upper trunk and lags in the lower trunk remain largely below zero, suggesting a poorly connected drainage system where tremor is generated largely by spatially distributed surface runoff. The correlation between modeled surface runoff and bay turbidity largely disappears during this time, with two oscillations with ~10 d periods in the turbidity record that are of similar amplitude to earlier variations but appear unforced by surface runoff on the glacier. Meanwhile, there are large variations in daily GNSS velocities (between 2 and 10 m d⁻¹) from 1 July to 11 August. The periods of speedup coincide with increases in meltwater supply and in seismic tremor. Broadly, our data seem to reflect a disrupted and relatively low efficiency system during the final weeks of the surge, with variability in ice velocities driven by changes in water supply from distributed surface runoff. Meanwhile, the proxy record for sub-aqueous frontal discharge suggests that considerable volumes of water are gradually released from sub- or englacial storage.

Figure 8. Glacier behavior during the 2021 melt season. (a) Modeled surface runoff and average daily water surface radiance measured with Sentinel 3 OLCI. Correlation coefficients calculated between 10 d rolling windows with a 12 h shift, of both time-series interpolated on a 12 h interval. (b) Median seismic power measured at SE7 and lags between tremor power time-series measured between SE15–SE7 (purple dots) and SE7–SW2 (red dots) station pairs. (c) Glacier surface velocity measured with GPS and satellite imagery feature tracking.

In early August, more sustained high coherence returns between the stations framing the lower trunk, with lags of ~5 h. From 1 August to 11 August, ice velocities undergo a decreasing trend from ~10 to ~3 m d⁻¹ and the turbidity in front of the terminus is gradually increasing, without significant increases in surface runoff. A rainstorm from 11 August to 14 August 2021 drives extreme surface runoff, which coincides with a peak in glaciohydraulic tremor and in ice velocities. The rising limbs in the runoff and tremor signals coincide with a peak in ice velocities of ~16 m d⁻¹, which is directly followed by a slowdown to velocities below 5 m d⁻¹. The spike in surface runoff, glaciohydraulic tremor and ice velocity fits well with previous suggestions of the surge termination process reported in Kamb and others (1985) and Kamb (1987), where a switch from a linked cavity system to a channelized drainage system drives surge termination. However, we do not observe any change in the time-lags between tremor time-series for either the upper or lower trunk. This is remarkable as the time-lags reliably reflect changes in the drainage system at other moments in our record and it seems likely that we would observe some change if there was a sudden and widespread switching between drainage system configurations from 11 to 14 August. Additionally, we do not observe a sudden spike in turbidity that would reflect a sudden release of large water volumes. The latter might be explained by more gradual prior water release during July 2021, which rather echoes behavior observed on Svalbard (Murray and others, 2000). Such a gradual water release, combined with the slowdown in ice velocities from 1 to 11 August and the 5 h time-lags existing from 1 August on-wards, suggest that any switching in drainage systems leading to surge termination was more gradual on Sít’ Kusá than on Variegated Glacier. We suggest some component of efficient drainage was established by 1 August in the lower trunk. Subsequently, the rainstorm-driven influx of water overwhelmed this existing efficient drainage and led to widespread high basal water pressure and low effective pressure, driving the peak in ice velocities. The 12–13 August speedup is superimposed on a trend of gradual slowdown and seems modulated by temporary overloading of the efficient component of the subglacial drainage system. Such a mechanism is common and well explained on non-surge-type glaciers (e.g. Iken and Bindschadler, 1986; Anderson and others, 2004; Bartholomaus and others, 2008; Schoof, 2010; Labedz and others, 2022; Hart and others, 2022).

During the remainder of the melt season, ice velocities are generally lower and less variable (Fig. 8c). Peaks in surface runoff result in sharp increases in glaciohydraulic tremor with muted responses in ice velocity variations (Fig. 8). Such behavior closely resembles that observed during the late-melt season on non-surging alpine glaciers (Nanni and others, 2020; Labedz and others, 2022). Time lags for the lower trunk remain around 5 h until 26 August, then decrease to become slightly negative by 31 August (Fig. 8b). Finally, the surface runoff and bay turbidity signals are highly coherent and vary in phase, as shown by the consistently high correlation between the two time-series (Fig. 8a). We suggest this points to a drainage system in the lower trunk through which water circulates quickly and where surface runoff is the dominant supply. The now ‘low volume’ system that lacks significant input from sub- or englacial water contributions might be gradually closing under ice overburden pressure by late August 2021, as evidenced by the lost coherence in the upper trunk and slightly negative lags in the lower trunk.

Our data point to gradual and non-monotonic changes in the drainage system during the 2021 melt season that combine toward surge termination. The glacier gradually releases water throughout the summer and the 11–14 August rainstorm seems to coincide with the final release and emptying of any sub- or englacial reserves. Throughout the summer, changes in surface runoff drive changes in ice velocity on the timescale of days, which are superimposed on the trend of surge slowdown. Such variability superimposed on longer trends is not exclusive to Sít’ Kusá (Benn and others, 2022). This variability shows the drainage system is evolving simultaneously at multiple time-scales and conflicts with the notion of a single switch marking termination.

7. Synthesis

Our observations are broadly consistent with existing theoretical mechanisms of glacier surging, but draw attention to some key points concerning the evolution of the subglacial drainage system during the surge. On Sít’ Kusá there is a component of the subglacial hydrological system that extends along the main glacier trunk, is connected to the glacier surface and is intermittently efficient during the active phase. This component is connected to a large enough area of the glacier base for its evolution to occasionally force high amplitude variability in ice velocities during the active phase. Similar high variability in surface velocities is observed during other surges when the temporal resolution of observations is high enough (Beaud and others, 2022; Benn and others, 2022) and has been linked to variability in basal water pressure (Kamb and others, 1985). Our data show that this variability can occur as suggested by Kamb and others (1985) during periods of disrupted, inefficient drainage (e.g. July 2021) but also through the overwhelming of existing efficient, channelized drainage (e.g. October 2020, August 2021). The second overwhelming mechanism has been observed repeatedly on non-surging glaciers (e.g. Anderson and others, 2004; Bartholomew and others, 2012; Cowton and others, 2013) and its presence during an active surge shows there are consistent mechanisms driving flow velocities throughout the surge cycle.

On the timescales of surge evolution, we observe variability in the efficiency of the subglacial system and possibly short-term changes in its configuration but no change in how it fundamentally behaves or clear evidence of a single widespread switching between configurations. Surge termination during the 2021 summer is a gradual process rather than a sudden switch in behavior. This echoes conclusions drawn for the slow surging Trapridge Glacier by Frappé and Clarke (2007), where authors hint at an out of equilibrium but fundamentally unchanged drainage system where subsequent changes in efficiency accumulate into the observed surge behavior. Benn and others (2022) suggest surge-type glaciers might be characterized by a basal water surplus that is too high to be accommodated by a ‘slow’ system and too low to transition to a ‘fast’ system. Reciprocally, such a situation would likely lead to the semi-efficient, or intermittently efficient, drainage system we seem to observe.

Finally, we note that there seems to be a shift in the availability of sub- or en-glacially stored water that occurs with surge termination. We suggest the surge-time drainage system differentiates itself primarily by an overabundance of available water that is independent of surface runoff, as shown by the continued water flow during the 2020–2021 winter that contrasts with the 2021–2022 winter. This notion is in line with the observation that overwinter storage of surface runoff allows springtime speedups on Sít’ Kusá (Liu and others, 2024), as well as with earlier work pointing to a potential role for gradual water accumulation in driving glacier surging (Humphrey and Raymond, 1994; Lingle and Fatland, 2003; Abe and Furuya, 2015). Furthermore, our observations align with the notion that a basal enthalpy surplus, which translates to a basal water surplus on temperate glaciers, accumulates during quiescence and surge onset. This basal water surplus then gradually dissipates during the later stages of the surge (Benn and others, 2019, 2022).

8. Conclusion

Although observing the subglacial environment remains a challenging task, the time-series presented in this study provide a window into the evolution of subglacial drainage over time, without direct access to the glacier bed. In particular, we note the novel application of wavelet coherence analysis to leverage phase lags in median seismic power toward observing subglacial water and pressure pulse migration. While there are limitations outlined in this study, the approach allows for continuous and remotely sensed monitoring of changes in the critical (seismically loud) components of the subglacial drainage system.

Our observations show that glacier surges can be resilient to continuous, observable and intermittently efficient drainage. The channelized components of subglacial drainage can intermittently restrict subglacial water flow, modulating surge dynamics. These observations conflict with a theory of ‘hard switching’ between fully efficient and inefficient drainage systems (Kamb, 1987). Nevertheless, they strengthen the broader applicability of the hydrologically regulated surge mechanism as they provide avenues through which a drainage system that continues to undergo seasonal and sub-seasonal changes can drive multi-year surge dynamics (Benn and others, 2022).

The evolution of the subglacial drainage system during the observed part of the surge cycle seems to express itself on a spectrum of efficiency. It underlines complexity and variability on sub-seasonal, seasonal and multi-year timescales that interfere to produce the spectacular glacier dynamics on Sít’ Kusá. We did not observe hydrological features that set Sít’ Kusá apart from glaciers with steady-state flow behavior. It seems worthwhile to consider whether similar multi-annual velocity variations might be present at all glaciers, just with lower amplitudes. Such a perspective seems consistent with the emerging notion that many glaciers not identified strictly as ‘surging’ have complex multi-year velocity patterns (Herreid and Truffer, 2016). This perhaps suggests that glacier surges are simply the most spectacular, and easiest to detect, expressions of hydrologically driven periodic velocity variations common to many glaciers.