Brief Communication: Contending Estimates of 2003–2008 Glacier Mass Balance over the Pamir–karakoram–himalaya

We present glacier thickness changes over the entire Pamir–Karakoram–Himalaya arc based on ICESat satellite altimetry data for 2003–2008. We highlight the importance of C-band penetration for studies based on the SRTM elevation model. This penetration seems to be of potentially larger magnitude and variability than previously assumed. The most negative rate of region-wide glacier elevation change (< −1 m yr −1) is observed for the eastern Nyainqêntanglha Shan. Conversely, glaciers of the western Kunlun Shan are slightly gaining volume, and Pamir and Karakoram seem to be on the western edge of this mass-gain anomaly rather than its centre. For the Ganges, Indus and Brahmaputra basins, the glacier mass change reaches −24 ± 2 Gt yr −1 , about 10 % of the current glacier contribution to sea-level rise. For selected catchments, we estimate glacier imbalance contributions to river runoff from a few percent to greater than 10 %.


Introduction and methods
Region-wide measurements of glacier volume or mass change are limited for the Pamir-Karakoram-Himalaya region, leaving room for speculation about the glacier response to climate change and its hydrological significance.Glacier mass change in high mountain Asia (or some part of it) have been obtained by (i) extrapolating the few existing in situ mass balance series (Cogley, 2011;Bolch et al., 2012;Yao et al., 2012), (ii) space gravimetry (Jacob et al., 2012;Gardner et al., 2013), (iii) laser altimetry (Kääb et al., 2012;Gardner et al., 2013;Neckel et al., 2014) and (iv) the differencing of digital elevation models (Gardelle et al., 2013).Between these studies that narrow down the range of uncertainties for core parts of this remote mountain region, significant inconsistencies remain.
The aims of this study are (i) to provide a new consistent regional-scale data set from the ICESat autumn laser campaigns (2003)(2004)(2005)(2006)(2007)(2008) by extending Kääb et al. (2012) to completely cover the study region by Gardelle et al. (2013) and several major river basins, (ii) to compare the results to other previous estimates of the Pamir-Karakoram-Himalaya glacier volume change and (iii) to roughly evaluate the contribution of glacier mass change to river run-off.
We follow the methods explained in Kääb et al. (2012) with a considerable extension towards the eastern Nyainqêntanglha Shan, the Pamir and part of the Tibetan Plateau (Fig. 1).In short, ICESat footprints are intersected with the February 2000 SRTM DEM and overlaid on the most snow-free multispectral Landsat images over ∼ 2000-2013 to manually classify footprints into three classes: glaciers, non-glaciers and water.Glacier elevation difference trends are then estimated regionally and at a 1 • × 1 • geographic grid by fitting a robust linear temporal trend to the time series of elevation differences between the SRTM DEM and individual ICESat footprint elevations.Trends are derived from autumn ICESat campaigns only (2009 ICESat winter campaigns excluded), because combined autumn and winter trends are sensitive to temporal variations in accumulation amount and timing, potentially introducing bias (see Supplement of Kääb et al., 2012).We confirm that our trends are not due to sampling bias of ICESat elevations by comparing ICESat elevation histograms with glacier hypsometry.The resulting elevation difference trends for all our zones are given in Table 1.

