Large-Scale Monitoring of Glacier Surges by Integrating High-Temporal- and -Spatial-Resolution Satellite Observations: A Case Study in the Karakoram

: Glacier surges have been increasingly reported from the mountain and high-latitude cryosphere. They represent active glaciological processes that affect the evolution of natural landscapes, and they possibly lead to catastrophic consequences, such as ice collapse, which threatens the downstream communities. Identifying and monitoring surge-type glaciers has been challenging due to the irregularity of the behavior and limitations on the spatiotemporal coverage of remote-sensing observations. With a focus on the Karakoram region, with concentrated surge-type glaciers, we present a new method to efﬁciently detect glacier-surging activities by integrating the high temporal resolution of MODIS imagery and the long-term archived medium spatial resolution of Landsat imagery. This method ﬁrst detects the location and initial time of glacier surges by trend analysis (trend and breakpoint) from MODIS data, which is implemented by the Breaks for Additive Seasonal and Trend (BFAST) tool. The initial location and time information is then validated with the detailed surging features, such as the terminus-position changes from Landsat, and the thickness-change patterns from surface-elevation-change maps. Our method identiﬁed 74 surging events during 2000–2020 in the Karakoram, including three tributary-glacier surges, and seven newly detected surge-type glaciers. The surge-type glaciers tend to have longer lengths and smaller mean slopes compared with nonsurge-type glaciers. A comparison with previous studies demonstrated the method efﬁciency for detecting the surging of large-scale and mesoscale glaciers, with limitations on small and narrow glaciers due to the spatial-resolution limitation of MODIS images. For the 38 surge-type nondebris-covered glaciers, we provide details of the surging, which depict the high variability (heavy-tailed distribution) in the surging parameters in the region, and the concentration of the surge initiation during 2008–2010 and 2013–2015. The updated glacier-surging information solidiﬁes the basis for a further investigation of the surging processes at polythermal glaciers, and for an improved assessment of the glacier-mass balance and monitoring of glacier hazards.


