A regionally resolved inventory of High Mountain Asia surge-type glaciers, derived from a multi-factor remote sensing approach

Knowledge about the occurrence and characteristics of surge-type glaciers is crucial due to the impact of surging on glacier melt and glacier related hazards. One of the "super-clusters" of surge-type glaciers is the mountains of Asia. However, no consistent region-wide inventory of surge-type glaciers in High Mountain Asia exists. We present a regionally resolved inventory of surge-type glaciers based on their behaviour across High Mountain Asia between 2000 and 2018. We identify surge-type behaviour from surface velocity, elevation and feature change patterns using a multi-factor remote sensing approach 5 that combines yearly ITS_LIVE velocity data, DEM differences and very-high resolution imagery (Bing Maps, Google Earth). Out of the ≈ 95000 glaciers in HMA, we identified 666 that show diagnostic surge-type glacier behaviour between 2000 and 2018, which are mainly found in the Karakoram (223) and the Pamir regions (223). The total area covered by the 666 surge-type glaciers represents 19.5% of the glacierized area in Randolph Glacier Inventory (RGI) V6.0 polygons in HMA. We validate 107 previously identified glaciers as surge-type and newly identify 491 glaciers. We finally discuss the possibility of self-organized 10 criticality in glacier surges. Across all regions of HMA, the surge-affected area within glacier complexes displays a significant power law dependency with glacier length.


