Surface Velocity Analysis of Surge Region of Karayaylak Glacier from 2014 to 2020 in the Pamir Plateau

The west branch of Karayaylak Glacier (eastern Pamir Plateau) surged in May 2015, significantly impacting on local socio-economic development. This event was also of great significance for studies of surging glaciers. Using Sentinel-1 imagery analyzed by offset tracking, based on normalized cross-correlation (NCC), and with the support of the Google Earth Engine (GEE) platform, we quantified the ice surface velocity of the west branch and terminus of Karayaylak Glacier from 13 October 2014 to 17 October 2020. Sentinel-1 images were acquired at intervals of 12 or 24 days. We also used a three-dimensional (3-D) laser scanner to measure the velocity of 3 ablation stakes and 56 feature points in the study region from 15 August to 6 October 2015, for the purpose of accuracy assessment. We set up an automatic meteorological station to record the air temperature in the same period and combined this with data from Tashkurgan meteorological station from 1957 to 2015. Analysis of this dataset provided insights into the glacier surge mechanism, with the following conclusions. (1) Surface velocity of the west branch and terminus of Karayaylak Glacier increased sharply after October 2014. The velocity then dropped significantly in the two months after the surge, and stayed at low values for nearly a year. After 2017, the velocity was slightly higher than in the previous period. (2) The surge event occurred from 11 April to 17 May 2015; the average surface velocity in this phase attained 2395 m a−1 with a maximum velocity of 4265 m a−1 at the west branch terminus. (3) From 2017 to 2020, the velocity showed periodic annual changes. (4) Based on the meteorological data analysis, we conclude that this surge resulted from the interaction between thermal and hydrological control mechanisms. Simultaneously, we demonstrate the high potential of the GEE platform and Sentinel-1 data to extract glacier surface velocity.


Introduction
The Tibetan Plateau and surrounding mountains represent the world's largest glacierized area outside of Alaska and the Arctic, Antarctic, and Greenland [1,2]. Research on the Pamir-Karakoram-Himalaya mountain glaciers has included studies of glacier length and area (the latter exceeds 70,000 km 2 [3]), debris cover [4,5], mass balance [6][7][8][9] and glacier impact on climate change and sea-level [1,10]. This region is also regarded as one of the most active surge zones [11,12]. Glacier surges are usually defined as quasi-periodic advances or increases in flow velocity unrelated to external triggers [13,14]. Surging glaciers continuous time series, due to cloudy conditions, especially in high mountain and cold regions. Radar remote-sensing technology can penetrate clouds and is, therefore, useful in all weather conditions. It has also been proven advantageous in rugged mountain areas with complex topography [52,53].
Offset-tracking procedures are usually adopted in deformation monitoring with longer time intervals, or glacier analysis with rapid changes [54][55][56]. Strozzi implemented two complementary offset-tracking methods in 2002: intensity and coherence tracking, and compared their respective results to DInSAR, finding offset tracking can be an alternative when the loss of coherence limits DInSAR [53]. In the NCC-based offset-tracking algorithm for displacement measurement, the search window size is a key factor affecting the magnitude of displacement. Huang et al. proposed the average velocity gradient (AVG) method to optimize window size for the locally adaptive flow field [57]. Many research works have been devoted to glacier surface velocity monitoring using the offset-tracking method [58,59]. Wendt [25]. The above research is all based on early SAR images that have short time series which do not continue beyond the surge phase. Ground-based measurements of glacier surface velocity are also absent from these studies, which means there is no guarantee of accuracy.
Sentinel-1/2 images are available from 2014 onwards and have become well suited to glacier research in recent years because of their high spatio-temporal resolution [61,62]. We also note that the studies noted above used Windows-based software such as GAMMA and COSI-Corr, which are time-consuming and labor-intensive for large SAR data processing tasks, thus restricting the acquisition of long-term, continuous time series [63]. The emergence of cloud data processing platforms in recent years has largely solved this issue: for example, the Pleiades supercomputer provided by National Aeronautics and Space Administration (NASA) Earth Exchange (NEX) [64], Amazon Web Services (AWS) [65] and Google Earth Engine (GEE). GEE combines a multi-petabyte catalog of satellite imagery and geospatial datasets with planetary-scale analysis capabilities [63,66], and has been widely applied to the remote-sensing [67][68][69][70] (including glacier remote sensing [70]) domain due to its powerful functions and convenient acquisition.
Given the advantages of the GEE platform, this study derived a continuous surface velocity time series of the surge region of Karayaylak Glacier, from 2014 to 2020, using offset-tracking with repeat Sentinel-1 images. This paper continues with a brief presentation of the study region. The method and data section illustrates the techniques used in the study and provides details about remote-sensing data and ground-based measurement data. The results and uncertainty analysis are presented in the subsequent section. Finally, after presenting and discussing the data from the automatic meteorological station and Tashkurgan meteorological station; this paper finishes with some conclusions about the surge mechanism.

