Assessing the supraglacier extent of Karakoram glaciers located in Hunza River Basin, Pakistan

B

band

CPs

change points

Glaciers are a significant source of fresh water, a sensitive indicator of climate change, and are retreating and thinning because of global warming (Kraaijenbrink et al. 2017; Tielidze et al. 2020). Recent studies on glacio-hydrological modeling, glacial area change, and mass balance response to climate change (Wu et al. 2018; Wang et al. 2019; Wouters et al. 2019) highlight impacts on the socioeconomic development of downstream riparian (Immerzeel et al. 2010). However, debris accumulation on mountainous slopes and glaciers as a result of climate-induced rock disintegration has been found at Karakoram Highway (KKH) (Khan et al. 2015; Cook et al. 2020) that affects the heat exchange at the surface affecting the mass balance of glaciers in terms of melt rate and mass loss. Moreover, the thickness of the debris covers slows down the melting, thus delaying its response to climate change. Therefore, the aerial assessment of supraglacial debris cover (SDC) is important for studies related to hydrologic modeling and glacier dynamics.

Assessment of debris cover has a unique thermal characteristic because of physical attributes including color, reflectance, and particle size which result in the melting process beneath the ice (Nicholson & Benn 2006). The similar reflectance of supraglacial debris with non-glaciated slopes and the discontinuous repository of remote sensing images without clouds and shadows make it challenging to track long-term variation in SDC (Paul et al. 2004). Previous research adopted various methods to distinguish the ice from debris cover by indices such as the normalized difference water index (NDWI), the normalized difference vegetation index (NDVI), the normalized difference snow index (NDSI), the spectral band ratio thresholds (NIR/SWIR), or a combination of these (Alifu et al. 2016; Mölg et al. 2018). Since these methods cannot accurately classify the debris cover (Robson et al. 2015) and urge to use of other parameters related to digital elevation models (DEMs) (Frey & Paul 2012; Patel et al. 2019), thermal characteristics from the infrared band to improve the accuracy of debris cover classification (LIPPL et al. 2018; Singh & Goyal 2018). Moreover, automatic algorithms based on machine learning approaches including support vector machine (SVM), random forest (RF), and classification and regression tree (CART) with remote sensing images and multiple DEM data are used by various researchers to map glacier facies (Racoviteanu & Williams 2012; Shukla & Yousuf 2017; Zhang et al. 2019; Yousuf et al. 2020) and debris-covered glaciers (Zhang et al. 2019; Xie et al. 2020). However, powerful computational systems are required to extract the large scale and time series of land information from high spatiotemporal resolution images and debris cover glacier mapping based on near-infrared (NIR) or visible (VIS) bands is challenging because of low detection efficiency (Bhambri et al. 2011; Frey et al. 2012). Therefore, field based survey coupled with remote sensing techniques is significant for detecting the debris-covered glacier (Vincent et al. 2013).

Satellite images are an important source for measuring the area of glaciers extent. However, field measurements are strongly recommended, e.g. the extraction of debris thicknesses to see their response to hydrology and water resources (Xie et al. 2020) and the selection of training samples for the validation of debris-covered glaciers (Zhang et al. 2019). The extracted glaciated extent area using the remote sensing techniques is different for different satellite products and often leads to erroneous results when used in climate models for prediction of climate change impact on glaciers and stream flow estimation which emphasizes the need for bias correction in the satellite products (Kääb et al. 2014). Moreover, the classification of snow on glaciers is difficult with band ratios, as the shape of the spectral curves of ice and snow is very similar and thus results in the same values for both facies. A robust and widely used approach for mapping SC on glaciers is based on the use of a threshold on either single-band reflectance values, preferably in the NIR to avoid saturated values over snow with Landsat data (Veettil & Kamp 2017), or albedo products that are readily available (Brun et al. 2015; Stephens et al. 2015). In previous studies, different researchers have conducted outlines of glacier data using different satellites, and agreeing on accuracy is a great challenge (Albert 2002; Friedl et al. 2020). However, the measured quantities of glaciers are not commonly verified by ground truthing (Eisen et al. 2008; Friedl et al. 2020), which results in insignificant and low-precision data which leads to overestimation or underestimation of the flows of glaciated rivers (Shea & Immerzeel 2016). Hence, there is a need to estimate glacier extent using remote sensing and GIS followed by ground verification of SDCs (Vincent et al. 2013).

