Abrupt and sustained increase in glacier surface velocity

Abrupt and large increases in glacier flow velocity, sustained for at least several months, are assumed to indicate surging (see Section 3.2.2 for more information), which we typically identify here as accelerations in the flow regime of a glacier. Although many existing definitions involve a surface velocity of at least ten times that observed during the quiescent phase, here we follow the reasoning of Guillet and others Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch(2022) and use a lower criterion of two times. In addition, we define a speed-up event as a surge candidate if the velocity is maintained above the threshold for at least four consecutive months (120 days), as it prevents the false identification of seasonal speed-ups as surge events. In the next section, we provide more information on how this and the subsequent proposed thresholds are used to detect surges.

Substantial and spatially concentrated thickening near the glacier terminus

We generally consider a glacier to be a candidate for surging if it presents substantial and spatially concentrated changes in surface elevation (over 1-10 years) at lower elevations (near the glacier terminus; see Section 3.2.3 for more details), as this deviates from recent trends in global glacier thinning and retreat (Hugonnet and others, Reference Hugonnet2021). Thus, our objective is to identify glaciers that exhibited a substantial and widespread surface elevation gain (dynamical thickening) over the receiving zone. We consider clear statistical breakpoints (defined by a unitless changepoint score) and positive trends in surface elevation time series over glaciers to be the result of surge-induced dynamical thickening.

Abrupt changes in glacier surface crevassing

Surges typically result in intense and widespread surface crevassing and changes in the crevasse patterns at the glacier surface (Truffer and others, Reference Truffer, Haeberli and Whiteman2021; Guillet and others, Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022; Kääb and others, Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi2023). Such changes are identifiable via proxy as abrupt changepoints in glacier SAR backscatter trend time series (Leclercq and others, Reference Leclercq, Kääb and Altena2021; Kääb and others, Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi2023). Again, we consider clear statistical breakpoints (defined by a unitless changepoint score) and positive changes in trend of SAR backscatter time series over glaciers to be indicative of intense surface crevassing associated with a surge.

In the present surge identification scheme, we require at least two of the aforementioned criteria to be met in order to verify a detected event as a surge. As an example, an event flagged as a surge candidate through extended positive surface elevation but for which no signal is detected in either surface velocity or crevassing changes will not be classified as a surge. Adopting a multi-criteria approach makes the framework conservative, but this is necessary as it reduces the occurrence of false positives, such as significant accelerations in ice flow resulting from dynamic adjustment of tidewater glaciers to major calving events or frontal collapses (De Angelis and Skvarca, Reference De Angelis and Skvarca2003; Benn and others, Reference Benn, Warren and Mottram2007; Benn and Evans, Reference Benn and Evans2014). Although surges often result in an advance of the glacier terminus, not all surge-type glaciers show a terminal advance during the active phase (Murray and others, Reference Murray, Dowdeswell, Drewry and Frearson1998; Benn and Evans, Reference Benn and Evans2014; Guillet and others, Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022), and thus we do not use glacier advance as a criterion for surge identification.

3.1. Data

3.2. Detection of surge events

3.2.1. Data sampling strategy

Based on the spatial resolution of the ITS_LIVE and DEM products (≈120m), we limited our analysis to RGI glacier polygons with flowlines long enough to be sampled at ten equally spaced vertices, i.e longer than 20 pixels and 2.4 km. This results in a total of 32133 glaciers with 54178 flowlines, and thus 541780 points to automatically evaluate for surging (Figure 1).

Figure 1. Example of the data sampling strategy applied to Khurdopin glacier, Karakoram. (a) The map shows the RGI7.0 outline of the glacier (blue area) as well as the different flowlines (greyed lines), with the main flowline highlighted in red. Black squares are vertices at which each dataset is sampled. (b) ITS_LIVE surface velocity estimates. (c, d) surface elevation change, and (e, f) SAR backscatter. The scatterpoints in c and e represent the pre-Gaussian process regression time series. The mean of each Gaussian process regression is presented a solid blue line, while the shaded blue area is the 95% credible interval. (d, f) Changepoint detection scores for the surface elevation change detection (d) and SAR backscatter (f) frameworks, see section 3.2.3. Note the variable x-axis ranges as the datasets span different periods. The synchronously detected surge in all three datasets is shaded in yellow.

3.2.2. Surface velocity

Automating the detection of glacier surge events within a time series of glacier surface velocity data faces challenges when relying solely on conventional signal processing techniques and time series analysis methodologies. Multiple sensors are used to prepare the ITS_LIVE velocity time series, which improves spatio-temporal sampling for all glaciers, but data quality issues arise due to the limitations of individual sensors, including cloud cover and illumination conditions for optical instruments. Such constraints lead to velocity data that have data gaps in time and space, with variable uncertainty, which presents difficulties for conventional time series analysis methods meant for regularly sampled data. In addition, glacier surges demonstrate intricate dynamics, showing different velocity patterns with diverse durations, magnitudes and temporal changes. Therefore, we chose a statistical modelling approach to automatically detect surges in the irregular surface velocity time series.

We apply the following framework independently to each point $s(x, y,z)$ along the flowline at which the velocity time series is sampled. We only consider surface velocity estimates computed with a temporal baseline of 6 to 110 days. We then normalise the surface velocity time series $U_s(t)$ using the mode of surface velocity magnitude for the full 2000–24 time period. We assume that the observed surface velocity $U_s(t)$ is the result of two distinct sources: the baseline surface velocity of the glacier $U_b(t)$ and an ‘excess’ surface velocity term $U_e(t)$ which potentially includes dynamic ice flow instabilities:

(1)\begin{equation} U_s(s, t) = U_b(s, t) + U_e(s, t) \end{equation}

Our surge identification scheme relies on estimating the excess velocity $U_e(t)$ at a given time t by modelling $U_b(t)$.

Several approaches to model time-varying glacier surface velocity have been proposed. However, these efforts rely on a periodic structure of the velocity signal (sine or cosine), often use higher-order polynomials to interpolate missing values and/or require continuous temporal data coverage (e.g. Vijay and others, Reference Vijay, Khan, Kusk, Solgaard, Moon and Bjørk2019; Greene and others, Reference Greene, Gardner and Andrews2020; Charrier and others, Reference Charrier, Yan, Trouvé, Koeniguer, Mouginot and Millan2022). Here, we model $U_b(s, t)$ for each point s of the flowline as a Gaussian process:

(2)\begin{equation} U_b(t) \sim \mathcal{GP}(\mu_{U}(t), k(t, t^{\prime})) \end{equation}

where $\mu_{U}(t)$ is the time-dependent mean and $k(t, t^{\prime})$ is the kernel function that captures the temporal covariance between data points.

Gaussian process mean

The Gaussian process proposed here acts as a surrogate model, or emulator, for what can be considered the baseline surface velocity of each glacier. Thus, its mean represents the most likely baseline velocity value at any given time t and is derived directly from the data. Since the distributions of glacier surface velocity values are non-symmetric and present heavy positive tails, the median value for the full 2000–24 period is not a robust estimate of the most probable baseline velocity value. Here we use the mode of the distribution of the surface velocity values as the mean for the Gaussian process. To calculate the mode of the distribution, we first used a Gaussian kernel to estimate the probability density function of the whole sample of glacier surface velocities. The mode was then estimated as the minimum of the negative probability density function. It is thus assumed that surges are represented by a relatively small proportion of the distribution and hence that the glacier is, more often than not, in quiescence during the 2000–24 period.

Gaussian process covariance

In practice, we model the changes in glacier velocity at a given point $s(x,y,z)$ and time t as a stochastically driven damped simple harmonic oscillator (SHO). This allows us to more intuitively capture potential variations in the frequency of baseline glacier surface velocity from one year to another. We use the celerite2 (Foreman-Mackey and others, Reference Foreman-Mackey, Agol, Ambikasaran and Angus2017) Python implementation of Gaussian process models with SHO kernel terms defined by the authors through the power spectral density (ω):

(3)\begin{equation} S(\omega)=\sqrt{\frac{2}{\pi}} \frac{S_0 \omega_0^4}{\left(\omega^2-\omega_0^2\right)^2+\omega_0^2 \omega^2 / Q^2} \end{equation}

with ω ₀ being the undamped angular frequency and Q the quality factor describing the resonance of the harmonic oscillator. celerite2 proposes an alternative parametrisation of the SHO term and further allows the user to define more intuitive kernel hyperparameters of the form:

(4)\begin{align} & \rho=2 \pi / \omega_0, \text{the undamped period of the oscillator (in days)} \end{align}

(5)\begin{align} & \tau=2 Q / \omega_0, \text{the damping timescale of the process (in days)} \end{align}

(6)\begin{align} & \sigma=\sqrt{S_0 \omega_0 Q}, \text{the standard deviation of the process (normalised)} \end{align}

We thus define the covariance function as a sum of two SHO terms capturing velocity variations through an additive model as annual and sub-annual (typically 4 months) periodic deviations:

(7)\begin{align} k(t, t^{\prime}) = & \underbrace{\mathrm{SHO}(\rho=365, \tau=365*2, \sigma=\mathrm{IPR}(U_s))}_{\mathrm{Annual~variations~term}}+ \\& \nonumber \underbrace{\mathrm{SHO}(\rho=365/4, \tau=365, \sigma=\mathrm{IPR}(U_s))}_{\mathrm{Sub-annual~variations~term}} \end{align}

where $\mathrm{IPR}$ represents the interpercentile range (between the 5th and 95th percentiles) of the entire available sample of glacier surface velocity measurements, ${U_s}$.

Using our fully-specified Gaussian process model(Equations 2 and 7), we can now derive an estimate of $U_b(s, t)$ through Gaussian Process regression for each sampling point s, and hence rewrite Equation 1 to estimate the excess velocity at a given time t (Figure 2):

(8)\begin{equation} U_e(s, t) = U_s(s, t) - U_b(s, t) \end{equation}

Figure 2. Time series of ITS_LIVE glacier surface velocity estimates (black dots) for a vertex along the main flowline of selected glaciers. For each glacier the top plot presents ITS_LIVE glacier surface velocity estimate (black dots) and the computed baseline velocity (solid coloured curves, mean of the Gaussian process prediction), the shaded area is the 90% confidence interval. The bottom one presents the excess velocity estimates, i.e. the residuals from the Gaussian process regression (black dots), and automatically identified surge events (shaded regions). Columbia Glacier is not a surging glacier and is presently used to demonstrate the proficiency of the Gaussian process in emulating glacier surface velocities.

3.2.3. Surface elevation changes and SAR backscatter

