The satellite observed glacier mass changes over the Upper Indus Basin during 2000–2012

Decadal glacier thickness changes over the Upper Indus Basin in the Jammu and Kashmir Himalaya were estimated using the TanDEM-X and SRTM-C Digital Elevation Models (DEMs) from 2000 to 2012. In the study area 12,243 glaciers having 19,727 ± 1,054 km2 area have thinned on an average of − 0.35 ± 0.33 m a−1 during the observation period. The highest thinning of − 1.69 ± 0.60 m a−1 was observed in the Pir Panjal while as the marginal thinning of − 0.11 ± 0.32 m a−1 was observed for the glaciers in the Karakoram. The observed glacier thickness changes indicated a strong influence of the topographic parameters. Higher thickness reduction was observed on the glaciers situated at lower altitudes (− 1.40 ± 0.53 m a−1) and with shallower slopes (− 1.52 ± 0.40 m a−1). Significantly higher negative thickness changes were observed from the glaciers situated on the southern slopes (− 0.55 ± 0.37 m a−1). The thickness loss was higher on the debris-covered glaciers (− 0.50 ± 0.38 m a−1) than on the clean glaciers (− 0.32 ± 0.33 m a−1). The cumulative glacier mass loss of − 70.32 ± 66.69 Gt was observed during the observation period, which, if continued, would significantly affect the sustainability of water resources in the basin.

www.nature.com/scientificreports/ and 25% 24 of the glacier thickness changes. Salerno et al. 21 on the other hand, indicated that the slope of a glacier tongue is the main topographic parameter controlling the glacier thickness changes. The influence of debris cover on the Himalayan glaciers is not very clear. Several geodetic mass balance studies have reported similar thinning rates for both debris-free and debris-covered glaciers 28 contradicting the reports of the slowdown of glacier melting under thick debris-cover 29 . The present study investigated the changes in glacier thickness based on two publicly available SAR-based Digital Elevation Models (DEMs) from 2000 to 2012 over the UIB in the Jammu and Kashmir Himalaya, India (Fig. 1). The study area, comprising of Jammu, Kashmir and Ladakh regions, is spread over an area of 222,236 km 2 and nearly 11% of the geographical area is covered by glaciers (glacier count = 15,064 spread over 24,022 km 2 ) 30 . Majority of the glaciers in the UIB are valley glaciers except in the Ladakh mountain range where cirque glaciers are more dominant 31 . The study area is often divided in six mountain ranges: Pir Panjal range, (PPR), Greater Himalaya range (GHR), Shamaswari range (SR), Zanaskar range (ZR), Ladakh range (LR) and the Karakoram range (KKR) each with distinct climatic and topographic characteristics. The climatology of the region is largely dominated by the western disturbances (WDs) compared to the monsoons which are predominant in the Indo-Gangetic plains 32 . The WDs, mostly originating from the Mediterranean Sea, are the major sources of precipitation particularly during winters in the study region 33 . Glaciers in the western Himalaya encompassing the study region receive 60-70% of their annual accumulation from the WDs 34 .
Data set. The Shuttle Radar Topography Mission (SRTM) operated between 11 and 22, February, 2000 and acquired Interferometric Synthetic Aperture Radar (InSAR) data from the space in two frequencies (X-band and C-band). However, compared to the C-band, the X-band system has limited coverage 35 . The 1-arc sec SRTM/X-SAR (X-SRTM) DEM is available from the German Aerospace Centre (DLR) and the non-void filled 3-arc sec SRTM/SIR-C DEM is distributed by the United States Geological Survey (USGS).
The TanDEM-X mission (TerraSAR-X-Add-on for Digital Elevation Measurements) with its twin satellites was launched in 2010 by the DLR. TanDEM-X and TerraSAR-X, flying in helix formation, record the backscattered signal with the satellite overpass lag as small as 3 s. The very small temporal baseline makes the data suitable for interferometric processing 36 . We used the TanDEM-X DEM composite data from 2010 to 2015 period and is The glacier outlines, used in this study, are based on the globally complete Randolph Glacier Inventory (RGI) version 6.0 30 . Furthermore, the Supra-glacial Debris Cover Dataset v1.0, a global supraglacial debris-cover dataset was used for glacier debris-cover (DC) characterization 37 . Keeping in view the limitations of the debris-cover dataset, we manually corrected the discrepancies in the data for glaciers with debris-cover fraction > 19% using Landsat data. We also used the MODIS LST (MOD11A2) to characterize the climatic settings of different mountain ranges in the study region. The details of the dataset used in this study are given in Supplementary Table S2.

