Spatio-temporal changes in the six major glaciers of the Chitral River basin (Hindukush Region of Pakistan) between 2001 and 2018

Glaciers in the northern Pakistan are a distinctive source of freshwater for the irrigation, drinking and industrial water supplies of the people living in those regions and downstream. These glaciers are under a direct global warming impact as indicated in many previous studies. In this study, we estimated the glacier dynamics in terms of Equilibrium Line Altitude (ELA), mass balance and the snout position variation using remote sensing data between 2001 and 2018. Six glaciers, having area ≥ 20 km2 each, situated in the Chitral region (Hindukush Mountains) were investigated in this study. Digital Elevation Model (DEM) and available cloud-free continuous series of Landsat and Sentinel the entire study area was a retreat of -231 ± 140 m. No obvious relationship was found between the glacier variation trends and the available gauged climatic data possibly due to the presence of debris cover in ablation zones of all the studied glaciers which provides insulation and reduces the immediate climatic effects.


Introduction
Mountain glaciers are the earth's freshwater ecosystems and essential components of global hydrological cycles. Water resources are directly affected by the dynamics of glacier-melt runoff in the high mountain Asia. These areas including Himalaya, Karakoram and Hindukush (HKH) Mountains, have the major glacier deposits outside Polar Regions (Rowan et al. 2018). These glaciers are key water suppliers of major rivers in Asia (Ashraf et al. 2012;Immerzeel et al. 2009;Naeem et al. 2016). The climate of high mountain Asia imparts indirect effect on water availability in rivers downstream. Climatology of the HKH varies such as southern and eastern parts are mostly under the influence of monsoon while northern and western parts are influenced mostly by the westerlies (Hasson et al. 2017;Fowler and Archer 2006).
Glaciers are also considered as the prime climate indicators and helpful element for detecting rapid response towards climate changes. Glacier dynamics in terms of glacier parameters (mass balance, snout position, and Equilibrium Line Altitude (ELA)) depends on local and regional climate variability. Therefore, monitoring these glacial parameters over a certain period of time are the key elements to analyze the impact of climatic disturbances on the health of glaciers (Oerlemans 2005). Some researchers (Gardner et al. 2013;Jacob et al. 2012) have reported a significant glacier contribution to sea level rise which is stressing the need to better monitor the glacier mass change. Difficult and challenging access to high mountain Asian glaciers (highest glacier concentration outside Polar Regions) leads the way to adopt remote sensing techniques for glaciers monitoring (Philip and Ravindran 1998) in this region. Ali et al. (2017) suggested that correct estimation about snow and ice within a basin facilitate the mass balance assessments based on the observations of accumulation area ratio obtained through remote sensing. Rabatel et al. (2012) indicated that remote sensing is a valuable tool to estimate the ELA of a glacier. ASTER and SRTM DEMs were used by Bolch and Kamp (2005) in a study for glacier delineation and they used different techniques including Landsat Thematic Mapper TM4/TM5-ratio images, multispectral and landform analysis for glacier mapping approach. Band ratio method and Normalized Difference Snow Index (NDSI) method are compared for extracting glacial parameters in a study done by Wang et al. (2017) and their results showed that the band ratio method is better for glacier boundary extraction.
Many studies have been carried out to estimate the glacier size fluctuations in the HKH region using different datasets and approaches. Kääb et al. (2012) used ICESat data to estimate the mass balance of glaciers in the HKH region from year 2003 to 2008. They indicated that the specific mass balance for entire HKH region was −0.21 ± 0.05 m.a −1 water equivalent (w.e). Mass balance was extended between 2000 and 2016 by Brun et al. (2017a) for >90% of the glacierized area of High Mountain Asia using ASTER DEMs'. The mean mass change for HMA was −0.18 ± 0.04 m w.e. a −1 in the studied period, which is less negative than most of the previous estimates. Debris covered glaciers of Himalaya has increased mass loss and surface elevation lowering rather than recession of snout position (Bolch et al. 2011;Rowan et al. 2017;Garg et al. 2017). Berthier et al. (2007) used remote sensing data to supervise changes in glacier elevation and mass balances in Western Himalaya, India by comparing a DEM of 2004 from SPOT5 satellite to the 2000 SRTM DEM and results suggested that glaciers are experiencing rapid loss in ice accumulated area with specific mass balance of −0.7 to −0.8 m.a −1 w.e. Ghosh (2017) successfully investigated retreat of Gangotri Glacier using multispectral and multi temporal Landsat images.
The HKH glaciers need continuous assessment due to their regional importance in terms of major source of water availability and for future water management planning. There are fewer studies reporting the status of Hindukush glaciers due to lack of ground observatories in this rugged terrain . Previous studies lack the updated information about status of HKH glaciers and it is necessary to compare recent dynamics of glaciers with past changes (Muhammad et al. 2019a;2019b) to get the outlook of prevailing climate and environmental effects. Aim of the current study was to estimate the glacier changes on a long time-period in the poorly studied Hindukush region and its relation to regional climatic factors. In this current study, we present the long term (17-years; base year 2001) changes in different glacial parameters (in terms of ELA, snout position and surface elevation change) of major glaciers situated in the Chitral River basin (Hindukush region of Pakistan). Advanced Land Observing Satellite (ALOS) DEM, Landsat and Sentinel satellite data sets are used to carry out the study with a continuous series of satellite images from year 2001 to 2018 excluding 2011 and 2012.

