Region-wide glacier mass balances over the Pamir-Karakoram-Himalaya during 1999–2011

The recent evolution of Pamir-Karakoram-Himalaya (PKH) glaciers, widely acknowledged as valuable high-altitude as well as mid-latitude climatic indicators, remains poorly known. To estimate the region-wide glacier mass balance for 9 study sites spread from the Pamir to the Hengduan Shan (eastern Himalaya), we compared the 2000 Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) to recent (2008-2011) DEMs derived from SPOT5 stereo imagery. During the last decade, the region-wide glacier mass balances were contrasted with moderate mass losses in the eastern and central Himalaya (-0.22 ± 0.12 m w.e./yr to -0.33 ± 0.14 m w.e./yr) and larger losses in the western Himalaya (-0.45 ± 0.13 m w.e./yr). Recently reported slight mass gain or balanced mass budget of glaciers in the central Karakoram is confirmed for a larger area (+0.10 ± 0.16 m w.e./yr) and also observed for glaciers in the western Pamir (+0.14 ± 0.13 m w.e./yr). Thus, the "Karakoram anomaly" should be renamed the "Pamir-Karakoram anomaly", at least for the last decade. The overall mass balance of PKH glaciers, -0.14 ± 0.08 m w.e./yr, is two to three times less negative than the global average for glaciers distinct from the Greenland and Antarctic ice sheets. Together with recent studies using ICESat and GRACE data, DEM differencing confirms a contrasted pattern of glacier mass change in the PKH during the first decade of the 21st century.


Introduction
The Pamir-Karakoram-Himalaya (PKH) mountain ranges are covered by more than 70 000 km 2 of glaciers . Assessing glacier evolution over such a large and remote region is challenging, but nevertheless required to characterize the impacts of climate change in the region (e.g. , to assess glacial contribution to water resources (e.g. Immerzeel et al., 2010;Kaser et al., 2010) and global sea level rise (e.g. Kääb et al., 2012;Gardner et al., 2013) and, ultimately, to reliably project glacier response to 21st century climate changes (e.g. Radic and Hock, 2011;Marzeion et al., 2012).
Length changes have been measured for about 100 glaciers in the PKH and revealed predominant retreat of glacier fronts since the mid-19th century, except in the Karakoram (Scherler et al., 2011a;Bhambri et al., 2012;. The majority of glaciers lost area over the past decades and in most cases, the rates of area loss have been increasing in recent years (e.g. . But changes in glacier length and area should be treated with care when used to evaluate the impact of climate change on glaciers, especially in the presence of surge-type and/or debris-covered glaciers, which are common in the PKH (Yde and Paasche, 2010).
Rather, mass balance is the most relevant variable to assess glacier responses to climate variability (Oerlemans, 2001;Vincent, 2002). Mass balance estimates, however, remain scarce in the PKH and may not adequately sample the wide range of climates and glacier responses in the region. Five Published by Copernicus Publications on behalf of the European Geosciences Union. different methods of mass balance measurements have been used in the PKH: i. Glaciological (or traditional) measurements are relatively short-termed, in general less than 10 yr, mainly because of difficult access to generally remote glaciers. In addition, measurement series are often biased towards small-to-medium-sized and debris-free glaciers for logistical reasons (e.g. Dobhal et al., 2008;Fujita and Nuimura, 2011;Azam et al., 2012;Yao et al., 2012). Vincent et al. (2013) discussed in detail the inadequate spatial and temporal sampling of those glaciological measurements in the western Himalaya. Gardner et al. (2013) have shown that, during [2003][2004][2005][2006][2007][2008][2009], the extrapolation of those few and shortterm glaciological mass balances to entire High Mountain Asia lead to an overestimation of the mass loss by about a factor of 3.
ii. The accumulation area ratio (AAR) method relies heavily on the availability of glaciological mass balances during various years in order to establish an empirical relationship between annual mass balance and AAR (e.g. Rabatel et al., 2005). The extrapolation to unmeasured glaciers is not straightforward. This method has been applied by Kulkarni et al. (2011) to 19 glaciers in north-west India.
iii. The hydrological method has only been used once in the PKH to our knowledge, providing 5 yr of mass balance for Siachen Glacier, Karakoram (Bhutiyani, 1999), and is very difficult to implement due to a lack of accurate precipitation and runoff measurements at high altitude over the PKH Immerzeel et al., 2012).
iv. The geodetic method, based on the comparison of topographic data (DEM or laser altimetry), can be applied at regional scales in remote areas using satellite data such as from ICESat, SPOT5 or SRTM. This method enables an increase in the number and diversity of glacier measurements (e.g. in terms of glacier size, elevation, aspect, percentage of debris cover) Bolch et al., 2011;Kääb et al., 2012;Nuimura et al., 2012;Pieczonka et al., 2013). An important, remaining limitation is the difficulty to convert volume changes to mass changes (e.g. Huss, 2013).
v. The gravimetric method uses the data from the Gravity Recovery and Climate Experiment (GRACE) mission launched in 2002. Gravity fields need to be corrected for the influence of the glacial isostatic adjustment (GIA) and the hydrological changes off-glacier (e.g. Zhang et al., 2013a) before the glacier mass change can be extracted. This method initially led to contrasted glacier mass budgets in High Mountain Asia (Matsuo and Heki, 2010;Jacob et al., 2012). The two most recent GRACE mass budgets for High Mountain Asia agree within their 2σ error bars (−14 ± 17 Gt yr −1 and −23 ± 24 Gt yr −1 in Gardner et al., 2013).
The ICESat-and GRACE-based methods are appropriate to assess the regional mass budget for large glacierized regions such as the PKH but also have their specific limitations. The ICESat sampling is restricted to satellite orbits that are separated by tens of kilometers at PKH latitudes and, thus, does not allow estimating the mass balance of individual glaciers. Glacier changes can only be investigated over the lifetime of ICESat, 2003ICESat, -2009. The gravimetric method has a coarse spatial resolution (typically ∼ 400 km) so that individual glaciers and small hydrological basins cannot be resolved.
The aim of our study is to present a homogeneous survey of regional glacier mass balances along the PKH between 1999 and 2011 using the geodetic method. This is achieved by differencing two digital elevation models (DEMs) over 9 study sites, selected to be representative of the PKH climatic and glaciological diversity. This allows for an estimation of the spatial pattern of glacier elevation changes along the 3000 km-long group of mountain ranges, as well as the influence of debris cover on thinning rates. Finally, we extrapolate these measurements to provide a new estimate of the mass budget of PKH glaciers, and their contribution to regional hydrology and sea level rise during the first decade of the 21st century.
Compared to Kääb et al. (2012), the length of the study period has been doubled and the study area has been significantly extended towards the east (Hengduan Shan, China) and the west (Pamir, Tajikistan). We also provide a nearly exhaustive coverage of glacier elevation changes for our 9 study sites, as opposed to the sampling by ICESat along satellite tracks only. However, contrary to the ICESat studies Gardner et al., 2013) our method is limited by the availability of SPOT5 DEMs on selected sites and does not provide insight into the seasonal and annual evolution of glacier mass balance.

