Heterogeneous mass balance of selected Glaciers in the Hindu Kush, Karakoram, and Himalaya between 2000 and 2018

ABSTRACT

Glaciers in the Hindu Kush-Karakoram-Himalayas (HKH) are a major source of freshwater in the agriculture-dependent economy of Pakistan. In recent decades, mountain glaciers have been threatened by global warming. In this study, we estimated the Equilibrium Line Altitude (ELA) and geodetic mass balance of fifteen representative glaciers from Hunza (Karakoram), Chitral (Hindu Kush), and Astore (Himalaya) using Landsat satellite images and Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER) DEMs between 2000 and 2018. The climatic trends (temperature and precipitation) in the 19 years (1995–2013) time series of the three regions were assessed using the non-parametric Mann-Kendall test and Sen’s Slope. The ELA of the observed glaciers except for two Atrak (Chitral) and Gulkin (Hunza) is shifted upward, while the mass balance indicates a heterogeneous pattern ranging from −0.23 ± 0.05 to −0.01 ± 0.015 m w.e a⁻¹, −0.13 ± 0.05 to+0.17 ± 0.11 m w.e a⁻¹ and −0.03 ± 0.02 to+0.23 ± 0.09 m w.e a⁻¹ in Chitral, Hunza, and Astore basin, respectively. A significant or slight increasing or decreasing trend was found in the mean annual temperature, mean total precipitation, and mean discharge of the studied basins. This study will be a good contribution to understanding snow and glacier dynamics in HKH, the factors that influence them, and their interaction with the environment.

KEYWORDS:

• Equilibrium Line Altitude
• Hindukush-Karakoram-Himalayas
• Geodetic method
• Glacier dynamics
• Mass balance

Introduction

The glaciers in the HKH region of the Upper Indus Basin (UIB), Pakistan are a major source of freshwater for agriculture, hydropower, industries, and domestic use (Anwar & Iqbal, 2018; C. Gul et al., 2017; Hayat et al., 2019). These water reserves are threatened by climate change with variable rates in different regions (C. Gul et al., 2017; You et al., 2017). It is challenging to assess glacier changes in such inaccessible and vast regions, yet it is pertinent to study their contribution to the freshwater supply and the impacts of changing climate on them (Alifu et al., 2015; Armstrong et al., 2019; Gardelle et al., 2013; Kääb et al., 2012; Muhammad et al., 2019). The glacier retreat or advance not only has impacts on local and regional scales but also works on a global scale significantly contributing to sea levels. Therefore, the data from the observation of major mountain ranges are pertinent in improving the understanding of glacier changes, which is an important tool to gauge climatological and hydrological processes and hazards. Remote sensing data sets are used for monitoring glacier change indicators, such as Snow Line Altitude (SLA), which is a good proxy for ELA, elevation changes, and mass balance (Alifu et al., 2015; C. Gul et al., 2017; Racoviteanu et al., 2008, 2010).

In recent decades, HKH glaciers have lost mass due to global warming (Brun et al., 2017; Sandhu et al., 2018; Wu et al., 2020). Concurrently, the frequency of natural hazards has increased, such as glacier surges and glacial lake outburst floods (Iturrizaga, 2011; Paul et al., 2017; Salerno et al., 2017). However, HKH glaciers are reported to have a heterogeneous response to climate (Berthier & Brun, 2019; Brun et al., 2017; Wu et al., 2020). There have been few recent studies on individual glaciers in HKH, and most of them are based on the regional scale (Paul et al., 2017; Sandhu et al., 2018). This study aims to extend glacier observations to the recent periods 2000 and 2018 in the HKH region and assess their relation to regional climatic factors. We have selected five representative glaciers from each comparable basin in HKH for glacier changes (in terms of ELA, surface elevation changes and mass balance). Landsat data sets and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Digital Elevation Models (DEMs) were used to estimate the glacier changes with the satellite images from the years 2000 and 2018.

Materials and methods

The study area and a flow chart of the methodology used for estimating ELA and geodetic mass balances are given in Figure 1.