Geographical Settings
Six glaciers (Figure 1) situated in the Chitral River basin of Hindukush region (north-west of Pakistan) were selected for this study. These glaciers are located between 7130′ to 7230′ E Longitude, and 36 to 37 N Latitude on the western portion of the Chitral River basin ( Figure  1). Glacier melt from the watershed vicinity contributes significantly to the flow of River Chitral as it runs down the valley. Highest mountain peak of the Hindukush region (Tirich Mir, 7708 m asl; Figure 1) is also situated in this river basin (Owen et al. 2002). There are many small and big glaciers in the Chitral River basin but we only selected six (mostly on western edge of the river basin) of those having an area greater than 20 km 2 . Key characteristics of the Chitral River basin are given in one of our previous study by Ahmad et al. (2018). Name or codes, geographical location, area and elevation ranges (on the basis of DEM and Randolph Glacier Inventory (RGI, version 6; RGI 2017)) of the six studied glaciers are presented in Figure 1 and Table 1.
Digital elevation models derived from ALOS DEM are superimposed by contour lines and the hypsometric curves of the six studied glaciers are presented in Figure 2 and Figure 3, respectively. Glacier ID05 has the largest area among all the 7studied glaciers covering 90 km 2 followed by glacier ID04, ID03, ID02, ID06 and ID01. Glacier  ID04 has the highest mean elevation followed by the glaciers ID06, ID05, ID02, ID03 and ID01.

Remotely-sensed satellite data
Main remotely sensed satellite data sets used in this study are the DEM and the multispectral LANDSAT and Sentinel images. DEM tiles of ALOS PALSAR were downloaded from the Alaska

