Ice volume estimates for the Himalaya – Karakoram region : evaluating di ff erent methods

H. Frey, H. Machguth, M. Huss, C. Huggel, S. Bajracharya, T. Bolch, A. Kulkarni, A. Linsbauer, N. Salzmann, and M. Stoffel Department of Geography, University of Zurich, Zurich, Switzerland Marine Geology and Glaciology, Geological Survey of Denmark and Greenland, GEUS, Copenhagen K, Denmark Department of Geosciences, University of Fribourg, Fribourg, Switzerland International Centre for Integrated Mountain Development, ICIMOD, Kathmandu, Nepal Divecha Center for Climate Change, Indian Institute of Science, Bangalore, India Climatic Change and Climate Impacts, Institute for Environmental Sciences, University of Geneva, Geneva, Switzerland


Introduction
The Himalaya-Karakoram (HK) region has the largest glacier coverage outside the polar regions, but the understanding of these glaciers is relatively sparse, due to their remoteness, the harsh topography, and complex political situations, all of which complicate physical access (Bolch et al., 2012).Glaciers influence the runoff regime of Indus, Ganges, and Brahmaputra, which all have their sources in this region, and affect thus more than 700 million people (Immerzeel et al., 2010;Kaser et al., 2010).A misinformation in the 4th Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR4) about the future evolution of glaciers in this region and the follow- ing debates in media (Cogley et al., 2010;IPCC, 2010;Schiermeier, 2010) emphasized both the concern as well as the poor understanding of glaciers and glacier changes in this region.Recent studies focusing on mass changes of this region revealed the complex and regionally heterogeneous behavior of HK glaciers (Fujita and Nuimura, 2011;Bolch et al., 2012;Gardelle et al., 2012).Furthermore, different findings based on various methodologies underline the intricacy of understanding complex processes in such a large region (Jacob et al., 2012;Kääb et al., 2012;Gardner et al., 2013).In this study 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), and the basis for sea-level rise projections.
Various approaches are in use to estimate glacier volumes, such as volume-area (V-A) relations (e.g., Chen and Ohmura, 1990;Bahr et al., 1997); slope-dependent ice-thickness estimations (Haeberli and Hoelzle, 1995); and more recently, a variety of spatially distributed ice-thickness models have been presented (e.g., Clarke et al., 2009;Farinotti et al., 2009;Huss and Farinotti, 2012;Li et al., 2012;Paul and Linsbauer, 2012).Such knowledge of the ice-thickness distribution is a requirement for state-of-the-art physically-based models to quantify the contribution of glaciers to runoff and sea level rise.
Until recently, estimates of global glacier volumes had to rely on extrapolations of glacier size distributions from regions with good data coverage to regions where glacier data is sparse (Meier and Bahr, 1996;Raper and Braithwaite, 2005;Radić and Hock, 2010).Only with the publication of the Randolph Glacier Inventory (RGI) (Arendt et al., 2012) a globally complete dataset 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); Huss and Farinotti (2012) used the RGI to apply a model to calculate spatially distributed ice thicknesses.Available volume estimates for the HK region exhibit a large variety; but deviating definitions and delineations of regions and different regional grouping hamper comparability.Bolch et al. (2012), using the same glacier inventory dataset as the present study, highlight that estimations of ice volumes in the Himalayas are highly uncertain and range from ∼ 2300 km 3 to ∼ 6500 km 3 .The IPCC AR4 even reports the total ice volume of the HK region to be 12 000 km 3 (Cruz et al., 2007).Ohmura (2009) gives a value of 3800-4850 km 3 for Pakistan, India, Nepal and Bhutan; Cogley (2011) estimates the mass of all glaciers in the HK to be 3600-7200 km 3 , depending on the chosen glacier inventory.Huss and Farinotti (2012)  .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.As input data, the digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) is used, in combination with the most recent and to our knowledge most accurate glacier outline dataset available (Bolch et al., 2012).Instead of developing a new volume estimation methodology, we compare results of the different existing approaches, investigate the influence of different uncertainty factors, and discuss the potentials and limitations of the different methods.