Study area
Mass balances are investigated for 9 study sites spread along the PKH mountain ranges to capture the climatic and glaciological variability of the region. The location of each study site is displayed in Fig. 1, as well as the extent of the corresponding sub-regions used to extrapolate the results to the whole PKH range (see Sect. 3.4). PKH glacier meltwater contributes to five major rivers of Asia (Brahmaputra, Ganges, Indus, Tarim and Amu Darya), whose catchment areas are also presented in Fig. 1 PKH range could not be achieved because of funding restrictions (SPOT5 is a commercial satellite), the busy acquisition schedule of this satellite and cloud coverage. A further constrain is that the optical stereo images used to compute the DEMs should be acquired as close as possible to the end of the glacier ablation season. In practice, a time window of 2-3 months is allowed for the acquisitions, which is short given the 26-day orbital cycle of SPOT5. Thus, a maximum of 3 to 4 acquisition attempts are possible each year, and only if the satellite is not programmed for higher priority targets (e.g. crisis situations). Despite these difficulties, our SPOT5 DEMs sampled in total 22 200 km 2 of glaciers, about 30 % of the total glacierized area in the PKH . The eastern-most sites (Hengduan Shan, Bhutan, Everest and West Nepal) are strongly influenced by the Indian and south-east Asian summer monsoons. Ablation and accumulation of monsoon-type glaciers occur during the summer season (e.g. Fujita, 2008). Towards the other end of the PKH, in the north-west (Pamir, Hindu Kush and Karakoram), the climate is dominated by the westerlies and glaciers are of winter-accumulation type. The Spiti Lahaul site lies in a transition zone, influenced both by the monsoon and the westerlies. The topography also plays a strong role in the moisture transfer, as it dries out the southerly air flow and can prevent the air mass from travelling further north, which results in a northward decrease of precipitation (Bookhagen and Burbank, 2010).
Monsoon-type glaciers are expected to be more sensitive than winter-accumulation type glaciers to a change in the rain/snow limit driven by an increase in temperature (Fujita, 2008). Indeed, a summer warming would increase the alti-tude of the rain/snow limit and reduce snow accumulation on glaciers, as well as decrease their surface albedo (less fresh snow with a high albedo) and thus enhance ablation.
The Karakoram and the Pamir are known to host numerous surge-type glaciers (Barrand and Murray, 2006;Hewitt, 2007;Kotlyakov et al., 2008;Gardelle et al., 2012a), characterized by the alternation between short active phases involving rapid mass transfer from high to low elevations, and longer quiescent phases of low mass fluxes. Copland et al. (2011) reported a doubling in the number of glacier surges after 1990 in the Karakoram, with a total of 90 surgetype glaciers observed in the region since the late 1960s. In the Pamir, an inventory from 1991 revealed 215 glaciers with signs of repeated surging (looped moraines, fast advance of the front) and 51 additional actively surging glaciers (Kotlyakov et al., 2008).
Many PKH glaciers are heavily debris-covered in their ablation areas, because of steep rock walls that surround them, and intense avalanche activity (e.g. . The size of these debris zones is highly variable and debris thickness can range from a few centimeters (dust or sand) to several meters (Mihalcea et al., 2008;Zhang et al., 2013b). The mean debris cover is about 10 % of the total glacier area in the Karakoram and Himalaya (Bolch et al., 2012), and culminate at 36 % in the Central Himalaya (Scherler et al., 2011a).
In the eastern PKH, glacier termini are often connected to pro-glacial lakes storing meltwater behind frontal moraines or dead-ice dams, whereas in the western PKH, pro-glacial lakes are less numerous and glacier ablation areas are often scattered with supra-glacial lakes , resulting from successive coalescence of small melting ponds.

Digital elevation models
For each study site, we used the Shuttle Radar Topographic Mission (SRTM) version 4 DEM as the reference topography (Rabus et al., 2003). This data set, acquired in February 2000, comes along with a mask to identify the pixels which have been interpolated (both were downloaded from http://srtm.csi.cgiar.org/). On each site, the SRTM DEM is subtracted from a more recent DEM, built from a stereo pair acquired by the HRS sensor onboard the SPOT5 satellite (Korona et al., 2009), except for the Hindu Kush study site where the recent DEM has been created from a stereo pair acquired by the HRG sensor, also onboard SPOT5 . The date of each SPOT5 DEM acquisition differs depending on the study site (Table 1) but 7 out of 9 were acquired between October 2010 and November 2011. Each SPOT5-HRS DEM is also provided with a correlation mask and an ortho-image generated from the rear HRS image. SPOT5-HRS DEMs were produced by the French Mapping Agency (IGN) using two sets of correlation parameters defined during the SPIRIT (SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies) international polar year project (Korona et al., 2009). Thus, two versions of the DEM are delivered with their respective mask, one optimized (v1) for a gentle topography and the second one (v2) for a more rugged relief. Earlier work has shown that v2 has a slightly larger percentage of interpolated pixels but, for the non-interpolated pixels, smaller errors (Korona et al., 2009). This version (v2) is preferred here. Given that both DEMs are calculated from the same image pair, little difference is expected between v1 and v2. However, locally, we identified some large, bull eye, elevation differences between the two versions of the DEM, both off and on glaciers. Comparison with the SRTM DEM shows that both versions of the DEM are erroneous for those areas and thus we preferred to exclude from further analysis all pixels for which the absolute elevation difference between the two versions of the DEM is larger than 5 m. About 20 % of the valid DEM pixels are excluded through that procedure.
The SRTM DEM and mask, originally at a 3 arcsecond resolution (∼ 90 m), are resampled bilinearly to 40 m (UTM projection, WGS84 ellipsoid) to match the posting and projection of the SPOT5 DEM. Altitudes are defined above the EGM96 geoid for both DEMs.

Glacier mask
Where available, glacier outlines were downloaded from the GLIMS database (Raup et al., 2007;Supplement) and manually corrected if needed. Otherwise, the glacier mask was derived from Landsat-TM and Landsat-ETM+ images. Path/row and dates of Landsat acquisition are given in Table 1 for each study site. A threshold varying, between 0.5 and 0.7 depending on the study site, was applied to the Normalized Difference Snow Index (TM2-TM5)/(TM2+TM5) to detect clean ice and snow automatically, whereas debriscovered parts were digitized manually by visual interpretation (e.g. Racoviteanu et al., 2009). The total glacier area and debris-covered part can be found for each study site in Table 1.
Most inventories are derived from Landsat images with minimal snow cover from around the year 2000, prior to the 2003 failure of the Scan Line Corrector of the ETM+ sensor onboard Landsat 7 that resulted in image striping. Glacier areas are expected to have experienced minor changes since then, except for glacier fronts that retreated due to the rapid The Cryosphere, 7, 1263-1286, 2013 www.the-cryosphere.net/7/1263/2013/ growth of connected lakes and for glaciers that surged and thus advanced. For the few glaciers that advanced, their front position was updated manually based on the SPOT5 orthoimage.
For each study site, we favored our user-defined glacier inventories instead of the global Randolph Glacier Inventory (RGI v2.0, Arendt et al., 2012) because the latter is heterogeneous and, for some study sites, did not match our accuracy requirements for the adjustments and corrections of the DEMs.
We also estimated the altitude of the equilibrium line (ELA) as the snowline on the ∼ 2000 Landsat data with minimal snow cover to separate ablation and accumulation areas. The ELA is assumed to be identical to the altitude of the snowline at the end of the ablation season (e.g. Rabatel et al., 2005), after the end of the monsoons or before the first snowfalls. Therefore, we manually digitized snowlines for more than 30 glaciers for each study site and computed their mean elevations. Our ELAs are slightly different from those given by Scherler et al. (2011a) and Kääb et al. (2012) probably due to the sensitivity of the ELA to the choice of the image (Table 2). Ideally, the ELA should have been estimated from images acquired at the end of the ablation season during various years covering our whole study period. This would be a highly time-consuming task that was not accomplished here. However, we stress that, in our study, the ELA is only used to compare elevation changes on the ablation/accumulation areas (Sect. 4.1) and to estimate the sensitivity of our mass balances to the choice of the volume-tomass conversion factor (Sect. 3.5). Our preferred mass balance estimate uses a unique volume-to-mass conversion factor for the whole glacier and is, thus, independent from the ELA.

Adjustments and corrections of DEM biases
Different DEM processing steps are followed to extract unbiased glacier elevation changes for each study site (e.g. Nuth and Kääb, 2011;Gardelle et al., 2012b).

Planimetric adjustment of the DEMs
A horizontal shift between two DEMs results in an aspectdependent bias of elevation differences (Nuth and Kääb, 2011). Here, the horizontal shift is determined by minimizing the root mean square error of elevation differences (SPOT5-SRTM) off glaciers, where the terrain is assumed to be stable over the study period . The SRTM DEM is then resampled (cubic convolution) according to the shift.

Along/across track corrections (drift of the satellite orbit)
For DEMs derived from stereo imagery, the satellite acquisition geometry can induce a bias along and/or across the track direction . We estimate the azimuth of each SPOT5 ground track and use it to rotate the coordinate system accordingly (Nuth and Kääb, 2011). Then, we compute the elevation differences along and across the satellite track on stable areas and when necessary correct the bias using a 5th order polynomial fit.

Curvature correction
As suggested by Paul (2008) and Gardelle et al. (2012b), the difference in original spatial resolutions of the two DEMs can lead to biases related to altitude in mountainous areas. When the curvature, the first derivative of the slope, is high (sharp peaks or ridges), the altitude tends to be underestimated by the coarse DEM, unable to reproduce high-frequency www.the-cryosphere.net/7/1263/2013/ The Cryosphere, 7, 1263-1286, 2013 slope variations. Thus, this apparent "elevation bias" is corrected using the relation between elevation differences and the maximum terrain curvature estimated over stable areas off glaciers (Gardelle et al., 2012b). Note that the units in Fig. 1b in Gardelle et al. (2012b) should be 10 −2 m −1 , instead of m −1 (A. Gardner, personal communication, 2012). However, this does not impact the validity of the correction.

Radar penetration correction
The SRTM DEM, acquired in C-band, potentially underestimates glacier elevations, since at this wavelength (∼ 5.6 cm), the penetration of the radar signal into snow and ice can reach several meters (e.g. Rignot et al., 2001). Here, we estimate this penetration for each study site by differencing the SRTM C-band (5.7 GHz) and X-band (9.7 GHz) DEMs. The latter was acquired simultaneously with the SRTM X-band at higher resolution (30 m) but with a narrower swath, and thus has a coverage restricted to selected swaths only. Since the penetration of the X-band is expected to be low compared to the C-band (an hypothesis that still needs to be confirmed; Ulaby et al., 1986), we consider the elevation difference (SRTM C-band − SRTM X-band ) to be a first approximation of C-band radar penetration over glaciers (Gardelle et al., 2012b). The correction is calculated as a function of altitude and applied to each pixel separately. For the Hindu Kush, Spiti Lahaul and West Nepal sites, due to the lack of SRTM X-band data, we were unable to determine the value of the penetration. Thus, we applied the correction computed over the nearest study site, i.e. Karakoram for Spiti Lahaul and Hindu Kush, and Everest for West Nepal. The mean values of the SRTM C-band penetrations over glaciers are given for each study site in Table 3.

Seasonality correction
SPOT5 DEMs were acquired between late October and early January, depending on the study site. Since the SRTM DEM is from mid-February, we need to account for possible mass changes during this 1-to-5-month period in order to estimate glacier mass balance over an integer number of years (Gardelle et al., 2012a). For the eastern sites (from West Nepal to Hengduan Shan) where glaciers are of summeraccumulation types (Fujita, 2008), the correction was set to 0. This choice is confirmed by recent field observations that show no winter accumulation on Mera and Pokalde glaciers in Nepal . For the Karakoram study sites, we used the winter accumulation rate +0.13 m w.e. per month measured on Biafo Glacier (Wake, 1989). For the other western study sites (the Pamir, Hindu Kush and Spiti Lahaul sites), the correction is derived from the mean of winter mass balances of 35 glaciers, all in the Northern Hemisphere, +0.89 m w.e. yr −1 over 2000-2005, corresponding to +0.15 m w.e. per winter month (Ohmura, 2011). The latter value is slightly lower than the average winter mass balance Table 3. Average SRTM C-band penetration estimates (in m) in February 2000 over glaciers in this study and from .

Mean elevation changes and mass balance calculation
Before averaging elevation changes, we exclude all interpolated pixels in the SRTM or SPOT5 DEMs, as well as unexpected elevation changes exceeding ±150 m for study sites (Pamir and Karakoram) that include surge-type glaciers or ±80 m for other sites. Those thresholds were chosen after careful inspection of the elevation difference maps and histograms in order to identify the largest realistic glacier elevation changes. Glaciers that are truncated at the edges of the DEM are also excluded from mass balance calculations as they may bias the final results if, for example, only their accumulation or ablation area is sampled. Although they are truncated, Siachen (Karakoram East) and Fedtchenko (Pamir) glaciers were retained because of their wide coverage and because their accumulation and ablation areas were correctly sampled in the DEMs. Elevation changes are then analyzed for 100 m altitude bins. Within each altitude band, we average the elevation changes after excluding pixels where absolute elevation differences differ by more than three standard deviations from the mean (Berthier et al., 2004;Gardner et al., 2012). This is an efficient way to exclude outliers, based on the assumption that elevation changes should be similar at a given altitude within each study site. This assumption does not hold for regions containing actively surging glaciers. Thus, for the Pamir and Karakoram study sites, surge-type glaciers are treated separately, i.e. elevation changes are averaged by 100 m altitude bands for each surge-type glacier individually.
Where no elevation change is available for a pixel, we assign to it the value of the mean elevation change of the The Cryosphere, 7, 1263-1286, 2013 www.the-cryosphere.net/7/1263/2013/ altitude band it belongs to, in order to assess the mass balance over the whole glacier area. The conversion of elevation changes to mass balance requires knowledge of the density of the material that has been lost or gained and its evolution during the study period. Given the lack of repeat measurement of density profiles over the entire snow/firn/ice column in the PKH, we applied a constant density conversion factor of 850 kg m −3 for all our study sites (Sapiano et al., 1998;Huss, 2013).
The mass balance of each surge-type glacier is computed separately and added (area-weighted) to the mass balance of the rest of the glaciers to estimate the mass balance of the whole study site.
For study sites with a high concentration of growing pro-glacial lakes, such as Bhutan, Everest and West Nepal , our mass losses are slightly underestimated because they do not take into account the glacier ice that has been replaced by water during the expansion of the lake. Subaqueous ice losses are only estimated for Lhotse Shar/Imja Glacier (Everest area, Fujita et al., 2009) for the sake of comparison to previous studious.
The region-wide mass balance of PKH glaciers, as defined in Sect. 2, is computed by extrapolating the mass balance from each study site to a wider sub-region ( Fig. 1) that is assumed to experience the same glacier mass changes because of its similar climate Kääb et al., 2012). The total glacier area for each sub-region is computed from the RGI v2.0 .
The periods over which we calculate mass balance are slightly different from one site to another (Table 1) depending on the year of acquisition of the SPOT5 DEM (two in 2008, two in 2010 and five in 2011). In the following, we will assume that all mass balances are representative of the period 1999-2011, in order to estimate a region-wide mass budget of PKH glaciers, neglecting thus part of the inter-annual variability.

Accuracy assessment
The error E h i , to be expected for a pixel elevation change h i , is assumed to equal the standard deviation σ h of the mean elevation change h of the altitude band it belongs to. The value of σ h can range from ±4 m to ± 20 m depending on the study site and elevation. This is rather conservative as the value of σ h contains also the intrinsic natural variability of elevation changes within the altitude band (i.e. it contains both noise and real geophysical signal).
The error E h of the mean elevation change h in each altitude band is then calculated according to standard principles of error propagation: N eff represents the number of independent measurements in the altitude band. N eff will be lower than N tot , the total num-ber of elevation change measurements h i in the altitude band, since the latter are correlated spatially (Bretherton et al., 1999): where PS is the pixel size (here 40 m), d is the distance of spatial autocorrelation, determined using Moran's I autocorrelation index on elevation differences off glaciers. On the average over the 9 study sites, d = 492 ± 72 m. The error on the penetration correction was assumed to be systematic due to the not fully understood physics involved and equal to ± 1.5 m. This error was obtained by comparing the SRTM C-band penetration estimated using two independent methods (i) SRTM C-band -SRTM X-band and (ii) difference between ICESat and SRTM when ICESat elevation change trends during 2003-2008 are extrapolated to February 2000 . The maximum difference of the penetration inferred from the two methods is 1.0 m (Table 3). A 50 % additional error margin was added to account for the fact that this comparison could only be made on 2 of our study sites.
Given the slender observational support for the seasonality correction (Section 3.3), we assumed its uncertainty to be ±100 % on the western sites, i.e. ±0.15 m w.e. per winter month. The same absolute uncertainty (±0.15 m w.e. per winter month) was used for the eastern sites where no seasonality correction is applied but ablation and/or accumulation may occur in winter .
All errors are summed quadratically within each altitude band, and then summed for all altitude bands assuming, conservatively, that they are 100 % dependent. The error on the SRTM penetration correction clearly dominates our error budget for volume change.
During the conversion from volume to mass, we assume a ±60 kg m −3 error on the density conversion factor (Huss, 2013). This ±7 % error is similar to the one estimated by Bolch et al. (2011). In a sensitivity test, the volume to mass conversion is also performed using two alternative density scenarios : (i) a density of 900 kg m −3 is assumed everywhere and (ii) 600 kg m −3 in the accumulation area and 900 kg m −3 in the ablation area. These two scenarios take, among other things, into account that elevation changes could be due to changes in glacier dynamics or due to surface mass balance. We then calculate the maximum absolute difference between the mass balances calculated using our preferred scenario (850 ± 60 kg m −3 everywhere) and the mass balances calculated using the two others scenarios. For the eastern sites (West Nepal to Hengduan Shan), the maximum absolute difference is 0.03 m w.e. yr −1 . For the western sites (Pamir to Spiti Lahaul), it reaches 0.06 m w.e. yr −1 . Those uncertainties, due to the choice of a given density scenario, remain negligible compared to other sources of errors.
We also include an error due to extrapolation of our mass balance estimate from our study site (the extent of each www.the-cryosphere.net/7/1263/2013/ The Cryosphere, 7, 1263-1286, 2013 SPOT5 DEM) to the unmeasured glaciers within each subregion. To estimate this regional extrapolation error, we compare the ICESat elevation change trends (computed as in Kääb et al., 2012) within 3 • × 3 • cells centered at the study sites (Table 5) with the ICESat-derived trends for the entire sub-regions. The absolute difference between both trends is less than 0.02-0.03 m yr −1 for Bhutan, West Nepal, Spiti Lahaul and Karakoram, i.e. negligible. For Everest, the elevation trend for the entire sub-region is 0.12 m yr −1 less negative than for the 3 • × 3 • cells around the SPOT5 DEM center; for Hindu Kush the trend for the entire sub-region is 0.17 m yr −1 more negative than for the 3 • × 3 • cell due to the variability of elevation changes within this sub-region ). The regional extrapolation error is obtained as the area-weighted mean absolute difference between the mass balances derived from ICESat for those 3 • × 3 • cells and the whole sub-region and equals ±0.04 m w.e. yr −1 . Finally, for sub-regions and total PKH estimates of mass change, we include an error on the glacier area computed from the RGI v2.0 . The inventory error is specific to each study site and evaluated based on the comparison of our user-defined inventory on each study site with the RGI v2.0 (Table 1). We note the remarkable accuracy of the RGI for all but one of our study sites. The relative errors are generally of a few percents and up to 12 % for West Nepal. The differences between our and existing inventories are probably due to the difficulty of delimitating debris-covered glacier parts (e.g. Frey et al., 2012; and accumulation areas. The Hengduan Shan study site is an exception. There, the RGI is based on the Chinese Glacier Inventory (Shi et al., 2009;Arendt et al., 2012) and overestimates the ice-covered area by 88 % compared to our Landsat-based inventory.
Unless stated otherwise, all error bars from this study and previously published studies correspond to one standard error.

Overview of the elevation changes
The maps of glacier elevation changes are given in Figs. 2-10 for our nine study sites with, as inset, the distribution of the elevation differences off glaciers. The full maps of elevation differences off glaciers are shown in the Supplement (Figs. S2-S10). Those off glacier maps indicate a local random noise of about ±5-10 m. However, at length scales of a few kilometers, the spatial variations in elevation differences are small, typically less than 1 m.
For the four eastern study sites, from the Hengduan Shan to the West Nepal (Fig. 1), the elevation changes averaged over the accumulation areas are small (between 0.03 ± 0.14 m yr −1 to 0.13 ± 0.15 m yr −1 ) whereas clear surface lowering is observed for all four study sites in their ablation areas, ranging from −0.50 ± 0.14 m yr −1 (Bhutan) to −0.81 ± 0.13 m yr −1 (Hengduan Shan). Over the Bhutan, Everest and West Nepal study sites, glaciers in contact with a proglacial lake often experienced larger thinning rates.
Further west, the Spiti Lahaul is the only site where surface lowering is measured both in the accumulation area (−0.34 ± 0.16 m yr −1 ) and the ablation area (−0.70 ± 0.14 m yr −1 ). Clear surface lowering is also observed in the ablation areas of the Hindu Kush glaciers (−0.50 ± 0.18 m yr −1 ).
In the Karakoram (east and west, Figs. 8 and 9) and in the Pamir (Fig. 10), the pattern of elevation changes is highly heterogeneous due to the occurrence of numerous glacier surges or similar flow instabilities. The elevation changes in the accumulation areas of those three study sites are slightly positive (between +0.22 ± 0.14 m yr −1 and +0.30 ± 0.18 m yr −1 ). In the ablation areas, small surface lowering is observed for the two Karakoram study sites (average of the two sites: −0.33 ± 0.16 m yr −1 ), whereas no elevation change is found in Pamir (−0.02 ± 0.13 m yr −1 ).
Note that those mean rates of elevation changes are sensitive to the choice of the ELA, assumed to equal the snowline altitude in Landsat images from around year 2000.

Influence of a debris cover
To evaluate thinning rates over debris-covered and debrisfree ice, we perform a histogram adjustment in order to compare data sets with similar altitude distributions. This is needed because debris-covered parts tend to be concentrated at lower elevations. The adjustment consists in randomly excluding pixels over debris-free ice from each elevation band, to match a scaled version of the debris-covered pixel distribution . Elevation changes with altitude can then be compared for clean and debris-covered ice (Fig. 11).
The thinning is higher for debris-covered ice in the Everest area and on the lowermost parts of glaciers in the Pamir. But on average, in the Pamir, western Karakoram and Spiti Lahaul, the lowering of debris-free and debris-covered ice is similar. The thinning is stronger over clean ice in the Hengduan Shan, Bhutan, West Nepal, Hindu Kush and, although to a lesser extent, in the eastern Karakoram.

Mass balance over the nine study sites
Glacier mass changes have been heterogeneous over the PKH since 1999. Slight mass gains or balanced mass budgets are found in the north-west, in the Pamir (+0.14 ± 0.13 m w.e. yr −1 ) and for the two study sites in the central Karakoram (+0.09 ± 0.18 m w.e. yr −1 and +0.11 ± 0.14 m w.e. yr −1 ). Mass loss is moderate in the adjacent Hindu Kush with a mass balance of −0.12 ± 0.16 m w.e. yr −1 . Glaciers in all our eastern study sites (two in Nepal, Bhutan and Hengduan Shan) experienced homogeneous mass losses with mass balances ranging The Cryosphere, 7, 1263-1286, 2013 www.the-cryosphere.net/7/1263/2013/  from −0.33 ± 0.14 m w.e. yr −1 to −0.22 ± 0.12 m w.e. yr −1 . The most negative mass balance is measured in Spiti Lahaul (western Himalaya), at −0.45 ± 0.13 m w.e. yr −1 . (Table 4). However, within a study site, mass balances can be highly variable from one glacier to another. In the Everest area, two large glaciers (71 km 2 each) experienced distinctly different mass balance with −0.16 ± 0.16 m w.e. yr −1 for the southflowing Ngozumpa Glacier and −0.59 ± 0.16 m w.e. yr −1 for the north-flowing Rongbuk Glacier (Fig. 4). In Bhutan (Fig. 3), mass loss is higher south (−0.25 ± 0.13 m w.e. yr −1 ) than north (−0.14 ± 0.12 m w.e. yr −1 ) of the main Himalayan ridge, though not at a statistically significant level.

Glacier surges
In the Pamir and the Karakoram, the spatial pattern of elevation changes of many glaciers is typical of surge events or flow instabilities. They can be identified because of their very high thinning and thickening rates that can both reach up to 16 m yr −1 (Figs. 8-10).

SRTM penetration
Alternative estimates of SRTM C-band penetration for the PKH region were obtained by Kääb et al.  glaciers are of summer-accumulation type. Further studies would be necessary to investigate the possible reasons behind the relatively low SRTM C-band penetration for the Pamir (1.8 m) that do not follow the east/west gradient mentioned above. Both our and Kääb et al.'s (2012) estimates agree within their error bars and seem thus to provide robust first-order estimates for the SRTM C-band radar penetration in the region. However, it remains to be verified if a penetration depth computed at the regional scale is appropriate for a single small glacier whose altitude distribution is different from the one of whole study site (e.g. Abramov Glacier in the far north of our Pamir study site).

Thinning over debris-covered ice
PKH glaciers are heavily debris-covered in their ablation areas, because of steep rock walls that surround them, and intense avalanche activity. Twelve percent (12 %) of the glacier area in our study sites are covered with debris (up to 22 % in the Everest area, Table 1), which is in agreement with previous estimates for the Himalaya and the Karakoram, 13 % in Kääb et al. (2012) and 10 % in . The debris layer is expected to modify the ablation of the under-lying ice. It will increase ablation if its thickness is just a few centimeters (dominance of the albedo decrease) or decrease ablation it if the layer is thick enough to protect the ice from solar radiation (Mattson et al., 1993;Mihalcea et al., 2006;Benn et al., 2012;Lejeune et al., 2013).
However, recent studies have reported similar overall thinning rates over debris-covered and debris-free ice in the PKH Gardelle et al., 2012a), or even higher lowering rates over debris-covered parts in the Everest area (Nuimura et al., 2012). Our findings are consistent with Nuimura et al. (2012) for the Everest study site, with a stronger thinning (rate of elevation change of −0.97 m yr −1 ) in debris-covered areas than over clean ice (rate of elevation change of −0.57 m yr −1 ) . In the Pamir, Karakoram and Spiti Lahaul study sites, these rates are similar, while in the Hindu Kush, West Nepal, Bhutan and Hengduan Shan study sites, thinning is greater over clean ice, in agreement with the commonly assumed protective effect of debris (Fig. 4). Those differences between our 9 study sites suggest a complex relationship between debris cover and glacier wastage.
As local glacier thickness changes are a result of both mass balance processes and a dynamic component, these unexpected high thinning rates over debris-covered parts The Cryosphere, 7, 1263-1286, 2013 www.the-cryosphere.net/7/1263/2013/ Fig. 8. Same as Fig. 2 but for the Karakoram east study site. Surge-type glaciers in an active surge phase (resp. quiescent phase) are identified with a solid triangle (resp. solid circle).
are potentially the result of glacier-wide processes. Over a glacier tongue, the ablation can be enhanced by the presence of steep exposed ice cliffs or supraglacial lakes and ponds (Sakai et al., 2000(Sakai et al., , 2002. In addition, it is also possible that the glacier tongues are mostly covered by thin debris layers, such that the albedo effect dominates the insulating effect, resulting in an increased tongue-wide ablation. The prevalence of thin debris was recently suggested for debris-covered glaciers around the Mount Gongga in the south-eastern Tibetan Plateau (Zhang et al., 2013b). Finally, ice-flow rates could be very low at debris-covered parts, as shown by Quincey et al. (2007) in the Everest area, Kääb (2005) for Bhutan, and Scherler et al. (2011b) in the Hindu Kush-Karakoram-Himalaya. Thus, all three factors are likely to increase the overall thinning under debris, despite the undoubted protective effect of a continuous thick debris cover at a local scale. Future investigations combining surface velocity fields and maps of elevation changes could allow evaluating more closely the relative role of ablation and ice fluxes in thinning rates over PKH glacier tongues (e.g. Berthier and Vincent, 2012;Nuth et al., 2012). Measurements of the spatial distribution of debris thickness over the PKH glacier tongues are also needed. Note that in the case of debris-covered tongues, the elevation change measurement represents the glacier thickness change but also the possible debris thickness evolution. The latter remains poorly known in the PKH or in any other mountain range. But based on the debris discharges given by Bishop et al. (1995) for Batura Glacier (Pakistan), 48 to 90 × 10 3 m 3 yr −1 , the thickening of the debris can be estimated between 2.7 and 5.5 mm yr −1 , which is negligible given the magnitude of the glacier elevation changes.

Comparison to previous mass balance estimates
Throughout the PKH, there is only a single peer-reviewed glaciological record (on Chhota Shigri Glacier) covering most of our study period . Our geodetic mass balance estimate for the Spiti Lahaul study site has been recently compared to those glaciological measurements on Chhota Shigri Glacier . Mass balance records reported in Yao et al. (2012) for the Tibetan Plateau and the surroundings mountain ranges only start during glaciological year 2005/2006. This is why, here, we do not compare our mass balance estimates to glaciological records and limit our comparison to others regional estimates obtained using the geodetic or gravimetric methods. We stress that the comparison of our mass balance measurements to published estimates is complicated by the fact that, generally, they do not cover the same time period. Little is known about the inter-annual variability in the PKH during the 21st century due to the limited number of long glaciological time series. The 9 yr mass balance record on Chhota Shigri Glacier reveals a relatively high inter-annual variability of the annual mass balance (standard deviation: ± 0.74 m w.e. yr −1 during , which is similar to the one obtained between 1968 and 1998 on Abramov Glacier (standard deviation: 0.69 m w.e. yr −1 , WGMS, 2012). Thus, inter-annual variability can explain part of the difference between two cumulative mass balance estimates over 5-10 yr even if they are offset by only one or two years.
Based on ICESat repeat track measurements and the SRTM DEM, Kääb et al. (2012)  For that same area (i.e., excluding the Pamir and the Hengduan Shan study sites), we found a mass balance of −0.15 ± 0.07 m w.e. yr −1 . Thus, our result agrees with Kääb et al. (2012) within the error bars, although our study period is longer (1999-2011) and our measurement principle and extrapolation approach are different. In addition, we have calculated updated mass balances from Kääb et al. (2012), i.e. averages over 3 × 3 degree cells centered over our study sites, and find that they are consistent, within error bars, with our mass balance estimates (Table 5). Similar results are found when our mass balances are compared to a spatially extended ICESat analysis for (Gardner et al., 2013. Mass losses measured with ICESat Gardner et al., 2013) tend to be slightly larger than ours. This could be due to the difference in time periods. Also, our smaller mass losses could be partly explained if we systematically underestimated the poorly constrained pene- In the present study, we found an average mass balance of −0.41 ± 0.21 m yr −1 over the ten glaciers mentioned above (Table 5, Table 4. Glacier mass budget of the nine sub-regions, to which the results of the study sites (mass balances) have been extrapolated. For each of them, the total glacierized area as well as the percentage of the glacier area over which the mass balance was computed (i.e. study sites) are given. The total glacier area is derived from the RGI v2.0 . The error for the mass balance in each sub-region is the root of sum of squares (RSS) of the mass balance error of each study site and the regional extrapolation error. The error for the mass budget in each sub-region includes additionally the error from the total glacier area in the RGI v2.0. +0.14 ± 0.14 +1.3 ± 1.3

Sub-region
Total/Area-72251 30 −0.14 ± 0.08 −10.1 ± 5.5 weighted average  11. Elevation changes (SPOT5-SRTM) in the ablation area for clean ice (blue circles) and debris-covered ice (black circles) for each study site. For the Karakoram (east and west) and the Pamir study sites, only non-surge-type glaciers are considered. Standard error of the elevation differences are typically less than ±2 m with each 100 m elevation band (not shown for the sake of clarity). Note, elevation changes over separate sections of a glacier cannot be treated as mass changes due to the disregard of glacier dynamics. estimate over this short time period. The differences in survey periods and in the glacier outlines may also explain part of these discrepancies. In particular, the delimitation of the accumulation area is notoriously difficult and Bolch et al. (2011) did not include some of the steep high altitude slopes, which we included due to their potential avalanche contribution. For the 10 glaciers mentioned above, our mass balances would become 0.12 m w.e. yr −1 more negative if we used the exact same outlines as Bolch et al. In Bhutan, our results indicate a stronger mass loss for southern glaciers than for northern ones. Interestingly, in this area, Kääb (2005) measured high speeds (up to 200 m yr −1 ) for glaciers north of the Himalayan main ridge, and rather low velocities over southern glacier tongues, by using repeat ASTER data. This is consistent with the more negative mass balance that we report for the southern glaciers in Bhutan, i.e. for the slow flowing, presumably downwasting, glaciers. Berthier et al. (2007) reported a mass balance between −0.69 and −0.85 m w.e. yr −1 for an ice-covered area of 915 km 2 around Chhota Shigri Glacier between 1999 (SRTM DEM) and 2004 (SPOT5-HRG DEM). For the same area, we found a mass balance of −0.45 ± 0.13 m w.e. yr −1 over 1999-2011. The difference in the study period may explain part of the differences. Also, the methodology by Berthier et al. (2007) did not explicitly take into account the SRTM radar penetration over glaciers, and also corrected an apparent elevation-dependent bias according to glacier elevations, which is now known as a resolution issue and is best corrected according to the terrain maximum curvature, (Gardelle et al., 2012b). More details and discussions about these differences can be found in Vincent et al. (2013).
Importantly, the results of the eastern Karakoram study site confirm the balanced or slightly positive mass budget reported by Gardelle et al. (2012a)  sites show similar mass balances over slightly different periods: +0.11 ± 0.14 m yr −1 over 1999-2010 in the east and +0.09 ± 0.18 m yr −1 over 1999-2008 in the west. This confirms a "Karakoram anomaly" that was first suggested based on field observations (Hewitt, 2005). We report a mass balance of +0.10 ± 0.13 m yr −1 for Baltoro Glacier (see location in Fig. 8, size = 304 km 2 ) between 1999 and 2010, which is consistent with its gradual speed-up reported since 2003 (Quincey et al., 2009). However, given the completely different methods used in our and the latter study, we cannot conclude that this demonstrates a real shift from negative to positive mass balance. The mass budget of Siachen Glacier (1145 km 2 , sometimes referred to as the highest battlefield in the world) is also balanced or slightly positive (+0.14 ± 0.14 m w.e. yr −1 ). Thus, the alarming statement by Rasul (2012) that this glacier retreated by 5.9 km and lost 17 % of its mass during 1989-2009 is very likely erroneous. Again, within error bars, our regional mass balance in the Karakoram agrees with the mass balance of −0.03 ± 0.04 m yr −1 found by Kääb et al. (2012Kääb et al. ( ) over 2003Kääb et al. ( -2008 and with the mass balance of −0.10 ± 0.18 m w.e. yr −1 (2-sigma error level) obtained by Gardner et al. (2013) for a region including both the Karakoram and the Hindu Kush.
In the northernmost part of the Pamir, in the Alay mountain range, the mass balance of Abramov Glacier has been measured in the field for 30 yr between 1968and 1998(WGMS, 2012, with a mean value of −0.46 m yr −1 . Here, we report a mass balance of −0.03 ± 0.14 m yr −1 for this glacier over 1999-2011. Our near zero mass budget estimate for this glacier disagrees with negative mass balances reconstructed for the same period with a model run using bias corrected NCEP/NCAR re-analysis data and calibrated with field data (Barundun et al., 2013). An underestimate of our SRTM C-Band penetration in the Pamir as a whole and in particular for the small Abramov Glacier may contribute to these differences.
For the whole Pamir study site, our mass budget is balanced or slightly positive (+0.14 ± 0.13 m w.e. yr −1 ). In the eastern Pamir, the mass balance of Muztag Ata Glacier was +0.25 m w.e. yr −1 between 2005 and 2010 (Yao et al., 2012). Although Muztag Ata Glacier is located outside of our Pamir study site and his glaciological records only cover the second half of our study period, this measurement is consistent with our remotely sensed mass balance. Thus the "Karakoram anomaly" should probably be renamed the "Pamir-Karakoram anomaly".
This slightly positive or balanced mass budget measured in the Pamir seems to contradict previous findings by Heid and Kääb (2012). They reported rapidly decreasing glacier speeds in this region between two snapshots of annual velocities in 2000/2001 and 2009/2010, an observation that would fit with negative mass balances (Span and Kuhn, 2003;Azam et al., 2012). We note, though, that the relationship between glacier mass balance and ice fluxes is not straightforward, and might in particular involve a time delay between a change in mass balance and the related dynamic response (Heid and Kääb, 2012). Thus, the decrease in speed could well be due to negative mass balances before or partially before 1999-2011. We acknowledge that this discrepancy needs to be investigated further, in particular by examining continuous changes in annual glacier speed through the first decade of the 21st Century. Indeed, Heid and Kääb (2012) measured differences between two annual periods which may not represent mean decadal velocity changes. Jacob et al. (2012), using satellite gravimetry (GRACE), measured a mass budget of −5 ± 6 Gt yr −1 over 2003-2010 for the "Karakoram-Himalaya", their region 8c, which corresponds to our study area excluding the Pamir, Karakoram and Hindu Kush sub-regions but including some glaciers north of the Himalayan ridge already on the Tibetan Plateau. For this subset of our PKH study area, we measured a mass budget of −12.6 ± 4.2 Gt yr −1 . The two estimates overlap within their error bars, especially if our error is doubled to match the two standard error level used by Jacob et al. (2012). The estimate of Jacob et al. (2012) over our complete study region (PKH), i.e. −6 ± 8 Gt yr −1 (sum of their regions 8b and 8c) also includes mass changes for the Kunlun Shan sub-region for which we did not measure glacier mass balances. Therefore the comparison with our PKH value (−10.1 ± 5.5 Gt yr −1 ) is not straightforward but suggests a reasonable agreement, especially if we consider the slight mass gain (1.5 ± 1.7 Gt yr −1 ) measured with ICESat in the West Kunlun (Gardner et al., 2013). , 7, 1263-1286, 2013 www.the-cryosphere.net/7/1263/2013/ Table 6. Annual average of the glacier seasonal and decadal mass loss contributions to Indus, Ganges and Brahmaputra discharges. Basin area and seasonal contributions are given by Kaser et al. (2010), glacier area in each basin is derived from the RGI v2.0 . Numbers in parentheses indicate the percentage of glacierized area within the river basin. Seasonal contributions were modeled by Kaser et al. (2010)

Reasons for the heterogeneous pattern of mass balance
The PKH glacier mass balance is slightly negative (−0.14 ± 0.08 m w.e. yr −1 ), but changes are not homogeneous from one sub-region to another and seem to coincide with the climatic settings of the mountain range. For the eastern sites (from the Hengduan Shan to West Nepal), which are under the influence of the Indian and south-east Asian summer monsoons, the mass balances are moderately negative (∼ −0.26 m w.e. yr −1 ). In the northwest of the PKH, at the Pamir and Karakoram sites, where glacier climate is characterized by the westerlies in winter, mass balances are balanced or slightly positive (∼ +0.10 m w.e. yr −1 ). In between these two influences, the glaciers of the Spiti Lahaul site, in a monsoon-arid transition zone, experienced the strongest mass loss (−0.45 ± 0.13 m w.e. yr −1 ). This heterogeneous pattern of mass balances seems to be related to trends in precipitation throughout the PKH. Archer and Fowler (2004) reported an increase in winter precipitation over the Upper Indus Basin since 1961, a trend confirmed by Yao et al. (2012) over the Pamir and the Karakoram, based on the analysis of the Global Precipitation Climatology Project data set (GPCP; Adler, et al., 2003). This precipitation increase is a source of greater accumulation for glaciers and thus a possible explanation for their slightly positive mass budgets. In addition, in the Upper Indus Basin, the mean summer temperature has been decreasing since 1961 (Fowler and Archer, 2006), which is consistent with the findings of Shekhar et al. (2010) in the Karakoram, who reported a decrease in both maximum and minimum temperature since 1988. These increases in winter precipitation and summer temperature likely translate in greater accumulation and reduced ablation, changes which are consistent with the balanced or slightly positive glacier mass budget.
On the other hand, in the east of the PKH, the Indian summer monsoon has been weakening since the 1950s (Bollasina et al., 2011), which could have reduced the amount of snowfall over the eastern study sites, and thus resulted in negative mass balances. In addition, maximum and minimum temperatures increased since 1988 in the western Himalaya (Shekhar et al., 2010), and maximum temperatures increased in the central Himalaya (Nepal) since the mid-1970s (Shrestha et al., 1999). Thus both precipitation and temperature trends in the Himalaya likely contributed to the negative mass balances measured over the corresponding study sites (from Spiti Lahaul to Hengduan Shan, Fig. 1).
However, gridded data sets, such as GPCP, are not necessarily representative of the climate at high altitude. Indeed, Hewitt (2005) reported a 5-to-10-fold increase in precipitation between the glacier front and the accumulation area in the Karakoram that is not captured by reanalysis data, as recently confirmed by Immerzeel et al. (2012). Thus, direct measurements of accumulation on High Mountain Asia glaciers (e.g. Azam et al., 2013) and high-altitude weather stations are eagerly needed to validate the large scale gridded data set (e.g. reanalysis), test their ability to describe the specific climate near to PKH glaciers and assess the relative role played by temperature and precipitation changes.

Contribution to water resources
Our glacier mass balance can be averaged over large river basins to estimate the mean contribution to river runoff due to a decade of glacier mass loss. We assume thereby only direct river runoff, without loss terms. For the three rivers for which we cover the entire glacierized area of their catchments within our sub-regions (Fig. 1), these contributions to the river discharge are equal to 103 m 3 s −1 (Indus), 103 m 3 s −1 (Ganges) and 147 m 3 s −1 (Brahmaputra) for the period 1999-2011 (Table 6).
Rivers of PKH being key water resources, measures of their discharges are highly sensitive and difficult to access for political reasons. Papa et al. (2012) built recent time series of discharges for Ganges and Brahmaputra by using altimetry measurements. The locations of the corresponding discharge estimates in Bangladesh are given in Fig. 1. We can therefore www.the-cryosphere.net/7/1263/2013/ The Cryosphere, 7, 1263-1286, 2013 evaluate the relative contribution of glacier decadal mass loss (Table 6) to annual river run-off at 0.9 % for the Ganges at Hardinge (basin area of ∼ 1 020 000 km 2 ) and 0.7 % for Brahmaputra at Bahadurabad (basin area of ∼ 530 000 km 2 ) between 1999 and 2011. Given the discharge measured at the Guddu station by the Surface Water Hydrology Project of the Water and Power Development Authority (Pakistan, see location on Fig. 1) between 1999 and 2003, we can estimate the contribution of glacier mass loss at 4.8 % for the Indus River at this station. The seasonal contributions of glaciers to river run-off (i.e. the part of the annual precipitation that experiences seasonally delayed release from the glaciers) directly results from seasonal variations of meteorological quantities, in contrast to the glacier mass loss contribution which results from the long-term climatic change. Both runoff components directly affect the amount of water available in a river at a given location and at a given time, so that their comparison could be of interest. Even if the seasonal timing of the glacier mass loss contribution is unclear, it could in general peak at a similar time of the year than the seasonal glacier contribution so that both components rather enlarge each other than cancel. Seasonal glacier contributions have been modeled by Kaser et al. (2010). The latter study assumed glaciers to be in equilibrium with the present climate and thus prescribed glacier mass balances equal to 0. Thus, their seasonally delayed glacier contributions do not include the part due to decadal glacier mass loss and is complementary to our estimates (Table 6). For the eastern basins under the influence of the Indian summer monsoon (Ganges and Brahmaputra), the contribution of glacier mass loss is higher than the seasonal contribution. In other words, for the first decade of the 21st century, the mass loss of glaciers is on average more important for water availability in those two rivers than their seasonal water storage effect. The situation is different for the Indus Basin where the glacier melt season is also the dry season, which results in a relatively higher seasonally delayed glacier contribution to the river discharge. There, the seasonally delayed water release from the glaciers and the decadal glacier mass loss contribute on average equally to the river discharge. Both contributions are strongly (and equally) dependent on the location of the discharge measurement and will be significantly larger (in relative term) in more heavily glacierized and smaller mountain catchments located in the upper part of these major river basins (Bookhagen and Burbank, 2010).

Conclusion
In this study, we assessed the spatial pattern of glacier mass balance over the PKH, by comparing the SRTM and SPOT5 DEMs over 9 well-distributed study sites, between 1999 and 2011. We found slightly positive or balanced mass budgets in the western part (Pamir, Karakoram), moderate mass loss in the east (Nepal, Bhutan, Hengduan Shan), and the most negative mass balance in the monsoon-arid transition zone (Spiti Lahaul in the western Himalaya). By extrapolating these values to climatically homogeneous sub-regions, we estimate the PKH overall mass balance to have been −0.14 ± 0.08 m w.e. yr −1 between 1999 and 2011, which corresponds to a contribution to sea level rise of 0.028 ± 0.015 mm yr −1 . This is about 4 % of the global glacier mass loss estimated by Gardner et al. (2013), though over a slightly shorter time period (2003)(2004)(2005)(2006)(2007)(2008)(2009). The largest source of uncertainties in those geodetic mass balance estimates is the poorly constrained penetration of the C-Band SRTM radar signal into snow and ice.
Our technique of DEM differencing allows mapping in detail the spatial pattern of glacier elevation changes. Those maps reveal many surge-type behaviors both in the Pamir and the Karakoram. It also appears that the glacier-wide mass balance does not seem to be considerably affected by the mass transfer from surges over the ∼ 10 yr time period studied. Similar lowering rates over debris-covered and clean ice are observed for four study sites, despite the commonly assumed insulating effect of debris covers. This suggests, for the scale of entire glacier tongues, a spatial mixture of reduced ablation under intact thick debris covers and enhanced ablation by supraglacial lakes or ice cliffs or due to thin debris layer, and very low glacier velocities in ablation areas that reduce the ice advection downstream.
The regionally heterogeneous pattern of ice wastage is in broad agreement with contrasted trends in large-scale precipitation and temperature. However, glaciological (e.g. accumulation) and meteorological records are needed at high altitude in the vicinity of glaciers to evaluate more closely the local influence of atmospheric variables on glacier mass balance and fully understand the reasons behind those contrasted mass budgets in the PKH.