Introduction
Glacier surges are internally triggered, quasi-periodic oscillations of a glacier's dynamical behaviour, alternating between slow and fast flow (Meier and Post, 1969;Raymond, 1987;Sharp, 1988;Truffer et al., 2021). While periods of slow flow (or 15 quiescence) usually last for tens to hundreds of years, periods of fast flow (surge) are shorter (months to 15 years), with flow velocity reaching up to 10-1000 times its standard order of magnitude (Meier and Post, 1969;Murray et al., 2003;Pritchard et al., 2005;Mansell et al., 2012). During surges, a substantial volume of ice is transferred from a reservoir zone down-glacier to a receiving zone. This mass transfer typically leads to drastic thinning in the upper reaches of the glacier and thickening at lower elevations, often causing a substantial advance of the glacier terminus (see Sund et al. (2009) for example). 20 Glacier surges can occur on both polythermal and temperate glaciers. In the former case, the bed can oscillate between cold and warm states over the course of a surge cycle, whereas in the latter the bed remains warm throughout. These two cases have 1 been used as the basis for a two-fold classification of surge-type glaciers ("Svalbard-type" vs. "Alaskan-type") with distinct instability mechanisms ("thermal switch" and "hydrological switch", respectively) (Murray et al., 2003;Jiskoot et al., 2001).
This binary view, however, reflects neither the wide diversity of surging behaviour (see Quincey et al. (2015) for example) 25 nor the consistent associations between surging and climate, glacier geometry, and substrate characteristics. Benn et al. (2019) have argued that the full spectrum of glacier dynamic behaviour can be understood in terms of basal enthalpy balance. That is, stable steady states (non-surging behaviour) are only possible if enthalpy (thermal energy and water) generated at glacier beds from geothermal heat, friction and other sources can be evacuated from the system at the same rate. Imbalances between enthalpy inputs and outputs result in dynamic instabilities, including surges. Some combinations of regional (climatic) and local 30 (topographic, geologic) factors encourage instability, consistent with observations (Sevestre and Benn, 2015). Many aspects of enthalpy balance theory need to be worked out in detail (e.g. interactions between ice motion, friction and basal drainage), and predictions of the theory need to be tested against local and regional datasets.
Many of these studies have focused on the modern satellite era and highlighted the non-uniform distribution of surge-type glaciers in HMA, giving a more accurate representation of the prevalence of surge-type glaciers compared to ground based efforts (Kotlyakov et al., 2008;Imran and Ahmad, 2021). The only existing inventory of surge-type glaciers in HMA originates 45 from the study of Sevestre and Benn (2015). Their inventory however relies on a compilation of scientific literature, albeit considering inconsistent diagnostic criteria, identification methods and study periods .
The cyclical advance of surge-type glaciers can result in repeated and widespread glacier hazard (i.e. outburst of glacierdammed lakes, gravitational instabilities etc.) formation (Gardner and Hewitt, 1990;Ding et al., 2018;Bhambri et al., 2019;Truffer et al., 2021). Hazards associated with surge-type glaciers impact both the local environment in the path of the advancing 50 glacier terminus (Shangguan et al., 2016) and communities further downstream, in the path of meltwater produced by glaciers at the head of mountain catchments. The source of many large Glacial Lake Outburst Floods (GLOFs) in the Karakoram can be traced back to lakes which have formed behind the advancing terminus of a surge-type glacier (Hewitt and Liu, 2010;Bhambri et al., 2019). The timing and magnitude of the drainage of such lakes is difficult to predict because of the active nature of the lake's ice dam, but several GLOFs may occur from the same lake during the active phase of the damming glacier (Round 55 et al., 2017;Gao et al., 2021;Bazai et al., 2021). The recent surge and associated glacial lake formation alongside Shishpare Glacier (Bhambri et al., 2020;Muhammad et al., 2021) is just one example of the threat posed by surge-related hazards to high mountain infrastructure such as the International Karakoram Highway -a route which provides a vital link to mountain communities in southwest China/Tibet and Pakistan (Ding et al., 2018).
The impact of ice redistribution through a glacier surge's cycle on its overall rate of mass loss or gain (mass balance) over 60 multi-decadal timescales is still disputed. The flux of ice from a surge-type glacier's reservoir zone to its receiving zone results in the redistribution of a large ice volume to lower altitude, where it is more prone to ablation (Sund et al., 2009;King et al., 2021). Studies have documented increased melt rate in the expanded ablation zone, outpacing the rate of quiescent phase ice build-up in the years following surge cessation, and leading to markedly more negative mass balance (Aðalgeirsdóttir et al., 2005;Kochtitzky et al., 2019;King et al., 2021). Similarly, synchronous surging by a number of glaciers in the Ak-65 Shirak region of the Northern Tien Shan resulted in enhanced ice loss in the 1980s and 1990s when no significant change in temperature or precipitation otherwise occurred to drive enhanced ablation (Bhattacharya et al., 2021). The surge-related fluctuation of individual glacier mass balance does not appear pronounced enough to influence regional ice loss rates over shorter time scales amongst the larger surge-type glacier clusters, such as the Karakoram (Gardelle et al., 2013;Berthier and Brun, 2019). A detailed inventory of surge-type glaciers is clearly important for a region such as HMA, where glacier hazard 70 mitigation and glacier meltwater yield will be a high priority as the temperatures of the region are projected to further rise in coming decades (Kraaijenbrink et al., 2017;Bolch et al., 2019;Lalande et al., 2021).
The objectives of this paper are twofold. First, we aim to assess the occurrence of surge-type glaciers across HMA using remotely sensed data spanning the period 2000-2018. Then, we aim to clarify what are the main geometric characteristics and controls for surging in HMA, as well as providing a clearer image of the spectrum of possible surge behavior, in light of the 75 enthalpy balance theory (Benn et al., 2019).

Data and methods for surge-type glacier identification
We propose to identify surge-type glaciers from distinct widely used criteria (Dowdeswell and Williams, 1997;Copland et al., 2003;Grant et al., 2009), directly contrasting with regional trends of glacier mass loss (Brun et al., 2017;Shean et al., 2020;Bhattacharya et al., 2021), slow down (Dehecq et al., 2019) and retreat . First, substantial and spatially concen-80 trated surface elevation changes (over 1-10 years), either in the reservoir zone or at lower elevations (near the glacier terminus), are assumed to be typical of surge-type glacier ice mass redistribution. Then, substantial variations in a glacier's velocity field over a similar time period is taken as indicative of surging. Finally, we consider chaotic glacier-wide crevasse patterns mixing longitudinal and transverse crevasses at low altitudes (typically close to the front) as diagnostic of active glacier surges. While terminus advance is often a consequence of glacier surges, not all surge-type glaciers display terminal advance during the active 85 phase and it is therefore not considered as a discriminating criterion for surge identification (Paul et al., 2017;Steiner et al., 2018).

