Surging Dynamics of Glaciers in the Hunza Valley under an Equilibrium Mass State since 1990

: Previous studies have shown that glacier changes were heterogeneous in the western Karakoram, with the coexistence of retreating, advancing, and surging glaciers. However, it remains unclear that the mechanisms driving these changes. Based on the Shuttle Radar Topography Mission (SRTM) DEM and TerraSAR-X / TanDEM-X images (2014), this study presents glacier surface height changes in the Hunza Basin of the western Karakoram, employing the method of di ﬀ erential synthetic aperture radar interferometry (DInSAR). A slight negative glacier mass balance was observed in the Hunza Basin during 2000–2014. Surge-type glaciers would not have an obvious e ﬀ ect on overall mass balance in regional assessments over long-time scales. Further, glacier surface velocities in the Hunza Basin were estimated from Landsat images for the period of 1990–2018 by utilizing published data sets and Landsat images. Compared to the annual glacier surface velocities, 22 surge events were observed in seven surge-type glaciers in the Hunza Basin. Glacier ﬂow can be attributed to thermally and hydrologically control, and the geomorphological characteristics of di ﬀ erent individuals. This study gives us a new insight into the situation of the “Karakoram anomaly” under the background of glacier mass loss in the high mountains of Asia (HMA).


Introduction
Glaciers in the high mountains of Asia (HMA) are the sources of many major rivers and lakes in Asia [1]. As sensitive climate indicators, continuous long-term monitoring on glacier is increasingly important for the interpretation of global changes [2]. Because of global warming in recent decades, the mass loss in the HMA from 2003 to 2009 accounts for~10% to the global glacier mass loss, which contribute 3% to the sea level rise [3], and debris coverage in glacier increased [4]. Meanwhile, the frequency of glacial hazards increased, such as glacial lake outburst floods (GLOF) and glacier surge [5][6][7]. Despite in-situ observations being encouraged, the great labor cost and rugged terrain has limited its application in the HMA. Based on geodetic methods, previous studies have shown that glacier mass balance was heterogeneous in the HMA [8][9][10]. A significant impact can be made by local parameters on glacier responses to climate change; therefore, there may have a heterogeneous pattern of glacier changes [11].
In contrast to the Nyainqentanglha, Tien Shan, and Himalayas, which were experiencing pronounced negative glacier mass balances [8,12], a complex pattern of glacier changes had been evident for the Karakoram, the Western Kunlun, and the Pamirs in recent decades [2,8,13]. Due to the occurrence of surges and complex glacier behaviors [14], these results of mass gain are difficult to interpret. Furthermore, many glaciers in the western Karakoram are covered by debris [13]. Influenced by supraglacial debris, glacier termini can reach to lower elevations and react differently to warming [13,15]. However, it remains unclear whether the debris cover has any impact on glacier surges or other complex behaviors. Hence, the details of the "Karakoram anomaly" remain controversial and deserves more attention.
Glacier surface height changes and velocities are particularly crucial parameters that determine glacier discharge. Measurements of these parameters could help us to better understand the glacier responses to climate change [16,17]. Predecessors in glacier surface height changes and velocities in western Karakoram have already done a lot of significant research work and achieved fruitful results [8][9][10]13,[18][19][20][21], and the relation between the two crucial parameters is fairly precise [22]. For complex patterns of glacier changes in western Karakoram, several mechanisms are proposed, while there are knowledge gaps and uncertainty [22]. In our study, the decadal glacier mass balance in the Hunza Basin of the western Karakoram was estimated from the differential synthetic aperture radar interferometry (DInSAR) method on the Shuttle Radar Topography Mission (SRTM) DEM and TerraSAR-X/TanDEM-X (TSX/TDX) images (2014). Additionally, glacier surface velocities in the Hunza Basin were derived from Landsat images for the period of 1990-2013 by utilizing cross-correlation feature tracking techniques. Based on the analysis of mass balance and annual variability of velocity, the spatiotemporal variabilities of glacier were investigated in the western Karakoram.

