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

Abstract

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 identified a confirmed 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⁷ m³ to 5.71 ± 0.69 × 10⁷ m³. Although it was impossible to obtain integrated velocity profiles 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.

1. 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], STGs throughout HMA have received significant attention in the past decades, for example, the Karakoram [8], Pamirs [9], West Kunlun Shan [10], Himalayas [11], and Tien Shan areas [12]. Although STGs account for only a small part of the HMA’s glacier population, they can provide valuable insights into the triggering mechanisms of surge-type behavior and fast glacier flow [13]. Moreover, the cyclical terminal advance of STGs can result in repeated glacial lake outburst flood hazard formation [14]. Furthermore, understanding the evolution of STGs is crucial to separate internal glacier dynamics from climate change signals. This further highlights the necessity to understand HMA glacier dynamics [15]. Therefore, exploring the surge mechanisms is crucial to preventing such glacial hazards, and, furthermore, it is also an effective way to reveal glacier movement responses to climate change in the HMA.

In the current literature, two main hypotheses have been used to explain glacier instability mechanisms. In the first case, surging occurs when the subglacial hydrological system switches to an ineffective cavity system, leading to increased pore water pressures on the glacier bed, accelerating fast basal sliding [16]. Such events are mostly reported in Alaska [17]. The second hypothesis, the thermal switch model, suggests that when a glacier thickens by accumulation, the increasing sequestration of geothermal heat warms its bed to the pressure-melting point; the ensuing meltwater production then weakens the bed to trigger a sliding motion [18]. This fast flow would have the effect of thinning the glacier and augmenting conductive cooling of the bed, leading to freezing that would terminate the surge [19]. This event type is mostly recognized in Svalbard, Yukon, and East Greenland [20]. In this scenario, the GSs can initiate during any season and tend to have long surge and quiescent periods [21]. In addition, Sevestre and Benn [15] proposed a general theory of GSs. Benn et al. [22] systematically refined the theory and developed a lumped parameter model, which parameterized key thermodynamic and hydrological processes.

STGs are unevenly distributed in the HMA [15]. Previous studies have focused on glacier-scale and regional-scale observations. In Tien Shan, Mukherjee and others [23] identified 39 potential STGs and 5 tributary surges using the available literature and remotely sensed images from 1960 to 2014. Unfortunately, surge mechanism studies are limited in this region. In contrast, STGs in the Pamir region have been described in detail, which indicates a hydrologically regulated surge mechanism at some STGs by observations of repeated short-lived surge cycles [24]. Some studies have also suggested that both thermal and hydrological controls are significant in the Pamir [25]. Surging has been widely studied at either a regional scale or the glacier scale across the Karakoram, e.g., the Hunza Basin [26], Khurdopin glacier [27], and Kyagar glacier [28]. Recently, Bhambri and others [8] produced an inventory of surging in the Karakoram, identifying 221 STGs using a wide variety of observations spanning the period 1840–2017. Reports of surging indicate no single specific trigger mechanism for these glaciers [29]. Yasuda and Furuya [10] identified 17 STGs in the West Kunlun Shan. Chudley and Willis [30] identified nine more STGs. Their findings indicate that the hydrological processes derived from surface meltwater played an important role in sustaining the active surging phase for up to a year. GSs in the Himalayas are isolated and have only been recorded on a glacier scale. Lovell and others [11] found that while at Sabche Glacier in the central Himalayas, surges are likely to be facilitated by subglacial processes, the subglacial topography of the glacier appears to modulate both the timing and duration of each surge [11]. STGs have also been well-documented in the Geladandong Mountain [31,32], Xinqingfeng, Malan Ice Caps [33], and Zangser Kangri Ice Field areas [34]. Recently, terminus advance has been reported for some glaciers in the TJP region, the QLMs; however, a newer list of confirmed surges for the entire area has not been undertaken [35].

In this paper, we use observations of ice surface change to identify STGs in the TJP, representing the first detailed description of STGs in the QLMs. We compare their characteristics to STGs elsewhere in the HMA region and discuss possible mechanisms 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.

2. 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² with an ice volume of ~84.48 km³ 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].

3. Data and Methodology

3.1. 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://earthexplorer.usgs.gov, accessed on 10 September 2021), and Sentinel-2 (https://scihub.copernicus.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. Optical remote sensing images were used in this study.

Date	Sources	Path/Row	Pixel Size (m)	Amount	Application
1973–1977	Landsat MSS	145/033	60 × 60	4	Identification of glacier outlines
1986–2011	Landsat TM	135/033	30 × 30	138	Identification of glacier outlines Visualization of glacier surface morphology Extraction glacier surface flow velocity
1999–2003	Landsat ETM+	135/033	15 × 15	17	Identification of glacier outlines
2013–2020	Landsat OLI	135/033	15 × 15	35	Identification of glacier outlines Extraction glacier surface flow velocity
2000–2020	ASTER	-	15 × 15	30	Identification of glacier outlines Visualization of glacier surface morphology Extraction glacier surface flow velocity
2015–2020	Sentinel-2	-	10 × 10	108	Accuracy assessment of glacier length Visualization of glacier surface morphology Extraction glacier surface flow velocity
2000	SRTM-C		30	1	Extraction of glacier centerlines

). 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.