Study Region
The Pamir Plateau (38 • N-41 • N; 73 • E-76 • E) is located in Central Asia, bordering the Taklamakan Desert to the east and the Karakoram region to the south, and lies partly in Kyrgyzstan, Tajikistan, Afghanistan, and China. About 10,234 glaciers are distributed across the Pamir Plateau, with a total area of about 10,500 km 2 [3]. Glacier resources in the Pamir Mountains, like in the other main mountain ranges in Central Asia, play an important role in modifying river run-off characteristics, particularly in sustaining summer river flows [71,72]. The part of the Pamir Plateau in Chinese territory is known as the Eastern Pamir Plateau (Figure 1a). The average elevation of the Eastern Pamir Plateau is more than 3000 m a.s.l. Muztagh Ata (7509 m a.s.l.) and Kongur Tagh (7649 m a.s.l.) are the two main summits in the region.  (Figure 1b), covers 115.2 km 2 (of which, 25.6 km 2 is debris-covered), and has a full length of 20 km [73]. The glacier consists of east and west branches, and lies at an altitude of 2800-7649 m a.s.l. Observations of Karayaylak Glacier's mass balance based on a digital elevation model (DEM) showed a loss of 0.32 ± 0.19 m water equivalent (w.e.)a −1 (a means one year) between 1971-1976, but a gain of 0.37 ± 0.25 m w.e. a −1 between 1999-2013/14 [74]. The west branch of this glacier surged in May 2015. This research focuses on the velocity of the west branch and the terminus (Figure 1c). The centre-line of the west branch and terminus are shown in Figure 1c. This branch has a total length of 10.65 km and is divided into 1 km subregions.

Data
We employed multisource satellite data as well as ground-based measurements. Images acquired by Sentinel-1 were used for the surface velocity retrieval. Landsat-8 and Sentinel-2 images were employed to analyze spatial changes of surge velocity and the glacier extent. Additionally, the GoLIVE dataset and the second Chinese Glacier Inventory provided indispensable references. We conducted three periods of fieldwork at Karayaylak Glacier after the surge in May 2015, during which times the surface velocity was obtained by 3D laser scanner and meteorological data were collected from an automatic meteorological station, for accuracy assessment and surge mechanism analysis.

Remote-Sensing Data
We selected 292 Interferometric Wide swath (IW), VV polarization Sentinel-1 Level-1 ground range-detected (GRD) images [61]. Sentinel-1 GRD images covering Karayaylak Glacier are available from 8 October 2014. As shown in Table 1, we filtered images according to "relative orbit number" and "slice number" under the premise of restricting the study region ( Figure 1). Further research on velocity in the surge phase required a single image covering the study region and a combination of two scenes to shorten the interval to 12 days from 30 March to 6 October 2015. This paper utilizes cloudless Landsat 8 images covering the study region to analyze glacier morphological variations during the surge phase. Due to the limitation of image quality, the Sentinel-2 image on 18 September 2016, was selected to extract the margin of the west branch and terminus of Karayaylak glacier. The Global Land Ice Velocity Extraction from Landsat 8 (GoLIVE), distributed via the National Snow and Ice Data Center, was used as a reference for glacier velocity and for choosing parameters in the velocity extraction. This dataset generates a normalized cross-correlation surface to determine the best match using Landsat 8 panchromatic images, thus yielding global glacier velocities, with accuracy range of~1 m/day to~0.02 m/day ± 5 m and a spatial resolution of 300 m [51].

