DERIVATION OF SUPRAGLACIAL DEBRIS COVER BY MACHINE LEARNING ALGORITHMS ON THE GEE PLATFORM: A CASE STUDY OF GLACIERS IN THE HUNZA VALLEY

Calculating the spatial-temporal distribution of supraglacial debris cover on glaciers is essential for understanding mass balance processes, glacier lake outburst floods, hydrological predictions, and glacier fluctuations that have attracted attention in recent years. However, due to the reflectance of supraglacial debris is similar to that of non-glacier slopes, mapping supraglacial debris cover based on optical remote sensing remains challenging. In this paper, we used NDSI and machine learning algorithm to extract debris cover on glaciers in Hunza Valley, Pakistan. Our result showed that the RF model has the best classification accuracy with kappa coefficient of 0.94 and overall accuracy of 96%. The debris-covered area increased by 21.31% from 1990 to 2019 (394.76 km ~ 478.88 km) in the study area. Results and the method are of significance in the assessment of meltwater modeling for glaciers with debris cover.  Corresponding author

In the past, attempts at boundary identification of glaciers with debris cover were made with the assistance of individual parameters such as the normalized difference vegetation index (NDVI), normalized difference snow index (NDSI), normalized difference water index (NDWI), and spectral band ratio thresholds (e.g., near-infrared/short-wave infrared (NIR/SWIR)) or their combination from optical remote sensing images (Alifu et al., 2016;Bolch et al., 2010;Mölg et al., 2018). The aim was to distinguish between clean ice and debris-covered ice. These methods can robustly delineate clean ice or snow, but they cannot accurately and automatically classify debris-covered glacier ice as distinct from clean ice and the surrounding land surface (Robson et al., 2015). This has stimulated studies on the use of other parameters such as geomorphic parameters derived from a digital elevation model (DEM) (Patel et al., 2019;Paul et al., 2004) and thermal characteristics from the infrared band (Singh and Goyal, 2018), as well as utilizing the coherence change between two successive synthetic aperture radar (SAR) images (Janke et al., 2015;Robson et al., 2015). However, complex pre-processing and severe terrain noise of SAR data make large-scale applications difficult.
In recent years, machine-learning-based classification methods have been applied to identify glacial outlines (Racoviteanu and Williams, 2012;Zhang et al., 2019). Machine learning has advantages for extracting land surface information from remote sensing images (Maxwell et al., 2018) because mining largescale and time-series land information from massive remote sensing data is a computationally intensive task and requires powerful computing platforms for analysis. Fortunately, some geospatial cloud-computing platforms are emerging to meet this demand, such as Google Earth Engine (GEE), Amazon Web Services (AWS), Earth Server (ES), and Earth Observation Data Centre (EODC) (Guo et al., 2020). Among these, GEE has obvious advantages because it is an open-source cloud-based platform for planetary-scale geospatial analysis that integrates mainstream free satellite data such as the Landsat archive, Sentinel series imagery, and other terrain products and climate data (Gorelick et al., 2017). It can quickly pull information from massive satellite-image data on various high-impact societal issues such as forest resources (Hazel et al., 2016), water resources (Pekel et al., 2016), land use classification (Dong et al., 2016), and other fields.
The objective of this study is to develop an automatic algorithm for identifying debris-covered ice and mapping its spatiotemporal distribution by combining glacier inventory data and remotely sensed images based on the GEE geospatial analysis platform. The study is focused on glaciers in the Hunza Valley in the Karakorum Mountains of Pakistan. Based on the traditional knowledge-based approach, Otsu's method was utilized to optimize thresholds of band-based indices and machine learning algorithms such as random forest (RF), support vector machine (SVM), and classification and regression tree (CART) that are used to classify supraglacial features, including debris-covered glacier ice, based on the spatial resolution of the images used. Raw spectral information, band ratios, and color-to-grayscale conversion from Landsat 5/8 optical satellite imagery and the topographical components derived from SRTM DEM products were extracted as feature variables in the machine learning model. The same scheme was used to generate a time-series of the debris-covered area and bare ice area in the study area. Finally, the results were comprehensively analyzed together with other data derived at the same time.

STUDY AREA
The Hunza Valley is an area measuring ~11000 km 2 located in western Karakoram, Northern Pakistan (latitude 36°00'15"~ 37°05'23" N, longitude 74°02'57" ~ 75°46'48" E) ( Figure 1). The topography across Hunza Valley is characterized by large altitudinal variations from 1,341 m to 7,831 m above sea level(a.s.l.). The valley is home to glaciers with a total area of ~3600 km 2 that accounts for ~33% of the basin area based on the RGI 6.0 dataset (RGI Consortium, 2017). Most of the glaciers (e.g., Hispar, Batura, Barpu) are debris-covered and in a state of surging and advancing (Bhambri et al., 2017). Debriscovered glaciers are potential factors driving glacial lake outburst floods (GLOFs) (Bhambri et al., 2019), which represent a major threat to local people, their properties, and infrastructure such as Karakoram Highway (KKH) in the Hunza Valley. Climatologically, the study area is arid to semi-arid and lies in the subtropical climate zone, with significant variations in precipitation and temperature (Immerzeel et al., 2012). Based on the MODIS 1 km LST Daily products, the mean land surface temperature (LST) for the entire region is -12.9℃ in January and 20.1℃ in July. Precipitation is mainly controlled by Indian monsoons and the westerlies, and average annual precipitation is between 180 mm and 690 mm (Qureshi et al., 2017). Snow cover occupies approximately 80% of the basin's land surface in the winter and decreases to 30% in the summer (Tahir et al., 2011). Types of land cover in the basin include forest (0.36%), shrubs (16.12%), farmland (0.76%), and barren land (82.78%). The main soil types include Leptosols (LP), rock outcrop soil (RK), and glaciated soil (GG); further, the major soil component in the region is highly active clay, followed by rock outcrops and glacial soil (Ali et al., 2018;Garee et al., 2017). Surging glaciers are marked with triangles (Bhambri et al., 2017).

Data and Processing
We used cloud-free imageries of Landsat 5 Thematic Mapper (TM) and Landsat 8 Operational Land Imager (OLI) acquired in each ablation season and are available on GEE. We used Sentinel-2A/B Multi-Spectral Instrument images with a spatial resolution of 10m and a recycle time of 5 days for validation purpose. The 1 arc-second SRTM DEM V3 was used to derive topographical parameters and elevation gradient, like slope, aspect, and plan curvature. RGI v6.0 glacier inventory for the region was used to constrain the spatial extend of model runs for debris cover extractions which have been modified based on Sentinel-2 images. The processing procedure can be summarised as follows: the selection of images excluding seasonal snow cover was done by confining images acquired during the ablation season (e.g., during July, or days 200 to 270 ± 10 days). Then, the pixel digital number (DN) values were converted to atmospheric top-level reflectance (TOA) (Chander et al., 2009). Third, images (pixels) with the least clouds were selected by applying a simple Landsat cloud score algorithm available on GEE, generating a composite image with no-cloud pixels of scenes in the region, and computing per-band percentile values (25%) from the accepted pixels. Whole image data pre-processing and the subsequent classification process were implemented by coding on the GEE platform.

Debris-covered ice extraction:
In this work, we used three algorithms to delineate boundaries for glaciers with debris cover. The first method is a multi-criterion-based approach that uses thresholds of NDSI < 0.4 (Dozier, 1989) to distinguish clean ice/snow from debris-covered ice and NDVI < 0.1 to exclude non-glaciers. However, a fixed NDSI threshold of 0.4 may not be applicable for all regions. We then developed an optimization method based on Otsu algorithms (Otsu) to optimize the NDSI threshold to better distinguish clean ice/snow and debris-covered ice in different areas. An Otsu algorithm is an automatic non-parametric and unsupervised method for thresholding that is used to automatically detect targets in computer vision and image processing fields (Ng, 2006). Otsu is a global threshold method, and its principles are the following: Assume that the grey value of an image is 1 ~ N, divide it into two groups at k value, G 0 = [1~k] and G 1 = [k + 1~N], and calculate the probability of the two groups, ω 0 and ω 1 , the average values for each group (μ 0 and μ 1 ), and the entire image (μ). Then, the variance of the two groups can be calculated by the following equation: 2 (k) = ω 0 (μ 0 − μ) 2 + ω 1 (μ 1 − μ) 2 where, 2 (k) is a threshold selection function. By changing the k value within 1 ~ N, the k value at which 2 (k) is maximized is the required threshold.
The third algorithm is a machine-learning algorithm that includes SVM (Suykens and Vandewalle, 1999), RF (Liaw and Wiener, 2002), and CART (Breiman et al., 1984). In this study, default parameters were used, but 500 trees were set for the RF classification, and the kernel function of the radial basis function (RBF) was applied in the SVM model. It is obvious that a single spectrum cannot fully solve the problem of the similarity of ice covered with debris to the surrounding terrain. We generated 14 feature variables: original spectrum (b1, b2, b3, b4, b5, b6, b7), band ratios (nir/swir1), NDVI [(nirred)/(nir+red)], NDSI [(green-swir1)/(green+swir1)], NDWI [(green-nir)/ (green+nir)], and luminance and geomorphic parameters (slope and aspect). Training data were visually sampled based on high-spatial-resolution Sentinel-2 and Google Earth images. Samples in the regions of interest (ROIs) were divided into clean ice or snow (ice/snow), debris-covered ice (debris), bare land, and others (e.g., vegetation, villages, rivers and lakes, and shadows) according to land cover types in the Hunza Valley (Ali et al., 2018). For example, for 2019, 1,024 sample points were selected, including 373 debris, 356 ice/snow, 270 bare land, and 205 other points. The spatial distribution of ROIs for 2019 is shown in Figure 1B.

Post-classification processing:
The classified results were processed by applying slope threshold and RGI6.0 + outlines and removing the 'salt-and-pepper effect' and glacier ice and debris-covered areas measuring less than 0.05 km 2 . (i) Slope threshold This is a key parameter for delineating glacier areas with debris cover. Some early studies proposed various thresholds for slope. For example, a slope < 24° can be used to distinguish debris-covered glaciers and the surrounding terrain (Paul et al., 2004). Some proposed smaller values such as a slope < 12° (Alifu et al., 2015) or a slope < 14~16° (Robson et al., 2015). We hypothesize that the slope threshold shows spatial heterogeneity for various glacierised mountains.
To obtain a typical value for this region, we selected 713 sampling points uniformly for glacier and non-glacier debris areas of the Hunza Valley. Our statistical results show that slopes of less than 25° dominate the glacier area (99% of the glacier debris area) with a mean slope of 6.71° (Figure 2A, B). Therefore, a slope threshold of 25°, which is consistent with Paul's suggestion, was used to confine the glacier debris area from the surrounding debris area (Paul et al., 2004).

Evaluation of classification accuracy:
To evaluate the accuracy of classification models and the authenticity of the estimation results, we used two methods. For machine learning, a cross-validation method was used by dividing the total sample into two parts, with 70% of the sample points from each class randomly selected to train the model and 30% of sample points withheld as a validation dataset. Using the validation data, a confusion matrix was generated to assess the accuracy of predictions across class and overall accuracy through the Kappa coefficient and overall accuracy (OA). Another accuracy evaluation method using the 'round robin' method is based on multiple manual digitizations, as proposed by Paul and colleagues (Paul et al., 2017). Based on high-resolution Sentinel-2 and Google Earth images, we selected three glaciers (Kukki Jerab, Virjerab, and Yashkuk Yaz) and performed manual digital digitization five times and used the mean value for evaluating automatically derived extent and standard deviation for digitization accuracy. Figure 3A shows that results from the machine learning algorithms generally have high accuracy, with Kappa coefficients ranging from 0.82 to 0.95. In comparison, classification accuracy for the SVM model is slightly lower than for the RF model, whose Kappa coefficient is 0.9 and overall accuracy is 96.02%. Table 1 shows the debris-covered area on three debris-covered glaciers derived by machine learning algorithms and by manual digitization. Manual digitization for each glacier by five professionals shows decreased standard deviations for the largest to the small glaciers, specifically Yashkuk Yaz (26.77 ± 2.19 km 2 ), Kukki Jerab (15.19 ± 0.80 km 2 ), and Virjerab (21.56 ± 0.35 km 2 ) glaciers. The linear fit between the mean of the digitizations and automatic estimates indicates that all the machine learning algorithms can produce results close to reality, but the RF model is slightly better than the other two with reference to the determination R 2 coefficient of 0.98 for RF, 0.97 for the multi-index method, and 0.93 for the Otsu-based model ( Figure 3B). Table 2 shows the confusion matrix of model prediction versus validation data using the RF model and 2019 remote sensing data. The accuracy of RF model prediction for different land types is more than 90%, especially for clean ice/snow with user and producer accuracy of 100%. Considering the accuracy of the RF model for satellite data acquired at different times, the Kappa coefficient of the RF model was within a range of 0.92-0.95 (Figure 4). In addition, to verify the accuracy and reliability of RF model estimation, we also compared the results with those of others ( Figure 4). For example, the manual digitization results of Mölg and colleagues (Mölg et al., 2018) are an area of 583.56 km 2 , which is much larger than the debris cover area of this study in 2010 (440.10 km 2 ). This is consistent with the overestimation described by the authors. Another result automatically extracted by Scherler and colleagues (Scherler et al., 2018) is roughly consistent with the results from this study, with a residual error of 25~33 km 2 .

Mapping Supraglacial Debris Cover
The glacial debris cover area for six periods (1990, 2000, 2010, 2013, 2016, and 2019) in the Hunza Valley was delineated based on the RF model. The total area of glaciers with debris cover in the study area showed a clear increase from 1990 (394.76 km 2 ) to 2019 (478.88 km 2 ) ( Figure 4). This increasing trend has been seen for other major glaciers during the same period as shown in Figure 5. This means glaciers in the region have been becoming dirty since 1990.
The results for bare ice and supraglacial debris cover area are based on the RF model (shown in Figure 6A) and on images acquired in 2019 for all glaciers in the Hunza Valley. More than 78% of the clean ice/snow area lies higher than 5,500 m a.s.l., while about 80% of the debris-covered ice is distributed between 4,000 and 5,000 m a.s.l.. The median elevations were approximately 5,365 m for clean ice/snow and 4,075 m for debris-covered ice. The median elevation of glaciers in Hunza is 5,230 m a.s.l., which is sometimes referred to as the equilibrium-line altitude of glaciers in the Hunza basin (Qureshi et al., 2017). This means that supraglacial debris is extensively distributed in the lower part of the ablation area of glaciers. Most glaciers have a North (N) and Northeast (NE) aspect, accounting for 38.7% of the glacier area, and only a few glaciers have a West (W) aspect ( Figure 6B).
Overall, the debris-covered area of each glacier has also increased over time but varies by altitude gradient. In the lowaltitude area (i.e., < 3,200 m a.s.l.), the debris-covered area steadily changes without a significant increase, with an average variance of about 0.1 km 2 but increases significantly in the midaltitude range (3,200 ~ 4,200 m a.s.l.). When the elevation is higher than 4,200 m a.s.l., the debris-covered area increases with large variance, especially in areas above 5,000 m. Seasonal snow is a key factor that causes large variance in the debris cover of these areas. The essence of debris expansion is the conversion of clean ice into impure ice, which changes its reflectance on remote sensing images, causing it to be identified as the debris class. The main material source of these moraines is impure ice formed by dust accumulation. Moreover, due to deglaciation, rocks in the ice body are exposed and accumulate, the bedrock on both sides of a glacier disintegrates under the influence of external forces, or debris such as rolled clastic rock accumulates on the ice surface because of climate change.

Uncertainties in Supraglacial Debris Mapping
Uncertainty analysis, including uncertainty from measurement errors, models, and scale effects, is an important step in validating the authenticity of remotely sensed products (Wu et al., 2014). Aiming to determine the extent of debris cover using remote sensing data, we analyzed sources of uncertainty in this paper. This mainly includes the following aspects: i) Errors from remote sensing data, including distortions generated by the sensor during image acquisition, and image quality and mixed pixel effects caused by spatial resolution. Satellite data directly determine the quality of the classification results, especially in high mountain areas, which are affected by clouds and steep terrain. In this paper, to exclude the impact of fresh snow, we selected images from the melting season with heavy cloud cover. This resulted in a lack of high-quality images for some years, requiring us to use images from previous and subsequent years. ii)Errors from ground observation data, such as the counts and spatial distribution of the sample points in this paper for machine learning classification. iii) Errors from classification methods, mainly including model types, selected feature variables, value selection for model parameters or thresholds, and accuracy evaluation method. iv) The complex surface types of debris-covered glaciers, such as stagnant ice in glacier tongue areas, cliffs and ponds on debris-covered area, large moraines, and dirty ice.
Generally, the methods for delineating debris include manual digitization and multi-rule-based and supervised classification. From the nature of different algorithms, manual digitization is time-consuming, and the interpretation accuracy varies depending on the experience of the interpreter. The multi-rulebased method can quickly produce binarised results, with simple post-classification processing. However, the rule-based method generally determines the threshold based on experience, and this threshold may not be suitable for all glacier areas. Fortunately, the Otsu-based method developed in this paper addresses this problem and automatically calculates the threshold, which is more flexible than the empirical threshold. Supervised classification is complicated and requires a huge amount of calculations, but when it is combined with geospatial big data analysis platforms such as GEE, it has great potential for classification. From the perspective of the classification effect, all methods are highly robust for identifying clean ice or snow, but there is a slight difference in the identification of the debris-covered layer, which is mainly reflected in the glacier tongue and the interface area between clean ice and debriscovered ice. All methods have inaccurate identification at the glacier tongue with thicker debris or vegetation (Vezzola et al., 2016). At the interface of clean ice and debris-covered area, the effects of machine learning classification are the best, followed by the Otsu-based method. On the whole, the Otsu-based method is fast and can combine RGI data with appropriate manual corrections, making it highly suitable for large-scale glacier inventory. Further, because the pixel-based classification considers only the spectral characteristics of a single pixel, it is easy to identify the ponds and cliffs on debris-covered glaciers as water or clean ice. In addition, the results with salt-andpepper effects often require post-classification processing, especially for machine learning algorithms. Therefore, objectbased classification based on image features such as shape or topological features has been developed with higher accuracy than pixel-based methods (Kraaijenbrink et al., 2016;Robson et al., 2015). For example, Robson and colleagues (Robson et al., 2015) showed that compared with the pixel-based method, an object-based method can improve the classification of debriscovered glaciers by 10%.

New Approaches for Future Work
The biggest challenges in debris-covered ice delineation are the complex surface features of glacial tongue areas (i.e., stagnant ice) and mixed pixels caused by the coexistence of debris and clean ice on upper glacier areas. Due to high-resolution optical remote sensing, images have more detailed characteristics, making the first approach conducive to high-precision detection of targets. Extraction of multiple feature variables, such as texture features, combined with the powerful capabilities of deep learning, is a method capable of extracting these complex surface features.
Another approach makes full use of the motion characteristics of glaciers to identify debris-covered ice. Ice flows, which are a major feature distinguishing glaciers from other natural bodies, control the basic process of glacial change. Observations of glacier ice flow are of great significance to reveal glaciers' changing rules and predict their future changes. Ice flow can be monitored by remote sensing methods such as visual tracking, cross-correlation (Fahnestock et al., 2016), and interferometric radar (Paul et al., 2013). Previous studies have used glacier velocity data to analyze glacial motion mechanisms, test glacial flow theories, and identify surging or advancing glaciers (Dehecq et al., 2018;Quincey et al., 2011), but no study has used glacier mapping. Therefore, to verify the possibility of using glacial surface motion characteristics to distinguish glaciers and non-glacier areas, we conducted a test using ITS_LIVE surface velocity products, which are derived by cross-correlation feature tracking (Gardner et al., 2019). The surface velocity is generated from Landsat images, and the resampling spatial resolution is 240 m. Figure 7A shows the mean surface velocity from 2011 to 2018. The surface velocity in non-glacial areas is between 0 ~ 1 m/yr. As shown in Figure  7B and 7C, we extracted areas with a surface velocity greater than 5 m/yr and 10 m/yr and overlaid the debris cover estimated in this paper. This showed good edge matches at coarse resolution. After improving the spatial and temporal resolution of surface velocity, this motion feature can solve the problem of mixed pixels between glaciers and lateral moraines, especially for large debris-covered glaciers.