DEM corrections and elevation changes.
In order to remove the horizontal and vertical offset between the two DEMs, the universal co-registration algorithm was used 38 (details in the Supplementary Section 1). Prior to their use for thickness change estimation, the DEMs were also corrected for radar penetration bias. The co-registered DEMs were differenced to generate the elevation difference (dH/dT) map over the glaciated terrain at pixel level. RGI6.0 glacier outlines 30 were used to calculate the mean glacier elevation changes between 2000 and 2012 and the volume changes thereof. Using the density conversion factor of 850 kg m −3 , the volume changes were then converted into glacier mass changes 39 (details in the Supplementary Section 2). It is pertinent to mention that only the glaciers with dH/dT coverage > 30% were considered in the present study. Furthermore, the uncertainties in the glacier mass changes owing to the uncertainty in DEM differencing, radar signal penetration, uncertainty due to void fill, glacier outlines and mass conversion have been addressed separately and discussed in the Supplementary Section 3. The glacier thickness change estimates are based on the elevation differences in the DEMs obtained ~ 12 year apart. Though, the SRTM DEM was obtained over a shorter period of time (11-20 February, 2000) but the timestamp of each TANDEM-X acquisition is not same and is spread over a wider period. This has the potential to add to the uncertainty of glacier thickness and mass changes which has been considered in the uncertainty analysis in the present study 40 . The bias correction and uncertainty analysis is discussed in detail in the Supplementary Sections 3 and 4.

Debris categorization and glacier topographical parameters.
To assess the influence of supraglacial debris cover, we used two criteria to differentiate between the clean and debris-covered glaciers: one proposed by Brun et al. 9 using > 19% debris-cover fraction threshold (Criterion 1) for identifying debris and non-debris glaciers and the other proposed by Ali et al. 41 which categorizes glaciers into three categories: clean glaciers with the debris-cover fraction < 25%; sparsely debris-covered glaciers with debris-cover fraction between ≥ 25 and ≤ 50% and debris-covered glaciers with debris-cover fraction ≥ 50% (Criterion 2). Prior to the use for analysis, the debris-cover dataset was corrected manually for any discrepancy using Landsat satellite images, however, we only corrected the glaciers with debris-cover fraction > 19%. The topographic parameters like elevation, slope and aspect for each glacier in the study area were extracted from the TanDEX-X DEM in ArcGIS environment (details in the Supplementary Section 4).