Ground-Based Measurement Data
The dynamic behavior of the surge in 2015 was further observed during three field visits to the Karayaylak Glacier from 16 to 20 May, 14 to 20 August, and 4 to 7 October, 2015. The purpose of the first visit was to obtain first-hand field material after the surge; the second and third visits were used to support the research program. The specific observation aims were to (1) locate the terminus and margin of the surge region of Karayaylak Glacier, using the real-time kinematic (RTK) method based on the Global Navigation Satellite System (GNSS); (2) obtain the glacier flow speed by scanning the surge region and stagnant ice distributed in lateral moraines using a 3D Laser scanner; (3) install ablation stakes on the glacier surface to monitor glacier ablation; (4) collect continuous observations of air temperature and precipitation at Karayaylak Glacier, by setting up an automatic meteorological station (Figure 1c) at the central glacier (altitude 3398 m a.s.l.). The present study uses the field-measured ice flow speed and the air temperatures from the automatic meteorological station.
The 3D laser radar technique quantitatively measures glacier velocity and the characteristics of ice surface landforms such as drumlins, moraines and debris by establishing a high-resolution DEM and quantifying its changes. The specific method is as follows.
(1) Perform coordinate system correction, data merging, and point cloud data filtering to obtain a 3D point cloud data from two field scans; (2) obtain a high-resolution DEM of processed scanning data through the RiSCAN PRO 1.81 software. We can obtain the surface morphological characteristics by analyzing the same region in two DEMs, usually with a high spatial resolution (0.1 m). Here we used a Riegl VZ-6000, which is the latest V-series 3D laser scanner launched by RIEGL Laser Measurement Systems in Austria. The accuracy of this type of survey has been reported to be within 15 mm and 10 mm for repeated surveys provide that the degree of conformity of a measured quantity to its actual (true) value [75].
After detailed exploration of Karayaylak Glacier and its surrounding terrain, we set up a fixed GPS control point (38 • 43 9.077 N, 75 • 15 38.065 E, 3235m a.s.l.) on a lateral moraine at the intersection of the west and east branches. This was used as a datum point for the 3D laser scanner (Figure 1c). We also installed three ablation stakes (Figure 1c) at the edge of glacier terminus to monitor glacier velocity. We conducted the scan using the 3D laser scanner at the datum point on 8 August 2015, and carried out a repeat scan on 6 October. We selected 56 feature points from the 3D point cloud to calculate the velocity, thus obtaining the daily average velocity of 3 stakes and 56 feature points in total.
The high-performance portable automatic meteorological station installed in the field was a Davis Vantage Pro (Davis Instruments, USA). The instrument has an air temperature range of −40 • C to 65 • C with a resolution of 0.1 • C and accuracy of 0.5 • C [76].

