Introduction
Fast glacier flow through basal sliding is ubiquitously linked to the presence of liquid water at the glacier bed (e.g. Iken and Bindschadler, Reference Iken and Bindschadler1986; Meier and others, Reference Meier1994; Truffer and Echelmeyer, Reference Truffer and Echelmeyer2003). The relationship between sliding velocity and basal water pressure has been demonstrated for transient speedup events on timescales of days (e.g. Kamb and Engelhardt, Reference Kamb and Engelhardt1987; Kamb and others, Reference Kamb, Engelhardt, Fahnestock, Humphrey, Meier and Stone1994; Nienow and others, Reference Nienow2005; Bartholomaus and others, Reference Bartholomaus, Anderson and Anderson2011; Bartholomew and others, Reference Bartholomew, Nienow, Sole, Mair, Cowton and King2012) and for seasonal velocity patterns (e.g. Moon and others, Reference Moon2014; Nanni and others, Reference Nanni, Scherler, Ayoub, Millan, Herman and Avouac2023). High basal water pressure and water volume at the glacier bed is also a central component of the community’s understanding of glacier surging (Kamb and others, Reference Kamb1985; Fowler, Reference Fowler1987; Björnsson, Reference Björnsson1998; Lingle and Fatland, Reference Lingle and Fatland2003; Murray and others, Reference Murray, Strozzi, Luckman, Jiskoot and Christakos2003; Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a), but the exact role of basal water volume in driving multi-year oscillations in ice velocity remains unobserved and poorly understood.
Glacier surges are semi-periodic, high-amplitude (2–3 orders of magnitude) accelerations in ice velocity that redistribute ice from higher to lower elevations (Meier and Post, Reference Meier and Post1969). Surge phases alternate with longer quiescent phases, during which ice velocities are low and ice flow is insufficient to balance the climatic mass balance gradient. Most explanations for glacier surging suggest some form of disruption of the basal drainage system to cause a surplus in basal water and a decrease in effective pressure, which is necessary for fast flow. Earlier theories tend to place this disruption during the surge phase: examples of surge drivers include the retention of basal water above a barrier of cold-based ice (Fowler and others, Reference Fowler, Murray and Ng2001; Murray and others, Reference Murray2000) or a linked cavity system that is resistant to growth of its orifices (Kamb, Reference Kamb1987). However, the recent advent of high-quality satellite-based ice velocity time-series has shown that seasonal velocity patterns often overlay the longer-term surge pattern (Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022; Benn and others, Reference Benn, Hewitt and Luckman2022; Liu and others, Reference Liu, Gendreau, Enderlin and Aberle2024b). Late summer slowdown, present in many surge velocity records, implies that basal water pressures can decrease during the surge and that basal water at least periodically escapes the system (e.g. Sundal and others, Reference Sundal, Shepherd, Nienow, Hanna, Palmer and Huybrechts2011). Field observations also suggest that channelized drainage can exist during surges (Terleth and others, Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024b). Newer models of glacier surging place a drainage disruption in the quiescent phase rather than in the surge phase (Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a, Reference Benn, Hewitt and Luckman2022). These models allow basal water to accumulate (and effective pressure to decrease) over multiple years, with positive feedback loops driving surge initiation. In this framework, the disruption in basal drainage does not need to be complete as long as there is a net annual gain in basal water during the quiescent phase and the response to change in water supply depends on the current volume of basal water (Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a; Terleth and others, Reference Terleth, Bartholomaus, Enderlin, Mikesell and Liu2024a).
The idea that subglacial water accumulates within the glacier implies that there exists some form of perennial englacial or subglacial water reservoir (e.g. Stenborg, Reference Stenborg1970; Tangborn and others, Reference Tangborn, Krimmel and Meier1975; Meier and others, Reference Meier1994; Murray and others, Reference Murray2000; Jansson and others, Reference Jansson, Hock and Schneider2003). Many Alaskan glaciers undergo downglacier propagating speed-ups during winter or early spring, which are likely driven by water stored overwinter (Abe and Furuya, Reference Abe and Furuya2015; Burgess and others, Reference Burgess, Forster and Larsen2013; Liu and others, Reference Liu, Gendreau, Enderlin and Aberle2024b). Basal water storage has also been suggested to play a role in the context of surging. There were large volumes of water stored in the reservoir zone of Variegated Glacier prior to its 1982–83 surge, and these high volumes of basal water have been hypothesized to be a required condition for sustaining surge velocities (Humphrey and Raymond, Reference Humphrey and Raymond1994). Similarly, Lingle and Fatland (Reference Lingle and Fatland2003) propose a mechanism where the transfer of englacial water to the glacier bed during winter drives winter surge initiations in Alaska. Finally, Zhan (Reference Zhan2019) presents observational evidence of an increase in basal water volume held in transverse basal crevasses during a surge of Bering Glacier. However, it is uncertain whether this water is a driver or a consequence of the surge on Bering Glacier, because no multi-year changes in water storage were apparent during the three years prior to surge onset.
The evolution of sub- and englacial water storage through the different stages of the surge cycle, and how it interacts with ice dynamics, remains unclear. Gradual water accumulation and multi-year hysteresis in the drainage system, predicted in the enthalpy model, have not yet been observed. Here, we use seismic noise interferometry to investigate if, when, and where basal water accumulates and diminishes during the 2020–21 surge and subsequent quiescent phase at Sít’ Kusá glacier in Alaska. We discuss possible mechanisms that could explain the periodicity in basal water volume fluctuations, its role in driving surge velocities, and basal water sourcing.
Methods
Field site and instrumentation
Sít’ Kusá (also briefly referred to as Turner Glacier) is a 
$\sim$30 km long and 
$\sim$2 km wide glacier located at the head of Disenchantment Bay, on Tlingit land in Wrangell-St. Elias National Park. Sít’ Kusá has the world’s shortest observed surge cycle, with a 
$\sim$1.5 year surge phase repeating every 
$\sim$5-8 years (Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021). The lowermost 
$\sim$7 kilometers of the main trunk receive mass during surge phases, thickening up to 
$\sim$100 meters (receiving zone in Fig. 1). The reservoir area extends from 7 km above the terminus to the headwall of a southern tributary and to the bottom of an icefall 19.5 km from the terminus in a northern tributary Nolan and others (Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021), reservoir zone in Fig. 1). During quiescent phases, glacier speedups (from 0 to 
$\sim$10 m/day) initiate in November-December in the northern tributary (initiation zone in Fig. 1) and propagate downglacier (Liu and others, Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a). Each subsequent year of the quiescent phase, these speedups propagate a further distance toward the terminus, until one speedup escalates into a full-scale surge (see Liu and others (Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a) for a more detailed description). To observe changes in glacier hydrology, we maintained a network of broadband seismometers around Sít’ Kusá from September 2020 to September 2025. These in situ observations captured the last of a spring 2020-summer 2021 surge, a 3.5 yr quiescent period and a major speedup (or minor surge onset; see discussion section) in spring 2025. Here, we focus on analysis from nine 3-component seismic stations with relatively complete data recordings (Fig. 1). Further details on the data acquisition strategy are found in Supplementary Table S.1.
Map showing the locations of the seismic stations. Contours (50 m) are from the US3DEP dataset. Glacier extent from RGI database and kilometer markers from Liu and others (Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a). Background imagery is a Sentinel-2 cloud-free natural color median composite acquired between 1 July and 30 September 2023. The inset map shows the Sít’ Kusá location in Alaska. Coordinate reference system (CRS): UTM 7N, EPSG:32606. The datum is WGS 84.
Figure 1 Long description
Glacier extent shown in a purple overlay of the glacier, the surge initiation zone in blue, and the receiving zone in green. Along with elevation contours. There are blue markers showing distance from the glacier terminus along the centerline each 2km. All the locations of the seismic stations are shown. There are four seismic stations on the south west side of the terminus called SW and a number reflecting km from the terminus, and five stations on the north east side called SE and a number reflecting km from the terminus.
Ambient noise analysis
Changes in seismic surface wave velocities can occur due to changes in the medium through which they propagate (e.g. Shapiro and Campillo, Reference Shapiro and Campillo2004; Curtis and others, Reference Curtis, Gerstoft, Sato, Snieder and Wapenaar2006). Passive seismic interferometry, or ambient noise correlation, uses the time correlation between seismic signals recorded simultaneously at two different locations to estimate the Green’s function (i.e. a description of how the ground moves between two stations as the result of a seismic impulse) under the assumption that the seismic noise sources are dispersed randomly around the recording locations (e.g. Curtis and others, Reference Curtis, Gerstoft, Sato, Snieder and Wapenaar2006; Larose and others, Reference Larose2015). Using the Python package seisgo (Jiang and Denolle, Reference Jiang and Denolle2020; Yang and others, Reference Yang, Bryan, Okubo, Jiang, Clements and Denolle2023), we pre-processed each channel of the raw 3-component (HHZ, HHE and HHN) seismic data by de-meaning, de-trending, whitening in the frequency domain and applying a 0.01–5 Hz band-pass filter. We then cross-correlated each receiver pair in 2-hour windows with 1-hour overlap and stacked the correlation functions within 72-hour segments (Fig. 2a,b). These stacks were corrected for clock drift following Hable and others (Reference Hable, Sigloch, Barruol, Stähler and Hadziioannou2018) (Supplementary Fig. S.2). Clock-drift corrections were primarily needed for the SE15-SW14 station pair, where they were applied to 
$\sim$22% of the data. We then computed the 
$dv/v$ values (i.e. the relative change in seismic wave velocity) within a 20-second window that started 6 seconds after the main surface wave arrival (
$\sim$9 seconds from the center of the correlation function) on both the causal and the a-causal side of the correlation function, so as to isolate the coda waves. We estimated 
$dv/v$ using the stretching method with 72-hour mean stacks relative to a reference correlation function derived from the mean correlation function of the entire data record (e.g. Curtis and others, (Reference Curtis, Gerstoft, Sato, Snieder and Wapenaar2006; Fig. 2)). The relative time shift (
$dt$) necessary to achieve the highest correlation coefficient between the two correlation functions yields the relative velocity change (assuming 
$dt/t=-dv/v$) (Guillemot and others, Reference Guillemot2020). We report 
$dv/v$ as a percentage and the unitless maximum correlation between each 72-hour segment and the reference correlogram (Fig. 2c). Because surface waves of different frequencies sample different depths, the frequency band can be chosen to target changes specifically near the glacier-bed interface (e.g. Mordret and others, Reference Mordret, Mikesell, Harig, Lipovsky and Prieto2016). Following Zhan (Reference Zhan2019), we consider changes in 
$dv/v$ within the 0.7–1.5 Hz frequency range to focus our analysis on that portion of the Green’s function with maximum sensitivity to the till layer at the bed of Sít’ Kusá in an effort to specifically monitor changes in water volume near the glacier bed (see supplementary information S.3). The 
$dv/v$ values for the three channels are averaged for each point in time to create one timeseries for each station pair (see supplementary information S.4).
(a) Time-dependent cross-correlation between 2-hour segments of pre-processed seismic data recorded simultaneously at SW7 and SE9, and band-pass filtered between 0.01 and 5 Hz. Green shading indicates time windows within which we compute the 
$dv/v$ values. (b) Reference correlation function obtained by median stacking all the cross-correlation functions. (c) Unitless maximum correlation value obtained through the stretching method between each 72-hour segment and the reference correlation stack (black line) and 
$dv/v$ values as percentage obtained from the time shift for which the maximum correlation occurred (blue dots). More opaque markers indicate stronger correlation values.
Changes in the pore pressure within the loose till and more consolidated sediment layers can occur as a result of changes in the overburden pressure driven by loading from seasonal snow and/or changes in glacier thickness (e.g. Wang and others, Reference Wang, Brenguier, Campillo, Lecointre, Takeda and Aoki2017; Guillemot and others, Reference Guillemot2020; Donaldson and others, Reference Donaldson, Winder, Caudron and White2019; Mordret and others, Reference Mordret, Mikesell, Harig, Lipovsky and Prieto2016). We follow methods outlined in Wang and others (Reference Wang, Brenguier, Campillo, Lecointre, Takeda and Aoki2017) to account for the variability in 
$dv/v$ driven by seasonal snow loading. We use a spatially averaged estimate of snow thickness modeled with the Energy Balance Firn Model (EBFM - van Pelt and Oerlemans (Reference van Pelt and Oerlemans2012)) forced with ERA-5 re-analysis data (Hersbach and others, Reference Hersbach2020) to construct a linear relationship with the 
$dv/v$ signal recorded at a set of cross-glacier and glacier-adjacent station pairs. Validation of the modeled snow thickness is described in the following section. The slope and intercept of this linear model are used to construct a synthetic 
$dv/v$ estimate based on the snow load on the glacier between each cross-glacier station pair. This estimate aims to represent the contribution of snow loading to the 
$dv/v$ signal, and is subtracted from the observed 
$dv/v$ values (see supplementary information S.5).
Additional timeseries
We complement the 
$dv/v$ measurements with additional timeseries that give insight into the evolution of the subglacial drainage system. We quantify the amplitude of glaciohydraulic tremor (GHT, the seismic tremor produced by fast subglacial water flow) as the median power spectral density between 3 and 10 Hz in 1-hour windows. The processing of this signal is described in Bartholomaus and others (Reference Bartholomaus, Amundson, Walter, O’Neel, West and Larsen2015); the timeseries of GHT presented here are a temporal extension of those presented in Terleth and others (Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024b).
We compute an estimate of the glacier surface runoff and of snow water equivalent on the glacier with the EBFM surface energy balance model. We force the model with climatic data down-scaled from the ERA-5 reanalysis dataset (Hersbach and others, Reference Hersbach2020). This approach was previously applied and validated through a comparison with observations of sub-seasonal surface elevation changes in Terleth and others (Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024b). Finally, we derive ice-velocity estimates by analysis of satellite image pairs with CautoRIFT (Liu and others, Reference Liu, Gendreau, Enderlin and Aberle2024b). CautoRIFT relies on the cross correlation of features in panchromatic image timeseries from Landsat 8 and 9, and the near infrared band of Sentinel-2 images to estimate ice velocities. CautoRIFT allows customization of cross correlation chip sizes as well as cross-platform imagery analysis, yielding estimates with high spatio-temporal resolution. The velocities we show here are a temporal extension of those presented in Liu and others (Reference Liu, Gendreau, Enderlin and Aberle2024b).
Results
Receiving zone - km 3
In the receiving zone, ice velocities are fast during the surge (
$ \gt $25 m/d, Liu and others (Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a)) and become near stagnant (
$ \lt $ 1 m/d) after October 2021 (Fig. 3f, i.e. lower panels). GHT amplitudes recorded at SW2 increase over the course of the 2020–21 winter to 
$ \gt 3.5 \times 10^{-8}$ m/s, then fall rapidly during the latter half of the 2021 melt-season (Fig. 3e). Low-amplitude seasonal variations in GHT persist throughout the rest of the record, with summer maxima of 
$\sim$1
$ \times 10^{-8}$ m/s and winter minima of 
$\sim$0.5
$ \times 10^{-8}$ m/s (Fig. 3e). The lower reaches of our seismic instrument network were plagued by animal encounters (Tape and others, Reference Tape, Heath, Baker, Dalton, Aderhold and West2019) and were only maintained until June 2024, leading to somewhat sparse observations of the receiving area.
Reservoir zone - km 8
The region 8 km from the terminus is part of the very lower end of the reservoir zone, with the dynamic balance line located 
$\sim$7 km from the terminus (Fig. 1; Liu and others (Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a)). Ice velocities peak at 
$\sim$25 m/d during the surge in May 2021 (approximately 15 months after surge initiation) and decrease to nearly-stagnant (
$ \lt $1 m/d) by late October 2021. Velocities remain nearly-stagnant until at least early April 2025, when they start increasing to reach 
$\sim$15 m/d by June 2025 coincident with the arrival of the 2025 speedup at this location. After peaking in June, ice velocities decrease to 
$\sim$10 m/day by August 2025 (Fig. 3d).
Evolution of Sít’ Kusá’s hydrology in three different sections of the glacier. (a),(c),(e) Show 
$dv/v$ values (i.e. the relative change in seismic wave velocity) on the left-hand y-axis, with the y-axis flipped so that slower seismic waves, indicating higher subglacial water content, plot negative up. Transparency of the markers is scaled with the correlation coefficient associated with the 
$dv/v$ measurement. The dotted 
$dv/v$ trend lines are 30-day rolling means of the plotted 
$dv/v$ values. The rolling means include both station pairs when two are plotted with overlapping timestamps (panels c and e). Right-hand y-axes show modeled surface runoff and glaciohydraulic tremor. (b),(d),(f) Show the evolution of ice velocities at each kilometer marker, estimated with CautoRIFT.
GHT amplitude recorded at SE9 is relatively steady 
$\sim$1.3
$ \times 10^{-8}$ m/s during the 2020–21 winter, and starts increasing in early May 2021, slightly prior to the onset of the melt-season (Fig. 3c). The SE9 GHT record then rises to 
$\sim$3 
$\times 10^{-8}$ m/s by July 2021. This increase relative to winter levels is higher than the increase relative to winter levels observed, further up-glacier, at SE15. Throughout the rest of the record, average winter and summer GHT amplitudes are consistent at SE9 (
$\sim$0.3
$ \times 10^{-8}$ m/s and 
$\sim$3
$ \times 10^{-8}$ m/s, respectively). There is no change in GHT amplitude or timing between 2022–23 and 2024–25. The increase in GHT at SE9 during the spring of 2025 seems to largely follow the melt signal. We record GHT values of up to 
$\sim$2.4
$ \times 10^{-8}$ m/s in July 2025, which is only slightly higher than values recorded in July of previous years (
$\sim$2.1
$ \times 10^{-8}$ m/s on average), but the record ends in July 2025 and we cannot determine whether these elevated values are a short term transient or a longer trend.
We calculate the relative change in seismic wave velocity (
$dv/v$) values for both SW7-SE9 and SW8-SE9 station pairs, to more robustly interpret changes within the same region of the glacier and fill in data outages. Because they exhibit similar behavior, we describe them jointly. Both timeseries are relatively stable during the 2020–21 surge winter and through the 2021 melt season, then increase during the 2021–22 fall and winter. Over the next three years of record, both station pairs show an annual pattern of rapid 
$dv/v$ decrease during the late melt-season (August-October) with an increase during winter and spring (November-July). Long-term trends are not apparent in the 
$dv/v$ record after April 2022, and we do not discern any changes in the partial 
$dv/v$ record preceding the onset of the 2025 speedup in ice velocities at km 8.
Reservoir zone - km 15
Fifteen kilometers from the terminus (km 15) is approximately the middle of the reservoir zone (Fig. 1). Ice velocities in this area are consistently 
$\sim$10 m/day during the 2020–21 surge phase and drop to 
$\sim$0.1 m/day with the termination of the 2020–21 surge (Fig. 3b). Ice velocities at km 15 are near stagnant (
$ \lt $1 m/day) until October 2023 (the end of the 2023 melt season) when they increase to 
$\sim$1.5 m/day. Ice velocities then return to near stagnation (
$ \lt $1 m/day) during November 2023. Velocities steadily increase from 
$\sim$1 m/day to 
$\sim$4.5 m/day between late January 2024 and mid-June 2024 before decreasing rapidly during the latter part of the melt season (
$ \lt $1 m/day by early August 2024). The seasonal speedup during spring 2024 is of the type described in Liu and others (Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a). Ice velocities are near stagnation (
$ \lt $1 m/day) during the fall of 2024 and begin increasing in the early months of 2025. We observe velocities of 
$\sim$5 m/day in February 2025, 
$\sim$10 m/day in April, and a peak at 
$\sim$18 m/day in June 2025. Velocities at km 15 then decrease to 
$\sim$8 m/day by August 2025 (Fig. 4b).
Detailed view of hydrology and ice velocity at km 15 during a period of inter-annual ramping-up of ice flow speeds. (a) 
$dv/v$ values with y-axis flipped so that negative values plot up. (b) Ice velocities estimated with CautoRIFT shown in black, modeled surface runoff shown in red and glaciohydraulic tremor recorded at SE15 shown in purple.
Glaciohydraulic tremor recorded at SE15 is stable and consistently high during the 2020–21 surge, with an increase of 
$\sim$0.3 
$\times 10^{-8}$ m/s in April 2021 that precedes the onset of the melt-season by about one month. GHT values at SE15 decrease over the course of the 2021 melt season, then remain consistent over the next two years at 
$\sim$ 
$0.5 \times 10^{-8}$ m/s during the winters and at 
$\sim$ 
$1.9 \times 10^{-8}$ m/s during the melt seasons. The timing of GHT amplitude changes mirrors the timing of melt changes. Beginning in January 2024, GHT amplitudes decouple from melt variations, with SE15 GHT amplitude increasing non-linearly between January 2024 and August 2024 (Fig. 4b). Values reach 
$0.9 \times 10^{-8}$ m/s by April 2024, one month before the onset of the melt-season, and continue increasing to 
$\sim$ 
$ 3 \times 10^{-8}$ m/s by early August 2024. This summer tremor amplitude is nearly twice the melt season maxima of the two previous summers and is similar to values recorded during the 2021 mid-surge early melt-season. GHT values then rapidly decrease 1 month prior to melt decreases, and reach 
$\sim$ 
$1.7 \times 10^{-8}$ m/s by early September 2024, a level that is more consistent with values at the same time-of-year during previous melt seasons. By October 2024, SE15 GHT values decrease with the end of the melt season to reach average quiescent-phase wintertime levels. However, GHT values start increasing again in November 2024 to reach 
$\sim$ 
$1.3 \times 10^{-8}$ m/s in March 2025, the end of our SE15 record. The tremor level recorded in March 2025 is near-double the level recorded in March 2024 and triple the levels recorded in March 2022 and 2023.
During the 2020–21 surge phase, the 
$dv/v$ values are consistently negative for station pair SW14-SE15 (Fig. 3a). The values gradually decrease between October 2020 and late August 2021, when station SW14 stopped recording data. The SW14 station was repaired in August 2022, re-establishing the 
$dv/v$ record for SW14-SE15. The outage period saw a large increase in seismic wave velocities since the 
$dv/v$ values observed in August 2022 are around 0.75%, the highest in our record. From August 2022 to March 2025, the 
$dv/v$ values at km 15 exhibit an annually repeating pattern of relatively gradual decrease outside of the melt season and rapid increase during the melt season. As such, annual peaks in 
$dv/v$ occur in September and annual minima occur in May (Fig. 3a). The annual pattern overlays a multi-year decrease in 
$dv/v$ values between 2022 and 2025. The September maxima decrease from 
$\sim$0.75% in 2022 to 
$\sim$0.35% in 2024, and the May minima decrease from 
$\sim$-0.25% in 2023 to 
$\sim$-0.5% in 2024. Our record ends in March 2025 with 
$dv/v$ values around 
$-0.25$%.
Discussion
Interpretation of 
$dv/v$ records
Seismic interferometry on glaciers is an emerging field. The method has been employed to monitor a wide range of glacier dynamics (e.g. Mordret and others, Reference Mordret, Mikesell, Harig, Lipovsky and Prieto2016; Guillemot and others, Reference Guillemot2020; van Ginkel and others, Reference van Ginkel, Walter, Lindner, Halló, Huss and Fäh2025) but it can be challenging to isolate individual processes. Our confidence that the 
$dv/v$ signal reflects changes in hydraulic properties near the glacier bed stems in part from the isolation of other processes that could drive the observed variability in the 
$dv/v$ signal. The 0.75–1.5 Hz frequency window over which we evaluate 
$dv/v$ is sensitive to temporal changes in medium properties specifically at or near the glacier bed (S.3, Zhan (Reference Zhan2019)). This limits the direct impact on the 
$dv/v$ timeseries of changes in surface crevassing or water volumes at or near the glacier surface. However, changes at the surface such as seasonal snow loading can affect seismic velocities in the subsurface and thus affect 
$dv/v$ timeseries (e.g. Donaldson and others, Reference Donaldson, Winder, Caudron and White2019). The timeseries shown in Fig. 3 have undergone an empirically derived correction to account for the effect of snow loading, following methods established in Wang and others (Reference Wang, Brenguier, Campillo, Lecointre, Takeda and Aoki2017) (S.5). The correction is relatively minor: for cross-glacier station pairs the average absolute correction is 0.19% while the 5th-95th percentile range in the corrected 
$dv/v$ timeseries is 1.70%. The remaining seasonal variability in the corrected 
$dv/v$ signal (Fig. 3 and Supplementary Fig. S5) is out of phase with the single-side changes in 
$dv/v$, assumed to result from snow-loading. Thus, we do not anticipate that seasonal snow loading remains a significant factor in the 
$dv/v$ timeseries plotted in Fig. 3.
While there are significant variations in ice thickness throughout the surge cycle on Sít’ Kusá that could result in overburden pressure changes at and around the glacier bed, the sign of the 
$dv/v$ series is inconsistent with an influence from changes in ice thickness. Glacier thickness decreased by 
$\sim$80 m between June 2020 and July 2021 at km 15, while increasing by a similar amount at km 3 during the same time-period (Liu and others, Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a). None of the 
$dv/v$ records show large amplitude, long term changes during this period. The only long term trend in 
$dv/v$ values is a gradual decrease between 2022 and 2025 at km 15. This area is thickening during the quiescent phase, meaning the overburden pressure on the glacier bed is increasing, which should compress the till and increase the surface wave velocities rather than giving the observed velocity decrease (Mordret and others, Reference Mordret, Mikesell, Harig, Lipovsky and Prieto2016).
Finally, changes in 
$dv/v$ values clearly coincide with changes in the glaciohydraulic tremor signals. As such, we interpret the variability in 
$dv/v$ values shown in Fig. 3 as, at least to a large extent, driven by changes in the volume of water stored near the glacier bed. We do not distinguish between water storage in a till layer, in basal crevasses, or in subglacial cavities (S.3; Zhan (Reference Zhan2019)). There are observations of water being held in all these configurations (Hart and others, Reference Hart, Rose, Clayton and Martinez2015; Zhan Reference Zhan2019; Guillemot and others, Reference Guillemot2024) and all seem consistent with our observations and with the following discussion.
Drainage system configuration
In this section, we use our observations to make inferences about the configuration of the subglacial drainage system over the course of our record. We begin with a few remarks that will support our analysis. Distributed subglacial drainage systems can store more water than centralized, arborescent, channelized systems (e.g. Iken and others, Reference Iken, Röthlisberger, Flotron and Haeberli1983; Raymond and others, Reference Raymond, Benedict, Harrison, Echelmeyer and Sturm1995; Bartholomaus and others, Reference Bartholomaus, Anderson and Anderson2008). As such, hydrologically driven variation in our 
$dv/v$ measurements likely reflects change in distributed subglacial water rather than changes in efficient drainage channels. Water residence time in the distributed components of the subglacial hydraulic system is generally long, and in the channelized components of the system is generally short (e.g. Raymond and others, Reference Raymond, Benedict, Harrison, Echelmeyer and Sturm1995), even though a binary distinction has been contested in recent work (Hart and others, Reference Hart, Young, Baurley, Robson and Martinez2022; Hoffman and others, Reference Hoffman2016; Mejía and others, Reference Mejía2021; Terleth and others, Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024b). The surface of Sít’ Kusá was extremely crevassed and had very rough micro-topography throughout the study period. We observed no supraglacial runoff during our repeated visits to the glacier between 2020 and 2025. Thus, rather than into discrete moulins as is commonly found (e.g. Hoffman and others, Reference Hoffman2018; Scholzen and others, Reference Scholzen, Schuler and Gilbert2021), the delivery of surface runoff to the glacier bed on Sít’ Kusá must have occurred in a largely distributed manner during the study period (Terleth and others, Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024b). In the paragraphs that follow, we interpret our geophysical observations to infer the changing configuration of the subglacial drainage system at the three locations described above.
Receiving zone (km 3). In the receiving zone, 
$dv/v$ values that are mostly limited to within 
$\pm 0.2\%$, with no apparent inter-annual change and indistinct seasonality lead us to infer that the volume of water storage varies little over the course of our observations (Fig. 3e). We suggest that most surface runoff is routed through the lower part of the glacier quickly. Glaciohydraulic tremor is high during the surge but shows little seasonal variation afterward. The low GHT values in summer could reflect a system where changes in water discharge are accommodated by changes in hydraulic radius rather than changes in pressure gradient, reflecting an efficient, low pressure, drainage system (e.g. Gimbert and others, Reference Gimbert, Tsai, Amundson, Bartholomaus and Walter2016; Nanni and others, Reference Nanni2020). We suggest the receiving zone during the quiescent phase is predominantly underlain by an extensive and well connected channelized drainage network that efficiently drains most of the glacier area, leaving little water storage capacity. This network likely persists in winter, maintained by intermittent melting at low elevations and stored water release from further up-glacier (see next item). Alternatively, the observed behavior could reflect unvarying water storage in a mostly disconnected system (Rada and Schoof Reference Rada and Schoof2018), but we would expect such a configuration to be disrupted by the surge, and we see little change in storage even in the presence of high ice velocities (Fig. 3e).
Lower reservoir zone (km 8). Here, just above the dynamic balance line, changes in 
$dv/v$ lead us to infer a gradual decrease in stored water during the winter following surge termination (2021–22), after which water storage varies seasonally until the end of our record (Fig. 3c). The observed 
$dv/v$ pattern leads us to suggest that the bed beneath the lower reservoir zone accumulates water in the latter portion of the melt season, which is then released over the following winter. Notably, the GHT is lower and less coupled to the runoff timeseries from June to August of the 2022 and 2023 melt seasons than during September to October (Fig. 3c). The period of discrepancy between the signals coincides with the periods of rapid decrease in 
$dv/v$, hinting that a share of surface runoff enters distributed storage rather than the seismically loud drainage system. We suggest this area of the glacier has a ‘slow’ but well connected distributed system connected to a channelized, arborescent system. Such a subglacial drainage configuration would explain the rapid ‘filling’ during the melt season by distributed supply from surface runoff while slow rates of discharge through the system would explain the delayed draining (e.g. Hock and others, Reference Hock, Iken and Wangler1999). A model that seems consistent with our observations at km 8 is outlined in Hart and others (Reference Hart, Young, Baurley, Robson and Martinez2022), where subglacial water can exchange between a wide anastomosing drainage system and storage within subglacial ponds and the till layer. In this conceptual model, subglacial water is forced from the channel system into storage during the melt season. In the autumn, the channel system becomes less connected, and stored water becomes trapped in more isolated pockets. Water is gradually released from storage into the drainage system during brief speedup and uplift events in winter. Such speedup events are not evident in our velocity record (Fig. 3d) but we cannot conclude whether they take place on Sít’ Kusá because they could be too brief or localized to be captured by the satellite-based velocity tracking.
Mid-reservoir zone (km 15). In the middle of the reservoir zone, our data suggest a drainage system exhibiting hysteresis on a multi-year timescale. GHT amplitude and ice velocity evolution are markedly different in 2024 and 2025 relative to 2023, seemingly driven by a larger volume of basal water (the increasingly negative 
$dv/v$ of Fig. 3a,b) that accumulated between 2022 and 2025. We identify two time scales for the Sít’ Kusá drainage system: (1) an annual cycle of increased GHT in response to surface runoff along with variability in water volume and (2) a multi-year cycle that is marked by high amplitude drainage events, large changes in water volume and surges.
The annual cycle (1) of increased GHT during the melt season also manifests as a reduction in water volume during the melt seasons and an increase in water volume during the winters. Notably, this seasonality in water storage is opposite to the one observed at km 8. During the melt season, water storage at km 15 begins to fall (e.g. late June 2023), while water storage at km 8 rises 1–2 months later (e.g. August 2023). The contrasting behavior of the drainage system at km 15 relative to km 8 could be driven by low basal hydraulic conductivity, especially during winter. If any efficient channelized drainage collapses at the end of the melt season and water flow through remaining distributed drainage elements is strongly restricted, any further water supply to the bed, from delayed transfer of distributed englacial water (Lingle and Fatland, Reference Lingle and Fatland2003) or water transfers from further up-glacier (e.g. Terleth and others, Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024b), will remain trapped in a distributed and poorly connected drainage system during the winter. This would gradually build up basal water volume and cause the observed decrease in 
$dv/v$. The advent of high surface runoff rates at the start of the melt-season re-establishes some form of ‘fast’ drainage, as evidenced by the increases in GHT during the summers (Fig. 3a). This drainage network would free the water ‘trapped’ overwinter from the parts of the bed it connects to, lowering the overall basal water volume and driving the observed increase in 
$dv/v$ (Fig. 3a). As melt decreases later in the summer, the efficiency and extent of the drainage network decreases and stored water volume begins to increase again (decreases in 
$dv/v$ last from 
$\sim$September until 
$\sim$June of the subsequent year).
Ice velocities at km 15 emerge from stagnation after the 2023 summer, showing seasonal variability modulated by the subglacial drainage system (Fig. 4). A brief speedup in late October 2023 coincides with the end of the melt season (Fig. 4b), characteristic of the temporary rise in basal water pressure caused by the collapse of a channelized drainage network (e.g. Hock and others, Reference Hock, Iken and Wangler1999; Nanni and others, Reference Nanni, Scherler, Ayoub, Millan, Herman and Avouac2023; Narayanan and others, Reference Narayanan2025; Thota and others, Reference Thota, Vijay, Sommers, Banerjee, Mey and Motagh2026). The late-winter and spring of 2024 are marked by a gradual speed-up that ends abruptly in June 2024 with the onset of the melt season and the partial drainage of water stored at the bed (Fig. 4b). Such a pattern of ‘winter speedups’ has been observed on many Alaskan glaciers (e.g. Abe and Furuya Reference Abe and Furuya2015; Burgess and others, Reference Burgess, Forster and Larsen2013). We speculate that those glaciers might be subject to drainage system behavior similar to what we observe at km 15 on Sít’ Kusá, while Alaskan glaciers exhibiting summertime speedups (e.g. Bartholomaus and others, Reference Bartholomaus, Anderson and Anderson2011) are subject to a drainage regime closer to what we observe at km 8.
The multi-year cycle in water storage (2) modulates ice velocities at the time scale of the surge cycle. This cycle is marked by a gradual multi-year buildup of basal water during quiescent phases, interspersed with major drainage events at surge terminations. We suggest this periodicity stems from areas of the glacier bed with a poorly connected distributed drainage system that must remain largely inaccessible to the seasonally varying channelized drainage component described in (1), so that they can build water volume over multiple years. This is evident in our data with a net increase in water volume (decreasing 
$dv/v$) from summer to summer between 2022 and 2024 (Fig. 3a). This water could be sourced through mechanisms operating seasonally, such as transfers of englacial water to the bed or water transfer from up-glacier. Additionally, the gradual increase in water volume could be driven by water sources with a lower flux that would still have a significant impact over multiple years, such as exchange with groundwater reservoirs (Rada and Schoof Reference Rada and Schoof2018) and basal melting (e.g. Benn and others, Reference Benn, Jones, Luckman, Fürst, Hewitt and Sommer2019b; Sommers and others, Reference Sommers2023).
Insights into the surge mechanism
Sít’ Kusá historically follows a 5–8 year surge cycle, with 3.5–6.5 year quiescent phases that are marked by intensifying seasonal speedups (Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021; Liu and others, Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a). This pattern continues in the ice velocity records in Fig. 3. In the years following the termination of the 2020–21 surge, we see both an acceleration of the upper part of the reservoir zone (15 km) over the 2023–24 winter and spring that does not extend to 8 km, and then a more vigorous speed-up over the 2024–25 winter and spring. The spatial extent of the 2025 speedup extends well past the dynamic balance line at km 7, reaching km 3 in September 2025. Also in September 2025, the terminus of Sít’ Kusá began to advance. At the time of this writing, in June 2026, it has progressed by up to 1 km, although its rate of advance appears to have slowed. Meanwhile, our ongoing monitoring of satellite velocity products shows that ice velocities slowed dramatically at the end of the 2025 melt season and velocity at the terminus (1 km) has so far never exceeded 5 m/d (Supplementary Fig. S6). This more limited speedup distinguishes the 2025 event from the 2020–21 surge when ice velocities remained elevated through one winter and terminus speeds reached 30 m/d (Liu and others, Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a). We leave to future analysis whether the 2025 increase in ice velocities should be classified as an early and weak surge (as described previously for Variegated Glacier in Eisen and others (Reference Eisen, Harrison, Raymond, Echelmeyer, Bender and Gorda2005) and can occur at Sít’ Kusá (Nolan and others, Reference Nolan, Kochtitzky, Enderlin, McNabb and Kreutz2021)), or whether it is limited to a particularly strong seasonal speedup.
We interpret seasonal speedups as following similar mechanisms as the full-scale surges on Sít’ Kusá (Liu and others, Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024a; Terleth and others, Reference Terleth, Bartholomaus, Enderlin, Mikesell and Liu2024a). As such, we let the multi-modal behavior of the drainage system during our period of observation inform our understanding of surge dynamics on Sít’ Kusá. There is a net accumulation of basal water at km 15 between October 2022 and July 2024 (Fig. 3a), in what we suggested are the distributed and disconnected elements of the basal drainage system (Fig. 5a,b). Furthermore, the accumulation of basal water precedes winter-spring ice velocity increases at km 15, seemingly priming the glacier to speed up. These observations support the notion that gradual enthalpy accumulation at the glacier bed, along with increasing ice thickness, promotes surging (Sevestre and Benn Reference Sevestre and Benn2015; Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a, Reference Benn, Hewitt and Luckman2022). We have not investigated whether non-surge type glaciers see similar multi-year basal water accumulation; however, a more rigorous comparison with observations made on glaciers that do not surge would help clarify similarities and differences in basal water evolution (Raymond and others, Reference Raymond, Benedict, Harrison, Echelmeyer and Sturm1995). We hypothesize that other Alaskan glaciers that undergo spring speedups but do not surge would show the seasonal component of the 
$dv/v$ signal in Fig. 3a, but not the multi-year cycle.
Conceptual sketch of a surge mechanism driven by multi-year basal water accumulation and a positive feedback between distributed drainage connectivity and ice velocity. We suggest that the annual speedups on Sít’ Kusá result from a similar mechanism at shorter time scale, without the complete drainage of the base. This sketch is in part inspired by Figure 7 in (Iken and Truffer Reference Iken and Truffer1997).
Figure 5 Long description
There are four panels, panel a shows the early quiescent phase, with a top-down view of an arborescent drainage network with a few isolated water pockets. The on figure caption reads: Low basal water volume. Basal water builds gradually in disconnected pockets from distirbuted supply. Panel two shows the same arborescent network and more disconnected water pockets. The caption says: High basal water volume in a poorly connected distributed system. Glacier begins acceleration when basal water pockets are widepsread enough. Panel C shows the active phase. It is again the same arborescent network but the water pockets are widepread and connected. The caption reads: high basal water volume in a more connected distributed system. Strong velocity response to additional water supply. High ice velocities further increase connectivity. The final panel D shows surge termination, with a draiange network and few cavities connected to the drainage network. The caption reads: widespread connection between channelized and distributed components during the surge allows the basal water to drain.
During the 2023–24 and 2024–25 winters, the GHT signal recorded at SE15 increases several months prior to the onset of the melt season (Fig. 4b). The GHT signal indicates that stored water is mobilizing internally, coincident with an increase in ice velocities (Fig. 4b). The increase in subglacial water flow is only present at km 15 and does not lead to increased water flow at downglacier locations (Fig. 3). This suggests that water motion at these times is a relatively local phenomenon, such as an increase in connectivity and subsequent water pressure equilibration among regions of the bed (e.g. Doyle and others, Reference Doyle2022), rather than a widespread drainage event. We speculate that speedups on Sít’ Kusá initiate when basal water in distributed and disconnected regions of the bed reaches some threshold in total volume or in areal extent. The value of such a threshold would also depend on the driving stress (e.g. Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a); however, we do not explore that relationship here. When the threshold is reached, the sliding velocity and connectivity of the bed increase, as reflected by the localized increases in surface speed and in GHT. This would be consistent with existing literature showing that increases in overall basal water volume and ice velocity can lead to increased connectivity (e.g. Iken and others, Reference Iken, Röthlisberger, Flotron and Haeberli1983; Walder, Reference Walder1986; Engelhardt and Kamb, Reference Engelhardt and Kamb1997; Hock and others, Reference Hock, Iken and Wangler1999; Hart and others, Reference Hart, Young, Baurley, Robson and Martinez2022; Doyle and others, Reference Doyle2022; Rada Giacaman and Schoof, Reference Rada Giacaman and Schoof2023). Moreover, a connected distributed drainage system can promote a more sensitive dynamic response to subglacial water pressure variations than a disconnected distributed system (Iken and Truffer, Reference Iken and Truffer1997). It seems plausible that the increased connectivity of the Sít’ Kusá drainage system during the 2023–24 and 2024–25 winters further promotes the increases in ice velocities. As such, the reorganization of the drainage system could provide a threshold, or at least a positive feedback mechanism, explaining the rapid increase in ice velocities under seemingly small changes in forcing from ice thickness or basal water volume (Fig. 5c). As the 2024 melt season starts, GHT amplitudes recorded at SE15 are nearly double the amplitudes recorded during previous summers, while ice velocities drop to near stagnation (
$ \lt $1 m/day). These coincident changes could occur as the result of the previously stored subglacial water now being available to be drained through the ‘fast’ channelized drainage component. We suggest this ‘freeing’ of water stored in previously disconnected regions of the bed also occurs, to a larger degree, during full-scale surges (Fig. 5d). This mechanism would explain the strong decrease in stored water during the 2021 surge termination. It is also consistent with the sustained high GHT levels during the 2020–21 surge and could provide a source for the high rates of frontal discharge often observed at surge terminations (e.g. Kamb and others, Reference Kamb1985; Harrison and Post, Reference Harrison and Post2003; Benn and others, Reference Benn, Jones, Luckman, Fürst, Hewitt and Sommer2019b). While speculative, this mechanism of disconnected to connected distributed basal drainage (Fig. 5) could be interesting to explore further in the context of surge initiation. Early models of hydrologically driven surging do not contain an in-depth explanation for surge initiation (Kamb and others, Reference Kamb1985; Lingle and Fatland, Reference Lingle and Fatland2003), and the enthalpy model relies on the crevassing feedback for a mechanism driving rapid (
$ \lt $1 year) ice velocity increase (Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a; Terleth and others, Reference Terleth, Bartholomaus, Enderlin, Mikesell and Liu2024a). Although observational evidence supports a crevassing feedback on Svalbard (e.g. Dunse and others, Reference Dunse, Schellenberger, Hagen, Kääb, Schuler and Reijmer2015; Benn and others, Reference Benn, Jones, Luckman, Fürst, Hewitt and Sommer2019b; Nanni and others, Reference Nanni2025), the mechanism seems unlikely to operate on persistently and extensively crevassed glaciers such as Sít’ Kusá. More broadly, our hydrological interpretation assigns an important role to the spatial configuration of the subglacial drainage system and to the evolution over time of connectivity within that system that is not yet fully represented in the community’s understanding of the surge mechanism (Lovell and others, Reference Lovell2026). We suggest that the spatio-temporal evolution of subglacial drainage during all stages of the surge cycle merits further exploration through both modeling and observational work.
Conclusion
In this work, we have applied seismic noise interferometry to monitor the evolution of water stored at or near the bed of Sít’ Kusá during major changes in ice flow behavior. From our observations, we infer that the receiving zone has relatively little capacity for water storage during the quiescent phase and is well drained throughout the observed part of the surge cycle. Basal water storage in the middle of the reservoir zone (15 km) increases during winter, when poor connectivity promotes the retention of water in a distributed system. Basal water volume in the same area decreases during summer when increased connectivity and efficiency, driven by high rates of surface runoff, promote water drainage from the distributed system. This pattern is consistent with the winter speedups and summer slowdowns widely observed on Alaskan glaciers (Burgess and others, Reference Burgess, Forster and Larsen2013). We also identify multi-year basal water accumulation in the reservoir zone at km 15 over the course of the quiescent phase, eventually supporting the onset of winter speedups and a strong increase in subglacial water flow in 2024 and 2025. Summer increases in basal connectivity might explain partial drainage of the glacier bed around km 15, which seem to terminate the 2024 and 2025 speed ups. A similar mechanism at a larger scale likely terminates the full scale surge in 2021. We suggest changes in the connectivity of the subglacial drainage system could prove important drivers of the dynamic behavior of the glacier during both seasonal speedups and surges. Our results are broadly consistent with predictions made by enthalpy-based models of glacier surging. Future work that provides comparison of multi-year water storage change between surging and non-surging glaciers, and that considers spatial variability in drainage conditions, is likely to prove valuable in untangling the relation between subglacial water and ice dynamics.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10206.
Data availability statement
Seismic data are archived and freely available via two networks at the Earthscope/IRIS DMC Mikesell and others (Reference Mikesell, Bartholomaus, Liu, Terleth and Enderlin2020); Bartholomaus and others (Reference Bartholomaus, Mikesell and Enderlin2024), or as the virtual network _SIT_KUSA (https://ds.iris.edu/mda/_SIT_KUSA/). Code for analysis with the seisgo package and for making the principal figures is repositoried at https://github.com/yoramterleth/NOISE. Code for seismic data processing toward glaciohydraulic tremor timeseries is available at https://doi.org/10.5281/zenodo.8102681. The timeseries presented in the paper are available at https://doi.org/10.5281/zenodo.20761017.
Acknowledgements
This work was funded through the NSF Arctic Natural Sciences program, via awards 1954021 and 1954006. T.D.M. acknowledges basic funding support to NGI from the Research Council of Norway. We thank Kate Bollen, Thomas Otheim, Chris Miele, Dakota Pyles, Jason Amundson, Jake Anderson, Galen Dossin, Susan Detweiler, Parker Wilkerson, Keeya Beausoleil, Hans Munich and Tanya Hutchins at Coastal Airline Services and Kelix DeZagottis at Coastal Helicopters for their roles in the data collection effort. Some seismic and GNSS instruments were provided by the EarthScope Consortium. The facilities of EarthScope Consortium are supported by the National Science Foundation’s Seismological Facility for the Advancement of Geoscience (SAGE) Award under Cooperative Support Agreement EAR-1851048 and Geodetic Facility for the Advancement of Geoscience (GAGE) Award under NSF Cooperative Agreement EAR-1724794. We would like to thank Martin Truffer and Elizabeth Cassel for helpful discussions leading to this manuscript, as well as two anonymous reviewers and Scientific Editor B. Kulessa and Associate Chief Editor H. Jiskoot for providing constructive feedback.
Open access