Figure 1. (A) Location map showing the studied glaciers (as red polygons) and basin boundaries (black outlines). (b) Methodology flow chart used for estimating ELA and geodetic mass balance.

Study area

Fifteen glaciers, each having an area >5 km², were selected from Chitral (Hindu Kush range), Hunza (Karakoram range), and Astore (Himalayan Range) basins (five from each basin) to assess the glacier dynamics in the study area (Figure 1(a)). These glaciers were selected based on the availability of minimum or cloud-free data. The physical attributes of the studied glaciers are given in Table 1.

Table 1. Attributes of the studied glaciers.

Basin	Glacier	Area (km^2)	Longitude (Decimal degrees)	Latitude (Decimal degrees)	Elevation range (m a.s.l.)	Mean elevation (m a.s.l.)
Chitral	Atrak	91	71.9	36.4	3716–7094	4984
	Kotgaz	49	72.2	36.7	3629–6970	5009
	Maski-o-Shayaz	44	72.1	36.6	3879–6969	5622
	Noroghikun	31	72.4	36.7	3661–6318	5129
	Lower Tirich	22	71.9	36.3	3854–7420	4921
Hunza	Hispar	499	75.28	36.05	3110–7794	4854
	Passu	62	74.71	36.46	2578–7487	5095
	Yazgil	133	75.3	36.24	3234–7709	5910
	Mulungutti	98	75.2	36.33	2911–7799	5664
	Gulkin	30	74.78	36.42	2556–7295	5364
Astore	Bazhin	13	74.63	36.22	3344–6633	4048
	Ghughuel	11	74.53	35.08	3997–6109	4929
	Dodhar	5.6	74.59	35.1	3808–5237	4360
	Sachen	11	74.74	35.32	3404–5530	3894
	Hunza	6	74.55	35.11	3890–5892	4336

The Chitral basin (34.4–36.5 N, 70.5–75.8 E) is located in the Hindu Kush range and covers an area of 12,400 km². Approximately 1903 km² (16%) of the Chitral basin is covered by glaciers (Naeem et al., 2015).

Table 1. Attributes of the studied glaciers.

The Hunza basin (35.7–37.1 N, 72.5–75.7 E) is in the Karakoram range and is one of the most glacierized regions in the Upper Indus Basin. It covers an area of ~13715 km², 4053 km² (~30%) of which is covered with glaciers (Muhammad and Tian, 2016). In recent years, many glaciers in the basin have surged, such as Hispar, Yazgil, Khurdopin, Mulungutti, and Shishper Glacier (Wu et al., 2020).

Astore basin (35.1–37.0 N, 74.0–75.1 E) located in the western Himalayas covers an area of 3990 km² with 252 km² (6%) of its area covered with glaciers (Muhammad et al., 2019). Approximately 16% of its glacierized area is covered with debris, and the remaining area has many small cirque glaciers (Muhammad et al., 2019b).

Data sets