Surface Velocities Retrieval
Glacier velocity was calculated using offset-tracking based on SAR intensity images. The offset-tracking method obtains the displacement between two SAR images, according to the offset of feature points in the search areas of a reference image and search image in a specific pattern-matching algorithm. The NCC [77] algorithm is widely used due to its simplicity and reliability. First, NCC calculates the correlation coefficient between a template in the first image (reference image) and subset (template) of the same size in the second image (search image). It then considers the pixel which produces the maximum correlation coefficient as the best matching in this template [78]. The correlation coefficient CC is calculated by formula (1): where (i, j) indicate the position in the search area, (k, l) the position in the reference area, r the pixel value of the reference template, s the pixel value of the search template, µ r the average pixel value of the reference template and µ s the average pixel value of the search template. It can be seen that the size of the search window is a crucial factor in NCC. The search window can contain one feature point and its neighboring pixels. This entity should become unique; thus, the ideal situation would be to find a match for every pixel independently, which is not practically feasible [78]. The search windows should be big enough to ensure that texture (and not noise) is matched, but small enough to limit displacement gradients within the window. A small window size is useful in some cases, for example shear zones or where glaciers flow over obstacles, and for small glaciers [50]. Therefore, this section tests different window sizes to identify the optimal matching. The matching result is considered to be optimized when assumed correct matches are obtained over most glacierized areas, without increasing the window size. Meanwhile, the greatest distance a pixel may shift is another parameter needed for the displacement estimation of GEE. Thus, the GoLIVE data from 2014-2020 were downloaded and converted to Tagged Image File Format (TIFF) in batches (using the Python programming language), to extract the maximum displacement of the glacier in different periods. The maximum shift and window size are considered as two input parameters to calculate the displacement. We used the Google Earth Engine to (1) transfer and mosaic images, (2) call the "Crosscorrelation" function for NCC performance, and (3) separate glaciers and export the tables and images along with the displacements. The above process was integrated into one GEE file, which automatically calculates displacement estimates for 140 image pairs in total. The processing methodology of GEE is presented in Figure 2. First, we filtered Sentinel-1 GRD images by region and period, then mosaicked images with "relative orbit number" and "slice number" of "107, 1" and "107, 2", respectively, from 30 March to 10 June 2015. Second, we calculated the displacement of the Sentinel-1 GRD images covering the study region by the "Cross-correlation" algorithm, which is provided by GEE and requires the input parameters of window size and the most significant distance a pixel may shift. Due to the considerable differences in displacements during the surge, quiescence and build-up phases, we ascertained the above parameters separately for various glacier stages (Table 2), with GoLIVE data as a reference. The "Cross-correlation" function outputs an image composed of four bands of information. The first three bands are distances: the horizontal and vertical displacements, and the Euclidean distance for each pixel in the first image to the pixel which has the highest corresponding correlation coefficient in the second image. The fourth band is the value of the correlation coefficient for that pixel. In the third step, we separated the study area into 10 bins at 1 km resolution along the centre-line of the west branch and terminus of Karayaylak Glacier. We then exported the mean displacement within each bin and the "Cross-correlation" function's output image. The margin of Karayaylak glacier is based on the second Chinese Glacier Inventory [73] and manual correction using Sentinel-2 images.

Accuracy Assessment of Surface Velocities
The accuracy of the velocity retrieval depends on the Google Earth Engine data processing and cross-correlation quality. We assess the correlation quality in this section. The accuracy of the cross correlation depends on the image quality, glacier surface characteristics and its spatial representation, co-registration error, and geocoding [58]. These influences are difficult to quantify or measure. Similar to previous studies [25,59,60], we calculated the displacement of a stable area adjacent to the glacier to estimate uncertainties in the velocity derived from different image orbits with varying time intervals. After distinguishing the non-glacier area using the Randolph Glacier Inventory 6.0 (RGI 6.0) [79], we produced 1000 randomly distributed samples over stable ground (Figure 3), then calculated the mean velocity and standard deviation (1σ) of these samples (Table 3).     [19] suggested that the complete surge phase of Karayaylak was from~13 April 2015 to~16 June 2015, while we considered the period with most intense movement as a more specific surge phase. However, this definition does not affect the consistency between the two studies. By cross-correlation calculations of Landsat-8 images from 3 October 2014 and 13 April 2015, Shangguan and others assessed the preceding active state of glacier. The consecutive velocities from 13 October 2014 to 13 April 2015, as calculated in this study, showed the steady acceleration before the surge. This is consistent with Shangguan et al., and here we extend the study period to 2020 to reveal the relatively active state of the glacier.  Figure 5b, the velocities accelerated continuously and showed a total mean increase of 874.03 m a −1 . However, the velocities of the glacier terminus did not increase. In contrast, the lower section and terminus of the west branch, 4-6 km up-glacier from the main glacier terminus, exhibited significant acceleration. Simultaneously, as shown in Figures 5b and 6a-6d, the west branch glacier terminus speeds were much higher than those of the main glacier terminus in the surge build-up phase, which indicates that this surge event originated in the west branch. (2) A continuous time series of peak velocities is displayed in Figure 5a. The maximum velocities in the study period occurred from 11 April to 10 June 2015, representing an average increase in speed of 2889.66 m a −1 . As shown in Figures 5a and 6f-6i, we found the greatest movement at the terminus of the glacier's west branch. The velocities at 3-4 km up-glacier from the terminus increased further, owing to ice advance from the west branch during the surge. It is noteworthy that the west branch terminus (bin 6) reached the highest velocity during the entire study period, of 4335.40 m a −1 , while the middle and upper parts of the west branch (bins 6-10) were still accelerating. This is in contrast to the other four periods of deceleration at bins 6-10. (3) The velocities after the surge are shown by the dotted line in Figure 5b. Velocity of the west branch showed a considerable decline, coincident with the decline in main glacier terminus velocity. The entire west branch of the Karayaylak glacier gradually stabilized, and its average velocities ranged from 216.13 to 468.46 m a −1 .

