Susceptibility Analysis of Glacier Debris Flow by Investigating the Changes in Glaciers Based on Remote Sensing: A Case Study

: Glacier debris flow is one of the most critical categories of geological hazards in high-mountain regions. To reduce its potential negative effects, it needs to investigate the susceptibility of glacier debris flow. However, when evaluating the susceptibility of glacier debris flow, most research work considered the impact of existing glacier area, while ignoring the impact of changes in glacier ablation volume. In this paper, we considered the impact of the changes in the glacier ablation volume to investigate the susceptibility of glacier debris flow. We proposed to evaluate the susceptibility analysis in G217 gullies with frequent glacial debris flow on the Duku highway, Xinjiang Province. Specifically, by using the simple band ratio method with the manual correction to identify glacier outlines, we identified the ablation zone by comparing the glacier boundary in 2000 with that in 2015. We then calculated ablation volume by changes in glacier elevation and ablation area from 2000 to 2015. Finally, we used the volume of glacier melting in different watersheds as the main factor to evaluate the susceptibility based on the improved geomorphic information entropy (GIE) method. We found that, overall, the improved GIE method with a correction coefficient based on the glacier ablation volume is better than the previous method. Deglaciation can be adapted to analyze glacier debris flow susceptibility based on glaciology and geomorphology. Our presented work can be applied to other similar glacial debris flow events in high-mountain regions.


Introduction
Debris flow is one of the extremely dangerous and threatening geo-hazards in high-mountain regions, which induce severely destructive outcomes to human life [1][2][3]. Glacier debris flow is a type of debris flow characterized by suddenness and high mobility in mountainous areas due to the loss of glaciers caused by climate change [4,5].
To minimize the loss induced by glacier debris flow, much research work has been conducted to investigate the formation mechanism of glacial debris flow in recent years. This work has concluded that glacial melting is one of the main reasons for the eruption of glacial debris flow [4]. It is necessary to evaluate the dangerousness and stability of glacial debris flow [6].
In order to investigate the typical characteristics and significant features of glaciers, much research work based on remote sensing has been conducted. Compared with conventional field-based survey, remote sensing-based methods have the advantage of large-scale coverage and convenient and real-time acquisition to monitor dynamic changes of glaciers [7]. For example, Bayr, K.J. and D.K. Hall [8] observed the glaciers recession and the delineation of the ice in the eastern Austrian Alps by using remote sensing data. Kulkarni, A.V. [9] used Indian remote sensing satellite data to map the glacial area, and obtained deglaciation of an entire glacier. In addition, remote sensing images can be used to provide observations of debris flow gullies and watershed morphology, which are necessary to assess the susceptibility and prediction of glacier debris flow. For example, A. Kääb1. [10] obtained geo-hazards susceptibility analysis in high-mountain areas through object classification and monitoring techniques based on remote sensing data. Elkadiri R. [11] developed a geographical information system (GIS) and evaluated the susceptibility of debris flows based on remote sensing data in the Jazan province of Saudi Arabia.
Much research work focused on the formation conditions of the glacial debris flow has been conducted. For example, Cheng, Z.L. [12] discovered that the loss of glacier area due to climate change can trigger debris flow in southeastern Tibet. Deng, M.F. [13] found that compared with rainfall-type debris flow, the occurrence of glacial debris flow is caused by the retreat of glaciers and the melting of internal ice particles and the changes in glacier can provide large volumes of active sediment.
Other research work focused on the susceptibility and forecasting of glacial debris flow has also been conducted. Risk assessment of debris flows is increasingly being considered an important method for dealing with disaster mitigation. For example, Golovko, D. [14] applied the susceptibility and forecasting mapping of debris flows to prevent regional debris flow. However, much research work has only considered rainfall conditions. For example, Cannon, S.H. [15] used the duration of rainfall to predict debris flow. Chen, N.S. [16] investigated the essential rainfall features for debris flows in the Wenchuan earthquake area. Rickenmann, D. [17] selected starting zone, critical rain conditions triggers, sediments source, potential event scale, and event frequency as the main factors to evaluate the occurrence of debris flow. Allen, S.K. [18] proposed a systematic analysis method by combining hydrological factors with topographical factors to forecast Kedarnath disaster.
Compared with rainstorm debris flow [19,20], little research work has been focused on glacial melting types of debris flow. Furthermore, when studying the factors that trigger the type of glacier melting debris flow, the triggering conditions are mostly considered as the two factors of temperature and rainfall without the factor of glaciers [21]. Little research work has been used the existing glacier area as a factor to evaluate while not considered the impact of the glacier changes [22]. Moraine deposits caused by glacier melting provide abundant provenance for debris flow motion, and glacier melting water and rainwater become the source of water for debris flow simultaneously [23].
To address the above issues, in this paper, deglaciations and icefalls were considered to be one of the main trigger factors. We present a case study analyzing the susceptibility by investigating the changes in glacier. Specifically, we first identified the ablation zone by comparing glacier boundaries. Four methods were used for classification, which include the simple band ratio method, the normalized difference snow index method (NDSI), the supervised classification of maximum likelihood method, and the supervised classification of neural networks. We chose the first method with the manual correction to identify glacier outlines by comparing the four methods. We then calculated the glacier area where it melting in the study area and computed glacier elevation changes by digital elevation model (DEM) differencing from 2000 to 2015 to calculate glacier changes of volume. Finally, we selected a target area for susceptibility analysis by the modified geomorphic information entropy (GIE) method.
Our contributions in this paper can be summarized as follows: (1) We identified the ablation zone by comparing glacier boundaries in 2000 and glacier boundaries in 2015 in the study area near Duku highway; (2) We calculated ablation volume by changes in glacier elevation and ablation area from 2000 to 2015; (3) We used the volume of glacier melting in different basins for glacier debris flow susceptibility based on the improved GIE method.
The rest of this paper is organized as follows. Section 2 will introduce the materials and methods in detail. Section 3 will describe the results by investigating the changes in ablation zone based on remote-sensing images and evaluate the susceptibility analysis for the glacier debris flow based on the modified GIE method with correction factor by glacier melting volume. Section 4 will discuss the advantages and shortcomings of our work, and also point out our future work. Finally, Section 5 concludes our work. High-mountain regions are vulnerable to permafrost melting and glacier recession caused by global warming [5,24]. The Tien Shan Mountains area is the main high-mountain distribution area of glaciers in Xinjiang Province, China (Figure 1a). It is necessary to study the susceptibility of glacial debris flow in typical areas of Tien Shan Mountains.