The overall aim of this study is to assess the selected glacier coverage through the well-proven application of NDSI using multiple satellite images (Landsat-8, Sentinel-2, and MODIS) and to outline the SDCs through field survey. The effectiveness of NDSI, derived from the visible and NIR spectral bands, has been widely documented for its ability to isolate snow and ice from other land cover types (Shukla et al. 2010). This study further compares the drawn SDC boundary with the GLIMS/RGI glacier outline inventory for its future applications regarding mass balance and hydrological modeling studies.

The Hunza is one of the main sub-basins of the Upper Indus Basin (UIB), and it contributes about 12% of the total flow of the UIB upstream of the Tarbela reservoir (Shrestha et al. 2015). The area of the Hunza catchment is about 13,733 km² and is covered by 25% glaciers. The basin hosts some extensive glacier systems, including Pissan (14.411 km²), Gulmit (14.621 km²), Gulkin (30.544 km²), Pasu (62.194 km²), Khurdopin (111 km²), Virjerab (112 km²), Batura (238 km²), Hispar (339 km²), and a few others (Arendt et al. 2017; Nazeer et al. 2022). The majority of the glaciers are valley and debris covered (Bhambri et al. 2017) and along the KKH. Glaciers in the area feed the Hunza River, as the main source of water for the area. The selected glaciers and others have retreated, as documented after evaluating the temporal glacier inventories by Anwar & Iqbal (2018), Shafique et al. (2018), and Zhang et al. (2019). Moreover, the climatic data specifically temperature is rapidly increasing resulting in the melting of glaciers that can cause the shrinkage of fresh water and glacier lake outburst floods (GLOFs) that impose a threat to downstream living communities and major infrastructure such as the KKH. Since these glaciers serve as essential water supplies for towns downstream, underscoring the urgent need for accurate hydrological modeling and monitoring that combines state-of-the-art satellite imaging with on-the-ground confirmation to efficiently manage water resources in the face of climate change.

The satellite images of Landsat-8, Sentinel-2, and MODIS were downloaded for the HRB for the month of September 2022, with cloud coverage of less than 15% from https://earthexplorer.usgs.gov/, https://search.earthdata.nasa.gov, and https://nsidc.org/data/MODIS. MODIS images were originally downloaded in hdf format and converted into Geotiff with a projection of UTM using HegTool Ver 2.15 Build 9.8. Table 2 provides information on salient features of used satellite images including the Universal Transverse Mercator (UTM) zone, date, resolution, and cloud cover. The two tiles for the MODIS image located in 42N and 43N UTM zones were used to cover the study area where zone refers to a particular geographic region on the earth's surface. Furthermore, these images were downloaded after the field visit of three glaciers so that a more precise glacier extent could be extracted.

The precise mapping of SDC glaciers is impeded by the presence of seasonal snow cover. As a result, the glacier mapping was conducted only in areas devoid of seasonal snow cover. Based on the available data and visual evidence, it is apparent that the process of snow or glacier melting persists in the Karakoram Mountain range until the end of August, followed by a subsequent decline in temperature throughout the month of September and snow accumulation again starting in October (Rittger et al. 2013; Atif et al. 2015). The identification of the lowest snow-free area, which refers to the region that remains covered by glaciers throughout the year, in the Karakoram mountains occurs around September. Subsequently, an analysis was conducted on the alterations seen in each of the three glaciers. The area was surveyed via a supervised classification technique that employed visible near-infrared (VNIR) and SWIR bands to delineate and analyze the various characteristics of the glacier features. A threshold range of <0.400 was established in order to delineate regions without snow (Atif et al. 2015).