Glacier Surface Changes in Surge Phase
We selected four high-quality Landsat 8 images to investigate the morphological evolution of the glacier surface and glacier extent in the surge phase ( Figure 8). We did not detect the glacier tongue's advance. On the contrary, the west tributary terminus presented a continuous advance from 13 April to 11 July 2015, as indicated by the curve representing the terminus position in Figure 8. However, this position does not reflect the period of highest velocity, when the ice flow direction shifted towards that of the main glacier as the west branch became obstructed by the east branch. The west branch also turns northeast at 6 km from the terminus (altitude around 3450 m a.s.l.), two arrows in Figure 8 indicate the two main ice streams at this turn during the advance phase of the surge. Meanwhile, as shown in Figure 6e-j, this position also presented the highest velocity during the surge time.

The Uncertainty Analysis Based on the Field Data
We selected two Sentinel-1 GRD images of the study region, with dates of 21 August 2015 and 15 October 2015, corresponding to the "relativeOrbitNumber_start "and "sli-ceNumber "of "107 and 6 "and "34 and 7 ", respectively, to analyze the correlation between the ground-based measurements and the remote-sensing results in the same period and position. Velocities in this image pair were computed with a window size and maximum displacement of 7 and 35 pixels, thus providing the average daily velocities at 3 ablation stakes and 56 feature points (Figure 9a). As shown in Figure 9b, the NCC results represent the pixel offset distances in the horizontal and vertical directions, but the field feature points were not located in the fast-moving part, which caused low calculated velocities in the field investigation period. Thus, although many similar velocities were extracted by remote sensing, these were not matched by field observations. Nevertheless, this does not affect the matching between remote-sensing retrievals and ground-based measurements. The average absolute velocity difference between remote-sensing and field data was 0.10 m/d. The correlation is indicated in Figure 9a, demonstrating that the method used here can be used for glacier velocity monitoring with high accuracy.

Accuracy of Glacier Surface Velocity
The errors in this study include the Google Earth Engine processing error and crosscorrelation error. The correlation error was assessed in Section 4.2. However, GEE processing errors are poorly defined. To maximize confidence in our results, the ground-based measurements were compared with the remote sensing retrievals. The 3D terrestrial laser scanning (TLS) uncertainties related to the instability of the TLS [80] and registration error of scan positions [81]. The instability of laser scanning is derived from small displacement of scan positions and the vibration of 3D laser scanner. Thus, we established each scan position on stable rock surfaces using reinforced concrete. Our experience was conducted with the windless weather conditions to avoid the 3D laser scanner's vibration. The standard deviation of errors (∂ MSA ) from the set of residuals obtained from registering the point cloud can be considered as an indication of registration quality [82]. The value of ∂ MSA in this study is 0.07 m. The good correlation between these two methods certified the high accuracy of velocity extraction based on the GEE platform, which provides a powerful computing tool in glaciology.
Sund et al. [83] reported that the surge occurred in three stages. Stage 1 was characterized by an initial surface lowering and acceleration in a section of the accumulation area. Stage 2 is defined as a larger mass displacement, with continuous surface lowering in the accumulation basin and down-glacier mass transport towards the ablation area. In stage 3, the pronounced acceleration extended throughout the glacier, and the terminus advanced. Not all glacier surges show these 3 typical stages. In an internal surge, the glacier surge terminates within another tributary of the glacier or its main trunk. This type of glacier surge is mainly found in compound valley glaciers [19]. In this study, rapid advance of the glacier's main terminus was not detected during the study period. However, the west branch of Karayaylak glacier accelerated strongly and its ice tongue advanced. Therefore, the surge was restricted within the west branch of the glacier, in contrast to typical stage 3 characteristics. Stages 1 and 2 were identified by glacier thickness and velocity changes, when the glacier was in in the active phase. This paper presents velocities during Stages 1 and 2 (October 2014 to July 2015). From July 2015 to the end of 2016, velocities decreased to their lowest values of the study period. After 2017, slight acceleration and an annual peak during April to September each year indicated that Karayaylak Glacier had returned to its active state again after the short quiescence. This activity should be closely monitored.

