Glacier branch lines and glacier ice thickness estimation for debris-covered glaciers in the Central Tien Shan

Abstract

Information about the ice volume stored in glaciers is of high importance for sustainable water management in many arid regions of Central Asia. Several methods to estimate the ice volume exist. However, none of them take the specific characteristics of flat terminus debris-covered glaciers into account. We present a method for deriving spatially-distributed ice thickness for debris-covered dendritic glaciers, which are common not only in Central Tien Shan but also in several other mountain ranges in High Asia. The method relies on automatically generated branch lines, observed surface velocities and surface topographic parameters as basic input. Branch lines were generated using Thiessen polygons and Dijkstra's path algorithm. Ice thicknesses for four debris-covered glaciers – South Inylchek, Kaindy, Tomur and Koxkar glaciers – have been estimated along the branch line network solving the equation of laminar flow. For Koxkar and South Inylchek glaciers, respectively, maximum thicknesses of ~250 and 380 m were estimated. These results differ by ~50 m compared with GPR measurements with an uncertainty for the debris-covered parts of ~40%. Based on geodetic mass balances, we estimate that the investigated glaciers lost between 6 and 28% of their volume from 1975 to the early 2000s.

INTRODUCTION

Large compound-valley and dendritic glaciers with numerous tributaries, debris-covered termini and complex shape are common in the Central Tien Shan (e.g. Tomur and Inylchek Glacier). Such glaciers are likely polythermal (Osmonov and others, 2013), and have been subject to several investigations regarding glacier area and glacier volume changes (Aizen and others, 2007; Kriegel and others, 2013; Osmonov and others, 2013; Pieczonka and Bolch, 2015; Shangguan and others, 2015).

Glacier ice volume has been frequently estimated with volume-area-scaling approaches that need only the glacier area as basic input (Bahr and others, 1997; Grinsted, 2013). This is accompanied by high uncertainties and no information is given regarding the spatial distribution of the ice volume (Frey and others, 2014). Few studies provide spatially-distributed ice thickness information for selected glaciers in the Central Tien Shan (Li and others, 2012; Hagg and others, 2013; Petrakov and others, 2014) or for the entire region as part of worldwide glacier volume estimations (Huss and Farinotti, 2012) using modeling approaches. These estimations are of uttermost importance for forecast models regarding water availability for sustainable water management along the northern margin of the Taklamakan desert.

Remote sensing analysis and modeling approaches provide the possibility to overcome problems of limited field accessibility in remote and rugged areas such as the Tien Shan. Several studies take advantage of remote sensing in combination with existing or estimated mass-balance data to derive glacier ice thickness estimations. Farinotti and others (2009), whose approach was further developed by Huss and Farinotti (2012), estimated ice thickness by inverting the estimated ice volume flux along the glacier relying on mass turnover, principles of ice flow dynamics and the shallow-ice approximation. McNabb and others (2012) derived ice thickness values by means of an inverse approach by solving the continuity equation with surface mass-balance rates, values of surface elevation change and surface velocities as a basic input. Other authors relate ice thickness to local topographic parameters assuming decreasing ice thickness with increasing slope angle (GlabTop-Model) (Linsbauer and others, 2012; Paul and Linsbauer, 2012) extended by Li and others (2012) to account for side drag on the valley walls. Using the same principle Frey and others (2014) calculated ice thickness values for randomly selected glacier pixels by considering local surface topography. The last two approaches need a good estimation for the basal shear stress which is estimated by the elevation range covered by the glacier system (Haeberli and Hoelzle, 1995). Clarke and others (2013) estimated the ice thickness for a sample of glaciers in Western Canada using ice extent, surface topography, surface mass balance and rate of surface elevation change as a basic input. By incorporating glacier surface velocity information, Gantayat and others (2014) derived ice thickness values for a clean-ice glacier based on the equation of laminar flow. The applicability of the presented approaches to debris-covered glaciers – with flat surface slopes and mostly low glacier velocities at their tongues – is questionable, in particular, when no slope threshold is used (Farinotti and others, 2009).

An intercomparison of different approaches for ice thickness estimation, incorporating also the above-mentioned approaches, has been done by Farinotti and others (2017) showing a mean deviation of ~10% between measured ice thicknesses and a composite solution of modeled ice thicknesses. Besides ensemble methods, which are rather more promising than one individual approach, they emphasize the importance of better accounting for uncertainties in the input data.

Glacier ice thickness is often calculated on a pixel basis (Frey and others, 2014; Gantayat and others, 2014) or along a branch line network, which covers all tributaries of a glacier and merges at confluences (Farinotti and others, 2009; Linsbauer and others, 2012). Several studies derived glacier flow lines or branch lines as surface flow lines following hydrological criteria (Schiefer and others, 2008), approximated by the glacier centerline by generating a grid-based least-cost-path (Kienholz and others, 2014), by generating a Triangulated Irregular Network (TIN) based max-min-connection of the highest and lowest point of a glacier system (Le Bris and Paul, 2013; James and Carrivick, 2016), or by constructing a grid-based path under consideration of maximum surface slope and maximum distance to the glacier margin (Machguth and Huss, 2014); the latter approaches with the focus on automatic glacier length determination. The approaches of James and Carrivick (2016), Machguth and Huss (2014) and Kienholz and others (2014) are also suitable to generate entire branch line networks. A graph-based approach has by now only been developed by Le Moine and Gsell (2015); however, focusing on glaciers of simple shape only.

The main objective of the presented study is, therefore, to extend the approach of Le Moine and Gsell (2015) to glaciers of complex shape and to estimate spatially-distributed glacier ice thickness for flat terminus debris-covered glaciers, where available GPR thickness measurements (Macheret and others, 1993) show less ice thickness than estimated by the widely used GlabTop model (Linsbauer and others, 2012).

STUDY SITE

The study is focusing on four large dendritic glaciers located in the Aksu catchment (Central Tien Shan) in the border region between Xinjiang/China and Kyrgyzstan (Fig. 1a, Table 1). Debris-covered glaciers in the Central Tien Shan cover a much larger elevation range compared with clean ice glaciers (Fig. 2). Tomur and Koxkar glaciers are located south of Jengish Chokusu (Chinese: Tomur Feng, Russian: Pik Pobeda) comprising an elevation range of 3755 and 2645 m. South Inylchek and Kaindy glaciers, in the Kyrgyz part of the Jengish Chokusu/Khan Tengri massif, have an elevation range of 4210 and 2390 m, respectively. Pieczonka and Bolch (2015) show that in the period between 1975 and 1999 glaciers in the Jengish Chokusu area experienced a significant mass loss despite their debris-cover but this mass loss was accompanied by relatively little area loss. Measurements of debris-thickness in the Central Tien Shan are sparse. For Koxkar Glacier, for instance, debris covers ~80% of the ablation area and debris-thickness increases from 0 m at 3900 m a.s.l. to over 2 m near the terminus at ~3060 m a.s.l. (Chen and Ding, 2009; Wu and others, 2013). For the same glacier, Juen and others (2014) found the strongest ablation in the elevation zone between 3750 and 3850 m a.s.l., where the thickness of debris diminishes from 10 cm to almost 0 cm.

Fig. 1. (a) Study region with debris-covered South Inylchek, Kaindy, Tomur and Koxkar Glacier. (b) Glacier velocities for 2002/03 (South Inylchek Glacier), 2010/11 (Kaindy Glacier) and 2013/14 (Tomur, Koxkar Glacier). The inset shows the location of the study region.

Fig. 2. (a) Glacier hypsometry in 100 m elevation bands (minimum and maximum elevations are depicted as red dashed lines). The percentage of each elevation level of the overall glacier area is shown as gray bars. (b) Glacier elevation range of debris-free and debris-covered glaciers in the Aksu catchment in dependence upon glacier size.

Table 1. Characteristics of the investigated glaciers

* Median elevation of the glacier.

DATA

Basic input datasets for glacier branch line and glacier ice thickness calculations are glacier outlines in vector format, a digital terrain model (DTM) covering the entire glacierized area, and a glacier velocity raster. As input datasets, we used Terra ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer), Landsat TM (Thematic Mapper) and Landsat OLI (Operational Land Imager) images for glacier surface velocity determination.

The gap-filled SRTM3 DTM provided by the Consultative Group on International Agricultural Research (CGIAR) (Jarvis and others, 2008) with a resolution of ~90 m was used for branch line generation and ice thickness estimation in combination with 2002/03, 2010/11, 2013/14 velocity information (Section 3.3, Table S1), because the original version of the SRTM3 dataset suffers from radar-related data gaps, in particular, in the accumulation regions of the investigated glaciers. Moreover, it was the only available DTM covering all glaciers with a suitable resolution and quality. Additionally, a 1975 KH9 DTM generated by Pieczonka and Bolch (2015) was used to calculate the volume changes and to relate them to the estimated overall volumes. To match the spatial resolution of the velocity raster and to preserve the information of the velocity dataset all DTMs have been bicubically resampled to 60 m. The use of cubic convolution maintains the elevation structure of the original DTM (Falorni and others, 2005).

Glacier outlines reflecting the state of the glaciers of ~2008 have been taken from the Aksu glacier inventory (Osmonov and others, 2013), which were also used for the mass-balance study of Pieczonka and Bolch (2015). The temporal mismatch of the glacier outlines in comparison with the other input data can be neglected, because the observed changes in area between 1975 and 2008 are very small or within the uncertainty range (Koxkar: −0.03 ± 0.08% a⁻¹, Tomur: −0.02 ± 0.08% a⁻¹, Kaindy: −0.14 ± 0.14% a⁻¹, South Inylchek: −0.07 ± 0.05% a⁻¹ [Pieczonka and Bolch, 2015]).

All datasets are projected to UTM (Universal Transverse Mercator) zone 44N with WGS84 as reference ellipsoid.

Glacier velocities

Glacier velocities were derived for the periods 2002/03 (South Inylchek Glacier), 2010/11 (Kaindy Glacier) and 2013/14 (Koxkar and Tomur Glacier). We used Terra ASTER Level 1A (South Inylchek), Landsat TM (Kaindy) and Landsat OLI (Koxkar and Tomur) Level 1 T data with a spatial resolution of 15 m (ASTER NIR), 30 m (TM NIR) and 15 m (OLI Pan), respectively (Table S1 in the auxiliary material). ASTER-L1A data were first orthorectified during pre-processing. The images were acquired in the late summer of the respective year, which roughly coincides with the end of the ablation season.

Velocities have been calculated with ENVI's add-on COSI-Corr (Co-registration of Optically Sensed Images and Correlation) (Leprince and others, 2007), which uses feature tracking to measure object displacements. The displacement vectors were filtered to remove outliers, potential mismatches and displacement anomalies in relation to the main flow direction (Scherler and others, 2008). In the last step, a 3 × 3 low-pass filter was applied to remove high-frequency noise. After outlier removal the uncertainty was estimated as the RMSE of all displacement vectors over nonglacierized terrain with 3.5 m a⁻¹ for 2002/03, 4.7 m a⁻¹ for 2010/11 and 3.7 m a⁻¹ for 2013/14. The spatial resolution of the final velocity raster is 60 m.

The resulting velocities show that most of the glaciers exhibit comparatively slow velocities of ~50 m a⁻¹ in the upper part of the tongues and further decreasing velocities to the termini (Fig. 1b) (Kröhnert and others, 2014; Shangguan and others, 2015). The maximum velocity of more than 150 m a⁻¹ (0.4 m d⁻¹) has been measured for the main trunk of South Inylchek Glacier. In contrast, a maximum velocity of 15 m a⁻¹ was measured on the debris-covered tongues of Tomur, Kaindy and Koxkar glaciers. These results are in line with published velocities. Using TerraSAR-X data Neelmeijer and others (2014) found velocities of ~0.05–0.1 m d⁻¹ (18–36 m a⁻¹ considering constant flow) for the lower part and 0.4 m d⁻¹ (~150 m a⁻¹) for the active upper part of the ablation region of South Inylchek Glacier for the period 2009/10. Based on ALOS/PALSAR, Li and others (2014) measured four year (2007/08, 2008/09, 2009/10, 2010/11) average surface velocities along the centerline for the active part of the ablation regions of ~0.33 m d⁻¹ (our study: 0.32 m d⁻¹ for 2013/14) for South Inylchek Glacier, 0.10 m d⁻¹ (here 0.094 m d⁻¹ for 2013/14) for Tomur Glacier and 0.07 m d⁻¹ (here: 0.075 m d⁻¹ for 2013/14) for Koxkar Glacier. Investigations by Xu and others (2011) for Koxkar Glacier could be confirmed by our results. For 2006/07 they found maximum surface velocities for Koxkar Glacier of ~120 m a⁻¹ observed in an elevation of ~4500 m a.s.l. which is close to our results of 125 m a⁻¹ for the 2013/14 period in the same elevation zone.

GPR data

For model testing, estimates of glacier thickness were extracted from previously published research. In June 2008 five transverse profiles comprising 136 single points located between 3100 and 3900 m a.s.l. were surveyed on Koxkar Glacier (Fig. 1a) using a Pulse EKKO-Pro penetrating radar system with a central frequency of 50 MHz and a bistatic antenna configuration. The GPR profiles were surveyed in common offset mode with a fixed distance between transmitter and receiver of 4 m. For thickness determination radiowave velocities of 0.167 and 0.118 m ns⁻¹ for glacier ice and supraglacial debris were used (Wu and others, 2013).

For South Inylchek Glacier thickness measurements from summer 1990 are available in Macheret and others (1993). They measured eight transverse profiles close to Merzbacher Lake and four transverse profiles and a longitudinal profile along the main trunk of South Inylchek Glacier using a high-frequency 700 MHz impulse radar TGU with an analog recording system from an oscilloscope on 36 mm film assuming a wave velocity in the ice of ~0.168 m ns⁻¹. All profiles, whose location is shown in Macheret and others (1993), were compiled from point measurements with a horizontal distance between 50 and 100 m.

The total vertical error of the radar systems can be estimated by the sum of the vertical resolution of the radar (approximated by one-quarter of the wavelength), picking errors of the ice-bedrock reflection (~25 ns) and errors due to variations in the wave velocity (±5%) (Lambrecht and others, 2014). Thus, the overall error for 100 m ice thickness is estimated at 9.8 m for the EKKO-Pro and 9.1 m for the TGU radar (≈10%).

In this study, we used the intersection points of the GPR-profiles with the glacier branch line, eight for South Inylchek Glacier and five for Koxkar Glacier, as reference points (Fig. 1a). All other points on Koxkar Glacier were used for validation and directly compared with our thickness estimates. Finally, the radar thickness data were corrected using the glacier surface elevation change rates for the 1975–99 period given by Pieczonka and Bolch (2015), which are comparable with those for the 1976–2009 period (Pieczonka and others, 2013), to match the date of the velocity measurements.

METHODS

The geometry of a glacier strongly depends on the topography and local climatic conditions. In the present study, the focus is on dendritic debris-mantled glaciers with complex glacier geometry. The utilized algorithms for glacier branch line and glacier ice thickness determination are written in Python in combination with ArcGIS and R. For the handling of raster (e.g. DTM and glacier velocity raster) and vector datasets (e.g. glacier outlines) the libraries of GDAL and OGR were used.

Glacier branch line

Glacier branch lines are important for glacier length determination (Le Bris and Paul, 2013; Machguth and Huss, 2014) and for some studies they are the basis for ice thickness determination, assuming that maximum ice thickness is concentrated along the branch line (Paul and Linsbauer, 2012; James and Carrivick, 2016).

We present a vector-based least-cost-path approach based on graph-theory to derive complete glacier branch line networks as input for glacier thickness estimates (Fig. S1 and S2). Basically, our construction of glacier branch lines relies on Thiessen polygons (Ladak and Martinez, 1996). Thiessen polygons are generated based on evenly distributed points on opposite edges of glacier polygons. These points were derived by splitting the glacier outline at regular intervals. Le Moine and Gsell (2015) only used points along the glacier outline as basic input. However, glaciers with large accumulation regions show a certain degree of simplification in their outlines, in particular, in higher regions where rock outcrops and steep headwalls are included as part of the glacier. We, therefore, generated an additional slope limited outline considering only parts of the glacier with slopes ⩽30°. Drawing on the example of South Inylchek Glacier this encompasses ~70% of the entire glacier system. The final set of input points is compiled from evenly distributed points along the original glacier outline and the slope-limited outline with a point distance of 180 m (Fig. S2).

All points are triangulated and the edges of the Thiessen polygons are derived by bisecting the sides of the triangles perpendicularly. All edges that are located inside the glacier polygon build the skeleton of the glacier which is modeled as an unweighted and undirected spatial graph of a set of nodes and links using the Python library NetworkX (Hagberg and others, 2008). In the branch line network, dangle nodes are defined as unconnected nodes of a dangling line segment (arc). These nodes are used to identify the respective glacier terminus and glacier heads where each branch line runs from. Glacier terminus is in this study defined as the dangle node with the lowest elevation.

The identification of glacier heads is not only dependent on their elevation or their distance to the end of the glacier tongue, but it is also rather necessary to spatially delimit the upper boundary of a glacier branch by means of a reliable geometric definition. Glacier heads are in a first step localized by identifying complex glacier outline segments describing a cumulative angle ⩾90° (Fig. S1 and S2). To facilitate the clear identification of complex segments, glacier outlines are smoothed by applying the Polynomial Approximation with Exponential Kernel algorithm (Bodansky and others, 2002) with a smoothing tolerance of 300 m. The value is a compromise between the level of detail and clarity of feature identification. As a result of the smoothing, we assume that there is only one distinct complex line feature and, consequently, one head per glacier branch.

All dangle nodes on arcs intersecting with complex segments are selected as potential glacier heads and tagged with the ID of the complex line feature. The dangle node with the highest elevation for each ID is finally stored as glacier head.

In the case of simple valley glaciers, there is only one branch line from the head to the terminus. Therefore, a simple shortest path algorithm would be sufficient for branch line derivation. Dendritic and compound-valley glaciers, however, consist of multiple glacier branches. Those branches are often connected in their accumulation regions resulting in more than one path connection between a glacier head and terminus. We, therefore, introduce an edge weight to the graph relying on the elevation difference covered by the edge and generated an oriented glacier branch line graph using breadth-first-search. Breadth-first-search is used to search graphs starting at a specific node (=glacier terminus) and examining all neighboring nodes at all subsequent levels until the end node (=glacier head) is found (Eppstein, 2007). The edge weights in the oriented graph are calculated using Eqn (1):

(1)$${\rm weight} = \; 1/e^{(z_{{\rm End}}-z_{{\rm Start}})/10}$$

where z _Start is the elevation of the from-node and z _End the elevation of the to-node. The weight increases exponentially when an edge runs uphill and decreases when the edge is oriented downhill.

The final branch line is then derived by finding a Dijkstra path between a glacier head and the glacier terminus. Dijkstra's algorithm returns a single-source shortest path running from the head to the terminus with minimum weight (Dijkstra, 1959). In general, Dijkstra's algorithm repeatedly selects from the set of unvisited nodes all adjacent nodes to the source or to the currently visited node and calculates the distance to the source. If the length of the new path from the source to the current node is shorter than the previously recorded distance, then the new distance is assigned as the shortest one. Originally, the output is the length of the shortest path and not the path itself. The path, that is a list of nodes, is retained by using Dijkstra's path algorithm from NetworkX.

Glacier ice thickness

We calculated ice thickness along the central branch line, based on the shallow-ice approximation, by applying the equation of laminar flow ignoring any longitudinal stress gradients (Gantayat and others, 2014):

(2)$$h = \; \displaystyle{{(n + 1)(1-z)\; u_{\rm s}} \over {2\; A\; \tau ^n}}$$

with u _s as surface velocity, n as creep parameter, A as flow rate factor, τ as basal shear stress, and z u _s as the basal velocity (Table 2) with z as the relative contribution of basal sliding to surface velocity. The flow rate factor A, primarily depending on ice temperature and water content is specified by the englacial ice temperature at the median elevation (≈ELA). We estimated the englacial ice temperature based on the mean annual air temperature (MAAT) in combination with a constant temperature offset of 7°C between MAAT and ice temperature (Huss and Farinotti, 2012). MAAT has been calculated by using the average temperature at Tien Shan Station (>3600 m a.s.l.) for the period 2000–10 (Williams and Konovalov, 2008, updated) in combination with a constant lapse rate of 0.6° per 100 m (Giese and others, 2007). This gives an englacial ice temperature of −4.6°C at the median elevation (4668 m a.s.l.) of South Inylchek Glacier and −2.7°C at the median elevation of Koxkar Glacier (4335 m a.s.l.). Accordingly, A was assumed to be 1.3E-24 Pa⁻³ s⁻¹ (Paterson, 1994) (Table 2). The creep parameter n is approximately constant and usually assumed to be 3 for glaciers (Paterson, 1994; Farinotti and others, 2009).

Table 2. Parameters used for ice thickness estimation

Basal shear stress

In many studies, glaciers are assumed to have a uniform basal shear stress between ~0.5 and 1.5 bar if the glacier is flowing over bedrock. Haeberli and Hoelzle (1995) used an empirical relation between the altitudinal range of a glacier and the basal shear stress. Carrivick and others (2016), on the other hand, calculated the mean basal shear stress using the mean of the shear stresses of pre-defined elevation bands where shear stress is related to area and slope of the respective elevation zone. In reality, local basal shear stress may deviate from this assumption due to the valley shape and the effects of basal sliding. For Greenland, Banks and Pelletier (2008) showed that basal shear stresses are heterogeneous with values from 0.5 bar to more than 2 bar. Driving stresses exceeding 1 bar are reported by Marshall and others (2011) for large glaciers and steep valley glaciers in western Canada. Iverson and others (2003) found that for debris-laden ice the spatially averaged shear stress on a rock-bed was between 1.5 and 3 bar. In general, high values of basal shear stress are expected in areas with extending flow, low values in an area with compressing or decelerating flow. In this study, the basal shear stress is calculated by Eqn (3):

(3)$$\tau = h\; \rho \; g\; f\; {\rm sin}(\alpha )$$

with ρ as ice density, g as gravitational acceleration, f as shape factor and α as surface slope.

Substituting τ in Eqn (2) by (3) gives the basic equation for thickness estimation:

(4)$$h = \; \root {n + 1} \of {\displaystyle{{(u_{\rm s}-z\; u_{\rm s})(n + 1)} \over {2\; A\; {(\rho \; g\; f\; {\rm sin}(\alpha ))}^n}}} $$

Slope α is calculated from the resampled SRTM DTM using the algorithm of Zevenbergen and Thorne (1987) which considers two cells for each gradient calculation. This approach has been proven to be robust and consistent with regard to different terrain geometries (Gonga-Saholiariliva and others, 2011). Nevertheless, SRTM slope angles are known to be overestimated in flatter topography and underestimated in steep terrain (Guth, 2006). To reduce overestimations in flat terrain as a consequence of noise, we applied a 7 × 7 low-pass filter. In contrast to the velocity dataset, where a 3 × 3 low-pass filter was applied, a larger filter size seems to be favorable as the surface is assumed to be rougher than the glacier bedrock, also with regard due to the occurrence of large ice cliffs in the flat glacier parts. Thus, particularly in the debris-covered regions with high surface roughness, the steep components are more generalized.

The slope is averaged over a vertical distance of 50 m resulting in a horizontal averaging distance between 2650 m on flat parts and 160 m in steep areas which are assumed to be several times the local ice thickness (Kamb and Echelmeyer, 1986). Glacier surface velocities are spatially averaged in the same manner.

For small slope angles Eqn (4) tends to overestimate ice thickness. Therefore, slope thresholds of 5 and 1.7° are used by Farinotti and others (2009) and Carrivick and others (2016), respectively, to avoid ice thickness overestimations in flat regions. In general, the chosen threshold should at least consider the uncertainty of the DTM. Guth (2006) showed that the slope error is highly dependent upon the steepness of the terrain. He found that SRTM in comparison with NED (National Elevation Dataset) of the US Geological Survey is too steep in low relief terrain and too smooth in steep terrain with an exponential trend. He proposed a 5% (=2.86°) slope cut-off to exclude regions where SRTM is expected to be too noisy. De Vente and others (2009) found that SRTM DTMs show more underestimations of slope gradients than overestimations. Based on the results of Guth (2006), we approximated the slope error dα with respect to the local terrain steepness α (in radians; Eqn (5)).

(5)$${\rm d}\alpha = 0.0319\; \log (\alpha ) + 0.071$$

The lower limit is defined where the slope angle α is larger than the slope error |dα| and was finally chosen with 2.1°.

For glacier velocities, the lower boundary was set in accordance with the respective uncertainty estimate du _s (2002/03: 3.5 m a⁻¹, 2010/11: 4.7 m a⁻¹, 2013/14: 3.7 m a⁻¹).

Basal sliding

Besides the basal shear stress, the contribution of the basal sliding velocity to the glacier surface velocity is one of the main unknowns in Eqn (2). Glaciers in the Central Tien Shan are assumed to be mostly polythermal (Osmonov and others, 2013). Temperate glaciers are, in general, characterized by a more uniform velocity distribution across the glacier, whereas in case of cold glaciers internal deformation is the dominating process with a velocity distribution across the glacier of parabolic shape. Velocity investigations on the debris-covered Fedchenko Glacier showed a contribution of basal sliding of ~10% to the surface velocity at the end of the ablation season (Lambrecht and others, 2014).

To infer the role of basal sliding summer (17 July 2010–17 October 2010) and winter (17 October 2010–04 March 2011) surface velocities given by Li and others (2014) for Koxkar and Inylchek glaciers were taken into account. For Koxkar Glacier the summer and winter velocities are 10% higher and 24% lower, respectively, than the annual mean velocities. According to the seasonal surface velocities given by Li and others (2014) a first estimate of the contribution of the basal sliding velocities is 33% for Koxkar Glacier (summer velocity 8 cm d⁻¹, winter velocity 6 cm d⁻¹) and 9% for South Inylchek Glacier assuming that there is no sliding during winter and no significant seasonal variation in ice deformation (Copland and others, 2003). This first estimate could not be proven based on 12 GPR-points, mainly distributed on the flat terminus of South Inylchek (7 GPR-points) and Koxkar (5 GPR-points) glaciers (Fig. 1), and Eqn (2) resulting in a contribution between 86% (z = 0.86) and 99% (z = 0.99) with regard to the glacier surface velocity at the flat terminus (Fig. 3).

Fig. 3. Empirical function to estimate the contribution of basal sliding to the glacier surface velocity (z) based on the ratio of the total glacier area and the glacier area upstream of a particular location assuming no sliding in the upper parts.

Rapid sliding at low gravitational driving stress, in winter and summer, as a consequence of increased meltwater production and converging tributary glaciers seems to play a major role in the case of South Inylchek Glacier; however, the influence of basal sliding usually varies spatially and decreases towards higher elevations. The point with a contribution of 86% was located at the lower ablation area of Koxkar Glacier where some artifacts in the velocity field are visible (Fig. 1b). Finally, we estimated the contribution of the sliding velocity for all four glaciers relying on the estimates of South Inylchek and Koxkar glaciers by using a sigmoid function in dependence upon the ratio of the total glacier area and the glacier area upstream of a particular location assuming no sliding in the upper parts where cold ice is assumed to predominate (Fig. 3).

For comparison, glacier ice volumes were also estimated using volume-area-scaling with the parameters defined by Su and others (1984) (denoted as Volume_Su) and the GlabTop approach (denoted as Volume_Lins) (Linsbauer and others, 2012; Paul and Linsbauer, 2012). For model calibration Su and others (1984) used gravimetric and radar ice thickness measurements of a sample of Tien Shan glaciers. Additionally, mean ice thickness, relying on a model that incorporates ice flow mechanics, were provided by Huss and Farinotti (2012) and compared with our results.

GlabTop is based on Eqn (3) only (where ρ is the ice density, g the acceleration due to gravity, α the surface slope) assuming an elevation range dependent constant basal shear stress τ along the central branch line. Because Koxkar, Tomur, Kaindy and South Inylchek glaciers cover an elevation range of more than 2000 m τ was 1.5 bar according to the empirical relation given by Linsbauer and others (2012).

Drag on the margins, relevant for high mountain glaciers which are usually bounded laterally, was considered by the empirical shape factor f (Table 2). For ice thickness measurements given in Macheret and others (1993) for South Inylchek Glacier and in Xie and others (2007) for Koxkar Glacier the shape factor varies between 0.65 and 0.85. We used the median value of 0.75 for all glaciers as we are compensating for ice thickness overestimations in flat terrain by a corresponding slope threshold.

Spatial interpolation

Based on cross-profiles just below the end of each glacier we found that a parabola with an exponent i of 2.5 for Koxkar, Tomur and Kaindy glaciers and 6.5 for South Inylchek Glacier is the best approximation of the glacier cross-section. The form of the parabola can be written as:

(6)$$y = h-a\; x^i$$

where h is the ice thickness at the branch line, a is a scaling factor, y is the ice thickness along the cross-section and x the distance between the branch line and the glacier outline. As input 192, 722, 340 and 788 profiles with an average horizontal distance of ~300 m were generated for Koxkar, Tomur, Kaindy and Southern Inylchek glaciers, respectively, with maximum ice thickness on the branch line and zero at the glacier margins. For spatial interpolation, the Ordinary Kriging algorithm implemented in ArcGIS 10 was utilized. To avoid unrealistic peaks, typical for distance-based interpolation schemes, the point density along the main branch line has been increased to 25 m by linear interpolation because lower point density than the final resolution of the thickness map was the main issue mentioned by Seroussi and others (2011) when using spatial interpolation schemes like kriging.

Glacier ice thicknesses for all glaciers under investigation have finally been estimated using Eqn (4) with the resampled SRTM DTM and 2002/03 (South Inylchek), 2010/11 (Kaindy), and 2013/14 (Koxkar and Tomur) glacier surface velocities as input. Because the layer of debris is comparatively thin in comparison with the estimated ice thickness, we assume that the specified volumes coincide with the total amount of stored ice.

Sensitivity and uncertainty

We focused our sensitivity analysis on exploring the influence of the model parameters on the estimated ice thickness. The sensitivity of Eqn (4) is investigated for 13 GPR points (Fig. 1a) by varying the flow rate factor A, the shape factor f, the ice density ρ and the contribution of basal sliding to the surface velocity z.

The parameter sensitivity S is determined as the ratio of the variation of the ice thickness as model output O to the variation of the relevant input I parameter (Eqn (7)). The variation is calculated by the ratio of the standard deviation (STD) to the mean value (AVG) of the respective parameter (Misra and Rose, 1996). This analysis is 1-D assuming independent input variables and neglecting parameter interactions.

(7)$$S = \displaystyle{{({\rm ST}{\rm D}_{\rm O})/({\rm AV}{\rm G}_{\rm O})} \over {({\rm ST}{\rm D}_{\rm I})/({\rm AV}{\rm G}_{\rm I})}}$$

The STD of the input parameters was prescribed by generating a sequence of ten values uniformly spread over a range bounded by the expected minimum and maximum value of the relevant parameter (Table 3).

Table 3. Sensitivity of glacier ice thickness to four different input parameters

The range of the flow rate factor A is specified by the englacial ice temperature at the median elevation (≈ELA) and at the glacier terminus. The lower boundary of the flow rate factor A is 9.3E-25 s⁻¹ Pa⁻³ (ice temperature −4.6°C). The upper boundary is the flow rate factor for temperate conditions of 2.4E-24 s⁻¹ Pa⁻³ (Cuffey and Paterson, 2010). The range of the factor z is given by the lower estimates assuming no sliding during winter and the maximum sliding velocity at GPR points derived by Eqn (4).

S = 1 indicates a linear relation between input and output items. S > 1 is an indication for greater sensitivity, whereas S < 1 indicates less sensitivity in comparison with a direct linear relation.

From Table 3 it can be seen that the ice thickness estimates are less sensitive to the flow rate factor A. With regard to the tested range of the relative contribution of basal sliding to surface velocity (z) ice thickness is sensitive to the contribution of the basal velocity when the glacier surface velocity is low (>0.7). The sensitivity of z decreases with increasing surface velocity (<0.5).

The uncertainty of the ice thickness values is estimated with regard to thickness measurements on Koxkar Glacier (yellow crosses in Fig. 1a) and a longitudinal GPR-profile on South Inylchek Glacier (Fig. 4b) not used for calibration. All measurements are located on the main trunk of the respective glacier.

Fig. 4. (a) Glacier branch lines for Koxkar, Tomur, Kaindy and South Inylchek glaciers, (b) Glacier ice thickness for South Inylchek, Kaindy, Tomur and Koxkar glaciers. I1 indicates the position of the longitudinal profile measured by Macheret and others (1993).

RESULTS

Glacier branch line

Relying on graph-theory, we derived branch lines for each of the four glaciers. The proposed algorithm successfully handles problems related to DTM artifacts introduced as a result of gap-filling, as well as inaccuracies in the glacier outlines. The resulting network is visually complete and contains all glacier branches and tributary glaciers (Fig. 3a). We estimated a maximum glacier length of 26.0, 40.7, 26.7 and 65.1 km for Koxkar, Tomur, Kaindy and South Inylchek glaciers based on the Aksu glacier inventory. In terms of the course of the branch line, the proposed method produces consistent branch line networks useful for glacier ice thickness estimations.

Glacier ice thickness

Total volume for the four investigated glaciers was found to be ~99 km³, corresponding to an average ice thickness of ~112 m (Fig. 4, Table 4).

Table 4. Glacier ice thickness and glacier ice volumes for investigated glaciers.

$\overline {h_{\rm D}} {\rm \;} $ Mean ice thickness for the debris-covered glacier part,

$\bar{h}$ Mean ice thickness for the entire glacier,

Max h _D Maximum ice thickness for the debris-covered glacier part,

Glacier ice thickness is highest at the debris-covered glacier tongues with an average thickness of ~72 m at Koxkar Glacier and 136 m at South Inylchek Glacier (Table 4). The mean ice thickness for the entire glacier is highest for Tomur and South Inylchek Glacier, while Kaindy Glacier reveals a lower mean thickness of ~77 m. For most of the glaciers, the majority of the ice is concentrated close to the median elevation. About 60% of ice volume is stored in the ablation regions. Using GlabTop (Eqn (3)) the amount of ice stored in the ablation region increases to 70–80% (Fig. 5).

Fig. 5. Percentage of ice volume stored in the accumulation and ablation region of the investigated glaciers (grey area). The black line shows the distribution according to Eqn (3) (GlabTop). The red line shows the percentage of glacier area per elevation.

For a sample of glaciers south of Tomur Peak Pieczonka and others (2013) showed that the specific mass budget for the 1975–2009 period (−0.35 ± 0.15 m w.e. a⁻¹) is approximately consistent with the 1975–99 mass budget (−0.42 ± 0.23 m w.e. a⁻¹). Thus, we estimated the relative volume change for the four investigated glaciers taking the 1975–99 geodetic mass balances given by Pieczonka and Bolch (2015). For Koxkar and Tomur Glacier, we found a relative volume loss of 17 and 28% for the 1975–2013 period. For Kaindy and South Inylchek Glacier the volume loss was 14% (1975–2010) and 6% (1975–2002).

DISCUSSION

Glacier branch line

The presented algorithm for branch line generation produces complete branch line networks for large dendritic glaciers of complex shape in a fully automated way and needs no user intervention. Shortest paths branch lines, under consideration of minimizing the cumulative elevation-costs of a branch line link between pairs of nodes, have the advantage that all paths in a glacier system have the same endpoint. So, it is easily possible to extract the longest branch line or the average length of all branch line paths in a glacier system for glacier length analysis.

The shape and completeness of the branch line network strongly depend on the quality of the underlying glacier outline which also affects the clearness of the identification of glacier heads used as branch line origin. The presented approach is resistant against convex cross-sections, typical for lower parts of the glacier tongue, which are causing deflections from the center position towards the glacier margin. Depending on the quality of the DTM (amount of data gaps, void filling) approaches following hydrological criteria would likely perform better in the higher reaches. For branch lines following hydrological criteria it has been shown that maximum path length is expected to be 10–15% longer than centerline based branch line lengths (Schiefer and others, 2008).

Glacier length estimates are already available as an attribute in the widely used Randolph Glacier Inventory (RGI) derived by the approach of Machguth and Huss (2014) which combines the steepest descent and greatest distance from the glacier margin. Because our branch lines are based on the Aksu glacier inventory we generated glacier branch lines also for the RGI 5.0 glacier outlines for validation.

To measure the similarity of the two branch lines connecting the highest and lowest points ($L_{(E_{{\rm min}}-E_{{\rm max}})}$ and L _LeBris (Le Bris and Paul, 2013)), we created buffers of different size around L _LeBris with distances of 30, 60 and 90 m and calculated the amount of overlap. For Kaindy Glacier 86.8% of the branch line $L_{(E_{{\rm min}}-E_{{\rm max}})}$ overlaps with the 90 m buffer of L _LeBris, whereas the overlap is only 52% for Tomur Glacier (Table 5). This is also evident in the smaller glacier length difference of 0.4 km compared with 3.1 km. However, lines connecting the highest and lowest point of a glacier system do not necessarily represent the longest branch line (Fig. S3). For Kaindy Glacier the length L _max of the longest branch line is ~12% longer than the branch line between the highest glacier head and the glacier terminus. The approach of Le Bris and Paul (2013) may, therefore, underestimate the glacier length. Compared with grid-based methods (L _RGI) our graph-based approach gives comparable length estimates with a ratio of ~1.0 (Table 5).

Table 5. Branch line characteristics in comparison with branch line lengths derived by the approach of Le Bris and Paul (2013) and to maximum lengths specifications in the RGI (Randolph Glacier Inventory) Version 5.0 (Arendt and others, 2015)

$L_{(E_{{\rm min}}-E_{{\rm max}})}$ Length of branch line connecting the highest and lowest point,

L _LeBris Length of branch line based on Le Bris and Paul (2013),

L _max Length of the longest branch line,

L _RGI Glacier length specified in RGI Version 5.0.

Glacier ice thickness

The presented approach is focused on debris-covered glaciers of complex shape; however, not accounting for mass conservation. The thickness estimations neglect the layer of debris on-top as the layer of debris is thin in comparison with the glacier thickness. The method for ice thickness estimation follows the idea of Gantayat and others (2017) with the difference of assuming spatially varying basal velocities and avoiding time-consuming branch line digitization.

In general, volume and area distribution show the same pattern (Fig. 5). For Kaindy Glacier most of the volume is concentrated around 4500 m where most of the glacier area is located. For South Inylchek Glacier most of the volume is located around and below the ELA, whereas the accumulation area represents the largest part of the glacier. Here, uncertain glacier velocities due to the snow coverage in the accumulation area might be one reason for the comparatively low amount of ice volume stored in the higher parts of the glacier.

Comparison with other approaches

The estimated mean ice thickness for the four investigated glaciers is lower compared with the mean ice thicknesses estimated with the scaling approach of Su and others (1984) and also lower compared to the results of the flux- and stress-driven approaches (Table 6). South Inylchek Glacier, in particular, shows much lower mean thickness.

Table 6. Glacier ice thickness in comparison with other studies.

Bold is used to emphasize the results from this study.

$\bar{h}$ Mean ice thickness,

$\overline {h_{{\rm Huss}}} $ Mean ice thickness (Huss and Farinotti, 2012),

$\overline {h_{{\rm Lins}}} $ Mean ice thickness based on the approach of Paul and Linsbauer (2012),

$\overline {h_{{\rm Su}}} $ Mean ice thickness [$\bar{h} = -11.32 + 53.21{\rm *Are}{\rm a}^{0.3}$, Su and others, 1984],

V _Lins Ice volume based on the approach of Paul and Linsbauer (2012),

V _Su Ice volume based on volume-area-scaling (Su and others, 1984).

We also compared our estimated glacier volumes with glacier volume estimations derived with the approach of Paul and Linsbauer (2012) (GlabTop – Eqn (3)). As a result, we found GlabTop glacier volumes to be between 34 and 41% (34% Koxkar, 34% Tomur, 41% Kaindy, 34% South Inylchek) higher than our results. From Figure 6a it is clearly visible that at the lower parts, where all of the available GPR points are located, GlabTop significantly overestimates ice thickness. An underestimation of GlabTop in the accumulation regions and an overestimation of ice thickness in the lowermost parts has also been observed by Petrakov and others (2016) for the Ak-Shirak massif. GlabTop was developed for Swiss glaciers and its suitability has been shown for alpine glaciers as well as for non-debris-covered polythermal glaciers in the Himalaya (Frey and others, 2014). The differences in the total volume between Volume_Lins and our results are assumed to be attributed to an underestimation of the average basal shear stress for dendritic debris-covered glaciers questioning the empirically set upper-bound value of 1.5 bar by Linsbauer and others (2012).

Fig. 6. (a) Koxkar GPR ice thickness measurements vs. thickness estimates based on Eqn (4) (This study) and Eqn (3) (GlabTop). (b) Longitudinal profile along South Inylchek Glacier (Profile I1, Fig. 4b) showing thickness measurements from Macheret and others (1993) and thickness estimates from this study.

Area-related scaling approaches have frequently been used. For example, Su and others (1984), Hagg and others (2013) and Petrakov and others (2016) estimated volume of Tien Shan glaciers; however, scaling relations with locally different regression parameters might adequately apply only to larger ensembles of glaciers with accurate outlines but not to single glaciers where uncertainties of ~50–60% have been reported by several authors (Marshall and others, 2011; Grinsted, 2013). The volumes based on the relation proposed by Su and others (1984) are 2.0–2.9 times higher than our volume estimations and 1.3–1.9 times higher than volume estimations using GlabTop (Table 6). This is comparable with Petrakov and others (2016) who found that calibrated volume-area-scaling approaches overestimate glacier volume in comparison with GlabTop by a factor of 1.5–1.7.

Our thickness estimates along the main trunk of South Inylchek Glacier are higher than the observations of Macheret and others (1993) who measured a thickness along the glacier tongue of ~150–250 m that increases to ~300–320 m adjacent to Merzbacher Lake (Fig. 6b). In comparison with Avsiuk and Kotlyakov (1967), who, based on seismosoundings, specified the ice thickness of South Inylchek to 100–120 at the lower parts of the tongue up to 400 m in the main part of the glacier with a mean ice thickness of 200–300 m, our estimates are lower with a mean thickness at the debris-covered glacier tongue of ~136 m. This can mainly be attributed to the time gap between both studies and a mean elevation change rate of −1.0 to −1.4 m a⁻¹ (Pieczonka and Bolch, 2015). Modeling results of Huss and Farinotti (2012) show an overall mean thickness for the entire South Inylchek Glacier system of 233 m (Table 6), which is probably slightly too high taking into account that the mean ice thickness along the main trunk of the glacier is ~200 m with regard to GPR measurements of Macheret and others (1993) and that ice thickness, in general, decreases towards the steep accumulation regions.

In general, the spatial distribution of ice thickness is comparable with that of Fedchenko Glacier, a dendritic glacier whose tongue is also debris-covered, with an ice thickness of 300 m at the lower part with gentle slopes between 1.5° and 2.5° up to 800 m in the middle section (Avsiuk and Kotlyakov, 1967; Lambrecht and others, 2014).

Uncertainty and sensitivity

Drawing on the example of Koxkar and South Inylchek glaciers, where thickness measurements for uncertainty estimation are available, the uncertainty for the estimated ice thickness values is high with approximately 54% for Koxkar Glacier and 32% for South Inylchek Glacier. Here it must be taken into account, that, on the one hand, the available GPR points are located only on the main trunk of South Inylchek and Koxkar Glacier without any information for the higher reaches. On the other hand, GPR points have sparse spatial coverage and a different spatial footprint compared with our estimates. In general, the uncertainty is comparable with other studies, taking the often high uncertainty ranges into account. ITMIX showed that in the case of glaciers only, shear-based and simple velocity-based approaches have average deviations in the order of 4 ± 72% and −16 ± 46% (Farinotti and others, 2017).

For Sary-Tor Glacier in the Ak-Shirak massif, Petrakov and others (2014) found a good agreement between radio-echo sounding derived ice thicknesses and thickness estimates by applying the GlabTop model with a difference of only 1% when using a shape factor of 0.6. For Koxkar Glacier the RMSE between the GPR thickness measurements and the ice thickness estimates based on Eqn (3) is 97.2 m (160% of the mean GPR thickness). Using equation Eqn (4) the RMSE is 57.1 m (95% of the mean GPR thickness; Fig. 6a).

To investigate the sensitivity of our ice thickness estimates to different input parameters, we analyzed the thickness and velocity variations along the central branch line of Tomur Glacier. We, therefore, measured surface velocities for the main trunk of Tomur Glacier using 2002/03 ASTER data as well and generated thickness maps using different DTM- and velocity-dataset combinations (see below). The mean elevation difference for the debris-covered part using DTM differencing is −36.4 and −23.3 m for the 1975–99 and 1999–2009 period which is higher than the average thickness change between ice thickness estimates of −3.5 m (SRTM + 2013/14 surface velocities – SRTM + 2002/03 surface velocities) and 0.9 m (SRTM + 2013/14 surface velocities – KH9 + 2013/14 surface velocities). The average slope for the debris-covered glacier tongue is slightly higher for 1999 (3.4°) compared with 1975 (3.1°). Decreasing surface velocities since the 1970s are, therefore, assumed to be the main reason for the discrepancy between multi-temporal DTM and ice thickness differencing for the 1975–99 period.

To approximate the sensitivity of the flow law parameter A, mean ice thickness for each glacier is estimated using the upper (2.4 * 10⁻²⁴ Pa⁻³ s⁻¹) and lower boundary (9.3 * 10⁻²⁵ Pa⁻³ s⁻¹) values of the flow law parameter. The variation of the mean ice thickness for Koxkar Glacier is between 76 and 96 m, for South Inylchek mean ice thickness varies between 114 m for temperate conditions and 144 m for cold ice conditions. The results show that variations in the flow factor have only an effect of ~9–15% on the estimated mean ice thickness.

For the basal velocity, Gantayat and others (2014) argue that ice thickness estimates are very less sensitive to changes in the contribution of the basal velocity to the surface velocity. Assuming basal velocity to be 90 or 85% of surface velocity would change the estimated ice thickness in the case of South Inylchek Glacier by more than 10% compared with <3% for Gangotri Glacier (Gantayat and others, 2014).

Relative mass loss

For the investigated glaciers in this study, we estimated a relative mass loss of ~6–28% since 1975, which is similar to the estimated total glacier mass loss of 27 ± 15% for the 1961–2012 period given by Farinotti and others (2015) for the entire Tien Shan; however, compared with the outer ranges the Central Tien Shan shows moderate loss only (Farinotti and others, 2015). To estimate the mass loss for the Central Tien Shan, we used the area and volume specifications from Avsiuk and Kotlyakov (1967) for the Khan Tengri-Tomur Peak region of 1517 km² and 379 km³ and the specific mass loss rate given by Farinotti and others (2015) of −0.21 ± 0.33 × 10³ kg m⁻² a⁻¹. The mass loss for the 1961–2012 period is then 18 km³ (density of ice 900 kg m⁻³) resulting in a mass loss of ~5 ± 7%. It has to be mentioned that the estimates of Avsiuk and Kotlyakov (1967) are based on the assumption that glaciers of the same type (hanging, valley, dendritic glaciers) have comparable mean ice thicknesses. A more detailed analysis would be required to prove and to understand the reason why glaciers in the Chinese part of the Central Tien Shan seem to be losing more mass than the regional average.

Overall, it is shown that glacier ice thickness can be estimated relying on GPR measurements for calibration, surface topography and glacier velocity information without additional data like observed or modeled mass balances. It is, therefore, promising for regions with limited field accessibility or uncertain input data. However, simplified model assumptions and uncertainties in the input parameters and datasets, particularly in the flow rate factor, the underlying DTM and DTM derivatives, and the velocity data lead to thickness uncertainties of 30–50% for the debris-covered parts of South Inylchek and Koxkar Glacier.

Approximately 60% of ice volume of a dendritic glacier system is concentrated on the glacier tongue, with the highest amount of volume around the median glacier elevation which is comparable with Linsbauer and others (2012) – for the Swiss Alps – and Frey and others (2014) – for the Himalayan-Karakoram region – who found the highest amount of ice stored around the glacier mean elevation and significant amount of ice in lower regions. The glacier ablation regions have been particularly influenced by glacier mass loss during the last decades (Pieczonka and Bolch, 2015). For the Tomur and the Koxkar Glacier, for instance, the surface elevation change rate between 1975 and 1999 in the debris-covered parts was negative and twice as high as for the entire glacier. These regions also exhibit the highest sensitivity regarding future climate change. Several authors reported rising temperatures in the Central Tien Shan during the last decades (Giese and others, 2007; Shangguan and others, 2009). Recent predictions from the IPCC AR5 suggest a MAAT increase of 1.0–1.3°C in Central Asia for the period 2016–35 in comparison with the period 1986–2005 (IPCC, 2013). The rise of mean winter temperature (December–February) is expected to be higher than for mean summer temperature (June–August). Hence, a significant increase of glacier downwasting can be expected in the long term causing glacier volume to decrease more rapidly than in the last decades.

CONCLUSIONS

In this study, we presented methods for glacier branch line and glacier ice thickness estimations for debris-covered glaciers of complex shape in the Jengish Chokusu/Khan Tengri massif (Central Tien Shan) using glacier outlines, glacier surface velocities and DTMs as a basic input.

The proposed method for branch line generation can handle DTM artifacts which may prevent the derivation of branch lines following strict hydrological criteria. We generated branch line networks considering all glacier branches and rock outcrops.

Glacier ice thickness has been calculated based on the equation of laminar flow with SRTM DTM and 2002/03, 2010/11 and 2013/14 glacier surface velocities as input. Under consideration of available field-based ice thickness measurements, we derived an empirical relation to approximate the contribution of the basal velocity to the glacier surface velocity. GPR measurements on Koxkar and South Inylchek Glacier differ by ~50 m on average, but the differences are lower compared with the well-known GlabTop approach.

The resulting maximum ice thicknesses vary between ~250 m for Koxkar Glacier and 380 m for South Inylchek Glacier. Uncertainties in ice thickness on debris-covered parts are between 30 and 50%. These uncertainties can be attributed to the uncertain contribution of basal sliding to the glacier surface velocity. Using geodetic glacier mass balances from existing studies, we estimated that the investigated glaciers lost between 6 and 28% of volume since 1975.

Because most of the ice volume is concentrated in regions below the equilibrium line, glaciers in the study region are highly sensitive to climate change. Due to rising temperatures glacier volume is likely to decrease more rapidly than in the last decades.

AUTHOR CONTRIBUTION

T.B. and T.P. designed the study and discussed the methodology and the results. T.P. generated glacier branch lines and estimated glacier ice thickness. M.K. calculated glacier velocities. J.P. contributed to the evaluation and validation of the results. L.S. provided validation data of Koxkar Glacier. T.P. wrote the draft of the article. All authors contributed to its final form.