We aim to identify glaciers that exhibited a substantial and widespread reorganisation of their surface as clear identifiers of a surge. This entails drastic increases in surface elevation gain (thickening) in their receiving zone in contrast to the expected glacier thinning trends (Hugonnet and others, Reference Hugonnet2021), or drastic changes in surface crevassing, as visible on backscatter time series (Leclercq and others, Reference Leclercq, Kääb and Altena2021; Kääb and others, Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi2023).

In practice, we are specifically aiming at detecting sharp transitions in trends in surface elevation change rates, as well as backscatter, as shown on Figure 1. Surface elevation change and SAR backscatter time series present similar structures in that they can be modelled as a trend term and one, or several, periodic terms (shown for elevation change in Hugonnet and others, Reference Hugonnet2021). However, they both present individual challenges which prevent the sole use of standard signal processing methodologies. Surface elevation change time series are plagued by their overall sparsity and irregular sampling, and SAR backscatter time series contain seasonal variations (largely due to precipitation and melt) with greater amplitudes than the typical change in trend from a surge that we wish to detect. Both existing data gaps, and high-amplitude periodic oscillations hamper the use of conventional trend change detection methods. We therefore set up a multistep framework to detect positive break points from both datasets. Similarly to section 3.2.2, the following framework is identical for each vertex $s(x, y,z)$ along the flowline at which the surface elevation and backscatter time series are sampled.

Gaussian process regression and signal resampling

We begin by modelling each time series y(t) as a realisation of a Gaussian process:

(9)\begin{equation} y(t) \sim \mathcal{GP}\left(\mu_{y}(t), k(t, t')\right) \end{equation}

We again rely on celerite2’s SHOTerm (Eq. 3) to capture both the periodic structure of the signal, as well as longer-term (typically close to 5 years) trend:

(10)\begin{align} k(t, t^{\prime}) =\ & \underbrace{\mathrm{SHO}(\rho=365, \tau=365*\theta_{1}, \sigma=\mathrm{IPR}(y))}_{\mathrm{Annual~variations~term}}+ \\ & \nonumber \underbrace{\mathrm{SHO}(\rho=365*\theta_{2}, \tau=365, \sigma=\mathrm{IPR}(y))}_{\mathrm{Long-term~trend~term}} \end{align}

where IPR represents the interpercentile range (between the 5th and 95th percentiles) of the entire available sample of either glacier surface elevation or backscatter measurements, y. θ ₁ and θ ₂ are parameters used to capture the wide diversity of surface elevation and backscatter patterns observed globally. They are fitted to each individual time series (from each flowline vertex) through maximum likelihood estimation using an L-BFGS scheme. Examples of realisations of individual Gaussian processes used for regression and interpolation of surface elevation and backscatter time series are given in panels b and d of Figure 1.

After fitting the Gaussian process model, each signal is resampled uniformly to ensure consistent spacing for subsequent analyses. This step minimises the risk of aliasing and ensures that our changepoint detection operates on an evenly sampled input signal.

Wavelet decomposition

In order to identify changes in signal trend, each resampled GP estimate needs to be decomposed into multiple time scales, accounting for trend and seasonality. In practice, we here want to isolate the trend in either measurement by filtering out any high-amplitude seasonal signal. Since the resampled GP signal is not purely periodic, we perform signal decomposition using discrete wavelet transform with Daubechies 4 (db4) wavelets. We choose Daubechies wavelets as they are effective for capturing both transient features and structural changes in signals. The wavelet transform represents the signal y(t) as a sum of approximations and details (Mallat, Reference Mallat1999):

(11)\begin{equation} y(t) = \sum_{j,k} c_{j,k} {\psi}_{j,k}(t), \end{equation}

where $ {\psi}_{j,k}(t) $ are wavelet basis functions indexed by scale j and position k, and $ c_{j,k} $ are the corresponding wavelet coefficients. An example of wavelet deseasonalisation is given in panels b and d of Figure 1.

Penalised changepoint detection

We consider $y(t=k)$ to be a changepoint if at time $ k \in \mathcal{T} = \{t_1, \dots, t_K \}_{K \leq T} $, there is a change in the mean value of the wavelet-decomposed GP estimate. We note $\mathcal{T}$ a partition of the input signal y(t) with K + 1 subsequences, where K is the number of changepoints. In our case, K is unknown but is assumed to be in the range [0, 1] and is found through a changepoint detection procedure.

A changepoint detection procedure aims to find the optimal segmentation $\hat{\mathcal{T}}$ by minimising the quantitative criterion $ V(\mathcal{T}, y) = \sum_{k=0}^{K} c\left(y_{t_k \ldots t_{k+1}}\right) $, with $c\left(\cdot\right)$ a cost function measuring, for each point $y(t=k)$ a changepoint partitioning precision; i.e. how well the signal is divided into segments such that each segment exhibits consistent statistical properties, with changepoints marking significant shifts in data characteristics. Finding $\hat{\mathcal{T}}$ is a discrete optimisation problem, with, in our case, an unknown number of K changepoints and written as follows:

(12)\begin{equation} \min_{\mathcal{T}} V(\mathcal{T}, y) + \text{pen}(\mathcal{T}) = \min_{\mathcal{T}} \tilde{V}(\mathcal{T}, y) \end{equation}

with $\text{pen}(\mathcal{T})$ constraining the number of detected changepoints by effectively penalising models with a higher number of changepoints. The proposed changepoint detection procedure is implemented using ruptures (Truong and others, Reference Truong, Oudre and Vayatis2020) and made up of 3 main parts:

• • (i) A search method or detection algorithm, to estimate $\hat{\mathcal{T}}$: we use the window search detection method, consisting of computing the discrepancy between two adjacent sliding windows along y(t) (Truong and others, Reference Truong, Oudre and Vayatis2020).

• • (ii) A constraint on the number of changepoints. In practice, we find that a penalisation parameter of 70 ensures the detection of a minimal number of changepoints in most cases.

• • (iii) A cost function $c\left(\cdot\right)$. We use an ensemble of three different cost functions: a rank-based cost function (Lung-Yut-Fong and others, Reference Lung-Yut-Fong, Lévy-Leduc and Cappé2011), a Mahalanobis-type metric (Truong and others, Reference Truong, Oudre and Vayatis2019), and a kernelised mean change (Arlot and others, Reference Arlot, Celisse and Harchaoui2019).

Similarly to Katser and others Reference Katser, Kozitsin, Lobachev and Maksimov(2021) we then rely on the aggregation of all models, combining different cost functions, to produce an aggregated cost; a method also known as ‘mixture of experts’ by the machine learning community. For each model and cost function, a score quantifies how well a particular segmentation fits the data. Scores are then scaled between 0 and 1 aggregated by taking the pointwise minimum. Since we are only interested in increases in the trend of the surface elevation change and SAR backscatter signals, we set the score of changepoints detected with a negative trend gradient to 0. Finally, vertices with a combined changepoint score greater than 0.6 are considered candidates for surging.

3.2.4. Combination of criteria and derivation of a surge inventory

Each of the aforementioned detection steps and data types leads to a different subset of candidate surge events (Figure 3). The subsets of potential surge-type glaciers are then compared based on their RGI v7.0 identification number. Sub-inventories are then combined as follows: glaciers detected by both our surface elevation change and surface velocity frameworks, or both the surface elevation change and SAR backscatter frameworks are considered to be surging (Figure 3). Since surge-induced changes in SAR backscatter are the result of changes in surface velocity, we do not consider the combination between our surface velocity and SAR backscatter frameworks to be indicative of surging, as other phenomena could be falsely identified as surges.

Figure 3. Dataflow diagram overview of the surge event detection scheme, showing the flow of information between input and outputs (blue blocks), internally used data (white blocks), and processing steps (rounded blocks). Green blocks represent the data-specific surge detection thresholds described in Section 2. Blue and red arrows represent the dataflow of surface elevation and SAR backscatter measurements respectively. The blocks defined as AND and OR are logical gates, stipulating that surges have to be detected by either i) the surface elevation change and velocities detection schemes or ii) the surface elevation change and backscatter detection schemes.