Surge Mechanism Analysis
To analyze the surge mechanism of Karayaylak Glacier, we first correlated air temperature at the automatic meteorological station with that at surrounding meteorological stations from 15 August to 6 October 2015. This indicated that the Tashkurgan station had the best match with the automatic station (R = 0.8444). Therefore, we obtained the 1957 to 2015 air temperature at Tashkurgan station ( Figure 10) from the China Meteorological Data Network [84] and used this to analyze the glacier surge mechanism. The glacier surge mechanism includes thermal and hydrological controls. In the thermal control mechanism, rising air temperatures lead to increasing ice temperature at the glacier bed. The resulting increase in subglacial meltwater triggers the surge [14,85]. This thermal trigger only operates at the interface between cold and warm subglacial conditions, normally in the upper regions of tributaries in dendritic glaciers. Therefore, this process mostly leads to internal surges and is directly related to climate warming [85]. In the hydrological control mechanism, excess meltwater formed by warming temperatures or other reasons cannot drain out quickly due to changes or obstructions in the subglacial drainage system [86,87].
The above thermal control is consistent with the Karayaylak Glacier surge and many glaciers in the Karakoram [25]. Thus, we considered this as the primary mechanism leading to the surge of Karayaylak Glacier. However, the thermal control cannot fully explain the surge initiation. As shown in Figure 10, the rapid air temperature increase at Tashkurgan after 1998 caused basal ice temperature to reach the melting point, suggesting that the surge was enabled by conditions at that time. However, the thermal control mechanism theory does not consider what factors trigger the surge. Besides, surges based on the thermal control mechanism may occur at any time of the year, while surges based on the hydrological control mechanism generally occur in winter [86]. The timing suggests this surge was more likely to have been triggered by the hydrological mechanism. Here, our detailed analysis suggest that the thermal control mechanism established optimal surge conditions, but the hydrological control mechanism triggered the surge; therefore, both processes were involved, as described in the following section.
The "upper cold, bottom warm" polythermal glaciers [87] distributed in the Karakoram, Pamir Plateau and Tianshan include the intensively-studied Batura glacier [88], and Qingbingtan No. 72 glacier, Tuomuer Peak [89] which is less than 600 km straight-line distance from Karayaylak Glacier. Although this study lacks englacial ice temperature data, on the basis of the very wide elevation difference (7649 to 2780 m a.s.l.) we speculated that Karayaylak Glacier belongs to the "upper cold bottom warm" type. Karayaylak Glacier exhibits distinct cold and warm ice areas, as indicated by the respective absence or presence of well-developed drainage tunnels. In the context of climate warming and rising temperatures, we infer that the cold-warm boundary has extended rapidly up-glacier, leading to increasing meltwater production. This excess meltwater will be evacuated quickly by the existing subglacial drainage system in the warm zone. However, as the warm zone expands to areas without such tunnels, the hydraulic pressure will rise sharply as meltwater accumulates. This dynamic state, coupled with meltwater lubrication, is most likely to result in the initial glacier sliding, which triggers a glacier surge through positive feedback [87].
Analysis of Tashkurgan air temperature from October 2014 to July 2015 indicated the anomalous warming from February to April in 2015 (Figure 11), when maximum daily air temperature in February 2015 reached 4.3 • C. Here we focus on 12 consecutive days of high-temperatures from March 21 to April 1 (Figure 11), when the average temperature (7.1 • C) was 7.2 • C higher than the mean for this period. Using an air temperature lapse rate of 0.7 • C (100 m) −1 [90], this daily mean temperature corresponds to −0.8 • C at 4300 m altitude, with a daily mean maximum temperature of 5.3 • C. Twelve days at these temperatures will have caused intensive ice surface ablation and melt water production; due to poorly-developed supraglacial and englacial drainage, this meltwater will have entered the ice through cracks and the glacier margin. This was most likely the main trigger of the surge. The possibility that the surge was triggered below the boundary between cold and warm ice is discussed further as follows. Before the surge, meltwater from the warm ice area below the cold-warm transition can drain easily though the existing drainage system, and will not trigger a surge [91]. However, if drainage tunnels are blocked for some reason, hydraulic pressure will rapidly increase, and eventually trigger a surge [85]. In theory, four reasons can cause changes to drainage tunnels: (1) increased pressure in the accumulation zone increases shear stress in basal ice; (2) subglacial temperature changes modify the temperature distribution in the ice and bedrock; (3) increased meltwater erodes drainage tunnels and causes them to collapse; (4) the ice bed comprises soft sediments which support only a weak drainage system that is easily altered.
At first, we thought that the first factor (increased shear stress in basal ice) was the main trigger, because Karayaylak Glacier's mass balance increased by 0.37 ± 0.25 m w.e. a −1 from 1999 to 2013/2014 [74]. However, according to the meteorological data for Tashkurgan, we found the precipitation from January to April of recent years was lower than the long-term average, especially during 2014 and 2015. Our inference would then depend on inaccurate mass balance data due to a lack of snow depth datasets. Thus, the shear stress trigger is lacking in supporting data. While we note that the accumulation rate on the upper glacier remains uncertain due to a lack of precipitation observations from the upper glacier, the first factor can be eliminated, and the other three factors causing changing subglacial drainage remain plausible. Therefore, the possibility that the surge was triggered by the hydrological control mechanism below the cold-warm transition is very high. Besides the "hydraulic effect" of meltwater, the "lubrication effect" cannot be ignored in these processes [91].
Our comprehensive analysis has indicated that the Karayaylak Glacier surge was likely triggered in the upper part of the glacier, but quickly propagated to cause a surge in the entire glacier branch. Continuous climate warming very likely caused the increasing meltwater and created favorable conditions for the surge. After October 2014, the glacier terminus and west branch velocity nearly doubled, indicating the onset of basal sliding or subglacial sediment deformation, prior to the main trigger. Warm air temperatures from 21 March to 1 April provided this trigger. This clearly shows that the surge mechanism is not a single process: instead, it comprises coupled thermal and hydrological controls, which play different roles at different times and places. The two mechanisms interact to produce superimposed effects.