Study Area
The Hunza Basin (35 • 42 -37 • 06 N, 72 • 30 -75 • 47 E), located in the western Karakoram, is one of the glaciated sub-catchment of the Upper Indus Basin (Figure 1). The typical characteristics of the study area are steep valley walls and narrow valleys [23]. The Hunza River is contributed by more than 14 small and medium sized tributaries, such as Khynjerab, Hispar, Hassanabad, Mulungutti, and Verjerab.
The Hunza Basin, located in subtropical climate zone, is affected by the Westerlies and the Indian Summer Monsoon [24]. Mean annual temperature ranges from 2.8 to 6.5 • C in the valleys at the elevation of 2810-3669 m a.s.l. [24], and mean daily temperature ranges from −9.0 • C in January to 14.6 • C in July based on observations at Ziarat (3669 m a.s.l.) [25]. Due to the limitation of high mountains, the influence of monsoon weakens north-westward [26]. While the Westerlies brings intensity precipitation in both summer and winter, about 2/3 of high-altitude snowfall in the study area result from the Westerlies [27]. The snow cover in winter is about 80% and decreases to 30% in summer [28].
Based on the Randolph Glacier Inventory, 2209 glaciers developed in the Hunza Basin, with glacier area of 4403.5 km 2 [29]. Due to great vertical height difference and abundant precipitation, many large valley glaciers developed in the Hunza Basin. The equilibrium-line altitude (ELA) in study area ranges from 4500 m to 5500 m a.s.l. [9,11,30]. The ablation areas of most glaciers are covered by debris, and about 90% of glacier area lies in the 5000-6000 m elevation range [31]. Using CORONA and Landsat imagery, glacier area in the Hunza Basin decreased by 1.36% from 1973-2014 [24]. Many glaciers surged in recent years, such as Hispar Glacier, Khurdopin Glacier, and Shishper Glacier [14,32,33]. Figure 1. Regional overview of the study area in the Hunza Basin, western Karakoram. Background: Landsat Operational Land Imager (OLI) image (Bands 6, 5, and 3) mosaicked from LC81490342013202LGN00, LC81490352013202LGN00, LC81500342014260LGN01, and LC81500352013209LGN00. Light blue areas represent bare ice, red-brown areas represent bare soil and rocks, green areas represent grass or other vegetation, and yellow areas represent debris-covered glaciers.
The Hunza Basin, located in subtropical climate zone, is affected by the Westerlies and the Indian Summer Monsoon [24]. Mean annual temperature ranges from 2.8 to 6.5 °C in the valleys at the elevation of 2810-3669 m a.s.l. [24], and mean daily temperature ranges from −9.0 °C in January to 14.6 °C in July based on observations at Ziarat (3669 m a.s.l.) [25]. Due to the limitation of high mountains, the influence of monsoon weakens north-westward [26]. While the Westerlies brings intensity precipitation in both summer and winter, about 2/3 of high-altitude snowfall in the study area result from the Westerlies [27]. The snow cover in winter is about 80% and decreases to 30% in summer [28].
Based on the Randolph Glacier Inventory, 2209 glaciers developed in the Hunza Basin, with glacier area of 4403.5 km 2 [29]. Due to great vertical height difference and abundant precipitation, many large valley glaciers developed in the Hunza Basin. The equilibrium-line altitude (ELA) in study area ranges from 4500 m to 5500 m a.s.l. [9,11,30]. The ablation areas of most glaciers are covered by debris, and about 90% of glacier area lies in the 5000-6000 m elevation range [31]. Using CORONA and Landsat imagery, glacier area in the Hunza Basin decreased by 1.36% from 1973-2014 [24]. Many glaciers surged in recent years, such as Hispar Glacier, Khurdopin Glacier, and Shishper Glacier [14,32,33]. Figure 1. Regional overview of the study area in the Hunza Basin, western Karakoram. Background: Landsat Operational Land Imager (OLI) image (Bands 6, 5, and 3) mosaicked from LC81490342013202LGN00, LC81490352013202LGN00, LC81500342014260LGN01, and LC81500352013209LGN00. Light blue areas represent bare ice, red-brown areas represent bare soil and rocks, green areas represent grass or other vegetation, and yellow areas represent debris-covered glaciers.