CONCLUSIONS AND OUTLOOK
The purpose of this study was to develop an automatic method for debris-covered ice mapping and further explore the change in the area of debris cover of glaciers over time. We coded three algorithms for debris cover extent mapping based on the GEE geospatial big data analysis platform, and we tested them with data on glaciers in the Hunza Valley using Landsat satellite imagery. The main conclusions are i) The RF has the best accuracy among the classification algorithms, with a Kappa coefficient of 0.94, an overall accuracy of 96.02%, and determination coefficient (R 2 ) of the linear fit between the manual digitisation of 0.98. However, after analysing the performance and practicability of the methods, we found that the Otsu-based method has great potential for large-scale glacier inventory when combined with manual corrections. ii) The supraglacial debris cover of the Hunza Valley in 1990Valley in , 2000Valley in , 2010Valley in , 2013Valley in , 2016Valley in , and 2019 was estimated based on RF classifier and Landsat images. The results show that the total debris cover area was 394.76 km 2 in 1990 and 478.88 km 2 in 2019, indicating that the debris cover on glaciers has increased over the past 30 years.
Identification of debris is a fundamental but challenging field for researching glacier change and water resources. Debris affects the melting rate of glaciers, with thin debris accelerating melting and thick debris slowing melting. Due to the heterogeneity of debris thickness, debris-covered glaciers experience differential ablation, which makes the surface moraine gather or disperse in space. Such gathering or dispersing, in turn, accelerates differential ablation. This interaction results in the formation of supraglacial ponds and ice cliffs, and these affect glacier-derived runoff to a great extent. Moreover, these cliffs and ponds are not only factors that influence the hydrological process but also sources of glacial disasters such as GLOFs, which affect production and living safety in the downstream areas. With the effects of global warming, clean glaciers in the high mountains of Asia are increasingly changing into debris-covered glaciers, which leads to a change of discharge in rivers mainly supplied by glacial meltwater. Consequently, the water resources in basins and future trends will be greatly affected. Glacier meltwater runoff is an important source of freshwater for human production and life. Research by Immerzeel et al. shows that owing to the retreat of mountain glaciers, reduction of snow cover, and widespread shrinkage and disappearance of alpine lakes, a quarter of the global population is at risk of more severe water supply problems (Immerzeel et al., 2019). Overall, studies on the identification of surface features of debris-covered glaciers (e.g., debris thickness, glacial ponds, and cliffs) and their spatial distribution and seasonal variation are meant to indicate the response of debris-covered glaciers to climate change and the effect on hydrology and water resources.