Conclusions
This study retrieved the surface velocities of the west branch and terminus of Karayaylak Glacier from 13 October 2014 to 17 October 2020, using Sentinel-1 GRD imagery and the GEE platform. The glacier dynamic and remote-sensing retrieval methods are focused on in this paper.
The west branch of Karayaylak Glacier had presented a substantial acceleration in October 2014. Velocities in May 2015 reached a maximum of 4335.40 m a −1 at the west branch terminus which showed strenuous movement with its advance of about 3 km. Thus, we concluded that this glacier tributary surge happened in May 2015 accordingly. Analysis of data at the automatic meteorological station and Tashkurgan meteorological station suggested that this surge resulted from coupled thermal and hydrological controls, with both playing different roles at a different times and altitudes. The velocities decreased sharply within two months after the surge and then remained slow for nearly a year. Since 2017, velocities have increased slightly and have shown the annual periodic variability, which revealed the current relative dynamic state of Karayaylak Glacier and should be closely monitored.
The accuracy assessment of field measurement and velocity retrieval was presented in this paper. Moreover, the high correlation between the field data and the remote sensingderived velocities supported this work's reliability. This demonstrated the high-precision monitoring method of the glacier surge based on the Sentinel-1 imagery and NCC, indicated in particular the great potential of the GEE platform to perform cloud-based computing in glaciology. Data Availability Statement: Surface velocities retrieval results are available upon request by email to the author (2019212427@nwnu.edu.cn).