The Karakoram Anomaly: Validation through Remote Sensing Data, Prospects and Implications

Abstract

Millions of people rely on river water originating from snow- and ice-melt from basins in the Hindukush-Karakoram-Himalayas (HKH). One such basin is the Upper Indus Basin (UIB), where the snow- and ice-melt contribution can be more than 80%. Being the origin of some of the world’s largest alpine glaciers, this basin could be highly susceptible to global warming and climate change. Field observations and geodetic measurements suggest that in the Karakoram Mountains, glaciers are either stable or have expanded since 1990, in sharp contrast to glacier retreats that are prevalently observed in the Himalayas and adjoining high-altitude terrains of Central Asia. Decreased summer temperature and discharge in the rivers originating from this region are cited as supporting evidence for this somewhat anomalous phenomenon. This study used remote sensing data during the summer months (July–September) for the period 2000 to 2017. Equilibrium line altitudes (ELAs) for July, August and September have been estimated. ELA trends for July and September were found statistically insignificant. The August ELA declined by 128 m during 2000–2017 at a rate of 7.1 m/year, testifying to the Karakoram Anomaly concomitant with stable to mass gaining glaciers in the Hunza Basin (western Karakoram). Stable glaciers may store fresh water for longer and provide sustainable river water flows in the near to far future. However, these glaciers are also causing low flows of the river during summer months. The Tarbela reservoir reached three times its lowest storage level during June 2019, and it was argued this was due to the low melt of glaciers in the Karakoram region. Therefore, using remote sensing data to monitor the glaciers’ health concomitant with sustainable water resources development and management in the HKH region is urgently needed.

1. Introduction

The study of the changes in the complex alpine glacier system (e.g., glaciers in the Hindukush, Karakoram, and Himalaya) is crucial for water resources development and management. Monitoring alpine glaciers imparts vast information about the response of glaciers to climate change together with the future prospects and implications for prosperity. According to the available studies, the rise in the global average surface temperature has caused glaciers to retreat and shift towards lower altitudes. Thus, glaciers are considered as a key and warning symbol of the escalating temperature due to changes in the climate [1,2].

The Indus River Basin (IRB) is one of the major basins of a giant and complex mountainous range of the HKH where more than 80% of the stream flow is generated by snow-glacier melts [3]. Various studies affirmed IRB as the most vulnerable and threatened region with a significant glacial loss in the Hindukush and Himalayas see e.g., [3]. Based on different glacier inventories, various existing studies estimated the glacier-covered area varying from 36,845 km² to 50,750 km² in the HKH region [4], The Karakoram range, roughly, contains about half of this glacierized area [4], and plays a pivotal role in water supply and energy generation.

Unfortunately, over the past few decades, retreat and negative glacier mass balance have been observed and reported in the central and eastern parts of the HKH region [5,6,7,8]. On the contrary, a few studies concluded that the Karakoram glaciers show diverse behavior and exhibit the thickening and growing of glaciers with a stable/positive mass balance [3,9,10]. This state of stable/positive mass balance in the highly glacierized central parts of the Karakoram is referred to as the “Karakoram Anomaly” [11,12].

Stable glaciers in the Karakoram may bring sustainable water availability to the Indus River; see, e.g., [13]. However, glacier surges and lakes will become common, and will pose a threat to infrastructure in the region; see, e.g., [14,15]. Therefore, there is a dire need for continuous monitoring of the Karakoram glaciers. The rugged terrain and harsh environment of the region restricts in-situ observations and monitoring.

One of the monitoring tools is the estimation of the equilibrium line altitude (ELA). ELA refers to an average altitude on alpine glaciers that represents a divide between glaciers lying in accumulation and ablation zones. This altitude is significantly influenced by the local/regional temperature and precipitation patterns. Thus, this altitude is considered as a fundamental and key measure of a “Glacier-Mass-Balance” in a particular glacier system. The accumulation zone, usually sited on the upper part of the glaciers, signifies “Winter-Glacier-Mass-Balance”, where glaciers experience thickening and growing owing to different processes including precipitation fall, avalanches and wind-blown snow. On the contrary, the lower part of the glaciers indicated by “Summer-Glacier-Mass-Balance” refers to an ablation zone, where thinning and shrinkage of glaciers can be observed. The quantitative measure of location of this altitude and its long-term fluctuations over the ablation period is of significant importance to consider the hydrological processes within glacier systems. If the ELA rises towards high altitudes, the glaciers indicate negative mass balance, and vice versa. However, none of the available studies provides a time-series ELA variation for the HKH region.