Results
Glacier thickness and mass changes. The investigation showed that the glaciers in the UIB were thinning at an average rate of − 0.35 ± 0.33 m a −1 (Table 1) amounting to the glacier mass loss of 297.5 ± 280.5 kg m −2 a −1 . The mean glacier-wide thickness changes in the UIB varies from ~ − 5.0 to ~ 16.0 m a −1 (Fig. 2), however, it is pertinent to mention here that ~ 64% (count) of the glaciers, accounting for ~ 73.3% of the glacier area in the UIB, fall in the 0 to − 2.0 m a −1 thickness change category ( Supplementary Fig. S3). The thickness change of − 0.11 ± 0.32 m a −1 is significantly lower in the KKR. The highest average glacier thickness loss of − 1.69 ± 0.60 m a −1 was observed in the PPR. The average glacier thickness loss of − 1.28 ± 0.46 m a −1 and − 1.12 ± 0.40 m a −1 was observed in the SR and GHR respectively. The glaciers in the ZR and LR have, on average, thinned − 1.17 ± 0.41 m a −1 and − 0.46 ± 0.26 m a −1 respectively (Table 1). Table 1. Changes in glacier thickness across different mountain ranges of the study region. Only the glaciers with < 30% voids corresponding to an area of 19,727 km 2 equivalent to ~ 82% of the total glacier area (24,022 km 2 ) were considered for the analysis.  www.nature.com/scientificreports/ Topographic influences. The investigation revealed a strong correlation between the mean glacier elevation and thickness change (R = 0.64). The influence of elevation on the glacier thickness change is evident from the fact that the highest glacier thickness loss (− 0.86 ± 0.56 m a −1 ) was observed in the glaciers located at mean elevations < 5,000 m a.s.l. and the rate of thickness loss is significantly lower (− 0.42 ± 0.33 m a −1 ) in the glaciers located at mean elevations > 5,000 m a.s.l. (Supplementary Fig. S6). It is pertinent to mention that ~ 25% and ~ 75% of the total glacier area is distributed at elevations < 5,000 m a.s.l. and > 5,000 m a.s.l. respectively. The analysis also indicated a strong positive correlation (R = 0.87) between the mean glacier slope and thickness change, i.e., the thickness loss is less for glaciers situated on steeper slopes. The glaciers with the mean slope < 10° have thinned at the rate of − 1.52 ± 0.40 m a −1 while as, mass gain of 0.28 ± 0.40 m a −1 was observed in the glaciers with the mean slope > 30° ( Supplementary Fig. S7).  Table S6). The significant impact of the debris-cover on the glacier thickness change is also evident when using the Criterion 1 for the debris-cover categorization. The clean and debris-covered glaciers based on the criterion 1, comprising ~ 92% and ~ 8% of the glacier area, have witnessed the thickness change of − 0.35 ± 0.33 m a −1 and − 0.53 ± 0.35 m a −1 respectively (Supplementary Table S6).
There is a positive, though weak, relationship (R = 0.25) between the glacier size and thickness changes. In general, the smaller glaciers experience higher thickness losses compared to the larger glaciers 42 . The larger glaciers with the size > 50 km 2 , comprising ~ 36% of the glaciated area under consideration in the UIB, have shown the lowest thickness changes (− 0.32 ± 0.26 m a −1 ). On the contrary, the glaciers with area < 50 km 2 comprising ~ 64% of the total glacier area have thinned, on an average, − 0.43 ± 0.26 m a −1 during the observation period (Supplementary Table S7 Table S8). The relatively higher elevation changes observed over the PPR and the lowest elevation changes in KKR are in agreement with the prevalent climatic regimes in the two ranges (Table S8). Similarly, the elevation change pattern observed in the GHR, SR, ZR and LR is also in line with the prevalent climatic settings observed across these mountain ranges. The details of the average summer and winter temperatures observed over different mountain ranges in the study region from 2000 to 2014 are provided in the Supplementary Table S8.
Validation of the results. Field observations of mass balance over the study region are very scarce. In fact only 4 glaciers over the study region have been studied for mass balances 43 . These observations are, however, very short (one year mass balance for the Kolahoi and Shishram, 2 years for the Rulung and 8 years for the Nehnar glacier respectively) and restricted to 1980s and as such were not considered for the validation in the present study. The elevation changes observed in the present study were, however, validated against the elevation changes based on the High Mountain Asia (HMA) DEM with specific date-stamps and SRTM DEM differencing 12 . The HMA-SRTM elevation changes for the selected glaciers fall within the uncertainty bars of the TanDEM-X-SRTM derived elevations changes observed in this study. The spatially distributed elevation changes of these glaciers along with the GLIMS ID are provided in Supplementary Fig. S10.