For the estimation of the ELA, Landsat 7 Enhanced Thematic Mapper (ETM+) images with a pixel size of 30 m were used. The scenes were selected for a two-month window (August–September) which corresponds to the end of the melting season in the study area. Due to the considerable cloud cover, Landsat 7 ETM+ images with minimum cloud cover (preferably less than 10%) for the years 2000 and 2018 were acquired from the United States Geological Survey (USGS)’s online portal (https://earthexplorer.usgs.gov/). Details of these Landsat images are given in Supplementary material (Table i).

The ASTER (AST14DMO) DEMs with a 30 m spatial resolution taken at the end of the ablation period years 2000 and 2018 were downloaded from NASA’s online portal (https://search.earthdata.nasa.gov/search/) to estimate mass balance. The Randolph Glacier Inventory version 6 (RGI6.0) was used for the extraction of glacier outlines (RGI, 2017).

The daily meteorological data (temperature and precipitation) of the study area over 19 years (1995–2013) were obtained from the Pakistan Meteorological Department (PMD) and Water and Power Development Authority (WAPDA).

Methodology

A flow chart of the methodology used for estimating ELA and geodetic mass balances is given in Figure 1(b).

Preprocessing of Landsat images

Landsat satellite data was hampered on 31^st May 2003, when Scan-line corrector onboard the Landsat 7 ETM+ suffered a failure, which resulted in data gaps (around 22% of pixels missing) in each Landsat image (J. Gul et al., 2020). However, USGS provides the gap mask files (compressed in GeoTIFF format) for each band in the downloaded data set. QGIS 3.5 was used to fill the gaps in downloaded images for the year 2018. After gap filling and correcting the atmospheric effect of the Landsat images, the Landsat image of the year 2000 was considered as a reference image for the co-registration of the corrected images. To detect and adjust geomatic errors, the acquired Landsat images were orthorectified using a reference image.

Estimation of mass balance

The DEM of Difference technique was used for the estimation of the surface elevation change in the studied glaciers. After computing the DEM differences, it was divided by the time separation (18 years) between the DEMs to attain maps of elevation change rate (dh/dt) (Dussaillant et al., 2018; Neckel et al., 2013; Vincent et al., 2013). The elevation change of each glacier with its altitude was calculated using 100 m elevation bins. Elevation changes were converted to volume changes by multiplying the mean elevation change rate by the area of each elevation bin. Assuming that only ice is lost or gained, an average ice density of 850 kg/m³ was used for the conversion of volume to mass for all glaciers (Brun et al., 2017; Rankl & Braun, 2016; Wu et al., 2020).

Generation of DEMs from ASTER data to estimate the change in surface elevation of glaciers

The ASTER DEMs were corrected for horizontal and along/cross-track shifts (Gardelle et al., 2013; Nuth & Kääb, 2011). The offsets were determined on stable terrain by minimizing the standard deviation of the elevation change (dh). Stable terrain was considered to be all pixels with valid elevation values excluding water, snow, ice, and clouds. In each 100 m elevation bin, pixels with values greater than three standard deviations from the mean were not considered for elevation change estimation (Berthier & Brun, 2019; Gardelle et al., 2013; Muhammad et al., 2019; Wu et al., 2020).

Accuracy assessment of the Landsat images

The accuracy assessment of the classified Landsat images of 2000 and 2018 of the basins was evaluated using an error matrix (Table 3). False band composite (unclassified Landsat images of 2000 and 2018) and high-resolution Google Earth images of the same year were used as reference images.

Table 3. Evaluation of the supervised classification of Landsat images in comparison to reference images based on an error matrix.

2.3.5. Extraction of glacier outlines

Initial outlines obtained from the Global Land Ice Measurements from Space (GLIMS) database were manually improved using the false-color composite (RGB-543) Landsat images.

Estimation of ELA from Landsat images

The SLA, estimated at the end of the melting season, is concurrent with the ELA. For the delineation of the accumulation area, supervised classification was used on the Landsat data for the years 2000 and 2018. The classified images were reclassified, and accumulation and ablation areas were extracted for each glacier (Benn et al., 2005; Huss, 2013; Owen & Benn, 2005). Once the accumulation area was determined, the Accumulation Area Ratio (AAR) was calculated using the total area of the glacier. The ELA was then determined using the hypsometric curve (area–elevation relationship, determined from the DEM in the years 2000 and 2018) to find the corresponding elevation against the AAR fraction and digitized using contours in ArcGIS 10.8.

Results

Estimation of ELA

The AAR and ELA values for the studied glaciers between 2000 and 2018 are shown in Table 2. The digitized estimated ELAs of the studied glaciers in Chitral, Hunza, and Astore basins for 2000 and 2018 are displayed in Supplementary material Figures I, ii, and iii, respectively.

Table 2. Estimated AAR and ELA of the studied glaciers for years 2000 and 2018 in all three basins.

Basin	Glacier name	AAR		ELA (m a.s.l.)
		2000	2018	2000	2018
Chitral	*Atrak	0.60	0.57	4780	4980
	*Kotgaz	0.53	0.54	5000	4960
	*Maski-o -Shayaz	0.75	0.64	4580	4920
	*Noroghikun	0.62	0.60	4920	5000
	*Lower Tirich	0.48	0.45	4980	5020
Hunza	*Hispar	0.70	0.52	4500	4860
	Passu	0.66	0.64	4580	4660
	Yazgil	0.82	0.79	4500	4700
	Mulungutti	0.79	0.77	4820	4910
	*Gulkin	0.68	0.72	4350	3980
Astore	*Bazhin	0.58	0.53	3890	3960
	*Ghughuel	0.65	0.60	4380	4500
	*Dodhar	0.36	0.34	4480	5095
	*Sachen	0.32	0.26	4200	4380
	*Hunza	0.40	0.38	4420	5270

Note: *Debris covered glaciers.

In the Chitral basin, an upward shift of ELA was observed for Atrak from 4780 m a.s.l. to 4980 m a.s.l. Maski-o-Shayaz experienced the highest upward shift from 4580 m a.s.l. to 4920 m a.s.l. among the studied glaciers of the Chitral basin, Noroghikun and Lower Tirich also had an upward shift in ELA from 4920 and 4980 m to 5000 m and 5020 m, respectively. However, only Kotgaz glacier among the studied glaciers experienced a lowering of ELA from 5000 m to 4960 m.

In the Hunza basin, an upward shift of ELA was observed for Hispar, Passu, Yazgil, and Mulungutti glaciers. Among these glaciers, the Hispar glacier experienced the highest upward shift from 4500 m to 4860 m.

In the Astore basin, the studied glaciers show an upward shift in their ELA from 2000 to 2018. Comparatively, among the studied glaciers in the Astore basin, the Hunza glacier experienced an upward shift from 4420 m to 5270 m.

Estimation of mass balance

The elevation changes and mass balance of glaciers in the Chitral, Hunza, and Astore basins are shown in Figures 2-4. In the Chitral basin, Maski-o-Shayaz experienced the most significant mass loss of −0.23 ± 0.05 m w.e a⁻¹ consistent with a significant upward ELA shift during the same observation period. The Kotgaz glacier has a large accumulation zone in comparison with other glaciers studied. It has undergone a considerable mass loss. Noroghikun has experienced mass loss in the ablation zone. Lower Tirich has ice mixed debris in the ablation zone, while no major contributing tributaries are in the accumulation zone. It has also experienced mass loss.

Figure 2. Elevation and mass balance changes from 2000 to 2018 for the studied glaciers in the Chitral basin.

In the Hunza basin, mass balance varied heterogeneously mostly due to glacier surges (Figure 3). Overall, the surface elevation increased as a function of altitude and varied considerably depending on altitude. Passu has experienced a significantly heterogeneous pattern with surface lowering in the ablation area (below ~4900 m) and thickening in the accumulation areas (above ~4900 m) which is comparable to Muhammad et al. (2019). Yazgil glacier in the Hunza basin experienced the most significant mass loss of −0.13 ± 0.05 m w.e a⁻¹. Yazgil and Hisper glaciers are reported to be surging (Paul et al., 2017; Wu et al., 2020).

Figure 3. Elevation and mass balance changes from 2000 to 2018 for the studied glaciers in the Hunza basin.

The mass balance of glaciers in the Astore basin also varied heterogeneously (Figure 4). Similar estimates of heterogenous mass balance in the Astore basin were reported by Brun et al. (2017). The glaciers in Astore are comparatively smaller than those of the Chitral and Hunza basins (Table 1) and are reported to be greatly influenced by topography and local climate changes (Kumar et al. 2019). A significant positive mass balance was estimated for Bazhin and Dodhar glaciers. This can be due to the regional climatic similarities with the neighboring Karakoram Range Muhammad et al. (2019). Bazhin is a debris-covered glacier and has experienced a surge event in the recent past (Muhammad et al., 2019). Dodhar is a heavily debris-covered glacier and has heterogeneously experienced a surface decrease and increase in the ablation and accumulation area with a positive net mass balance. Sachen glacier is stable during the observed period. The Hunza glacier is almost stable during the observed period. Among the observed glaciers in the basin, the Ghughuel glacier has slightly lost mass in the observation period (below 4300 m).

Figure 4. Elevation and mass balance changes from 2000 to 2018 for the studied glaciers in the Astore basin.

Estimation of climatological trends

Trend analysis of climate data was performed using 19 years (1995–2013) time-series data of Chitral, Hunza, and Astore basins as shown in Figure 5. In the Chitral basin, the Mann Kendall (MK) trends and Sen’s slope values of+0.41 (p < 0.05) and+0.18 (°C/yr) suggest a significant increase in mean annual temperature. MK trend test of total precipitation indicated a slight negative trend with a value of −0.192 (p > 0.05), while Sen’s slope value of −0.56 (mm/yr) indicates a slight decreasing trend (Figure 5a).

Figure 5. Trends in mean temperature and total precipitation of the a) Chitral, b) Hunza and c) Astore from 1995–2013.