This study, therefore, utilizes remotely sensed images to estimate the ELA, and uses it for glacier health monitoring from 2000 to 2017 in the Hunza Basin. The temporal trend analysis of the ELA in such basins can significantly contribute to understanding glacier response to climate change. The variability in the historical ELA can also be used to assess the possible future implications and sources of sustainable water resources.

The Study Area

The Hunza River basin, a highly glacierized sub-basin of the HKH, runs through high mountains of the Central Karakoram. The basin lies within an altitudinal range of 1415 m to 7809 m with a mean altitude of 4570.61 m. The total surface area of the study area is 13,733 km^2, stretching between 35°54′18.29″ N to 37°5′47.86″ N (latitude) and 74°0′37.76″ E to 75°47′34.72″ E (longitude) (See Figure 1). According to the Randolph Glacier Inventory (RGI v 6.0), the total glacier area covered by the study area is approximately 4290.44 km², making 31.27% of the overall study area. The RGI shows that the total number of glaciers in the Hunza River basin is about 2329. The statistics show that 78.08% of the glacier area in the Hunza River basin is occupied by large glaciers s greater than 5 km² in size; 16.83% of the glaciers are in sizes ranging from 1 km² to 5 km², 8.73% of the glaciers have sizes varying from 0.05 km² to 1 km² and 0.36% of the glaciers are less than 0.05 km². About 85% of the Hunza River basin remains covered in seasonal snow in the winter season, which reduces to 30% in the summer season. For measuring precipitation, there are three-gauge stations installed by WAPDA at different elevations, i.e., 4730 m (Khunjerab), 3669 m (Ziarat) and 2858 m (Naltar). These are shown in Figure 1b. According to the record for 1999–2008, the average precipitation at these stations is 165 mm/year, 292 mm/year and 660 mm/year, respectively [16]. One rain gauge station at an elevation of 2156 m (Hunza) was installed by PMD in 2007. The available stations data do not show a direct increase with altitude due to the complex nature of the HKH region. The precipitation variation depends on various factors in the Karakoram, such as altitude, latitude, longitude, orographic effects, and moisture from the westerlies and monsoon. The Hunza River carries an average annual flow of 323 m³/s recorded over the period from 1966 to 2008 at Dainyor Bridge [17]. The summer flows, which originate from snow and glacier melt, are about 70% to 80% of the annual flows.

2. Materials and Methods

2.1. Materials

In this study, we used a variety of satellite datasets including a Moderate Resolution Imaging Spectro-radiometer (MODIS) 8-day snow cover product (500 m spatial resolution), Landsat imageries (30 m spatial resolution), an RGI v.6.0 (Randolph Glacier Inventory), and MODIS 8-day Surface Reflectance (250 m spatial resolution).

For appropriate hydrological-climatic studies there is a dire need for accurate evaluation of snow- and ice-covered extents. For this purpose, the United States’ Space Agency National Aeronautics and Space Administration (NASA) employed two satellites, Terra and Aqua, launched in 1999 and 2002, respectively, which are operated by a MODIS instrument. Both Terra and Aqua satellites managed to take an image of the same part of the Earth at 10:30 am and 1: 30 pm, respectively, and hence provided two values of snow-covered area for the same site per day. In this study, snow-cover product from MODIS/Terra (MOD10A2) was utilized as the principal data source on the basis of the study conducted by Wang et al. [18] for China, in which the MODIS/Terra showed better results as compared to MODIS/Aqua for the least cloud coverage. Besides this, many authors have successfully used the MODIS/Terra snow cover product in the Upper Indus Basin (including the Hunza Basin) or parts of the Upper Indus Basin such as [1,17,19,20,21,22]. Therefore, we downloaded the MODIS 8-day snow cover product for 18 years from 2000 to 2017 followed by the mosaicking of two tiles (h23v05 and h24v05 shown in Figure 2) to swathe the entire Hunza River basin.