Study region and data
In the past, comparisons of ice volumes calculated for the HK region were hindered by the use of different definitions of the Karakoram and the Himalayas (Gurung, 1999;Shroder, 2011).Here we divide the HK region into four sub-regions, i. same study region and the identical separation between Karakoram and Himalayas, but didn't 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, and are as well the same than used by Bolch et al. (2012).Glacier outlines were mapped by the International Centre for Integrated Mountain Development (ICIMOD) for the southern part of the Karakoram range, parts of the Indian Himalayas (eastern Uttarakhand and Sikkim), Nepal, and Bhutan (Bajracharya and Shresta, 2011).For the western Himalayas, outlines from Frey et al. (2012), and for the Shyok basin (eastern Karakoram) from Bhambri et al. (2013), are used (Fig. 1).If no alternative datasets were available, outlines for Chinese glaciers were taken from the Chinese glacier inventory (Shi et al., 2009), as available from the GLIMS database.For many regions, these outlines are also contained by the RGI, but in most parts of the Karakoram, as well as in the central and eastern Himalayas (outside Nepal and China) they are more recent and -based on visual inspections -more accurate than the outlines contained by the RGI, as in these regions data from ICIMOD is used, which is not publicly available and therefore not represented in the RGI.
Required topographic parameters such as minimum and maximum elevation as well as mean slope were obtained by fusing the outlines with the void-filled SRTM version from CGIAR (Farr et al., 2007;Jarvis et al., 2008;Frey and Paul, 2012).Table 1 gives an overview of total glacier areas and mean glacier elevations per region.

Methods
In the following, the different ice volume estimation approaches applied in this study are described.For the existing methodologies the descriptions are restricted to short summaries, whereas the presentation of the two ice-thickness distribution models are more extensive, including newly developed modifications.

Volume-area scaling
V-A scaling has been the most frequently used approach for ice volume estimations so far.Ice volume is calculated as a function of area, as large glaciers tend to have large volumes.The reason of the popularity of V-A scaling is its simple and rapid application; once the scaling parameters are determined, glacier volumes can be easily calculated for all glaciers with a known area.Equation (1) shows the general form of V-A scaling: with V representing the glacier volume, A the glacier area and c and γ the two scaling parameters.
Here we apply three sets of scaling parameters which have been used by Cogley (2011) for the same study region: (i) Chen and Ohmura (1990), who used a set of 63 glaciers where volumes were interpolated based on radio-echo soundings to determine the related scaling parameters; (ii) Bahr et al. (1997), who confirmed these parameter values in a theoretical study; and (iii) Arendt et al. (2006)

Slope-dependent thickness estimations
Haeberli and Hoelzle (1995) presented an approach for estimating glacier volume based on the average surface slope.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, which are available for many regions in the world (cf. the Global Land Ice Measurements from Space initiative GLIMS; e.g., Raup et al., 2007) or the Randolph Glacier Inventory (RGI; Arendt  et al., 2012).An application of this approach to modern glacier data can be found, for instance, in Hoelzle et al. (2007) or Salzmann et al. (2013).Haeberli and Hoelzle (1995) calculate the volume V based on the following equation: where A is the glacier area and h f the mean ice thickness.h f is defined as where τ is the average basal shear stress along the central flowline; f a shape factor (chosen constant as 0.8); g the gravitational acceleration; and α the surface slope.
To account for the extrapolation from the mean ice-thickness along the central flowline (h f ) to the mean ice-thickness of the entire glacier (h F ), in accordance with the assumption of a semi-elliptic cross sectional geometry, a multiplication with (π/4) is added to Eq. (3): The parameterization of τ with the elevation range ∆H is based on reconstructed late-Pleistocene Alpine glaciers (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 normally 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 and the elevation range ∆H with: By using zonal statistics 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 ∆H and l (α l ) is different -generally smaller -than 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 −10 , respectively; and mean differences of the three size classes are all significantly different from each other and from 0 (at a 95 % confidence level).

Modeling of ice-thickness distributions
Based on the approach of Haeberli and Hoelzle (1995), Paul and Linsbauer (2012) developed a model for assessing the spatial distribution of ice thickness by estimating the glacier depths at several points along so-called glacier branch lines and interpolating between these points and the glacier margins.Farinotti et al. (2009) proposed an alternative approach to calculate glacier ice thickness distribution that is based on the principles of ice flow dynamics.Local ice thickness is inverted from surface topography by calculating ice balance fluxes through cross profiles along the glacier and applying Glen's (1955)

flow law.
A major disadvantage of these spatially distributed approaches is the timeconsuming preparation of the input-data, as central flow lines (Linsbauer et al., 2012) or catchment areas for each glacier branch (Farinotti et al., 2009) need to be digitized manually.Hence, we present a modified version of the approach by Paul and Linsbauer (2012), where the ice-thickness is estimated at randomly selected points on the glaciers (cf.Sect.3.3.1).For similar reasons Huss and Farinotti (2012)  oped the method by Farinotti et al. (2009) and applied it to all glaciers worldwide.The only input requirements for both approaches are glacier outline data and a DEM.

GlabTop2
The first approach to calculate distributed glacier thickness applied here is almost identical to the GlabTop model (Linsbauer et al., 2012;Paul and Linsbauer, 2012) but circumvents 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 the glacier margins where ice thickness is known to be zero.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. Figure 3 schematically illustrates the model.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).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 number 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.This procedure avoids very small slope values with corresponding extremely high ice thicknesses and thus makes a slope cutoff (i.e. a minimum local slope considered) obsolete.The mean surface slope of all glacier cells in the buffer is used to calculate local ice thickness according to Eq. ( 3) and the result is assigned to the corresponding random cell.From ice thickness at all random cells and an ice thickness value h ga assigned to all glacier adjacent cells ice 4821 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 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 is explained in the following.Parameters from Eq. ( 3) are set identical to Linsbauer et al. (2012) and Haeberli and Hoelzle (1995): for all glaciers f is set to 0.8 and τ is calculated as a function of vertical glacier extent (∆H) according to Eq. ( 5).With the approach from Huss and Farinotti (2012) (see below), these parameters can be calculated for each glacier individually, which in turn can be used for a comparison.Linsbauer et al. (2012) calculate the mean slope α over 50 m elevation intervals along the branch lines to ensure averaging of α over a reference distance of 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 r and h ga govern the shape of the glacier cross section.The latter has been measured only for a few glaciers but the cross sections of formerly glaciated valleys are generally described by a parabolic function (e.g., Graf, 1970).In the original model by Linsbauer et al. (2012) the modeled cross section is controlled by deciding manually how many branch lines run parallel and where to draw them (cf.Paul and Linsbauer, 2012).Here, random points are chosen automatically and their number and spatial distribution across the glacier-width determine the modeled cross-sectional profile of ice thickness.A high density of random points r reduces the influence of the glacier adjacent cells (set to ice thickness h ga ) in the interpolation and vice versa.Thereby the cross-section also changes from more v-shaped at low r to a parabolic profile with very steep sidewalls at high r .Both h ga and r were calibrated by adjusting calculated cross sections to theoretical parabolic glacier cross sections and set to h ga = 15 m and r = 0.3 (i.e., in each model run, ice thickness is calculated for 30 % of all inner glacier .n (number of model runs) was set to 3; however, this parameter controls the smoothness of the modeled glacier bed topography and has only a minor influence on the resulting total ice volume of a glacier.To facilitate reading, this approach is called GlabTop2 in the following.

HF-method
The second approach to calculate ice thickness distribution applied here is based on the model proposed by Huss and Farinotti (2012).The general idea is based on the ITEM model (Farinotti et al., 2009) and relies on the calculation of local ice thickness from ice volume fluxes.The original model has been significantly modified to be applicable using glacier inventory data only and to various glacier types in 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-method, 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 of the glacier.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), and reduced surface melt rates for debris-covered glacier tongues are accounted for.Based on the estimated surface mass balance distribution, ice volume fluxes along the glacier are calculated.Using an integrated form of Glen's (1955)  The model parameter values used in the present study correspond to those of Huss and Farinotti (2012), and details about all model approaches are given therein.For every glacier, the HF-method 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 point along the glacier.
In this study, hence, six volume calculations using four different approaches have been performed to estimate the ice volume of all HK glaciers: V-A scaling with the parameters from Chen and Ohmura (1990), Bahr et al. (1997), andArendt et al. (2006); slope-dependent thickness estimations following Haeberli and Hoelzle (1995); and the ice-thickness distributions models GlabTop2 (based on Paul and Linsbauer, 2012) and the HF-method from Huss and Farinotti (2012).

Uncertainty and sensitivity assessment
There are general sources of uncertainties related to all approaches, including for example the level of accuracy of glacier outlines and the DEM used.These input data uncertainties affect the results of each method, depending on how they propagate through the model.To assess the influence of these uncertainties on the resulting glacier volumes, we made a series of sensitivity tests.As the different approaches use different derivatives of the input data (e.g., 2-D glacier outlines or only glacier area; average slope over the entire glacier or local slopes within variable buffer), different 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.Furthermore, inaccuracies in the glacier outlines of the glacier inventory as a source for uncertainty in derived glacier areas 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.By considering the older reference epoch of Cogley's (2011) inventory, a modification of glacier areas of ±5 % is assumed as an upper boundary for uncertainty in the input glacier area.Thus, glacier areas used for V-A calculations were increased and decreased by +5 % and −5 % for the sensitivity analysis.
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, and two alternative parameterizations of τ with DH were performed; a high shear-stress version with an upper limit for τ max (i.e., the basal shear-stress for all glaciers with ∆H > 1600 m) of 1.8 bar, and a low shear-stress version with a lower limit with τ max = 1.2 bar.The parameterization of τ for glaciers with ∆H ≤ 1600 m have been adapted accordingly (Fig. 4).For uncertainty and sensitivity analyses of the HF model, see Huss and Farinotti (2012).

Evaluation and validation
Since direct measurement of glacier volume is not possible, glacier volumes or mean glacier thicknesses obtained from V-A approaches and slope-dependent thickness estimations can only be evaluated on an imperfect basis.Reference values build on the interpretation of ground penetrating radar (GPR) measurements and interpolations between measured profiles.Considering the uncertainties related to interpretation and analysis of GPR data, the comparison of estimated vs. "measured" total glacier volumes is rather an evaluation than a validation.However, to our knowledge only two published measurement-based estimates of total glacier volumes exist in the entire HK region (Gergan et al. (1999) for Dokriani and Ma et al. (2010) for Kangwure), hence it was not possible to perform such an evaluation.Furthermore, V-A relations are designed to estimate the volume of a larger glacier ensemble (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.This is a major advantage and reduces the uncertainties of the reference values.In the entire HK region, however, only a very limited number of ice-thickness measurements are available 4825  (Bolch et al., 2012).To evaluate our estimates, we compared the results of the icethickness distribution model to published point measurements of ice thicknesses at 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 found in the literature.

Results
Volumes were calculated for all 28 100 glaciers of the glacier dataset using six approaches: three V-A relations, the adapted approach from Haeberli and Hoelzle (1995), and two ice-thickness distribution models.Figure 5 shows the total volume as well as resulting ice volumes for all sub-regions.The calculated volumes for the entire HK region vary between 2955 and 6455 km 3 , depending on the method applied.By far the largest glacier volumes exist in the Karakoram region, amounting to about 50-60 % of the total ice volume in the HK region.The smallest volumes are found in the Eastern Himalaya (194-408 km 3 ), while the Western and Central Himalayas show comparable glacier volumes (504-1128 km 3 each).All V-A approaches lead to larger volumes in all sub-regions than the other three approaches.The largest differences are found between the V-A relation using the scaling parameters of Arendt et al. (2006), and the slope-dependent thickness estimation and the two ice-thickness distribution models, with differences of up to 118 % for the entire HK region and up to 139 % for the Karakoram.
The slope-dependent thickness estimation and the two ice-thickness distribution models, however, lead to comparable results, with total absolute differences of < 13 %.For the sub-regions, largest deviations are observed for the Karakoram, where the largest ice volumes are 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 from GlabTop2 and the HF-method.The general pattern of the ice-thickness distribution is similar in both models, however local differences occur (see Sect. 4.2 and Fig. 7 for a comparison of the model results with measured ice thicknesses).

Uncertainty and sensitivity assessment
Table 3 shows that varying the glacier areas by +5 % and −5 % alters the resulting total ice volumes calculated by means of V-A relations by +6.9 % and −6.8 %, respectively, on average.This is in agreement with other uncertainty assessments of V-A approaches (Slangen and van de Wal, 2011).Modifying the parameters f and τ for the slope-dependent thickness approach and the ice-thickness distribution with GlabTop2 has a high and a low ice-volume combination (high: f = 0.7 and τ max = 1.8 bar; low: f = 0.9 and τ max = 1.2 bar).These parameter modifications change the results of the slope-dependent thickness approach by +41 % and −27 %, and the ice-thickness distribution with GlabTop2 by +36 % and −26 % (Table 4).Huss and Farinotti (2012) found for the HF-model an overall sensitivity of total glacier volumes in the HK-region of ±10 % for typical uncertainties in the model parameters and uncertainties in the input data.

Evaluation and validation
Ice thicknesses modeled with GlabTop2 and the HF-model were compared to 86 local point ice thicknesses, derived from GPR measurements at 6 glaciers (Fig. 7).The average differences between the models and the measurements are for GlabTop2 −25.7 m and for the HF-model −19.0 m; the RMSE of all validation points are 63.7 m for GlabTop2 and 60.7 m for the HF-method.The negative mean differences indicate an underestimation of the ice thicknesses by both models, which to a certain extent 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 ice-thickness measurements with model results on a 90 m grid, which can cause large differences in particular at the steep margins of glacier beds.
Locally large differences are observed; however, in some of these cases the two model approaches (GlabTop2 and HF) agree well, but exhibit large difference to measured ice-thicknesses (e.g., validation at Khumbu glacier, see Fig. 6).The average difference between GlabTop2 and the HF-model for the 86 validation points is −6.6 m or −4.1 % with a standard deviation of 41 m.The good agreement of these two independent methods strengthens confidence in the two approaches.

Discussion
Our results of HK glacier volumes from GlapTop2 and the modified slope-dependent thickness estimation are lower than the results from the V-A relations used here, but also smaller than the estimates from other studies using such approaches (cf.Cogley, 2011;Radić and Hock, 2010;Marzeion et al., 2012;Grinsted, 2013;Radić et al., 2013).However, the results are in good agreement with volumes calculated with the HF-model and validation with local ice-thickness measurements revealed in general a good agreement.
Most assessments of the Sea Level Equivalent (SLE) stored in glaciers so far rely on statistically based volume estimations (e.g., Raper and Braithwaite, 2005;Radić and Hock, 2010).Our results point to a systematic overestimation of glacier volumes in the HK-region when using V-A relations.Differences to locally better adapted ice thickness assessment methods can be highly important and be up 100 %.Based on our model results, the SLE of the HK glaciers is 16.1 mm for the V-A relation by Arendt et al. (2006), and to 7.4 mm for GlabTop2 (assuming a density of 900 kg m −3 for the entire glacier volume and a total ocean area of 362 km×106 km).It is thus likely that the global total potential contribution of glaciers and ice caps to sea-level rise might be smaller than previously assumed, as it was also found by Huss and Farinotti (2012) at the global scale.Potential sources for this difference are (i) the input data, i.e. the glacier inventory used for the calculation; (ii) conceptual aspects of the different methodologies; and (iii) the values used for the required model parameters.These points are addressed in the following.

Input data
The glacier inventory used for this study is the most detailed dataset currently available for this region.Previous studies had to rely on extrapolation of glacier size classes (Raper and Braithwaite, 2005;Radić and Hock, 2010) while others used glacier data of older date (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).The separation of individual glacier entities has as well a large influence on volume estimates, in particular on results from V-A relations (see conceptual aspects below).
The quality of the glacier inventory input data is also of major importance for potential ice volume estimations in other regions of the world.The RGI, the most recent global-level inventory for instance, is of heterogeneous quality (Arendt et al., 2012) and especially the separation of glacier complexes into individual glacier entities is still a challenge and source of errors in some regions (e.g., Huss and Farinotti, 2012;Grinsted, 2013).Nevertheless, extrapolations of glacier size classes from one region to another are not required anymore.

Conceptual aspects
Estimating a glacier's volume based on its area has three major shortcomings: (i) mean ice thickness and area have a weak correlation.By plotting the two variables against each other, a scatter of about one order of magnitude is evident (Fig. 8).The high correlation of volume and area comes from the self-correlation in the relationship of area and volume, as volume is the product of area and mean thickness.(ii) The correlation of glacier area and thickness is rather weak (cf.Fig. 8), and modified by the surface slope of the glacier.(iii) Large uncertainties are introduced as (a) the "measured" volumes of the glaciers are based on the interpolation of GPR measurements, which in general are biased towards flat and crevasse-free glacier parts (Binder et al., 2009;Fischer, 2009); and (b) the scaling parameters are determined on the basis of only a few hundred glaciers with measurements at most.There is a variety of possible combinations for the scaling parameters c and γ, regardless whether they are calculated with a set of measured or idealized glaciers (e.g., Adhikari and Marshall, 2012;Grinsted, 2013;Farinotti and Huss, 2013).And (iv) due to non-linearity of the scaling relationships, the separation of glacier complexes into individual glaciers has major impact on resulting ice volumes: in the extreme case of a glacier complex not separated at all, an overestimation of the volume by 70-80 % can result (Huss and Farinotti, 2012;Grinsted, 2013).

Model parameters
Besides modeled regional glacier volumes, Huss and Farinotti (2012) also provide region-specific thickness-area scaling parameters according to the best fit to their calculated glacier volumes (see their Fig. 7).Applied to our inventory dataset, 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 Himalayas, respectively, calculated with the parameters for the "South Asia West" region; and 511 km 3 and 188 km 3 for the Central and Eastern Himalayas, respectively, calculated with the parameters for "South Asia East").This corroborates the hypothesis that the overestimations of the V-A approaches are mainly caused by the lack of measurement-based volume estimates of glaciers in the HK region, which is obvious as well from the large volumes resulting from the scaling parameters of Arendt et al. (2006), which were calibrated with (maritime) Alaskan glaciers.
In Table 5, volume estimates are given following the V-A approaches from Marzeion et al. (2012), Grinsted (2013), andRadić et al. (2013), applied to the glacier inventory used in this study.Results following Marzeion et al. (2012) are virtually identical to the results according to Bahr et al. (1997), as their approach was directly used.Results obtained by the V-A approach of Grinsted (2013) correspond for most sub-regions to the lower-bound range of the three V-A approaches applied here (between the volumes following Chen and Ohmura (1990) and Bahr et al. (1997), except for the Karakoram, where the volume estimate is in the range of the results from the slope-dependent ice-thickness estimations and the ice-thickness distribution models.Results obtained by the approach from Radić et al. (2013) exceed the other estimates and are only surpassed by the estimates following Arendt et al. (2006).
For the slope-dependent thickness estimations following Haeberli and Hoelzle (1995), here, mean slope had to be obtained in a different manner than originally proposed, as glacier length is not available from the glacier data set used here.Glacier length information is, however, missing in many modern glacier inventories (Paul et al., 2009), due to the lack of a suitable methodology that automatically determines glacier length without manual corrections (e.g., Le Bris and Paul, 2013).If α DEM (the slope from the DEM averaged per glacier) is directly used for Eq. ( 4) without the suggested correction (see Sect. 3.2 and Fig. 2), the resulting total ice volume is only 49 % of the volume calculated with corrected slope values (55 % of the GlabTop2 estimate), since alpha DEM is larger than the average slope along the central flowline which has to be used in the parameterization of Haeberli and Hoelzle (1995) 3 ), as different correction factors had been applied, which were based on fewer measurements.Obviously, the corrections applied here are also a generalization and thus introduce considerable uncertainties.Hence, we assume that results from this approach could be improved, if (automatically derived) glacier length values were available for all glaciers in the inventory.
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 of τ, which is used for glaciers with elevation extents of more than 1.6 km.As the HF-model calculates the average basal shear stress for each glacier, these τvalues can be used for an independent comparison to the τ-parameterization 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 HF method, but if weighted by the glacier area 4.1 % higher than in the HF-model.In Fig. 9 τ-values from the HF-model are plotted for each region.
According to this comparison, the parameterization of Haeberli and Hoelzle (1995) overestimates τ for HK glaciers with ∆H between about 0.5 km and 2 km, but slightly underestimates it for glaciers with ∆H > 2.5 km.The values of τ max used for the uncertainty assessment cover almost all τ-values as calculated with the HF-method, in particular those of glaciers larger than 100 km 2 .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 model runs with the modified parameters (τ max = 1.8 bar and 1.2 bar, and f = 0.7 and 0.9, respectively) constitute upper and lower bound estimates (Table 4).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 as well to this approach.Glacier volumes were calculated for all glaciers in the HK region, based on recent and up-to-date dataset of glacier outlines, containing more than 28 100 glaciers covering an area of 40 775 km 2 , using three V-A scaling relationships, a slope-dependent mean ice thickness estimation approach, and two models for assessing the ice-thickness distribution.Results from the latter three approaches are in good agreement and indicate a total ice volume of 2955-3360 km 3 (1683-2122 km 3 in the Karakoram, 504-543 km 3 in the W-Himalayas, 512-560 km 3 in the C-Himalayas, and 194-215 km 3 in the E-Himalayas), corresponding to a mean ice thickness of 72.5-82.5 m.These volumes are lower than estimates based on V-A relations, which are by average 50 % larger (ranging from 15 % to over 100 %, depending on the scaling parameters).These differences are caused by the lack of measurement-based volume estimates of HK glacier as well as by conceptual aspects of the different methods.The results are important in the context of improved estimates of water stored in the Himalaya region, both in view of fresh water resources and sea level rise.
Comparisons of ice-thickness measurements with modeled ice thicknesses in general showed good agreement, and thus increase the confidence in the results from the distributed ice-thickness models.Although they require more computational effort, approaches that take the three-dimensional shape of the glacier surface into account are considered as superior, since they are based on ice-mechanical and ice-dynamical considerations on the one hand, but also because they offer the possibility for validation with ice-thickness measurements at points or in GPR profiles.Furthermore, knowledge about ice-thickness distribution is important for several other fields of glaciology, hydrology or natural hazards.Due to the good agreement with local ice-thickness measurements, we have higher confidence in the results of ice-thickness distributions models such as GlabTop2 and the HF-method.By looking at these advantages and the numerous potential for further applications, we recommend using distributed ice-thickness modeling approaches for ice volume estimations ranging from a single-glacier up to the global scale, as they solely require digital glacier outlines and terrain information, which today are both globally available.Previous estimates of regional or global ice volumes, and sea level equivalent may result in systematic deviations of up to 100 % or even more, and thus may need thorough revision.
To further validate and improve ice-thickness distribution models, more ice-thickness measurement, for instance from GPR is needed, in particular in regions with sparse data coverage, such as the HK region (Bolch et al., 2012).New technologies, such as airborne GPR (Blindow et al., 2011) offer new possibilities for such measurements, which are less biased towards flat and easy accessible glacier parts and, hence, can serve as a basis for improved estimations of total glacier volumes as well.) calculated with the slope-dependent thickness estimation approach and GlabTop2, using modified parameters: f = 0.7 and τ max = 1.8 bar for the large ice-volume combination; and f = 0.9 and τ max = 1.2 bar for the small ice-volume combination.Parameterizations of τ as shown in Fig. 4, relative change to ice volumes calculated with original parameterization given in brackets.f = 0.7 and τ max = 1.8 bar f = 0.9 and τ max = 1.Table 5. 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.34 ig. 2. Differences between α l (mean slope from arctan(∆H/l)) and α DEM (mean slope from EM) for selected test glaciers, including the linear regression lines that are used as correction actors for glaciers > 20 km 2 (blue), glaciers between 5 km 2 and 20 km 2 (purple), and glaciers 5 km 2 (orange).Glaciers were randomly selected, but the 20 largest glaciers of the inventory re included.
35  are converted to a raster matching the DEM cells (red outline).Cells are discriminated as inner glacier cells (light blue), glacier marginal cells (powder-blue), glacier adjacent cells (yellow), and non-glacier cells (white).Pink cells represent randomly selected cells (r) for which local ice thickness is calculated; the blue square symbolizes the buffer of variable size, which is enlarged (dashed blue square), until an elevation extent of h min is reached within the buffer.are converted to a raster matching the DEM cells (red outline).Cells are discriminated as inner glacier cells (light blue), glacier marginal cells (powder-blue), glacier adjacent cells (yellow), and non-glacier cells (white).Pink cells represent randomly selected cells (r ) for which local ice thickness is calculated; the blue square symbolizes the buffer of variable size, which is enlarged (dashed blue square), until an elevation extent of h min is reached within the buffer.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

4815
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | calculate a volume of glaciers in South Asia East and South Asia West of 4552 ± 239 km 3 , whereas for the same region Marzeion et al. (2012) obtained a volume of 5350 ± 282 e. Karakoram, Western, Central, and Eastern Himalayas.Definitions of sub-regions extents and subdivisions are given in Bolch et al. (2012) (see also Fig. 1); Cogley (2011) used the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

4817 Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

) 4819 Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | further devel-Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | cells) 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.4823 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

4827
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

4829
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Karakoram and Western Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | . The slope-dependent 4831 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ice thickness estimations by Bolch et al. (2012) are considerably smaller (total volume of the HK region 2330 km

4833
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Study region and sources of the glacier inventory used for this study as well as by Bolch et al. (2012).

Fig. 1 .
Fig. 1.Study region and sources of the glacier inventory used for this study as well as by Bolch et al. (2012).

Fig. 2 .
Fig.2.Differences between α l (mean slope from arctan(∆H/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 km 2 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.

Fig. 3 .
Fig.3.Schematic sketch of the functionality of GlabTop2: glacier polygons (blue curved line) are converted to a raster matching the DEM cells (red outline).Cells are discriminated as inner glacier cells (light blue), glacier marginal cells (powder-blue), glacier adjacent cells (yellow), and non-glacier cells (white).Pink cells represent randomly selected cells (r) for which local ice thickness is calculated; the blue square symbolizes the buffer of variable size, which is enlarged (dashed blue square), until an elevation extent of h min is reached within the buffer.

Fig. 3 .
Fig.3.Schematic sketch of the functionality of GlabTop2: glacier polygons (blue curved line) are converted to a raster matching the DEM cells (red outline).Cells are discriminated as inner glacier cells (light blue), glacier marginal cells (powder-blue), glacier adjacent cells (yellow), and non-glacier cells (white).Pink cells represent randomly selected cells (r ) for which local ice thickness is calculated; the blue square symbolizes the buffer of variable size, which is enlarged (dashed blue square), until an elevation extent of h min is reached within the buffer.

Table 3 .
Total ice volumes (km 3 ) calculated with V-A relations, using glacier areas modified by ±5 %; relative difference to calculations with original areas in brackets.