Glacier Outlines
Glacier outlines could be acquired from the Randolph Glacier Inventory [29]. For the Hunza Basin, these glacier outlines acquired from Landsat images during the period 1998 -2002. To analyze the changes of surface elevation and velocities, glacier extents in different periods should be considered, especially in the Karakoram region, where more than 50% of glaciers advanced. Hence, glacier mapping in 2014 was delineated in this study.
A semi-automated approach was employed to acquire glacier outlines in 2014 from Landsat OLI images. This approach was applied at the cloud-computing platform Google Earth Engine (GEE) [34]. First, a pixel-based approach was used to select the suitable images [34]. To minimize the effect of snow cover, all images during January 2013-December 2015 were selected from the ablation season from July to October. Using a cloud score algorithm, pixels of cloud cover in the images were masked [35]. Next, to filter the remaining clouds and shadows, the 30th to 70th percentile range of the entire image stack were used to calculate a mean value for each pixel [34]. Then, ice/snow body was identified from resulting image by the normalized difference snow index (NDSI) [36]. The NDSI [(λ green − λ SWIR1 )/(λ green + λ SWIR1 )] includes Landsat OLI band 3 (green) and band 6 (SWIR1). Based on the NDSI image and threshold value, ice/snow body was distinguished from other land cover types. For isolated ice/snow patches smaller than 0.01 km 2 , they were eliminated by a 3 × 3 median filter [37].
Finally, to discriminate clean ice and seasonal snow, the derived glacier outlines were checked manually. The distribution of debris cover is an important parameter in the western Karakoram; therefore, the same method (NDSI from Landsat images in the GEE platform) but different threshold value (NDSI < 0.4) were used to derive debris-covered ice surfaces.
Comparing delineated outlines with ground truth or reference data derived from high-resolution air photos is a good way to assess the uncertainty of glacier mapping [38]. A previous study indicated that the average offsets of glacier outlines acquired from Landsat image were ±10 m for clean ice and ±30 m for debris-covered areas, respectively [39]. Based on these results, an uncertainty of ±4.6% were calculated for glacier outlines in 2014.

Glacier Height Changes
Glacier height changes (∆H) were calculated from seven TSX/TDX acquisitions and the SRTM DEM. The SRTM C-band DEM with one arc second resolution in WGS84/EGM96 was used in this study. Flying in close orbit formation, the TerraSAR-X and TanDEM-X act as a flexible single-pass SAR interferometer [40]. Acquired in bistatic InSAR stripmap mode, seven pairs of X-band bistatic TSX/TDX acquisitions were used to derive glacier height change ( Figure 1, Tables 1 and 2). In the ground range and azimuth direction, the resolutions of TSX/TDX image were both~2.5 m.  Based on the method of DInSAR, the TSX/TDX acquisitions were processed by interferometric processing software [12]. Glacier height changes can be acquired from bistatic interferograms, while topographic phases and flat earth, contained in bistatic interferograms, should be removed first. The flat earth and topographic phases were simulated by orbital information from SRTM DEM and TSX/TDX acquisitions, and then removed from the bistatic interferogram and leave the topographic residual phase, which called differential interferogram. The topographic residual phase, mainly caused by the change of glacier height, can be converted to height changes directly [41][42][43].
Using an intensity cross-correlation algorithm, a refinement of the lookup table was acquired from a horizontal co-registration between the SRTM DEM and the TSX/TDX acquisition [44]. The refinement of the lookup table and precise orbit information of the TSX/TDX scene were used to eliminate the horizontal offsets between the TSX/TDX acquisition and SRTM DEM. Using GAMMA's minimum cost flow algorithm, differential interferogram was unwrapped and converted to absolute differential heights [45].
Because of the poor quality of coherence for topographic measurement, foreshortening, layover, and shadow regions of the TSX/TDX interferogram cannot be used. The unusable regions were masked out. Assuming that elevation changes is similar at the same altitude [9], the changes of glacier height were calculated in each 100 m altitude bands [46], and "no data" pixels were replaced with the mean ∆H of the corresponding altitude band [47].
The penetration depth of microwaves into the snow, firn, and ice can range from 0 m to 10 m [48]. Hence, glacier surface elevation derived from InSAR might be lower than the actual glacier surface. Due to the advantage of high resolution, large coverage, and freely availability, the SRTM DEMs are being used in many studies; therefore, the penetration depth of SRTM is widely discussed [9,10,12,41,46]. The published value from Lin et al. [46]-a penetration depth of 2.41 ± 0.17 m for the Karakoram region-was applied in this study. Considered that the observed time of TSX/TDX were closed to the performed time of SRTM, and the carrier frequencies of the SRTM X-band and TSX/TDX are almost the same, Lin et al. [46] presumed that no penetration difference exists between the SRTM X-band DEM and TSX/TDX. Hence, the elevation difference between SRTM X-band and C-band DEM can be regarded as the penetration depth when the changes of glacier height calculated from TSX/TDX and SRTM C-band DEM.
The changes of glacier height were converted into mass balance by the ice/firn/snow density of 850 ± 60 kg/m 3 [49]. Assuming that heights remains unchanged in non-glacier areas for a long time, the mean elevation differences (MED) and standard deviation (SD) in non-glacier areas can be used to estimate the errors of the changes of glacier heights. In order to reduce the error of autocorrelation, a distance of 500 m was chosen for difference map derived from SRTM DEM and TSX/TDX. Hence, the overall errors of the changes of glacier height can be calculated by where N is the pixel number with effective measurements. Overall, the error of mass balance was calculated by the errors of glacier area, height changes, and ice density uncertainty.