3.2. 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 sub-pixel 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⁻¹) was acquired using Equation (1):

$$GSV=GSD\times T_y/T_i$$

(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⁻¹ or 3 m d⁻¹ 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]:

$$e_{off}=\sqrt{S{E}^2+ME{D}^2}$$

(2)

where 𝑒_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:

$$SE = STDV/\sqrt{N_{eff}}$$

(3)

In this case, STDV is the standard deviation of the average displacement in the non-glacier region, and N_eff represents the number of effective pixels from which autocorrelation is removed.

$$N_{eff}=N_{total}\times PS/2D$$

(4)

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⁻¹ to 6.63 m a⁻¹ and 0.04 m d⁻¹ to 2.56 m d⁻¹, respectively.

3.3. 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 MMASTER [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. Statistics of vertical errors between the TOPO, SRTM, and ASTER DEM. MED is mean elevation difference, STDV = standard deviation, N = number of considered pixels, SE = standard error, and δ = the overall error of the derived surface elevation change.

Item	MED (m)	STDV (m)	N	SE (m)	δ (m)
SRTM-TOPO DEM	−1.06	8.35	7756	0.09	1.07
ASD02-SRTM	−2.87	6.05	12700	0.05	2.87
ASD07-ASD02	−0.64	6.46	12700	0.06	0.64
ASD12-ASD07	2.66	5.42	12700	0.05	2.66
ASD18-ASD12	2.57	6.22	12700	0.06	2.26

) can be calculated by:

$$\delta =\sqrt{ME{D}^2+{\left(\frac{SD}{\sqrt{N}}\right)}^2}$$

(5)

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.

3.4. 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.

4. Results

4.1. 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², 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. Basic attributes of the STGs in this study. All attributes are obtained from the Randolph Glacier Inventory 6.0.

No.	GLIMS_ID	Area (km2)	Lengthmax (km)	Max/Min/Med Elevation (m)	Mean Slope (°)	Aspect	Debris Cover (%)
NE1	G097793E38464N	14.95	7.34	5582/4411/5016	17.8	NE	1.73
NE2	G097776E38483N	5.89	6.51	5727/4516/5063	14.5	NE	2.39
NE3	G097759E38506N	5.14	6.02	5716/4371/4986	21.3	NE	5.78
NE4	G097729E38514N	6.94	5.48	5750/4557/5064	22.4	NE	1.23
NE5	G097694E38534N	6.10	4.05	5522/4602/5011	20.1	NE	0.56
SW1	G097699E38498N	4.45	5.50	5760/4623/5278	20.7	SW	0.84
SW2	G097731E38477N	14.10	7.85	5755/4557/5297	17.2	SW	1.40

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. Detailed information of each surge-type glacier in this study, including evidence of each surge event, their initiating and terminating year, their duration, their area changed, and their length changed in the active phase. The result was obtained from Landsat, ASTER, and Sentinel-2 images between 1973 and 2020.

No.	Evidence of Surge Events	Surge Initiating Year	Surge Terminating Year	Surge Duration (years)	Area Changed (km2)	Length Changed (m)
NE1	Terminus advanced from 2002 to 2004; Velocities during the active phase typically reach at least an order of magnitude higher than during the passive, quiescent phase.	2002	2006	5	0.64 ± 0.78	718.41 ± 18.35
NE2	Changes in ice crevasses; looped moraines; Clear thickened reservoir area and thinned receiving area.	2000	-	>20	−0.32 ± 0.24	−103.83 ± 10.62
NE3	Terminus advanced; Obvious thinned reservoir area and thickened receiving area during study period.	Before 1986	2000	>15	0.48 ± 0.32	539.88 ± 12.35
NE4	Terminus advanced since 2015; Obvious thinned reservoir area and thickened receiving area.	2013	-	>8	0.05 ± 0.28	98.97 ± 22.35
NE5	Terminus advanced since 2012; Obvious thinned reservoir area and thickened receiving area.	2012	-	>9	−0.03 ± 0.37	114.65 ± 14.21
SW1	Surface features show clear movement; Clear thickened reservoir area and thinned receiving area.	2012	-	>9	−0.08 ± 0.23	−78.14 ± 2.80
SW2	Changes in ice crevasses; Clear thickened reservoir area and thinned receiving area.	2012	-	>9	−0.09 ± 0.49	−304.13 ± 17.60

). 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² (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.

4.2. 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⁷ m³. 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. Mean surface elevation changes of the reservoir area and receiving area and transfer of ice volume for studied glaciers.

No.	GLIMS_ID	Data	Mean ∆H of Reservoir Area (m)	Mean ∆H of Receiving Area (m)	Transfer of Ice Volume (×107 m3)
NE1	G097793E38464N	2002–2007	−22.24 ± 0.64	53.10 ± 0.64	5.71 ± 0.69
NE2	G097776E38483N	2012–2018	−6.41 ± 2.26	1.93 ± 2.26	0.11 ± 0.13
NE3	G097759E38506N	1966–2000	−29.93 ± 1.07	45.07 ± 1.07	2.71 ± 0.06
NE4	G097729E38514N	2012–2018	−13.76 ± 2.26	10.76 ± 2.26	0.34 ± 0.07
NE5	G097694E38534N	2012–2018	−5.57 ± 2.26	5.96 ± 2.26	0.32 ± 0.12
SW1	G097699E38498N	2012–2018	13.44 ± 2.26	−4.38 ± 2.26	-
SW2	G097731E38477N	2012–2018	7.45 ± 2.26	−10.02 ± 2.26	-

). 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⁷ m³ 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).