Discussion
Topographical parameters, to a large extent, explain the thickness change variability observed across the six mountain ranges in the UIB. This is due to the fact that there is a marked difference in the topographical variables across the ranges which is strongly related to glacier response time and sensitivity of mass balance to climate change 18,21 . The glaciers situated on the gentle slopes have relatively longer response time, leading to the slow glacier dynamics 44 . They do react immediately to warming by retreating, however, they remain in equilibrium by making dynamic adjustments to climate forcing and therefore remain out of balance for a longer time 45 . The relationship between the glacier slope and the observed glacier thickness change is strong (R = 0.76).
The glacier aspect has a profound effect on the glacier melting 20,46 . In the Himalaya, the north facing slopes generally receive less solar radiations compared to the south facing slopes 47 . This explains the robust relationship (R = 0.86) observed between the aspect and the glacier thickness changes in all the six ranges of the basin. Relatively higher melting rates observed for the glaciers situated on the warmer southern slopes have been previously reported in the Himalaya 47 .
The glaciers situated at lower altitudes are more sensitive to the rising temperatures, as these glaciers tend to have larger mass turnover to balance the relatively higher melting in their lower reaches 48 . This explains the influence of glacier elevation on the observed ice thickness changes as evident from the high correlation coefficient (R = 0.78). Similar relationship has been previously reported in the Himalayan region 9,49 . Contrary to the common sense, the glaciers with mean elevation > 6000 m a.s.l. have shown higher ice thickness losses compared to the glaciers situated at lower altitudes between 5,000 and 6,000 m a.s.l.. This asymmetry is explained by the fact that nearly ~ 35% of the total glacier area above 6,000 m a.s.l. is south-oriented compared to ~ 23% of the south-oriented glaciers situated between 5,000 and 6,000 m a.s.l. elevation. The shallow mean slope of the glaciers situated between 4,000 and 5,000 m a.s.l. elevation explains the slightly higher observed ice thickness recession compared to that of the glaciers located at elevations below 4,000 m a.s.l. (Supplementary Fig. S6).
The higher ice thinning (− 1.69 ± 0.60 m a −1 ) observed in the PPR correlates with its lower mean altitude compared to the other mountain ranges in the UIB ( Table 1). The relatively higher glacier thinning observed in the SR despite its higher mean elevation compared to that of the glaciers in the GHR is explained by the presence of the relatively more supra-glacial debris of the glaciers in the SR range ( Table 1).
The variable glacier thickness changes observed in the climatologically somewhat similar LR and ZR ranges 50 situated in the cold desert region of Ladakh are also attributed to the differential topographic and morphological settings of the glaciers in these two ranges. The glaciers in ZR have relatively lower mean altitudes (5,032 m a.s.l.), relatively more area distributed on the southern slopes (~ 26%) and relatively higher (~ 13%) debris-cover (Table 1). All these parameters favour enhanced glacier melting and therefore explain the relatively more thinning of the ZR glaciers compared to the glaciers in the LR. Almost zero glacier thickness changes observed in the KKR are in line with the previous studies in the mountain range 8,10,[13][14][15] . The stability of the glaciers in the KKR, despite the fact that the morphological and topographical variables (Table 1) in general seem to favour enhanced melting, corroborate with the concept of the preponderance of climatic influence on the KKR glacier dynamics [51][52][53] .
The influence of the glacier morphological and topographic variables on glacier mass balance has been previously reported by Huss 27 , who concluded that 35% of the glacier mass balance variability is explained by the median glacier elevation, mean glacier tongue slope and glacier area. Similarly, Rabatel et al. 24 reported that 25% of the glacier mass balance variability is explained by glacier median elevation and mean slope of the glacier tongue.
Among the non-climatic parameters, supra-glacier debris is considered as one of the important parameter influencing glacier mass balance 21,54 . Supra-glacial debris cover alters glacier surface energy balance acting as a barrier between the atmosphere and the ice. It can lead to the reduction of melt, but in case of thinly debriscovered glaciers or the glaciers with patchy supra-glacial debris, melt rates are enhanced compared to the bare ice 55 . Differential thinning of debris-covered and clean glaciers is reported in several other studies 26,56 , however, the picture is not very clear as the debris influence on the melting of glacier varies as a function of the debriscover extent, distribution and thickness 42 . Some studies have suggested reduced melting of the debris-covered glaciers 56,57 , and contrarily, a few other researchers have reported enhanced melting from the debris-covered glaciers 58,59 . Even some studies have reported similar thinning rates from both debris-covered and debris-free glaciers indicating no overall influence of the debris-cover on mass balance 28  www.nature.com/scientificreports/ or less in tune with the studies suggesting relatively higher thinning rates from the debris-covered glaciers which is corroborated by the higher mass balance sensitivity of the debris-covered glaciers to the rise in temperature 60 . However, the influence of debris-cover on the glacier thickness changes is not uniform; out of the 1,011 debriscovered glaciers (based on the criterion 1), 275 glaciers showed stability or even a slight gain in thickness, while as 736 glaciers showed the negative thickness changes. The observed variability in the glacier thickness changes from the debris-covered glaciers is explained by the variability of topographic variables; the glacier with positive or no-change in the ice thickness are situated at higher mean elevation (5,431 m a.s.l.) and have steeper slope (~ 29°) compared to the glaciers showing the negative thickness changes which have a mean elevation of 5,228 m a.s.l. and ~ 24° slope ( Supplementary Fig. S9). Larger glaciers are laggard in responding to climate perturbations because of their long response time, with the consequent lower retreat rates 42 . The results indicate a slight positive correlation between the glacier size and the glacier thickness changes, however, the relationship is not significant or even uniform across the different glacier size categories (Supplementary Table S7). The observed variation in the glacier thickness for different glacier-size categories is largely controlled by the topographic variables. For instance, the glaciers with area < 1 km 2 have thinned relatively less due to their location at the higher mean elevations (5,292 m a.s.l.). The observed lower thinning rates of the smaller glaciers is corroborated by the fact that the smaller glaciers across the mountain ranges in the UIB are situated at higher altitudes, mostly in cirques and below the rock cliffs where the accumulation is significantly pushed up by snow avalanches and wind drift 61,62 , resulting in the slowdown of the glacier thinning. Generally, the lower mean elevation, shallower slopes, mean southern aspects and relatively higher debris-cover fraction of the glaciers (Supplementary Table S7 The differential elevation changes observed across different mountain ranges were also found in good agreement with the prevalent climatic regimes over the mountain ranges (Supplementary Table S8). It is pertinent to mention here that Shekhar et al. 63 based on the analysis of 18 meteorological stations distributed over different mountain ranges of the study region reported an increase in both the maximum and minimum temperatures except for the KKR. The study reported an increase of 0.8, 2.0 and 1.0 °C in the maximum temperature in the PPR, SR and GHR respectively from 1988 to 2008. The minimum temperatures during the same period increased by 0.6, 1.0 and 3.4 °C in the PPR, SR and GHR mountain ranges respectively. However, the study reported a decrease of around 1.6 and 3.0 °C in the maximum and minimum temperatures respectively over the KKR. A decrease of ~ 280, 80 and 440 cm in the precipitation was also reported over the PPR, SR and GHR respectively during the same period 63 65 , unlike the DEMs of same resolution used in this study, explains the deviation of the two ice-thickness change estimates. The observed deviation may also be due to the different DEM resolutions and processing techniques used in the two studies. In the case of Vijay and Braun 13 , though the data and the period of observation are same, but the ice-thickness estimates are based on a small subset (2,308) of the large number of glaciers used in this study (12,243). The glacier thickness changes estimated in this study agree well with the previous studies showing stability or even mass gain in the Karakoram and negative changes in the rest of the Himalaya 8,9,13,15 . The observed ice thickness changes of − 0.11 ± 0.32 m a −1 over the KKR agree reasonably well with the estimates of + 0.12 ± 0.19 m a −1 (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) reported by Gardelle et al. 15 ; − 0.03 ± 0.04 m a −1 by Kääb et al. 65 ; and − 0.19 ± 0.22 m a −1 by Vijay and Braun 17 . The results are however relatively higher, though with significant intra basin variations, when compared with the average glacier thickness change estimate of − 0.21 ± 0.05 m a −1 for the Hindu Kush-Karakoram-Himalaya (HKKH) 65,66 . conclusion In this study, the glacier thickness changes over the UIB in the Jammu and Kashmir Himalaya were quantified using the SRTM-C and TanDEM-X Digital Elevation Models (DEMs) during the period from 2000 to 2012. The study concluded that the glaciers in the region have thinned at the rate of − 0.35 ± 0.33 m a −1 which amounts to the glacier stored water loss of 70.32 ± 66.69 Gt during 12 year observation period. The variability in the observed glacier thickness changes across the six mountain ranges in the region is explained well by the variability of topographic parameters across the ranges. However, further investigations aimed at understanding the range-wise glacier response to climatology would provide further insights into the differing glacier behaviour and response observed over the topographically complex mountainous UIB. The glacier thickness changes across different mountain ranges of the data scarce Jammu and Kashmir Himalaya presented in this study is vital for determining the sustainability of water resources in the UIB.

Data availability
The dataset will be available from the corresponding upon request.