In spite of having a coarse spatial resolution of 500 m, the MODIS snow-cover product compared to Landsat images is a widely used dataset on account of the availability of long-term series daily and 8-day data.

Similarly, MODIS/Terra Surface Reflectance 8-day L3 Global 250 m SIN Grid V005 (MOD09Q1 with spatial resolution of 250 m) for 2000 to 2017 were downloaded from (https://reverb.echo.nasa.gov/reverb/#utf8=%E2%9C%93&spatial_map=satellite&spatial_type=rectangle) followed by mosaicking the two tiles (h23v05 and h24v05) to cover the study area (accessed during 1–30 January 2019).

In addition, the Randolph Glacier Inventory (RGI) version 6.0 data was downloaded using the website: http://www.glims.org/download/ (accessed during 1–30 January 2019).

2.2. Methods

A step-by-step methodology of the current study is provided in the following schematic diagram (see Figure 3). The overall analysis can be divided into two parts: (i). Pre-processing the MODIS snow-cover data, and (ii). Post-processing the MODIS snow-cover data. Each step is explained in the following paragraphs:

• After downloading 430 snow-cover images, and mosaicking the same for the study area, No-Data Cells were searched, using a methodology explained by Mukhopadhyay [23]. However, none of the images contain No-Data Cells for the study area.
• Misclassification of snow-ice due to the presence of cloud cover in MODIS images is the prime limitation in the use of MODIS snow cover data for estimating snow- and ice-covered areas. Since the HKH region is persistently covered in clouds, it is imperative to rectify the MODIS images from cloud coverage to curtail the effect of cloud cover and to avoid under-estimation and/or over-estimation errors. The technique used in the current study to reduce cloud coverage from the MODIS data results in removing 100% of cloud cover over the glacier area. Therefore, the cloud cover in MODIS images was identified by masking the cloud pixels over the glacier region by using the global glacier inventory RGI v 6.0. One such example is shown in Figure 4a in which the original MODIS image (dated 22 September 2003) is shown with initial cloud coverage and in Figure 4b a processed image after reduction in the cloud cover is shown. A total of 191 out of 215 MODIS images show cloud-cover pixels less than 4% of the total pixels while 24 images possess cloud cover ranging from 4% to 14% over the snow-covered area. It is noteworthy that underestimation of the glacier-cover area may adversely affect the ELA estimates. Therefore, the exclusion of cloud cover over the glacier area improves the accuracy of the ELA estimation. However, since the main focal area for determining the ELA is the glacier area, the cloud cover in the surrounding areas or over the snow area does not influence the results.

Figure 4. MODIS image of 22 September 2003. (a) shows the original cloud coverage in that particular image. (b) shows the processed MODIS of same date after reduction in cloud cover. In (c) it is clearly visible that clouds are removed over the glacier area.

Figure 4. MODIS image of 22 September 2003. (a) shows the original cloud coverage in that particular image. (b) shows the processed MODIS of same date after reduction in cloud cover. In (c) it is clearly visible that clouds are removed over the glacier area.

• One of the key issues associated with the MODIS product is that it cannot identify small glaciers due to its coarse resolution and that glaciers with a size less than 0.01 km² may be under-estimated by MODIS data [20]. Additionally, MODIS data cannot differentiate between snow and ice extents. For this purpose, the RGI data was used to reinstate the small glaciers in the MODIS data. The re-instatement of small glaciers in the Hunza River basin shows a variation in the small glaciers over the ablation period (July–September) ranging from 0.5% to a maximum 7.2% of the total glaciers in the study area. The smallest percentage of reinstatement could be attributed to the wet year with a maximum snowfall and low melt, whereas the maximum reinstatement could be associated with dry years with minimum snow-fall or excessive snow-melt in a particular year.
• MODIS/Terra Surface Reflectance 8-day L3 Global 250 m SIN Grid V005 (MOD09Q1 with spatial resolution of 250 m), a very powerful and useful albedo data particularly for separation of snow and ice, shows variation in albedo for snow and ice surfaces which influences the melt rate significantly. The lower and upper limits of surface reflectance provided in the albedo data for the study area range from −100 to 16,000′; however, no such threshold is available to separate the snow and ice from this data. Khan et al. [20] separated snow and ice for some of the sub-basins in the UIB using a single time series dataset Landsat images. The Albedo images together with the pre-processed MODIS images of the same date as were used for the Landsat images by Khan et al. [20] were selected for determining the threshold. Therefore, by using a trial-and-error technique, the threshold value was continuously adjusted until the same/approximated snow/ice-covered area was obtained as was estimated by Khan et al. [20]. The threshold limit of albedo (surface reflectance) for separating snow and ice over the glacier area was found to be 3241 with a value of less than 3241 indicating an ice-covered area and the value above 3241 providing a snow-covered area.
• The summer months (July, August and September) show maximum exposed ice (see Figure 5), and hence were selected for estimation of the ELA. The separated snow- and ice-covered images were used to estimate the ELA, using SRTM 90 m DEM. To acquire monthly ELAs, weekly ELAs were first estimated, and then averaged for the particular month.

Figure 5. Processed MODIS images for each month from the year 2011 are shown. It is clearly visible that a major part of the ice is exposed during the ablation period (July, August and September) due to the melting of seasonal snow.

Figure 5. Processed MODIS images for each month from the year 2011 are shown. It is clearly visible that a major part of the ice is exposed during the ablation period (July, August and September) due to the melting of seasonal snow.

• The temporal trend analysis of the ELA during the ablation period, i.e., the months of July, August and September, was then carried out using Mann–Kendall and Sen’s slope non-parametric tests. The brief description and reasons for selection of these tests are provided below:

2.2.1. Non-Parametric Tests

Parametric statistical tests usually require some strict assumptions to be satisfied for their correct application (e.g., normality in data distribution, sample/data independency, linearity in data, etc) [24]. However, climatic data due to continuing changes in climatic conditions may infringe these assumptions. Thus, distribution-free statistical tests also known as non-parametric tests are often preferred over parametric tests for time-series trend analysis of hydro-climatic variables. They are less sensitive to the outliers/missing values present in the data and are also suitable for small-size sample data. The Mann–Kendall test (M–Kt) is the simplest and most widely used non-parametric test for identifying the existence/significance of trends in hydro-climatic variables. For the detection of climate changes and their relationship with hydro-climatic variables, it is important not merely to determine the significance of the trend but also to estimate its magnitude. Thus, the M–Kt complemented with the Sen’s slope test (S-St) can be used to estimate trend magnitude.

The existing studies show that several authors, around the world, have employed the M–Kt in trend analysis of various hydro-meteorological variables and sediment loads, e.g., [25] used M–Kt for evaluating temperature and precipitation patterns in Eastern Slovakia. Latif et al. [26] used M–Kt together with S-St to detect the trends and changes in seasonal and annual precipitation in the UIB. Many researchers have applied the M–Kt along with Sen’s slope pattern in the UIB. Trends in the monthly, seasonal and annual precipitation for the Amhara regional state of Ethiopia were analyzed by Gedefaw et al. [27]. Gumus et al. [28] carried out trend analysis of various hydro-meteorological variables including temperature, total precipitation, humidity and wind speed in Sanliurfa, Turkey. Salami et al. [29] examined trends in various hydro-meteorological variables including temperature, rainfall, sea level rise, relative humidity and wind speed in the coastal region of Lagos, a city in West Africa.

2.2.2. Mann–Kendall Test, M-Kt

M–Kt is the commonly used and robust statistical method presented by Mann [30] and Kendall [31] to investigate the existence of significant trends in time-series datasets of hydro-climatic variables. This method is not restricted in/by the use of specific types of data distribution as in the case of parametric test, and thus can be used for non-normal distributed datasets. The efficiency of the test is validated by comparing it with other statistical tests [32]. The significance level (SL) indicated by “p-value” is used to accept or reject the “Null-hypothesis (H_o)” of the M–Kt. If p-value > SL, the trend is significant and p-value < SL, the trend is insignificant. In this study, trend analysis is carried out using two significance levels, i.e., 5% and 10%.

The M–Kt statistic “S_MK” is determined as:

$$S_{MK}=\sum _{P=1}^{N-1}\sum _{Q=P+1}^N\mathrm{sgn}(X_Q-X_P)$$

(1)

Here, “N” represents the number of observations in the data sample. X_P and X_Q are the consecutive values in the data sample. Taking (X_Q − X_P) = Y, the value of “sgn Y” can be calculated as:

$$\mathrm{sgn}\mathrm{Y}=\left\{\begin{array}{lll}1& if& \mathrm{Y}>1\\ 0& if& \mathrm{Y}=1\\ -1& if& \mathrm{Y}<1\end{array}\right.$$

(2)

The positive value of S_MK shows a rising trend while the negative value indicates a declining trend.

2.2.3. Sen’s Slope Test, SSt

SSt (non-parametric test) is coupled with MKt to determine the slope of a trend line identified in a hydro-meteorological time series datasets. This method is based on the assumption of a linear trend present in the dataset. The Sen’s slope (SS) is calculated by:

$$SS=Median\left. [\frac{X_Q-X_P}{Q-P}\right]$$

(3)

where

• time Q > time P
• X_P = data value taken at time P
• X_Q = data value taken at time Q

The Sen’s slope (SS) provides the information on change in any hydro-meteorological variable with respect to time, which helps in estimating the total change occurring over the analysis period.

3. Results and Discussion

The results of this study suggest that an average ELA for the Hunza River basin is 5022 m ± 103 m during 2000 to 2017. Similar ELA values for the study area have also been estimated using different techniques and datasets. Butt et al. [33] estimated this altitude between 5000 m to 5500 m for the HKH region by using the techniques of supervised and unsupervised classification of Landsat 7 images coupling with ASTER DEM (30 m resolution). Shrestha et al. [16] suggested a mean ELA of 5050 m for the Hunza River basin using hypsographic analysis. Scherler et al. [34] estimated the ELA for the Hunza Basin is 5140 m, using Landsat images. Hewitt [15,35] reported the ELA for the Hunza Basin is 4800 m to 5200 m, using in-situ observations for 15 glaciers in the Hunza Basin in 1959. Khan et al. [20] estimated an average value of ELA of 5000 m based on the average of the snow and clean-ice extents using Landsat images. In another study conducted by Racoviteanu et al. [36] an ELA of 4917 m to 5336 m was reported for the period 2000–2016 for the Hunza River basin by taking the average elevation of snow and ice-extents. Thus, the estimated ELAs in the current study and available studies are consistent and in good agreement considering differences in data, time and methods.

The estimated ELAs for July, August, and September during 2000 to 2017 have been used in time series trend analysis. The fundamental statistical parameters based on M–Kt and SSt i.e., Tau, p-value and SS are provided in

Table 1. Trend analysis results for monthly ELAs of the Hunza River basin.

Basin Name	Statistical Parameters of M–Kt and SSt	July	August	September
Hunza	Kendall’s tau (Ƭ)	0.111	−0.294	−0.163
	SMK	17	−45.000	−25.000
	p-value (Two-tailed test)	0.55	0.096	0.369
	Null Hypothesis (HO) at SL = 5%	Accepted	Accepted	Accepted
	Null Hypothesis (HO) at SL = 10%	Accepted	Rejected	Accepted
	Sen’s slope, SS	1.687	−7.106	−2.846

Note: H_o = No trend exists in time series data, and H₁ = Trend exists in time series data.

and shown in Figure 6.

The statistical parameters in Table 1 show that a positive trend exists in the time series dataset for the month of July; p-value 0.55 indicates that the trend is not statistically significant. On the contrary, August and September show a contrasting pattern of ELA with a decreasing trend, although the trend for September is statistically insignificant. Unlike July and September, the ELA trend for August is statistically significant at SL 10% (p-value is 0.096 < 0.1). The results of the August ELA trend show a recession of ELA at 7.1 m/year during 2000–2017. Thus, the ELA in August receded by about 128 m during the study period.

The decreasing trend in the ELA during August supports the “Karakoram Anomaly”, which refers to the observations that are in sharp contrast to those made in all other parts of the world where glaciers are shedding mass [9,12,37]. The “Karakoram Anomaly” has also been reported and agreed on by many researchers, based on their research on glaciers advancement, positive precipitation trends, negative temperature trends, positive mass balance studies, trend analysis of snow-cover areas, and declining flows trends during the summer months.

Of the available studies, the climate data have been assessed by many researchers. The decreasing trend in the ELA in August is in good agreement with the fact of observing significant positive trends in precipitation and negative trends in mean temperature in the summer by Fowler and Archer [38,39]. They conducted linear regression analyses of seasonal and annual temperature records from five meteorological stations within the UIB with data spanning the period 1961–1999. From these analyses, they concluded that winter and annual temperatures are rising in line with global warming but the summer temperature, particularly the minimum summer temperature, is falling. Khattak et al. [40] analyzed precipitation and temperature data for various stations in the UIB (including the Hunza Basin) during 1966–2005, and found similar results to Fowler and Archer [38,39].

Similarly, Hussain et al. [22] discovered that in high elevated areas of the Hunza River basin, pre-monsoon and monsoon seasons are complemented by decreasing trends in temperature. Trydtle et al. [41] also argued for an increase in precipitation in recent decades (the 1970s onwards) in the Hindukush–Karakoram mountain basins, using oxygen isotope measurements in tree rings. Atif et al. [1] found a positive trend in solid precipitation and a declining trend in temperature in the Hunza River basin. Similar trends have been reported by Reggiani et al. [42] for 1979–2014.

In addition, due to the scarcity of climatic stations at high altitudes in the Himalaya–Karakoram–Tibetan Plateau (HK–TP) region, long-term observations of atmospheric temperature acquired from the Microwave Sounding Unit (MSU) have also been used by several researchers to determine tropospheric temperature trends in recent decades [43,44]. In general, these data are for the period 1979–2008, and indicate a statistically significant rise in temperature for the winter to spring months (December to May), with relatively neutral or falling temperatures at a statistically insignificant rate for the summer months (June to November) over the Karakoram region [44].

The climate datasets analysis has further been supplemented by remotely sensed datasets. Forsythe et al. [45] noticed a decline in night-time cloud cover in the Hunza Basin, and argued that it is responsible for the declining summer minimum temperature, as cloud cover provides a radiative effect that can influence surface temperature. Tahir et al. [17] noticed an increasing trend in the snow cover in the Hunza watershed during 2000–2009. Hewitt [9,14,35,46] reported more than 34 surges in the Karakoram region since 1985. Rankl et al. [47] also noticed surging glaciers in the Hunza watershed during 1994–1999. Bolch et al. [5] provided advancement and retreat details of the glaciers in the HKH region. According to their study, most of the glaciers in the Karakoram (Hunza watershed) have been surging or have remained stable.

Similar to the above studies, the “Karakoram Anomaly” has also been attested by glacier mass balance studies. Glacier mass balance studies show slightly negative to positive mass balance in the western Karakoram (including the Hunza watershed) at a rate of + 0.11 ± 0.22 m/yr water equivalent (w.e), based on 850 kg/m³ density for ice/snow, during 1999–2008 [6], + 0.09 ± 0.18 m.w.e/yr during 2000–2008, using 850 kg/m³ density for ice/snow density [42], while Kaab et al. [7,8] noticed a −0.06 ± 0.04 m.w.e/yr (using 900 kg/m³ density for ice/snow) during 2003 to 2009. Muhammad et al. [48] also show similar results as shown by Kaab et al. [7,8] for the western Karakoram (Hunza basin) during the same time period.

Many researchers have also argued that positive mass balance and a declining ELA will decrease summer flows from the Hunza Basin. Fowler and Archer [38,39] noticed and argued that the Hunza Basin’s flow varies by about 16% per 1 °C rise/decline in temperature. They also noticed decreasing summer flow trends during 1966–1999. The decreasing summer flow trends (June, July, August, and September) have also been reported by Mukhopadhyay [49] during 1966 to 2010. The maximum decrease has been noticed during the month of August (16.3 m³/s per year). During summer months (June to September) the Hunza Basin contributes about 85% of its annual flows to the Tarbela Dam [49,50]. Therefore, a decline in flows from the Hunza Basin affects the reservoir operations of the Tarbela Dam. The Tarbela reservoir reached three times its lowest operation level during June 2019 due to low glacial melt in the month. Therefore, a reduction in summer flows warrants a review of the operation rules of the Tarbela Dam.

It is also noteworthy that the state and fate of the Karakoram glaciers have a profound influence on the river flows in the Upper Indus Basin (UIB), and will play a vital role in flows in a future changing climate. The projected climate in the Hunza Basin (Karakoram) is expected to be warm and wet. On a longer-term basis (to 2100) the warming climate may further enhance the temperature in the region [13,51]. According to Ridley et al. [52] and Wiltshire [53] it is argued that western disturbances will bring more and more frequent precipitation in winter to the Karakoram mountain range, which may maintain a positive glacier mass balance.

On one hand, the stable glaciers will store and provide sustainable water during future periods (as argued by Immerzeel et al. [51]), whereas on the other hand the summer flows decline, if continued, will require adequate water resources planning and management. In addition to this, since the 1860s, the Karakoram glaciers have experienced more than 34 surges, and these have become more frequent since 1985 [14,46], while Campbell [54] has identified 165 potentially dangerous glacier lakes across the whole HKH region, that could generate devastating floods. Thus, the advancement of glaciers in the Karakoram will bring both prosperity and threats in the near and far future periods, and require adequate planning and policy-making.

Thus, in sum, the current study findings are in good agreement and consistent with the available literature and field observations related to the historic ”Karakoram Anomaly”, and can be used in future glaciers’ health monitoring together with precise and accurate water resources planning and management.

4. Conclusions and Recommendations

The most glaciated watersheds of the UIB straddle the Karakoram Mountains and contribute a significant portion to the discharge volumes at the outlet of the basin (see [43] for various estimates). During the summer months (June to September), the Hunza Basin contributes more than 85% of its annual flows to the Tarbela Dam reservoir. The summer flows of the Hunza Basin have been reported to be declining during the summer months. The decrease in summer flows have been argued to be 16% per 1 °C decline. The main causes of the decreasing flows remain a decline in summer temperatures and an increase in winter snowfalls, resulting in positive glacier mass balance (the “Karakoram Anomaly”).

The available studies argue for an increase in snowfall in the Karakoram region together with an increase in temperatures. The future increase in snowfall and temperature will alter the state and fate of the glaciers. With the advent of technologies and accessibility to datasets, monitoring of the Karakoram glaciers using remotely sensed datasets and robust methods is a priority.

The current study, therefore, successfully used the remotely sensed MODIS data for the glacial health assessment, through estimation of ELAs during summer months in the Hunza Basin. Key findings of the current study are summarized below. Future sustainable, precise and accurate water resources planning and management depends on the state and fate of the glaciers in the HKH region. Therefore, the method developed and validated in the current study is recommended to be used in other parts of the HKH region.

• The MODIS snow-cover and albedo data can be successfully used in the monitoring of glaciers’ health, mass balance, and future water resources potential in the Karakoram and nearby similar mountain regions.
• The current study provides estimates for the period 2000–2017. The latest and future years’ data need to be assessed in future studies.
• The ELA trend for the month of August shows a decline in the ELA of 128 m at a rate of 7.1 m/year during 2000–2017.
• The results of this study are in good agreement with the available studies and successfully validate the “Karakoram Anomaly”.
• The current study results can be efficiently and effectively utilized in the Tarbela reservoir operation, management, future planning and policy-making.
• The current study argues that stable glaciers in the Karakoram may provide sustainable water resources concomitant with low flows in the current period and GLOFs in the future.
• To evaluate the sensitivity of the Hunza River basin’s ELA tow continuing climate changes, it is recommended to study and determine the correlation of the ELA, temperature, precipitation, flow and snow-ice extents.