Satellite
Facility's Vertex Data Portal (http://vertex.daac.asf.alaska.edu/). ALOS DEM has a resolution of approximately 12.5 m and was used for watershed delineation in this study. DEM data was also used to study the topographical characteristics (hypsometry, contour lines and elevation variation) of studied glaciers which were then used for the ELA estimation.
The Thematic Mapper (TM) instrument of Landsat is available since 1982 that responds the reflected radiation from the surface of earth in the (MSI) which has important characteristics for glaciers investigations Paul et al. 2016). Sentinel data can be downloaded from the online portal of Copernicus and European Space Agency (ESA).
Available images of Landsat (TM, ETM and OLI) and Sentinel 2 (S2) satellites for the study area from the year 2001 to 2018 with less than 10% cloud cover and minimum snow cover season were obtained from the United States Geological Survey (USGS) Earth Explorer website (http://earthexplorer.usgs.gov/) and ESA online portal, respectively. Scene IDs, date of acquisition, sensor name and resolution of downloaded images are presented in Table 2. These images were used for the glacier dynamics after the standard pretreatment procedure such as atmospheric correction, conversion of digital numbers to reflectance values and co-registration.

In-situ data
Mean temperature (Tavg) and total precipitation (PT) data of Chitral weather station (latitude 3551 and longitude 7150; altitude: 1498 m asl) is obtained from Pakistan Meteorological Department (PMD). PMD maintains the sequential data of climate stations all over Pakistan. The data has been used for climatic trends assessment to relate it with the obtained glacier dynamics results. Climatic variables data available for approximately 24-years (19902013) was used for trend analysis. Climatic data used is from lower elevation as compared to studied glaciers, therefore, it may have a magnitude difference of Tavg and PT but the trends may still be similar as those on higher elevation due to close proximity of the station.

Extraction of the studied glaciers and their topographical characteristics from the ALOS DEM
Shapefile of glacier boundaries obtained from Randolph Glacier Inventory (RGI) version 6.0 (RGI 2017) was used as a mask to extract the area of six studied glaciers ( Figure 1) from the ALOS DEM of Chitral River basin to study the topographic characteristics such as elevation, contour lines and hypsometry as shown in Figure 2 and 3. Different researchers (Porter 1975;Thomson et al. 2016) have reported contouring method with varied contour intervals for glacier dynamics studies. Contour lines were mapped ( Figure 2) for the studied glaciers with a contour interval of 400 m from the extracted DEM of each glacier. The contour lines were used further for the analysis of change in the ELA for the studied glaciers.

Pre-treatment of the LANDSAT and Sentinel satellite data
Scan Line corrector (SLC) onboard LANDSAT satellite sensor (ETM) that work for the forward motion of LANDSAT 7, suffered major failure in 2003 (Maxwell et al. 2007). The gaps of LANDSAT 7 (ETM) images due to SLC failure are ranging up to 14 pixels near the boundaries with 22% loss of data (Zhang et al. 2007). QGIS 2.18.13 software has been used to fill the gaps of LANDSAT 7 images. The data set downloaded from USGS has the raster band files as well as the gap data files of each band image that are used for gap filling.
After the acquisition of data and gap filling (for LANDSAT 7), the atmospheric correction has been done for LANDSAT and Sentinel images to reduce the effects caused by sensor, solar radiations, atmospheric and topographic factors according to the method used by Fraser (1992) and Azabdaftari and Sunar (2016). For the atmospheric correction, firstly, the data in digital numbers has been converted to radiance data using Eq.(1). Radiance data thus obtained was further converted to reflectance data using Eq.(2). Reflectance data sometimes contain minor negative values which are false and generally forced to be positive by using Eq. (3).
Rλ= reflectance (unitless ratio), Lλ= radiance calculated in Eq. (1) These atmospherically corrected images were then co-registered using the image of the year-2001 as a reference image used for further analysis to study the glacier dynamics parameters (ELA, mass balance and snout position).

Estimation of the ELA using band spectral indices, band ratio and band color composite methods
If measured at the end of the melt season, the snowline altitude (SLA) is concurrent with the ELA. The equilibrium line marks the position where accumulation of snow in accumulation zone is exactly balanced by the ablation taking place from the ablation zone for the period of one year. The ELA has been estimated for each selected glacier using LANDSAT, Sentinel images and DEM data. Mir and Majeed (2016) suggested that the manual digitization is more accurate than automated methods for glaciers with debris cover. Different spectral techniques applied on images in this study aided to manually delineate ELA position on glaciers. Normalized difference snow index (NDSI), Normalized difference glacier index (NDGI), band ratio and false color composite (NIR, Red, Green) techniques were used to enhance the contrast between the accumulation and ablation area over glaciers. Band ratio method and false color composite technique were observed to be better for glacier dynamics studies (Wang et al. 2017). Normalized difference debris index (NDDI) was used to analyze the debris covered portion on glaciers to better analyze the actual accumulation area of glaciers. Formulas used for NDSI, NDGI and NDDI are presented in the Table 3. The results of all above mentioned techniques are compared and analyzed to obtain ELA of each glacier. By using the above mentioned techniques with suitable threshold values and contour lines obtained from DEM, we have digitized the equilibrium line position on each glacier for each available image. The mean ELA of each glacier are calculated as the average of three observations taken on the ELA for each year.

Estimation of glaciers elevation change
Estimation of glaciers elevation change was derived from ASTER DEM differenced data available at PANGAEA (Brun et al. 2017b) using two ASTER DEM's from year 2000 and 2016. Geodetic DEM differencing method provides reliable mass balance estimates for areas of few kilometers (Bhattacharya et al. 2016). Surface elevation changes derived from DEM differenced images (Brun et al. 2017b) were converted into mass balance using an average ice density of 850 kg/m 3 as also suggested by Huss (2013) for longterm mass balance studies. The method, data, and ice density used to calculate the mass balance in this study is similar as by Muhammad et al. (2019a).

Estimation of the snout position
Snout position analysis is significant to assess glaciers' retreat or advancing behavior in response to climate change. The analysis is important to estimate future contribution in stream flow as result of the glacier melt. Rate of melting at the terminus can be influenced by debris cover. Excessive debris prevents melting and protects a glacier from recession (Kulkarni et al. 2005). It is comparatively easy to map clean glaciers than debris covered glaciers in which the mapping is inhibited by the presence of debris. Therefore, different feature analysis techniques has been used in previous studies to map the glacier snouts having debris accumulation i.e. the presence of moraine dammed lakes, ice-wall shadow, texture of glacier surface and periglacial area, braided streams, moisture content of supraglacial debris, and lateral moraines Shukla and Qadir 2016). In this study, snout position for each glacier has been analyzed comparing all the images obtained from different techniques (NDSI, NDGI, NDDI, Band ratio and Band composite) over entire study period. The change in the position of snout has been measured using the measure tool in ArcGIS (version 10.5) that calculated the advance and retreat in meters by keeping the reference snout position as observed in year-2001. The advance has been taken as positive and retreat as negative measures from reference year-2001. Manually digitized snout positions of each glacier have been saved as separate shapefiles. Presence of water streams coming out from the snout area of glaciers, presence of debris and side wall calves and texture facilitated in snout position mapping.

Climatic trends assessment
The climate data (mean temperature and total precipitation) taken from PMD has been used to obtain linear trends of local climate on the study area for 24 years because the trends of climate are generally in correlation with glacier dynamics. Long-term patterns of climatic factors were investigated to build up an understanding of prevailing and future climatic trend and their link with observed glacier changes. The trends were estimated on the mean annual temperature (Tavg) and the total annual precipitation (PT). XLSTAT data analysis software available as Microsoft Excel plugin (Chang et al., 2017) was used for the statistical applications. Nonparametric Mann-Kendall (MK) trend test (Hirsch and Slack 1984;Kendall and Gibbons 1990;Yue et al. 2002) and Theil-Sen slope estimator (also known as Sen's slope) (Gilbert 1987) for periodic data were used to estimate the trends (significance level = 5%) in the time-series of available mean temperature, and total precipitation. Linear trend line equation is also presented for all the trend analyses in this study. Sen's slope method has been reported in numerous researches for trends analysis in hydro- climate variables (Teegavarapu and Nayak 2017).

Dynamics of the Equilibrium Line Altitude (ELA)
Above or below shifting of annual ELA from referenced ELA position (2001 for current study) is a key indicator of the state of health of a glacier. Application of techniques as mentioned in section 3.3.3 aided in enhancement of ELA between accumulation and ablation area which was then digitized manually for each studied glacier. The best results are obtained using the combination of band composite (NIR, Red, Green), NDSI, NDGI and band ratio (NIR/SWIR) methods. The range of the applied threshold values obtained for all the spectral indices and band methods are given in the Table 4.

Dynamics of glaciers surface elevation
Glaciers surface elevation change was carried out as detailed in section 3.3.4. Results of the glaciers surface elevation change from 2001-2016 for each glacier using data provided by Brun et al. (2017a), similarly as derived by Muhammad et al. (2019a) are presented in Figure 6. An average geodetic mass balance obtained for the entire investigated study area comprising of six glaciers is -0.110.07 m w.e. a -1 . We obtained a mass balance value of −0.14 m, -0.17 m, -0.19 m, −0.06 m, −0.08 m and −0.02 m w.e. a −1 for glaciers ID01, ID02, ID03, ID04, ID05 and ID06, respectively ( Figure 6).
It has been observed that the overall surface elevation loss of studied glaciers is slightly negative. For the glaciers ID01, ID04, ID05 and ID06, the observed elevation change during the studied period (2001-2016) is −0.14 m, -0.061 m, -0.08 m and −0.02 m, respectively, that shows the receding status of these glaciers. Glacier ID06 has very low value of elevation change and may be termed as slightly receding. For the glaciers ID02 and ID03, the surface elevation change values are calculated as -0.17 m and -0.19 m, respectively, which shows comparatively significant surface elevation lowering as compared to other studied glaciers (Muhammad et al. 2019b).

Dynamics of the snout position
Snout position of studied glaciers has been digitized using treated LANDSAT and Sentinel band composite images after suitable threshold values (Table 4)

General trends in the glaciological parameters (ELA, mass balance and snout position)
General trends obtained after the investigation of different glaciological parameters in this study are summarized in Table 5. Trends among the entire studied glacier are almost consistent i.e. glaciers are receding at snout and there is an upward shift of ELA reducing the accumulation area (Table 5). Sarikaya et al. (2012) reported the reason of observed glacier fluctuations (as observed in this study) in the Hindukush region as more  strongly influenced by precipitation than by temperature variations, alongside the significant dependence on microclimatic conditions governed by topography. There is a slight discrepancy amongthe glaciers health on the basis of glaciers surface elevation change. All the studied glaciers have negative surface elevation change. GlacierID02 and glacierID03 have the highest surface elevation loss as assessed through DEM differencing method. ELA shift and snout retreat has also been observed as maximum in glacierID02. GlacierID04 has shown advancing behavior in snout position which is consistent with the observed surface elevation change that calculated as less negative as compared to other studied glaciers.

Climate trends assessment
A trend analysis was conducted on the 24 years time-series  of climate data (mean temperature, Tavg and total precipitation, PT) to assess general climate pattern prevailing in the study area as shown in Figure 9.
Although the estimated trends were nonsignificant for both the variables but they have slight tendency towards a decrease. Mann Kendall trend values of −0.087 and −0.188, and the Sen's slope values of −0.011 o C.a -1 and −4.94 mm.a -1 were found for the mean annual temperature and the total annual precipitation, respectively. Trends of annual precipitation and mean temperature from Chitral basin may be regarded as stable during the study period. A seasonal trend analysis by the Ahmad et al. (2018) has revealed a significant increase in the mean temperature which may be the reason of gradual snout retreat, loss of mass balance and ELA rise in glaciers.

Discussion
Our observed range of ELA values for reference year is in agreement with the values given by Haserodt (1989) for the Chitral region. Figure 5 is demonstrating the ELA fluctuation on studied glaciers during study period which shows the overall increasing trend (upward shifting) of ELA on all glaciers as shown by positive slope values. Result of our conducted study shows that  This is consistent with the estimated ELA by Scherler et al. (2011a) and Gardelle et al. (2013) in Hindukush region.
The slight difference in ELA values reported by other researchers may be due to the different study period. Our studied glaciers show varied trends of ELA among themselves which may be explained by the statement of Vashisht et al. (2016) who suggested that different behaviors of ELA shift observed in same area may be due to the micro climatic conditions prevailing in that area. The studied glaciers that have experienced regular upward shift in the ELA may be due to the increased avalanche events triggered by steep topography (Azam et al. 2018) and the regular earthquake events in the area. As avalanching, rockfalls and landslides from adjacent rocks (triggered by freeze-thaw process or earthquakes) contribute to the glacier surface debris (Nagai et al. 2013;Markus and Christian 2012;Miles et al. 2017) which cause rapid melting in some cases.
The mean mass balance observed (-0.110.07 m w.e. a -1 ) value is in agreement with the mass balance value (-0.12±0.07) calculated by Brun et al. (2017a) for the Hindukush region using a different time series (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) of DEMs. Results showed that the mass balance values in Hindukush region are slightly less negative than those found in the Himalayas by other researchers. For example, Berthier et al. (2007) reported the range of specific mass balance as -1.4 m to +0.1 m w.e. a -1 for the glacier in western Himalaya but their findings suggest the overall negative mass balance of glaciers in the region which resembles our results of negative mass balance for most of studied glaciers. Pellicciotti et al. (2017) reported a mass balance for debris covered glaciers of Langtang Himal, Nepal as -0.320.18 m w.e. a -1 ). The decreasing mass of studied glaciers indicates that the trend may be similar in the entire Hindukush region but it needs to be explored in detail.
The snout analysis of glacierID05 has revealed that the glacier has experienced slight retreat as compared to other studied glaciers while glacierID04 has shown a slight advance in snout position during study period. There may have been the effects of local climatic conditions and topography of this glacier. The glacierID04 has the second largest surface area among all the six studied glaciers so the impact of climatic changes on the snout may appear after a long time (not significant in our study period). Most of the ablation areas including the snout in the studied glaciers are debris covered (Figure 7) which results in a delayed response to climatic variability. Many previous studies have also witnessed that the glaciers having thick debris cover have unpredictable response towards temperature increase due to the insulation properties of debris (Salerno et al. 2017;Vincent et al. 2016). The observed snout position fluctuation may also be due to the geomorphic parameters in the area (Kulkarni and Karyakarte 2014;Benn et al. 2012). Scherler et al. (2011b), Bahr et al. (1998) and Mir et al. (2017) stated that the glaciers having debriscover lose most of their mass by decreasing surface elevation rather than snout retreat and small glaciers are generally more prone to climatic changes because of their shorter response time to climatic factors, similar to the trends observed in some of our studied glaciers.

Conclusions
The current study presents analysis of glacier dynamics between 2001 and 2018 using remote sensing data for six glaciers (each with an area ≥20 km 2 each) situated in the Chitral river basin of Hindukush region. Glacier variations were studied based on standard glaciological parameters such as ELA, snout position and the surface elevation change. Following can be concluded from the results obtained in this study.
 ELA has shifted upward from its reference position in 2001 during the data period of 19 years at an average rate of ~16 m.a -1 for the entire study area indicating a loss in the accumulation area. All the studied glaciers on average show an upward shift of ELA with varying rates during the study period.
 There is overall loss observed in surface elevation of studied glaciers during the data period. Glacier ID04 and ID06 show a less negative surface elevation loss but the general trend over entire studied area is negative.
 There is a general retreat in the snout positions with a mean retreat of ~130 m for the entire studied glacier area during the study period.
Variation in responses of glaciers parameters may be the result of local climatic condition of each glacier as well as difference in topographic features and location. All the studied glaciers are indicating a generally gradual upward shifting of ELA, mass balance loss and the retreat of snout position. A set of threshold values found during the estimation of ELA and snout position using NDGI, NDDI and NDSI during this study may be further used for this region and similar data for future glaciological studies. No significant trends in climatic data of the studied area were found but previous studies on seasonal data suggested a warming trend in this area which may explain the reason of spatial changes in glacier parameters. The overall trends may be similar on other glaciers of the Chitral River basin but needs more detailed field and remote sensing studies in the region.

Acknowledgment
Financial support for this research work by the National Natural Science Foundation of China (NSFC) and ICIMOD (Grant no. 41761144075) is highly acknowledged. Authors thank the United States Geological Survey, European Space Agency and Copernicus for providing the Landsat and Sentinel satellite data free of cost for this study. The authors extend their thanks to the Water and Power Development Authority (WAPDA) and the Pakistan Meteorological Department (PMD) for contributing their hydrological and meteorological data, respectively. We are grateful to the anonymous reviewers whose comments and suggestions significantly improved the manuscript.

Electronic supplementary material:
Supplementary material (Appendix 1) is available in the online version of this article at https://doi.org/10.1007/s11629-019-5728-9 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by/ 4.0/, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons License, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons License and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this License, visit http://creativecommons. org/licenses/by/4.0/.