2.1 Glacier surface elevation changes
We considered two contrasting patterns of glacier surface elevation change to be diagnostic of surge-type glaciers. We aimed to identify glaciers which exhibited substantial and widespread surface elevation gain (thickening) over their receiving zone, 90 which is also commonly accompanied by surface elevation decrease (thinning) over a glacier's reservoir zone due to the flux of ice between the two areas ( Figure 1 (A)). We consider this pattern of elevation change to be diagnostic of the active phase of a glacier's surge cycle. We also aimed to identify glaciers which exhibited substantial and widespread reservoir zone thickening and concomitant terminal thinning, which we attribute to quiescent phase ice build-up occurring alongside post-surge phase ice loss at lower elevations (Figure 1 (B)). To identify the mass displacement associated with the active phase of a glacier's surge 95 cycle, or substantial ice mass build up associated with the quiescent phase of the surge cycle, we examined multi-temporal datasets of surface elevation change (dH) over HMA glaciers. We primarily used the dH data generated by Shean et al. (2020), which covers the period 2000-2018, to identify the elevation change patterns described above. To aid the identification of surge-like behaviour which may have occurred at the beginning of the period covered by Shean et al. (2020), the signal of which may have been obscured by subsequent long-term thinning, we also examined dH data generated by Hugonnet et al. 100 (2021) over the periods 2000-2004 and 2005-2009, and by Brun et al. (2017) over the period 2000-2016. Distinctive surface elevation patterns were identified manually and corroborated by multiple users.

Glacier surface velocity
The NASA MeaSUREs ITS_LIVE project (Gardner et al., 2019) provides measurements of glacier surface velocity at monthly to yearly temporal resolution over all major land ice regions, with a resolution of 240m. Glacier surface velocities are estimated 105 over the 1985-2018 period from Landsat 4,5,7 and 8 images using the auto-RIFT feature tracking algorithm (Gardner et al., 2018). As a consequence of unequal data quality and scarcity in the years covered by the earlier Landsat data archives, we only considered yearly glacier surface velocity derived from Landsat 7 and 8 imagery, between 2000 and 2018, which also matches the period covered by surface elevation change datasets.

