remote sensing The Evolution of the Glacier Surges in the Tuanjie Peak, the Qilian Mountains

: Glacier surges (GSs) are a manifestation of glacier instability and one of the most striking phenomena in the mountain cryosphere. Here, we utilize optical images acquired between 1973 and 2021 to map changes in glacier surface velocity and morphology and characterize differences in surface elevation using multi-source DEMs in the Tuanjie Peak (TJP), located in the Qilian Mountains (QLMs). These data provide valuable insights into the recent dynamic evolution of glaciers and hint at how they might evolve in the next few years. We identiﬁed a conﬁrmed surge-type glacier (STG), three likely STGs, and three possible STGs. Our observations show that TJP GSs are generally long-term, although they are shorter in some cases. During the active phase, all glaciers exhibit thickened reservoir areas and thinned receiving areas, or vice-versa. The ice volume transfer was between 0.11 ± 0.13 × 10 7 m 3 to 5.71 ± 0.69 × 10 7 m 3 . Although it was impossible to obtain integrated velocity proﬁles throughout the glacier surge process due to the limitations of available satellite imagery, our recent observations show that winter velocities were much higher than summer velocities, suggesting an obvious correlation between surge dynamics and glacial hydrology. However, the initiation and termination phase of GSs in this region was slow, which is similar to Svalbard-type STGs. We hypothesize that both thermal and hydrological controls are crucial. Moreover, we suggest that the regional warming trend may potentially increase glacier instability and the possibility of surge occurrence in this region.


Introduction
Glaciers are recognized as important variables in the global climate system [1]. Most glaciers in High Mountain Asia (HMA) have retreated significantly against the backdrop of recent climate warming [2]. In contrast, some glaciers have remained stable, advanced, or even surged during the same period within HMA sub-regions [3]. GSs are quasi-periodic oscillations of a glacier's dynamical behavior, during which flow velocity reaches up to 10-1000 times its standard order of magnitude [4]. During surges, a substantial volume of ice is rapidly transferred from the reservoir to the receiving area, leading to dramatic changes in the morphology of the glacier surface [5]. On the other hand, the quiescent periods usually last for tens to hundreds of years, during which mass accumulates in the upper glacier area [6]. With the progress of remote sensing observation technology [7],

Study Site
The QLMs are located on the northeastern edge of the Qinghai-Tibet Plateau (QTP) and are composed of a series of parallel mountains and valleys, which run roughly northwest (Figure 1), starting from Wushaoling in the east, ending at Dangjinshan Mouth in the west, bordering the Qaidam Basin to the south and the Hexi Corridor in the north, with a total length of around 800 km and a width of around 300 km [36]. The QLMs belong to the plateau continental climate zone; the western section is controlled by the Westerly circulation, while the eastern section is affected by the southeast monsoon and the southwest monsoon passing over the QTP [37]. According to the Second Chinese Glacier Inventory (SCGI), the QLMs contained 2684 glaciers covering a total area of 1597.81 ± 70.30 km 2 with an ice volume of~84.48 km 3 between 2005 and 2010 [38]. The terrain elevation rises gradually from northeast to southwest, with the highest elevations at TJP. The TJP region (36.5 • -39 • N, 93.5 • -103 • E) is located on the northwestern fringe of the QLMs (Figure 1). The peak elevation is 5826 m a.s.l., with an average firn line at about 4862 m a.s.l. during late summer. The average annual precipitation and temperature at 4000 m a.s.l. are approximately 300 mm and −6.3 • C, respectively [35].
Remote Sens. 2022, 14, x FOR PEER REVIEW 3 of 20 controlling their behavior. Overall, the results of this study have important implications for understanding glacial dynamics and will aid in the interpretation of heterogeneity in HMA GSs.

Study Site
The QLMs are located on the northeastern edge of the Qinghai-Tibet Plateau (QTP) and are composed of a series of parallel mountains and valleys, which run roughly northwest (Figure 1), starting from Wushaoling in the east, ending at Dangjinshan Mouth in the west, bordering the Qaidam Basin to the south and the Hexi Corridor in the north, with a total length of around 800 km and a width of around 300 km [36]. The QLMs belong to the plateau continental climate zone; the western section is controlled by the Westerly circulation, while the eastern section is affected by the southeast monsoon and the southwest monsoon passing over the QTP [37]. According to the Second Chinese Glacier Inventory (SCGI), the QLMs contained 2684 glaciers covering a total area of 1597.81 ± 70.30 km 2 with an ice volume of ~84.48 km 3 between 2005 and 2010 [38]. The terrain elevation rises gradually from northeast to southwest, with the highest elevations at TJP. The TJP region (36.5°-39°N, 93.5°-103°E) is located on the northwestern fringe of the QLMs (Figure 1). The peak elevation is 5826 m a.s.l., with an average firn line at about 4862 m a.s.l. during late summer. The average annual precipitation and temperature at 4000 m a.s.l. are approximately 300 mm and −6.3 °C, respectively [35].