MK trend and Sen’s slope values of+0.68 (p < 0.05) and+0.2 (°C/yr) of mean temperature suggests a significant positive trend with a value of+0.68 (p < 0.05) in the Hunza basin (Figure 5b). MK trend test and Sen’s slope value of total precipitation suggest a significant increase with a value of+0.32 (p < 0.05) and+0.78 (mm/yr), respectively (Figure 5b).

The MK trend test and Sen’s slope implied a significant increase in mean annual temperature with a value of+0.36 (p < 0.05) and+0.23 (°C/yr) in the Astore basin, respectively (Figure 5c). These tests indicated a significant decrease with values of −0.38 (p < 0.05), and −0.78 (mm/yr) in total annual precipitation (Figure 5c). These values suggest contrasting trends in mean annual temperature and total annual precipitation in the Astore basin.

Discussion

The variability in the mass balances using the geodetic method is likely due to the complex terrain and environment of glaciers in the region. Some glaciers are thickening in the ablation zones, possibly due to thick debris cover. Most surface lowering is observed in the middle elevations of the glaciers, comparable to the ground observations (Muhammad et al., 2019). Overall, there is no significant difference in surface lowering of debris-covered and clean ice observed at the same elevation. A comparison of the estimated mass balances in this study and previously published studies in the same region are given in Supplementary material Table iii.