Thickening in the Karakoram and western Kunlun Shan
A first striking feature in the regional map of elevation difference trends (Fig. 1) is glacier thickness gain in the west-ern Kunlun Shan (∼ +0.1 m yr −1 ), agreeing with in situ mass balance and length change measurements (Yao et al., 2012).There is a southwest-to-northeast gradient from considerably negative glacier mass balances in Hindu Kush and Spiti-Lahaul to positive values in the Pamir-Karakoramwestern-Kunlun-Shan region (Fig. 1).This suggests the socalled Karakoram glacier mass-balance anomaly (Hewitt, 2011;Gardelle et al., 2012a), or Pamir-Karakoram anomaly (Gardelle et al., 2013), is rather the edge or southwest limit of an anomaly centred more to the northeast over the western Kunlun Shan, or Tarim Basin.The anomaly seems thus indeed the result of a larger-scale meteorological or climatic feature, and peculiarities of the Karakoram topography or glaciers (e.g.surge type, hypsometry, avalanche contribution; Hewitt, 2011) do not necessarily play a decisive role.Combined, the results by Gardner et al. (2013), Neckel et al. (2014) and the glacier elevation change pattern of Fig. 1 suggest the centre of the anomaly could be located over the Tibetan Plateau.Direct precipitation measurements in this region are scarce and thus trends are uncertain.Satellite-retrieved precipitation and gauge data (Global Precipitation Climatology Project) suggest an increase of precipitation over the study region north of Karakoram and east of Pamir (Yao et al., 2012).Chinese measurements show increased precipitation over the Tibetan Plateau (C.-Y. Xu, personal communication, 2014), andTao et al. (2014) suggest wetter conditions over the Tarim Basin since the mid 1980s.A number of abnormally wet years occurred during the early 21st century over the Tarim Basin and the Tibetan Plateau (Becker et al., 2013), in particular for the hydrological years 2003/2004 and 2005/2006.A recent climate modelling study proposes that stable or increasing snowfalls characterise the Karakoram anomaly on a background of increasing air temperatures (Kapnick et al., 2014).Despite the available studies and data, further research seems necessary to consolidate the precipitation and temperature trends and the reason behind the slight glacier volume gains.

Massive thinning in eastern Nyainqêntanglha Shan and Spiti-Lahaul
The other striking feature in Fig. 1 is the massive glacier thickness loss in the eastern Nyainqêntanglha Shan (between −1 and −1.7 m yr −1 ), also consistent with the large negative mass balances and frontal retreats in this zone (Yao et al., 2012).The glaciers of eastern Nyainqêntanglha Shan have the smallest total elevational range in our study region, indicating a large sensitivity to fluctuations in the equilibrium line altitude (Pelto, 2010;Loibl et al., 2014).The few available in situ mass balance measurements in the area suggest that the equilibrium line was over the vertical limits of the monitored glaciers in the late 2000s, and precipitation in this zone shows the strongest long-term decrease over the entire Pamir-Karakoram-Himalaya region (Yao et al., 2012;Becker et al., 2013).A similar pattern of glacier shrinkage, though less distinct, is found at the western end of the Great Himalaya Range within our Spiti-Lahaul zone and forms the cluster of second-largest thickness loss rates in this study (−0.5 to −0.7 m yr −1 ).Also here, Landsat data indicate that firn lines have risen in several years towards high glacier el-evations, resulting in very small accumulation areas or even their complete loss.
The 2003-2008 glacier thickness changes in the other study zones are all similar, on the order of ∼ −0.4 to −0.5 m yr −1 (Table 1), with more negative values in the Bhutan zone at the transition between the eastern Nyaiqêntanglha and Everest zones.We note that glaciers dominated by the summer monsoon (i.e.east of the Spiti-Lahaul) all show thickness losses (summer-accumulation type glaciers; Fujita, 2008; Kapnick et al., 2014;Maussion et al., 2014).Eastern Nyaiqêntanglha Shan, the zone with strongest glacier thickness loss, receives most accumulation during March-May (spring-accumulation type; Maussion et al., 2014).The glaciers with considerable winter accumulation under influence of the westerlies show a more mixed picture, with stable or growing thicknesses in the Karakoram and western Kunlun Shan but thickness losses for instance in the Hindu Kush.