Glacier Surface Velocities
Based on the auto-RIFT feature tracking, the Inter-Mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) project extracted annual mean surface velocities for major glacier-covered regions from Landsat imageries [50,51]. Due to the low quality of data sources, the ITS_LIVE data product has large uncertainty in the earlier product years. Therefore, glacier surface velocities in the Hunza Basin from 1990-2018 were used in this study.
In addition, seasonal surface velocities in the period of 2013-2019 were estimated from Landsat OLI images using cross-correlation feature tracking techniques. This application was processed on Co-Registration of Optically Sensed Images and Correlation (COSI-Corr), an IDL module for the remote sensing platform ENVI [52]. Previous studies have proved that cross-correlation feature tracking on COSI-Corr is an effective method for measuring glacier surface velocity [53][54][55].
First, the 15-m panchromatic band of the Landsat OLI images was co-registered to a common reference dataset. Then, horizontal ground displacements were estimated from sub-pixel correlation of co-registered images. With search windows of 64 × 64 pixels (pattern size) and 128 × 128 pixels (search area), image correlation was acquired iteratively [20,55]. The east-west and north-south displacements and signal-to-noise ratio (SNR) were yielded by image correlation. Glacier velocity were obtained by the combination of horizontal ground displacements. In order to improve the accuracy of surface displacement, image matches that SNR < 0.9 and extreme values were removed [20]. Finally, the surface displacement converts into annual glacier surface velocity by using the following equation: where GSV is the glacier surface velocity, GSD is the glacier surface displacement, T i is the time interval, and T y is a constant of 365 or 366.
As with the glacier height changes, the uncertainties of glacier velocity were estimated by the residual value in non-glacier areas. Gaussian low path filtering was used to smooth the velocity map in order to minimize the effect of autocorrelation [54]. Thus, image pairs separated by~365 days were used to measure annual glacier surface velocity. The uncertainties of annual and seasonal glacier surface velocity depend on the time separation between the images, ranging from 2-20 m a −1 and 0.02-1.00 m d −1 , respectively.

Mass Balance
The average surface thinning of glacier in the Hunza Basin was −0.61 ± 0.07 m from 2000 to 2014 ( Figure 2a and Table 3). Glaciers with an area of 2878.2 km 2 experienced a mean mass balance of −0.04 ± 0.02 m w.e. a −1 . Overall, mass balance increased with altitude, and it varied significantly in different elevation ranges ( Figure 3). Glaciers experienced a negative mass balance of −0.21 ± 0.02 m w.e. a −1 below 5100 m a.s.l., while it turned to a positive mass balance of 0.11 ± 0.02 m w.e. a −1 above 5100 m a.s.l. during 2000-2014. Even though glaciers experienced a negative mass balance below 5100 m a.s.l., a significant thickening was observed in the 3800-4100 m a.s.l.; the main reason for this phenomenon is glacier surging in Hispar Glacier.
The ablation area of many glaciers in the Hunza Basin are covered by debris. The total debris-covered area of 587.8 km 2 account for~13.3% of the total glacier area. Approximately 98.3% of the debris-cover is found in the 2700-5200 m a.s.l. range. Compared with clean ice in the same elevation range in all non-surging glaciers, negative mass balance was greater on the debris-covered area during 2000-2014 (−0.46 ± 0.15 m w.e. a −1 on debris-covered area vs. −0.08 ± 0.02 m w.e. a −1 on clean ice). This indicated that the debris cover has an acceleration effect on surface melting overall, whereas the debris cover has different impact in different elevation ranges (Figure 3).
Glaciers in the study area experienced retreating, advancing, surging, or quiescent phases. After surging, heterogeneous mass balance was observed in these glaciers during 2000-2014. According to the classification of surge and non-surge glaciers, surge glaciers experienced a negative mass balance of −0.06 ± 0.02 m w.e. a −1 , it was larger than the negative mass balance of non-surg glaciers (−0.03 ± 0.02 m w.e. a −1 ). All the non-surge glaciers experienced a negative mass balance, while mass balance of surge-type glaciers was heterogeneous (Table 3)    The ablation area of many glaciers in the Hunza Basin are covered by debris. The total debriscovered area of 587.8 km 2 account for ~13.3% of the total glacier area. Approximately 98.3% of the debris-cover is found in the 2700-5200 m a.s.l. range. Compared with clean ice in the same elevation range in all non-surging glaciers, negative mass balance was greater on the debris-covered area during 2000-2014 (−0.46 ± 0.15 m w.e. a −1 on debris-covered area vs. −0.08 ± 0.02 m w.e. a −1 on clean ice). This indicated that the debris cover has an acceleration effect on surface melting overall, whereas the debris cover has different impact in different elevation ranges (Figure 3).
Glaciers in the study area experienced retreating, advancing, surging, or quiescent phases. After

