Estimating the volume of glaciers in the Himalayan-Karakoram region using different methods

: Ice volume estimates are crucial for assessing water reserves stored in glaciers. Due to its large glacier coverage, such estimates are of particular interest for the Himalayan–Karakoram (HK) region. In this study, different existing methodologies are used to estimate the ice reserves: three area–volume relations, one slope-dependent volume estimation method, and two ice-thickness distribution models are applied to a recent, detailed, and complete glacier inventory of the HK region, spanning over the period 2000–2010 and revealing an ice coverage of 40 775 km². An uncertainty and sensitivity assessment is performed to investigate the influence of the observed glacier area and important model parameters on the resulting total ice volume. Results of the two ice-thickness distribution models are validated with local ice-thickness measurements at six glaciers. The resulting ice volumes for the entire HK region range from 2955 to 4737 km³, depending on the approach. This range is lower than most previous estimates. Results from the ice thickness distribution models and the slope-dependent thickness estimations agree well with measured local ice thicknesses. However, total volume estimates from area-related relations are larger than those from other approaches. The study provides evidence on the significant effect of the selected method on results and underlines the importance of a careful and critical evaluation. Abstract. Ice volume estimates are crucial for assessing water reserves stored in glaciers. Due to its large glacier coverage, such estimates are of particular interest for the Himalayan–Karakoram (HK) region. In this study, different existing methodologies are used to estimate the ice reserves: three area–volume relations, one slope-dependent volume estimation method, and two ice-thickness distribution models are applied to a recent, detailed, and complete glacier inventory of the HK region, spanning over the period 2000–2010 and revealing an ice coverage of 40 775 km 2 . An uncertainty and sensitivity assessment is performed to investigate the inﬂuence of the observed glacier area and important model parameters on the resulting total ice volume. Results of the two ice-thickness distribution models are validated with local ice-thickness measurements at six glaciers. The resulting ice volumes for the entire HK region range from 2955 to 4737 km 3 , depending on the approach. This range is lower than most previous estimates. Results from the ice thickness distribution models and the slope-dependent thickness estimations agree well with measured local ice thicknesses. However, total volume estimates from area-related relations are larger than those from other approaches. The study provides evidence on the signiﬁcant effect of the selected method on results and underlines the importance of a careful and critical evaluation.


Introduction
The Himalayan-Karakoram (HK) region has among the largest glacier coverage outside the polar regions but knowledge of the dimensions of these glaciers and their behavior in response to climate change is still limited, due to their remoteness, the harsh topography, the complex political situation, and the associated difficult physical access (e.g., Bolch et al., 2012). Glaciers influence the runoff regime of major river systems, including Indus, Ganges, and Brahmaputra, which all have their sources in the HK region, and affect thus more than 700 million people (Immerzeel et al., 2010;Kaser et al., 2010). Recent studies focusing on mass changes revealed the complex and regionally heterogeneous behavior of HK glaciers (Fujita and Nuimura, 2011;Bolch et al., 2012;Gardelle et al., 2013). Furthermore, applications of modern methodologies, such as measurements of gravity field anomalies (Jacob et al., 2012) or laser altimetry , and combinations of them with field data (Gardner et al., 2013) lead to deviating findings, underlining the difficulty of measuring complex processes in such a large region. Here, we focus on determining the ice volume for the HK region, which is a basic parameter required for estimations of future glacier evolution (e.g., Le Meur et al., 2007), runoff projections (Huss et al., 2008;Gabbi et al., 2012), and estimates of future sea-level rise contributions.  Table 4 for glacier names and volumes.
For estimates of global glacier volumes, the inclusion of regions with only sparse glacier data had until recently to rely on extrapolations of glacier size distributions from regions with good data coverage (Meier and Bahr, 1996;Raper and Braithwaite, 2005;Radič and Hock, 2010). Only with the release of the Randolph Glacier Inventory (RGI) in 2012 (Pfeffer et al., 2014) a globally complete data set of glacier coverage has become available and can now be used for assessing glacier volumes without relying on data extrapolation. Volume estimations of all glaciers and ice caps of the world based on the RGI using V-A relations are given by Marzeion et al. (2012), Grinsted (2013), and Radič et al. (2013). A first estimate of the ice thickness distribution of all glaciers around the globe, based on the RGI, was presented by Huss and Farinotti (2012).
Available volume estimates for the HK region indicate large differences. However, but inconsistent delineations of regions and differing regional grouping hamper a comparative analysis. Using the same glacier inventory data set as the present study, Bolch et al. (2012) highlighted that estimates of ice volumes in the Himalayas range from 2300 to 6500 km 3 , depending on the estimation approach. Ohmura (2009) gives a value of 3800-4850 km 3 for Pakistan, India, Nepal, and Bhutan, whereas Cogley (2011) estimates the volume of all glaciers in the HK to be 3600-7200 km 3 , depending on the chosen glacier inventory. Huss and Farinotti (2012) calculate a volume of glaciers in the RGI regions (cf. Pfeffer et al., 2014) "South Asia East" and "South Asia West" (which cover also some regions adjacent to the Karakoram and Himalayas) of 4552 ± 239 km 3 , whereas for the same region Marzeion et al. (2012) obtained a volume of 5350 ± 245 km 3 , Grinsted (2013) of 6017 km 3 , and Radič et al. (2013) of 5595-6327 km 3 . Currently we lack a more systematic analysis and a comparison of ice-volume estimation methods based on a consistent data basis.
Here, we present ice volume estimations for the HK region using three different V-A relations, a slope-dependent ice thickness estimation method, and two ice-thickness distribution models. These methods are applied to the same base data, comprising the digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) and a recent glacier inventory, that to our knowledge is the most accurate glacier outline data set available for this region (Bolch et al., 2012). Results from the different approaches constrain the range within which the ice volume of HK glaciers is expected. Also the influence of different uncertainties is investigated and results of the two ice-thickness distribution models are validated with the sparse ice-thickness measurements available.

Study region and data
Here, the HK region is divided into four sub-regions, i.e., Karakoram, western, central, and eastern Himalayas (cf. Gurung, 1999;. Definitions of the extent of sub-regions are taken from Bolch et al. (2012) and visualized in Fig. 1. Cogley (2011) used the same study region and an almost identical separation between Karakoram and Himalayas, but did not further divide the Himalayas.
The glacier outlines used for the calculations have been compiled by remote sensing, based on satellite scenes acquired between 2000 and 2010. Glacier outlines were mapped by the International Centre for Integrated Mountain Development (ICIMOD) for the southern part of the Table 1. Number, area, and mean elevation of glaciers as derived with the void-filled SRTM DEM from CGIAR; see also Supplement in Bolch et al. (2012). Karakoram range, parts of the Indian Himalayas (eastern Uttarakhand and Sikkim), Nepal, and Bhutan (Bajracharya and Shresta, 2011). For the western Himalayas, outlines from , and for the Shyok basin (eastern Karakoram) from Bhambri et al. (2013), are used ( Fig. 1). For regions where no specific data sets were available, outlines for Chinese glaciers were taken from the Chinese glacier inventory (Shi et al., 2010), as available from the GLIMS database. For many regions, these outlines are also contained in the RGI, but differences exist in most parts of the Karakoram, as well as in the central and eastern Himalayas outside Nepal and China: in these regions data from ICIMOD is used (Bajracharya and Shresta, 2011), which only recently became available via GLIMS and is not represented in the RGI (version 3.2). Based on visual inspections with satellite scenes and high-resolution imagery from GoogleEarth, outlines of our inventory were found to be more accurate than the outlines contained by the RGI 3.2 in these regions. Glacier polygons smaller than 0.05 km 2 were removed as they are subject to high uncertainties and do not add much to the total volume. This inventory was already used by Bolch et al. (2012). Required topographic parameters such as minimum and maximum elevation as well as mean slope were obtained by combining the outlines with the void-filled SRTM version from the Consultative Group on International Agricultural Research (CGIAR) (Farr et al., 2007;Jarvis et al., 2008;. Table 1 provides an overview of total glacier area and mean glacier elevation per region.

Methods
In the following, the different ice volume estimation approaches applied in this study are described. For the published methodologies the descriptions are restricted to short summaries, whereas the presentation of the two ice-thickness distribution models is more extensive including newly developed modifications. Further details on the GlabTop2 approach are given in the Appendix.

Area-related thickness estimations
V-A scaling has been the most frequently used approach for ice volume estimations so far. Ice volume is calculated as a function of surface area, as large glaciers generally tend to be thicker. Area-related scaling techniques have been extensively applied for two reasons: first, their application is simple and fast. Once the scaling parameters are determined, glacier volumes can be quickly calculated for all glaciers with a known area. Secondly, area data had been measured and compiled long before digital terrain information became available. The general form of V-A scaling relation is with V representing the glacier volume, A the glacier area, and c and two scaling parameters. In order to facilitate comparisons with results from other methods and ice thickness measurements, Eq. (1) can be translated into the thickness-area relation where H is the mean ice thickness and = 1. Here we use three sets of scaling parameters, which have been applied by Cogley (2011) for the same study region: (i) Chen and Ohmura (1990), who used measurement from 63 glaciers to determine the related scaling parameters; (ii) Bahr et al. (1997), who derived the parameters in a theoretical study; and (iii) LIGG et al. (1988), who established a thickness-area relation based on ice-thickness measurements on 15 glaciers (Su et al., 1984). Although not calibrated with measurements from HK glaciers, the latter is the only available area-related parameterization that exists for High Asian glaciers, and it is also used for all ICIMOD glacier inventories (e.g., Mool et al., 2006;Bajracharya and Shresta, 2011). The applied scaling parameters are given in Table 2. Comparing our results with those from Cogley (2011) allows examining the influence of the different glacier inventory data used. Haeberli and Hoelzle (1995)  and vertical glacier extent. As this parameterization scheme has been designed to estimate glaciological parameters with detailed tabular glacier inventory data, it can be directly applied using modern glacier inventories, such as available from the Global Land Ice Measurements from Space initiative GLIMS (Raup et al., 2007) or the RGI. An application of this approach to modern glacier data can be found, for instance, in Hoelzle et al. (2007), Svoboda (2009), or Salzmann et al. (2013). Haeberli and Hoelzle (1995) calculate the volume based on mean ice thickness along the central flowline h f :

Slope-dependent thickness estimations
where ⌧ is the average basal shear stress along the central flowline, f a shape factor, g the gravitational acceleration, and ↵ the mean surface slope along the central flowline. f is chosen constant as 0.8, which is considered as the typical value for valley glaciers (Paterson, 1994). To account for the extrapolation from the mean ice-thickness along the central flowline (hf ) to the mean ice-thickness of the entire glacier (h F ), in accordance with the assumptions of a semi-elliptic cross sectional geometry and a non-branched glacier, a multiplication with (⇡/4) is added to Eq. (3): The parameterization of ⌧ (Pa) with the elevation range 1H (km) is based on reconstructed late Pleistocene glaciers of the European Alps (cf. Fig. 1 in Haeberli and Hoelzle, 1995): When applied to modern remotely sensed glacier inventories, a challenge arises from the definition of the surface slope ↵: glacier length is included in older tabular glacier inventories, such as the World Glacier Inventory (WGI), which was used by Haeberli and Hoelzle (1995). The surface slope ↵ can then be calculated with the glacier length l (hereafter called ↵ L ) and the elevation range 1H with In combination with a DEM, mean slope can also be directly derived for each glacier without knowing its length. However, mean surface slope derived form 1H and l (↵ L ) is different -generally smaller -from mean surface slope averaged over all DEM cells of a glacier (↵ DEM ). In order to determine a correction factor for ↵ DEM , glacier length has been manually determined for 50 glaciers of the study region, including the 10 largest glaciers. Resulting ↵ L values have then been compared to the ↵ DEM values (Fig. 2). Differences between ↵ DEM and ↵ L appear to be nearly constant for different glacier size classes. Based on this comparison, the correction factors to be applied to ↵ DEM have been determined as Figure 2. Differences between ↵ L (mean slope from arctan (1H /l)) and ↵ DEM (mean slope from DEM) for selected test glaciers, including the linear regression lines that are used as correction factors for glaciers > 20 km 2 (blue), glaciers between 5 and 20 km 2 (purple) and glaciers < 5 km 2 (orange). Glaciers were randomly selected, but the 20 largest glaciers of the inventory are included (note the logarithmic scale of the x axis).
10, 5, and 2.5 for glaciers with an area > 20 km 2 , 5-20 km 2 , < 5 km 2 , respectively. The standard deviations for these slope differences are 2.51, 4.16, and 2.02 for glaciers with areas > 20 km 2 , 5-20 km 2 , and < 5 km 2 , respectively, and mean differences of the three size classes are all significantly different from each other and from 0 (at a 95 % confidence level). Linsbauer et al. (2009) developed GlabTop (Glacier bed Topography), a model for assessing the spatial distribution of ice thickness by estimating the glacier depths at several points along so-called glacier branch lines. Ice thicknesses at these base points are calculated using Eqs. (3) and (5), ice thickness distribution is then derived by interpolating between these points and the glacier margins (for technical details see also Paul and Linsbauer, 2012). Here we apply a new version of this model, which is similar to GlabTop, but avoids the laborious process of manually drawing branch lines. Instead, ice thickness is calculated for an automated selection of randomly picked DEM cells within the glacierized areas. Ice-thickness distribution for all glacier cells is then interpolated from the ice thickness at the random cells and from ice thickness at the glacier margins, known to be zero. The calculation of ice thickness is grid-based and requires a DEM and the glacier mask as input (Fig. 3). Detailed explanations of how the model works and how the model parameters have been determined are given in Appendix A. To facilitate reading, this approach is called GlabTop2 in the following.

GlabTop2
The ice-thickness calculation with GlabTop2 requires estimating the parameters ⌧ (basal shear stress) and the shape factor f . Like in the slope-dependent ice thickness estimation approach, ⌧ is parameterized with the vertical glacier extent 1H (Eq. 5) and f is set to 0.8 for all glaciers (according to Paterson, 1994). With the approach from Huss and Farinotti (2012) (see Sect. 3.3.2 below), these parameters can be calculated for each glacier individually, which thus can be used for comparison.

HF-model
The second approach to calculate ice thickness distribution applied here is the model proposed by Huss and Farinotti (2012). The general idea is based on the ITEM model presented by Farinotti et al. (2009) and uses the principles of ice flow dynamics for the calculation of local ice thickness from ice volume fluxes. Local ice thickness is inverted from surface topography by calculating ice balance fluxes through cross profiles along the glacier and applying Glen's flow law (Glen, 1955). Huss and Farinotti (2012) modified the original model so that the labor-intensive digitalization of catchment areas for each glacier branch is obsolete, only glacier inventory data and a DEM are required as input data. The model can be applied to various glacier types and to different climatic regions. It was validated with thickness measurements from almost all glacierized mountain ranges around the globe (Huss and Farinotti, 2012). Hereafter, this ice thickness model is referred to as HF-model, and its structure is shortly summarized.
First, the SRTM DEM is interpolated to a metric grid with a cell size of 25-200 m depending on glacier size. Glacier hypsometry in 10 m elevation bands is derived for each glacier individually and glacier surface characteristics (mean slope, width, length) for each band are evaluated. All calculations are performed using this simplified 2-D shape. Apparent mass balance gradients for the ablation and the accumulation area (see Farinotti et al., 2009) are estimated for each glacier individually by taking into account continentality (Huss and Farinotti, 2012). Reduced surface melt rates for debris-covered glacier tongues are considered. Based on the estimated surface mass balance distribution, ice volume fluxes along the glacier are calculated. Using an integrated form of Glen's flow law, ice flux is converted to local ice thickness. The variations in the valley shape factor f , and the basal shear stress ⌧ in the longitudinal glacier profile are explicitly included in the model, and simple parameterizations account for the temperature-dependence of the flow rate factor, and the variability in basal sliding. Finally, calculated mean elevation band thickness is extrapolated to a regular grid by considering local surface slope, and the distance from the glacier margin. The model parameter values used here correspond to those of Huss and Farinotti (2012), none of the parameters were recalibrated for this study. Full details are given in Huss and Farinotti (2012). For every glacier, the HF-model provides a fully distributed ice thickness map that is directly comparable to the results of the GlabTop2 model. Furthermore, the shape factor and the basal shear stress are calculated for every elevation band.

Uncertainty and sensitivity assessment
There are general sources of uncertainties related to all approaches, including for example the correctness of glacier outlines and the DEM used. To assess the influence of these uncertainties on the resulting glacier volumes, we performed a series of sensitivity tests. As the approaches need different input data (e.g., 2-D glacier outlines or only glacier area; average slope over the entire glacier, in the surrounding of a cell, or of an elevation band), different uncertainty analyses and sensitivity tests have been performed for each method.
The influence of the scaling parameters on the results from V-A scaling is examined by comparing the results from the three applied scaling parameter sets (Table 2). Furthermore, uncertainties in the inventory-based glacier area have been examined. In the raw version of the glacier inventory used by Cogley (2011) for the same region, total glacier area is 43 178 km 2 , i.e., 5.9 % larger than in the more recent inventory used here. Considering the older reference period of Cogley's (2011) inventory, a modification of glacier areas of ±5 % is assumed as an upper bound for uncertainty in the input glacier area, which corresponds to the findings of the mapping accuracy analysis by Paul et al. (2013). Larger area differences between different glacier inventories (e.g., Nuimura et al., 2014) are normally caused by differing definitions of glaciers used for the mapping. Thus, the area of each glacier used for V-A calculations was modified by ±5 % for the sensitivity analysis.
Another important factor influencing the results of arearelated scaling methods is the chosen algorithm used for separating individual glaciers. However, for ice-divides no reference data exists, because their location depends on the purpose of the inventory and it is hardly possible to quantify the errors of glacier separations (Racoviteanu et al., 2009). Therefore, this factor is not considered for the uncertainty and sensitivity assessment. Nonetheless, the effect of using glacier complexes instead of individual glaciers can be evaluated, an issue that is a major source of uncertainty of area-related volume estimations (Huss and Farinotti, 2012;Grinsted, 2013). Many of the currently existing RGI-based estimates of global glacier volume used previous RGI versions that for many regions contained only glacier complexes rather than individual glaciers. Area-related volume calculations resulted in 70-80 % larger glacier volumes when glacier complexes instead of individual glaciers were used (Huss and Farinotti, 2012;Grinsted, 2013). Compared to this, with the HF-model resulting volume differences are only 4 % in the European Alps and 7 % in "Arctic Canada South" (Huss and Farinotti, 2012). This is, because uncertainties in glacier areas only scale linearly with these models, whereas they propagate exponentially for area-related scaling approaches. In the meantime, with the update to version 3.2, separations of glaciers in the RGI have been improved (Pfeffer et al., 2014).
For sensitivity tests of both the slope-dependent thickness estimations and the GlabTop2 model, the two scaling parameters f and ⌧ were modified. f was altered by ±0.1, which according to Paterson (1994) is a reasonable range for valley glaciers, assuming a semi-elliptic cross section. Two alternative parameterizations of ⌧ with 1H were performed; a high shear-stress version with an upper limit for ⌧ max (i.e., the basal shear-stress for all glaciers with 1H > 1600 m) of 180 kPa, and a low shear-stress version with a lower limit of ⌧ max = 120 kPa. The parameterizations of ⌧ for glaciers with 1H  1600 m have been adapted accordingly (Fig. 4). The range of maximum basal shear stresses corresponds to a change of ±20 % to the original assumption based on reconstructed late-Pleistocene glaciers, and agrees with independently derived ⌧ values from the HF-model.
For the uncertainty assessment of the HF-model, we rely on the analysis done by Huss and Farinotti (2012). They detected three main sources of uncertainty in their approach: (1) the estimation uncertainty of parameters and model simplifications, (2) the uncertainty in the DEM, and (3) the uncertainty in glacier shapes. By varying each of the five most important parameters within a physically reasonable range, they found an ice volume uncertainty for (1) of ±7.8 % for the HK region. The influence of erroneous DEM data (2) is estimated to be only ±1 % in these regions, but for individual glaciers, the influence of DEM errors can be considerably higher. Uncertainties in glacier shapes (3) are due to errors in georeferencing and a resulting overestimation of the slope, and the lack of the separation of individual glaciers in the  Haeberli and Hoelzle (1995), the modified parameterizations (⌧ max ± 30 kPa) are shown in red and orange.
RGI. This results in an estimated uncertainty of ±4-5 % in South Asia, but might be smaller for the application of the HF-model in the present study, since the glacier data set used here is considered superior to the RGI in some parts of the HK region (see Sect. 2). The total uncertainty of the ice volume estimate from the HF-model is estimated to be ±8-10 % for the HK region. See Huss and Farinotti (2012) for more details of this uncertainty and sensitivity analysis.

Evaluation and validation
V-A approaches and slope-dependent thickness estimations result in a single number of the total glacier volume or the mean ice thickness, respectively. Ground penetrating radar (GPR) measurements can thus only be used indirectly for validation, since interpolations between the GPR profiles are required. Considering the uncertainties related to interpretation and analysis of GPR data in combination with the often poor spatial coverage of GPR profiles on the entire glacier area, the comparison of estimated and 'measured' total glacier volumes is rather an evaluation than a validation. To our knowledge only two published measurement-based estimates of total glacier volume exist in the entire HK region: Gergan et al. (1999) for Dokriani glacier, and Ma et al. (2010) for Kangwure glacier, where GPR profiles were extrapolated to estimate the total glacier volume. Hence, it was not possible to perform such an evaluation. Furthermore, V-A scaling relations are designed to estimate the volume of a larger glacier ensemble (e.g., Farinotti and Huss, 2013) but are not suitable to assess the volume of individual glaciers, which further hampers their comparison with measurements.  Results from spatially distributed ice-thickness models, in contrast, can be directly validated with local ice-thickness measurements from GPR profiles. In the entire HK region, however, only a very few ice-thickness measurements are available (Bolch et al., 2012). We compared the results of the ice-thickness distribution model to published point measurements of ice thicknesses from one glacier in the Karakoram (Baltoro, measured by Marussi, 1964), one in the western Himalayas (Chhota Shigri, Azam et al., 2012), and four glaciers in the central Himalayas (Dokriani, Gergan et al., 1999;Kangwure, Ma et al., 2010;and Khumbu and Lirung, both from Gades et al., 2000). For the eastern Himalayas, no ice-thickness measurements suitable for such a validation were available. Figure 5 shows the total ice volume as derived by the six approaches for the HK region and for all sub-regions. Table 3 provides the related mean ice thicknesses and corresponding sea level equivalents (SLEs) (assuming a mean glacier density of 900 kg m 3 ). The calculated volumes for the entire HK region vary between 2955 and 4737 km 3 , depending on the method applied. By far the largest glacier volumes exist in the Karakoram region, corresponding to about 50-60 % of the total ice volume in the HK region. Also mean ice thicknesses are significantly higher in the Karakoram (94-158 m) than in the Himalayas (54-83 m). The smallest volumes are found in the eastern Himalayas (194-322 km 3 ), while the western and central Himalayas show similar glacier volumes (504-704 and 512-883 km 3 , respectively). SLEs estimates for the entire HK region range from 7.4 mm (GlabTop2) to 11.8 mm (V-A scaling with LIGG et al. (1988) parameters). Volumes and mean thicknesses for the largest glaciers of all sub-regions are given in Table 4. All applied V-A relations lead to larger volumes in all sub-regions compared to the other three approaches. The largest differences in the Himalayas (western, central, and eastern) are found between the V-A relations and the slope-dependent thickness estimation, in the Karakoram between the V-A relations and the two ice-thickness distribution models. Differences between results after LIGG et al. (1988) and results from GlabTop2 are of up to 60 % for the entire HK region and up to 67 % for the Karakoram.

Results
The slope-dependent thickness estimation and the two icethickness distribution models lead to similar results, with total absolute differences of less than 13 %. The most important deviations are found for the Karakoram, where the largest ice volume is stored and which contains the largest individual glaciers of the region. Figure 6 shows the modeled ice thickness distribution of Chhota Shigri in the western Himalayas as modeled by GlabTop2 and the HF-model. The general pattern of the ice-thickness distribution and the location of maximum ice-thicknesses are similar in both models, however local differences occur (see Sect. 4.2 for a comparison of the model results with measured ice thicknesses). For the entire glacier, these differences average out and both models yield similar mean thicknesses.
Results from GlabTop2 and the HF-model can furthermore be used to derive volumetric glacier hypsometry (Fig. 7). In the entire HK region, the elevation bands between 5100 and 5600 m a.s.l., contain the largest amounts of ice. In the east (central and eastern Himalayas) maximum volumes are at slightly higher elevations than in the west (Karakoram and western Himalayas), corresponding to mean glacier elevation (cf . Table 1), probably because the eastern glaciers are at lower latitudes. Such information can be valuable for an implementation of glaciers in regional climate models (RCMs) (Kotlarski et al., 2010). Ratios between area and volume hypsometric curves indicate ice thicknesses, for instance relatively large ice thicknesses at lower elevations in the Karakoram.

Uncertainty and sensitivity assessment
Varying the glacier areas by ±5 % alters the resulting total ice volumes calculated by means of V-A relations by ±6.9 % on average (Table 5). Modifying the parameters f by ±0.1 and ⌧ max by ±30 kPa for the slope-dependent thickness approach and the ice-thickness distribution with GlabTop2 has a combination yielding to high ice volume (f = 0.7 and ⌧ max = 180 kPa) and a combination yielding to a low ice volume (f = 0.9 and ⌧ max = 120 kPa). These parameter modifications change the total HK ice volume calculated with the slope-dependent thickness approach by +41 and 27 %, and the total volume calculated with GlabTop2 by +36 and 26 % (Table 6).

Evaluation and validation
Ice thicknesses modeled by GlabTop2 and the HF-model were directly compared to 86 ice-thickness values derived from GPR measurements at six glaciers (Fig. 8). For Chhota Shigri direct GPR data were available, for the other glacier, ice-thickness values were inferred from published GPR profiles. The average difference between the modeled and the measured ice thickness of all validation points is 25.7 m for GlabTop2 and 19.0 m for the HF-model; but average of relative differences are +6.9 and +6.5 %, respectively. The root mean square error (RMSE) of all validation points are 63.7 m for GlabTop2 and 60.7 m for the HF-model (the average ice thickness of the measurements is 180.5 m). The negative mean differences on the one hand indicate an underestimation of the ice thicknesses by both models, the positive relative differences, on the other hand, indicate that small ice thicknesses are rather overestimated. The latter to some degree is caused by glacier changes between the dates of the measurements and the acquisition dates of the DEM and the glacier outlines used by the models. In addition, errors and artifacts in the input data and simplifications and parameterizations in the models might also account for the differences. But also uncertainties related to the measurements (resolution, interpretation, and spatial reference of GPR data) and their digitalization influence the results of the validation. Another factor leading to differences between measured and modeled ice thicknesses is the comparison of local icethickness measurements with model results on a 90 m grid, which can cause large differences in particular in steep regions and yield to an overestimation of the ice thickness close to the glacier margin.
Locally large differences are observed between modeled and measured ice thickness; however, in some of these cases the two model approaches (GlabTop2 and HF) agree well among each other (e.g., Khumbu Glacier, Fig. 8). At the 86 validation points, GlabTop2 calculates ice thicknesses that on average are 6.6 m (4.1 %) smaller than calculated with the HF-model. The standard deviation of these differences is 41 m. The good agreement of these two independent methods strengthens confidence in the two approaches.

Region-specific glacier characteristics
Our results of HK glacier volumes from GlapTop2, the HFmodel and the modified slope-dependent thickness estimation are lower than the results from the V-A relations. On the other hand, total ice volumes from the slope-dependent thickness and the two ice-thickness distribution models are comparable to each other and, with one exception in the western Himalayas, always smaller than all V-A based estimates. In addition, validation with local ice-thickness mea-  surements revealed in general a good agreement with results from GlabTop2 and the HF-model. Comparability of the results form the different methods applied is hampered by the diverse nature of the approaches and the consequent different uncertainty analyses. The range of parameter modifications for GlabTop2 aims at representing the variability of f and ⌧ for different individual glaciers and within each glacier. Results from this analysis can thus be interpreted as extreme upper and lower bound estimates for the entire region. The sensitivity of V-A scaling on uncertainties in the input glacier areas is much smaller and in agreement with other uncertainty assessments of V-A approaches (Slangen and van de Wal, 2011). However, larger differences are caused by the choice of the scaling parameters: sub-region volumes calculated with parameters after LIGG et al. (1988) are on average 34 % larger than calculated with parameters after Chen and Ohmura (1990). Area-related scaling requires only two parameters that are based on available measurements. Adaption to local conditions thus depends on the availability of region-specific measurements, suitable for parameters determination, whereas slope-dependent approaches consider local conditions by directly taking into account local topography. In addition to the three area-related scaling relations presented so far, volume estimates following the V-A approaches from Marzeion et al. (2012), Grinsted (2013), andRadič et al. (2013) are given in Table 7. Results following Marzeion et al. (2012) are virtually identical to the results according to Bahr et al. (1997), as the same scaling parameters were used. Results from the V-A relation applied by Grinsted (2013) correspond for most sub-regions to the lowerbound range of the three V-A approaches applied here (between the volumes following Chen andOhmura, 1990 andBahr et al., 1997), except for the Karakoram, where the volume estimate is in the range of the results from the slope- Table 7. Ice-volume estimates (km 3 ) of the HK region according to the V-A approaches from Marzeion et al. (2012), Grinsted (2013), andRadič et al. (2013) applied to the glacier inventory used in this study.

Marzeion Grinsted
Radič et al. Our results point to a systematic overestimation of glacier volumes in the HK-region when using V-A relations (see Fig. 5 and Table 3). Differences to locally better-adapted ice thickness assessment methods can be considerable and in some cases resulting ice volumes for the sub-regions deviate by a factor of 1.6 or more. This difference is probably caused by the relatively steep slopes and small mass turnover of HK glaciers compared to the global average. Pfeffer et al. (2014) applied the three scaling relations mentioned above (Table 7)  and . However, for mountainous regions such as the HK region, V-A scaling overestimate the volumes compared to the HF-model. Compared to area-related scaling, slope-dependent approaches have the advantage to take such topographic effects directly into account and to consider region-specific conditions. Many glaciers in the Karakoram region exhibit a surgetype behavior (Hewitt, 2005;Copland et al., 2011), which might explain the larger average ice thickness in this region. Surging glaciers are not treated specifically in this study, but one should investigate whether model parameters need to be adjusted for this glacier type. However, to do so, icethickness data of surge type glaciers in both the active and the quiescent phase would be required, and such data does not exist.

Shear-stress related considerations
Results of both the slope-dependent thickness estimations and GlabTop2 strongly depend on the parameterization of the average basal shear stress ⌧ , in particular the upper limit ⌧ max , which is used for glaciers with elevation extents of more than 1600 m. As the HF-model calculates the average basal shear stress for each glacier, these values can be used for an independent comparison to the ⌧ parameteriza-tion used for GlabTop2 and the slope-dependent thickness estimation: averaged by the number of glaciers, ⌧ values used for GlabTop2 are 7.7 % lower than as derived using the HFmodel, but 4.1 % higher than in the HF-model if weighted by the glacier area (Fig. 9).
According to this comparison, the parameterization of Haeberli and Hoelzle (1995) overestimates ⌧ for HK glaciers with 1H between about 500 and 2000 m, but underestimates it for glaciers with 1H > 2500 m. Most ⌧ values calculated with the HF-model, in particular those for glaciers larger than 100 km 2 , are below 180 kPa, the value that was used for ⌧ max in the high shear-stress parameterization of the sensitivity analysis. However, an independent assessment of the plausibility of average basal shear-stress values is not possible without a representative set of measured glaciers. Taking these considerations into account together with the validation, we conclude that the results from the GlabTop2 modelruns with the modified parameters (⌧ max = 180 and 120 kPa, and f = 0.7 and 0.9, respectively) constitute upper and lower bound estimates (Table 6). On the other hand it also indicates that calculated volumes can exhibit relatively large differences on the level of individual glaciers. Since the same equation (Eq. 3) is used for the slope-dependent ice thickness estimations, these findings apply to that approach as well. , 8, 2313-2333, 2014 www.the-cryosphere.net/8/2313/2014/ Figure 10. Differences of mean ice-thicknesses from all approaches, plotted against mean slope. All HK glaciers larger than 10 km 2 are considered, they hold about 75 % of the entire glacier volume. For V-A scaling, results from the Bahr et al. (1997) approach are used because these numbers represent about the average from the different V-A approaches applied here.

The Cryosphere
Results from the slope-dependent thickness estimation method can be calibrated to give the same average ice thicknesses as obtained from V-A scaling by changing the value of ⌧ max and the related parameterization of tau (cf. Fig. 4). Only the average ice thicknesses from the approach of Chen and Ohmura (1990) can be reproduced with ⌧ max values of less than 180 kPa. To get the same average ice thickness as calculated with the scaling parameters from Bahr et al. (1997), ⌧ max values of 190 kPa are required for the Karakoram, and 170, 220, and 200 kPa for the western, central, and eastern Himalayas, respectively. Related to results from the LIGG et al. (1988) scaling parameters, these values are 200, 190, 250, and 230 kPa, respectively. For practical reasons there is a lack of average basal shear stress measurements, but the calculated mean shear stresses calculated independently by the HF-model suggest, that ⌧ values of more than 200 kPa occur only in exceptional cases (Fig. 9).

Slope-and area-related considerations
In Fig. 10, the differences of average ice thicknesses between the approaches are plotted against the mean slope of the glacier (↵ DEM ). In order to reduce the sample to a reasonable size only glaciers larger than 10 km 2 are considered. Despite representing only 2 % (571) of the total number of glaciers in the HK region, they cover about 50 % (18 860 km 2 ) of the area and hold about 75 % of the entire volume. For the V-A scaling derived thicknesses, results from the scaling parameters proposed by Bahr et al. (1997) are used for this figure, since the results from this approach are closest to the average of all three V-A scaling relations used here. Using the other approaches (not shown) gives different numbers, but the characteristics remain the same. Related to V-A scaling the following can be observed (Fig. 10 ing results in higher ice-thicknesses than the other three approaches. The first observation has strong implications, since the large glaciers strongly influence the total volume of a region. For the two largest glaciers, Siachen and Baltoro (black dots in Fig. 10), the average thickness based on V-A scaling is 166 m (+61 %) to 243 m (+57 %) greater compared to the other approaches (see also Table 4). Depending on the difference to the other three approaches, this corresponds to a volume difference of 287 to 336 km 3 only from these two glaciers. In the case of Siachen and Baltoro the difference can likely be attributed to the strongly branched character of these glaciers which cannot be adequately considered by V-A scaling approaches. Chen and Ohmura (1990) suggested estimating the volume of glaciers larger than 20 km 2 individually, because of the exponential increase of the absolute volume error with increasing glacier area. This recommendation is supported by this evaluation. The second observation indicates that the data basis for the calibration of V-A scaling parameters might be biased towards gently sloping glaciers or flat glacier parts.
Comparisons of the slope-dependent thickness estimations to the other results exhibit a strong slope dependency of the ice-thickness differences but no dependence on area. Glaciers with average surface slopes of less than 15-20 get -independently of their size -higher thicknesses with this approach compared to the results from all other methods. GlabTop2 and the HF-model, finally, lead to similar ice thicknesses. With three exceptions, all differences of average ice thicknesses from the two models are of less than 50 m.

Input data
The glacier inventory used for this study is one of the most detailed data sets currently available for this region. Previous studies had to rely older glacier data (e.g., Cogley, 2011). In many of these cases glacier accumulation areas were too large, especially if the inventories are based on the digitalization of topographic maps. Smaller glacier areas in the new glacier inventory (and resulting lower ice volume estimates) are thus rather due to improved mapping accuracy than real glacier shrinkage (Bolch et al., 2012).
Obviously, the conversion from ↵ DEM to ↵ L applied for the slope-dependent thickness estimation introduces considerable uncertainty since this approach is very sensitive to the mean slope input. Neglecting the difference of these slope values and assuming ↵ DEM = ↵ L yields to a total HK ice volume of only 1637 km 3 , i.e., less than half of the volume calculated with the correction (3360 km 3 ). Bolch et al. (2012) calculated a total volume of 2330 km 3 , using a less elaborated correction factor for ↵ DEM . Algorithms for automated glacier length calculations became available in the past years (Le Kienholz et al., 2014) but still require manual corrections; only Machguth and Huss (2014) recently presented a method for straightforward calculation of glacier lengths. Applying this approach to our glacier inventory would allow us to directly using the slope-dependent thickness estimation as suggested by Haeberli and Hoelzle (1995), without the application of slope correction factors.
Due to limited availability of ice-thickness measurements in the HK region in general, it is not possible to obtain HK-specific parameters for area-related scaling approaches. However, Huss and Farinotti (2012) provide region-specific thickness-area scaling parameters according to the best fit to their modeled glacier volumes (see their Fig. 7). Applied to our inventory data set, their thickness-area relations lead to considerably smaller ice volumes, which are comparable to our results from the slope-dependent thickness estimations and GlapTop2. The total volume for the HK region based on these scaling parameters from Huss and Farinotti (2012) is 3113 km 3 (1935 and 479 km 3 for the Karakoram and western Himalayas, respectively, calculated with the parameters for the "South Asia West" region, and 511 and 188 km 3 for the central and eastern Himalayas, respectively, calculated with the parameters for "South Asia East"). Although it is selfevident that these volumes correspond well to the estimates from the HF-model, they are independent from GlabTop2 and the slope-dependent approach. The good agreement of these methods corroborates the hypothesis that the overestimations of the area-related scaling approaches are mainly caused by the lack of measurement-based volume estimates, a finding that applies to most if not all glacierized mountain regions of the world. Therefore, volume estimations for mountain glaciers should preferably take slope directly into account since these glaciers typically have slopes steeper than the global average.

Conclusions and perspectives
Glacier volumes were calculated for all glaciers in the HK region, based on a glacier inventory containing more than 28 000 glaciers and covering an area of 40 775 km 2 , using V-A scaling, a slope-dependent mean ice thickness estimation approach, and two models for assessing the ice-thickness distribution. Results range from 2955-4737 km 3 . Mean icethicknesses are significantly higher in the Karakoram (94-158 m) than in the Himalayas (54-83 m). Area-related scaling approaches lead to higher volumes than the slopedependent thickness estimations and the two ice-thickness distribution models. The latter three approaches are in good agreement and indicate a total ice volume of 2955-3360 km 3 . These volume estimates are lower than most of the previous assessments of HK glacier volumes. Presumably, the main reason for these differences is the superior quality of the glacier inventory used here.
Comparison of mean ice thicknesses from area-related scaling and ice-thickness distribution models showed that V-A scaling yields substantially higher thicknesses for large, as well as for relatively steep glaciers. Validation of modeled ice thicknesses with direct observations generally showed good agreement. Although spatial ice-thickness distribution models require more computational effort, they provide a significant advantage since they take the three-dimensional shape of the glacier surface into account and offer the possibility for direct validation with ice-thickness measurements. Furthermore, errors in the input glacier area or separation of individual glaciers have a smaller effect on the results of ice-thickness models, since uncertainties related to area do not propagate exponentially as for the scaling approaches. Besides total glacier volume, knowledge about ice-thickness distribution is also important for several other fields of glaciology, such as hydrology, regional climate modeling, or the assessment of glacier hazards.
The generally steep surface slopes of HK glaciers, and the lack of measurement-based volume estimates of HK glaciers -that would be required for a more accurate estimate of scaling parameters -probably cause differences between V-A scaling and the other approaches that consider the surface geometry. The results of this study highlight the uncertainties related to estimates of fresh water reserves stored in the HK glaciers and their potential contribution to sea-level rise.

A1 Concept
An initial approximation of ice thickness is calculated for an automated selection of randomly picked DEM cells within the glacierized area. Ice thickness distribution for all glacier cells is then interpolated from (i) the ice thickness guesses at the random cells and (ii) the glacier margins where ice thickness is known to be zero. The procedure is explained and discussed in detail below. The calculation of ice thickness is grid-based and requires a DEM and the glacier mask as input. In a first step, all groups of glaciers sharing common borders, i.e., glacier complexes, are assigned a unified ID. All following steps are performed for one ID (i.e., all cells of a glacier complex) at a time, disregarding all cells of differing IDs. A second mask is generated where a code is assigned to all non-glacier cells directly adjacent to the glacier cells (called "glacier-adjacent cells", see Fig. 3 for a schematic illustration of the model). A different code is assigned to all glacier cells being located directly at the glacier margin (called "marginal glacier cells"). From the remaining glacier cells ("inner glacier cells") a set of random cells are drawn whereas their number corresponds to a predefined percentage (r) of the inner glacier cells. An initial buffer of 3 ⇥ 3 cells is then laid around each random cell and each individual buffer is enlarged until the difference in elevation between the lowest and the highest DEM cell falling within the buffer is equal or greater h min . Thereby all glacier cells in the buffer (marginal and non-marginal) are considered. The mean surface slope of all glacier cells in the buffer is used to calculate an initial guess of local ice thickness according to Eq. (3) and the result is assigned to the corresponding random cell, to which the buffer has been applied. Extending every buffer to a minimum elevation difference of h min avoids in most cases (in particular in mountainous terrain) very small slope values with corresponding extremely high ice thicknesses and thus makes a slope cutoff (i.e., a minimum local slope considered) redundant. From the ice thickness guesses at all random cells and an ice thickness value h ga assigned to all glacier adjacent cells, ice thickness is interpolated to all glacier cells using inverse distance weighting (IDW). The ice thickness calculation for each ID is repeated n times with different sets of random points and the n ice thickness distributions are averaged into a final estimate of ice thickness distribution.
The repeated interpolation is crucial to GlabTop2 (see the Discussion below). Only the final estimate after repeated interpolation is meant to represent ice thickness distribution and can be compared against theoretical cross sections (Figs. A1 and A2), other modeling approaches (Fig. 6), or observations (Fig. 8).

A2 Determination of model parameters
The ice-thickness calculation with GlabTop2 requires estimating the parameters ⌧ and f as well as the model parameters h min , r, h ga and n. The chosen values and their influence on the output are explained in the following with special emphasis on parameters exclusive to GlabTop2. As mentioned in the main text, ⌧ and f are set identical to Haeberli and Hoelzle (1995) and the original GlabTop1 version from Paul and Linsbauer (2012): ⌧ is calculated as a function of vertical glacier extent (1H ) according to Eq. (5), and f is set to 0.8 for all glaciers.
In GlabTop1 the mean slope ↵ is calculated over 50 m elevation intervals along the branch lines. This ensures averaging of ↵ over a horizontal reference distance which is approximately one order of magnitude larger than local ice thickness (cf. Kamb and Echelmeyer, 1986). Consequently h min = 50 m was applied here as well. The parameters n, r and h ga are not included in GlabTop1 and need calibration.
We calibrate GlabTop2 against idealized valley cross sections derived from fitting a power law to measured cross sections in formerly glaciated terrain (e.g., Svensson, 1959;Graf, 1970). By calibrating GlabTop2 against idealized shapes we avoid using the few ice thickness measurements in the HK regions for the purpose of model calibration. Instead the data are used for model evaluation and model comparison (e.g., Fig. 8). In the following we first describe the calculation of the idealized cross sections, after that the actual model calibration.
To calculate idealized glacier cross sections we measured glacier width W and surface slope ↵ from DEM data and maps of a number of real glaciers. Thereby ↵ was determined along the central flow line, over a distance of one kilometer to several kilometers for larger glaciers. along the central flow line (i.e., in the middle of the idealized cross section) was determined according to Eq. (3). The so-called form-ratio F R (Graf, 1970) was calculated according to F R = h f /W . Idealized glacier cross sections are finally calculated using a power law of the form y = a x b (Svensson, 1959), where y is the elevation of the valley side and x the horizontal distance (always measured positive) from the valley center. The factor a is used to scale the cross sections to their respective F R . The form-factor b defines the profile of the cross section, from v-shaped at b = 1, parabolic at b = 2 to flat and wide valleys with steep side walls at b 2. Because there is no generally valid form-factor, we chose to apply three different values for b (cf. Figs. A1 and A2), representing approximately the range given in the literature (cf. Svensson, 1959;Graf, 1970;Pattyn and Decleir, 1995). It needs to be kept in mind that real glacier cross sections can vary much more, depending on a variety of factors including e.g., geological properties and characteristics of present or past glaciations. In the original model by Linsbauer et al. (2012) the modeled cross sections are controlled by deciding manually how many branch lines run parallel and where to draw them (cf. Paul and Linsbauer, 2012). In GlabTop2 random points are chosen automatically and cross sectional glacier profiles are the result of the interpolation between ice thickness values of the marginal cells and guessed ice thickness at the random glacier cells (here and in the following: "guessed" refers to ice thickness values calculated in one of the n model runs while "calculated ice thickness" denoted the final output resulting from averaging the values of the n model runs). Apart from the initial guesses of ice thickness at the random points (Eq. 3), cross sections are controlled by the three model parameters n, h ga and r. The number of model runs n was set to 3 to achieve a smoother glacier bed compared to the output of a single model run (see the discussion in the following section). The choice of n has only a minor influence on the resulting total ice volume of a glacier. A minimum ice thickness in very narrow glacier sections is guaranteed by choosing an ice thickness of h ga = 15 m for the marginal cells. As for n, the influence on total ice volume is limited.
The frequency r of the random sampling controls the ratio of the number of random cells to the given number of marginal cells. The ratio of the two plays a crucial role when interpolating ice thickness from the random cells and the marginal cells to all grid cells of the glacier. If the interpolation involves a low ratio of random to marginal cells (rewww.the-cryosphere.net/8/2313/2014/ The Cryosphere, 8, 2313-2333, 2014 sulting from choosing a low r value) then the marginal cells with h f = h ga receive more weight, and the result is a shallow glacier cross section as illustrated in Fig. A1 for the case r = 0.1. If the interpolation involves a high ratio of random to marginal cells (this is the case when the chosen value of r is large), then the random grid cells with hf hga receive more weight. Hence a cross section with very steep side walls results as shown for the case r = 0.9 in Fig. A1. For these reasons the parameter r needs adjustment to achieve reasonable cross sections of valley glacier tongues. Different values for r were tested to achieve cross sections similar to the idealized ones. Best agreement resulted at r = 0.3. In Fig. A2 GlabTop2 cross sections using the above settings are compared to the idealized cross sections. The agreement is not perfect; the cross sections of wide glaciers resemble power functions with b values clearly exceeding 2.5. A better agreement could be achieved using a lower r value, but the drawback would be a general underestimation of ice thickness in narrower glacier sections.

A3 Discussion of GlabTop2
GlabTop2 involves two major changes from GlabTop1, namely branch lines become obsolete and surface slope is derived not along lines, but as the average surface slope of all grid cells within a certain buffer. These changes allow a fully automated calculation, but also result in certain conflicts with theory that are discussed below. According to Kamb and Echelmeyer (1986) Eq. (3) is intended to calculate ice thickness from surface slope measured along a central flow line. We deviate from this concept by using Eq. (3) with surface slope averaged over two-dimensional areas. The reasons for this modification are firstly that at the time of writing no data set of central flow lines was available. Secondly, an alternative experiment with automatically calculated local flow directions and slope values derived therefrom resulted in spurious variations in ice thickness because assessing local flow direction is often hampered by DEM irregularities. Glacier surfaces are usually of a rather smooth character and further experiments have shown that averaging slope over a certain surface area avoids most of the spurious variations in assessed ice thickness. By averaging surface slope over a certain area we also implicitly make the assumption that ice flow is isotropic within the area. This assumption is not valid, but the DEM data that we use need smoothing that is here achieved by directly averaging surface slope.
As explained above, only a certain fraction r of all inner glacier cells is sampled to achieve reasonable cross sections of ice thickness distribution. In contrast to manually drawn branch lines in GlabTop1, random points are sampled without control on their locations. Consequently there is a risk that the envisaged ratio of random to margin points is not given everywhere. Since we need to operate with a low ratio of r = 0.3, random points are often clustered at one location while elsewhere only a few or none exist. Directly interpolating over these points would result in a "correct" mean ice thickness but the spatial distribution of ice thickness guesses is dominated by the random clustering of the random points. Hence we repeat the calculations n times with different random sampling and each time interpolate ice thickness distribution. Eventually we average all ice thickness distributions and thereby strongly reduce the effects of the inevitable nonuniform distribution of random points.
It is shown in Fig. A2 that glacier beds tend to become rather flat for wide glacier tongues. The issue reflects the simple fact that the ratio of glacier-outline length to glacier-surface area decreases with increasing glacier area. Since r is independent of glacier size, the random cells gain more weight in the interpolation on larger glaciers. The effect on calculated ice volume outside idealized examples (Fig. A2) is difficult to quantify. For instance the comparison of GlabTop2 to the HF-method in Fig. 6 indicates no tendency towards flatter valley bottoms and steeper side walls in GlabTop2 despite of the glacier being rather wide. Furthermore Fig. 8 provides no evidence of a general overestimation of ice thickness where glaciers are thick and thus also likely wide. Then again only few measurements are available and one might need to analyze in more detail how the data are distributed within the glacier perimeters.