Materials and
Duku highway is located in the middle part of the Tien Shan Mountains, which has a north-south trend leading-from Dushanzi to Qiaoerma (Figure 1b). The Duku highway has a special terrain with steep slopes and more than 280 km of road with altitudes above 2000 m (Figure 1c). This highway is vulnerable to natural disasters due to natural factors such as snowfall and glaciers in Tien Shan Mountains.

Details of the Target Area
The glacial debris flow in the Tien Shan Mountains is controlled by temperature changes in the high-mountain area. With the condition of continuous high temperature, the infiltration of glacial meltwater with the subsequent high-intensity storm easily causes glaciers or glacier-rainstorm mixed mudslides in the watersheds.
We calculated glacier changes in the study area, and then selected a section with frequent debris flow from K580+000 to K680+000 in the Duku highway as a target area. The target area is a concentrated distribution area of modern glaciers with a watershed area of approximately 79.724 km 2 ( Figure 2). This area has the characteristics of the long gullies, steep slopes, and bended routes due to melting glaciers and most of the debris flow gullies in this area. The glacier tongue is dominated by unsteady hanging glaciers, which readily cause ice collapse in the case of continuous high temperatures. Collapsing ice objects are also easy to wash quickly into the gully. As a result, melting glacial water cuts out mud-rock flow gully bed from narrow to deep, which provides good gully bed conditions [25]. Melting glaciers and ice collapse are some of the major factors that induce debris flow. For example, the K636 mud-rock flow watershed has blocked the Lapat River several times and formed a mud-rock flow barrier lake, which caused the valley where the accumulation fan is on the upper reaches of the main river has changed from the original canyon to the wide valley ( Figure 2a). With the loose residual materials in the gullies caused by glacier ablation, it was essential to obtain the debris flow risk assessment of this section.