Comparison to previous thickness change studies
The following comparison to other studies uses average glacier thickness changes rather than total mass changes in order to minimise effects from different delineations of study zones, glacier cover areas and density assumptions.From Hindu Kush and Karakoram in the west to Nepal in the east, results of all studies agree within their errors (Table 1).Results are most sensitive to zone delineation in the Hindu Kush, reflecting the strong spatial variability of glacier thickness change rates in this area (Fig. 1) and presumably also locally heterogeneous glacier behaviour (Sarikaya et al., 2012; see also below for Pamir).
Significant differences among the results of all studies are found over eastern Nyainqêntanglha Shan.Our results and those from Neckel et al. (2014) agree within the errors but not with Gardner et al. (2013), although all three studies are based on ICESat.While our study and Neckel et al. (2014) use ICESat footprint classifications from contemporary satellite images, Gardner et al. (2013) use Randolph Glacier Inventory outlines (RGI version 2.0; Pfeffer et al., 2014), which contain considerable errors of commission and omission in this zone (see Table 1 in Gardelle et al., 2013).Repeating our analysis with footprint classifications based on the Randolph Glacier Inventory results in less negative elevation difference trends on glaciers (∼ 20 % less negative) due to inclusion of non-glacier footprints.Contrarily, the elevation difference trends on land (−0.10 ± 0.06 m yr −1 when using our own footprint classification) become more negative if ICESat footprints are classified using RGI due to inclusion of glacier footprints (−0.16 ± 0.07 m yr −1 ).The remaining discrepancy is presumably due to the fact that the ICESatbased results of Gardner et al. (2013) Gardner et al. (2013) seem to agree, but we believe this might be a coincidence.Previously, we argued why the Gardner et al. (2013) results might be less negative.Additionally, the results in Gardelle et al. (2013) rely crucially on an estimate of SRTM C-band penetration.Over any glacier globally, the SRTM radar waves will typically have penetrated into the snow and ice, with potential largest penetration depths through snow and firn and smallest through ice (Kääb et al., 2012;Dall et al., 2001;Rignot et al., 2001).As a consequence, SRTM glacier elevations do not, in general, reflect real mid-February 2000 glacier surface elevations but some lower horizon, the elevation of which depends on, among others, the dielectric properties and structure of the penetrated glacier volume during the SRTM campaign.For elevation difference studies where one of the elevation data sets is the SRTM, its penetration depth needs to be estimated for correction, and biases in this estimate translate directly into offsets in thickness change.Gardelle et al. (2013) Kääb et al., 2012).For eastern Nyainqêntanglha Shan this analysis indicates an average penetration of 8-10 m (7-9 m if based on the winter trends that might alternatively be assumed to reflect February conditions), much more than the 1.7 m assumed in Gardelle et al. (2013), while the corresponding off-glacier penetration is not discernible from zero.Clearly, our penetration depth lies at the high end but remains within the range of possible C-band phase-centre penetrations (Kääb et al., 2012;Dall et al., 2001;Rignot et al., 2001) In the Pamir, our results are more negative than Gardner et al. ( 2013) and in particular Gardelle et al. (2013).
As above, we suggest that our manual classification of ICE-Sat footprints versus the Randolph Glacier Inventory contributed to the difference between this study and Gardner et al. (2013) (remark: Gardelle et al. (2013) used their own inventory).Also, the difference between our study and Gardner et al. (2013) is reduced when only the results from their Method B (similar to ours) are considered.Gardelle et al. (2013) find glacier thickness changes of +0.16 ± 0.15 m yr −1 over the Pamir, whereas the present study suggests −0.48 ± 0.08 m yr −1 .Again, we find larger SRTM C-band penetration of 5-6 m compared to 1.8 m (Gardelle et al., 2013).Applying the average C-band penetration from the present study again reconciles the results of both studies.However, comparison of both studies in Pamir is complicated by a number of glacier surges (Gardelle et al., 2013) in connection with particularly sparse ICESat glacier coverage.Superimposing ICESat tracks over Landsat images and the elevation change map of Gardelle et al. (2013) reveals that they cross areas of either strongly positive or negative elevation change zones from surge waves.The ICESat trends thus become biased depending on where they sample surges, and the total ICESat sample size over Pamir is not large enough to compensate for these effects.The different observation periods for both studies (2000-2011 versus 2003-2008) may also have considerable impact due to surge activities and climate interannual variability (Yi and Sun, 2013).