3.2.5. Validation and quality control

Validating the semi-automated detection of glacier surges is not straightforward due to existing inventories being derived from inconsistent methods and data sets. The most comprehensive inventory of surge-type glaciers is the one presented by Sevestre and Benn Reference Sevestre and Benn(2015). However, we argue that, due to the mismatch in the 2000–24 period of our inventory and the 1861–2013 period of the Sevestre and Benn Reference Sevestre and Benn(2015) inventory, as well as the conceptual differences in the identification of surges, a detailed comparison between the two inventories would be inconclusive. The most recent global inventory of glacier surges was presented by Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) and relied on the visual interpretation of Sentinel-1 backscatter signals, following the methods described in Leclercq and others Reference Leclercq, Kääb and Altena(2021).

We rely on the inventory of Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) as a validation data set and assess the capability of our detection framework quantitatively through the use of the Jaccard index, also known as the Intersection over Union, a similarity measure commonly used in object detection and expressed as follows:

(13)\begin{equation} J(A, B)=\frac{|A \cap B|}{|A \cup B|}=\frac{|A \cap B|}{|A|+|B|-|A \cap B|} \end{equation}

where A and B represent two different sets of glaciers: A, a validation set of surge-type glaciers identified and B, the glaciers with semi-automatically detected surge events identified in this study. $|A \,\cap\, B|$ thus represents the total surface area (in km²) covered by the glaciers identified in both A and B, while $|A \cup B|$ is the surface area covered by the sum of A and B. The Jaccard index is thus bounded between 0 and 1, with a value of 0 when the sets are completely dissimilar and 1 if the sets are identical.

In addition to the proposed validation experiment, each identified glacier has been manually checked following the criteria proposed in Section 2. In the final product, potential misidentifications have been culled.

4.1. Global distribution of glaciers with active surges between 2000 and 2024

Of the ≈ 32000 glaciers analysed, we detected 246 glaciers with active surges over the 2000–24 period (Figure 4). Most glaciers surged only once, while five glaciers surged two or more times between 2000 and 2024, and one glacier surged four times (Sit’Kusa/Turner glacier, Alaska). Compared to version 7 of the RGI, in Alaska and the Yukon, we identify one surge on a glacier that was previously classified as ‘No evidence’ (Nadina Glacier), as well as 4 previously classified as ‘Probable’ surge-type (Martin River, Marvine, Valerie and Ferris glaciers). In Svalbard, we note a surge of Lilliehöökbreen, previously classified as having no evidence of surge-type behaviour. We further identified nascent surges at Nordsysselbreen, Aavatsmarkbreen and Doktorbreen, all classified as ‘Probable’ surge type. Similarly, we identify a surge on Borebreen, classified as ‘Possible’ surge-type. We finally report the onset of surges at Deltabreen and Seftrombreen which, according to Sevestre and Benn Reference Sevestre and Benn(2015), have not surged in the 20th century, effectively displaying quiescence for the past 165 and 128 years, respectively. All other identified surges in our global sample affect previously identified surge-type glaciers.