Remote Sensing Image Data
In our study, we used Landsat series data with a spatial resolution of 30 m for the period from 2000 to 2020, which included thematic mapper (TM) and operational line imager (OLI). We also used DEMs, which included SRTM DEM and TanDEM-X to observe glacier changes.
Landsat data have been used to map glacier outlines for many years [26]. For example, Andreassen, L.M. [27] introduced the method of extracting glacier extents through the Landsat thematic mapper (TM) with aerial images and digital maps. Dozier, J. [28] distinguished snow from other surface cover and cloud layer by the Landsat TM. Frey, H. [29] used data involving the enhanced thematic mapper plus (ETM+) images from 2000 to 2002, coherence images from ALOS PALSAR image pairs, the SRTM DEM and the ASTER global DEM (GDEM) for glacier mapping.
Cloud and snow coverage are the main reasons for the suitability of remote sensing data for glacier observation and change monitoring analysis. We selected data with minimal shadow-cover, low-coverage of snow and clouds in this research. The details of the remote sensing data are listed in Table 1. Table 1.
Details of the remote sensing data used in the study (Landsat data source: http://earthexplorer.usgs.gov/ (accessed on 23 June 2020)).

Images Resolution (m) Snow-Covered Applicability
Landsat TM Much research work has been conducted to investigate the change in glacier elevation and the loss of glaciers based on the DEM data. For example, Gardelle, J. [30] measured glacier elevation and calculated the changes in volume from DEM differencing. In this study, DEM was used for dividing the sub-basin by extracting the topographical parameters and computing the change in glacier elevation. The details of the DEM data are listed in Table 2. In this paper, the processing procedure for susceptibility assessment of glacier debris flow consisted of three main steps: (1) The first step is to identify the glacier ablation zone, where glaciers were present in 2000, but disappeared in 2015. We extract the ablation zone by comparing the glacier boundaries between 2000 and 2015; (2) The second step is to calculate the changes in glacier charateristics, which included their areas and elevations. We then calculated the volume of the abaltion zone by multiplying the area and elevation changes from DEM differencing using an integral method; (3) The third step is to convert the volume of glacier melting into a correction factor to be applied to the debris flow susceptibility analysis. These above steps are the basis of the susceptibility analysis of debris flows. The entire workflow of the susceptibility analysis is illustrated in Figure 3.

Step 1: Identification of Ablation Zone
Due to the gradual warming of the global climate, the area and volume of glaciers have decreased considerably, and even disappeared. Glacial melting is the main method of glacial material depletion. The methods of glacier melting include ice surface melting, ice melting, and sub-ice melting, and ice surface melting is the main method. In this paper, we focused on the impact of glacier changes in the ablation zone on the occurrence of glacial debris flow.
(1) Pre-processing of the Remote Sensing Images Due to the interference of solar radiation, clouds, topography, etc., it is necessary to preprocess remote sensing images to avoid the loss in the process. These steps of remote sensing image preprocessing included radiometric calibration, atmospheric correction, geometric correction, and image sharpening of Gram-Schmidt (GS) pan sharpening. These steps can reduce the impact, like noise and improve the quality of images before the classification operation.
(2) Main Processing: Identification of Glaciers Boundaries Satellite data have been commonly used for glacier outline identification. For example, by stereoscopic interpretation with aerial photo and Landsat TM data, Aniya, M. [31] identified glacier boundaries and glacier watersheds. Gjermundsen, E.F. [32] calculated the changes of glacier area and extracted the 2002 glacier boundaries in the central Southern Alps, New Zealand. Accurate glacier boundaries are quite important to the study of glacier changes. When glacier extents are regarded as a primary factor for change analysis, it is essential to know how accurate they are. Therefore, it is necessary to obtain with the manual correction of the original outlines (e.g., for debris cover, shadow) [32]. Creating an automated glacier-mapping method that is compared with digital boundaries can be calculated for selected glaciers in the study area. The identification of glacier boundaries based on satellite data was inspired by related work. Among the methods for glacier classification, the common method is the simple band ratio method, which has a fast speed, high precision, and great robustness. A number of studies on glaciers are based on this method. For example, Andreassen, L.M. [28] obtained glacier mapping by threshold ratio images (TM3/TM5). Le Bris, R. [33] used this method based on an application of a threshold (TM3/TM5) to automatically map accurate glacier outlines through the manual correction of misclassified lakes in the study area. This method calculates the reflectance ratio between two bands. The basic principle is to use different reflectivities for different wavebands to distinguish glaciers from other ground objects. The accuracy of this method depends on the selection of thresholds. There are other commonly-used methods, such as NDSI, the maximum likelihood method, and the neural network method. NDSI is computed as (Green-SWIR)/(Green+SWIR) for extracting glacier extents. The maximum likelihood method and the neural network method generally belong to the algorithms of supervised classification, which classify new sample pixels by establishing a typical training sample [34]. As illustrated in Figure 4, it shows different glacier boundaries based on the four methods. The classification results of these methods are as follows. Due to the NDSI method has low recognition of moraine-covered glaciers [35,36], it is easy to identify bare land as glaciers, and the results of glacier classification are fragmented. The maximum likelihood method misclassifies rivers as glaciers [37]. In addition, it is easy to misclassify snow and bare land as glaciers, causing the salt-and-pepper phenomenon. The maximum likelihood method is easy to misclassify snow as glaciers, and the classification results are also fragmented. The neural network method is more accurate than the maximum likelihood, while the training time of this method for samples is relatively long, and in this case, the snow, bare land, and glaciers cannot be well-identified. Although the band ratio method creates some cavities inside the glacier with a too small threshold, it can be solved by adjusting multiple times on the selection of thresholds according to the actual situation. Among the four methods, the simple band ratio method is the most suitable in this study. The results based on Red/SWIR for classification are better than those based on NIR/SWIR after testing. (3) Post-processing of the Glaciers Boundaries In the results of automatic classification of glaciers based on remote sensing images, there are usually many small patches or holes. Thus, we need to manually correct interpretation of glacier outlines after the classification operation is implemented. These patches can easily be selected and deleted in the layer by manual editing to improve the accuracy of the glacier outlines. Furthermore, the automatic classification method of glaciers can lead to the misclassification of glaciers (e.g., snow cover, rivers, and clouds). The proper correction of glacier boundaries can use high-resolution images, such as Google Earth images, Landsat series data by overlay of three individual bands, and even DEM data. These high-resolution data can obtain a classification result of the snow and cloud conditions, even shaded areas. After manual correction, we compared the glacier boundary in 2000 with the glacier boundary in 2015, and extracted the ablation region where the glacier existed in 2000 but disappeared completely in 2015. The ablation zone is illustrated in Figure 5.

Step 2: Calculation of the Changes in Glacier Ablation Zone
Temperature and rainfall are the two main control factors for development of glaciers, Temperature determines the ablation area, and rainfall determines the accumulation area [38]. We analyzed the average temperature and precipitation data in the Tien Shan Mountains area for the past 40 years. Although this area receives a great amount of precipitation, glaciers are more sensitive to temperature [39]. Falatkova, K. [40] predicted that the area of glaciers will decrease significantly by 2050 based on climate model scenarios. The replenishment of solid precipitation to the glacier cannot offset the material loss caused by glacier melting. The glacier recession will be accompanied by increases in glacier run-off and the formation of new lakes. As a result, most of the glaciers in the area are shrinking and melting is accelerated. In this step, we calculated the changes in glacial ablation zone through the three aspects of area, elevation, and volume.
(1) Determination of Ablation Area Glaciers are sensitive to climate change, even that of very slight temperature changes. With a background of a globally warming climate, any excess energy can melt glaciers, leading to glacier shrinkage. Glacial flow transports material from the accumulation area to the ablation area, where it melts [31]. Rainfall and glacier melt are the main reasons for increasing surface run-off in response to climate change. In this study, we focus on ablation area, where all glaciers are transformed into water to trigger glacial debris flow. Based on the identified ablation zone boundary range, we determined the area of the glacier ablation.
(2) Calculation of Elevation Changes from DEM Differencing Measured elevation changes have been used to investigate volume and mass changes in glacier for many studies. Accurate DEM differencing is key to calculate glacier elevation changes. Therefore, it is necessary to process the registration of the DEMs before calculating the glacier elevation changes. Due to the long wavelength of the C-band, there is an obvious penetration effect on the glacier, and the penetration depth is between 0∼10 m, especially in the glacier accumulation area [41]. Therefore, we preprocessed the SRTM DEM before DEM registration, as a result, the C-band SRTM DEM data corrected by the X-band SRTM DEM can represent the glacier surface height [30,42]. In this paper, the DEMs must be projected to the same map system, such as the universal transverse Mercator (UTM) system. In addition, because of the different spatial resolutions, one of the DEMs needs to be resampled to ensure the accuracy of the results. The abnormal values of surface elevation changes are altered through the Kriging interpolation algorithm. The elevation value we used now was proposed by McVicar, T.R. [43]. The method of processing the changes in elevation was taken from previous research work. The most frequently used DEM co-registration method, which is simple, analytic and robust was proposed by Nuth, C. and A. Kaab [44], and is based on the elevation derivative of slope and aspect and elevation difference residuals to correct the registration errors (Equation (1)).
where dh is the elevation difference between the different DEMs, α and Ψ are the terrain slope and aspect in the region, a is the magnitude elevation difference bias, b is the direction of the shift vector, c is the mean bias between the DEMs divided by the mean slope tangent of the terrain. Slope and aspect can be calculated by employing GIS. Previous research work has shown that the offset obtained by this method may not be the most suitable offset, and this process needs to be iterated to reach the optimal solution. In this paper, we iterate this process until the solution is solved. The process is terminated when the standard deviation is less than 2% to obtain the best registration result. The bias between DEMs is effectively eliminated through the use of fitting and solving to obtain the correction parameter pair after the DEM is registered. The result of DEM co-registration is illustrated in Figure 6. The glacier elevation changes were calculated based on the mean elevation difference between the DEM of two acquisition periods (2000-2015).

(3) Calculation of Ablation Volume
Compared with the area of glaciers, we are more concerned about the changes in glacier volume. The volume of the ablation area is calculated by multiplying the area and elevation changes using an integral method. The combination of temporally consistent ablation area and elevation measurements is the crucial advantage of our method, which improves volume estimate accuracy. We improved the two-dimensional changes in the glacier to the three-dimensional changes in our study to obtain the susceptibility analysis of glacier debris flow based on glacier ablation.  2000) and TanDEM-X (2015) before (a) and after (b) co-registration. Histograms of the elevation differences before (c) and after(d) co-registration.