Surface Velocity
Through visual inspection of the successive satellite imagery and examining the spatiotemporal evolution of surface velocities of each glacier, 16 glaciers in the Hunza Basin were classified into surge Remote Sens. 2020, 12, 2922 9 of 21 and non-surge types. Their characteristics, including area, debris cover, elevation, and the length of centerline, are listed in Table 4.  Seven glaciers were identified as surge-type glaciers, with 22 discrete surge events during the investigated period (Table 5). All surge-type glaciers experienced surge events more than once over the study period.  Seven glaciers were identified as surge-type glaciers, with 22 discrete surge events during the investigated period (Table 5). All surge-type glaciers experienced surge events more than once over the study period.   For Yazgil Glacier, four surge events were captured during 1990-2018 ( Figure 5). From the terminus to accumulation area, there have two peak velocities in profiles of Yazgil (13-15 km and 19-22 km from the terminus). The first peak velocity was induced by the increase of slope, and there were no significant fluctuation in surface velocities above the altitude of first peak velocity over the study period. Below the altitude of first peak velocity, profiles of Yazgil revealed four surge events between 2 km and 11 km from the terminus, with 1997-1998, 2005-2006, 2010-2011, and 2018, respectively. These surge events lasted~2 years, and the peak velocities of surge events decreased gradually from 1990-2018.
Profiles of the two branches of Khurdopin Glacier revealed that two surge events occurred in the east branch, and four surge events occurred in the west branch ( Figure 5). In the east branch of Khurdopin, the two surge events lasted~3 years in the periods of 1998-2000 and 2016-2018 ( Figure 6). The interval of the surge events is~16 years, if it be taken as a surge cycle period, the east branch of Khurdopin could be about to enter a long-term quiescent phase. The four surge events in the west branch lasted~2 years, with~7 years interval. The latest surge event began in 2018; this indicated that the west branch is entering a short-term active surge phase. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 23   For Yazgil Glacier, four surge events were captured during 1990-2018 ( Figure 5). From the terminus to accumulation area, there have two peak velocities in profiles of Yazgil (13-15 km and 19-22 km from the terminus). The first peak velocity was induced by the increase of slope, and there were no significant fluctuation in surface velocities above the altitude of first peak velocity over the study period. Below the altitude of first peak velocity, profiles of Yazgil revealed four surge events between 2 km and 11 km from the terminus, with 1997-1998, 2005-2006, 2010-2011, and 2018,