Glacier Delineation and Glacier Centerline Extraction
We utilized data from the Landsat (https://glovis.usgs.gov, 20 October 2021), Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER, Figure 1. The distribution of STGs in the TJP. Flowlines shown here are used for velocity and elevation profiles. Black and white lines mark the distance from the terminus position in kilometers, plotted at 1 km intervals. Background image is from Landsat TM on 1 August 2006, projected on a WGS84 grid. The white border in the sub-picture indicates the range of the QLMs. We numbered each glacier by combining their orientation. NE and SW indicate the northeast slope and the southwest slope, respectively.

Glacier Delineation and Glacier Centerline Extraction
We utilized data from the Landsat (https://glovis.usgs.gov, accessed on 20 October 2021), Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER, https://eart hexplorer.usgs.gov, accessed on 10 September 2021), and Sentinel-2 (https://scihub.coperni cus.eu, accessed on 2 September 2021) optical satellite systems to map glacier outlines and extract glacier centerlines from 1973 to 2020. All the data are free and open access. We used images with less than 10% cloudiness over the study area (see Table 1). Sentinel-2 Level 1C products include orthoimages containing radiometrically corrected top reflectance values of the atmosphere. These orthoimages are geometrically corrected based on an accurate geometric model [39]. The processing of Landsat and ASTER products is similar to the processing of Sentinel-2 data for radiometric and geometric correction, orthorectification, and resampling to the map grid [40]. We took the SCGI as a reference used manually to measure the changes in glacier termini for each period using ArcMap 10.8. The assumed accuracy of glacier boundary extraction is mainly affected by sensor and image registration errors because interpretation errors and technical errors are difficult to quantify [41]. Thus, we only consider errors resulting from the spatial resolution of satellite remote sensing images in this study, e.g., [42]. To extract the glacial centerlines, we used an automatic method proposed by Zhang et al. [43] A Sentinel-2 image was used to evaluate the accuracy of these centerlines [44], the uncertainties were no more than ±23 m.

Velocity Mapping
Glacier surface velocity was quantified using feature tracking in Co-registration of Optically Sensed Images and Correlation (COSI-Corr) software [45]. This software has been successfully used for glacier velocity extraction [46]. In this study, the annual glacial velocity fields were derived from 1986 to 2020 utilizing Landsat, ASTER, and Sentinel images (See Table S1). To compare the interannual changes for the studied glaciers, the velocity profiles in different years were extracted along the glaciers' centerlines and combined. To investigate seasonal motion variability in greater detail, we incorporated Landsat 8 (OLI) and Sentinel-2 images from 2013 to 2021 to obtain the displacements of several STGs during their surge processes (See Table S2). The horizontal displacement was measured by a frequency subpixel correlator; using a multi-scale method, the correlation window size was changed from 64 × 64 to 32 × 32, sliding every pixel. The measured column and row displacements were combined to acquire the magnitude and direction of glacier surface displacement [25]. The glacier surface velocity (GSV) in meters per year (m a −1 ) was acquired using Equation (1): where GSD is the displacement measured between each image pair, T i is the interval between the image pair in days, and T y is a constant of 365. Due to the limitations of the clear surface features of optical satellite images, the detection of mountain glacier motion inevitably has the problem of decorrelation, e.g., variable snow-cover [26]. To reduce these impacts, we smoothed the displacement map using a Gaussian low pass filter with a kernel size of 3 × 3. We took values that differed by >20 m a −1 or 3 m d −1 relative to those in surrounding areas as errors and manually removed them from the resulting data. The cross-correlation algorithm was used to extract the glacier surface velocity, the errors of which include system error, image quality, and image registration error [47]. To quantify the error, assuming that the non-glacier region is stable and no displacement changes occur, the velocity error can be evaluated in the non-glacier region, according to the following formula [48]: where e off is the displacement error of the non-glacier region, MED is the average displacement, and SE is the standard deviation of the average displacement. The value of SE is calculated as: In this case, STDV is the standard deviation of the average displacement in the nonglacier region, and N eff represents the number of effective pixels from which autocorrelation is removed .
Here, PS is the pixel resolution, and D is the spatial autocorrelation distance to remove the influence of autocorrelation, which is generally 20 times the pixel resolution. Tables S1 and S2 show the mean displacement and corresponding standard deviation values in the non-glacier region. The accuracy of the velocity data depends on the time interval between images. The error ranges of the interannual and seasonal velocities are 1.03 m a −1 to 6.63 m a −1 and 0.04 m d −1 to 2.56 m d −1 , respectively.

Glacier Surface Elevation Changes
Surface elevation changes were documented using four ASTERs, one Shuttle Radar Topography Mission (SRTM), and four 1:50,000 topographic (TOPO) maps. We used MMAS-TER [49] to create ASTER DEMs with ∼10 m vertical uncertainty [39] and 30 m spatial resolution using stereo imagery data from 02 November 2002 (ASD02), 16 November 2007 (ASD07), 28 October 2012 (ASD12), and 10 March 2018 (ASD18). The SRTM mission collected near-global data at 30 m-resolution from 11 to 22 February 2000 using C-band radar [50]. The SRTM DEM is provided by the USGS and is non-void-filled in our study area. The vertical and horizontal accuracies can be better than 10 m [51]. Four TOPO maps were used at a scale of 1:50,000, compiled from aerial photos taken in 1966 by the Chinese Military Geodetic Service. The nominal vertical accuracy of TOPO is within 3-5 m for flat and hilly areas and 8-14 m for mountainsides and high mountain areas [52]. These data have previously been successfully applied to regional-scale glacier elevation change monitoring on the QTP e.g., [44]. We manually digitized the contours and then converted them into a raster DEM (TOPO DEM) with a 30 m grid cell size using the Thiessen polygon method in ArcMap 10.8 software. The surface elevation accuracy of the TOPO DEM, ASTER DEM, and SRTM DEMs for the TJP was tested by comparing their values over the non-glacier region with the Advanced Topographic Laser Altimeter System (ATLAS) product of the National Aeronautics and Space Administration's Ice, Cloud, and Land Elevation Satellite (http://nsidc.org/data/, accessed on 17 September 2021). The surface elevation difference over the non-glacier region between ATLAS and the TOPO DEM was 2.52 ± 4.71 m (mean value ± standard deviation) and 1.48 ± 2.06 m between ATLAS and SRTM. The surface elevation differences between the ATLAS and ASTER DEMs are 1.94 ± 3.31 m, 3.07 ± 5.68 m, 1.84 ± 2.45 m, and 2.52 ± 4.71 m, from oldest to most recent.
Multitemporal DEMs require co-registration to remove horizontal and vertical offsets before any difference analysis, as DEMs generated from different data sources may have inconsistent geolocations resulting from different postprocessing procedures and methodological limitations, etc. [53]. Here, we co-registered the TOPO and ASTER DEMs to the baseline SRTM DEM using the three-step correction framework of Nuth and Kääb [54]. See literature [55] for specific steps. Eventually, glacier surface elevation changes were obtained by using DEM differencing. It is difficult to evaluate DEM data errors through field measurements because of harsh glacier environments and poor accessibility [56]. Assuming that the heights in non-glacier areas remain essentially unchanged for long periods, the mean elevation differences (MED) and standard deviation (SD) in non-glacier areas can be used to estimate the glacier height change errors [26]. In order to reduce the autocorrelation error, a distance value of 600 m was chosen for the difference map derived from TOPO DEM, SRTM DEM, and ASTER DEM [57]. Therefore, the overall errors of the changes in glacier surface elevation (δ) ( Table 2) can be calculated by: where N is the number of pixels with effective measurements. Radar signals can penetrate into snow and ice up to several meters [52]. The penetration depth can range from 0 to 10 m depending on a variety of parameters [57]. The datasets of the SRTM C/X-Band radar penetration depth differences on the 1 • × 1 • grid of HMA glaciers (2000) were obtained from the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn/en/, accessed on 22 September 2021) in this study [58]. These datasets show that the difference of the SRTM C/X-band penetration depth in the HMA region is between 0 m and 7.64 m, with an average value of 2.40 ± 0.04 m. Differences between the SRTM-C and X-bands showed an average C-band penetration depth of 5.54 ± 0.91 m in the TJP region.

Identification of STGs
The most direct evidence of GSs is sudden acceleration. However, due to the limitation of data availability, it is necessary to consider the long sequence variation of the surface elevation and frontal position when identifying surges [4].
In this paper, we used the following criteria for the identification of STGs: I Velocities during the active phase typically reach at least an order of magnitude higher than during the passive, quiescent phase, which was classified as a confirmed STG. For example, the glacier NE1.
II Obvious thinned reservoir area and thickened receiving area accompanied by a sudden and rapid advance of the frontal position, out of synchrony with the behavior of neighboring glaciers, was classified as a likely STG. For example, the glaciers NE3, NE4, and NE5.
III Clear thickened reservoir area and thinned receiving area accompanied by a glaciological/geomorphological evidence: looped moraines, heavy surface crevassing, push moraines, eskers, etc., was classified as a possible STG. For example, the glacier, NE2, SW1, and SW2.

STGs Overview
We identified a confirmed STG, three likely STGs, and three possible STGs from the studied glaciers by extracting glacier surface velocity, analyzing glacier surface elevation changes, and manually inspecting glacier surface morphology. The average area of these glaciers is 8.22 km 2 , and the corresponding average maximum length and average slope of the glacier surface are 6.11 km and 19.14 • , respectively. The glaciers are distributed on the northeast and southwest slopes of TJP (Table 3 and Figure 1). The proportion of debris in the glacier areas does not exceed 6%. The identified surges generally exhibit two styles by whether the surge impacts the glacier terminus change or not. Four GSs led to a terminus advance (NE1, NE3, NE4, and NE5) between 1973 and 2020 (Table 4). According to the change of glacier length data shown in Table 4, glaciers NE1 and NE3 advanced more than 500 m, with surge durations varying from 5 years (NE1) to nearly 20 years (NE3). Notably, although the terminal of glacier NE5 advanced by 114.65 ± 14.21 m, the area of the glacier was reduced by −0.03 ± 0.37 km 2 ( Table 4). Three glaciers surged without affecting their termini (NE2, SW1, and SW2); in contrast, their termini are retreating. From the perspective of the initiation time of the GSs, one GS started before 2000 (NE3), while the others all occurred after 2000.

Glacier Surface Elevation and Ice Volume Changes
In total, five elevation change maps were generated from the TOPO, SRTM, and ASTER DEMs. Between 1966 and 2018, all glaciers exhibited thickened reservoir areas and thinned receiving areas (or thinned reservoir areas and thickened receiving areas). The glacier NE1 underwent the maximum transfer of ice volume, with a net volume change of 5.71 ± 0.69 × 10 7 m 3 . This glacier also experienced a mean surface elevation increase of up to 53.10 ± 0.64 m at the glacier terminus and a mean surface lowering of −22.24 ± 0.64 m in the reservoir area (Table 5). This elevation change coincided with the advance of the glacier terminus ( Table 4). The largest surface elevation decreased along the centerline occurred between 1.8 and 4.9 km from the glacier terminus, and the largest surface elevation increase occurred between 0 and 1.8 km (Figure 2). In comparison, glacier NE2 had the smallest ice volume shifts ( Table 5). The redistribution of surface mass on glaciers NE3 and NE4 both caused the glacier termini to advance ( Figure 2 and Table 4). The average ice surface elevation of glacier NE3 decreased by −29.93 ± 1.07 m at a distance between 0.9 km and 4.6 km from the end of the glacier, while the average thickness of the glacier terminus increased by 45.07 ± 1.07 m ( Table 5). Compared with NE3, glacier NE4 exhibited a smaller change in surface elevation, resulting in only 0.34 ± 0.07 × 10 7 m 3 of ice volume transfer ( Table 5). The average elevation change of the NE5 glacier surface is almost the same in the reservoir area and the receiving area (Table 5). Although the SW1 and SW2 glaciers did not undergo ice volume transfer, both glaciers showed a thinning of reservoir areas and a thickening of receiving areas between 2012 and 2018 ( Figure 2 and Table 5).

Ice Surface Velocities
The annual velocity of the glacier centerline profile provides quantitative information about the surface displacement during the GSs.  Figure 3c,e, initiation of both glaciers NE2 and NE4 was slow; their peak velocities eventually reached 40 ± 4.54 m a −1 and 46 ± 4.54 m a −1 , respectively, in 2019. By contrast, the surge initiation of glaciers NE5 and SW1 in 2012 was sudden (Figure 3f,g). Interestingly, the glacier terminus of SW1 was not affected. Changes in the glacier terminus of NE3 were indicative of surge behavior ( Table 4). The profiles of NE3 show an active terminal, although our results are limited to the decay of the surge (Figure 3d). The peak velocity of this glacier reached~17 ± 3.0 m a −1 m in 1989 and then gradually decreased.
The final sub-figure shows the dynamic characteristics of glacier SW2 between 1987 and 2020 ( Figure 3h). The profiles show a clear downstream acceleration process between 1987 and 2000. Moreover, the recent accumulation of ice mass in the glacier reservoir area allowed it to reach its peak velocity (~68 ± 4.62 m a −1 ) in 2020.

Ice Surface Velocities
The annual velocity of the glacier centerline profile provides quantitative information about the surface displacement during the GSs.    (Figure 4a). The ice crevasses in the middle of glacier NE2 showed clear and dramatic changes over relatively short periods (Figure 4b). The terminal morphology of glacier NE3 changed significantly between 25 July 1986 and 15 July 2000, with a distinct narrowing of the downstream area of the glacier observed (Figure 4c). Glaciers NE4 and NE5 were all partially debris-covered. In general, GSs destroy the original surface morphology and remold both the glaciers and debris during the active phase. Compared with glacier NE4, the lateral ice crevasse in glacier NE5 s reservoir area markedly increased (Figure 4d,e). Although the surges of glaciers SW1 and SW2 did not lead to terminus changes, they reshaped the surface structures of the glaciers. For example, Figure 4f shows clear movement of the ice cliff on the surface of the glacier SW1.

Glacier Surface Morphology Changes
Glacier surface morphology changes were analyzed using ASTER (  glacier NE4, the lateral ice crevasse in glacier NE5′s reservoir area markedly increased (Figure 4d,e). Although the surges of glaciers SW1 and SW2 did not lead to terminus changes, they reshaped the surface structures of the glaciers. For example, Figure 4f shows clear movement of the ice cliff on the surface of the glacier SW1.

Surge Characteristics
Several independent lines of evidence strongly indicate that the glaciers identified in this paper are STGs, including undoubtedly, likely, and possibly. Our data suggest that the seven glaciers have experienced a total of seven surge events during the last 34 years. Of these, four glaciers are currently in the active phase. Although there are gaps in the optical remote sensing image dataset between 1973 and 1986, we anticipate that glacier surge durations last for more than 10 years (except for glacier NE1) in the TJP of the QLMs.
Glacier surge characteristics are important for our understanding of HMA glacier instability [11]. The STGs in the 'Karakorum anomaly' region was the most widely distributed in the HMA. Among these, glaciers with relatively short (≤5 years) active phases are widespread throughout the Pamirs [9]. In the Karakoram, the active phases range from several months to over 15 years [8]. In contrast, the active period of GSs in the West Kunlun Mountains is more than 5 years [10]. The GSs in the inner QTP last for no more than 10 years (Table 6). However, the duration of the TJP GSs recorded in this study is similar to those of the Geladandong glaciers, ranging from 3 to 20 years (Table 6). This may also imply that STGs have a relatively long recurrence interval in the TJP area; unfortunately, we currently have insufficient information to constrain this hypothesis further. In general, the surge cycle is several decades in the HMA. However, the whole surge period is only observed for a few glaciers; thus, this impression is potentially biased by the observation window [9]. In addition, as shown in Table 5, the magnitude of the STGs' maximum ice surface velocities and terminal advance in TJP are smaller than those of other surge clusters. However, this may be related to the area of the glaciers. Notably, we do not consider detachments of low-angle mountain glaciers here, although this behavior is similar to GSs [59]. Table 6. Averaged geometric information (area, slope, maximum length, and maximum/minimum/ medium elevation) of glaciers in nine surging regions and the summary of the surge duration, period, maximum velocity (m a −1 ), and terminus advance (km) for STGs in different regions of HMA [3,[8][9][10][11]15,25,[32][33][34]. We use the symbol-when no values are specified. A previous study found that some glacier geometry variables (e.g., surface slope and glacier length) show statistically significant correlations with surging behavior [15]. However, to date, only three regions have been researched by glacier inventory, including Karakoram [8], Pamirs [9], and Tian Shan [12]. To advance the understanding of glacier surge characteristics and extend our predictive capabilities of regional-scale changes, using the outlines from RGI V6 [60], we analyzed differences in area, slope, range of elevation, and length between STGs and non-STGs in the HMA. In addition, we also comprehensively compared the surge features of the STGs in the active phase. Table 6 shows the mean results for all geometry attributes for the main surge regions. In all regions where both non-STGs and STGs coexist, the latter exhibits larger values. The largest STGs are found in the West Kunlun Mountains, with a mean area of 89.41 km 2 compared to 2.46 km 2 for non-STGs. The smallest size differences between non-STGs and STGs were found in the TJP, where STGs are still larger than non-STGs, but the mean difference in size between the two is <7.5 km 2 . According to Table 6, STGs are longer than non-STGs, similar to the pattern shown by glacier area. The mean slope is also an obvious difference between the two populations; however, in this case, an inverse relationship is shown, where STGs have shallower slopes than their non-STGs equivalents. Interestingly, we observed that the average slope of STGs in the inner QTP is less than 14 • , such as Puruogangri, Xinqingfeng and Malan, Geladandong and Muztag, etc. (Table 6) The TJP have an average slope of STGs of 19.14 • . All other main surge regions have slopes greater than 22 • due to the extremely steep and elevated topography [25]. The glacier elevation range is related to the glacier length. Overall, the STGs' average median elevation at TJP is only higher than those of the Karakoram and Tian Shan glaciers (Table 6).

Motion Patterns of STGs
The existence of seven STGs on TJP in the QLMs strongly proves that the HMA region is one of the most active surge regions in the world. Except for NE1, none of the glaciers detailed here have been previously identified as STGs. To better understand the movement patterns of STGs in this region, we analyzed the surface velocities of four glaciers (NE2, NE4, NE5, SW1, and SW2) from 2013 to 2021. Although the size and shape of these glaciers are different, we observe that their surge patterns are quite similar based on 131 datasets of seasonal surface velocity. Figure 5 shows the mean velocity of the glacier centerlines in the five likely and possibly STGs from 2013 to 2021. The velocity-time series from STGs indicates an obvious correlation between glacial hydrology and surge dynamics. We found that (1)

Motion Patterns of STGs
The existence of seven STGs on TJP in the QLMs strongly proves that the HMA region is one of the most active surge regions in the world. Except for NE1, none of the glaciers detailed here have been previously identified as STGs. To better understand the movement patterns of STGs in this region, we analyzed the surface velocities of four glaciers (NE2, NE4, NE5, SW1, and SW2) from 2013 to 2021. Although the size and shape of these glaciers are different, we observe that their surge patterns are quite similar based on 131 datasets of seasonal surface velocity. Figure 5 shows the mean velocity of the glacier centerlines in the five likely and possibly STGs from 2013 to 2021. The velocity-time series from STGs indicates an obvious correlation between glacial hydrology and surge dynamics. We found that (1) the glacier is in a relatively stable stage before the onset of the surge, where the mean velocities in summer are greater than those in winter; (2) a gradual increase in winter speeds in 2016 and 2017; (3) the peak velocity was reached in winter 2018; and (4) subsequent gradual velocity decreases. When observations with a sufficiently high time resolution are available, the velocity fluctuations during surges seem to be usually recorded. Previous studies have reported seasonal velocity changes during GSs, such as seasonal-scale observations, or even changes within hours, of the Variegated glacier in 1982-1983 [61]. These studies have further established an obvious connection between velocity fluctuations and glacier hydrology during surge events [61]. GSs in the Karakoram are widespread and universal [9]. Quincey et al. [29] observed that the velocity of the Karakoram STGs, such as the Shakesiga glacier, changed from slow to fast movement in winter. Interestingly, the GSs of Shispare and Khurdopin accelerated in late spring and early summer in the Karakoram, which may hint at the involvement of glacier meltwater [27]. Similarly, two-phase surges of peak velocity associated with surface meltwater production have also been found in Kyagar Glacier [28]. However, some researchers have argued that the emergence of GSs in the Karakoram largely depends on the high elevation, complex topography, and climatic background [9,29]. Another example of seasonal velocity changes during surges is When observations with a sufficiently high time resolution are available, the velocity fluctuations during surges seem to be usually recorded. Previous studies have reported seasonal velocity changes during GSs, such as seasonal-scale observations, or even changes within hours, of the Variegated glacier in 1982-1983 [61]. These studies have further established an obvious connection between velocity fluctuations and glacier hydrology during surge events [61]. GSs in the Karakoram are widespread and universal [9]. Quincey et al. [29] observed that the velocity of the Karakoram STGs, such as the Shakesiga glacier, changed from slow to fast movement in winter. Interestingly, the GSs of Shispare and Khurdopin accelerated in late spring and early summer in the Karakoram, which may hint at the involvement of glacier meltwater [27]. Similarly, two-phase surges of peak velocity associated with surface meltwater production have also been found in Kyagar Glacier [28]. However, some researchers have argued that the emergence of GSs in the Karakoram largely depends on the high elevation, complex topography, and climatic background [9,29]. Another example of seasonal velocity changes during surges is documented in the West Kunlun Mountains, where winter speed-ups are also recorded [30]. In contrast, a gradual increase in summer velocities leading up to the surge of Monomah Glacier, Central Kunlun Mountain Range, was reported by Guo and others [62]. Although we only obtained seasonal velocity changes from four glaciers because of the availability of satellite data, the remaining glaciers also showed signals of winter initiation; for example, on 25 April 2002, the NE1 glacier clearly showed a ring-shaped terminal connected with the east branch ( Figure 4).
Raymond [63] suggested that hydrologically controlled surges showed acceleration during winter and deceleration during summer. These underlying mechanisms would mean that the acceleration and deceleration processes are rapid. Thermally controlled surges are the opposite [19]. Even though our data provided a seasonal signal of the surge evolution, the surges clearly reached their peak velocities in the several years following initiation. In contrast to hydrologically regulated GSs recorded in Alaska (where the initiation phase lasted as little as several months), the rates recorded in our study seem to be slow, as the acceleration and deceleration phases of the GSs of TJP are relatively long-lasting ( Figure 3). The possible differences between these GSs may challenge the view that glaciers in different climate regions have different surge mechanisms, rather than suggest that a single mechanism has different manifestations [29]. An interesting result of the study is that using only Landsat 8 and Sentinel-2 data from 2013 to 2021 would have resulted in an imperfect map of the dynamics of STGs in the TJP, as these data would ignore the surge characteristics of other glaciers, especially glacier NE3. Therefore, we believe that the uncertainty in the mechanism and evolution of specific GSs in different regions of the world is partly due to the limitations of the observations.

Surge Trigger Mechanisms
Two principal, broadly recognized hypotheses of surge triggering mechanisms have been proposed: hydraulic and thermal [39]. The hydraulic mechanism consists of an event that causes the base water pressure to increase, reduces bed resistance, and causes a rapid slip to occur; this trigger type is expected to produce more sudden behavior [16]. The thermal trigger represents an event involving the thawing of a previously frozen glacier bed, enabling widespread sliding to take place, and is characterized by a slow build-up and termination [19]. Previous studies on GSs in western China have shown that the control process of each glacier changes due to its different geomorphological characteristics, hydrological conditions, and thermal conditions [25]. No single process has been reported to be responsible for their instability.
In the TJP, the observed winter speed-up and summer slowdown during the active surging phase can be readily interpreted in terms of causal mechanism, as these processes resemble hydrological rather than thermal regulation, although the long active phase distinguishes these events from hydrologically controlled surges in Alaska [64]. In addition, the thickening of the glacier accumulation area may lead to the occurrence of basal sliding (Table 5). Notably, this feature is more consistent with previous reports of the thermal regulation mechanism in Svalbard [5]. Due to limitations of the available satellite images from 1970 to 2012, the apparent seasonal variations of glaciers NE1 and NE3 could not be detected. We could also not provide evidence of thermal condition changes inside the glacier due to the harsh environment. However, based on variations in their annual velocity, the surge in glacier NE1 is sudden, which was highly similar to Alaska-type STGs. Interestingly, the duration of the glacier NE3 is more than 15 years. This may be similar to Trapridge Glacier, Yukon, Canada, belonging to a "slow surge" [65]. This slow phase may be due to the massive climate-driven accumulation deficit in recent decades that prevented its rapid phase from progressing [65]. Therefore, we suggest that there is no unified mechanism for GSs in the TJP. A previous study has also suggested that surges in the eastern Pamir [25], Karakoram [29], and West Kunlun Mountains [10] are triggered by a blend of hydrological and thermal regulation.
Furthermore, till deformation was considered an important process of surge generation [10]. For instance, the finding GSs of Gilbert and others [66] suggests the subglacial sediment played a key role in the two giant glacier collapses documented to date in Tibet. However, the glacier till exhibits plastic rheology, and its shear strength strongly depends on the high basal water pressure [19]. The high sensitivity of the ultimate shear strength of the substrate to pore water pressure makes the friction conditions of the substrate change drastically and continuously, which leads to this type of instability [66]. In other words, irrespective of the presence of a till layer beneath some STGs in the TJP area, the surge trigger mechanism should include a process that can provide high basal water pressures. The mechanisms of hydrological and thermal controls are absolutely different in terms of their basal water pressure generation [10]. The seasonal evolution of flow velocity in Alaska-type STGs proves that the basal water originates from surface meltwater in the hydrological regulation mechanism [67]. In contrast, in terms of the trigger processes for the thermal regulation mechanism, the basal water originates fully from the bottom of the ice based on observations from surges in polythermal glaciers in Svalbard [5]. The presence of seasonal changes in surface velocities shows the impact of surface meltwater in TJP. In summer, new crevasses are formed, which can transport surface meltwater to the glacier bed. Although the inflow of meltwater is increased, the increased void volume and cracks can store water, resulting in a drop in water pressure and a slower surface velocity in summer ( Figure 5).
Despite this mechanism appearing to be the opposite of the proverbial spring-summer acceleration [8], the base water pressure is suggested to be not only proportional to the total englacial water volume, but also inversely proportional to the macroscopic porosity [68]. Moreover, the inflow of surface meltwater may be smaller than the glaciers, which may prevent channel development. Since the channel may not exist or only develop very weakly, the summer meltwater may not be discharged as effectively during GSs [69]; therefore, a portion of the meltwater would remain until late fall to winter. This may imply that due to the reduced macroscopic porosity of the glacier surface, higher basal water pressure is reached in winter [10]. The TJP surge velocity speeds are much slower than those observed in the Karakoram GSs, which could be due to the smaller englacial and subglacial water volumes inherent to the relatively colder and drier climate [35]. On the other hand, the glaciers in the TJP are relatively small. Therefore, there is sufficient ice mass in the accumulation area to successfully repeat the seasonal cycle of the surge period [19], which may potentially explain the longer duration of the active period in this area.
Although the surge in climate control activities is controversial, the thickening of the reservoir area in the TJP must be associated with a meteorological forcing. Previous studies have found that the precipitation in this region showed minimal change before 2000; however, after 2000, the precipitation increased at a rate of 60.40 mm/10a [70]. This is consistent with the recorded increase in summer precipitation since the mid-1990s in the Inner Tibetan Plateau due to the weakening of the Westerlies [71]. Compared with other parts of the QLMs, the forcing mechanism must favor greater accumulation [72]. Over time, the accumulation zone will accept more ice mass-as a result, the bottom part of the glacier may reach the pressure melting point, which produces meltwater and causes glacier sliding to occur [73]. At the beginning of the base slip, the frictional heating caused by the ice movement is thought to produce additional meltwater [30], thus triggering GSs. This is consistent with the thermal regulation mechanism of the Svalbard-type GSs. In addition, warming has been observed at a mean annual rate of 0.46 • C/10a during the period 1979-2018 at 150 stations over the Tibetan Plateau [74]. The TJP area is no exception, where an overall temperature increase has been recorded at a rate of 0.39 • C/10a [70]. Regional warming may have an impact on the thermal conditions in and under the glacier, leading to changes in the thermal state of the glacier. For example, a recent modeling study found that Laohugou Glacier No. 12, located around 100 km away from TJP, is probably polythermal; however, it was previously considered fully cold [75]. Based on the above discussion, we believe that regional warming and increasing precipitation trends may increase the instability of glaciers in the region and may cause surges.

Conclusions
In this study, we report a confirmed surge-type glacier (STG), three likely STGs, and three possibly STGs in the TJP. Using a combination of manual digitization, feature tracking, and DEM differencing, we mapped the velocity, surface elevation, and morphology of the seven STGs between 1966 and 2021. STGs are distributed in the TJP region and thus should be included in mass balance studies based on changes in glacier areas or ice thicknesses, etc. Our data demonstrate that the TJP GSs are generally long-term, lasting more than eight years from initiation to termination, although they are shorter in some cases. These characteristics classify the recorded surges as Svalbard-type GSs. Such surges have been assumed to be controlled by the thermal regulation mechanism; however, the seasonal velocity pattern appears to be more strongly related to the hydrological control mechanism. Therefore, the dynamic evolution of TJP GSs does not fit clearly with the established paradigms for hydrologically driven or thermally driven surge mechanisms. We speculate that climate change may complicate the interpretation of the mechanisms of the region's GSs. To thoroughly understand surge mechanisms, the regular in situ monitoring of climate, glacier motions, thermal regime, drainage systems, and basal conditions is required in the future.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14040852/s1, Table S1: Landsat 5 TM (LT05), Landsat 7 ETM+ (LE07), Landsat 8 OLI (LC08), ASTER (AST) and Sentinel-2 (S2) scenes used to produce the annual velocity maps. Scene pair dates (DD/MM/YYYY) is also listed. Bias and uncertainty metrics for the magnitude components of the velocity are listed. MED is mean surface velocity difference in stable terrain area, STDV is standard deviation. The spatial autocorrelation distance (D) is 2400 m. The resolution (PS) of the interannual velocity is 120 m; Table S2: Landsat 8 OLI (LC08) and Sentinel-2 (S2) scenes used to produce the seasonal velocity maps. Scene pair dates (DD/MM/YYYY) is also listed. Bias and uncertainty metrics for the magnitude components of the velocity are listed. MED is mean surface velocity difference in stable terrain area, STDV is standard deviation. The spatial autocorrelation distance (D) is 600 m. The resolution (PS) of the interannual velocity is 30 m.
Author Contributions: Conceptualization, S.L., X.Y. and Y.G.; data curation, Y.G. and F.X.; methodology, Y.Z., K.W., F.X. and Y.G.; validation, M.Q. and Y.G.; visualization, Y.G. and M.S. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: The study did not report any data.