4.3. Ice Surface Velocities

The annual velocity of the glacier centerline profile provides quantitative information about the surface displacement during the GSs. Figure 3a–g shows the annual surface velocities of the seven STGs during one surge cycle. Note that only glacier NE1 has a complete velocity profile due to the limitations of available satellite imagery. Profile (a) reveals that glacier NE1 peaked between 2004 and 2006 (Figure 3a). The velocity peaks that occurred before 2005 exhibit a fluctuating propagation from upstream to downstream. The maximum value of ~8.0 ± 1.61 m a⁻¹ lies at 1.2 km from the glacier terminus in 2004, then subsequently reached a value of ~17 ± 5.62 m a⁻¹ in 2005 at a distance of 0 km to 1 km from the glacier terminus (Figure 3b). After 2006, the peak velocity decreased markedly over time. According to Figure 3c,e, initiation of both glaciers NE2 and NE4 was slow; their peak velocities eventually reached 40 ± 4.54 m a⁻¹ and 46 ± 4.54 m a⁻¹, 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⁻¹ 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⁻¹) in 2020.

4.4. Glacier Surface Morphology Changes

Glacier surface morphology changes were analyzed using ASTER (15 m), Landsat TM (30 m), and Sentinel-2 (10 m) data. The surfaces of the seven STGs underwent clear morphological changes between 1986 and 2020 in the active phase. As shown in Figure 4a, between 26 April 2002 (toward the beginning of the surge) and 25 April 2004 (midway through the surge), the terminus of glacier NE1 advanced. By 25 April 2007, the crevassing had propagated up-glacier to cover most of the glacier surface, with large compressional crevasses occurring at the glacier terminus and extensional crevasses appearing in the upper glacier area (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.

5. Discussion

5.1. 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. 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⁻¹), 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.

Region	Category	Sample	Area (km2)	Slope (°)	Lengthmax (km)	Max/Min/Med Elevation (m)	Surge Duration (yr)	Period (yr)	Velocitymax (m a−1)	Terminus Advance (km)
Tian Shan	Surge	23	65.59	22.38	12.82	5540/3447/4392	2 to 10	35 to 60	-	0.12 to 3.2
	No-surge	1934	1.59	28.58	1.47	4483/3888/4178
Pamirs	Surge	186	16.84	22.52	7.87	5809/3876/4817	<1 to ≤5	8 to 10	40 to 150	0.2 to 5.88
	No-surge	9597	0.57	27.06	1.08	5100/4650/4868
Karakoram	Surge	258	42.11	25.14	11.19	6425/4379/5384	<1 to >15	8 to 34	300 to 5200	0.11 to 6.1
	No-surge	13576	0.89	30.35	1.19	5665/5156/5417
West Kunlun Mountains	Surge	20	89.41	11.36	18.07	6632/5170/6102	>5	>42	200 to 1200	0.03 to 1.65
	No-surge	528	2.46	23.39	1.87	6206/5684/5976
Puruogangri	Surge	3	22.54	11.67	7.59	6300/5390/5879	-	-	-	-
	No-surge	56	5.89	18.37	3.11	6102/5574/5852
Xinqingfeng and Malan	Surge	8	32.43	13.81	11.07	6485/5031/5685	1 to 7	-	-	0.04 to 1.2
	No-surge	97	3.42	18.19	2.39	5812/5260/5557
Geladandong	Surge	12	24.15	12	9.03	6391/5366/5780	>3 to >19	-	30 to 90	0.17 to 1.8
	No-surge	207	2.68	18.89	2.12	5971/5513/5754
Muztag	Surge	5	55.31	11.04	12.94	6618/5199/5735	4	>40	-	0.55
	No-surge	209	1.84	22.03	1.49	5727/5342/5543
TJP	Surge	7	8.22	19.14	6.11	5687/4519/5102	5 to >20	-	20 to 68	0.11 to 0.72
	No-surge	71	1.09	24.02	1.57	5257/4751/5031

). 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].

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² compared to 2.46 km² 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². 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).

5.2. 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 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.

5.3. 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.

6. 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.