Glacier Mass Balance
Previous studies revealed both negative and positive glacier mass balances in the study area. Gardelle et al. [9] showed a slight mass gain of +0.09 ± 0.18 m w.e. a −1 using the geodetic method based on C-band SRTM and SPOT5 DEMs between 2000 and 2008. However, Kääb et al. [30] and Gardner et al. [3] observed slight mass loss in the Karakoram, with glacier mass changes of −0.06 ± 0.04 m and −0.10 ± 0.14 m w.e. a −1 derived from SRTM and ICESat altimetry measurements during 2003-2009, respectively. Lin et al. [46] recorded a near-stable mass balance of −0.020 ± 0.064 m w.e. a −1 during 2000-2014 in the western Karakoram. Brun et al. [8] presented a mean mass loss of −0.03 ± 0.14 m w.e. a −1 from 2000-2016 in Karakoram. In this study, a slight mass loss of −0.04 ± 0.02 m w.e. a −1 was measured in the Hunza Basin. Compared to previous studies, the result of Gardelle et al. [9] were slightly different than those of other studies. Given that large uncertainty existed in previous studies, it is concluded that glaciers in west Karakoram experienced a near-stable mass balance since 2000.
In the Hunza Basin, debris cover has an accelerating effect on surface melting overall, while it plays different roles in different elevation ranges, with insulating and accelerating effects below and above 3500 m a.s.l., respectively. Except for debris cover, supraglacial/proglacial lakes and ice cliffs may also affect surface melting (Figure 7). Many crevasses and supraglacial lakes were developed in the ablation area, and the surface meltwater and liquid precipitation can flow down into the glaciers by crevasses. The englacial melt funnels/holes collapsed due to the coexistence of thick debris and inner ablation, and then ice cliffs and supraglacial lakes formed easily. The very fine debris attached to ice cliff, can reduce the ice albedo and absorb more shortwave radiations [56]. Hence, because of the effect of heat conduction, supraglacial lakes and ice cliffs aggravated the glacier surface melting. For surge glaciers, large volumes of ice were stored at high elevation for many decades before rapidly discharging the mass down-glacier during surge events. In short time scales, the surface elevation of the ablation area thickened significantly. A significant thickening was observed in the 3800-4100 m, which was the main reason for this is glacier surging in the north branch of Hispar between 2009-2011 ( Figure 7). However, that ice came from high elevation and experienced a significant thinning at low elevation over long time scales. Therefore, surge-type glaciers in the Hunza Basin experienced a greater mass loss than non-surge-type glaciers. It is concluded that even though some surge-type glaciers may have fluctuations in mass balance over short-time scales, their inclusion will not have a significant influence on overall mass balance in regional assessments.

Glacier Surface Velocity
Previous studies have shown that glacier surges, controlled by hydrology, tend to decelerate For surge glaciers, large volumes of ice were stored at high elevation for many decades before rapidly discharging the mass down-glacier during surge events. In short time scales, the surface elevation of the ablation area thickened significantly. A significant thickening was observed in the 3800-4100 m, which was the main reason for this is glacier surging in the north branch of Hispar between 2009-2011 ( Figure 7). However, that ice came from high elevation and experienced a significant thinning at low elevation over long time scales. Therefore, surge-type glaciers in the Hunza Basin experienced a greater mass loss than non-surge-type glaciers. It is concluded that even though some surge-type glaciers may have fluctuations in mass balance over short-time scales, their inclusion will not have a significant influence on overall mass balance in regional assessments.