Figure 4. Global distribution of 246 glaciers with active surges between 2000 and 2024 identified with our methodology. Prominent but well-known clusters of surge-type glaciers are evident in High Mountain Asia and the Arctic Ring. N refers to the number of surge-type glaciers detected.

Glaciers with active surges during the 2000–24 period represent around 1% of all glaciers studied and 0.1% of all glaciers in RGI V7.0. However, the total surface area of these glaciers with active surges is approximately 38000 km², which corresponds to approximately 5.5% of the worldwide glacierised area outside of the polar ice sheets.

The spatial distribution of surge events is consistent with previous assessments. Approximately 167 (68%) glaciers with identified surges are located in High Mountain Asia, with 78 (32%) in the circum-Arctic region (the Arctic Ring). Within the Arctic Ring, Alaska-Yukon and Svalbard-Russian Arctic had the highest number of glaciers with active surges, with estimates of 28 and 24, respectively. Our method identified 16 active surge-type glaciers in Greenland, 12 in the Canadian Arctic and one in the Andes. No surge events were detected during the 2000–24 period in Iceland. A near or total absence of surges in Iceland is expected (Kääb and others, Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi2023), as most known surge-type glaciers in Iceland surged during the 1990s (Hannesdóttir and others, Reference Hannesdóttir2020) and are therefore likely to be quiescent (Björnsson and others, Reference Björnsson, Pálsson, Sigurdsson and Flowers2003) or senescent (Benn and others, Reference Benn, Hewitt and Luckman2023).

4.2. Climatic controls on the distribution of surge-type glaciers

Surge-type glaciers are found within a specific climatic envelope where they are accompanied by non-surge-type glaciers; beyond this envelope, only non-surge-type glaciers are typically observed. Figure 5 shows the distribution of non-surge-type and surge-type glaciers in relation to the median temperature and the median total annual precipitation from ERA5-Land data for the period between 1990 and 2024. The surge-type glacier climatic envelope typically encompasses a median annual temperature of -1 to $-20^{\circ}\,\mathrm{C}$ and a median total precipitation of 100 to 1200 mm a⁻¹ (Figure 5, panel a).

Figure 5. Median temperature and total precipitation for surge-type (orange) and non-surge-type glaciers (blue). (a, b, c) Due to the disparity in sample size between surge-type (246) and non-surge-type glaciers (more than 270000), the median temperature/median total precipitation relationship is represented using a kernel density estimate, all levels represent lines of probability, or density of the 2D distributions: 20%, 40%, 60%, 90% and 99%. (a) Median total annual precipitation versus median annual temperature. (b) Median total winter precipitation versus median winter temperature. (c) Median total summer precipitation versus median summer temperature. Distribution of surge-type glaciers from this study (red scatter plots) and using the glaciers identified bySevestre and Benn (Reference Sevestre and Benn2015); Falaschi and others (Reference Falaschi, Bolch, Lenzano, Tadono, Lo Vecchio and Lenzano2018); Guillet and others (Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022); Lovell and others (Reference Lovell, Carrivick, King, Yde, Boston and Małecki2023) and Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) (orange kernel density estimate).

Seasonal climate data reveal additional insights on surge-type glacier distributions (Figure 5, panels b and c). Surge-type glaciers typically occur over a broad range of median total winter (October–April for the Northern Hemisphere, June–September for the Southern Hemisphere) precipitation (20 to 1000 mm a⁻¹) and median winter temperature (-10 to $-30^{\circ}\,\mathrm{C}$), while they are still absent from either end of the ranges for the larger population of global glaciers. The occurrence of surge-type glaciers is firmly bound by median summer temperature, as surge-type glaciers are almost non-existent in regions where the median summer temperature exceeds $3^{\circ}\,\mathrm{C}$. Similarly, 90% of the surge-type glacier population lies within a median total summer precipitation range of 100-1000 mm a⁻¹.

The clearest climatic relationship used to distinguish surge-type glaciers is between the median summer temperature and the median total winter precipitation (Figure 6). Surge-type glaciers are absent beyond median summer temperatures higher than $3^{\circ}\,\mathrm{C}$ and less than $-10^{\circ}\,\mathrm{C}$. Similarly, they are seldom present in regions where the median total winter precipitation exceeds 1000 mm a⁻¹.

The characteristics of the climatic envelope established by our analyses are similar to those of Sevestre and Benn Reference Sevestre and Benn(2015) and in very strong agreement with the envelope described by Lovell and others Reference Lovell, Carrivick, King, Yde, Boston and Małecki(2023) using ERA5-Land. Using the glaciers identified by Sevestre and Benn (Reference Sevestre and Benn2015); Falaschi and others (Reference Falaschi, Bolch, Lenzano, Tadono, Lo Vecchio and Lenzano2018); Guillet and others (Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022); Lovell and others (Reference Lovell, Carrivick, King, Yde, Boston and Małecki2023) and Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023), we calculated an updated climate envelope using the more modern ERA5-Land products, which very strongly concurs with our results (Figures 6 and 5).