110
The ITS_LIVE yearly glacier surface velocity data are provided with an associated error map. From the collection of yearly glacier surface velocity and their respective error maps, we then form V 0 , the error-weighted mean velocity map for the study period, as follows: where V i is the glacier surface velocity dataset for year i and the weights w i are defined as: with i being the error map for year i.
The yearly map of velocity change (dV ) is computed as the difference between V 0 and V i for each year i. We then follow the method from Mouginot and Rignot (2015) and eliminate potential outliers using a 3x3 median filter. From these data, we aim to identify positive velocity anomaly with a magnitude commensurate to that of a glacier surge. We 120 thus assume the distribution of dV to be a 0-mean Gaussian distribution for glaciers with stable flow over the studied period ( Figure 2 (a)). Surging behaviour typically displays more complex patterns and one can expect the distribution of surface velocity variations to be either positively (during the active phase) or negatively (quiescent phase) heavy-tailed ( Figure 2 (b)).
In the present work, we aim to identify glaciers displaying strong positive heavy-tails resulting from active surges on particular years. We thus compute the range between the median (P 50 ) and the percentile 95 (P 95 ) for each glacier and year : where IPR i is the inter-percentile range for year i. From Equation 3, we propose to quantify surge magnitude as a surge index, defined as follows: where s i is the surge index for year i. k corresponds to a threshold for surge identification. Surge velocities are usually described 130 as 10 to 1000 times greater than "standard" glacier velocity (Cuffey and Paterson, 2010;Jiskoot, 2011). However, surges in HMA have been documented with surface velocities close to 4 times the standard velocity (Quincey et al., 2011(Quincey et al., , 2015Bhambri et al., 2017;King et al., 2021). Thus, we here work with k = 4.
In practice, a glacier is defined as surge-type when it presents at least one occurrence of s i >= 1 within the whole study period. Equation 4 allows to quantify the relationship between the observed velocity anomaly for each specific year, and V 0 .

135
We thus highlight variations in the magnitude of the active phase of surge events from one year to another.

Glacier surface features
Substantial changes occur in the stress and strain regime at the surface of a surge-type glacier during its active phase which result in the widespread modification of surface structures such as crevasse patterns (Jennings and Hambrey, 2021). Such changes in structural glaciology have previously been used to identify glacier surging (See Grant et al. (2009);Lovell et al. 140 (2018) for example).
We performed systematic visual checks for the presence and temporal variation of diagnostic crevasse patterns on Very High Resolution (VHR, spatial resolution finer than 5m) optical imagery. We used multi-temporal, VHR imagery in Google Earth (GE) and additional images in Bing Maps to confirm or rebut the identification of surge-type glaciers, alongside glacier surface elevation and velocity changes. We aimed to identify crevasse patterns indicative of longitudinal extension or lateral shear 145 margins, which are indicative of ice flow increases and ice mass redistribution in the down-glacier direction (Figure 3), as well as potholes on the surface of the glacier which suggest ice flow stagnation. The availability of VHR optical imagery allowed for the examination of structures in glacier reservoir and receiving zones, whereas moderate resolution optical imagery (Landsat, for example) typically only allows for the identification of larger glacier features such as looped moraines, which might not be prevalent on all surge-type glaciers 150
The RGI provides a snapshot of glacier extent near the beginning of the 21st century and so current glacier extent, particularly of glaciers which have since surged, may differ somewhat from the area estimates associated with the RGI. However, such differences are unlikely to cause erroneous comparisons at a regional level and we have not further modified RGI glacier In the present study, we manually digitized the maximal surface area impacted by a given surge within glacier complexes observed between 2000 and 2018. Surge-affected areas are delineated from the extent of substantial glacier surface elevation change caused by ice mass transfer during the active phase of a surge, using very-high resolution optical imagery from Google 165 Earth and Bing Maps to corroborate each delineation. We therefore provide an estimate of the maximum surface area affected by surges (n = 300/666), in addition to the positively-biased surface area of surge-type RGI polygons. Due to the relatively low resolution (size of one pixel on the ground) of the velocity data used in this study (240m, see Section 2.2), we follow the method proposed by Dehecq et al. (2019) and only consider RGI polygons with an area greater than 5 km 2 for velocity-based analyses.

Identification of surge-type glaciers
We carried out a preliminary manual search of surface elevation change datasets to highlight possible surge-type glaciers by identifying typical surge-induced surface elevation change patterns (see Section2.1) which contrasted directly with local thinning patterns evident on non-surge-type glaciers. In parallel, automatic surge detection was carried out from the ITS_LIVE 175 datasets on RGI polygons with an area greater than 5 km 2 . Two sub-inventories were thus created before being merged together.
Each glacier was then individually inspected for the presence of surface features typical of active surging (see Section 2.3). To be qualified as surging, individual glaciers must display at least 2 of the 3 proposed identification criteria of rapid changes in surface elevation, surface velocity, surface crevassing and potholes. By combining different criteria, we typically mitigate false identification of phenomena other than surges which may cause a similar signal in, for example, surface velocity data such 180 as lake-terminating glacier frontal speed-up. The multi-factor approach also mitigates problems such as remnant noise in the surface elevation and velocity datasets among others.

Results
Out of the ≈ 95000 glaciers in HMA, 666 have displayed surge-type behavior between 2000 and 2018 ( Figure 4). Traditional surge-type glacier clusters are clearly defined, as approximately two thirds of the identified surge-type glaciers are located in 185 the Karakoram (223) or in the Pamirs (223) ( Table 1). Smaller clusters are located in the Tien Shan and the Kunlun Shan, whilst isolated and less numerous examples are found across the Tibetan plateau and it's peripheral mountain ranges (Table 1).
From the 666 identified glaciers, 68 had been previously identified as "Observed Surging" in the RGI V6.0. We further confirm 107 glaciers, previously identified as "Probable" or "Possible" surge-type glaciers. Finally, we newly identify 491 190 glaciers, previously categorized as "No Evidence" and "Not Assigned". Compared to the Sevestre and Benn (2015)  RGI polygons of surge-type glaciers cover a total surface area of 19044 km 2 , which corresponds to ≈ 19.5% of the glacierized area in HMA. Unsurprisingly, surge-type glaciers account for the greatest portion of RGI glacier area in the Karakoram

Distribution and geometry of surge and non surge-type glaciers
Previous studies have demonstrated significant differences in glacier geometry (e.g. length, slope) between surge-type and 200 non surge-type glaciers (see Sevestre and Benn (2015) among others, for example). We here further document the impact of geometry on surge behaviour by studying glacier elevation range, slope, length area and aspect, using attributes derived from the RGI V6.0. Note that, due to the wide ranges covered by each attribute, we compare common logarithm values.
Throughout HMA, surge-type glaciers systematically present greater elevation range and area ( Figure 5). As glacier length, slope and elevation range are strongly correlated attributes, such a pattern is unsurprising. We further observe generally shal-

HIMAP region
Greater HIMAP region N STGs N Glaciers Area of STGs (km 2 ) Area of all glaciers (km 2 ) STG area (%) the median surge-type glacier slope indeed lies within one standard deviation of the median for non-surge-type glaciers. Slope alone is thus not a sufficient predictor for surge-type behavior in HMA. We however observe that the frequency of surge-type glaciers drastically increases with greater lengths and shallower slopes (Figure 6), as predicted by the enthalpy balance theory (Benn et al., 2019).

210
While in the Tien Shan, Pamirs and Kunlun Shan, aspect distributions do not significantly differ, we nonetheless note narrower distributions and increased number in south quadrants for surge-type glaciers, with ≈ 8% of glaciers flowing S in the Tien Shan, ≈ 7% flowing SW and S in the Pamirs and ≈ 9% flowing SE in the Kunlun Shan (≈ 2 % for non-surge-type glaciers in all three regions) (Figure 7). We note distinct variability in the aspect distribution of surge-type glaciers in the Tibetan Mountains (Figure 7), which show no dominant orientation. Inter-regional differences in glacier type (defined in the

Analysis of glacier geometry and its impact on surge-type behaviour
We here further document the potential impacts of geometry on surging behaviour, both at HMA and regional scales. More specifically, we study the relationships between surge-affected area in glacier complexes and median surge indices with glacier median elevation, elevation range, slope and length. We observe no particular relationship between the area affected by surge in glacier complexes and the median elevation or the slope (Figure 9). Similarly no evident correlation can be described between 225 the median surge index and the same geometric attributes. The common logarithm of the surge-affected area in glaciers complex  slope distribution in the Kunlun is uniform with no clearly identifiable mode. All identified correlations display a heavy scatter and more data is needed to consolidate these relationships.

The impact of surge-type behaviour on glacier mass balance
The mass balance measurements presented here originate from the studies of Shean et al. (2020)

Temporal variability and maximum velocities in surge active phase
We observe no major differences in the length of the active phase of surges captured by the velocity data in different regions across HMA, the distributions of which are presented in Figure 11  We note no clear difference in surge active phase duration between regions (Figure 11, c). The majority of surges typically last less than 3 years, with a median duration of 2.6 ± 0.1 years. However, the distributions are systematically heavy-tailed, with all regions displaying active phases lasting 18 years.
We finally study the distributions of duration and cumulative sum of IP R (see Equation 3) for each surge over HMA. We 270 use the latter as a crude proxy of the total energy dissipated during a given surge.
Both distributions are presented in Figure 12 and display the characteristic shape of power law (or Pareto)-like distributions.
Power laws are common in geophysical systems where energy accumulated over a significant period is released rapidly (also known as avalanching systems), and arise from a wide variety of physical phenomena (criticality and self organized criticality (SOC), among others. See e.g. Turcotte (1989), Sachs et al. (2012), Åström et al. (2014), and Corral and González (2019) for 275 more). To further study the hypothesis that the velocity anomaly during active phases and the duration of surges follow power law distributions, we use the widespread method of Clauset et al. (2009), relying on maximum likelihood estimation of power law distribution parameters (Pawitan, 2001;Bauke, 2007) and Kolmogorov-Smirnov distance minimization or log-likelihood ratio goodness-of-fit tests (See Park and Kim (1992); Press et al. (2007) for example). The proposed power law model fits the velocity anomaly during active phases better than exponential models, with log-likelihood ratio (R) significantly greater than 280 zero (see Table 3). On the other hand, the exponential model fits the duration of active phases of surges better than the best-fit power law model, with a log-likelihood ratio of -6.  Table 3. Table of power law parameters: α represents the power law exponent. σ is the standard deviation associated to α. R represents the log-likelihood ratio between power law and exponential distributions.
The fits between the estimated probability densities and the data are further presented in Figure 13. We find that the theoretical cumulative distribution function (CCDF) is in good agreement with the cumulative sum of IPR data. As mentioned above, the fitted distribution does not replicate active phase duration distribution as well. This is likely to be a consequence of our 285 heavily truncated observation period: we do not highlight surges with an active phase greater than 18 years and thus truncate the power law's heavy tail. Similarly the yearly sample rate of the ITS_LIVE data used prevents identification of active phase lasting less than a year, further truncating the potential power law.

Identification of surge-type glaciers
Our study combines different, individual criteria used in previous studies to identify glacier surging, each of which may be more or less reliable in the detection of glacier surging in different scenarios. For example, the identification of glacier surging through the monitoring of glacier terminus positions is effective when glacier terminus advance is anomalous in comparison to the local pattern of glacier retreat in response to climate conditions, or where glaciers are devoid of stagnant, debris-covered 315 tongues which may inhibit surge propagation towards glacier termini (Quincey et al., 2015). Glacier surging and mass gain related terminus advance may be confused in some scenarios (Lv et al., 2020) and additional observations are required to more accurately interpret glacier behaviour. In this study we qualify a glacier as surging if, and only if, it presented at least two of the three discriminatory criteria outlined in sections 2.1 to 2.3 (rapid surface elevation and/or velocity variations and/or intense surface crevassing and/or potholes) . We are therefore confident that our interpretations are robust and that the prevalence of 320 surging we document is an improvement on approaches using fewer surging criteria.

Identification of individual surges, temporal variability and IPR
In On the other hand, the distribution of cumulative sum of IPR is positively biased, as a result of the truncation of lower values.

330
Increasing the detection threshold does not affect the shape of the distributions.

Comparison with other inventories
Inventories of surge-type glaciers have been assembled for many sub-regions of HMA and particularly detailed studies have been carried out where hazards associated with surging have been widespread and recurring, such as the Karakoram and the Pamirs. A careful comparison of the similarities and differences between our inventory and those of previous studies is difficult 335 because of the contrast in the methods used to identify glacier surging, the criteria used to define glacier surging and because the periods examined to identify glacier surging vary substantially between studies. Sevestre and Benn (2015) proposed the only existing regional inventory of surge-type glaciers in HMA. The number of surge-type glaciers documented in our inventory significantly differs from that of Sevestre and Benn (2015), especially in the Pamirs. A further examination of the Sevestre and Benn (2015) inventory reveals that, out of the 827 surge-type glaciers 340 documented, 284 correspond to individual tributaries within glacier complexes which are not individualized in the present study. From the remaining 543 glaciers, 35 documented in the RGI V5.0 (on which Sevestre and Benn (2015) is based) do not exist in the RGI V6.0. Furthermore, we found that the proposed inventory and the one from Sevestre and Benn (2015) only share 83 identified surge-type glaciers in the Pamirs. This yields a difference of 390 surge-type glaciers between the two inventories. Upon further examination of the remaining glacier population we note a median glacier area of 1.6km 2 , a sixth of 345 the median area of the surge-type glacier population described in the present inventory (9.6km 2 ). Of those 390 in the Sevestre and Benn (2015) inventory, 30% present an area small than 1.0km 2 . Close examination of these glaciers (Ujsu Glacier and the glaciers in its direct vicinity such as Aldzhaylau and Rakzou glaciers) using the surface elevation and surface velocity change data over the period 2000-2018 did not yield any evidence of surge-type behavior. We rather observed constant glacier mass loss and recession, with no clear signal of instability-related velocity anomalies. Furthermore, no mention of such a high number 350 of surge-type glaciers in the Pamirs can be found in the literature used by Sevestre and Benn (2015). Kotlyakov et al. (2008) indeed refer to Osipova et al. (1998), mentioning "630 glaciers with indications of dynamic instability, 51 of them identified as surging", they later state that 55 surge-type glaciers had been documented up to 2006 in the Pamirs. Given these unclear, conflicting and interpretive records of unstable glacier dynamics, the small sizes of glaciers in that population compared to other surge-type glaciers in the present and other studies (see next section), the lack of evidence for unstable behavior over the 2018 suggest 223 surge-type glaciers in the same region. The limited difference between these results may suggest that few surges in the Pamirs exhibit low magnitude elevation changes (e.g. Figure 1) and require additional observations (surface velocity) to aid surge-type glacier identification. The low rate of thinning and ice mass loss in the Pamirs (see Shean et al. 365 (2020) for example) has also likely preserved the strong signal of elevation change associated with surging in multi-decadal geodetic studies.
In contrast, Rankl et al. (2014) identified 101 surge-type glaciers in the Karakoram from glacier surface velocity and glacier terminus observations over the period 1972-2012, which is substantially less than the 223 identified in this study using data covering the period 2000-2018. The main methodological difference between our study and Rankl et al. (2014) is our additional 370 use of surface elevation change data and the annual resolution of the ITS_LIVE glacier surface velocity fields, which were intermittent (1992, 1993, 2003 and 2006-2013) in the study of Rankl et al. (2014). Many short duration (median 2.3 years) Karakoram glacier surges ( Figure 12) may have occurred between the dates of velocity fields presented by Rankl et al. (2014) and the longer lasting mass transfer signal of the same surges has only been picked up by our subsequent analyses of glacier surface elevation change. Also focusing on the Karakoram, Bhambri et al. (2017) compiled a list of 221 surge-type glaciers, 375 2 fewer than our estimate, but using a wide variety of observations spanning the period 1840-2015. As a broad portion of the study period of Bhambri et al. (2017) covers the pre-satellite era, during which many short-duration or low magnitude surges could have occurred undocumented, we would suggest that our inventory of surge-type glaciers provides a higher time-resolved estimate of the contemporary prevalence of surge-type glaciers in the Karakoram. Mukherjee et al. (2017)   has resulted in a more complete picture of the prevalence of surge-type glaciers here. The comparisons above emphasise the benefits of a multi-factor approach to identify surge-type glaciers, particularly in regions where glacier recession and mass loss rates are high. The availability of multi-temporal, globally resolved remotely sensed datasets should facilitate the assembly of extensive, accurate surge-type glaciers inventories in the future.

Distribution and geometry of surge and non surge-type glaciers
We have studied the relationships between glacier geometry and surging behaviour. Our results are consistent with previous knowledge that surging glaciers systematically present greater areas, length and elevation range than non-surge-type glaciers (Clarke et al., 1986;Hamilton and Dowdeswell, 1996;Jiskoot et al., 1998Jiskoot et al., , 2000Jiskoot et al., , 2003Barrand and Murray, 2006;Grant et al., 2009;Sevestre and Benn, 2015).

395
Greater lengths are correlated with increased surge-affected area in glacier complexes, with up to 45% of total glacier area impacted in the Karakoram (Hispar Glacier). With surge-type glaciers being consistently in the higher-end of the glacier area spectrum, such considerations need to be further studied as they are likely to have a significant impact on glacier mass balance (King et al., 2021) and glacier-related hazards (Kääb et al. (2021) for example) at glacier scale.
While we observe lower median surface slope for surge-type glaciers, the results are not as clear-cut as for other studies 400 such as Sevestre and Benn (2015) for example. The slight, non-statistically significant correlation between slope and surging described in our study is consistent with previous results in Yukon Territory (Clarke, 1991), Spitsbergen (Hamilton, 1992;Hamilton and Dowdeswell, 1996) and Karakoram (Bhambri et al., 2017). We however demonstrated that, in HMA, longer glaciers with gentler surface slopes are more likely to be surge-type than short steep glaciers. Length and slope thus act as independent controls on surging behaviour, a result that was first predicted in the enthalpy balance theory proposed by Benn We have shown that the aspect distributions for surge-type glaciers and non-surge-type glaciers do not substantially differ, which is consistent with previous knowledge on the studied regions (Bhambri et al., 2017;Mukherjee et al., 2017;Goerlich et al., 2020). Aspect distributions however differ in the Tibetan Mountains. We interpreted this result as the consequence of ice cap outlets outnumbering valley/mountain glaciers in this area, with the aspect of the latter being controlled by topography and 410 geology, while the former present less topographically constrained flow. From this, we argue that, in the Tibetan Mountains, aspect has no control over surging behaviour.

Impact of surges on mass balance
In Section 3.3, we presented glacier mass changes computed by the studies of Shean et al. (2020) and Hugonnet et al. (2021) with an emphasis on the potential differences in mass balance between surge-type and non-surge-type glaciers. In line with 415 the results of Gardelle et al. (2013); Berthier and Brun (2019), we do not report any significant difference in mass budget between surge-type and non surge-type glaciers over multi-decadal timescales, for all the greater HIMAP regions. We however described a marked increase in mass loss following surge cessation for two glaciers in the Karakoram. These results are consistent with the findings of Kochtitzky et al. (2019); King et al. (2021) and Bhattacharya et al. (2021). While such a mass loss is clearly identifiable in a minimal time frame after surge termination, it is likely averaged out over multi-decadal 420 timescales. The regional distribution of mass changes over multi-decadal time scales is thus unlikely to be altered by surge-type behavior, unless a significant number of glaciers surge repeatedly during the studied period (see Bhattacharya et al. (2021) for more) or in case of extreme events. Surges will however locally lead to increased variability in a glacier's meltwater runoff and hamper the water resource availability, potentially jeopardizing human livelihoods in the glacier's vicinity. The impact of surging over more localized scales thus need to be further quantified, through the thorough study of a greater number of mass 425 balance measurements, covering the entire surge cycle of a significant glacier sample.

Temporal variability and velocities in surge active phase
In this study, we present distributions of surge active phase duration derived from yearly ITS_LIVE velocity data. The moments of the presented distributions (here the percentiles) are in line with previous studies focusing on each individual region (Quincey et al., 2011;Yasuda and Furuya, 2015;Quincey et al., 2015;Bhambri et al., 2017;Lv et al., 2019;Goerlich et al., 2020;Paul, 430 2020; Zhu et al., 2021). We nonetheless demonstrated that the differences between distinct HMA regions lie in the number of surges rather than the duration of individual surges. This result contrasts with the findings of Dowdeswell et al. (1991) stating that active phases in the Pamirs typically last for 1 to 2 years, and disputes the claims of Zhu et al. (2021) than Tien Shan covers a greater range of active phase duration than other regions in HMA.
We here further dispute the claim of Goerlich et al. (2020) that active phase duration in the Pamirs is "random". The presented 435 results indeed demonstrate that the distribution of surge active phase duration (both for HMA and the Pamirs) is more likely to be heavy-tailed or power law-like.

On the physics of glacier surges
Our results suggest potential power law-like and heavy-tailed distributions for a variety of surge-type glacier parameters.
Using standard methods for power law validation, we most notably highlighted that the cumulative sum of IPR over the 440 entire duration of surge active phases (as a crude proxy of energy dissipated during the active phase) potentially follows a power law distribution. The emergence of power-law distributed quantities in avalanching geophysical systems (systems that suddenly release energy slowly accumulated over a period) arises from a wide variety of physical reasons (see Corral and González (2019) for example), with the potential role of Self-Organized Criticality (Bak, 2013) in glacier surges being first discussed by Kavanaugh (2009). Given the importance of power law-generating physical phenomena in understanding 445 geophysical dynamical systems, we emphasize the need to validate the potential power-laws described in this study with 1) complementary data (Sentinel 1 velocity data for example) providing velocity and surge active phase duration estimates over a wider range of order of magnitudes and 2) adapting the existing glacier surge models proposed by Benn et al. (2019) or Thøgersen et al. (2019) to estimate the energy dissipated during the active phase of a surge.