Introduction
Glacier surge, which is an unusual phenomenon that is characterized by the abnormal acceleration of the glacier flow, can have catastrophic consequences, such as ice collapse and avalanches, inundating land and blocking river valleys, and forming ice-dammed lakes and the subsequent glacier lake outburst flood (GLOF) [1][2][3][4]. Surge-type glaciers display cyclic behavior between two distinct phases: the quiescent and active phases [5,6]. During the active phase, the glacier velocity accelerates from at least 10 to 100 times that of the quiescent phase, accompanied by the substantial downward transfer of the glacier mass, main basis for identifying surging glaciers, such as the shape and area of the glacier-tongue change [43][44][45]. An emerging method is to identify the occurrences of surging events by analyzing the ice-mass-distribution changes based on multi-temporal DEMs [2,[46][47][48]. With this method, surge-type glaciers and the area impacted by surging events can be clearly identified. However, this method usually cannot determine the time information (active period) of surging events. Recently, Guillet et al. (2022) [23] presented a comprehensive inventory of the surge-type glaciers across HMA over the last two decades (2000-2020) based on the combination of surface-velocity estimations, multitemporal DEM products, and the interpretation of the surge features (surface crevassing and potholes). Guo et al. (2022) [49] built a surge-type-glacier inventory for HMA (specifically, from 1970 to 2010) by using the glacier-surface-elevation change from multiple DEMs. However, one limitation of these studies is that the recent active-phase information about the surge-type glaciers is missing because the surface-elevation observations cannot offer accurate surging time information due to the coarse temporal resolution. These prior studies mostly used from medium-to high-spatial-resolution satellite images (e.g., from instruments including the Landsat, ASTER, Envisat, ALOS/PALSAR, TerraSAR-X), which are normally constrained in terms of the scene width, spatial coverage, and temporal resolution, and therefore are limited in detecting glacier surging events over extended regions with the necessary time information. Given that the occurrence of surging events is essentially hard to predict, the updating of the current inventory of surge-type glaciers, and the continued monitoring of surging events, remain challenging tasks for the community.
In this study, we propose a new method for deriving the location and temporal information of glacier-surging events over extended regions by analyzing the trends of the normalized difference snow index (NDSI) time series derived from the Moderate-Resolution Imaging Spectroradiometer (MODIS). With global coverage, high temporal resolution (daily), and a long-term timespan (since 2000), the MODIS imagery offers the potential for the regular monitoring of the Earth's surface changes over extensive regions. The rationality of the method lies in the considerable snow/ice mass redistribution caused by surging, which usually leads to land-cover changes, from bare land or debris to snow/ice (e.g., terminus advance, snow/ice-coverage change on glacier tongues), which will be reflected as an abrupt jump (breakpoint) or significant increase from the NDSI time series. We employed the Breaks for Additive Seasonal and Trend (BFAST) tool, which has been proposed [50] and widely used for time series analyses [51][52][53], to detect such breakpoints and trends in the NDVI time series. The first objective of this study was to develop an efficient approach for detecting the location and time information of glaciersurging events over extended regions. This information was further validated and used to guide the detailed retrieval of the surging features (terminus change/active period) by referring to high-spatial-resolution Landsat images and surface-elevation-change maps. The second objective of this study was to update the current database of surge-type glaciers and the surging characteristics (e.g., active period, the magnitude of terminus advance) in the 21st century in the Karakoram, with the aim of improving the references for studying and understanding the glacier-surging events in the mountainous cryosphere.

Study Area
The Karakoram is located in HMA (Figure 1), which hosts the largest glacier concentration outside the polar regions [54,55], with elevations varying between 3793 m and 8569 m. There are 10,937 glaciers within the region, covering a total area of 20,577.98 km 2 , which accounts for 21% of the total glacierized area in HMA, according to the Randolph Glacier Inventory 6.0 (RGI v6.0) [56]. With a cold semiarid climate and rugged terrain, this region is a hotspot of large-scale continental mountain glaciers and debris-covered glaciers [19,42,57]. In situ climatological data are rarely available in this region with a high altitude and harsh environments. Previous literature has documented that the climate in the Karakoram region is influenced by westerly winds in the winter and the Asian monsoon in the summer. Snowfall in the upper glacier accumulation zones (between 4800 and 6000 m above sea Remote Sens. 2022, 14, 4668 4 of 19 level) accounts for a major part of the precipitation in the region [9,37,58]. Glacier-surging events and surge-type glaciers have been widely found in the study region, according to the glacier surface features (e.g., heavily crevassed surface, collapsing or downwasting at the tongue, and looped moraines), surface-flow velocities, and surface-elevation changes [9,19,23,38,42]. In contrast with the significant glacier-mass loss in other parts of HMA in the global climate response [59][60][61], equilibrium or slight glacier-mass gain during the early 21st century was observed in the Karakoram [54]. A recent study suggests that the glacier melting in this region has been accelerating in recent years [26].
is a hotspot of large-scale continental mountain glaciers and debris-covered gla [19,42,57]. In situ climatological data are rarely available in this region with a high alti and harsh environments. Previous literature has documented that the climate in the rakoram region is influenced by westerly winds in the winter and the Asian monsoo the summer. Snowfall in the upper glacier accumulation zones (between 4800 and 60 above sea level) accounts for a major part of the precipitation in the region [9,37,58]. G ier-surging events and surge-type glaciers have been widely found in the study reg according to the glacier surface features (e.g., heavily crevassed surface, collapsin downwasting at the tongue, and looped moraines), surface-flow velocities, and sur elevation changes [9,19,23,38,42]. In contrast with the significant glacier-mass loss in o parts of HMA in the global climate response [59][60][61], equilibrium or slight glacier-m gain during the early 21st century was observed in the Karakoram [54]. A recent s suggests that the glacier melting in this region has been accelerating in recent years [

Materials
Our proposed method mainly relies on the freely accessible MODIS pro (https://modis.ornl.gov/globalsubset/, accessed on 18 March 2021), and on Landsat im archived on the platform of Google Earth Engine (GEE). The MODIS instrument opera on both the Terra and Aqua spacecrafts provides high-temporal-resolution observa that are beneficial for detecting changes from the regional to global scales. In this st we used the MODIS surface-reflectance product (MOD09A1), at a spatial resolution o m, and a temporal resolution of 8 days, spanning from February 2000 to December 2 This product represents a composition of the best observations within eight days, therefore contains less noise associated with atmospheric conditions and other distur factors compared with daily products (e.g., surface-reflectance daily global MOD0 product). The TM, ETM+, and OLI sensors aboard the Landsat mission provide longobservations of the Earth's surface at a moderate spatial resolution of 30 m, which i vorable for the confirmation and extraction of detailed surging features, such as gla

Materials
Our proposed method mainly relies on the freely accessible MODIS product (https: //modis.ornl.gov/globalsubset/, accessed on 18 March 2021), and on Landsat images archived on the platform of Google Earth Engine (GEE). The MODIS instrument operating on both the Terra and Aqua spacecrafts provides high-temporal-resolution observations that are beneficial for detecting changes from the regional to global scales. In this study, we used the MODIS surface-reflectance product (MOD09A1), at a spatial resolution of 500 m, and a temporal resolution of 8 days, spanning from February 2000 to December 2020. This product represents a composition of the best observations within eight days, and therefore contains less noise associated with atmospheric conditions and other disturbing factors compared with daily products (e.g., surface-reflectance daily global MOD09GQ product). The TM, ETM+, and OLI sensors aboard the Landsat mission provide long-term observations of the Earth's surface at a moderate spatial resolution of 30 m, which is favorable for the confirmation and extraction of detailed surging features, such as glacierterminus changes and the surging period (starting and ending of the terminus advance). Because of the Scan Line Corrector Failure of the Landsat 7 instruments, Landsat 7 images that are free from data gaps on the studied glaciers were used in this study. The Landsat images used were carefully selected (mostly acquired in June-September) to limit the influence of seasonal snow on the identification of the glacier extents.
The auxiliary datasets used in this study include the RGI v6.0 glacier inventory dataset [56], and the surface-elevation-change maps (hereafter, dh/dt maps) available at http://maps.theia-land.fr/theia-cartographic-layers.html (accessed on 8 January 2021) [26]. The glacier inventory dataset consists of various glacier metrics (area, length, mean elevation, mean slope, and aspect) that are used to identify the characteristics of surge-type glaciers. Hugonnet et al. (2021) [26] generated dh/dt maps for different periods (e.g., 2000-2005, 2005-2010, 2010-2015, 2015-2020, 2000-2020) on the basis of annual dh/dt maps derived by fitting the time series of the elevation measurements constructed from the stereo ASTER images. We used the dh/dt maps [26] to validate the pattern of the glaciersurface-elevation changes associated with glacier surging (e.g., substantial thickening on the receiving zone (lower reaches) of glaciers versus thinning on the reservoir parts (higher reaches)). This pattern is particularly important for identifying debris-covered-glacier surging, which generally shows no obvious terminus advance, or the terminus advance cannot be clearly revealed from the Landsat images due to the debris cover.

Methods
We propose the combination of the high temporal resolution of MODIS and the moderate spatial resolution of Landsat images to detect glacier surging and retrieve the associated characteristics (space, initiation date, duration period, terminus-advance length). The method ( Figure 2) consists of three main sections: (1) the preprocessing of MODIS images to obtain pixel-by-pixel NDSI time series of the study glaciers; (2) the identification of glacier surging by detecting the trend changes of each NDSI time series with the BFAST module; (3) the validation and fine-featured determination of the glacier surging by referring to Landsat images and surface-elevation-change maps (dh/dt maps) [26].
terminus changes and the surging period (starting and ending of the terminus advance). Because of the Scan Line Corrector Failure of the Landsat 7 instruments, Landsat 7 images that are free from data gaps on the studied glaciers were used in this study. The Landsat images used were carefully selected (mostly acquired in June-September) to limit the influence of seasonal snow on the identification of the glacier extents.
The auxiliary datasets used in this study include the RGI v6.0 glacier inventory dataset [56], and the surface-elevation-change maps (hereafter, dh/dt maps) available at http://maps.theia-land.fr/theia-cartographic-layers.html (accessed on 8 January 2021) [26]. The glacier inventory dataset consists of various glacier metrics (area, length, mean elevation, mean slope, and aspect) that are used to identify the characteristics of surge-type glaciers. Hugonnet et al. (2021) [26] generated dh/dt maps for different periods (e.g., 2000-2005, 2005-2010, 2010-2015, 2015-2020, 2000-2020) on the basis of annual dh/dt maps derived by fitting the time series of the elevation measurements constructed from the stereo ASTER images. We used the dh/dt maps [26] to validate the pattern of the glacier-surfaceelevation changes associated with glacier surging (e.g., substantial thickening on the receiving zone (lower reaches) of glaciers versus thinning on the reservoir parts (higher reaches)). This pattern is particularly important for identifying debris-covered-glacier surging, which generally shows no obvious terminus advance, or the terminus advance cannot be clearly revealed from the Landsat images due to the debris cover.

Methods
We propose the combination of the high temporal resolution of MODIS and the moderate spatial resolution of Landsat images to detect glacier surging and retrieve the associated characteristics (space, initiation date, duration period, terminus-advance length). The method ( Figure 2) consists of three main sections: (1) the preprocessing of MODIS images to obtain pixel-by-pixel NDSI time series of the study glaciers; (2) the identification of glacier surging by detecting the trend changes of each NDSI time series with the BFAST module; (3) the validation and fine-featured determination of the glacier surging by referring to Landsat images and surface-elevation-change maps (dh/dt maps) [26].

Construction of NDSI Time Series
The NDSI metric has been widely used to extract the snow/ice cover and glacier extent, and to monitor glacier changes [62][63][64]. The NDSI (Equation (1)) is calculated based on the Remote Sens. 2022, 14, 4668 6 of 19 differences in the green-band reflectance (ρ GREEN ) (band 4) and shortwave-infrared-band reflectance (ρ SWIR ) (band 6): The NDSI time series may contain extreme values (beyond the normal range of [−1, 1]) and data gaps resulting from the influence of extreme weather and instrument malfunction, which need to be processed by filtering and filling. We filtered out all the NDSI values beyond the normal range, and we applied NDSI image filtering using quantile statistics in a 3 × 3 sliding window. After the noise filtering, the percentage of data gaps for each image was mostly less than 5%, and the missing observations accounted for less than 5% of each time series sample. From the viewpoint of NDSI time series, the missing observations are randomly and discretely distributed, and therefore they can be estimated by the linear fitting of the two adjacent valid NDSI values. To identify the NDSI changes related to glacier surging, we restricted the NDSI time series analysis to pixels within a 5 km buffer of the glacier polygons for efficient and simplified processing.

Identification of Surge-Type Glaciers
Glacier surging, accompanied by substantial snow/ice transportation from higher reaches to lower reaches, and sometimes the advance of the glacier terminus, will present an abrupt or persistent increase in the NDSI after the events of surging. This abrupt or persistent increase in the NDSI usually occurs on the glacier terminus or on the middle parts of debris-covered glacier tongues, which shift to endurable snow/ice cover from bare lands/debris cover due to the substantial mass from the upper snow/ice reservoirs. We characterized this abrupt or persistent increase in the NDSI by analyzing the long-term time series of the NDSI by using the BFAST package implemented in 'R' [50]. The BFAST algorithm can effectively decompose time series data into three components, including trend, seasonal, and remaining noise, enabling the robust detection of the abrupt and gradual changes in the trend component [50]. This tool has been successfully applied to detect abrupt changes in the time series of variables caused by land-cover changes [51,53,65,66].
In the decomposition of NDSI time series, we set the two general parameters 'h' and 'break' as 0.15 and 1, respectively. The 'h' parameter refers to the minimum period length (related to the whole period length) for the identified different states. By setting it as 0.15, we restricted the minimum period length to three years (total study period of 20 years) to exclude short-term variations or the influence of noise in the time series. The 'break' parameter represents the maximum number of breakpoints detected from the trend component. We only extracted one breakpoint for each NDSI time series, which is the most significant breakpoint, corresponding to the most prominent abrupt change. Other parameters, such as the maximum iteration (2), seasonal model (harmonic), and level of significance (95%), were set to default values.
From the breakpoint-detection results, we identified two types of NDSI-trend changes pertinent to glacier surging. First, the abrupt increase in the NDSI-trend component indicates the possible land-cover shift from no snow/ice to snow/ice. When an abrupt change is detected, the trend component is split into two stages: the pre-change phase and post-change phase, which may separately display upward or downward trends, or no significant trends ( Figure 3). The date of the breakpoint corresponds to when an abrupt change (e.g., from bare land/debris-covered with a low NDSI to full snow/ice coverage with a high NDSI) roughly occurred. We applied the following criteria with three conditions to identify the grids associated with glacier surging. where NDSI post and NDSI pre refer to the NDSI values at the start point of the post-change phase and the endpoint of the pre-change phase, respectively. Their difference indicates the magnitude of the NDSI jump before and after the breakpoint, and a positive magnitude denotes an NDSI increase. T post denotes the trend of the NDSI time series in the postchange phase. NDSI pre−mean represents the mean NDSI value in the pre-change phase. The threshold of the first condition is to identify the pixels with a substantial increase in the mean NDSI level. Because glaciers normally keep stagnant or slowly recede due to melting in the post-surge period, we applied a threshold on the NDSI trend in the post-stage T post (−0.0006) to filter out cases with fast decreases in the snow/ice cover, which may be associated with extreme events, such as widespread and heavy snowfall. The thresholds used in the former two conditions were empirically determined by several rounds of testing. The third condition is to filter out cases in which snow/ice coverage dominated even before the abrupt change (glacier surging). The threshold of 0.4 has been widely used to the classify the snow/ice cover of MODIS observations [67].
where NDSI and NDSI refer to the NDSI values at the start point of the post-change phase and the endpoint of the pre-change phase, respectively. Their difference indicates the magnitude of the NDSI jump before and after the breakpoint, and a positive magnitude denotes an NDSI increase. T denotes the trend of the NDSI time series in the post-change phase. NDSI represents the mean NDSI value in the pre-change phase. The threshold of the first condition is to identify the pixels with a substantial increase in the mean NDSI level. Because glaciers normally keep stagnant or slowly recede due to melting in the post-surge period, we applied a threshold on the NDSI trend in the post-stage T (−0.0006) to filter out cases with fast decreases in the snow/ice cover, which may be associated with extreme events, such as widespread and heavy snowfall. The thresholds used in the former two conditions were empirically determined by several rounds of testing. The third condition is to filter out cases in which snow/ice coverage dominated even before the abrupt change (glacier surging). The threshold of 0.4 has been In the second type, we found that, in the case of non-abrupt changes, the trend component of the NDSI time series, with a trend > 0.0001, is likely associated with glacier surging (e.g., the long-term persistent or ongoing advance of the glacier terminus). The magnitude of the NDSI trend (>0.0001) is abnormally high when analyzing the distribution of the trend magnitude of all the grids, given that glaciers generally retreat or vary between retreat and advance in response to climate change. Therefore, we classified all the pixels showing a trend > 0.0001 in the NDSI-trend component as probably related to glacier surging.
The above two diagnostic criteria retrieved a set of individuals or clusters of pixels/grids of interests. We applied a spatial filter process to remove the individual pixels and retrieve grid clusters within or around the glacier outlines by overlaying the 2 km buffer zone of the glacier outlines. The grid clusters were linked with the potential surge-type glaciers, with the locations (grids) and time information (the date of the breakpoint for the first type, and the starting year 2000 for the second type) serving as the initial guidance on where and when glacier surging may have occurred. The high-resolution Landsat images were then used to validate and determine the exact terminus-position changes and the active phases of surging.
The effectiveness of NDSI time series analyses for revealing glacier-surging activities is exemplified in Figure 3. Pixel 1 (one of the surging pixel clusters) shows an obvious abrupt rise in the trend component, with the NDSI increasing abruptly from 0.3 to 0.7 at the breakpoint (in 2008). The glacier G1 (GLIMSID: G077970E34577N), which covers the location of Pixel 1, was confirmed as a surge-type glacier by analyzing the time series of the Landsat images. Similarly, Pixel 2 (one of the surging pixel clusters) on glacier G3 (GLIMSID: G078078E34602N), where the NDSI time series show a persistent increase from 2000 to 2020, was identified as a surge-type glacier that persistently advanced during the study period. By contrast, the NDSI trends on the nonsurging impacted areas are different, as exemplified by Pixel 3 (one of the MODIS pixels) and Pixel 4 (one of the MODIS pixels). Pixel 3, located on the boundary of glacier G7 (GLIMSID: G077736E34888N), shows a significant decline in the NDSI-trend component due to the glacier recession. No significant trend or abrupt change was identified for Pixel 4, which is located at the inner-glacier boundaries. Additionally, on glaciers G2 (GLIMSID: G078017E34593N), G4 (GLIMSID: G077894E34764N), G5 (GLIMSID: G077682E34933N), and G6 (GLIMSID: G077715E34929N), there are clusters of grids showing NDSI trends similar either to Pixel 1 or Pixel 2.

Validation of Glacier Surging and Retrieval of Surging Details
The MODIS-based NDSI analysis identified the grid clusters that are potentially impacted by glacier-surging activities (termed as surging-impacted areas). Two methods were employed to verify whether surging occurred on the glaciers, and to retrieve the associated details (terminus changes/active periods) of the surging. For the surgingimpacted grid clusters located around the clean-ice or partly debris-covered terminus, the glacier terminus positions from 2000 to 2020 were retrieved from the Landsat imagery. We used the Google Earth Engine Digitization Tool (GEEDiT) [12,68] to efficiently process the Landsat images for retrieving the glacier outlines. The GEEDiT tool is particularly efficient at selecting remotely sensed images according to the date, cloud volume, and other parameters with a standard browser, and it is able to manually extract the centerline and terminus positions of glaciers. The terminus-position change was measured as the distance between each terminus position and the earliest observed glacier terminus position along the downward centerline, according to the standard single-centerline approach [12,69]. Then, for each potential surging event, the time series of the terminus-position changes were derived. We excluded the label of 'surge-type' for these pixels and glaciers when no signs of surging characteristics were found, whereas for those exhibiting surging characteristics, we obtained fine details, such as the exact initiation year and duration period of the terminus advance, total distance, and rate of terminus advance, from the time series of the terminusposition changes. Note that our classification of the active-phase (duration) period is solely based on the changes in the terminus position rather than the velocity.
For grid clusters located within the glacier polygon or at the debris-covered glacier terminus, in which cases the terminus change cannot be clearly revealed from the Landsat images (due to debris cover), we utilized the dh/dt maps to examine whether there was a clear pattern of mass transfer from a reservoir area (higher reaches) to a receiving area (lower reaches). Such a substantial mass transfer, as a clear signal of surging, can be revealed from the spatially concentrated glacier-surface-elevation increase/decrease on the lower/higher reaches on the dh/dt maps. We chose dh/dt maps for which the timespan covers the starting time of the glacier surging to verify the mass-transfer pattern.

The Temporal and Spatial Distributions of Surge-Type Glaciers
Our method detected 74 places of surging activities associated with 71 surge-type glaciers in the Karakoram from 2000 to 2020 ( Figure 4). For some large glaciers, more than one occurrence of glacier surging was observed in different tributaries, termed as tributary glacier surges [70]. We identified tributary glacier surges on three glaciers (GLIMSId: G076483E35901N, G075281E36055N, and G075593E36006), which experienced surging at both the main termini and one of their tributaries. For example, the Hispar Glacier (GLIM-SID: G075281E36055N) surged at its tributary, Kunyang, in 2008, followed by surging at the tongue terminus in 2015. The identified 71 surge-type glaciers are mostly concentrated in the central and southeastern parts of the study area (Figure 4), with the total area accounting for 34% of the glacier area in the region. Among the 74 surging events, 38 occurred on the clean-ice-covered or partly debris-covered glacier tongues, which showed significant terminus advance. We term this type as clean-ice glacier surging to differentiate it from the surging at heavily or fully debris-covered glaciers, and the former type of glacier surging can be directly retrieved from optical images. We detected 33 surge-type glaciers with heavy debris cover, mainly located in the central Karakoram. The surging at these heavily debris-covered glaciers could be successfully detected as parts of the upper glacier tongues shifted from debris-covered to snow/ice due to the substantial snow/ice transfer from the upper-reservoir areas. Statistics of the years of surge initiation (a rough estimation from the MODIS-based NDSI analysis) reveal the concentration of the surging initiation during 2008-2010 and 2013-2015 (the upper-right inset in Figure 4). revealed from the spatially concentrated glacier-surface-elevation increase/decrease on the lower/higher reaches on the dh/dt maps. We chose dh/dt maps for which the timespan covers the starting time of the glacier surging to verify the mass-transfer pattern.

The Temporal and Spatial Distributions of Surge-Type Glaciers
Our method detected 74 places of surging activities associated with 71 surge-type glaciers in the Karakoram from 2000 to 2020 (Figure 4). For some large glaciers, more than one occurrence of glacier surging was observed in different tributaries, termed as tributary glacier surges [70]. We identified tributary glacier surges on three glaciers (GLIMSId :  G076483E35901N, G075281E36055N, and G075593E36006), which experienced surging at both the main termini and one of their tributaries. For example, the Hispar Glacier (GLIM-SID: G075281E36055N) surged at its tributary, Kunyang, in 2008, followed by surging at the tongue terminus in 2015. The identified 71 surge-type glaciers are mostly concentrated in the central and southeastern parts of the study area (Figure 4), with the total area accounting for 34% of the glacier area in the region. Among the 74 surging events, 38 occurred on the clean-ice-covered or partly debris-covered glacier tongues, which showed significant terminus advance. We term this type as clean-ice glacier surging to differentiate it from the surging at heavily or fully debris-covered glaciers, and the former type of glacier surging can be directly retrieved from optical images. We detected 33 surge-type glaciers with heavy debris cover, mainly located in the central Karakoram. The surging at these heavily debris-covered glaciers could be successfully detected as parts of the upper glacier tongues shifted from debris-covered to snow/ice due to the substantial snow/ice transfer from the upper-reservoir areas. Statistics of the years of surge initiation (a rough estimation from the MODIS-based NDSI analysis) reveal the concentration of the surging initiation during 2008-2010 and 2013-2015 (the upper-right inset in Figure 4).  Among the detected 74 surging activities, the majority (63) were identified by the breakpoint detection (breakpoint detected in the NDSI-trend component), with 33 and 30 cases on the clean-ice and debris-covered glacier tongues, respectively. The 11 remaining surges are related to the significantly increasing NDSI trends (hollow circles in Figure 4), with 5/6 cases of clean-ice/debris-covered glacier surging. Two styles of clean-ice surging are exemplified at glaciers G8 (GLIMSID: G076225E36119N) and G9 (GLIMSID: G075542E36235N) ( Figure 5). The first style is characterized by abrupt glacier advance in a short period. In this case, a cluster of pixels showing an evident breakpoint in the NDSI trend was identified on the terminus position of G8 (Figure 5a), providing a rough time for the surging occurrence (around 2009). The terminus-position changes reveal that fast terminus advance persisted during the active period of 2009-2014 (Figure 5c,d), in contrast with the high stability of the glacier-terminus position out of the active periods (before 2009 and after 2014). The second style is featured by persistent slow glacier advance during a long period (on the scale of decades), as exemplified by the glacier G9 (Figure 5g,h, persistent glacier advance since the early 2000s). This type of surging usually shows a stable snow/ice-coverage increase around the glacier terminus, with the clustering pixels presenting the significant increasing trends of the NDSI (Figure 5e,f). Among the detected 74 surging activities, the majority (63) were identified by the breakpoint detection (breakpoint detected in the NDSI-trend component), with 33 and 30 cases on the clean-ice and debris-covered glacier tongues, respectively. The 11 remaining surges are related to the significantly increasing NDSI trends (hollow circles in Figure 4), with 5/6 cases of clean-ice/debris-covered glacier surging. Two styles of clean-ice surging are exemplified at glaciers G8 (GLIMSID: G076225E36119N) and G9 (GLIMSID: G075542E36235N) ( Figure 5). The first style is characterized by abrupt glacier advance in a short period. In this case, a cluster of pixels showing an evident breakpoint in the NDSI trend was identified on the terminus position of G8 (Figure 5a), providing a rough time for the surging occurrence (around 2009). The terminus-position changes reveal that fast terminus advance persisted during the active period of 2009-2014 (Figure 5c,d), in contrast with the high stability of the glacier-terminus position out of the active periods (before 2009 and after 2014). The second style is featured by persistent slow glacier advance during a long period (on the scale of decades), as exemplified by the glacier G9 ( Figure  5g,h, persistent glacier advance since the early 2000s). This type of surging usually shows a stable snow/ice-coverage increase around the glacier terminus, with the clustering pixels presenting the significant increasing trends of the NDSI (Figure 5e,f). The surges on the fully debris-covered glacier tongues are exemplified in Figure 6. These surging activities were confirmed with the dh/dt maps, which show abnormal thickening on the lower reaches (receiving area), in contrast with significant lowering on the upper reaches (reservoir area) (e.g., Figure 6c,g). For the fully debris-covered glacier tongues, glacier surging may lead to terminus advance, or cause snow/ice-coverage changes in the middle parts of the glacier tributaries (around the snow altitudes covered by snow/ice), as exemplified by the glaciers G10 (GLIMSID: G076527E35548N) and G11 (GLIMSID: G075780E36040N), respectively. For G10, we detected pixels showing The surges on the fully debris-covered glacier tongues are exemplified in Figure 6. These surging activities were confirmed with the dh/dt maps, which show abnormal thickening on the lower reaches (receiving area), in contrast with significant lowering on the upper reaches (reservoir area) (e.g., Figure 6c,g). For the fully debris-covered glacier tongues, glacier surging may lead to terminus advance, or cause snow/ice-coverage changes in the middle parts of the glacier tributaries (around the snow altitudes covered by snow/ice), as exemplified by the glaciers G10 (GLIMSID: G076527E35548N) and G11 (GLIMSID: G075780E36040N), respectively. For G10, we detected pixels showing significant increasing trends of the NDSI around the glacier terminus (Figure 6a,b). The surfaceelevation change (dh/dt 2000-2020) reveals remarkable thickening on the glacier tongues (Figure 6c). By the visual interpretation of the Landsat images, G10 has significantly advanced (Figure 6d), although the terminus position cannot be retrieved directly from the spectral index of the optical images due to debris cover. Another glacier, G11, was typical of glaciers with no terminus advance during surging. On the upper glacier tongues of this glacier (Figure 6e), there are grids displaying evident breakpoints in the NDSI-trend component (Figure 6f). The breakpoints of these grids roughly estimated the surging initiation (2013) and active period (2013-2015). The dh/dt maps (2010-2020) revealed substantial thickening on the receiving zone (middle reaches of the glacier tongue) and thinning on the reservoir parts (higher reaches) (Figure 6g). It seems that this glacier surging did not impact the lowermost 10 km of debris-covered glacier tongues (Figure 6g,h), and therefore, there was no significant terminus advance.
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 19 significant increasing trends of the NDSI around the glacier terminus (Figure 6a,b). The surface-elevation change (dh/dt 2000-2020) reveals remarkable thickening on the glacier tongues (Figure 6c). By the visual interpretation of the Landsat images, G10 has significantly advanced (Figure 6d), although the terminus position cannot be retrieved directly from the spectral index of the optical images due to debris cover. Another glacier, G11, was typical of glaciers with no terminus advance during surging. On the upper glacier tongues of this glacier (Figure 6e), there are grids displaying evident breakpoints in the NDSI-trend component (Figure 6f). The breakpoints of these grids roughly estimated the surging initiation (2013) and active period (2013-2015). The dh/dt maps (2010-2020) revealed substantial thickening on the receiving zone (middle reaches of the glacier tongue) and thinning on the reservoir parts (higher reaches) (Figure 6g). It seems that this glacier surging did not impact the lowermost 10 km of debris-covered glacier tongues ( Figure  6g,h), and therefore, there was no significant terminus advance. For the 38 surging events on clean-ice or partly debris-covered glaciers, we further determined the details, including the active periods/surge durations (defined as the start year when the terminus changes accelerate, and the end year when the terminus changes slow down), the total terminus-advancing distance, and the terminus-advancing rate (a proxy of the ice velocities) according to the time series of the terminus-position changes (e.g., Figure 5d,h). Here, to better differentiate the different styles of surging rather than the underlying mechanisms, we classified the surge-type glaciers into the traditional Alaskan type and Svalbard type, according to the features of the terminus changes over time (Figure 7). The start and end years of the active period are flagged with black squares for each glacier (Figure 7). Note that for some of the ongoing surging glaciers (particularly the Svalbard type), the start year or end year of the active period is not discernible from the time series of the terminus changes, and they are flagged with black dashed squares. The majority of the glaciers surging on clean-ice (or partly debris-covered) glacier For the 38 surging events on clean-ice or partly debris-covered glaciers, we further determined the details, including the active periods/surge durations (defined as the start year when the terminus changes accelerate, and the end year when the terminus changes slow down), the total terminus-advancing distance, and the terminus-advancing rate (a proxy of the ice velocities) according to the time series of the terminus-position changes (e.g., Figure 5d,h). Here, to better differentiate the different styles of surging rather than the underlying mechanisms, we classified the surge-type glaciers into the traditional Alaskan type and Svalbard type, according to the features of the terminus changes over time (Figure 7). The start and end years of the active period are flagged with black squares for each glacier (Figure 7). Note that for some of the ongoing surging glaciers (particularly the Svalbard type), the start year or end year of the active period is not discernible from the time series of the terminus changes, and they are flagged with black dashed squares. The majority of the glaciers surging on clean-ice (or partly debris-covered) glacier tongues/termini (23 out of 38) can be categorized as the Alaskan type, and they experienced abrupt terminus advance during a short period (e.g., 1.8 km during 2009-2014, Figure 5c,d). We detected 15 examples of Svalbard-type surging, which exhibits gradual initiation or termination over a protracted period (over 10 years in our studied database).
The Alaskan-type surging glaciers exhibited a mean active period of 5 years, with a mean advancing rate of 374.68 m/year (varying from 100 m/year to over 1000 m/year for different glaciers). In contrast, the Svalbard-type surging generally showed a much more prolonged advancing period (15 years on average) and a lower advancing rate (generally less than 100 m/year, 65.27 m/year on average) (Figure 8a). According to the enthalpybalance theory, both the Alaskan-type and Svalbard-type surging events occur by the same mechanism. Our statistics on the relationship between the mean advancing rates and activation periods indicate that the two types generally follow the same power-law distribution, but at different spaces (Figure 8a), whereas the total advancing distances present no significant differences between the two types ( Figure 8b).
Remote Sens. 2022, 14, x FOR PEER REVIEW 12 tongues/termini (23 out of 38) can be categorized as the Alaskan type, and they exp enced abrupt terminus advance during a short period (e.g., 1.8 km during 2009-2014, ure 5c,d). We detected 15 examples of Svalbard-type surging, which exhibits gradual tiation or termination over a protracted period (over 10 years in our studied datab The Alaskan-type surging glaciers exhibited a mean active period of 5 years, with a m advancing rate of 374.68 m/year (varying from 100 m/year to over 1000 m/year for di ent glaciers). In contrast, the Svalbard-type surging generally showed a much more longed advancing period (15 years on average) and a lower advancing rate (generally than 100 m/year, 65.27 m/year on average) (Figure 8a). According to the enthalpy-bal theory, both the Alaskan-type and Svalbard-type surging events occur by the same m anism. Our statistics on the relationship between the mean advancing rates and activa periods indicate that the two types generally follow the same power-law distribution at different spaces (Figure 8a), whereas the total advancing distances present no sig cant differences between the two types ( Figure 8b).   tongues/termini (23 out of 38) can be categorized as the Alaskan type, and they experienced abrupt terminus advance during a short period (e.g., 1.8 km during 2009-2014, Figure 5c,d). We detected 15 examples of Svalbard-type surging, which exhibits gradual initiation or termination over a protracted period (over 10 years in our studied database). The Alaskan-type surging glaciers exhibited a mean active period of 5 years, with a mean advancing rate of 374.68 m/year (varying from 100 m/year to over 1000 m/year for different glaciers). In contrast, the Svalbard-type surging generally showed a much more prolonged advancing period (15 years on average) and a lower advancing rate (generally less than 100 m/year, 65.27 m/year on average) (Figure 8a). According to the enthalpy-balance theory, both the Alaskan-type and Svalbard-type surging events occur by the same mechanism. Our statistics on the relationship between the mean advancing rates and activation periods indicate that the two types generally follow the same power-law distribution, but at different spaces (Figure 8a), whereas the total advancing distances present no significant differences between the two types ( Figure 8b).   The surge parameters of the 38 surge-type glaciers (clean-ice or partly debris-covered) were analyzed to better depict the surging characteristics in the Karakoram. The surge duration, total terminus-advancing distances, and mean advancing rates follow the shape of heavy-tailed distributions (Figures 9 and 10). The surge duration varied from 2 to 20 years, with a median surge duration of 6.5 years (Figure 9a). The total terminus-advancing distances ranged from 152.12 m to 3566.98 m, with a median value of 858 m (Figure 9b). Accordingly, the mean terminus-advancing rates are mostly less than 300 m/year (with a median value of 136.76 m/year). We observed generally prominent glacier terminus advance (>300 m/year), mostly concentrated in the southeastern part ( Figure 10). The maximum terminus-advancing rate and the longest advancing distance occurred at glacier G077544E35215N (location flagged in Figure 10), which had advanced 3566.98 m at 1188.99 m/year in about three years.
The surge parameters of the 38 surge-type glaciers (clean-ice or partly debris ered) were analyzed to better depict the surging characteristics in the Karakoram. surge duration, total terminus-advancing distances, and mean advancing rates follow shape of heavy-tailed distributions (Figures 9 and 10). The surge duration varied fr to 20 years, with a median surge duration of 6.5 years (Figure 9a). The total term advancing distances ranged from 152.12 m to 3566.98 m, with a median value of 8 ( Figure 9b). Accordingly, the mean terminus-advancing rates are mostly less than m/year (with a median value of 136.76 m/year). We observed generally prominent gl terminus advance (>300 m/year), mostly concentrated in the southeastern part (Figur The maximum terminus-advancing rate and the longest advancing distance occurr glacier G077544E35215N (location flagged in Figure 10), which had advanced 3566. at 1188.99 m/year in about three years.   The surge parameters of the 38 surge-type glaciers (clean-ice or partly debris-covered) were analyzed to better depict the surging characteristics in the Karakoram. The surge duration, total terminus-advancing distances, and mean advancing rates follow the shape of heavy-tailed distributions (Figures 9 and 10). The surge duration varied from 2 to 20 years, with a median surge duration of 6.5 years (Figure 9a). The total terminusadvancing distances ranged from 152.12 m to 3566.98 m, with a median value of 858 m (Figure 9b). Accordingly, the mean terminus-advancing rates are mostly less than 300 m/year (with a median value of 136.76 m/year). We observed generally prominent glacier terminus advance (>300 m/year), mostly concentrated in the southeastern part ( Figure 10). The maximum terminus-advancing rate and the longest advancing distance occurred at glacier G077544E35215N (location flagged in Figure 10), which had advanced 3566.98 m at 1188.99 m/year in about three years.

Geometry Characteristics of Surge-Type Glaciers
Using the dataset of RGI v6.0, we further investigated the characteristics of the 71 surgetype glaciers in the fields of the elevation, aspect, slope, size, and length in the Karakoram. The scatter distribution plot (Figure 11a) reveals that the mean elevations of the surgetype glaciers are more concentrated in the range of 5000 m to 6000 m than those of the nonsurge-type glaciers. We found that the surge-type glaciers depend on a certain level of topographical drop. The elevation range (the difference between the maximum and minimum elevations of the glacier) of the surge-type glaciers is above 1000 m, with more than half (52%) of them being over 2000 m (Figure 11b), while most nonsurge-type glaciers have an elevation range less than 1000 m. The surge-type glaciers are mainly dominated by the northeast orientation, accounting for 31% of the identified surge-type glaciers, and they are evenly distributed in the other orientations (Figure 11c). For the slope, the surge-type glaciers are characterized by lower slopes (median: 25 • ) than those of the nonsurge-type glaciers (median: 31 • ) (Figure 11d). In terms of the area and length (Figure 11e), the surge-type glaciers are larger and longer than the nonsurge-type glaciers, with a mean area and length of 97 km 2 and 18 km, respectively. A few relatively small and short surge-type glaciers are distributed at higher altitudes. As surge-type glaciers have a bias towards greater lengths and low slopes, they mostly represent the lower-right corner of the triangles formed by the slope-length parameters of all the glacier samples (Figure 11f). These statistics are in line with the enthalpy-balance theory [35], which predicts that surge-type glaciers increase with greater lengths and shallower slopes.

Geometry Characteristics of Surge-Type Glaciers
Using the dataset of RGI v6.0, we further investigated the characteristics of the 71 surge-type glaciers in the fields of the elevation, aspect, slope, size, and length in the Karakoram. The scatter distribution plot (Figure 11a) reveals that the mean elevations of the surge-type glaciers are more concentrated in the range of 5000 m to 6000 m than those of the nonsurge-type glaciers. We found that the surge-type glaciers depend on a certain level of topographical drop. The elevation range (the difference between the maximum and minimum elevations of the glacier) of the surge-type glaciers is above 1000 m, with more than half (52%) of them being over 2000 m (Figure 11b), while most nonsurge-type glaciers have an elevation range less than 1000 m. The surge-type glaciers are mainly dominated by the northeast orientation, accounting for 31% of the identified surge-type glaciers, and they are evenly distributed in the other orientations (Figure 11c). For the slope, the surge-type glaciers are characterized by lower slopes (median: 25°) than those of the nonsurge-type glaciers (median: 31°) (Figure 11d). In terms of the area and length (Figure  11e), the surge-type glaciers are larger and longer than the nonsurge-type glaciers, with a mean area and length of 97 km 2 and 18 km, respectively. A few relatively small and short surge-type glaciers are distributed at higher altitudes. As surge-type glaciers have a bias towards greater lengths and low slopes, they mostly represent the lower-right corner of the triangles formed by the slope-length parameters of all the glacier samples (Figure 11f). These statistics are in line with the enthalpy-balance theory [35], which predicts that surge-type glaciers increase with greater lengths and shallower slopes.

Comparison and Limitation
Previous studies, mostly based on medium-resolution satellite observations, have reported extensive glacier-surging behavior from the late 1960s to 2020 in the Karakoram [9,19,20,23,34,49]. Guillet et al. (2022) [23] compiled a relatively complete inventory of surge-type glaciers for the HMA region from different sources of observations during the same study period. They identified a total of 207 surge-type glaciers in the study region. In comparison, our method detected 80% (by area) of the surge-type glaciers with sizes >10 km 2 . As debris-covered surge-type glaciers tend to have a relatively large size, our method successfully detected most of them (95% by area). Additionally, we identified seven new surge-type glaciers (with a total area of 1190 km 2 and a mean size of 170 km 2 ) that have not been reported before. By merging our results with previous studies [9,19,23,37], we found a total of 217 glaciers exhibiting surging behavior in the Karakoram, covering about half of the total glacier area (49%) in the region. In addition to the inclusion of the new surge-type glaciers, we updated the glacier-surging database with information about the most recent active periods and advancing rates. Such information is not provided in the database of Guillet et al. (2022) [23] because of the limitations of different methods in identifying surging events and the laborious work of compiling information from different studies.
The method we proposed for identifying and characterizing surge-type glaciers relies on analyses of the NDSI trends derived from the MODIS imagery. One main limitation of the method lies in the coarse spatial resolution of the MODIS imagery (500 m), which means that it is challenging to identify surging on narrow glacier tongues (<500 m) and for smallscale snow/ice redistributions (i.e., glacier advance less than 1000 m (two-grid size)). For the efficient detection and noise removal, we filtered the 'surge-related' pixel clusters that consisted of less than three pixels, which can further eliminate certain small-scale surging activities. This probably explains why a large number of small surge-type glaciers were not detected in this study. Additionally, the observation period of the method depends on the length of the NDSI time series. In cases of breakpoint detection, the time series analyses will not detect breakpoints located on the two ends ('h' times of the total length) of the period. For example, by setting the 'h' value as 0.15 in this study, we set the minimal segment size from the breakpoint in the trend component for about three years (0.15 × 20 years). Consequently, surge activities occurring between 2000 and 2003 and 2017 and 2020 could not be identified. By setting a small 'h' parameter, the target observation period can be extended slightly, but it will increase the sensitivity of the results to short-term disturbances, such as noise or seasonal snow-cover changes.
Our method takes advantage of the high temporal resolution of MODIS images and the relatively high spatial resolution of Landsat images for the efficient detection of surging activities over extended regions and in the long temporal span. The principle of the method is mainly related to the land-cover changes, or snow/ice-coverage increases, associated with glacier surging. In certain cases (e.g., surging on fully debris-covered glacier tongues, which causes no significant snow/ice-cover changes), the method may fail. However, in comparison with the previous inventories of surge-type glaciers [23], our results successfully detected the majority of the debris-covered surge-type glaciers. We also provide the terminus changes of the clean-ice-covered glacier surging and the active periods by using long-term Landsat images. Note that the active period defined here differs from that determined by the surface-velocity changes. However, it offers a reliable proxy on the velocity change, active period, and magnitude of surging.
Besides the location, active period, and terminus changes, other metrics on glacier surging, such as the velocity change, glacier surface elevation, and geometrical changes, are not documented in this study. Particularly, for the surging of debris-covered glaciers with no significant advance, velocity measurements are required to define the active period. Based on our inventoried surge-type glaciers and surging activities, these metrics can be further investigated with velocity data and geodetic surface-elevation observations.

Characteristics of Glacier Surging and Implications
The surge-type glaciers in the Karakoram tend to have longer glacier tongues (>10 km) and milder slopes (<30 degrees) than the non-surge glaciers, which is consistent with the geometry characteristics of global surge-type glaciers [23,34], and agrees with the implications of the recent enthalpy-balance theory [35]. We reported a median active period of the surge-type clean-ice glaciers of 6.5 years (ranging from 2 to 20 years), which is in line with the typical active period of from 4 to 9 years, as suggested in Paul (2020) [71], but we are missing documents about the active periods fewer than 2 years. The differences may partly arise from the fact that our results do not include a number of small-scale glaciersurging activities due to the spatial-resolution limitation. Note that our determination of an active period could be slightly different from the velocity-based definition and could exclude fully debris-covered surge-type glaciers. Our results confirm the high variability (heavy-tailed distributions) of the surging parameters, including the length of the active period, and the magnitudes (maximum advancing length, mean advancing rates) of the surging. However, we observed the years of surge initiation concentrated in 2008-2010 and 2013-2015. According to the enthalpy theory of Benn et al. (2019) [35], glacier surging arises for a wide variety of geophysical reasons, including the climate (precipitation, temperature) and bed properties (hydraulic conductivity), as well as the geometrical conditions (length, slope), while climatic factors may influence the time and/or magnitude of the surging. Our observations agree with the predictions and implications of this theory, and they suggest that the high concentration of surging events in the past decades may be related to climatic factors. Further studies with velocity observations and the glacier-surge model proposed by Benn et al. (2019) [35] are necessary to examine the energy dissipation during the active phase of a surge to better understand the glacier-surging evolution and the relationship with climate.

Conclusions
The identification of surging glaciers has been challenging due to the difficulties of monitoring long-term glacier behavior, which has led to the limited understanding of the distribution, characteristics, and mechanism of intensified glacier surging in the Karakoram, which is a hotspot of surge-type glaciers in High-Mountain Asia. Here, we present a new method of detecting glacier-surging activities over extended regions by detecting the trends and breakpoints of NDSI time series, which are available globally and almost daily from the MODIS imagery. We developed two detailed criteria (significant upward trends and valid breakpoints) to identify the surging-related signals from the NDSI-trend analysis implemented in the BFAST tool. The surging signals (time and location information) were further validated with the Landsat imagery and geodetic surface-elevation-change information. We identified 71 surge-type glaciers with 74 surging activities from 2000 to 2020 in Karakoram, about half (33) of which are debris-covered glaciers. The surge-type glaciers showed a high clustering in the central and southeastern parts of the study area, and they tend to have long glacier tongues and mild slopes, distributed at high altitudes (between 5000 masl and 6000 masl), and with considerable topographical drop. We further provide detailed information, including the terminus-advancing distances/rates, for the active phase during 2000 and 2020. These metrics reveal a long active period of surge-type glaciers (2-20 years) and the high diversity in the surging characteristics in this region. These geometry characteristics and different metrics of surging information benefit the understanding of the conditions of the glacier-surging developments, as well as future investigations of the surging characteristics in this region for the better assessment and monitoring of glaciological hazards. Due to the coarse spatial resolution of MODIS images, the ability to detect small surging glaciers with narrow tongues may be limited by this method. However, our method can efficiently detect mesoscale surge-type glaciers and derive the related surging information at the regional scale, and it thus has the potential to be applied in other glacierized areas for the near-real-time detection of the occurrence of glacier surging.