Figure 6. Median summer temperature and median total winter precipitation for surging and non-surge-type glaciers. Distribution of surge-type glaciers from this study (red scatter plots) and using the glaciers identified by Sevestre and Benn (Reference Sevestre and Benn2015); Falaschi and others (Reference Falaschi, Bolch, Lenzano, Tadono, Lo Vecchio and Lenzano2018); Guillet and others (Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022); Lovell and others (Reference Lovell, Carrivick, King, Yde, Boston and Małecki2023) and Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) (orange kernel density estimate). Lines represent different density levels: 30%, 50%, 70%, 90% and 99%. Also note the difference in median summer temperature between surge-type and non-surge type glaciers.

We finally focus on further validating the proposed climatic envelope. We here aim at avoiding sampling bias in climatic distributions, resulting from sample size differences between non-surge-type (around 270000 glaciers) and surge-type glaciers. For the latter, we use both the 246 surge-type glaciers from this study, as well as the composite inventory formed from the datasets of Sevestre and Benn (Reference Sevestre and Benn2015); Falaschi and others (Reference Falaschi, Bolch, Lenzano, Tadono, Lo Vecchio and Lenzano2018); Guillet and others (Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022); Lovell and others (Reference Lovell, Carrivick, King, Yde, Boston and Małecki2023) and Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023), containing close to 2200 glaciers. From each sample of glaciers, we obtain continuous estimate probability density functions (PDFs) on the number of glaciers in a given temperature/precipitation range by kernel density estimation, using a Gaussian kernel. Figure 7 represents 200000 samples drawn from each PDF and highlights the clustering of surge-type glaciers within the defined climatic envelope. The samples derived from the proposed inventory (Figure 7, panel a) show a greater variance compared to those from the composite inventory (Figure 7, panel b), which results from initial sample size disparity.

Figure 7. Median summer temperature and median total winter precipitation for (a) surge-type glaciers from this study, (b) the surge-type glaciers from by Sevestre and Benn (Reference Sevestre and Benn2015); Falaschi and others (Reference Falaschi, Bolch, Lenzano, Tadono, Lo Vecchio and Lenzano2018); Guillet and others (Reference Guillet, King, Lv, Ghuffar, Benn, Quincey and Bolch2022); Lovell and others (Reference Lovell, Carrivick, King, Yde, Boston and Małecki2023) and Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) and (c) non-surge-type glaciers from the RGI v7.0. Note the important clustering of surge-type glaciers within the defined climatic envelope in a and b. Non-surge type glaciers are more uniformly distributed over the temperature/precipitation spectrum.

4.3. Surge event statistics

The following section describes surge events statistics that are solely derived from velocity time series. Due to the data quality issues described in Section 4.6.1, velocity time series for glaciers in the polar regions lack the temporal resolution necessary to clearly identify surge onset and termination dates. Consequently, we do not consider the surges identified in the Canadian Arctic and the periphery of Greenland within the present section.

The peak surface velocity is defined here as the maximum velocity measured over all sampled ITS_LIVE time series. Similarly, the surge duration is taken as the longest time interval for which the glacier is detected to be actively surging, for all individual surges and across all flowline vertices. The distributions of the peak surge velocity and duration for 193 surge events are presented in Figure 8.

Figure 8. Empirical cumulative distributions of (a) peak surface velocity and (b) surge duration, for the full sample of detected surge events. Regional distributions of (c) peak surface velocity and (d) surge duration for Alaska-Yukon, High Mountain Asia and Svalbard-Russian Arctic. Vertical dashed lines are the median of each distribution.

4.4. Regional differences in surge dynamics: on the Svalbard vs Alaska-type classification

The hydrological switch (Alaska-type surges) hypothesis posits relatively short active phases, typically lasting between 1 and 3 years and reaching a peak velocity of 10 to 100 metres per day (Murray and others, Reference Murray, Strozzi, Luckman, Jiskoot and Christakos2003). Conversely, the thermal switch hypothesis supposes longer durations than Alaska-type surges, typically lasting up to more than 10 years and reaching peak velocity between 1 and 15 metres per day (Murray and others, Reference Murray, Strozzi, Luckman, Jiskoot and Christakos2003). We therefore test for the existence of statistical differences between the distributions of peak surface velocity and surge duration between different surge clusters, with a focus on Alaska-Yukon and Svalbard-Russian Arctic.

The distribution of peak velocity during surge events appears similar in both regions, with surge events reaching more than 20 m day⁻¹ of peak velocity and both distributions displaying a median close to 6 m day⁻¹, with an average spread of 5 m day⁻¹. More quantitatively, a two-sample Kolmogorov-Smirnov test comparing the distributions of peak surface velocity for surges in Alaska-Yukon and Svalbard-Russian Arctic return a statistic of 0.18 and a p-value of 0.81 indicating that the two distributions are not significantly different. Our results therefore do not statistically support a Svalbard-vs-Alaska-type classification. This finding corroborates the idea of dynamical unity between what was previously considered two ends of a behavioural spectrum resulting from distinct processes (Jiskoot, Reference Jiskoot, Singh, Singh and Haritashya2011; Sevestre and Benn, Reference Sevestre and Benn2015; Herreid and Truffer, Reference Herreid and Truffer2016; Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a).

A closer investigation of the differences between glaciers in Alaska-Yukon, Svalbard-Russian Arctic and High Mountain Asia further reveals the breadth of dynamical behaviours displayed by surges.