Step 3: Susceptibility Analysis of Glacier Debris Flow
Glacial debris flow is a serious issue in high-mountain areas that need to be prevented and mitigated by susceptibility assessments. Glacial geomorphology can be used for completing susceptibility assessments of glacier debris flows. For example, Huggel, C. [45] assessed the susceptibility of glacial disaster, such as ice avalanches and glacial debris flows based on glaciological, geomorphological, and hydraulic principles. Zaginaev, V. [46] found that glacier lake outburst floods and glacial debris flows in northern Tien Shan area are negatively correlated with the rates of glacier area and shrinkage ratios of moraine-glacier. These studies showed that the susceptibility methods based on glaciology and geomorphology can improve the assessment accuracy of glacier disasters. Our method of evaluating susceptibility was inspired by the related work.
(1) The Divisions of Sub-watersheds in the Target Area We divided sub-watersheds in the target area considering hydrological statistics and an analysis module based on ASTER data. The main process is swale filling, traffic analysis, confluence accumulation statistics, water flow length calculation, river network extraction, etc. The result of watershed division is illustrated in Figure 7. (

2) Brief Descriptions of the Geomorphic Information Entropy Method
Geomorphic condition is one of the three conditions necessary for the formation of debris flows, it is the background condition and potential energy condition of debris flows. The geomorphic information entropy method (GIE) of erosion watershed systems was proposed by Ai and Yue (1988) [47]. The GIE method depends on the area-elevation curve (Strahler) and combined the information entropy theory (Equation (2)). This theory has been validated to the evaluate of debris flows developed from traditional geomorphology. Some research work have been evaluated the susceptibility of debris flow based on the watershed erosion degree. For example, Wang, J. [47] used the geomorphic information method to evaluate the susceptibility in the strong earthquake area, Wenchuan. Li, Y.C. [48] used this method to obtain the susceptibility of each sub-basin in debris flow gullies.
where H is the GIE value, S is the area-elevation integral value, and f (x) is the fitting function. Due to glacier-covered and rich sediments, glacial geomorphology has particularities. From the view of thermodynamics, glaciers are the process of increasing entropy under high temperature conditions. Furthermore, glaciers increase the instability of the entire basin and the probability of glacial debris flows. From the perspective of geomorphological development, when glacier coverage is below a certain value in glacier-covered areas, the risk of glacial debris flow is positively associated with glacier coverage. However, beyond this value, the large area of glaciers has a protective effect on the evolution of watershed geomorphology, but most glacier-covered areas have low glacier coverage. Based on related research work, glacier coverage is 3.4% in the Tien Shan area in total. Therefore, the risk of glacial debris flow is positively associated with glacier coverage in the study area. According to this theory, a modified method was proposed by a glacier correction factor based on the existing glacier area (Equations (3) and (4)).
where i is the ratio of glacier area to basin area, G is the non-glacial coverage rate within the watershed. On this basis, we turn the glacier area into the volume of melting glacier to improve the accuracy of debris flow susceptibility. We divided the glacier volume level into 11 levels, which correspond to the correction factor by volumetric range. The values of the correction coefficients are as follows, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95. By calculating the value of GIE according to the improved correction coefficient, we evaluated the degree of development and risk of glacial debris flow.
Previous criteria were proposed by Ai and Yue in 1988 [47], which divided the development of different basins into three stages. In this paper, we appropriately adjusted our criteria for division. We divided the basin into five different stages: H 0.11, the development of debris flow is in infancy stage with strong erosion, the slope is steep with mainly convex slopes, glacier coverage is extremely high, the susceptibility of debris flow is very high; 0.11< H 0.20, the development of debris flow is in mature but bias infancy stage, erosion begins to weaken, the slope begins to transform into concave, glacier coverage is high, the susceptibility is high; 0.20 < H 0.30, the development of debris flow is in mature stage, erosion is obviously weak, the slope shape is further transformed into concave, the glacier coverage is moderate, the susceptibility of debris flow is moderate; 0.30 < H 0.40, the development of debris flow is in mature but bias old stage, erosion weakens clearly, the landform surfaces is gently undulating, the glacier coverage is low, the susceptibility of debris flow is low; H > 0.40, the development of debris flow is very old stage, most of the terrain glacier coverage is very low, the susceptibility is very low.
The GIE value under the glacier landform is obtained by calculating the correction factor. This GIE method can be used to divide the development stage and obtain the risk assessment of glacial debris flow based on the above-mentioned criteria.

