1. Introduction
Sít’ Tlein (“Big Glacier” in Tlingit, temporarily known as “Malaspina Glacier” (pers. comm. Yakutat Elders, 2023)), located in Southeast Alaska, is the largest glacier-complex in the world outside of the ice sheets (Windnagel and others, Reference Windnagel, Hock, Maussion, Paul, Rastner, Raup and Zemp2023). Its piedmont lobe holds approximately 700 
$\mathrm{km}^{3}$ of ice almost entirely in the ablation zone, and is separated from the Pacific Ocean by moraines of variable width (Figure 1, (Thompson and others, Reference Thompson2021)). The glacier originates in the Saint Elias mountain range and is exposed to strong climatic gradients with topographic rain-shadow effects (Whiteman, Reference Whiteman2000, Wendler and others, Reference Wendler, Gordon and Stuefer2017).
Radar depth sounding shows that the piedmont lobe’s bed is mainly below sea level (Tober and others, Reference Tober, Holt, Christoffersen, Truffer, Larsen, Brinkerhoff and Mooneyham2023) with deep subglacial troughs, which are sometimes related to geological faults (Cotton and others, Reference Cotton, Bruhn, Sauber, Burgess and Forster2014), some of which extend close to peripheral lakes on the outskirts of the glacier’s lobe. In 2021, our field campaign discovered seawater intrusion in one of the peripheral lakes (Thompson and others, Reference Thompson2021), and satellite images have shown an increase in peripheral lake area during the last 10 years. The small ocean separation together with the over-deepened bed make the piedmont lobe susceptible to further intrusion by warm seawater, with a potential for rapid glacier retreat (Motyka and Beget, Reference Motyka and Beget1996, Post and others, Reference Post, O’Neel, Motyka and Streveler2011, Brinkerhoff and others, Reference Brinkerhoff, Truffer and Aschwanden2017). This could make Sít’ Tlein Alaska’s greatest sea level rise contributor, with Alaska already projected to lead glacial sea level rise contribution outside of Greenland and Antarctica during this century (Edwards and others, Reference Edwards2021, Rounce and others, Reference Rounce, Hock, Maussion, Hugonnet, Kochtitzky, Huss, Berthier, Brinkerhoff, Compagno and Copland2023).
Sít’ Tlein’s flow is characterized by surging. The unusually high flow rates of surging glaciers, as shown by the heavy crevassing and accompanied by glacier advance, were recognized more than 100 years ago by Tarr and Martin Reference Tarr and Martin1914 along the coast of Southeast Alaska. They attributed surges to earthquakes, but this was not confirmed by subsequent events. Meier and Post Reference Meier and Post1969 first proposed a definition of surges as quasi-periodic events of one to several years duration with glacier velocities up to two orders of magnitude higher than quiescent flow. These events are generally separated by years to decades. Since then, global inventories of surge-type glaciers have been compiled (e.g. Guillet and others, Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022, Kääb and others, Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi2023), and Sevestre and Benn Reference Sevestre and Benn(2015) grouped those glaciers into preferred climatic zones. A detailed process study at Variegated Glacier (Kamb and others, Reference Kamb1985) highlighted the important role of pressurizedsubglacial water in enabling the high flow velocities, while other studies emphasized the role of frozen/temperate transitions at the glacier bed (e.g. Clarke and others, Reference Clarke, Collins and Thompson1984). These mechanisms were put into a consistent framework via the enthalpy model by Benn and others Reference Benn, Hewitt and Luckman(2022). Despite some progress (e.g. Thøgersen and others, Reference Thøgersen, Gilbert, Schuler and Malthe-Sorenssen2019, Minchew and Meyer, Reference Minchew and Meyer2020), physical mechanisms remain insufficiently understood to make predictive surge models possible. In particular, it is difficult to even agree on a clear definition for surging, and some authors have argued for a full spectrum of flow behaviors ranging from steady flow to total glacier collapse (Herreid and Truffer, Reference Herreid and Truffer2016, Kääb and others, Reference Kääb2021).
Understanding the glacier’s dynamics is crucial to predicting glacier evolution during the next few decades. This is difficult due to complicated surge-type flow patterns in Sít’ Tlein’s three main tributaries (Agassiz, Seward, and Marvine-Hayden), which converge into a single lobe. These surges have variable amounts of impact, from affecting most of the glacier surface to only a small portion of it. They are generally asynchronous, but there is at least one documented occurrence of simultaneous surging (Muskett and others, Reference Muskett, Lingle, Sauber, Post, Tangborn and Rabus2008). To date, no predictive model can forecast the timing and extent of these complex dynamical patterns.
Here, we will use a combination of published reports from recorded Indigenous oral history and early European explorers, scientific reports from the pre-satellite era, and a compilation of satellite derived data to compile a century-scale history of ice flow of the Sít’ Tlein and its tributaries. We adopt a technique new to glaciology, Stochastic Matrix Factorization, to fill gaps in the satellite-derived velocity record, and we discuss the implications of our results on the health of the glacier.
2. Methods
We use three types of data sources to study the flow and surge histories: pre-satellite (aerial pictures, expedition, and scientific reports), Landsat images, and ITS_LIVE (Gardner and others, Reference Gardner, Fahnestock and Scambos2022). Every source captures distinct aspects of ice flow and surges (extent, duration) at different spatiotemporal resolutions. Thus, we provide dataset-specific surge definitions. This is in line with previous work that required adjustments for surge criteria depending on the methodology used for classification (Jiskoot and others, Reference Jiskoot, Murray and Boyle2000, Copland and others, Reference Copland, Sharp and Dowdeswell2003, Sevestre and Benn, Reference Sevestre and Benn2015). Several authors have discussed the limitations of those surge definitions (Sevestre and Benn, Reference Sevestre and Benn2015, Herreid and Truffer, Reference Herreid and Truffer2016, Truffer and others, Reference Truffer, Haeberli and Whiteman2021, Benn and others, Reference Benn, Hewitt and Luckman2022, Guillet and others, Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022). Sít’ Tlein’s dynamics cover a broad, continuous range, making it challenging to identify distinct events without defined thresholds. To address this, we adopt threshold-based classifications to clearly define surges and facilitate the discussion. Our Landsat surge classification considers a spatial displacement threshold, while the ITS_LIVE surge classification involves both spatial and temporal dimensions, enabled by the dataset’s enhanced temporal resolution. We also calculate yearly fluxes at the throat (transition area from glacier trunk to glacier lobe), and near the terminus of each tributary. We then detail the Stochastic Matrix Factorization interpolation method and apply it to the ITS_LIVE dataset. Finally, we describe the novel data decomposition method ROCK-PCA and apply it to the interpolated velocity product. These methods allow an analysis of ice-flow dynamics for large datasets.
2.1. Pre-satellite reports
Sít’ Tlein is located on Tlingit land that has been inhabited for thousands of years (Schurr and others, Reference Schurr2012, Shorty, Reference Shorty2016). Some sources mention Sít’ Tlein’s impact on communities and the landscape in the early 19th century (Cruikshank, Reference Cruikshank2001). From the 1800s to the 1950s, mountaineering and scientific expeditions (Washburn and Sharp expeditions Tarr, Reference Tarr1907a, Tarr and Martin, Reference Tarr and Martin1914, Cruikshank, Reference Cruikshank2001), brought written and photographic evidence of the glacier’s state. We use those sources to determine whether the glacier was in quiescence or actively surging at the time of the historical encounters.
During quiescence, Sít’ Tlein’s lobe is flat and not heavily crevassed. A surge will generate chaotic surface features and lead to moving ice reaching the terminus. If a report mentions travelers crossing the ice on foot, we will consider it as local evidence of quiescence. However, if there is a broken ice surface, demolished trees, or changes at the terminus, we consider it as surge evidence.
2.2. Displacement from Landsat records
We quantify ice displacements using Landsat 1 to 9 cloud-free summer images merged into mosaics and cropped to Sít’ Tlein’s lobe. The time between each mosaic is centered around a year (figure S1). We compare the annual mosaics spanning from 1973 to 2023 to find large displacements indicative of surges.
We assess the three tributaries separately (Figure 1). We further distinguish displacements of the trunks (constrained laterally) from displacements in the lobes (less lateral confinement) between each mosaic. We use moraines as tracking features to manually assess flow displacement and identify potential surges (Herreid and Truffer, Reference Herreid and Truffer2016). We create an index (Figure 2) describing the type of displacement a tributary is experiencing. This index is composed of 2 numbers and 1 letter. The first number indicates displacement in the tributary’s trunk, the second in its lobe, with values ranging between 1 (no displacement), 2 (
$ \lt $500 m), and 3 (
$ \gt $500 m). A surge is expected to exceed a 500 m displacement threshold, while non-surging years might not quite reach it. This value was chosen after observing the entire record of yearly mosaics for the three tributaries. We choose the flux gates (Figure 1) of each tributary as a separation between lobes and trunks. To quantify the extent of each surge, we introduce the third index, a letter, either F for ‘Full’ (displacements reach further than half of the lobe) or P for ‘Partial’ (displacements reach less than half of the lobe). In most cases, an ‘F’ case will propagate to the terminus. An example of each case is shown in figure S8.
Figure 1. (a) Contextual map of Sít’ Tlein. Velocity transects are shown by the colored lines (blue, orange, purple, and pink), and the important landforms are outlined in yellow. The background color depicts the glacier’s bed elevation (Tober and others, Reference Tober, Holt, Christoffersen, Truffer, Larsen, Brinkerhoff and Mooneyham2023), elevations close to sea-level appear as a line in the legend due to a compression of the colorbar, although they are represented by a gradient of green on the figure. (b) The image shows the iconic folded moraines of Sít’ Tlein (Sentinel-2, August 4th 2024, retrieved from the Copernicus Browser, courtesy of ESA).
Figure 2. Surge classification based on Landsat displacements.
2.3. Velocities
2.3.1. ITS_LIVE velocity dataset
Ice velocities are quantified by analyzing 39 years of surface velocities (1984 to 2023) generated with the autonomous Repeat Image Feature Tracking (autoRIFT) algorithm (Gardner and others, Reference Gardner, Moholdt, Scambos, Fahnstock, Ligtenberg, Van Den Broeke and Nilsson2018, Lei and others, Reference Lei, Gardner and Agram2021) for the ITS_LIVE project, interpolated to a 5-day timestep. This global dataset uses optical and Synthetic Aperture Radar (SAR) images from various satellites to calculate pixel displacements between two images through a cross-correlation technique. The number of available images increased dramatically twice during the record: in 2014, when Landsat-8 products started becoming available, and in 2017 with Sentinel 2 products. The velocity dataset was downloaded as a netcdf file from the ITS_LIVE Amazon Web Service server (code adapted from https://github.com/nasa-jpl/ITS_LIVE); only slices with more than 40% of valid data were accepted, regardless of the image pair time separation. This dataset (datacube) is a stack in time of 2D matrices with its three axes representing time, northing, and easting (t, y, x). We group each entry within a 5-day window and calculate the median. We note that this approach does not take into account the finite time interval over which the velocity is measured (i.e., the interval of time between repeat images from which the velocities are derived (see Minchew and others (Reference Minchew, Simons, Riel and Milillo2017), Greene and others (Reference Greene, Gardner and Andrews2020), Charrier and others (Reference Charrier, Yan, Koeniguer, Leinss and Trouvé2021) for alternative approaches). Then, we interpolate missing values using a method based on Stochastic Matrix Factorization (Mnih and Salakhutdinov, Reference Mnih and Salakhutdinov2007), detailed in the next subsection.
We analyze three items: ice velocities along longitudinal transects on each tributary, ice-flow direction on the Seward tributary transect, and spatial patterns of surge velocities. The first item allows us to quickly identify main events in the dynamics of Sít’ Tlein by capturing the flow variations at the transition between the trunks and lobes. The second item is generated by creating an index that calculates the mean speed-weighted position along the Seward transect as a proxy of flow direction:
\begin{equation}
\langle x \rangle = \frac{\sum x_i v_i}{\sum v_i},
\end{equation}where vi is the speed at position xi along the transect. The purpose of analyzing this quantity is to assess subtle changes in directionality as the ice enters the main lobe.
Finally, the third item represents the spatial extent of the strongest surges. A location is considered surging if its monthly median speed exceeds four times the time serie’s median, for at least five months within a 12 months rolling window. This condition filters out seasonal speed-ups. Once the surge dates have been obtained, we determine which year has the highest surge velocities. The resulting map shows when each location was the most strongly impacted by a surge.
2.3.2. Ice flux
We trace two flux gates for each tributary (Figure 1): one where the ice exits the trunk and spills into the lobe (throat of the tributary), and one close to the terminus of the lobe. The near-terminus fluxes of Agassiz and Marvine glaciers are very small, and we will not discuss them further. The gates at the throats are curved approximately perpendicularly to the ice flow. This shape marks the onset of flow divergence (figure S2) in the lobe. The throat flux gates are close to the equilibrium line altitude (Brinkerhoff and others, Reference Brinkerhoff2025) and allow us to evaluate the amount of transferred ice from the accumulation area to the ablation area.
We use surface velocities and ice thickness to estimate ice fluxes through the flux gates. Ice thickness is obtained by differencing the complete Interferometric Synthetic Aperture Radar DEM (IFSAR, USGS 2013), and the bed elevation from Tober and others Reference Tober, Holt, Christoffersen, Truffer, Larsen, Brinkerhoff and Mooneyham2023. We make one main assumption to calculate the fluxes: the surface velocity is related to the depth-averaged velocity by a factor k: 
$\bar{v} = k v_\mathrm{surface}$. For motion of an infinite plate with a flow law exponent of n = 3, this factor is bounded between k = 0.8 (no slip) and k = 1 (sliding without internal deformation). Since we expect significant amounts of basal sliding (the ice is assumed temperate), we choose a value of 
$k=0.95 \pm 0.05$. Therefore, the flux of ice through the flux gate is approximated as:
\begin{align}
\int_F \bar{\mathbf{v}} H \cdot \mathbf{n} \mathrm{d} S \approx \sum_{i}^{N} \bar{\mathbf{v}}_i H_i \cdot \mathbf{n}_i \Delta_i = \sum_{i}^{N} H_i (v_x \Delta y - v_y \Delta x),
\end{align} with 
$\bar{\mathbf{v}} = (v_x, v_y)$ the depth-average velocity, H the ice thickness, n a vector normal to the flux gate, N the number of segments in the flux gate, and F the path along the flux gate. 
$\Delta_i$ is the length of the ith segment of the flux-gate. We calculate fluxes from the yearly velocity mosaics provided by ITS_LIVE to reduce temporal variability.
We obtain the uncertainty of the flux from standard error propagation. ITS_LIVE and the ice thickness datasets both provide error estimates.
Since the ice flux is strongly dominated by the the Seward trunk, we calculate normalized fluxes to compare temporal variations among tributaries using: 
$\frac{q-\mu}{\sigma}$. Here, q represents the flux, µ is the time-averaged flux, and σ is its standard deviation.
2.4. Stochastic matrix factorization interpolation
The ITS_LIVE dataset is a stack of 2D velocity fields along time, thus forming a 3D array (X). Each element of this array falls into one of three categories: on-ice, off-ice, or invalid (no data). We determine on-ice locations based on outlines from the RGI Consortium (2017) version 6.
The interpolation assumes that the flattened matrix Xflat, containing values at different spatial locations (rows) and times (columns), can be decomposed into limited sets of spatial and temporal orthogonal basis functions, U and V. These bases fully characterize the spatio-temporal variability of the dataset, with their product 
$X_{pred} = UV^T$ being a prediction of the spatiotemporal velocities without gaps. However, computing the singular value decomposition (SVD) of datasets with missing data is not possible, so we adopt an alternative method based on Stochastic Matrix Factorization (Mnih and Salakhutdinov, Reference Mnih and Salakhutdinov2007). This method uses spatiotemporal variance patterns to infer missing values in our largely sparse but organized dataset.
To apply the interpolation to the ITS_LIVE dataset, we make the following assumptions:
1. Our dataset X (or in its flattened form,
$\boldsymbol{\hat{X}_{flat}}$) has a low rank (hence its factorization into two low rank matrices, U and V)2. Each entry shares a similar normally distributed error σ 2
The Stochastic Matrix Factorization Interpolation revolves around minimizing the L2 norm:
\begin{equation}
\Vert\boldsymbol{\hat{X}_{flat}}-\boldsymbol{\widehat{UV}}\Vert^{2}.
\end{equation}This norm is the squared differences between our dataset X and the product of two matrices U and V, evaluated only for on-ice values (hence the hat notation on both matrices). While this minimization is evaluated only on valid entries, it applies to the entire dataset, such that missing entries are interpolated. The objective is to find the optimal U and V whose product is as close as possible to X.
The algorithm proceeds as follows:
1. Initialization
We flatten X in the spatial dimensions to obtain a 2D matrix Xflat of dimension (RC, T), with R and C the number of rows and columns, and T the number of timesteps. We mask the invalid entries and obtain
$\boldsymbol{\hat{X}}_{flat}$. At each time step, we compute the median off-ice velocity and discard entries where the absolute value exceeds the 95th percentile of the absolute time series. We then subtract the remaining median off-ice velocities from all on-ice pixels. We then standardize the data following:
$\boldsymbol{\bar{X}}_{flat} = \frac{\boldsymbol{\hat{X}_{flat}} - \mu}{std}$ where µ and std are the the mean and standard deviation of
$\boldsymbol{\hat{X}_{flat}}$. Thus
$\boldsymbol{\bar{X}}_{flat}$ is a standardized, masked, and flattened array. We solve a regularized minimization problem by finding two matrices U and V of dimensions (RC, M) and (M, T) that approximate X subject to regularization terms measuring spatial and temporal smoothness. The product UV is a flattened interpolation of X, with U and V approximating the SVD components of X if it had no missing values. The matrix UV contains the M spatial modes of X, where each column captures a portion of the dataset’s spatial variance. V represents the temporal strength of the modes; its dimensions are
$\left[C, T \right]$. We selected M = 10 modes, as higher-order modes beyond 10 typically capture noise rather than meaningful variance. This choice is supported by the SVD of an interpolated dataset using M = 100 modes (a conservative value ensuring over 99% of the variance is included), showing that the first 10 modes account for approximately 95% of the total variance. Finally, we define
$\boldsymbol{\widehat{UV}}$ as the matrix UV with masked invalid entries.2. Least Squares with minimization of L2 norm
To minimize expression 3, we use the steepest gradient descent fitting the matrices U and V. We use four regularization terms with associated regularization weights (λu, λv, λx, λt), determining the contribution of each term relative to the goodness of data fit:
(4)\begin{align}
E_{misfit} &= \overline{(\boldsymbol{\bar{X}_{flat}}-\boldsymbol{\widehat{UV}})^{2}}
\end{align}(5)\begin{align}
E_{reg} &= \frac{\lambda_{u}\Vert\textbf{U}\Vert_{Frob}^{2}}{(RC)^{2}} + \frac{\lambda_{v}\Vert\textbf{V}\Vert_{Frob}^{2}}{(RC)^{2}}
\end{align}(6)\begin{align}
E_{space} &= \frac{\lambda_x}{R} \left(\sum\partial U_{R}^{2} + \sum\partial U_{C}^{2}\right)
\end{align}(7)\begin{align}
E_{time} &= \frac{\lambda_t}{R} \sum\partial t^{2}
\end{align}The misfit is only computed for the non-masked entries of X. Norm regularization uses the squared Frobenius norm weighted by λu and λv (row and column spaces) to enforce U and V to have comparable magnitudes. Spatial smoothness is ensured via the squared norms discrete gradients
$\sum\partial U_{R}^{2}$ and
$\sum\partial U_{C}^{2}$. Finally, the temporal regularization promotes smoothness over time using the squared differences
$\sum\partial t^{2}$. Once these terms have been calculated, we compute a total regularized cost functional:
(8)\begin{equation}
E_{tot} = E_{misfit} + E_{reg} + E_{space} + E_{time}.
\end{equation}Minimizing this functional provides us with optimized matrices U and V.
3. Reconstruction of the 3D interpolated array
The optimized matrix product
$\textbf{U}\textbf{V}$ from the optimization represents an interpolated dataset without gaps. We call XC the unflattened and interpolated
$\textbf{U}\textbf{V}$ product matching the 3D shape of X. Most of the noise has been removed, assuming that the modes beyond the first 10 primarily capture the dataset’s noise.4. Quality assessment
Finally, we evaluate the quality of the interpolation by calculating the difference between each corresponding entry in both interpolated and initial datasets:
$\boldsymbol{\hat{X}}_{C} - \boldsymbol{\hat{X}}$. We fit a normal and a Student’s t distribution to the histogram of the differences. We also calculate the average and standard deviation of the normalized differences per velocity bin from 0 to 2000 m a−1. We finish by comparing the differences in flow direction between the two datasets.
2.5. Mode analysis of the variance
We extract the spatio-temporal patterns of the velocity dataset by using a technique derived from the Principal Component Analysis (PCA), called Rotated Complex Kernel—Principal Component Analysis (ROCK-PCA) (Bueso and others, Reference Bueso, Piles and Camps-Valls2020), which has the advantage of decomposing nonlinear modes (Cristianini and Scholkopf, Reference Cristianini and Scholkopf2002) and decoupling them in space and time (Horel, Reference Horel19841984, Esquivel and Messina, Reference Esquivel and Messina2008). We find that our dataset is quasi-stationary and non-linear, based on the Augmented Dickey-Fuller test (Fuller, Reference Fuller2009) and low Pearson coefficients from linear fits to the time-averaged velocities. The ROCK-PCA differs from a standard PCA in that it calculates complex modes, allowing for the derivation of phased differences for each mode, and the modes are not required to be orthogonal.
3. Results
3.1. Reports
We find evidence for advances, surges, retreats, and quiescence as early as 1837 to 1956, with several hiatuses of observations. We report the findings in their chronological order below.
• Late Holocene Sít’ Tlein retreat: The earliest evidence of retreat is shown by two sources. The Tlingit tribe in Yakutat uses ‘Sít’ Tlein’ (’Big Glacier’ in Tlingit) for both Malaspina and Hubbard glaciers. This indicates that they were once one single body of ice reaching to the coast (Cruikshank, Reference Cruikshank2001, Thompson and others, Reference Thompson, Loso, Mooneyham, Tober, Larsen and Holt2024). Old moraines located in Disenchantment Bay, disappearing under Sít’ Tlein’s current eastern side of the piedmont lobe (Barclay and others, Reference Barclay, Calkin and Wiles(2001), Calkin and others, Reference Calkin, Wiles and Barclay2001) show that the moraine was deposited by Hubbard Glacier. This means two things: the two glaciers were probably joined (which corroborates the Tlingit name), and Sít’ Tlein had retreated further than its contemporary position.
• Late Holocene/1837 Sít’ Tlein advance: The Vancouver expedition (1794 in Icy Bay) reported trees along the shore near Sít’ Tlein. In 1837, the Belcher expedition’s observations found no evidence of trees along the shore (Cruikshank, Reference Cruikshank2001), indicating a possible advance of Sít’ Tlein’s lobe.
• 1890/91 Sít’ Tlein quiescence: In 1890 and 1891, Russell crossed the ice from the Marvine Glacier to the West twice (Russell, Reference Russell1893, Tarr, Reference Tarr1907a). Other expeditions from the Duke of Abruzzi and Bryant both used the glacier because of its easy terrain (Tarr, Reference Tarr1907a), indicating quiescence along their routes.
• 1907 Marvine surge: Tarr (Reference Tarr1907a, Reference Tarr1907b) observed a chaotic glacier surface forbidding any crossing by foot. They recorded photographic evidence of the fractured and moving terminus of Marvine, and also of trees being destroyed by the ice. They observed from a distance that, according to a member of this expedition and also of 1897’s Abruzzi’s expedition, Seward’s throat looked “far more broken”.
• 1935 Marvine surge: After a hiatus of 30 years, Bradford Washburn (Washburn, Reference Washburn1935) provides aerial photographs of eastern Sít’ Tlein showing Marvine’s lobe undergoing a strong surge based on how broken up the surface is, reminiscent of the photographs from Tarr in 1906.
• 1948/53 Seward quiescence: Sharp (Reference Sharp1951, Reference Sharp1953, Reference Sharp1958) provide a U.S. Coast Guard aerial pictures mosaic from 1953, and an extensive description of fieldwork spanning 1948 to 1953, focused solely on the Seward tributary. They mention the lobe being accessible by foot, while the photographs do not show extensive crevasses.
• 1954 Seward surge: In the same publication, Sharp Reference Sharp(1958) dedicates a chapter on the sudden increase in activity of Seward’s lobe from 1954 to 1956. In contrast with their previous years’ fieldwork on the ice, they noted increased noise from crevasses and debris falling into said crevasses, along with an increase in crevasses size and quantity. They indicate that other parties observed such activity continuing in 1955 and 1956.
After Sharp’s work, we are not aware of any relevant records until the first Landsat images in 1972.
3.2. Landsat
Table 1 summarizes ice activity from 1973 to 2021 for each tributary based on the XXF/P classification explained in the methods.
Table 1. Classification of ice movement from Landsat mosaic spanning 1972-2021. Full surges have been highlighted in orange, partial surges in blue. 1 = no displacement, 2 = displacement 
$ \lt $ 500 m, 3 = displacement 
$ \gt $ 500 m. “F” is a displacement propagating past the first half of a lobe, “P” otherwise
Agassiz fully surged twice over the 50-year record. Most of the surges we identified do not propagate further than the middle of its lobe, but we record two occurrences (1986, 2002) reaching its terminus (identified as “33F” in the table). The 1986/87 surge was synchronous with the other tributaries, as the entire Sít’ Tlein lobe experienced unrest at this time. Agassiz’s flow oscillates slightly from southwest to northeast and produces banded moraines on its eastern part (figure S3). The ice flow forms a medial moraine marking the suture line between the Seward and Marvine lobes.
Marvine’s two primary feeder trunks result in more complex dynamics, whereas Seward and Agassiz each have only one trunk. We document four partial surges in Marvine’s trunks that stopped in the upper half of the lobe in 1974/75, 1984, 1986, and 2013. We classify the 1984 and 1986 surges as separate events because no significant displacements were observed between 1985 and 1986. 1987 is the only surge that moved ice to Marvine’s terminus, displacing it by almost 10 kilometers (figure S4). It activated both trunks and the lobe in its entirety, and occurred synchronously with Seward’s surge. Marvine’s lobe displays banded moraines indicative of a flow oscillation at the throat of the trunk. They are more pronounced than Agassiz’s but not as spectacular as Seward’s. Like Agassiz, the border with Seward’s lobe is marked by lower displacements and stacked moraines.
The Seward tributary surged five times in the 50-year record; all of its surges propagated past the upper half of the lobe. They took place in 1977, 1981/82, 1986/87, 2008/09, and 2020/21. With Seward’s trunk flowing consistently above 500 m of displacement per year, our surge threshold for Landsat displacements is frequently exceeded. Hence, our surge assessment for Seward is only valid for its lobe. The 1977 and 2021 surges bear a strong spatial resemblance in how local they are and propagate along a well-defined subglacial valley to Fountain Stream (Figure 1), displacing the terminus by roughly 500 m. Seward’s surges are peculiar in that they only activate the ice in specific parts of the lobe, leaving other areas completely unaffected. In 1986/87, the eastern part of the lobe spectacularly surged while the western part was not impacted. It coincided with Agassiz’s surge, and most importantly with Marvine’s surge. It moved the terminus by a few kilometers and almost entirely covered Malaspina Lake. The ice in the zone of influence of the surge was mobilized regardless of its surface cover (debris or clean ice). Surges of this magnitude can redefine entire portions of the lobe and change the boundaries of the glacier (figure S5). Outside of these periods, ice displacements are generally constrained to areas of clean ice and tend to slow by an order of magnitude over a short distance at the outskirts of the lobe, where the moraines stack. On a few occasions, the lobe shows movement below 500 m per year that stops a few hundred meters from the terminus.
3.3. ITS_LIVE
3.3.1. Velocity interpolation and modes
Figure 3. (a) Amplitude of the first 6 modes of the velocity. A non-zero amplitude represents changes in the variance. (b) Principal Components (modes in times) with explained variance of the dataset (in %) associated with their respective modes (color coded). High positive values mean the spatial pattern shown in (a) is at its maximum contrast. Negative values mean the spatial pattern is the opposite of the positive spatial pattern. Values around zero indicate that the spatial pattern is flat.
The ROCK-PCA decomposition provides spatiotemporal modes of the variance, each described by complementary spatial maps and an intensity time series. The analysis shows that 78.8% of the dataset’s variance is explained by its first 6 modes.
The spatial patterns of the modes show similarities in active areas of the tributaries. Mode 1 emphasizes western Seward’s lobe and the lower trunk. Mode 2 shows strong variance in the lower trunk, Seward’s lower lobe, and Marvine’s eastern trunk. Mode 3 brings out features in Marvine’s throat and lobe, in Seward’s eastern lobe, and in parts of Seward’s western lobe. Mode 4 distinguishes Marvine’s eastern trunk from the rest of the glacier. Mode 5 draws attention to Agassiz’s lobe and Seward’s lower trunk. Finally, Mode 6 highlights Seward’s lower trunk, nearby small areas, and part of Agassiz’s upper trunk.
Temporally, modes 1, 2, 4, and 5 display yearly cycles. Mode 1 has strong peaks in 1986/87 and 2021. Mode 2 presents a multi-year oscillation with a strong break after 2000. Mode 3 distinctively picks out the unusual pattern of the 1986/87 surge, which also leaks into modes 1, 5 and 6. Modes 4, 5, and 6 do not show clear patterns before 2014, the year after which annual variations appear. The appearance of better resolved annual cycles later in the time series could be a result of better temporal resolution in the recent ITS_LIVE products. Mode 4 also displays a drop in 2013/14. Mode 5 shows negative values for a multi-year period in 1985/86/87 and 2001/02/03. Finally, mode 6 does not display any noticeable pattern similar to the other modes.
We do not show the phase patterns of the ROCK-PCA here, as they are difficult to interpret in the absence of modes with very regular oscillations in time.
3.3.2. Ice flow history
The ITS_LIVE dataset is richer than Landsat’s in terms of data density. Whereas the observed surges in the Landsat mosaics record agree with those from the ITS_LIVE dataset, other events do not fall into a binary surge classification.
We retrieve ice velocities at transects located at the throat of each tributary’s trunk (Figure 4) to assess the temporal variations. For reference, the un-interpolated equivalent can be seen in figure S12. To quantify the propagation of surges down the lobes, we select velocities along each tributary’s central flow line. Only speeds exceeding 4 times the median (see Methods) are plotted (figure S6). Surges tend to propagate far down Seward and Agassiz lobes, while their trunks do not meet the surge requirement to appear on the plot. Marvine’s surges are mostly located in its main branch and rarely propagate to the lobe, except for the 1986/87 surge that reached the terminus.
We identify several surges based on Figures 4 and S2. In 1986 and 1987, all the tributaries experienced fast flow propagating far down their lobes. All of Seward’s lobe was impacted with a particularly strong response from the eastern part. The next surge occurred in 1999 in the western part of Seward’s lobe and propagated to the terminus (Fountain Stream, Figure 1). In 1999, Agassiz started showing a noticeable increase in peak velocities until its 2002 maximum. This signal propagated from the throat to the terminus (figure S6). In 2008/09, Seward surged again in both the eastern and western parts of the lobe. In 2013, Marvine surged, but the signal weakened at the throat, while prior to July, Seward was almost surging. Finally, in 2020/21, Seward surged only in the western part of the lobe to its terminus. This surge is well defined with two peaks (one in 2020, another in 2021), and a slowdown in between.
Figure 4. (a) Seward’s throat transect flow direction variations. The x-axis represents the speed-weighted mean flow position from west to east along the transect (equation 1). The vertical dashed line helps visualize the flow transitions from a southeastern or southwestern direction. (b) Velocities along transects spanning the three tributaries. The x-axis shows position along the transects (blue = Agassiz, orange = Seward, purple = Marvine), and the y-axis shows time (from top to bottom, 1984 to 2022). The vertical white lines separate the tributaries. The dashed black rectangles surround identified surges. (c) Velocities of three points at the center of the tributary’s throat’s transect. The dashed horizontal lines represent 4 times the median of the points’ time series.
Both Seward and Agassiz show a seasonal signal of summer speedup and winter slowdown (Figure 4). However, Seward is prone to surging without noticeable precursor signs, while Agassiz has a gradual buildup during the years preceding a surge. For Marvine, the velocities are very low aside from the 1986/87 and 2013 surges, and seasonal patterns or precursor signs of surges are difficult to observe.
The calculation of Seward’s speed-weighted mean flow direction (Figure 4a) shows that the 1986/87 surge resulted in the most eastward flow diversion of the whole record. The signal displays a noticeable change eastward in 2002/03. 2013 marks a change in the appearance of the time series, probably due to data density. Before 2013, some yearly oscillations are visible, but the poor data resolution makes it difficult to see a structure beyond the annual cycle. After 2013, the annual signal shows a slight jerk eastward in the fall of every year, then a change in direction towards the west until the next Fall.
We investigate when each part of the lobe changed the most throughout the record. For each pixel, we calculate the median monthly speeds and retrieve the date of the fastest month of the time series, under the condition that the year has 5 months or more of surging velocities (figure S7). We can thus identify which surge had the most impact at any given map location. The 1986/87 surge is the strongest on record for the eastern part of Sít’ Tlein’s lobe and low velocity areas. Agassiz’s strongest surge occurred in 2002, while western Seward’s lobe’s strongest surge happened in 2021. Marvine had significant surges both in 1986/87 and in 2013, although the most recent one was limited to the throat. We hypothesize that the 2013 signal at Seward’s terminus reflects erroneous ITS_LIVE velocities caused by thermokarst changes misinterpreted as displacement by the cross-correlation algorithm.
3.4. Summary of surge history
Table 2 summarizes identified surges, periods of confirmed quiescence, and observational ambiguities and data gaps. Data density increased markedly with the start of the satellite era, resulting in higher certainty of identifying surges. The 1980s stand out as an active decade, with surges occurring at all tributaries. Generally, surges repeat at about decadal intervals, but are not synchronized across tributaries.
Table 2. Table summarizing surges for the three tributaries in the past century. Purple represents surges witnessed by the expeditions, green represents surges shown by the Landsat mosaics only, and orange represents surges detected both through the ITS_LIVE dataset and Landsat mosaics. Grey rectangles and dots represent observations of quiescence, while black rectangles and black dots represent a lack of observations. Ambiguous data, in crimson red, means that one method indicates a surge while the other does not. It can also be due to a non-confident observation from a report
3.5. Flux gates
The average influx of ice through all three tributaries into the lobe is 5.1 km3a−1. Divided by the ablation area (2308 km2), this averages to 2.2 m a−1 of ice spread over the lobe. Seward accounts for 86% of the total influx to the piedmont, while Agassiz and Marvine supply 13.65% and 0.35% of ice, respectively.
Years of higher velocities (e.g. 1994, 2012, 2017) at Seward’s throat are not matched at the terminus. 2021 marks Seward’s throat’s maximum flux for the entire record. Although Seward’s terminus’s spike is the 2nd highest of the record in 2021, the year of highest flux through the transect is in 1987. Spikes in fluxes tend to reach the terminus with a 1-year delay from those observed at the flux gate at the throat. The 2021 and 2022 fluxes at Seward’s terminus are equivalent to 116% of the cumulated flux between 2010 and 2020, indicating that surging is essentially the only mechanism supplying ice to the near-terminus region. The highest flux of ice entering Sít’ Tlein’s lobe occurred in 2021, while the 1986 surge produced the 2nd highest peak (Figure 5a). Agassiz and Seward present a strong correlation in their variations, with a varying lag changing from 0 years to 3 years.
Figure 5. (a) Annual fluxes for 4 flux gates (Figure 1) with the calculated propagated errors (vertical red bars). The curve labeled ‘Total influx’ (black) represents the sum of Agassiz, Seward, and Marvine flux gates values for each year. The fluxes at Marvine Glacier throat and terminus, and Agassiz terminus, are too small to be shown on the figure. (b) Scaled fluxes of the Seward terminus and throat. The fluxes have been scaled according to the minimum and maximum values of their time series, for comparative reading. The grey areas spanning the two sub-figures represent identified surge events.
We estimate a total cumulative ice flux to Sít’ Tlein’s lobe of 
$193.5\pm18.5\,\mathrm{km}^{3}$ between 1985 and 2022, which corresponds to 28% of the total 691.6 
$\mathrm{km}^{3}$ of ice contained in the lobe (Tober and others, Reference Tober, Holt, Christoffersen, Truffer, Larsen, Brinkerhoff and Mooneyham2023). With the current yearly average ice influx of 5.1 
$\mathrm{km}^{3}$, it would take 132 years to restock the lobe at its current volume.
Supplementary figure S10 shows the scaled cumulative flux at each point along Seward terminus’s flux gate. It displays a highly asymmetrical profile with four different peaks, three of which are spatially close together. Most of the flux at the terminus reaches the western and central parts of the flux gate. The highest peak coincides with Fountain Stream (Figure 1), while the three grouped peaks represent the area north of Sít’Kagi Lagoon and Backdoor Lake. There is an increase in flux in the eastern part of the flux gate, which corresponds to the area around Malaspina Lake. The locations of these flux peaks are a reflection of the bed topography and show how bed troughs direct ice flow in the lobe.
3.6. Dataset quality
We compare the ITS_LIVE original and interpolated datasets, and assess the quality of the original dataset.
Sít’ Tlein’s piedmont lobe has better data coverage than the trunks (Figure 6a). Generally, fast ice has low data density, while the stagnant moraines have higher data density. The apparent patches are due to the footprint of the images used in ITS_LIVE. Data density from 1984 to 2000 is low (Figure 6b) compared to the twofold data increase in 1999 (launch of Landsat-7). 2013 to 2019 marks the biggest improvement in data density following the launches of Landsat-8 and Sentinel 1 & 2 missions. This explains the over-representation of the latest decade in the dataset.
Figure 6. Dataset quality in space and time. Both figures show how many values occur depending on their location in space, and which years are the most represented. (a) Fraction of valid data points for each pixel compared to the complete dataset (number of 2D slices with any valid entries) (b) Number of valid data pixels as a function of time.
For both velocity components, the spread of scaled differences between interpolated and initial datasets increases as velocities decrease (Figure 7a), while their averages are stable. Therefore, errors do not scale with velocities, and there is no magnitude bias introduced by the interpolation. The algorithm corrects outliers, as shown by the better fit of the Student’s t distribution to the histogram of differences (Figure 7b), compared to the normal distribution fit. Mean differences are 0 m a−1 for both velocity components, confirming the absence of bias from the interpolation. The main flow direction is southwesterly, consistent with glacier flow from its accumulation to its ablation areas (Figure 7c). Including lower velocities lowers the distribution width (turquoise and red distributions).
Figure 7. (a) Plot of scaled mean (solid lines) and standard deviation (faded areas) of the differences between initial and interpolated velocity components Vx (purple) and Vy (orange). Each value is plotted per velocity bins spaced every 100 m a−1. (b) Histograms of the differences between Vx (purple) and Vy (orange) reconstructed and initial values, with fitted normal (continuous lines) and Student’s t (dashed lines) distributions over the histograms. The distributions overlap almost completely, hence why only the orange distribution is visible. (c) Histograms of directions are shown for the reconstructed dataset (turquoise), the initial dataset (red), and the differences between them (gray). Results are shown for various velocity thresholds. If red and turquoise are indistinguishable, it means the two distributions are overlapping.
4. Discussion
Our results show that despite a shared regional setting, each of the main tributaries of Sít’ Tlein displays different dynamics. Albeit commonalities (seasonality, surges), they are generally not synchronized. We thus present the analysis for each tributary separately. We first review the interpolation method and modes analysis of the velocity dataset, then we detail each tributary’s specificity, and finish with a discussion on folded moraine formations and the role of surges in ice resupply at the terminus.
4.1. Modes analysis
Although the ROCK-PCA provides a phase metric for the modes, we exclude it from the analysis, as interpreting the phase of non-periodic modes is complicated and would offer limited insight for our study.
The ROCK-PCA decomposition reveals that Seward’s dynamics dominate the dataset’s variance, while surge signals spread through almost every mode. The fact that different surges are picked up by specific modes shows that their expression in the variance of the dataset differs. The modes’ spatial patterns align with surge footprints, but also reveal variance dipoles (i.e., regions opposed in magnitude of variance). Lobes and trunks represent one pair of data variance, eastern and western Seward’s lobe are another, as well as Seward and the rest of Sít’ Tlein. The Marvine glacier illustrates this, with decoupled trunks and lobe dynamics as seen in the Landsat record, and reflected in the spatial patterns of modes 2, 4, and 5.
The ROCK-PCA can guide the analysis of a dataset by streamlining the identification of its main characteristics. Even without prior knowledge, one could identify key patterns: Seward as the most dynamic glacier, each tributary exhibiting surges, decoupling between western and eastern Seward’s lobe, Marvine’s trunk and lobe behaving as distinct dynamic systems, and the glacier showing annual oscillations in flow dynamics. However, there is no guarantee that a derived mode can be linked to a specific physical process, and the mode analysis should always be done with careful consideration of other spatial and temporal knowledge of system behavior.
4.2. Interpolation quality assessment
The differences between the interpolated dataset and observations do not reveal any clear spatial or temporal structure, indicating a lack of systematic errors (figure S11). The interpolation corrects for outliers by prioritizing solutions that avoid large gradients (in space or time), as shown by the good fit of the Student’s t distribution, which is associated with fewer but more extreme outliers (Figure 7b). Penalizing large gradients through a strong regularization is an informed choice: glacier surface velocities can vary in space and time, but we assume that they do so smoothly. Finally, the constant histogram width regardless of the threshold suggests ITS_LIVE’s errors are independent of velocity magnitude.
4.3. Agassiz Glacier
According to the Landsat mosaics and ITS_LIVE datasets surge definitions, the Agassiz Glacier surged several times, most dramatically in 1986. Its flow is cyclical: it builds up over several years, surges for a year, then slows down. This sequence is unique when compared to other nearby glaciers. Trapridge Glacier’s slow surge is less dramatic and ends more gradually (Frappé and Clarke, Reference Frappé and Clarke2007), while Variegated Glacier shows a more abrupt surge onset despite a build-up phase (Raymond and Harrison, Reference Raymond and Harrison1988) similar to Agassiz’s. The binary surge classification fails in this case, and could be improved upon by considering the spatial extent of elevated velocities (Liu and others, Reference Liu, Enderlin, Bartholomaus, Terleth, Mikesell and Beaud2024).
Another striking feature is how Agassiz’s ice is blocked by Seward’s and forced to take a 90-degree turn to the west, past its throat. Figure 1 shows how a bed trough is aligned with the fastest flow, which is also observed for Seward’s tributary.
4.4. Marvine Glacier
Marvine is the slowest tributary, although it displays the widest range of flow dynamics. Its baseline velocities are almost stagnant in the lobe, while its surges exhibit a large range of magnitudes, scaling from spatially limited to redefining the glacier’s extent in 1986/87.
The 1906, 1935, and 1986/87 surges were spectacular in their reach, magnitude, and duration. They could be rare and related to synchronously occurring eastward surges of Seward. Two out of these three surges coincide with a Seward surge. The 1935 pictures only capture Marvine’s lobe and arguably the eastern portion of Seward’s. Without any cover of Seward’s upper lobe, we have a lack of evidence of either a surge or quiescence there.
Marvine’s main trunk is made of two branches. The eastern one provides most of the ice to the lobe (we use velocities as a proxy for ice supply). This branch sometimes surges beyond its throat and remains repeatedly active, unlike the more stagnant lobe and western branch (figures S8, S9). Surges in one branch don’t typically trigger surges in the other, but both branches show elevated velocities for several years before and during the 1986/87 surge. Marvine might enter a major surge (1906, 1935, 1986/87) only when both branches surge; otherwise, the ice supply might not be enough to activate the stagnant ice dam in the lobe.
4.5. Seward Glacier
Seward’s dynamics dominate the variance of velocity in the piedmont lobe (Figures 3 and 5). Its quasi-decadal surge cycle regularly impacts the texture of the piedmont lobe while heterogeneously mobilizing areas at the terminus. The pre-satellite data show that surges are not a new phenomenon, but only since Landsat images can we properly estimate their impact on the lobe. 1986/87 and 2020/21 surges are the most impactful, considering ice influx into the lobe, impacted area, or impact location. However, their effects seem limited in time, as shown by the post-1986/87 surge Malaspina Lake geometry recovery in barely a decade. Yet, the terminus’s flux profile (figure S10) shows ice resupply to preferential zones. The locations of higher flux coincide with the presence of deep troughs in the glacier bed (Figure 1).
In the past 39 years, Seward’s lobe flow has been highly directional (Figure 4a), which raises the question of the impacts of such behavior. We hypothesize that the direction a surge takes across Seward’s lobe may depend on how ice flow is directed at the time of surge initiation. The 1986 surge was mostly active between June and late September 1986, as shown by images from the same year. The flow direction at this period was directed towards the southeast. The 2020 surge started in spring, when the ice flow was in its south-westernmost direction. Both times, the flow direction followed the apparent seasonal cycle.
4.6. Ice flow interactions in a broader context
All three tributaries exhibit evidence of ice flow interactions at different scales. They range from intra-tributary to inter-tributary, including flow oscillation and glacier redirection (blockage/damming).
We hypothesize that routing of subglacial water causes flow redirection (intra-tributary oscillation). The ice in the throat is heavily crevassed, allowing for a fast transfer of meltwater from the surface of the glacier to its bed. Seasonal meltwater on glaciers induces seasonal drainage system changes, which can impact basal motion through the modification of effective pressure (Kamb, Reference Kamb1987, Truffer and others, Reference Truffer, Harrison and Echelmeyer2000, Hewitt, Reference Hewitt2013, Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013).
The second hypothesis, applicable to both intra- and inter-tributary settings, is that ice flow redirection can be guided by a gradient in compressive strain from the ice damming an area. We speculate that it creates the folded moraines by rerouting the ice flow yearly (see next subsection). This hypothesis can also be used for inter-tributary ice dynamics, as observed for other glaciers (Main and others, Reference Main, Copland, Flowers, Dow, Van Wychen, Samsonov and Kochtitzky2024, Polashenski, Reference Polashenski2024, Van Wychen and others, Reference Van Wychen, Jiskoot, Shannon and Gorwill2025). Seward’s lobe blocks and redirects ice from the other tributaries. We do not have sufficient insight to determine if Marvine or Seward caused significant drag on the other during the 1986/87 surge, although we can observe on satellite images that the two lobes are impacting each other’s shape while surging. Finally, although Seward and Agassiz have a strong ice flux correlation, their dynamics and accumulation area settings are very different. The inactive area between the two lobes is fairly large and does not show movement that could indicate a drag induced by Seward’s western lobe.
Sít’ Tlein does not have wavy moraines typically associated with alternating surges between tributaries (Kotlyakov and others, Reference Kotlyakov, Osipova and Tsvetkov2008, Young and others, Reference Young, Flowers, Jiskoot and Gibson2024). Other studies report interacting tributaries in dendritic glaciers (Kristensen and Benn, Reference Kristensen and Benn2012, Main and others, Reference Main, Copland, Flowers, Dow, Van Wychen, Samsonov and Kochtitzky2024, Young and others, Reference Young, Flowers, Jiskoot and Gibson2024), but Sít’ Tlein might not belong in this category (cf. Jiskoot and others, Reference Jiskoot, Fox and Van Wychen2017). Its size and highly variable flow-affected areas might prevent drag-induced surges most of the time, as illustrated by the 2020/21 surge confined to Seward. No spatially-conjoint surges have been observed between western Seward and Agassiz, possibly due to the slow debris-covered zone between them acting as a cushion that prevents surge propagation. The interface between Seward and Marvine is also made of a slow, debris-covered zone, although it is more active than its Agassiz-Seward counterpart.
4.7. Morphogenesis of folded moraines
Sít’ Tlein’s folded moraines are probably its most striking surface feature (Figures 1, S3, S4, S5, S8). They were previously thought to be generated by surges (Lingle and others, Reference Lingle, Fatland, Voronina, Ahlnaes and Troshina1997, Muskett and others, Reference Muskett, Lingle, Tangborn and Rabus2003), but the oscillation of flow direction suffices to generate them, as shown by the annual Landsat mosaics displaying moraines being folded even during quiescence. Small variations in moraine patterns are subsequently amplified by the laterally divergent flow as the ice spreads into the piedmont lobe. Ramberg Reference Ramberg(1964) used an analogue wax model to recreate the shape of folded moraines without introducing surge-like events. The mechanism inducing the flow variation direction is inherent to the hourglass shape of the glacier: two reservoirs connected by a narrow outlet. Minor variations in the velocity and composition of the material would be enough to create disturbances that induce flow direction oscillation. Surges enhance the process by exaggerating the moraines’ folds and movements, but they do not appear to be essential to their morphogenesis.
4.8. Surge implications for ice resupply
Parts of Seward’s terminus have been stagnant for many decades in places the ice is covered in debris and vegetation, including mature trees. The stacking of moraines is a reflection of ice resupply, which happens primarily through surges, and this stacking is heterogeneous throughout the terminus. While we identified seasonal flow direction changes, predicting surge direction also requires precise timing of surge onset. This information is important for assessing the location and pathways of Sít’ Tlein’s retreat.
Glacier surges most often resupply the area north of Sít’Kagi Lagoon and Backdoor Lake (Figure 1). This leaves other areas without ice supply and subject to thinning and retreat, particularly in the eastern part of the lobe. If a surge delivers 10 years worth of ice but occurs infrequently, retreat becomes inevitable as melt outpaces incoming flux. Although events like the 1986/87 eastern surge can reshape the extent of Sít’ Tlein, the longevity of the effects vary. Malaspina Lake recovered its original shape in a decade, while the terminus of Marvine still bears the marks of the surge.
Given the low elevation of Sít’ Tlein’s lobe, current resupply and mass balance would not permit its formation under today’s climate. The yearly average influx to Sít’ Tlein’s lobe is 2.2 m a−1 of ice resupply averaged over the lobe. Any melt exceeding this amount indicates a disequilibrium for the area. Sharp Reference Sharp(1951) measured ablation rates above 3 m over an entire summer, while we measured 2.9 m of melt between July and August 2021 (averaged over the lobe). The entire lobe is located in the ablation area, the summer melt (after snow melt) thus reflects the annual climatic balance. Our measurements started after snow melt and were terminated before snowfall, therefore, they somewhat underestimate the annual climatic balance. Nonetheless, those numbers exceed the average ice supply and show that the lobe is at a disequilibrium. Modeling work from Brinkerhoff and others Reference Brinkerhoff2025 and several mass balance studies over the area (Muskett and others, Reference Muskett, Lingle, Sauber, Post, Tangborn and Rabus2008, Larsen and others, Reference Larsen, Burgess, Arendt, O’Neel, Johnson and Kienholz2015, Hugonnet and others, Reference Hugonnet2021) are consistent with this assessment.
5. Conclusion
Sít’ Tlein has experienced surges for over a century. Although the record is incomplete, snapshots offered by reports and satellite mosaics allow us to ascertain that the dynamic patterns observed over the last 37 years, the period corresponding to satellite observations, have been occurring before that with a different reach and perhaps frequency.
We showed that approaching a vast dataset like ITS_LIVE requires adapted pre-processing and data analysis methods capable of handling non-linear patterns. The Probabilistic Matrix Factorization interpolation performed very well with such an organized data structure and corrects the original dataset by removing its outliers and performing an informed interpolation across data gaps. The ROCK-PCA enabled us to streamline the dataset analysis by highlighting key features of the spatiotemporal variance. It revealed an organized ice dynamics structure and paired dynamical systems.
We categorized 100 years of surges following methods specific to each dataset, by identifying surge regimes for each tributary, assessing their spatiotemporal distribution, and quantifying their impact on their respective areas. Velocity and flux analysis revealed the Seward tributary as the primary driver of ice dynamics in Sít’ Tlein and clarified that seasonal ice flow redirection explains the formation of banded moraines across the lobe. The surges, previously thought to be solely responsible for their formation (Muskett and others, Reference Muskett, Lingle, Sauber, Post, Tangborn and Rabus2008), might speed up and enhance the process. Our analysis is consistent with a lab experiment (Ramberg, Reference Ramberg1964) showing that moraine bands can occur without surges.
All three tributaries show repeated surges with different characteristics, and generally occurring at different times. However, they can coincide on occasion, as last happened in 1986/87. In the current climate, a surge is the only mechanism to activate the stagnant terminal ice that is buried under sediment or vegetation. Surges from the Seward Glacier have a localized influence on the lobe, and the main direction of flow appears to be related to the timing of surge onset.
Overall, on a decadal scale, surges fail to offset the imbalance induced by high surface melt. The challenging detail for models is to capture the spatial distribution of a given surge event, which can have a large temporary influence on a very localized part of the terminus.
Data availability statement
Code, data, and supplementary figures can be found at https://doi.org/10.5281/zenodo.15652862.
Acknowledgements
This study is part of a collaborative research project funded by the National Science Foundation (NSF) and the National Park Service (NPS). VDC, MT, MF, and CFL were supported by the NSF Grant #1929566. DB was supported by the NSF Grant #1929718. MSC, MD, BST, and JWH were supported by the NSF Grant #1929577). MGL was supported by Wrangell-St. Elias National Park and Preserve, and by the PMIS project 230360. We thank Gregoire Guillet, Emma Marshall, Kristin Timm, Anna Thompson, Tyler Kuehn, Russ Mitchell, Natalie Wagner, and Elvie Underwood for discussions, manuscript review, and advice during the production of this study. We gratefully acknowledge the constructive feedback provided by Alex Gardner, one anonymous reviewer, and Hester Jiskoot, Scientific Editor. Their insightful comments and suggestions significantly improved the quality and clarity of this manuscript.
Author contributions
MT, DJB, MF, JWH, ML and CFL developed the project and acquired financial support. CFL, JWH, BST, MSC, and MT acquired and processed bed altimetry data. BST and DJB developed and analyzed bed altimetry data and products. MF provided velocity data and guidance on the surface velocity product. DJB wrote the numerical interpolation method. VDC, MT, and MC developed surface velocity products analysis. VDC developed the paper concept and wrote the original draft. MT, DJB, MGL, and MD contributed to the revisions. All co-authors contributed to discussions and development that led to this manuscript.
Open access