First, it is worth noting that we expect the distributions described in Section 4.3 to be updated by future work as, out of the 31 surges detected in Svalbard-Russian Arctic, the onset and termination date could only be estimated for 18 of them. Most ongoing surges have already lasted longer than the estimated median surge duration in the Svalbard-Russian Arctic and are likely to further positively skew the distribution. As an example, the surge of the Nathorstbreen Glacier system, the onset of which is invisible on ITS_LIVE, is believed to have occurred around 2006 (Sund and others, Reference Sund, Lauknes and Eiken2014). Similarly, the surges of Austfonna Basin-3 and Scheelebreen, which started in 2012 and 2021, respectively, are still ongoing.

Second, to understand the regional differences in surge dynamics, we find added value in investigating the velocity signals from individual surges. Figure 9 shows the similarities between the surface velocity signals for Alaska-Yukon and High Mountain Asia, as well as how they differ from surges in the Svalbard-Russian Arctic. The acceleration phases of all the surges presented in Figure 9 first shows relatively low linear increase in surface flow velocity, before a rapid switch to a quasi-exponential regime roughly 500 days prior to reaching peak value (Figure 9 panels 13, 26, 39). All three clusters further seem to reach around 10% of maximum velocity value. Individual differences are apparent in the transition between slow linear flow acceleration and active surge phase. In the selected surges from Alaska-Yukon, few glaciers, apart from Donjek and Lowell (Figure 9 3, 8), show notable increases in surface flow earlier than 500 days prior to the peak. This dynamic is similar in Svalbard-Russian Arctic, where only Monacobreen ice cap (Figure 9, 18) shows a gradual increase in ice flow before a dramatic spike in velocity, marking the transition to an active surge phase.

Figure 9. Example of ITS_LIVE surface velocity times for surge events in (1-12) Alaska-Yukon, (14-25) Svalbard-Russian Arctic and (27-38) High Mountain Asia. (13, 26, 39) Regional stacks of surface velocity time series. Bold lines represent stack median and shaded areas cover the minimum to maximum range. Subtitles list the name/RGI identification number of each glacier as well as the sampling location of each surface velocity time series. The x-axes are centred on the reference date when peak velocity was reached for each event, and y-axes show the maximum-normalised surface velocity.

Notable differences between regions are further observable in the deceleration phase of individual surges. Surges in the Alaska-Yukon cluster show a median deceleration to 10% of the maximum surface velocity within 180 days of their peak. In addition, most glaciers in the sample are back to quiescent velocity 500 days after peak. In the Svalbard-Russian Arctic cluster, surging velocities are maintained at 50% of the maximum surface velocity for 180 days. Most deceleration phases last for more than 1000 days, with Negribreen and Sonklarbreen still close to 20% of peak surface velocity 2000 days after (Figure 9 15, 21). Surges in High Mountain Asia typically show a more symmetric velocity profile, centred around the date of peak velocity. We observe a median deceleration to around 30% of the maximum surface velocity at 180 days after peak, with most surges terminated by day 800. In addition, slowdown phases for surges in Svalbard-Russian Arctic show very strong seasonal fluctuations, that are mostly absent from in Alaska-Yukon and High Mountain Asia.

We interpret the regional differences in surge duration as different sensitivities of individual glaciers to the propagation of frictional instability at the glacier bed. Glaciers such as Kluane, Arnesenbreen, or RGI2000-v7.0-G-13-05693 (Figure 9 panels 5, 14, 30) show a very clear spike in surface velocity, without gradual flow acceleration before the transition point. This testifies to the rapid onset of transient feedbacks between drainage and friction (Benn and others, Reference Benn, Jones, Luckman, Fürst, Hewitt and Sommer2019b; Thø gersen and others, Reference Thøgersen, Gilbert, Bouchayer and Schuler2024). For surges with a more gradual onset such as Donjek, Osbornbreen and Gulyia ice cap (Figure 9 panels 3, 16, 27), lower enthalpy production likely first leads to the observable gradual increase in ice flow. The glacier then transitions to a velocity weakening regime (Thø gersen and others, Reference Thøgersen, Gilbert, Bouchayer and Schuler2024), as increased sliding velocity leads to lower basal friction triggering a sliding/frictional heating feedback (Benn and others, Reference Benn, Hewitt and Luckman2023), further enhancing glacier motion and propagating the surge.