The Observation of Glacier Ablation and Recession
We observed the glacier ablation zone in the study area and calculated glacier changes in different regions during the observation period from 2000 to 2015. Due to the severely fragmented glacier distribution near the Duku highway, glacier change measurements have been aggregated over regional subdivisions to show the scale of ablation more intuitively.
From the results, mostly negative mean elevation changes (<−0.8 m) are recorded in the study area with a small part of region having positive values at lower glacier elevations (Figure 8a). In many regions, most glacier elevation changes are negative throughout all altitudes, indicating that the glacier is thinning and downwasting, as small and fragmented glaciers are completely melting. In addition, the accumulation area is much smaller than the ablation area. Only a small portion of the glacier elevation changes in the accumulation area formed by the transportation of glacial meltwater are positive (>0.4 m).
From the results, we observed the region that has the largest melting area and the melting volume in the study area. It was showed that the largest ablation area is approximately 15.6 km 2 and the largest ablation volume is approximately 21.2 km 3 (Figure 8b,c). The possible reason is that the glaciers shrank drastically during the research period and there were many small ice bucket glaciers located in valleys and hanging glaciers on steep mountain slopes in this largest ablation area, especially hanging glaciers with short ice tongues and thin thicknesses from a few meters to tens of meters, which are vulnerable to regional climate change. With increasing ice temperature, the number of ice fissures increases, glacier fragmentation increases, and the melting surface increases.
We have obtained relevant records on the glacier changes according to related books and references [49][50][51][52]. For example, Schiefer, E. [49] quantified the changes in glacier volume in British Columbia, Canada for the period 1985-1999 and investigated regional glacier thinning rates. Tak, S. [50] estimated the surface mass loss in Parvati glacier using 19 years of satellite images. Romshoo, S.A. [52] observed glacier recession in the Kashmir Himalaya in India from 1980 to 2018 and the observed glacier loss is higher (0.77 ± 0.31 km 2 a −1 ) compared with the other Himalayan regions. According to related references, the results of glacial changes, which include glacier melting area, elevation change, and volume variation are all within a realistic and reasonable range. Related research can validate the rationality and accuracy of our investigation results.
In regions with larger area and volume of glacier melting, it is necessary to obtain the risk assessment of debris flow caused by a large amount of glacier melt water will flow into the river course to cause changes in surface run-off or form dammed lakes, especially with conditions of torrential rains and steep terrain. The observation of glacier changes in the entire study area was promoted to obtain the distribution of glaciers and glacier ablation in the target basins, which conduced to analysis the type determination and mechanism analysis of glacial debris flow, and facilitates the improvement of the accuracy of debris flow susceptibility analysis.