The glaciers in the Chitral basin are significantly debris-covered. Avalanches, landslides, and rockfalls from surrounding rocks (triggered by a freeze-thaw process) accumulate as debris on the glacier. J. Gul et al. (2020) have estimated the ELA, mass balance, and snout position of six glaciers and reported an increase in the ELA from 2001–to 2018 and a decrease in mass balance and retreat of snout positions. Mean annual temperature and total precipitation in the Chitral basin were found to be increasing over the study period. J. Gul et al. (2020) have reported similar climate trends in the Chitral basin and Ahmad et al. (2018) have reported a significant increase in the mean annual temperature and linked it with negative mass balance and upward shift of ELA. The mass balance of studied glaciers in the Hindu Kush, Karakoram, and the Himalayas in this study are presented in Figure 6.

Figure 6. Mass balance in the a) studied glaciers in Chitral (2000–2018). b) studied glaciers in Hunza (2000–2018). c) studied glaciers in Astore (2000–2018).

Among the studied glaciers in the Hunza basin, Hispar has a major debris-covered portion, which may have contributed to an upward shift in the ELA of the glacier. However, Hispar has repeatedly surged in the past few decades (Paul et al., 2017; Wu et al., 2020). The elevation change of the Hispar glacier is not continuous along the main trunk, with increasing and decreasing elevation change patterns. When the glacier surges, it sometimes deposits the disintegrated mass on the main trunk or tongue of a glacier, which might explain the discontinuous thinning and thickening along the main trunk. Only the ELA of the Gulkin glacier in the Hunza basin shifted downward from 4350 m to 3980 m, indicating an increase in the accumulation area over the period which is consistent with the ground-based elevation changes in the lower part of the glacier estimated by Muhammad et al. (2019). Although it has a debris-covered tongue, it may indicate the effect of other factors, such as the thickness of debris cover or topography. Consequently, small differences between the estimated ELAs in this study and other studies may be due to different observation periods and locations. Racoviteanu et al. (2019) reported an average ELA of 5,177 ± 108 m a.s.l. over 2000–2016 in the Hunza basin. While Minora et al. (2016) reported an ELA of 5200–5300 m a.s.l. in the Central Karakoram from 2001–to 2010. In addition, Kääb et al. (2012) reported the highest ELA of 5540 m a.s.l. in the Karakoram region. Mukhopadhyay and Khan (2016) and Racoviteanu et al. (2019) reported a strong east-to-west precipitation gradient in the Karakoram, which is one of the reasons for ELA and mass balance variability across the region. Our elevation changes estimate indicates a heterogeneous pattern.