Glacier mass changes and water resources
We assume an average density of 850 kg m −3 for all 2003-2008 volume changes to convert the thickness changes to water equivalent quantities (Huss (2013); see Kääb et al. (2012), for different density scenarios).The total glacier area is estimated using a simple cross product: we multiply the number of ICESat glacier footprints in each zone with the ratio between the total zone area and total number of ICESat footprints.Our method to estimate the total glacier areas is certainly open to discussion, but we prefer the above procedure over using areas from the Randolph Glacier Inventory because of the large deviations to our estimates, mainly for eastern Nyainqêntanglha and Pamir, from obviously outdated glacier outlines and voids in the Randolph Inventory (Nuimura et al., 2014).The uncertainty of water equivalent quantities includes the standard error of the elevation difference trend fit, the off-glacier trends, an error due to temporal offset of the ICESat autumn campaigns from maximum cumulative ablation conditions, an uncertainty of ±20 % for the glacier cover areas and an uncertainty of ±60 kg m −3 for density (Kääb et al., 2012;Huss, 2013).The effects of these individual sources of uncertainty, all converted to error in mass change, are combined through the root sum of squares to arrive at the total uncertainty.Note that water equivalent results from this study are not identical to Kääb et al. (2012), even when elevation difference trends agree, due to the simplified density assumption and the different glacier area estimates used.Given their fundamentally different approaches, it is challenging to discuss potential sources of disagreement between the two studies in the Himalaya and eastern Nyainqêntanglha Shan.Groundwater depletion (Rodell et al., 2009), glacier imbalance run-off into endorheic basins (Zhang et al., 2013) and errors and biases in the ICESat-derived trends as discussed above and in Kääb et al. (2012) are all likely explanations.Note that Gardner et al. (2013) offer a second, more negative gravimetric estimate for the entire combined high mountain Asia that is, however, not spatially resolved enough to compare to our results.The uncertainties of our results in this entire paragraph are given at 2σ confidence level to better agree with the uncertainty level in Jacob et al. (2012), whereas elsewhere in this contribution uncertainty is provided at 1σ confidence level.

River run-off
The glaciers of the Tarim Basin (only 40 % of its total glacier area is covered here, with notably Tien Shan missing) and the Amu Darya basin (all glacier areas covered) drain into en-dorheic basins and thus their mass changes do not contribute to sea-level changes (Table 2).The glacier mass changes in the Indus, Ganges and Brahmaputra basins from the present study contributed together ∼ 0.06 ± 0.01 mm yr −1 to eustatic sea-level rise, that is ∼ 10 % of the current sea-level contribution of 0.71 ± 0.08 mm yr −1 from glaciers outside the ice sheets (Gardner et al., 2013).
The discharge equivalent of these mass changes, that is the annual average glacier imbalance contribution to river runoff, is given in Table 2 for the major river basins covered.Note that computation of our discharge equivalents is a pure unit conversion from Gt yr −1 to m 3 s −1 , neglecting any hydrological processes and with the sole aim to roughly evaluate the relative importance of glacier mass changes for river flow in the catchments.
The Tarim Basin glaciers most likely stored water over 2003-2008 (+24 ± 33 m 3 s −1 discharge equivalent, DE).The glacier imbalance contribution to run-off is largest for Brahmaputra (−400 ± 60 m 3 s −1 DE), followed by the Indus (−220 m 3 s −1 DE), Ganges and Amu Darya (each −130 m 3 s −1 DE).Comparison of the discharge equivalent of glacier imbalance to measured river run-off is increasingly biased the further downstream the gauging stations are situated from the glaciers due to cumulative natural and manmade losses.It is also important to note that the available runoff data from literature and databases refer to various time periods, in parts considerably older than the ICESat period.Figure 2 illustrates thus only roughly the hydrological significance of the 2003-2008 glacier mass change in selected gauged catchments.For details on the gauging stations used and the uncertainty of the contributions see Supplement.As an example, the 2003-2008 glacier imbalance within the upper Indus basin at Besham Qila contributes ∼ 6 % to annual average river discharge (Fig. 2; Supplement), and we roughly estimate a very similar number for the Amu Darya (Supplement).For the upper Indus basin, the hydrological balance is under ongoing discussion (Reggiani and Rientjes, 2014), and we hope that our glacier mass change estimates can contribute towards balance closure and better understanding of spatial-temporal patterns of run-off or high-elevation precipitation amounts in the region (e.g.Immerzeel et al., 2012).
Figure 2. The percentage of discharge equivalent from annual glacier imbalance to measured average river run-off for selected catchments.Note that the actual numbers will be somewhat lower due to unaccounted water losses such as from evaporation or to groundwater.For details on the gauging stations used and the uncertainty of the contributions see the Supplement.