Glacier Surface Velocity
Previous studies have shown that glacier surges, controlled by hydrology, tend to decelerate during summer and accelerate during winter [57,58]. This indicated that a seasonal signal could be exist to hydrologically controlled surge front propagation, whereas thermally controlled glacier surges can initiate at any time because these surges rely on the conditions of the glacier bed. Although there is no enough data to analyze the thermal regime of surge events, cold-based ice is considered to predominate at lower and higher elevations, thick warm-based ice predominates at middle elevations [14], and this has been borne out in Aru Glacier [5]. Meanwhile, profiles of several surge glaciers (Batura, the east branch of Khurdopin, Virjerab) revealed that surface velocities accelerated gradually before surging ( Figure 5). This phenomenon meets the characteristic of thermally controlled surges. In order to detect if there has been a seasonal signal in surge events, seasonal variabilities of surface velocities were estimated from Landsat imageries with a high temporal resolution (  (Figure 9a). This suggests that these flow instabilities in Hispar can be attributed to a thermally controlled surge. However, surge events in Yazgil and the west branch of Khurdopin are similar to the surge event of Kelayayilake Glacier in Pamir, which was initiated in winter and terminated in the following summer [6]. It is reported that water discharge reduced significantly before and during the surge; therefore, the surge event of Kelayayilake can be regarded as a hydrologically controlled surge. Even though there have no hydrological data in the Hunza Basin to analyze the processes of water discharge, profiles of Yazgil revealed that surface velocities have a significant seasonal fluctuation during 2013-2019 ( Figure 9b). Apparently, surge events of Yazgil are partially hydrologically controlled.
Except for the thermal and hydrological conditions, the geomorphological characteristics of glaciers, such as surface features, the size of the accumulation area, the width of the outlet, an surface slope, have a significant impact on glacier surge [5,7,20,54]. Seven surge glaciers with 22 surge events were identified in the Hunza Basin over the study period, and all surge events occurred in the middle parts of glaciers. Due to the obstruction of debris cover and long glacier tongues, surging ice were inhibited to reach the terminus; therefore, the termini of these glaciers remained stable after the surging [54]. New crevasses appeared, surface features, i.e., ice cliffs and supraglacial lakes, were rearranged, and the remnant of a previous glacier surge was overridden; all of these resulted from the spatial difference of surface velocities [7,54]. In turn, ice cliffs and supraglacial lakes aggravate the surface melting and more meltwater flowed down into glacier via crevasses. Then, sufficient hydrological conditions can prompt glacier surge. Besides, the main branch of Hispar has a broad outlet, which make ice transport to downstream stable, resulting in a longer surge cycle period. While Yazgil Glacier has a narrow outlet and a large accumulation area, ice/snow accumulates quickly and is frequently transported downstream, resulting in a shorter surge cycle period. Although cycle periods of glacier surge vary widely between individual glaciers, in general, the narrower the outlet and the larger of the accumulation area, the shorter of the surge cycle period [54].
According to the surge cycle period during 1990-2018, some glaciers are about to enter a quiescent phase, whereas some glaciers will enter an active phase. The first surge event of Virjerab lasted~4 years from 2004-2007, and the second surge event began in 2018. The characteristics of profiles of surface velocities at the same position in the years before a surge were similar to each other, which indicated that it is approaching its next surge. For Hispar, Batura, and Khurdopin, longer surge cycle periods resulted from a broad outlet or a glacier surge having just terminated so that there has not much ice that could be transported downstream. This indicates that those glaciers are about to enter a quiescent phase. As mentioned above, flow instabilities can be attributed to thermal and hydrological control and the geomorphological characteristics of different individuals. These factors are influenced by climate change. In the context of global warming, the thermal and hydrological conditions change accordingly. For example, four surge events were captured in Yazgil during 1990-2018; with a relatively shorter surge cycle period, surges will occur more frequently in the future according to the cycle period. However, the peak velocities of surge events in Yazgil decreased gradually from 1990 to 2018. Due to atmospheric warming, Yazgil experienced a significant mass loss during 2000-2014; this could weaken the thermal and hydrological conditions, resulting in a gradual decrease of the peak velocities of surge events.   Although some other surge glaciers remain undiscovered, and such surge-type glaciers in the Hunza Basin represent a small part of all glaciers in HMA, the investigation of fast glacier flow, glacier surge and processes have a great contribution to glacier research [59]. Many surge events, rapid acceleration and deceleration, and decreased peak velocities of surge events give us a new insight into the "Karakoram anomaly" in the context of glacier mass loss in HMA.  Although some other surge glaciers remain undiscovered, and such surge-type glaciers in the Hunza Basin represent a small part of all glaciers in HMA, the investigation of fast glacier flow, glacier surge and processes have a great contribution to glacier research [59]. Many surge events, rapid acceleration and deceleration, and decreased peak velocities of surge events give us a new insight into the "Karakoram anomaly" in the context of glacier mass loss in HMA.

The Relationship between Mass Balance and Surface Velocity
All the non-surge-type glaciers in the Hunza Basin experienced a negative mass balance during 2000-2014. Except for Minapin Glacier, there have been no significant fluctuation in surface velocities of non-surge glaciers over three decades. Negative mass balance weakens the thermal condition but produces more meltwater, which can flow down into glaciers via crevasses and moulins and enhance the hydrological conditions. Therefore, glacier mass loss would not result in the decrease of surface velocity significantly over short time scales. For surge glaciers, there has been a rapid acceleration and deceleration in surface velocities, while mass balance was heterogeneous over the study period. Overall, mass balance increased with altitude, and glaciers experienced a significant thinning in the 2700-5100 m a.s.l. range due to the accelerating effect of debris cover, while it turned to mass gain above 5100 m a.s.l. during 2000-2014. This can cause an increase of the glacier gradient, resulting in an increase of driving stresses and ice velocities [60,61].
Our data indicated that even though some surge-type glaciers may have fluctuations in mass balance over short time scales, their inclusion will not have a significant influence on overall mass balance in regional assessments. Glaciers experienced a mass gain above 5100 m a.s.l., and then more ice was transported from up-glacier to terminus due to the rapid acceleration. Mass loss at higher elevation resulted from glacier surge; meanwhile, climate warming caused a quick mass loss at lower elevation. This suggests that glacier surge could result in a greater mass loss compared to non-surge-type glaciers over long time scales [8,9,21].

Climatic Considerations
The Hunza Basin climate is characterized by the westerlies, and the glaciers have experienced a near-stable mass balance since 2000. Previous studies indicate that this balanced pattern of mass balance is related to the changes in precipitation [9]. The Global Precipitation Climatology Project (GPCP) data set shows precipitation increasing over the Pamir and the Karakoram since 1970 [2]. In addition, the mean summer temperature in the Karakoram has been decreasing since 1961 [62]. These precipitation increases in winter and temperature decreases in summer result in greater accumulation and reduced ablation. Due to the effect of glacier surging, large volumes of ice rapidly discharged the mass down-glacier during surge events and experienced a significant thinning at low elevation over long time scales. Therefore, a nearly balanced glacier mass budget was present in the Hunza Basin.
In order to analyze the relationship between glacier surge and extreme precipitation, Climatic Research Unit TS V4.03 (CRU V4.03) data was employed in this study. The CRU V4.03 data is the gridded time-series dataset from 1901 to 2018 [63], which was constructed using the Climate Anomaly Method with observations of the whole world. The dataset provides monthly total precipitation with a spatial resolution of 0.5 • ×0.5 • . The change of winter precipitation derived from this data is shown in Figure 10. It is clear that winter precipitation increased slightly in the Hunza Basin during 1991-2018. This could provide an abundance of matter sources sufficient power for glacier surging. In addition, heavy precipitation in winter occurred in 1998,2005,2009, 2013, and 2017; meanwhile, several glacier surge events were initiated in the same year. Furthermore, heavy precipitation has become more frequent in recent decade, and these changes are consistent with the more frequent glacier surges. The IPCC AR5 suggests that the occurrence of heavy precipitation events might increase in the future [64], and this might induce more frequent glacier surges.

Conclusions
In this study, the differential synthetic aperture radar interferometry (DInSAR) method on the SRTM DEM and TerraSAR-X/TanDEM-X were employed to calculate glacier mass balance in the Hunza Basin of the western Karakoram. The temporal and spatial evolution of glacial motion in the Hunza Basin were acquired from published data sets and Landsat images.
Glaciers in the Hunza Basin experienced a slight negative mass balance during 2000-2014. The debris cover has an acceleration effect on surface melting overall, whereas the debris cover has different impacts at different elevation ranges. Even though some surge-type glaciers may have fluctuations in mass balance over short time scales, their inclusion will not have a significant influence on overall mass balance in regional assessments over long time scales.
By analyzing and quantifying glacial motion from 1990 to 2018, seven glaciers were identified as surge-type glaciers in the Hunza Basin, with 22 discrete surge events from the studied glaciers. Flow instabilities can be attributed to thermal and hydrological control and the geomorphological characteristics of different individuals.

Conclusions
In this study, the differential synthetic aperture radar interferometry (DInSAR) method on the SRTM DEM and TerraSAR-X/TanDEM-X were employed to calculate glacier mass balance in the Hunza Basin of the western Karakoram. The temporal and spatial evolution of glacial motion in the Hunza Basin were acquired from published data sets and Landsat images.
Glaciers in the Hunza Basin experienced a slight negative mass balance during 2000-2014. The debris cover has an acceleration effect on surface melting overall, whereas the debris cover has different impacts at different elevation ranges. Even though some surge-type glaciers may have fluctuations in mass balance over short time scales, their inclusion will not have a significant influence on overall mass balance in regional assessments over long time scales.
By analyzing and quantifying glacial motion from 1990 to 2018, seven glaciers were identified as surge-type glaciers in the Hunza Basin, with 22 discrete surge events from the studied glaciers. Flow instabilities can be attributed to thermal and hydrological control and the geomorphological characteristics of different individuals.