We further infer longer deceleration phases in the Svalbard-Russian Arctic, as well as the drastic seasonal speed-ups, to be correlated to the large proportion of tidewater surge-type glaciers in the region (Figure 10). Indeed, contrarily to High Mountain Asia where all surge-type glaciers are land-terminating, or Alaska-Yukon where only Sit’Kusa (Turner), Valerie and La Perouse glaciers are marine-terminating, all but one of the identified surge-type glaciers in Svalbard-Russian Arctic are tidewater glaciers. The surge detected at Vallåkrabreen, the only land-terminating glacier in our Svalbard-Russian Arctic sample, displays a velocity signal similar to surges affecting land-terminating glaciers in Alaska-Yukon and High Mountain Asia. We however want to note that 3 of the tidewater glaciers in our sample, Sit’Kusa (Alaska-Yukon), Arnesenbreen (Svalbard-Russian Arctic) and Sortebræ (East Greenland), present velocity signals closer to that of land-terminating glaciers. While the surge Sortebræ still displays a strong periodicity, Arnsenbreen and Sit’Kusa only multiple surges over the studied time period. We did no find any correlation between glacier geometry as given in version 7 of the RGI (mostly length, area, altitude range and slope) and either the maximum velocity or duration of surge events. We cannot control for differences in thermal regime between land-terminating and tidewater-glaciers, or in our sample. We thus posit that ice-ocean interactions at the front of surging tidewater glaciers lead to the observed longer active phases in the Svalbard-Russian Arctic cluster. More specifically, higher hydrostatic pressure at the glacier terminus likely hampers frontal drainage of the water at the glacier bed resulting in the inability of basal drainage systems to act as efficient enthalpy sinks (e.g. Terleth and others, Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024). This ultimately leads to a more gradual termination of surges (Benn and others, Reference Benn, Jones, Luckman, Fürst, Hewitt and Sommer2019b; Reference Benn, Hewitt and Luckman2023) and higher sensitivity to seasonal meltwater input. The interpretation of our results therefore corroborates the hypothesis that surge termination is controlled by how efficiently the drainage system can evacuate the basal enthalpy (Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019a; Ravier and others, Reference Ravier, Lelandais, Vérité and Bourgeois2023; Terleth and others, Reference Terleth, Bartholomaus, Liu, Beaud, Mikesell and Enderlin2024; Thø gersen and others, Reference Thøgersen, Gilbert, Bouchayer and Schuler2024).

Figure 10. Stacks of surface velocity time series for glaciers from all clusters, by terminus type. Bold lines represent stack median and shaded areas cover the minimum to maximum range. The x-axes are centred on the reference date when peak velocity was reached for each event, and y-axes show the maximum-normalised surface velocity. N specifies the number of time series used to generate each subplot. Time series were sampled across clusters, only accounting for terminus type. Note the strong periodic component of the velocity signal for tidewater glaciers, as well as the overall greater variance in surface flow velocities.

The presented differences in surge peak surface velocity and duration distributions between Alaska-Yukon, Svalbard-Russian Arctic and High Mountain Asia reflect variations in enthalpy budgets. While individual variations might originate from different enthalpy producers, such as changes in a glacier’s thermal regime or mass balance, poor basal drainage or rapid influx of meltwater, they do not represent any fundamental contrast in surge mechanism and rather form a wide spectrum of dynamical behaviours.

4.5. Evaluation of inter-inventory consistency

We now attempt to quantify the similarity between our catalogue of surge events and that of Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) for global detection as well as Koch and others Reference Koch, Seehaus, Friedl and Braun(2023), solely for Svalbard. It is important to note that we rely only on surges that Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) classify as certain.

For this experiment, we first focus on Svalbard, for which ITS_LIVE surface velocity estimates are only reliably available from 2013, effectively reducing the mismatch in the periods considered between the inventories from Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) and Koch and others Reference Koch, Seehaus, Friedl and Braun(2023) and the one proposed here. Between the proposed inventory and that of Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023), the Jaccard index estimate for Svalbard gives an acceptable similarity between the sets of glaciers with a result of 0.61. This is due to the presence of four relatively large glaciers in Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023), which are not captured by either the surface velocity or the SAR backscatter scheme. When comparing our results to the Svalbard inventory of Koch and others Reference Koch, Seehaus, Friedl and Braun(2023), the Jaccard index quantifies a good similarity with a result of 0.71. The difference between our inventory and that of Koch and others Reference Koch, Seehaus, Friedl and Braun(2023) first lies in the detection of a surge on Austfonna Basin-3, which, while surging, currently undergoes linear deceleration with marked seasonal accelerations. Second, Koch and others Reference Koch, Seehaus, Friedl and Braun(2023) describe a 2017 surge of Orsabreen, which is invisible in the three datasets we analysed, and not reported in Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023).

In High Mountain Asia, we choose to study the regional Jaccard index for the two main clusters of surge-type glaciers in High Mountain Asia: Pamir range/Tibetan Plateau (RGI region 13) and the Karakoram Range (region 14). In the Pamirs and Tibetan Plateau, we obtain a Jaccard index of 0.58, similar to that of Svalbard. The similarity is greater in the Karakoram, where the Jaccard index is 0.81. Finally, in Alaska (region 01), the similarity index is 0.67.

Upon further investigation, false-negatives originated from our surface velocity identification scheme led to several points of discussion. First, in some rare cases, the quiescent behaviour of the given glacier is obscured as the glacier is surging over the whole considered period, preventing the computation of reliable baseline behaviour. Second, the glacier is already decelerating at the beginning of the time series, and hence, there are no positive anomalies to be detected (Austfonna Basin-3 and Nathorstbreen glacier system in Svalbard, for example). Data availability and quality problems (RGI2000-v7.0-G-13-16640, RGI2000-v7.0-G-14-08450, etc.) obscure potential surges; a point further discussed in Section 4.6.1. Missed detections from the SAR backscatter and surface elevation change scheme are a direct consequence of our changepoint detection threshold, which is purposely conservative to avoid false positives. Finally, in Alaska, the difference between our inventory and Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) is a result of their non-detection of the surges of the three lobes forming the Sít’ Tlein (Malaspina - RGI2000-v7.0-G-01-15261) glacier system. We however want to mention that Kääb and others Reference Kääb, Bazilova, Leclercq, Mannerfelt and Strozzi(2023) identify the propagation on an instability on the Seward lobe of Sít’ Tlein glacier but classify it as uncertain.

4.6. Towards a fully automated detection of glacier surges: uncertainties, challenges and perspectives