Conclusions
From 2003 to 2008 ICESat-derived elevation difference trends over Pamir-Karakoram-Himalaya and from comparisons to geographically overlapping studies, we draw the following conclusions: -Glacier thickness loss over the study region is most pronounced for the eastern Nyainqêntanglha Shan, followed by the western end of the Great Himalaya Range.Glaciers in and around the western Kunlun Shan are in balance or even gaining volume, and Pamir and Karakoram seem to be on the western limit of this mass balance anomaly rather than its centre.This suggests it is a meteorological or climatic anomaly (rise in precipitation).However, the cause and duration of this regional glacier anomaly is not fully understood yet.
-Our glacier volume changes seem especially uncertain in Pamir and, to a lesser extent, Hindu Kush.The heterogeneous behaviour of individual glaciers in these two zones, for instance from glacier surges, may lead to biases when extrapolating elevation difference trends from particularly sparse ICESat tracks, or areas covered by differential DEMs, to the entire zones.
-Extrapolation of ICESat trends back in time to the SRTM acquisition date suggests a much larger potential magnitude and variability of SRTM C-band phasecentre penetration than often assumed.Given the crucial importance of radar penetration for glacier thickness change studies based on radar DEMs, such as the SRTM or the upcoming TanDEM-X, we recommend being critical of penetration assumptions used in previous studies and to investigate the issue more extensively and systematically (Langley et al., 2007;chapter 7 in Müller, 2011).The problem is complicated by the fact that radar penetration has to be known specifically for certain dates from the past.
-The glacier mass changes in the Tarim and Amu Darya basins of +0.7 ± 1.0 and −4.0 ± 0.8 Gt yr −1 do not contribute to sea-level rise.The combined Ganges, Indus and Brahmaputra basin glacier mass change is −23.7 ± 2.1 Gt yr −1 , almost 10 % of the glacier contribution to sea-level rise.
-Neglecting water losses downstream of the glaciers, the 2003-2008 glacier imbalances amount to ∼ 6 % of the annual discharge of Amu Darya and upper Indus where they leave the mountains.This is a considerable amount given the significance of the rivers for the Aral Sea (Amu Darya) and massive irrigation schemes and household use in these dry climate regions.Maximum glacier imbalance contributions to annual average river run-off of up to ∼ 17 % are found for the Shyok (Indus) and ∼ 10 % for Vakhsh (Amu Darya), while minimum contributions are only ∼ 1-3 % for the monsoon-type catchments in Nepal.
-Our results on glacier mass loss agree with those from satellite gravimetry (Jacob et al., 2012) over Pamir, western Kunlun Shan and Karakoram but significantly diverge over the Himalaya and eastern Nyainqêntanglha Shan.
It is important to note that our results only cover 5 years, 2003-2008, and it remains open to what extent those years are representative for longer periods, such as the 10 years covered by Gardelle et al. (2013).For short mass balance series, single anomalous years may have large impacts on The Cryosphere, 9, 557-564, 2015 trends.Our water equivalent results are also sensitive to density and glacier area assumptions.We find that glacier outlines and areas in the study region are still quite uncertain and invite the reader to use improved glacier area estimates for upscaling our results and own assumptions for the conversion of volume changes to mass changes.
The Supplement related to this article is available online at doi:10.5194/tc-9-557-2015-supplement.
Author contributions. A. Kääb designed the study, performed the data analysis and wrote the paper.D. Treichler, C. Nuth and E. Berthier contributed to data analysis, performed supporting analyses and edited the paper.

Figure 1 .
Figure 1.Study region and trends of elevation differences during 2003-2008.Data are shown on a 1 • grid with overlapping rectangular geographic averaging cells of 2 • × 2 • .Trends are based on autumn ICESat acquisitions.Only ICESat footprints over glaciers are indicated.The zones indicated by black outlines are equivalent to the ones of Gardelle et al. (2013) with the western Kunlun Shan-Tarim zone (dashed outline) being the only additional one.Trends for all cells (coloured data circles) are statistically significant except for the cells that are marked with grey centres.The uncertainty of the temporal trends per cell is indicated through circle sizes indirectly proportional to the standard error of trends at 68 % level.
used an average C-band penetration of 1.7 m for eastern Nyainqêntanglha Shan estimated from the difference of SRTM C-band and X-band DEMs(Gardelle et al., 2012b).Here, we extrapolate our ICESat elevation trends over 2003-2008 and their uncertainty back in time to the SRTM acquisition period in February 2000.Under the coarse assumption that the 2000-2003 trends equal to those of 2003-2008, the extrapolation of February 2000 should result in a zero elevation difference to ICESat since the SRTM DEM was used as elevation reference.Offsets in this elevation difference for February 2000 are mainly attributed to SRTM radar penetration into ice and snow (for method and discussion see . Sakai et al. (2014) suggest the highest accumulation rates of the entire study region occur in eastern Nyainqêntanglha Shan, together with Hindu Kush.Correction of the Gardelle et al. (2013) results by our present C-band penetration estimate completely reconciles their results with ours.Note, however, that extrapolation of our 2003-2008 elevation difference trend back to 2000 is based on the risky assumption that the 2000-2003 trend equals the 2003-2008 trend.For the Bhutan zone, Gardelle et al. (2013) estimated a C-band penetration for February 2000 of 2.4 m, whereas our extrapolation of ICESat trends suggests around 6 m, which again reconciles the results of both studies for this zone.

Table 1 .
Gardelle et al. (2013)erence trends over the Pamir-Karakoram-Himalaya from this and other studies.Note thatGardelle et al. (2013)cover the period 2000 to ∼ 2010, while the other studies cover 2003 to 2008/2009.Note also that the zones of this study andGardelle et al. (2013)coincide whereas the zones of the other do so only roughly, which can potentially explain parts of the disagreements.See text in Sects.3 and 4 for an explanation of how the glacier areas were estimated.
a Named Hengduan Shan in Gardelle et al. (2013); b two zones of Gardner et al. (2013) overlap with our zone and both their values are given.

Kääb et al.: Contending estimates of 2003-2008 glacier mass balance over the Pamir-Karakoram-Himalaya At
are averaged from three different methods.Their results based on autumn footprints only (method B, Gardner et al., 2013) suggest a thickness change rate of −0.86 m yr −1 , which is in closer agreement with our results.a first glance, eastern Nyainqêntanglha Shan results from Gardelle et al. (2013; zone called Hengduan Shan) and www.the-cryosphere.net/9/557/2015/TheCryosphere, 9, 557-564, 2015A.

Table 2 .
Glacier thickness and mass changes over the major river basins of the study area.The discharge equivalent is a unit conversion from mass change and neglects any losses such as by evaporation or to groundwater.The Tarim Basin is endorheic.Only parts of the glacier area (∼ 40 %) within the Tarim Basin are covered in this study.b Endorheic basin. a