The Influence on the Glacier Shrinkage and Downwasting
The altitude of the glacier has a substantial influence on the retreat and melting of the glacier. The high temperature at low altitude caused glacier to recede rapidly [53,54]. Through glacier classification based on the main altitude of the glaciers, we observed the altitudinal impact on glacier ablation. The result is illustrated in Figure 9. It was observed that the main area of glacial melting was found at altitudes varying between 3600∼4600 m during the study period. The highest ablation area of approximately 14.80 km 2 and greatest volume of approximately 21.80 km 3 were witnessed in the glaciers situated between 4000 and 4200 m altitude, whereas the lowest ablation area of approximately 0.03 km 2 and lowest volume of approximately 0.04 km 3 were observed in the glaciers situated above 4800 m altitude. The reason for the location of highest ablation region may be that more glaciers melt because there are many small glaciers, such as hanging glaciers, in this altitude range. When located above a certain altitude, glaciers rarely melt due to low temperature and precipitation inputs. In general, the altitudinal influences on the changes in volume and area are basically the same in high-altitude areas where glaciers melt. When the average altitude of the glacier distribution area is higher than 4800 m, the susceptibility of debris flow can be analyzed temporarily without considering the influence of glacier melting.  This aspect also has an important influence on glacier ablation. Glacial ablation is more heavily focused on southern aspects caused by the high intensity and long duration of solar radiation [55]. However, it was observed that north-oriented glaciers melted more than the south-oriented glaciers from 2000 to 2015 ( Figure 10). The result is illustrated in Table 3, and the reason for the analysis of this result may be the fragmented glaciers distributed in the north-facing direction. According to the glacier changes in various slope directions in the study area from 2000 to 2015, glacier ablation regions were mainly distributed on the north slope, west slope, northeast slope, and northwest slope. The glacier melting areas to the north and northwest melted highly, the west and northeast slopes were in the middle, and the glaciers on the south, southwest, east and southeast slopes melted less.