450
In this paper, we presented a new inventory of surge-type glaciers for High Mountain Asia. This inventory is based on a multifactor remote sensing approach, combining yearly velocity fields (ITS_LIVE), surface elevation change datasets published in prior studies and very high-resolution satellite imagery (Bing Maps and Google Earth) to identify surge-type glaciers between 2000 and 2018. Overall, we identified 666 surge-type glaciers across HMA over the studied period, confirming 107 for which surging behaviour was believed "Probable" or "Possible", and newly identifying 491. We further studied the geometry of surge-455 type glaciers compared to non-surge-type glaciers and found relationships coherent with previous studies which have focused on smaller sub-samples of surge-type glaciers. We established a relationship between surging behaviour and glacier size, the most important being that for glacier complexes, the area affected by surging is dependent on overall complex length following a power law relationship. The duration of the active phase of surge-type glaciers across HMA shows little variation in between greater HIMAP regions. Finally we defined the sum of velocity anomalies during a surge as a crude proxy for energy dissipated 460 over the duration of the active phase and studied its distribution. We showed that this distribution is heavy tailed and, using standard methods for power law validation, discussed that the sum of velocity anomalies and surge duration are potentially power law distributed. We would recommend that our inventory serves as a baseline dataset in the further study of surge-type glacier characteristics.