Due to the remoteness, vastness, and inaccessible nature of the mountain glaciers, satellite data are likely an effective tool for regular mapping of glaciers in a comprehensive and effective manner (Shukla et al. 2009; Bhambri et al. 2011). Over the years, a number of remote sensing techniques for automated mapping of glacier snow and ice by means of multispectral classification have been available. Commonly used techniques such as single-band ratios (Paul et al. 2011) and NDSI (Nagai et al. 2013) take advantage of the high brightness of snow and ice in visible wavelengths to separate them from darker areas such as rock, soil, or vegetation. However, the greatest difficulty in glacier mapping from remotely sensed data is the presence of debris-covered areas on glaciers (Bolch et al. 2008). Glacier areas covered by debris confound the existing techniques because of similar VIS and NIR spectral signatures to the surrounding terrain, lateral, and terminal moraines – outside of a glacier margin (Florath et al. 2022) even with additional information of topography, thermal data, surface texture, and geomorphometric features (Paul et al. 2004; Shukla et al. 2010; Bhambri et al. 2011; Pellicciotti et al. 2015) for neighborhood analysis (Florath et al. 2022).

However, these automatic methods have constraints to extracting the glaciers from SD areas precisely, depending on various factors including steepness of the terrain slope, spectral reflectance, cloudiness, type of snow or ice, and required field measurement. For example, Gjermundsen et al. (2011) highlighted that all the lakes in their study area were classified as glaciers through the application of NDSI due to the similar reflectance of turbid lakes and glacier ice. Also, automated mapping is not generalized and generates reproducible results (similar outlines are going to be generated with the same threshold values) (Paul et al. 2013). Moreover, the SDC area remains underestimated by 8–23% using the band ratio and NDSI methods, respectively (Pratibha & Kulkarni 2018).

Keeping in view the above limitation, the mapping for SDC was completed by the field survey for each selected glacier. The details of the field survey are placed in a subsequent section. The field survey consisted of marking the glacier snout points for each studied glacier. A total of 240 snout points covering latitude, longitude, and altitude were noted for demarcating the SDC areas. The details of the surveyed snout points for each studied glacier are provided in Supplementary Tables S1–S3.

A preliminary site assessment evaluated field conditions and topography, finalizing the study's implementation plan. Data were processed using advanced Global Navigation Satellite System (GNSS) software, but comprehensive field mapping of the glacier is not feasible due to challenging terrain and climatic conditions.

The field survey was conducted from 19 August to 25 August 2022, to measure the SDCs of three glacier sites as depicted in Figure 1. The mapping of the glacier terminus/snout was conducted in the field using a combination of long-range prismless total station and handheld GPS devices and collected a sequence of data points along the outside edge of the terminus.

In this study, the selection of the position for the Real-Time Kinematic Differential Global Positioning System (RTK DGPS) was determined by stability, safety, and the sky view. A control network with benchmarks and control points was implemented for optimal precision and minimal spatial discrepancies. The benchmark in close proximity to the project area was initially built in the WGS1984 coordinate system, which serves as the default system. GPS control points were established, observations were taken, and it was assured that at least four satellites with a Geometric Dilution of Precision (GDOP) value of 5 were available for all-time observations. The control points already established by GNSS were used to carry out the survey along the boundary of SDC. One of the established control points was utilized as a reference for setting the total station, while the other established control point was utilized for obtaining readings during the back station process. Points known as change points (CPs) were strategically placed at regular intervals to mitigate the impact of the curvature factor on the true ground heights. Spot heights were measured generally at close intervals, particularly at abrupt changes in the ground.

The GNSS instrument, Satlab SL800 provided by SATLAB Geosolutions (for converting the coordinate and projection system), and computer software SetSurv (to transfer control networks to layout through GIS/AutoCAD) were used. The accuracy of Satlab SL800 in RTK mode is established as horizontal/vertical ±10 mm/ ± 20 mm +1 ppm. The recorded field survey data through Total Stations and GNSS was transferred and digitized in a CAD environment to check the quality and gaps. Further, the observed data was also exported in CSV format and KML file for further processing in the GIS environment and presented in Supplementary Tables S1–S3.