Susceptibility Analysis of Glacier Debris Flow
Based on the calculated changes in ablation area and elevation, it is easy to extract the ablation volume in different watersheds, as illustrated in Table 4. Through our standard of dividing the ablation volume level, it is converted into the corresponding correction coefficient.
After calculating each elevation class area of each sub-basin by reclassifying the DEM files of the sub-watershed at 100 m intervals, we needed to process the data in Excel. We then fit the data by using a cubic polynomial and obtained the Strahler area-elevation curve and fitting function. Finally, the results of the Strahler area-elevation integral value and geomorphic information entropy are obtained through this process. We directly quoted the results of related studies. The range of the GIE values is from 0.09 to 0.37. The susceptibility level of each watershed based on the glacier coverage is illustrated in Table 5. The range of the modified GIE values is from 0.0896 to 0.333. By multiplying by the correction factor based on glacier ablation volume, the susceptibility level of each watershed by our division criteria is illustrated in Table 6. The range of the modified GIE values is from 0.0855 to 0.296. Based on the GIE method with the consideration of glacier coverage, the susceptibility classification result is illustrated in Figure 11a. The results showed that in the 13 watersheds, 1 sub-watershed exhibited infancy development and very high susceptibility; 5 sub-watersheds exhibited mature development but biased infancy development and high susceptibility; 6 sub-watersheds exhibited mature development and moderate susceptibility and 1 sub-watershed exhibited low susceptibility.
Based on the GIE method with the consideration of the glacier ablation volume, the susceptibility classification result is illustrated in Figure 11b. Specifically, in the 13 watersheds, 3 sub-watershed exhibited infancy development and very high susceptibility; 8 sub-watersheds exhibited mature development but biased infancy development and high susceptibility; and 2 sub-watersheds exhibited mature development and moderate susceptibility. The watershed area of very high susceptibility covers approximately 10.09 km 2 , which accounts for 12.65% of the total basin area. The watershed area of high susceptibility covers approximately 44.97 km 2 , which accounts for 56.4% of the total basin area. The watershed area of moderate susceptibility covers approximately 24.66 km 2 , which accounts for 30.9% of the total basin area. Comparative results of the two susceptibility maps show that several sub-watersheds move from moderate susceptibility to very high susceptibility. As illustrated in Figure 12 (e.g., K636 gully), the reason for the change of the susceptibility may be that the more glacier terminus is distributed in the K636 gully. When Moraine deposits, caused by glacier melting, provide abundant provenance for debris flow motion, and glacier melting water and rainwater become the source of water for debris flow simultaneously. The former susceptibility analysis only considers geomorphic features such as transport distance and elevation while not consider the impact of the glacier changes. After considering the factor of glacier melting, the susceptibility level of its sub-watersheds increased. Based on the field survey, combined with related work and statistics from highway maintenance departments, we collected the data from a total of 15 glacial debris flow disaster events in G217 gullies in the basin from 2001 to 2019, of which 2 were large-scale debris flows, 4 were medium-scale debris flows, and 9 were small-scale debris flows. It shows that debris flows events frequently occur in the region, which could damage major projects, traffic corridors, and urban construction. The statistics on the occurrence of glacial debris flows events based on the maximum scale in the region since 2001 are listed in Table 7. The distribution of glacial debris flow events in the G217 gullies is illustrated in Figure 13. The susceptibility analysis we proposed by investigating the changes in glacier ablation based on remote sensing data can improve the accuracy that risk assessment of the K636 and K637 watersheds deemed to be at risk. According to the actual number and scale of debris flows from some historical information (Figure 13), it can be judged that the susceptibility of the improved GIE method is more accurate.