We found heterogeneous glaciers’ status with a stable net mass balance in the Astore basin. Thick debris cover and climatic effects of monsoons and westerlies explain the conditions of the balanced glacier in the western Himalayas. Other studies have shown a rapid thinning of glaciers in the western Himalayas. Berthier et al. (2007), Gardelle et al. (2013) and Kääb et al. (2015) have reported a negative mass balance of −0.49 ± 0.12 and −0.53 ± 0.16 m w.e a⁻¹ during 2003–2008 and 1999–2011, respectively. There are local climatic differences within and across the region (Kayastha and Harrison, 2008). Muhammad and Tian (2016) have reported similar trends in the climate data in Astore from 2000–to 2015 as reported in this study.

The variability of glacier mass balances is a net result of climatic forcing, surface mass balance changes, and ice flow. The net effect of a specific mass balance may be small due to the cancelling effects of other drivers (Kumar et al., 2019). Whereas ELA was determined for the years 2000 and 2018, it is highly influenced by snow cover and temperature of the region. For this purpose, historical snow-cover and climate data should be accessed and processed for the basins of HKH to generate snow-cover data to find significant trends over a longer timescale to understand the interaction between climatic variables, ELA and mass balance.

Simultaneously, the inefficient distribution of water supply is the source of contention between the provinces of Pakistan, which can escalate due to the uncertain availability of water. The heterogeneous mass balance and the uncertainty over the snow and glacier melt will further escalate the situation, especially with respect to this region, which is prone to climate hazards, particularly floods. The socioeconomic conditions of the people living in the region and the intensity of the impact of climate change on them should be considered in developing sustainable livelihood plans involving Nature-Based Solutions.

Conclusion

The main conclusions abstracted from the results of this study are:

• An upward shift of ELA for 13 of the 15 studied glaciers was observed, indicating a decrease in the accumulation area. Mass balance estimated from DEM signified heterogeneous behavior of the studied glaciers in the basins. We found heterogeneous mass balances ranging from −0.23 ± 0.05 to −0.01 ± 0.015 m w.e a⁻¹, −0.13 ± 0.05 to+0.17 ± 0.11 m w.e a⁻¹ and −0.03 ± 0.02 to+0.23 ± 0.09 m w.e a⁻¹ in Chitral, Hunza and Astore basin, respectively.

• The mean annual temperature suggested a significant increase in trend in the studied basins. The total annual precipitation in the Chitral basin was slightly decreased, while there were significantly increasing trends in the Hunza and Astore basins. Variations in the response of glaciers in terms of ELA and mass balance are most likely due to differences in local climatic conditions, topographic features, and glacier morphology.

The heterogeneous behaviour of glacier demands more long term, consistent, detailed field and remote-sensing studies in the region. Future work should include the response of glaciers to climate change in the region with precise ground measurements and high-resolution remote sensing data. The unsustainable supply of freshwater in Pakistan requires a proactive water resource management by the relevant Government authorities, policymakers and stakeholders to ensure economic stability in the country.

Acknowledgments

This work was financially supported by the Grant number 42171129 of the National Natural Science Foundation of China (NSFC). ICIMOD gratefully acknowledges the support of its core donors: the Governments of Afghanistan, Australia, Austria, Bangladesh, Bhutan, China, India, Myanmar, Nepal, Norway, Pakistan, Sweden, and Switzerland. The views and interpretations in this publication are those of the authors and they are not necessarily attributable to their organizations.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/22797254.2023.2173086