Differences in sensor resolution, spectral band sensitivity, picture acquisition timing, and image processing methods affect the variability in glacier area estimations from the Landsat-8, Sentinel-2, and MODIS satellites. Since Landsat-8 has a better resolution (30 m) than MODIS (500 m), it can catch finer features, which may result in more accurate boundary delineations (Concha & Schott 2014; Azzoni et al. 2018). Sentinel-2 is in the middle of the two with its 10-m resolution (Azzoni et al. 2018). The sensitivity of each satellite to snow and ice is influenced by variations in spectral bands, especially when there are changes in debris cover and meltwater. Additionally, because glacier boundaries may be obscured by seasonal snow cover, the timing of image acquisition might have a considerable impact on the results. Discrepancies can also be caused by the processing methods, which may differ in how they address elements such as cloud cover and image calibration.

Debris-covered glaciers do, in fact, have distinct spectral properties that require specialized methods for precise identification and extraction. Azzoni et al. (2018) concluded that medium-resolution Landsat ETM+ are ineffective for evaluating the SDC in the Alpine context and Sentinel-2 proved only for preliminary mapping of debris features. Therefore, the studies emphasized the field survey to estimate the precise SDC area (Zhang et al. 2019; Xie et al. 2020).

Based on the results, it is apparent that none of the three satellites were able to detect the SDC. Consequently, a substantial variation was detected in the assessment of glacier extent. This discrepancy has the potential to introduce large errors in the prediction of present and future calibrations of flows. The HRB region is characterized by the presence of pristine glaciers at higher altitudes, exhibiting high reflectivity in their upper sections. In the intermediate sections, these glaciers are partially covered by thin layers of debris, while the lower parts, particularly in the snout area, are heavily covered by debris.

The study also reveals that a substantial proportion of glaciers are covered by debris and on average 20.326% of the SDC area remained undetected by three satellite images. Moreover, all of the glaciers exhibit partial debris coverage, particularly the Pissan glacier having the highest undetected SDC area of 33% by Sentinel-2 (Figure 7(b)). The findings of the current study align with prior research on glacier extent and snow and debris cover (SDC) in the Hunza basin where field measurements conducted in the basin have consistently documented a decline in glacier coverage (Shrestha & Nepal 2019). It is also evident that there is a persistent decline in glacier coverage over time, accompanied by an increase in SDC (Alifu et al. 2016). The utilization of the current SDC findings in modeling can facilitate the estimation of forthcoming changes in flow patterns and mass balance depletion caused by climate change.

This study utilized Landsat-8, Sentinel-2, and MODIS satellite products to estimate the NDSI of the Gulkin, Gulmit, and Pissan glaciers. Subsequently, ground verifications were conducted to assess the supraglacial portion of these three glaciers. The following conclusions are inferred from the results presented in this study:

• The less variability in quantitative glacier extent measurements by Landsat-8 compared with others shows a good agreement with other glacier extent datasets such as GLIMS, RGI, and aerial surveys.

• The glacier extent remained underestimated by the NDSI of Landsat-8, Sentinel-2, and MODIS because of the undetected significant SDC that ranges from 7.652 to 33.548% where MODIS produces the highest undetected SDC area.

• The presence of an additional supraglacial area of 4, 1.530, and 3.200 km² in the Gulkin, Gulmit, and Pissan watersheds, respectively, emphasizes the need to validate the accuracy of glacier extent prediction by the imagery data for water balance and hydrological model studies.

• The SDCs areas of 18.433, 23.980, and 11.642% by Landsat-8, Sentinel-2, and MODIS for Gulkin, Gulmit, and Pissan glaciers may present an average of 18.018% glacier area under debris cover in HRB.

In order to predict the impact of climate change on future flows as well as the loss of mass balance, the results of the present SDC can be used in the models. This study shows that Landsat-8 and Sentinel-2 data with GLIMS and RGI inventory with few changes can be used for SDC studies at basin scale. The Landsat-8 data exhibits better results compared with Sentinel-2 and MODIS platforms.

The authors are thankful to USGS, NASA, and NSIDC for providing the Landsat-8, Sentinel-2, and MODIS satellite images for the SRTM digital elevation model. The author also wants to acknowledge GLIMS and RGI for providing glacier extent data from the websites.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.