Discussion
In this paper, we used the improved GIE method to evaluate the susceptibility of glacier debris flows in high-mountain areas, which makes our susceptibility assessment more accurate and realistic. There are two specific contributions in our work. First, we used remote sensing data for glacier observation and change monitoring analysis. Remote sensing data with large scale coverage and real-time acquisition is more convenient and reliable in glacier mapping. Second, we used a correction coefficient based on glacier ablation volume to evaluate susceptibility with our modified division criteria. The results presented in this paper are more accurate than traditional GIE method.
However, there are several shortcomings in this paper. First, the major important factor in the susceptibility analysis and prediction of glacier debris flows is precipitation, such as average daily precipitation. We focused on changes in glaciers without considering precipitation. In this study, the volume of glacier melting over a period of time is selected instead of precipitation for the susceptibility of glacial debris flows. Second, the glaciers exhibit dynamic changes, i.e., melting while accumulating. The volume lost through glacier ablation does not completely turn into water, and may lead to the creation of new glaciers in the accumulation area, but for the prediction of geological disasters, we simulated the worst case in which all the melting of the glacier becomes the amount of water that triggers the occurrence of mudslides in the rainy season. Our presented work in this paper can be applied to glacier disasters in high-altitude areas.
In the future, we will consider some field campaigns in order to validate the data obtained. Our future work will focus on field survey records to improve the accuracy of our outcomes. The specific work is to design investigation and verification routes based on the topographical conditions, image interpretation, etc., we will discuss how to apply glacier changes as an indicator to other methods of debris flow susceptibility, such as machine learning methods. The changes in glaciers in high-mountain regions cause many geo-hazards. For example, floods are caused by the massive melting of glaciers and snow. Meltwater with a lot of sediment can block ditches and bridges, even silt roads. The weathered gravel begins to slide down along with the melting of the glacier covering its surface to trigger landslide activities. Therefore, our work will focus on how to use glaciology and geomorphology systematically to obtain more accurate prediction results of glacial disasters.

Conclusions
In this paper, by investigating the changes in glacier ablation and distribution based on remote sensing applications, we proposed analyzing the susceptibility for glacial debris flow in high-altitude regions. We presented a case study evaluating the susceptibility of G217 gullies in Tien Shan Mountains area. Through the comparison result of mapping susceptibility classification of the G217 gullies, the investigation results indicate that the improved GIE method with correction coefficient based on the glacier ablation volume is indeed better than the previous method, and combining glacier distribution and changes to the prediction of debris flow can improve the accuracy. In the future, we will develop machine learning method for analyzing debris flow susceptibility based on glacier changes, and explore a new prediction method of geo-hazards based on glaciology and geomorphology, especially that of glacier dynamics.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
The authors would like to thank the editor and the reviewers for their helpful comments.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: