Upward Expansion of Supra-Glacial Debris Cover in the Hunza Valley, Karakoram, During 1990 ∼ 2019

Introduction

Glaciers are not only a diminishing natural freshwater resource but also a sensitive indicator of global climate change (Yang, 1995; Kaab et al., 2012; Kraaijenbrink et al., 2017; Zemp et al., 2019). However, due to global warming, global glaciers are showing a trend of retreat and thinning. Recent studies on mass balance (Rowan et al., 2015; Lynch et al., 2016; Dehecq et al., 2018; Wu et al., 2018; Wang et al., 2019a; Wouters et al., 2019; Zemp et al., 2019), glacial area change (Paul et al., 2013; Patel et al., 2019; Reinthaler et al., 2019), surface velocity (Wang et al., 2018; Altena et al., 2019; Garg et al., 2019), glacio-hydrological modeling (Shrestha et al., 2015), and glaciers’ response to climate change (Scherler et al., 2011; Rowan et al., 2015) have revealed that changes in glaciers are exerting an impact on socioeconomic development in their downstream areas through changes in glacial water resources in High Mountain Asia (Immerzeel et al., 2010). Climate change may force geomorphological processes on high mountain slopes (Tipper et al., 2012; Cook et al., 2020) by accelerating the disintegration of rocks and increasing the accumulation of debris on glaciers and mountain slopes. As a result, the alpine glaciers widely distributed in the Pamir, Karakorum, Kunlun, Nyainqentangula, and Himalaya Mountains (Scherler et al., 2011; Khan et al., 2015; Zhang et al., 2016; Shukla and Garg, 2019) have widespread supra-glacial debris cover. Some studies have referred to glacier areas with sporadic debris cover as “dirty ice” (Robson et al., 2015; Fyffe et al., 2019a, b). The debris cover on glaciers affects the energy exchange at the ice surface and therefore has an impact on the surface mass balance of glaciers. It is therefore possible for glaciers to gain more heat from the environment, increasing the melt rate and causing the mass loss of glaciers.

Compared with clean ice or snow, the debris layer has a unique thermal process due to differences in physical properties such as reflectance, particle size, and color, which results in different ablation processes in the underlying ice (Østrem, 1959; Nicholson and Benn, 2006). Also, the differential ablation caused by an uneven distribution of debris thickness makes it easy for glaciers to form cliffs (Kindermann et al., 2008; Herreid and Pellicciotti, 2018) and ponds (Miles et al., 2016; Chand and Watanabe, 2019) in the ablation zone. Notably, these cliffs and ponds are not only factors that affect the hydrological process but also home to numerous glacial lakes, which can pose a serious threat to downstream communities and lead to catastrophic socioeconomic disasters in cases of glacial lake outburst flood (Benn et al., 2012; Dubey and Goyal, 2020). Therefore, obtaining information about the spatial distribution and temporal variation of supra-glacial debris cover would enhance the understanding of debris-covered glaciers and the glacial hydrological model. Some studies have shown that the area of debris cover is increasing as, overall, glaciers shrink and lose mass. For example, the supra-glacial debris cover in the Greater Caucasus increased from 48.3 ± 3.1 km² in 1986 to 79.0 ± 4.9 km² in 2014, based on Landsat and images from the years 1986, 2000, and 2014 (Tielidze et al., 2020). Furthermore, the glaciers in Pamir, Karakoram, and West Kunlun have moved forward or backward, either in equilibrium or showing an increase in mass (“Karakoram anomaly”) in recent years (Salerno et al., 2017; Farinotti et al., 2020; Gao et al., 2020), and it is of great interest to know whether this debris cover increase phenomenon exists for glaciers in anomaly areas, such as the Hunza Valley. Although some studies have mapped the spatial distribution of supra-glacial debris cover in this area (Khan et al., 2015, 2020; Mölg et al., 2018; Gao et al., 2020), they were limited to the extraction of spatial distribution information and lack any analysis on the temporal dynamic changes of supra-glacial debris.

Since the reflectance of supra-glacial debris is similar to that of non-glaciated slopes (Paul et al., 2004) and there is a lack of continuous, large-scale, high-quality cloud and shadow-free optical images, obtaining the long-term supra-glacial debris variations based on remote imagery is challenging. Previous studies have attempted to distinguish clean ice from debris-covered ice by using 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 (NIR)/short-wave infrared (SWIR)] or their combination from optical remote sensing images (Bolch et al., 2010; Alifu et al., 2015, 2016; Mölg et al., 2018). These methods can robustly delineate clean ice or snow, but they cannot accurately and automatically classify debris-covered 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 digital elevation models (DEMs) (Paul et al., 2004; Frey and Paul, 2012; Patel et al., 2019) 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; Yang et al., 2016; Lippl et al., 2018), and the recognition accuracy was improved. However, complex preprocessing and severe terrain noise from SAR data make large-scale applications difficult. In conclusion, especially for the long-time series involved in debris monitoring, optical images, with decades of continuous observations, such as the Landsat archive, are more practical.

In recent years, machine-learning-based classification methods have been applied to mapping glacier facies (Racoviteanu and Williams, 2012; Shukla and Yousuf, 2016; Zhang et al., 2019; Yousuf et al., 2020). Studies have shown that machine learning has advantages in extracting land-surface information from remote sensing images, which can effectively improve the accuracy of object recognition (Lary et al., 2017; Maxwell et al., 2018). However, mining large-scale and time series land information from high spatiotemporal resolution remote sensing data was found to be a computationally intensive task, requiring powerful computing platforms for analysis. Fortunately, some geospatial cloud-computing platforms are emerging that meet this demand, such as Google Earth Engine GEE, Amazon Web Services, Earth Server, and the Earth Observation Data Centre (Guo et al., 2020). Among these, GEE has 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 remove the parts of cloud, cloud shadow, and terrain shadow that affect each scene at the pixel scale, compared with the local image processing and analysis software, which takes the scene as a processing unit. The GEE platform has been widely used for various high-impact societal issues such as forest resources (Hazel et al., 2016), water resources (Pekel et al., 2016; Wang et al., 2019b), and land-use classification (Dong et al., 2016; Midekisa et al., 2017; Hao et al., 2019_ENREF_18) and has achieved remarkable results.

The aim of this study was to develop an automatic algorithm to identify debris-covered ice and map its spatiotemporal distribution in order to explore the dynamic process of supra-glacial debris cover by combining glacier inventory data and remotely sensed images on the GEE geospatial analysis platform. The study focused on glaciers in the Hunza Valley in the Karakorum Mountains of Pakistan. Otsu’s method was utilized to optimize thresholds of NDSI, and three machine learning algorithms [random forest (RF), support vector machine (SVM), and classification and regression tree (CART)] were used to classify supra-glacial features, including debris-covered ice, at the pixel level. Raw spectral information, band ratios, and color-to-grayscale conversion from Landsat 5/8 optical satellite imagery and the topographical components derived from DEM products were extracted as feature variables in the machine learning models. The same scheme was used to generate a time series of the debris-covered and clean ice areas in the study area. Finally, the results were comprehensively analyzed and discussed together with other data derived at the same time.

Study Area

The Hunza Valley is an area measuring ∼11,000 km², located in the western Karakoram, northern Pakistan (36°00′15′′ ∼ 37°05′23′′ N, 74°02′57′′ ∼ 75°46′48′′ E) (Figure 1). The topography across the Hunza Valley is characterized by large altitudinal variations, from 1,341 to 7,831 m above sea level (a.s.l.). The valley is home to glaciers with a total area of ∼3,600 km² (1,878 glaciers, with ∼12% of them being larger than 0.5 km²) that accounts for ∼33% of the basin area based on the Randolph Glacier Inventory (RGI) 6.0 dataset (RGI Consortium, 2017). Most of the glaciers (e.g., Hispar, Batura, and Barpu) are debris-covered and in a state of surging and advancing (Bhambri et al., 2017). Debris-covered glaciers are potential factors driving glacial lake outburst floods (Bhambri et al., 2019), which represent a major threat to local people, their property, and infrastructure such as the Karakoram Highway in the Hunza Valley. Climatologically, the study area is arid to semi-arid, situated in the subtropical climate zone and experiences significant variations in precipitation and temperature (Immerzeel et al., 2012). Based on the Moderate Resolution Imaging Spectroradiometer (MODIS) 1-km land-surface temperature daily products, the mean land-surface temperature for the entire region is −12.9°C in January and 20.1°C in July. Precipitation is mainly controlled by Indian monsoons and the westerlies, and the average annual precipitation is between 180 and 690 mm (Qureshi et al., 2017). Snow cover occupies approximately 80% of the basin’s land surface in the winter, decreasing to 30% in the summer (Tahir et al., 2011). The types of land cover in the basin include forest (0.4%), shrubland (16.1%), farmland (0.7%), and barren land (82.8%). The main soil types include Leptosols (type LP), rock outcrop soil (type RK), and glaciated soil (type GG). The primary soil component in the region is highly active clay, followed by rock outcrops and glacial soil (Garee et al., 2017; Ali et al., 2018).

Data and Methods

Landsat Imagery and Preprocessing

The Hunza Valley is covered by three WRS2 path/rows (149/034, 149/035, and 150/034) of Landsat images. We collected all the available standard level 1 terrain-corrected products of the Landsat raw scenes that were consistent in geometry and radiometry on the GEE platform, representing the years 1990–2019, including 10 Landsat TM images, four Landsat Enhanced Thematic Mapper Plus (ETM+) images, and 18 Landsat Operational Land Imager (OLI) images (see Table 1). The selection of images was constrained by the acquisition time during the ablation season (e.g., during July, or 200–270 days) when the images showed minimum snow cover, with little or no cloud cover. Mosaics of three path-row images within the same year in the research domain has been processed ready for debris extraction except for 2010 which were the mosaic from TM images acquired in 2008 (149/034), 2009 (149/035), and 2011 (150/034), respectively, due to wide snow cover of images in 2010. The Landsat 7 ETM+ data from 2000 only were used because of data gaps caused by scan line corrector failure since June 2003. Spectral bands with 30-m resolution were used in this study. All the data were publicly available, for free, on the United States Geological Survey (USGS) website.¹ We applied a standard top-of-atmosphere calibration, which is available on the GEE platform, to all the USGS Landsat Raw Scenes, converting the pixel digital number values to top-of-atmosphere reflectance (Chandera et al., 2009). Then, we assigned a cloud score to each pixel, using the Landsat simple cloud score algorithm,² which computes a simple cloud-likelihood score from 0 (not cloudy) to 100 (most cloudy), using a combination of brightness, temperature, and NDSI. We used a <10 threshold on the cloud score to mask cloudy pixels and took the per-band median values from the accepted pixels. Whole-image data preprocessing and the subsequent classification process were implemented by coding on the GEE platform.

Supra-Glacial Debris Extraction

Debris-Covered Ice Extraction

In this work, we used three algorithms to map the debris-covered outlines of glaciers. A detailed framework for debris-covered ice extraction is presented in Figure 2. The first algorithm was a machine-learning algorithm including SVM (Suykens and Vandewalle, 1999), RF (Liaw and Wiener, 2002), and CART (Breiman et al., 1984). Default parameters were used, but 500 trees were set for the RF classification model, and the kernel type of radial basis function that is suitable for the case of linear inseparability 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 (Band 1 ∼ 7); band ratios (NIR/SWIR1); NDVI ([NIR − Red]/[NIR + Red]); NDSI ([Green − SWIR1]/[Green + SWIR1]); NDWI ([Green − NIR]/[Green + NIR]); and luminance (0.3*Red + 0.59*Green + 0.11*Blue); and geomorphic parameters (slope and aspect). The training data were visually sampled, based on Landsat images at coincident times, combined with high-resolution images from Sentinel-2 and Google Maps. Samples in the regions of interest were divided into clean ice or snow, debris-covered ice (or supra-glacial debris), bare land, and other (e.g., vegetation, villages, rivers, lakes, and shadows), according to the land cover types pertaining to the Hunza Valley (Ali et al., 2018). For example, for 2019, 1,024 samples (shown in Figure 1) were selected, including 373 debris, 356 ice/snow, 270 bare land, and 205 other. The second method was a band ratio segmentation method (NDSI < 0.4) (Dozier, 1989) for eliminating the clean ice/snow part, under the restrictions of the RGI 6.0 outlines. However, a fixed NDSI threshold of 0.4 may not have been applicable to all periods. Therefore, we developed an optimization method based on Otsu algorithms to optimize the NDSI threshold to better distinguish between clean ice/snow and debris-covered ice in different periods. The 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). It is a global threshold method, and its principles are the following: assume that the gray value of an image is 1 ∼ N, divide it into two groups at value k, G₀ = [1∼k], and G₁ = [k + 1∼N], and calculate the probability of the two groups, ω₀ and ω₁, the average values for each group (μ₀ and μ₁), and the entire image (μ). Then, the variance of the two groups can be calculated by the following equation:

where σ²(k) is a threshold selection function. By changing the k-value in 1 ∼ N, the k-value at which σ²(k) is maximized is the required threshold.

Slope Delimitation and “Salt and Pepper Effect” Removal

Topographic parameters are a key factor in delineating glacier areas with debris cover. The spatial distribution of supra-glacial debris depends on the geomorphology and elevation gradient (Paul et al., 2004). Thus, parameters such as slope, aspect, and plan curvature derived from DEM were applied to improve the accuracy of classification of the debris cover and clean ice on the glacier surfaces. Some previous studies proposed various thresholds for the slope; for example, a slope < 24° (Paul et al., 2004) or smaller values [<12° (Alifu et al., 2015) or <14 ∼ 16° (Robson et al., 2015)] were used to distinguish debris-covered glaciers from the surrounding terrain. We hypothesized that the slope threshold would show spatial heterogeneity for various glacierized mountains. We used the slope gradient derived from the Shuttle Radar Topography Mission DEM void-filled version, known as “SRTM Plus,” at a resolution of 1 arc-second (approximately 30 m), which was available from the National Aeronautics Space Administration Jet Propulsion Laboratory. To obtain a typical slope threshold for this region, we evenly selected 713 slope samples from the glacier areas covered by debris in the Hunza Valley using visual interpretation. Our statistical results indicated that slopes of <25° dominated 99% of the debris-covered glacier areas with an average slope of 6.7° (Figures 3A,B). Therefore, a slope threshold of 25°, which is consistent with suggestions in Paul et al. (2004), was used to distinguish the glacier debris areas from the surrounding debris areas.

The “salt-and-pepper” effect is a common noise problem in pixel-based remotely sensed imagery classification, occurring when the same features (adjacent pixels) on the image are divided into different categories. Usually, median or morphological filtering is considered as a method for noise reduction (Serra and Vincent, 1992; Jassim, 2013). To remove this noise (also known as island pixels), we applied a morphological reducer filter to classified images using a custom kernel. We found that a 5 × 5 square kernel, with 50 connected pixels, was the best format for removing this effect, after a comparison of different kernel sizes (kernel radii of 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11, 13 × 13, 15 × 15, and 25 × 25) (Figure 3C) and types (square, circle, octagon, diamond, cross, and plus) and maximum number of connected pixels (5, 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, and 1000).

Definition of Glacier Outlines

At present, there are several datasets defining glacier outlines of the world including Global Land Ice Measurements from Space (GLIMS) (Raup et al., 2007), RGI (Pfeffer et al., 2014), Glacier Area Mapping for Discharge from the Asian Mountains (GAMDAM) (Nuimura et al., 2015), and the Second Glacier Inventory of China (Liu et al., 2015). In our study, the glacier outlines consisted of clean ice and supra-glacial debris outlines. We automatically extracted the clean ice outlines robustly using the classification algorithm, while relying on the existing RGI 6.0 outlines for the supra-glacial debris outlines, as in previous studies (Mölg et al., 2018; Scherler et al., 2018; Tielidze et al., 2020). For terminus-advancing or retreating glaciers, however, it is not reasonable to use a fixed glacier outline in time series of supra-glacial debris mapping. Therefore, we manually improved the RGI 6.0 outlines by visually interpreting glacier boundaries based on historical images of Landsat, Sentinel-2, and Google Maps, especially for the surging glaciers proposed by Bhambri et al. (2017) (Figure 1) and debris-covered glaciers. A total of six glacial outlines were modified based on the historical images of 1990, 1998, 2014, 2018, and 2019. We found that although the Karakoram contains a large number of surge-type glaciers (Quincey et al., 2011; Bhambri et al., 2017), there were no significant terminus-advancing or retreating glaciers in our study area. Thus, we assumed that the limiting effect of RGI outlines did not affect the dynamic analysis of supra-glacial debris. Nevertheless, there were voids or missing data in the glacier outlines mapped by RF modeling, possibly due to some darker pixels on the glaciers being misclassified, with so much attention being paid to these in the process during “salt and pepper effect” removal. Also, we found that the clean ice/snow automatically mapped by NDSI < 0.4 and Otsu thresholding contained some water pixels located outside the glacier area, so the results were masked using a 500-m buffer of the RGI 6.0 outlines.

Evaluation of Classification Accuracy

We used two methods to assess the accuracy of the classification results. For machine learning, a cross-validation method was used, where the total sample was divided into two parts, and 70% of the sample points from each class were randomly selected to train the model while 30% were withheld to form a validation dataset. Using the validation dataset, a confusion matrix was generated to assess the accuracy of the predictions across the classes, and the overall accuracy through the Kappa coefficient. The confusion matrix prediction versus the validation data based on the RF model and 2019 Landsat 8 OLI image are shown in Table 2. The prediction accuracy of the RF model for different land types exceeded 89%, especially for clean ice/snow with user and producer accuracy of 100%. Therefore, for this study, we used the RF model because it had the greatest accuracy compared with the SVM and CART models. Another evaluation tool is the “round robin” method proposed by Paul et al. (2017), which is based on multiple manual digitizations. Based on high-resolution Sentinel-2 and Google Earth images, we selected three glaciers (Kukki Jerab, Virjerab, and Yashkuk Yaz) and performed manual digitization five times, using the average value for evaluating the automatically derived extent and standard deviation for digitization accuracy. Table 3 shows the debris-covered area on three debris-covered glaciers derived by machine learning algorithms and by manual digitization.

Results

Changes in Supra-Glacial Debris Cover From 1990 to 2019

Overall, the supra-glacial debris-covered area in the Hunza Valley showed a trend of slowly increasing from 1990 to 2019. This increasing trend has been seen in other major glaciers during the same period, as indicated in Figure 4. 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.96 (0.94 ± 0.01), the overall accuracy was 95.5 ± 0.9, and the estimated debris coverage was 410.9 ± 25.9 km². The debris cover estimated by the Otsu thresholding method was 559.1 ± 23.4 km², with an optimized threshold of 0.53 ± 0.04. The NDSI < 0.4 method gives a value of 487.1 ± 16.9 km². The detailed estimates of supra-glacial debris cover are shown in Table 4 and Figure 5. The total area of supra-glacial debris cover in the Hunza Valley has expanded by about 8.1–21.3% from 1990 to 2019. In comparing the estimated values from the three methods, we found that those based on Otsu thresholding were the largest, followed by the NDSI < 0.4 method and the RF model, while the estimated standard deviation based on the RF model was the largest, with a value of 25.9 km². The results indicated that the supra-glacial debris cover in this region has been expanding slowly since the 1990s, due to an increasing number of clean ice areas being covered by impurities, such as rock and dust, which reduce the albedo, resulting in a debris-covered ice classification.

In order to further explore the variations in supra-glacial debris cover, we drew the spatial distribution of supra-glacial debris in the Hunza Valley from 1990 to 2019 using spatial overlay analysis, as shown in Figure 6A. The supra-glacial debris has migrated up-glacier, and the main area of change is located in the middle and upper regions of the glacier, close to the altitude of the equilibrium line (4,500–5,500 m a.s.l.). The variation in debris on the slopes showed a normally skewed distribution, with ∼60% of the variation in coverage being located in the area below a slope of 10° (Figure 6B). The maximum variation in supra-glacial debris was at an altitude of 5,000 m a.s.l. (Figure 6C) and with a northeast aspect. Also, We mapped the distribution of supra-glacial debris and clean ice/snow cover in the Hunza Valley at elevation gradient and aspect based on the estimation results by RF model on images acquired in 2019. 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 was 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 (Figure 7A). The median elevation of glaciers 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 supra-glacial 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 glaciers’ area, with only a few glaciers having a west-facing aspect (Figure 7B).

Status of Glaciers in the Hunza Valley

Over the past 30 years, the glacier area in the Hunza Valley has shown a downward trend. The glacier area estimated by the NDSI < 0.4, Otsu thresholding, and the RF model is 4,355.5 ± 210.2 km², 4,104.8 ± 167.6 km², and 3,810.8 ± 321.7 km², respectively. The clean ice/snow and supra-glacial debris areas of the Hunza Valley for the years 1990, 1998, 2000, 2010, 2013, 2014, 2015, 2017, 2018, and 2019 are included in Table 4 and Figure 5. The results presented in this study indicate a 9.311.6% decrease in clean ice/snow area in this region between 1990 and 2019. This decrease may be due to the retreat of glaciers and the overall expansion of debris-covered areas on glaciers (8.1–21.3%). According to the estimates from the RF model, the total glacier area of the Hunza Valley was 3,497.1 km² in 1990 and 3,286.6 km² in 2019. Overall, there was a 6.0% decrease in glacier area from 1990 to 2019.

Discussion

Uncertainty Analysis

Uncertainty, including that stemming from measurement errors, models, and scale effects, is an important step in validating the mapping results of supra-glacial debris. In this study, the quality of the satellite imagery was the most important factor, as this directly determined the classification results, especially in the high mountain areas with clouds and steep terrain. To exclude the impact of fresh snow, images from previous and subsequent years were used for some periods when there were few cloud-free images, such as in 2010, and this may inevitably have added to the uncertainty. Besides, the uncertainty mainly stemmed from the ground observational data (such as the counts and spatial distribution of the samples for the machine learning model), the classification methods, including the model types, the selected feature variables, model parameters or thresholds, and the postprocessing process. We analyzed this uncertainty from the following two aspects.

Uncertainties in Mapping Glacier Outlines

Uncertainty of supra-glacial debris extraction is a major source of uncertainty in mapping glacier outlines compared with the uncertainty from the extraction of the clean ice/snow part. The mapping uncertainty for the clean-ice areas was mainly affected by the seasonal snow cover. Seasonal snow cover has a greater impact on small glaciers, and the uncertainty may exceed 50% or even more (Mölg et al., 2018). From 2010 to 2014, the clean ice/snow area estimated by the three methods was on the high side, which may be caused by seasonal snow cover. In addition, we compared the glacier area estimated by Mölg et al. (2018), hereafter called the Mölg dataset, with our inventory to determine the main differences between them. The RF-based classification accuracy for clean ice/snow was 100% with an estimated glacier area of 3,810.8 ± 321.7 km², which is ∼6% larger than the Mölg dataset (3,617.1 km²). The glacier area estimated by NDSI < 0.4 and Otsu thresholding tended to be overestimated, being 738.4 km² and 487.8 km² higher than the Mölg dataset, respectively. Unfortunately, there were still voids or missing data in the glacier outlines mapped by the RF model, which also lead to the underestimation of glacier area in some years.

Uncertainties in Mapping Supra-Glacial Debris

Figure 8A shows that the results from the machine learning algorithms generally had high accuracy, with Kappa coefficients ranging from 0.82 to 0.95. In comparison, the classification accuracy for the SVM model was slightly lower than for the RF model, where the Kappa coefficient was 0.9 and the overall accuracy was 96.0%. The manual digitization for each glacier by five professionals showed decreasing standard deviations from the largest to the smallest glaciers, specifically Yashkuk Yaz (26.8 ± 2.2 km²), Kukki Jerab (15.2 ± 0.8 km²), and Virjerab (21.6 ± 0.4 km²) glaciers. The linear fit between the mean of the digitizations and the automatic estimates indicated that all the machine learning algorithms could produce results close to reality but that the RF model was better than the other two, in terms of the determination of an R² coefficient of 0.99 for the RF model, 0.87 for the NDSI <0.4 method and 0.86 for the Otsu thresholding algorithm (Figure 8B).

As shown in the results, there were differences in the supra-glacial debris cover area estimated by the three methods, with the Otsu thresholding tending to overestimate the debris coverage, while the RF model gave the smallest estimated value. To verify the reliability of these estimates, we compared the results with those from other works. For example, the supra-glacial debris area extracted from the 2007 and 2009 ALOS-1 PALSAR-1 coherence images in the Mölg dataset was 583.6 km², which was larger than the debris area estimated by our study for 2010–481.9 km² (NDSI < 0.4), 529.6 km² (Otsu thresholding) and 405.4 km² (RF model). Another comparative dataset is one obtained by extracting regions with NDSI < 0.4 and clipping them with RGI 6.0 outlines (Scherler et al., 2018), the same as method 2 of this study. In this work, the automatically extracted areas, using Landsat 8 composite images from 2013 to 2017 and Sentinel-2 composite images from 2015 to 2017, were 465.1 and 453.3 km², respectively. For the same period, the estimated results from the RF model in this study were comparable (408.9 ± 22.2 and 405.4 ± 13.2 km², respectively). However, the supra-glacial debris results for the NDSI < 0.4 and Otsu thresholding were higher than those of Scherler et al. (2018). It is difficult to evaluate the estimated results from this study in terms of previous results because they derive from different remote sensing data, and the methods used for data processing and classification are different; however, the supra-glacial debris-covered area falls within a reasonable range, in terms of magnitude.

Uncertainty from complex surfaces debris-covered glaciers, such as stagnant ice in glacier tongue areas, cliffs, and ponds on the surfaces of debris-covered areas, and boulders and dirty ice, have extensive effects on pixel-based classification. The supra-glacial ponds or lakes widely distributed on glaciers such as Batura, Hispar, Yashkuk Yaz, and Kukki Jerab frequently emerge or disappear, demonstrating obvious supra-glacial characteristics with time. We used salt-and-pepper-effect removal, as described above, to compare these effects. By comparison, we found that convolution kernel and maximum-connected pixel counts have to be properly decided in order to avoid the removal of correctly classified patches, or vice versa. In this study, a 5 × 5 kernel with 50 connected pixels was used, but this may not have been optimal. The problem of the salt-and-pepper effect is unavoidable in pixel-based classification. Although some studies have pointed out that the object-based classification method can effectively solve this problem and improve the accuracy (Rastner et al., 2014; Robson et al., 2015; Kraaijenbrink et al., 2016; Sahu and Gupta, 2018), its universality and operability on a large scale need to be further considered. Also, the stagnant ice at glacier tongues, such as the Batura and Hispar Glaciers, tend to have thick overlying debris layers with signs of vegetation (shrub) growth on their surfaces. Parts of these regions are classified as non-debris cover regions, which increases the uncertainty of estimating supra-glacial debris coverage, especially for RF models. Vezzola et al. (2016) showed that the number of debris-covered glaciers featuring supra-glacial trees is increasing globally as a response of the alpine environment to climate warming. They found supra-glacial trees present on the Miage Glacier when a debris thickness threshold (≥19 cm) was exceeded, and where there was a gentle slope (≤10°) and a low glacier surface velocity (≤7.0 m yr^–1) and where the vertical changes due to glacier dynamics were positive.

Additionally, although the image acquisition time was limited to 200–270 days, there may still be snow cover at high altitude, which is one of the significant uncertainties in mapping debris-covered ice. Due to the annual difference in temperature and precipitation, the melting limits of seasonal snow on the upper part of glaciers in summer are different each year, which may lead to an underestimation of the supra-glacial debris area. Similarly, the effects of shadows and cloud coverage inevitably pose challenges to the mapping of supra-glacial debris.

Characteristics and Possible Causes for Supra-Glacial Debris Cover Changes

Overall, the supra-glacial debris cover migrated up-glacier. Similar patterns of up-glacier migration have also been described for the Zmutt Glacier, Swiss Alps (Mölg et al., 2019) and the Greater Caucasus (Tielidze et al., 2020). Variations in supra-glacial debris cover were mainly observed in the middle and upper parts of the glacier, while the distribution of supra-glacial debris from the lower part to the glacier tongue remained stable. In the middle of the glacier, the area showing the most intensive variation was at the junction between the clean ice and the debris-covered ice. We found that the expansion of debris-cover over clean ice was slow. Due to the characteristics of glacier movements, the spatial distribution of supra-glacial debris cover in the downstream part of the glacier changes with time, while the area changes slightly. This redistribution of the surface moraine is a key factor in the diversification of glacier surface types, and it also changes the debris coverage. Recycled material from lateral moraines may be the source of this debris. As the glacier lowers and the permafrost melts, the lateral moraines become unstable, causing moraine material to fall onto the glacier (Nakawo et al., 1986). Debris eroded from the glacier bed may also be entangled in the ice. The degree of under-ice entrainment mainly depends on the thermal conditions at the base of the glacier and the substrate erodibility.

Previous studies have indicated that the main sources of debris are mass movements such as rockfalls, rock avalanches, debris-laden ice, and snow avalanches from the surrounding slopes (Hambrey et al., 2008). These sudden and sometimes massive debris migration/relocation events are common in alpine areas. For example, between September 1999 and June 2000, a landslide occurred in the upper region of the Batura Glacier, possibly as a result of slope failure. The accumulated rock fragments covered the glacier surface and moved downstream with the glacier ice, following a dynamic change in their shapes and sizes (Figure 9). Such phenomena play a key role in the increase and fluctuation of supra-glacial debris-covered areas.

Figure 10 shows the hypsometry of the supra-glacial debris coverage, and the five representative debris-covered glaciers, in 1990, 1998, 2010, 2014, and 2019. The zone at 4,500 ∼ 5,500 m a.s.l. displayed the greatest variation in supra-glacial debris (Figure 10A) and is within the altitude equilibrium line zone of glaciers in the Hunza Valley (Scherler et al., 2011; Kaab et al., 2012; Gardelle et al., 2013). These areas have steep topography and are prone to rockfalls and rock avalanches, in which there is an abundant release of supra-glacial debris. The amount of debris reaching a given glacier depends on the characteristics and extent of the catchment area, and especially its weathering and erosion rates (Haeberli et al., 2006). These are also affected by the lithology of the bedrock (Singh et al., 2011). Not all glaciers have experienced an increase in supra-glacial debris cover, with minor changes in supra-glacial debris having been recorded for the Batura, Parpik, and Gharesa glaciers (Figures 10B,C,E), which differ significantly from the obvious increases reported for the Kukki Jerab and Yashkuk Yaz glaciers (Figures 10D,F), especially between 1990 and 1998. In addition to subjective factors, the apparent increase in supra-glacial debris cover on Kukki Jerab and Yashkuk Yaz glaciers is mainly attributed to the evolution of ice ponds and cliffs and the upward expansion of surface moraines.

Impact of Climate Change

Thinning of glaciers and a warming atmosphere can lead to permafrost melting and slope instability at higher altitudes (Deline et al., 2015). Glacier changes in the Karakoram Mountains have been attributed to the dynamics of the Indian monsoon and the westerlies (Qureshi et al., 2017). To demonstrate the relationship between supra-glacial debris and climate change, we analyzed temperature, albedo, and precipitation datasets taken from MODIS products and Terra Climate Monthly³ and Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) products.⁴ The temporal changes in temperature, albedo, and precipitation derived from this data are shown in Figure 11. It is clear that warming and wetting have been a dominant phenomenon in the Hunza Valley in recent decades (Figures 11A,C,D). Temperature had a downward trend in 1990–1998; however, since 2000, the temperature derived from both MODIS and the Terra Climate products have shown a slight increase. Albedo on the ice surface is decreasing (Figure 11B)—-a sign of supra-glacial debris increase. The amount of absorbed solar radiation on a glacier’s surface increases where the ice is covered by supra-glacial debris or discontinuous debris (dirty ice); such debris causes an increase in glacier surface temperature and a decrease in albedo (Singh et al., 2011). As shown in Figures 11E,F, we also analyzed the relationships between temperature or precipitation and supra-glacial debris cover. Decreasing temperature and increasing precipitation have a positive relation with supra-glacial debris increase. This is particularly evident in the period 1990–1998. Global warming has accelerated the melting of glaciers resulting in changes in ice surface morphology and the redistribution of surface moraine material. A positive mass balance under decreasing temperatures and increasing precipitation may be caused by an increase in the production of material triggered by snow/ice avalanches or rockfalls.

Conclusion

The expansion of supra-glacial debris and decrease of glacier surface albedo have been reported widely in several parts of the word (Ming et al., 2012; Jiang et al., 2018; Fugazza et al., 2019; Tielidze et al., 2020). In this study, we mapped the dynamics of supra-glacial debris cover in the Hunza Valley, using 34 Landsat images (10 TM, four ETM+, and 18 OLI) acquired from 1990 to 2019. Firstly, an image composite method was applied, and sequential Landsat images with less cloud and snow cover at the pixel level were used to generate annual images of supra-glacial debris on glaciers. Then, the Otsu algorithm was utilized to optimize thresholds of NDSI segmentation, and the supra-glacial features (e.g., clean ice/snow and debris-covered ice) were classified on the state-of-the-art GEE cloud computing platform together with three machine learning algorithms (RF, SVM, and CART). All the training and validation datasets were derived from the visual inspection of freely available, high-spatial-resolution satellite imagery (Landsat, Sentinel-2, and Google Earth). Among these, the RF model produced the best classification accuracy, with a Kappa coefficient of 0.94 ± 0.01 and an overall accuracy of 95.6 ± 0.9%. Based on these estimates, we found that significant increases in areal supra-glacial debris cover had occurred during the last 30 years. The total area of supra-glacial debris cover in the Hunza Valley has expanded by about 8.1–21.3% from 1990 to 2019. It has migrated up-glacier, and the main area of change is located in the middle and upper regions of the glacier, close to the altitude of the equilibrium line (4,500–5,500 m a.s.l.). We also found that the interannual temperature decrease and precipitation increase had positive relations with the supra-glacial debris increase. This supra-glacial debris areal change information can also be used in mass balance, glacier hydrology, glacier hazard, and glacier response to climate change models.

The identification of debris is a fundamental, and yet still challenging field in studying glacier change and water resources. With the effects of global warming, clean glaciers in the high mountains of Asia are increasingly changing into debris-covered glaciers, which leads to changed discharge in rivers mainly supplied by glacial meltwater, resulting in the water resources in basins, and their future trends, likely to be greatly affected. To explore the response of debris-covered glaciers to climate change, and their effect on hydrology and water resources, future work should focus on extracting information about debris thickness and other surface features (e.g., glacial ponds and cliffs) using high-resolution aerial or satellite imagery and deep-learning techniques.

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Author Contributions

The study was designed by SL and FX. FX ran models and performed the data analysis. FX and SL performed the analysis and drafted the first manuscript. KW, YZ, YG, MQ, SD, MS, and AT discussed and improved the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This project has received funding from the National Natural Science Foundation of China (Grant Nos. 41761144075 and 41801031), the Second Tibetan Plateau Scientific Expedition and Research Program (STEP, Grant No. 2019QZKK0208), the Research Fund for Introducing Talents of Yunnan University (No. YJRC3201702), and the initiative for the promotion of Yunnan University in the development of a world-class university. FX was also supported by the Innovation Fund Designated for Graduate Students of Yunnan University (Grant No. 2019226).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the Google Earth Engine Science team for the freely available cloud-computing platform and USGS for Landsat imagery and SRTM DEM. We thank Sidou Zhang, Xinxin Qiang, and Ying Yi at the Institute of International Rivers and Eco-Security, Yunnan University, who involved in the conception of the manuscript. Finally, we thank the editor TZ and three reviewers for their constructive comments on the manuscript.

Footnotes

1. ^ https://www.usgs.gov/land-resources/nli/landsat
2. ^ https://developers.google.com/earth-engine/landsat
3. ^ http://www.climatologylab.org
4. ^ https://chc.ucsb.edu/data/chirps

References