Heterogeneous glacier thinning patterns over the last 40 years in LangtangHimal, Nepal

Abstract. This study presents volume and mass changes of seven (five partially debris-covered, two debris-free) glaciers in the upper Langtang catchment in Nepal. We use a digital elevation model (DEM) from 1974 stereo Hexagon satellite data and seven DEMs derived from 2006–2015 stereo or tri-stereo satellite imagery (e.g., SPOT6/7). The availability of multiple independent DEM differences allows the identification of a robust signal and narrowing down of the uncertainty about recent volume changes. The volume changes calculated over several multiyear periods between 2006 and 2015 consistently indicate that glacier thinning has accelerated with respect to the period 1974–2006. We calculate an ensemble-mean elevation change rate of –0.45 ± 0.18 m a−1 for 2006–2015, while for the period 1974–2006 we compute a rate of −0.24 ± 0.08 m a−1. However, the behavior of glaciers in the study area is heterogeneous, and the presence or absence of debris does not seem to be a good predictor for mass balance trends. Debris-covered tongues have nonlinear thinning profiles, and we show that recent accelerations in thinning correlate with the presence of supraglacial cliffs and lakes. At stagnating glacier areas near the glacier front, however, thinning rates decreased with time or remained constant. The April 2015 Nepal earthquake triggered large avalanches in the study catchment. Analysis of two post-earthquake DEMs revealed that the avalanche deposit volumes remaining 6 months after the earthquake are negligible in comparison to 2006–2015 elevation changes. However, the deposits compensate about 40 % the mass loss of debris-covered tongues of 1 average year.


Introduction
Atmospheric warming has caused widespread recent glacier thinning and retreat in the Himalayan region (Bolch et al., 2012).The impact of current and future glacier changes on Himalayan hydrology and downstream water supply strongly depends on the rate of such changes.However, planimetric and volumetric glacier changes are difficult to characterize due to limited data availability, and many recent studies have highlighted the spatially heterogeneous distribution of glacier wastage in the Himalayas (Fujita and Nuimura, 2011;Bolch et al., 2012;Kääb et al., 2012).Prominent examples of current-day regional differences in glacier evolution across the Hindu Kush-Karakoram-Himalaya (HKH) are the reported positive glacier mass balances in the Pamir and Karakoram.Glaciers in the rest of the HKH are thinning and receding (Bolch et al., 2012;Kääb et al., 2012;Gardelle et al., 2013).Across regions, differences in recent glacier evolution can often be associated with differences in climatic regimes (Fujita, 2008), particularly the varying influence of the south Asian monsoon and westerly disturbances (Yao et al., 2012).However, within the same climatic region the rate of glacier changes can also be heterogeneous (Scherler et al., 2011b).A main focus of current research is on the effect of supraglacial debris cover on glacier response to climate.Thick debris cover is a common feature in the HKH (Scherler et al., 2011b;Racoviteanu et al., 2015) and a homogenous layer of thick debris effectively reduces melt rates of the underlying ice (e.g., Östrem, 1959;Mattson et al., 1993).However, the characterization of debris-covered glacier re-sponse to climate is complicated by the frequent occurrence of ice cliffs and supraglacial lakes.At exposed cliffs, melt rates are much higher compared to the ice covered by a thick debris mantle (Sakai et al., 1998(Sakai et al., , 2002;;Immerzeel et al., 2014a;Steiner et al., 2015;Buri et al., 2016), and also at supraglacial ponds energy absorption is several times larger than that at the surrounding debris-covered surface (Sakai et al., 2000;Miles et al., 2016a).Recent large-scale geodetic studies based on remote sensing have provided evidence that the present-day surface lowering rates of some debriscovered areas in the HKH might be similar to those of debrisfree areas even within the same altitudinal range (Kääb et al., 2012;Nuimura et al., 2012;Gardelle et al., 2013), and surmise this could be due to enhanced melt from exposed ice cliffs and supraglacial lakes.Several detailed modeling studies, however, have provided evidence for a melt reducing effect of debris at the glacier scale (e.g., Juen et al., 2014;Ragettli et al., 2015) and have shown how supraglacial debris prolongs the response of the glacier to warming (Banerjee and Shankar, 2013;Rowan et al., 2015).Discrepancies between the different conclusions may be associated with glacier samples that are not comparable or with model uncertainties (particularly regarding the representation of the effect of supraglacial cliffs and lakes on total melt).Models can also provide actual melt rates while geodetic studies only provide glacier thinning rates, which are affected by glacier emergence velocity.
Programs to monitor debris-covered glaciers have been initiated in the Karakorum (e.g., Mayer et al., 2006;Mihalcea et al., 2006Mihalcea et al., , 2008) ) and in the Central Himalaya (e.g., Pratap et al., 2015;Ragettli et al., 2015).However, due to the logistical and financial constraints, long-term mass balance measurements are basically inexistent in the HKH.To document changes in debris-covered glacier thinning over time, declassified high-resolution reconnaissance satellite data available from the 1960s and 1970s are an important source of information (Bolch et al., 2008(Bolch et al., , 2011;;Maurer and Rupper, 2015).However, a common problem of previous multitemporal geodetic studies is the relatively low statistical significance of detected changes in glacier thinning over time.The uncertainties in digital elevation models (DEMs) derived from optical data and mass change measurements by DEM differentiation in the HKH arise from the difficult conditions for photogrammetric elevation analysis (due to extreme topography, surfaces with low contrast like bright snow cover or cast shadows).Radar-derived DEMs provide more accurate results in snow-covered areas but have even higher problems in steep terrain due to the side looking geometry.Unknown penetration of the radar beam into snow and ice is another shortcoming of DEMs derived from radar data.Uncertainties in volume loss estimates are therefore usually higher than identified acceleration in glacier thinning (Nuimura et al., 2012).The uncertainties are especially high over short periods.For long periods with much larger absolute eleva-  1. Monsoon snowcover frequency is based on Landsat 1999 to 2013 land cover classifications (Miles et al., 2016b).The 1974 glacier area (dotted lines) is shown for the seven studied glaciers only.
tion changes, the effect of errors in the DEMs weighs less and uncertainties in glacier volume changes are lower.
This study presents volume and mass changes of seven (five partially debris-covered, two debris-free) glacier in the upper Langtang catchment in Nepal.The aim of this study is to determine changes in thinning rates with high confidence by considering multiple independent DEM differences for short periods.For this we use seven DEMs derived from 2006-2015 stereo or tri-stereo satellite imagery and one DEM obtained from 1974 stereo Hexagon satellite data.We obtain an ensemble of multiannual elevation changes that provides a range of plausible values for the period between October 2006 and October 2015.We then assess if the elevation changes between different overlapping periods between 2006 and 2015 show similar characteristics.If this is the case, the ensemble of results can be used to identify statistically significant changes in thinning rates with respect to the longer period .Three main research questions are then addressed.First, we assess if overall thinning of glaciers in the region has accelerated.Second, we determine if spatial thinning patterns have changed over time.To explain changes in thinning rates we derive a number of glacier surface properties and glacier surface velocities.Third, we evaluate if there are major differences between the response of debris-covered and debris-free glaciers in the sample.Finally, we also look at the cryospheric impact of the April 2015 Nepal earthquake (7.8 magnitude, epicenter approximately 80 km west of the Langtang Valley).The earthquake devastated large parts of the Langtang catchment by triggering large avalanches and landslides (Kargel et al., 2016;Lacroix, 2016).Two post-earthquake DEMs from May and October 2015 are used to quantify the impact of the avalanche events on the mass balance of the debris-covered glacier tongues and assess its significance in comparison to multiannual volume changes.

Study site
We analyze the seven largest glaciers in the Langtang Valley (Langtang, Langshisha, Shalbachum, Lirung, Ghanna, Yala, Kimoshung), located in the monsoon-dominated Central Himalaya in Nepal, approximately 50 km north of Kathmandu and 100 km west of the Everest region.While Yala and Kimoshung glaciers are debris-free glaciers, all other studied glaciers have tongues that are almost entirely covered by supraglacial debris (Fig. 1).Langtang Glacier is the largest glacier in the valley with an area of 46.5 km 2 in 2006 (Table 1) and a total length of approximately 18 km.The smallest glacier is Ghanna Glacier with an area of 1.4 km 2 .
Critical debris thicknesses leading to a reduction of melt rates are exceeded over most parts of the debris-covered glacier area (Ragettli et al., 2015).Relatively thin debris appears only at the transition zone between accumulation and ablation area.At Lirung, Shalbachum, Ghanna and Langshisha glaciers the upper margins of debris-covered sections are located at the foot of steep cirques and icefalls, and transition zones are therefore very short.Ice cliffs and supraglacial ponds increase the heterogeneity of glacier surface characteristics in the Langtang Valley (Pellicciotti et al., 2015).
The ablation season of glaciers in the Langtang Valley lasts from April to September.The monsoon season (mid June-September) is at the same time the warmest and the wettest period of the year.Snow cover at the lower elevation of debris-covered glaciers is common only in winter (December-March).However, outside the monsoon period precipitation is limited and winters are rather dry (Collier and Immerzeel, 2015).

Satellite imagery
Multitemporal high-resolution data from different sensors are applied to assess glacier change in the upper Langtang catchment.Each type of remote-sensing data employed to calculate glacier elevation changes is listed below.Spatial and radiometric resolutions and base to height (b / h) ratios are provided in Table 2.
-Cartosat-1 is a remote-sensing satellite built by the Indian Space Research Organisation (Tiwari et al., 2008).We purchased radiometrically corrected alongtrack stereo imagery (processed at level "ortho-kit") of the upper Langtang catchment from October 2006 and November 2009.Cartosat-1 data have been previously used for DEM generation, e.g., in the Khumbu region in the Nepal Himalaya by Bolch et al. (2011) and Pieczonka et al. (2011).
We purchased a radiometrically calibrated along-track triplet mode scene from December 2010.
-SPOT6/7 (Système pour l'Observation de la Terre) along-track tri-stereo images were acquired upon request in April 2014 and May and October 2015.SPOT6 and 7 are the newest satellites of the SPOT series which have been frequently used for geodetic glacier mass balance studies (e.g., Berthier et al., 2007Berthier et al., , 2014;;Pieczonka et al., 2013).We acquired stereoscopic images in panchromatic mode corrected for radiometric and sensor distortions.Two of the three SPOT6/7 scenes used in this study were acquired in April/May, which means that limited amounts of winter snow are still present on the images.However, the imagery has a high spatial resolution (1.5 m) and high radiometric depth of 12 bit, which leads to good correlation results also over snowy parts.
-Overlapping pairs of high-resolution images acquired by the WorldView-2 and 3 satellites in February 2014 provide the basis of 8 m DEMs downloaded from http: //www.pgc.umn.edu/elevation(Noh and Howat, 2015).

DEM generation
The Hexagon DEM used here was generated for the study by Pellicciotti et al. (2015).We therefore refer to this study for further technical details regarding the Hexagon DEM.The SPOT6/7, Cartosat-1 and ALOS PRISM DEMs were generated for this study using the OrthoEngine module of PCI Geomatica 2015.We used the same parameters for DEM generation as proposed by Berthier et al. (2014) except setting the parameter "DEM detail" to "very high" instead of "low", which provided better results for the rugged debriscovered glacier surfaces., 2016).Because glacier motion and ablation have to be accounted for when using on-glacier dGPS points, we first generated a DEM from an across-track Pléiades stereo image pair from 1 and 9 November 2014 using the available dGPS points as ground-control points (GCPs).Glacier melt between 23 October and the acquisition dates of the Pléiades scenes is negligible due to the low temperatures during this period.The horizontal shift due to glacier motion during this period is less than the grid size of the Pléiades image (0.5 m) and is therefore also negligible.Subsequently, we determined 17 off-glacier GCPs on the basis of the Pléiades scene, which were then used to derive a DEM from the SPOT6 April 2014 tri-stereo scene.The Pléiades DEM itself is not used in the following to calculate glacier elevation changes since it covers only a small part of the catchment and since only low stereo matching scores were achieved at elevations higher than 4300 m a.s.l.due to snowfall onset between 1 and 9 November 2014.To guarantee high-quality GCPs, only pixels with correlation scores higher than 0.7 were considered for GCPs.Since the Pléiades scene covers only about one-fourth of the upper Langtang catchment, an additional 60 GCPs were determined on the basis of the April 2014 SPOT6 scene for the DEM extraction from the Cartosat-1, ALOS Prism and SPOT7 scenes.In addition to the GCPs, approximately 100 tie points for each scene were used to match stereo pairs before DEM extraction.The WorldView DEMs are 8 m posting DEMs produced using the Surface Extraction with TIN-based Search-space Minimization (SETSM) by Noh and Howat (2015).The WorldView DEMs rely on the satellite positioning model to locate the surface in space.The scenes from February 2015 which provide the basis of the two WorldView DEMs used in this study were acquired only 20 days apart (Table 2) and are adjacent to each other.The WorldView-2 DEM covers the western part of the study catchment and the WorldView-3 DEM the eastern part.Those DEMs were merged for this study and in the following are referred to as one single DEM representative of February 2015.
In addition to the DEMs discussed above, the 2000 SRTM (Shuttle Radar Topography Mission) 1 arcsec Global DEM (30 m spatial resolution) was used to calculate slopes and accumulation area ratios (AARs) of glaciers (Table 1) and to define 50 m altitude bands.However, the SRTM DEM was not used for DEM differencing because of the uncertainty re-garding the penetration depth of the radar signal into snow and ice (Gardelle et al., 2013;Kääb et al., 2015;Pellicciotti et al., 2015).Only DEMs extracted from optical stereo imagery are therefore employed to calculate elevation changes in this study.

Co-registration and DEM differencing
Co-registration of DEM pairs is applied in order to minimize the errors associated with shifts.Systematic errors in the elevation change maps due to tectonic uplift which could be relevant after the April 2015 Nepal earthquake are also corrected with the co-registration.For this purpose we exclude from each DEM the nonstable terrain such as glaciers and in general all off-glacier area at elevations higher than 5400 m a.s.l.(which is the estimated equilibrium line altitude (ELA) in the upper Langtang catchment; Ragettli et al., 2015).The correlation score maps, indicating which pixels have been matched successfully during the DEM extraction process, are used to exclude all DEM grid cells with a correlation score below 0.5.Then, horizontal shifts are determined by minimizing the aspect-dependent bias of elevation differences (Nuth and Kääb, 2011) between each DEM pair.Because of the slope dependency of the method all terrain below a slope of 10 • is excluded.The "older" DEM is then resampled (bilinear interpolation) according to the determined horizontal shift.In a second step the vertical DEM shifts and possible tilts are corrected using second-order trend surfaces fitted to all gently inclined (≤ 15 • ) stable terrain (Bolch et al., 2008;Pieczonka et al., 2011;Pieczonka and Bolch, 2015).
We resample all DEMs bilinearly to the grid size of the coarsest DEM (30 m) to reduce the effect of different resolutions.Elevation differences are calculated by subtracting the older from the younger DEM (such that glacier thickening values are positive) and are converted to elevation change rates by dividing by the number of ablation seasons between the acquisition dates.Seasonal effects on elevation change rates are neglected when discussing time intervals between DEMs of 4 years or longer, since elevation changes during the winter half-year are usually minor.On average, less than 20 % of annual precipitation occur during post-monsoon and winter (Immerzeel et al., 2014b) and less than 3 % of annual glacier ice melt (Ragettli et al., 2015).Area-average glacier elevation change rates are calculated using always the maximum glacier extent between two acquisition dates.

Processing of elevation change maps
Processing of the elevation change ( h / t) maps involves two main steps: (i) removal of pixel values identified as outliers and (ii) filling of gaps.

Outlier removal
The stereo matching score maps provided by PCI Geomatica are used to identify elevation data that can be considered for elevation change calculations.If the correlation score of a given DEM pixel is below 0.5, this indicates a poor matching score (Pieczonka et al. 2011) and therefore the corresponding h / t values are treated as "no data".Very unrealistic elevation change data (exceeding ± 150 m) are also excluded from the analysis.
We use the standard deviation (σ ) of observed elevation changes to identify h / t outliers.Outliers are defined separately for debris-covered glacier areas and debris-free glacier areas.For the latter we additionally distinguish between glacier area below and above the ELA (estimated at 5400 m a.s.l., see above).σ levels are thus calculated for each of the three area types in every h / t map.Below the ELA (both debris-free and debris-covered area), pixels are defined as outliers if h / t values differ from the average by > 3σ (e.g., Gardelle et al., 2013).This means that only very few data are classified as outliers, since 3 standard deviations account for 99.7 % of the sample (assuming the distribution is normal).The conservative outlier definitions are justified by the shallow slopes and high contrast, which also explains why stereo matching scores are generally higher below the ELA (Fig. 2c).Above the ELA, steep terrain or featureless snow surfaces lead to low DEM accuracy and therefore the outlier criteria should be more restrictive (e.g., Pieczonka et al., 2013;Pieczonka and Bolch, 2015).On debris-free glacier area above the ELA, pixels are therefore defined as outliers if h / t values differ from the average by > 1σ (which applies to approximately 32 % of the values if the distribution is normal).A stricter criterion for the accumulation area is also justified by the fact that it can be assumed that elevation changes in the accumulation areas over periods of several years are on average small (Schwitter and Raymond, 1993;Huss et al., 2010).Because we use different σ thresholds above and below the ELA we test the sensitivity of calculated glacier volume changes to a ±100 m ELA uncertainty.Furthermore, we test the sensitivity to different outlier definitions by comparing our results to the results obtained with a 2σ level applied to all area types.

Gap filling
On the glacier areas below the ELA, with only very few data gaps, missing data are replaced using inverse distance weighting.In the accumulation areas, in contrast, data gaps can extend over a wide elevation range if the terrain is steep or if the gaps are very large.Because of the elevation dependency of h / t values (e.g., Huss et al., 2010) only values from the same altitudinal range should be used to fill data gaps.We thus replace missing data in the accumulation areas by median h / t values per 50 m elevation band considering all available data for a given glacier (also from h / t maps representative of different periods).For this, we first calculate the mean elevation change rates per 50 m elevation band of each glacier and every h / t map and then determine the median of the ensemble.h / t maps that are re- jected from the ensemble (see Sect. 3.2.5 below) and in general all values representative of short periods ( t < 4 years) are not considered to calculate the ensemble-median values.

Uncertainty
Elevation change uncertainty estimates are based on the standard error E h calculated per elevation band (Gardelle et al., 2013).The standard error quantifies the effect of random errors on uncertainty according to the standard principles of error propagation: σ h, noglac is the standard deviation of the mean elevation change of non-glacierized terrain per elevation band, and N eff is the effective and N tot the total number of observations.PS is the pixel size (30 m) and d is the distance of spatial autocorrelation.d is equal to the range of the spherical semivariogram obtained by least-squares fit to the experimental, isotropic variogram of all off-glacier elevation differences (Wang and Kääb, 2015;Magnússon et al., 2016).The distance of spatial autocorrelation of the elevation change maps varies between 260 and 730 m with an average of 495 m.
To quantify the elevation change uncertainty of glacier area spanning several elevation bands, weighted averages of E h are calculated.E h of each individual elevation band is weighted by the glacier hypsometry.Elevation change uncertainties therefore vary for each individual glacier because of the different glacier area-elevation distributions.E h tends to increase with altitude (Figs. 3 and 4) due to steeper slopes, snow and deep shadows, which are factors that decrease the accuracy of DEMs derived from stereo data (e.g., Nuimura et al., 2011).Uncertainty estimates for each individual glacier therefore account for the spatially nonuniform distribution of uncertainty.Elevation change uncertainties of glaciers with a high accumulation area such as Kimoshung and Lirung glaciers (Table 1) are 50-100 % higher than those of other glaciers, in accordance with lower DEM matching scores (Fig. 2).The low uncertainty associated with debris-covered areas agrees with the 30-100 % lower off-glacier errors on shallow slopes (s < 18 • , 95th percentile of debris-covered glacier slopes) than on steeper slopes (s < 45 • , 95th percentile of glacier slopes; Fig. S1 in the Supplement).
The standard error can be interpreted as the 68 % confidence interval of the sample mean if the distribution is normal.Since we are conservatively assuming no error compensation across elevation bands the approximate confidence level in our uncertainty estimates per glacier is higher than 68 %.
This study aims at obtaining an ensemble of results about elevation change rates from the set of seven DEMs available for the period 2006-2015 and we thus calculate an ensem- For overall mass budget uncertainties we assume an ice density of 850 kg m −3 to convert the volume change into mass balance (Sapiano et al., 1998;Huss, 2013) and consider the elevation change rate uncertainties and an ice density uncertainty of 60 kg m −3 .

Ensemble selection
The possible two-fold combinations of all available DEMs are classified in two groups: maps that involve the Hexagon 1974 DEM and maps that represent only 21st century elevation changes (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).From the first group we only use the 1974-2006 h / t map, to strictly separate our two main study periods in 1974-2006 and 2006-2015.The redundancy of information between 2006 and 2015 allows extracting 12 maps of h / t with a time interval larger than 4 years.Since h / t uncertainties increase with shorter time intervals between DEMs (Fig. 5, Table 3) h / t maps with t < 4 years are not considered for the 2006-2015 ensemble.After careful evaluation of the DEMs in terms of h / t uncertainties and the 2015 Nepal earthquake impact we discard from the 2006-2015 ensemble also all h / t maps involving the ALOS PRISM 2010 and the SPOT7 May 2015 scenes.The 2006-2015 ensemble consists therefore of six maps of h / t (Fig. S2 in the Supplement).
The ALOS PRISM DEM is discarded because uncertainties associated with h / t maps involving this DEM are 30-100 % higher than if other DEMs are involved (Table 3).The ALOS-PRISM sensor has a radiometric resolution of 8 bit, which means that in comparison to a 12 bit image (SPOT6/7, Table 2), 2 4 = 16 times less information is provided per panchromatic image pixel.The image contrast is therefore lower, which decreases the accuracy of this DEM.
The May 2015 DEM is not considered for the 2006-2015 ensemble because of massive deposits of avalanched snow and ice as a consequence of the April 2015 Nepal earthquake, which strongly limits the representativeness of this DEM for the 2006-2015 period.The October 2015 SPOT7 DEM is still considered for the 2006-2015 ensemble because 6 months after the earthquake most avalanche material disappeared from the glacier (Sect.4.6).
Due to the incomplete representation of Langtang Glacier on the SPOT6 Apr 2014 scene (the scene does not cover the area north of 28 • 19 N; Fig. S2a and d in the Supplement), h / t maps involving this DEM are excluded when discussing ensemble results for Langtang Glacier.

Delineation of glaciers, debris-covered areas and supraglacial cliffs/lakes
The glacier outlines were manually delineated.We used the orthorectified satellite images with the least snow cover (the Cartosat-1 2006 and 2009 scenes) to delineate the accumulation areas and assumed no changes in the accumulation area over time.The tongues of the seven studied glaciers and debris extents were re-delineated for every year for which satellite images are available (1974, 2006, 2009, 2010, 2014 and 2015), using the corresponding orthorectified satellite images.A first operator delineated the outlines and a second operator provided feedback in order to improve delineation accuracy.To quantify the uncertainty in derived glacier area changes we consider a 0.5 pixel size delineation uncertainty (Paul et al., 2013).
The four largest glaciers in the valley were already delineated manually by Pellicciotti et al. (2015)    and 2000.However, we decided not to use those outlines because of the considerably higher resolution of the images that are available for this study and for consistency in the procedure applied for different outlines.We also re-delineated the catchment boundaries using the SRTM 30 m DEM and flow accumulation to accurately identify the ice divides between neighboring catchments.As a result, the calculated glacier areas (Table 1) changed considerably with respect to Pellicciotti et al. (2015).The 1974 glacier area of Langshisha Glacier changed by −40.4 % (Fig. S3 in the Supplement), mostly due to clipping with the catchment mask which reduced the extent of the accumulation areas.The 1974 areas of Langtang, Shalbachum and Lirung changed by −8.7, −9.5 and +8.0 %, respectively.
To identify glacier area associated with small glaciers in the catchment that are not discussed in this study we used the glacier outlines provided by the GAMDAM glacier inventory (Nuimura et al., 2015).Those areas were masked out from off-glacier terrain for the co-registration of the DEMs and stable terrain accuracy assessments.
Six quality-checked maps of supraglacial cliffs and lakes are used to characterize debris-covered glacier surfaces (Steiner et al., 2016).The cliff and lake inventories were generated based on the available satellite imagery for the Table 5. Sensitivity to outlier correction and equilibrium line altitude (ELA) definitions.2σ is the difference in results if a 2σ level is used to define outliers at all area types instead of a 3σ level above and a 1σ level below the ELA.The estimated ELA is 5400 m a.s.l.(Sugiyama et al., 2013;Ragettli et al., 2015).ELA ± 100 m represents the differences in results obtained with an ELA at 5500 m a.s.l. in comparison to results obtained with an ELA at 5300 m a.s.l.Sensitivity values that exceed uncertainty ranges as indicated in Table 4   period 2006-2015(October 2006, November 2009, December 2010, April 2014, May 2015and October 2015).As for the glacier outlines, cliff and lake outlines have been delineated by two independent operators.To further improve the accuracy of the inventories, a third operator used slope and elevation change maps to identify potential cliff and lake locations.The first two operators then used these indications to review the inventories.All outlines have been obtained by manual delineation on the basis of the orthorectified satellite images.
We calculated the fraction of pixels including lakes and cliffs per 50 m elevation band of each debris-covered tongue (excluding tributary branches, Fig. 6).In the following, we only discuss median 2006-2015 cliff and lake area fractions to minimize seasonal effects.Large avalanche cones, such as those present on Lirung and Langtang glaciers after the April 2015 earthquake, are masked out from the inventories before calculating median values.

Surface velocities
To assist with the interpretation of volumetric changes, we use glacier velocities determined with the COSI-Corr cross-correlation feature-tracking algorithm (Leprince et  The selected orthorectified images (5 m resolution) were adjusted according to the shifts determined by co-registration (Sect.3.2.2).Since the window size must be large enough to avoid correlating only noise but small enough not to degrade the output resolution (Dehecq et al., 2015), we tested several configurations.The best results for the COSI-Corr correlation analysis were achieved using multiscale window sizes of 128 down to 32 pixels, as also proposed by Scherler et al. (2008).To post-process the velocity data we removed pixels with x or y velocity values greater than 40 m a −1 , since these were identified as errors by manually measuring the surface displacement on the basis of the orthorectified images and prominent features.We then ran a median filter on the data to remove areas which show a local reversal in x or y directions.Missing values were then filled with the mean of the adjacent eight values.Finally, the velocity map was resampled to 30 m resolution with a bicubic algorithm.

Assessment of the April 2015 earthquake impact
We quantify the impact of the avalanche events after the April 2015 earthquake on volume changes of debris-covered tongues.For this purpose we use the April 2014-May 2015 h map to quantify the accumulated volumes less than 2 weeks after the earthquake and the April 2014-October 2015 h map to quantify the remaining volumes after one ablation season.To identify glacier area where avalanche material accumulated we consider all glacier grid cells with significant positive elevation changes ( h > 5 m).Approximately 7.9 % (1.9 km 2 ) of all debris-covered areas were affected by avalanches according to this definition.To account for preearthquake volume losses we first subtract from the DEM differences the elevation changes between April 2014 and April 2015 determined on the basis of the October 2006-February 2015 mean thinning rates.Note that we do not use the February 2015-May 2015 and the February 2015-October 2015 h maps to quantify avalanche debris volumes because the calculated uncertainties associated with these maps are up to 300 % higher than the uncertainties associated with differential DEMs involving the April 2014 scene (Table S1 in the Supplement).

Mean glacier surface elevation changes
The 2006-2015 ensemble consistently indicates an increase in mean glacier thinning rates in comparison to the period 1974-2006 (Fig. 7h).For 2006-2015 we calculate an ensemble-mean thinning rate of −0.45 ± 0.18 m a −1 , while for the period 1974-2006 we identify a thinning rate of −0.24 ± 0.08 m a −1 (Table 4).This corresponds to an increase in determined mean thinning rates by 0.21 m a −1 or 87.5 %.The error bounds associated with the two periods are overlapping at the extremes.However, error bounds are not overlapping at 80 % confidence levels (assuming normal distribution).Given the probability of less than 10 % for 1974-2006 and 2006-2015 thinning rates for being above or below this confidence interval, the estimated confidence level of accelerated thinning rates is higher than 99 %.
From the seven studied glaciers in the valley, the thinning rates of Langtang, Langshisha and Yala glaciers have accelerated at 99 % confidence levels (Fig. 7, Table 4).At Shalbachum Glacier the error bounds are overlapping but the estimated probability of thinning acceleration is more than 90 %.At Lirung and Kimoshung glaciers the mean elevation change rates have likely remained approximately constant (Table 4).Mean thinning rates of these glaciers increased by less than 0.10 m a −1 between 1974-2006 and 2006-2015.Also at Ghanna Glacier the 1974-2006 value and the 2006-2015 ensemble mean differ by only 0.05 m a −1 (Table 4).However, the scatter in the 2006-2015 values is such that no trend can be identified.Ghanna Glacier is the only glacier where the ensemble of values available for the period 2006-2015 did not narrow down the uncertainty associated with individual periods (Fig. 7).
Elevation change rates are also calculated separately for the five debris-covered tongues (Fig. 8, Table 4).An increase in thinning rates is evident on the Langtang, Langshisha, Shalbachum and Lirung tongues.Thinning rates increased between 15 % (Langtang tongue) and 68 % (Langshisha and Shalbachum tongues).Changes in thinning rates are not significant for Ghanna tongue, but five out of six members of the 2006-2015 ensemble decreased as compared to the previous period.
Of all debris-covered areas, the downwasting rates on Lirung tongue are the highest.This applies to both the period 1974-2006 (−1.03 ± 0.05 m a −1 , Table 4) and to the 2006-2015 ensemble mean (−1.67 ± 0.59 m a −1 , Table 4).The 2006-2015 ensemble uncertainty is very large on Lirung tongue (± 0.59 m a −1 ), which we believe is due to systematic errors in the 2009-2014 differential DEM that represents an outlier in the ensemble (Fig. 8).However, neither on Lirung nor on Langtang tongue (the two glaciers most   estimates are in most cases not significantly influenced by outlier definitions and ELA estimates.

Altitudinal distribution of elevation changes
The altitudinal distribution of mean elevation changes clearly show that the thinning patterns of all debris-covered tongues have changed over time (Figs. 9 and 10).Areas with clear increases in thinning rates can be identified for Langtang Glacier 5000-5150 m a.s.l.(25-100 % thinning rate increase), for Langshisha Glacier 4650-5100 m a.s.l.(25-260 %), for Shalbachum Glacier 4500-4800 m a.s.l.(25-180 %) and for Lirung Glacier 4300-4350 m a.s.l.(80-170 %).Thinning rates have remained mostly constant in the lower third of the elevation ranges of the tongues (Langtang, Shalbachum and Lirung glaciers).At Ghanna Glacier, thinning rates have recently declined near the glacier terminus at 4800-4850 m a.s.l.(60-90 % thinning rate decrease; Fig. 9e).This pattern contrasts with all other temporal patterns for debris-covered glacier areas.On Langshisha Glacier thinning patterns in 1974-2006 and 2006-2015 are different near the terminus (Fig. 9b).Here, the glacier tongue became very narrow in the last decade and ultimately a small part below 4500 m a.s.l.dis-connected from the main tongue (Fig. 1) between 2010 and 2014.The fragmentation of the tongue lead to mean thinning rates close to 0 at elevation bands where a substantial part of the glacier area disappeared.
Overall, the thinning profiles of 2006-2015 ensemble members show very similar characteristics (Figs. 9 and 10).The profiles diverge for the uppermost elevation bands of the tongues and in the accumulation areas.This agrees with the larger error that is attributed to higher elevations (Fig. 3).Above 5500 m a.s.l. it is impossible to separate uncertainty from actual differences in thinning rates.
On Yala Glacier there has been a 3-fold increase in thinning rates below 5400 m a.s.l, comparing 1974-2006 to the 2006-2015 ensemble results (Fig. 10d).Maximal thinning takes place at the terminus and then decreases nearly linearly with altitude until it reaches values close to 0. This is in clear contrast to the much less uniform patterns on debris-covered glaciers (Fig. 10a-c).On debris-covered glaciers, the elevation corresponding to the maximum thinning rates is different from glacier to glacier.On Shalbachum and Lirung glaciers the maximum is reached somewhere close to the upper end of the tongue (4650-4750 and 4300-4400 m a.s.l. as indicated in Fig. 9c and d, respectively), on Langtang and Ghanna glaciers more in the middle part (4950-5150 and 4900-5000 m a.s.l., indicated in, respectively, Fig. 9a and e) and on Langshisha Glacier closer to the terminus (4450-4700 m a.s.l., Fig. 9b).On the large debris-covered glaciers, areas of maximum thinning seem to have shifted and extended to higher elevations only at Langtang Glacier, where during the period 1974-2006 maximum thinning occurred between 4850 and 4950 m a.s.l.(Fig. 9a).On Shalbachum Glacier maximum thinning during the period 1974-2006 occurred slightly higher up at 4750-4800 m a.s.l.(Fig. 9c).
Note that the altitudinal h / t profiles (Figs. 9 and 10) always refer to the same position in space, since 50 m elevation bands were delimited only once on the basis of the SRTM 1 arcsec Global DEM.To account for the up-valley movement of on-glacier elevation bands over time due to surface lowering, profiles would have to be slightly shifted relative to each other.However, given the maximum thinning rates of 1-1.5 m a −1 in 1974-2006, the maximum relative adjustment of values in Figs. 9 and 10 would never exceed one 50 m elevation band.Accounting for the shifting of elevation bands over time would therefore not lead to different conclusions regarding changes in spatial h / t patterns.

Glacier area changes
Debris-free Yala Glacier experienced the strongest increase in relative annual area loss of all studied glaciers : −0.43 ± 0.05 % a −1 ; 2006-2015: −1.77 ± 0.16 % a −1 ; Table 6).During the same two time intervals debris-free Kimoshung Glacier shrank only at rates of 0.08 ± 0.01 and 0.05 ± 0.02 % a −1 , respectively.This represents significantly lower retreat rates for the second period than at Yala Glacier.The differences in area change rates are consistent with the identified differences in mean glacier surface elevation changes, where the two glaciers also represent opposite extremes (Sect.4.1).

Surface velocities and supraglacial cliff/lake areas
Approximately 10 % of all grid cells for the three largest debris-covered tongues (Langtang, Langshisha, Shalbachum) contain supraglacial cliff features ("cliff area" in Table 7).At Lirung and Ghanna tongues this value decreases to 8 and 3 %, respectively.For Ghanna tongue practically no supraglacial lakes could be identified, while at the other debris-covered tongues "lake area" is between 2.3 and 3.3 %.
The mean surface velocities of the tongues range between 1.6 m a −1 (Ghanna tongue) and 7 m a −1 (Langshisha tongue).The mean and the standard deviation of off-glacier surface velocities are 1.3 and 1.9 m a −1 , respectively.At Ghanna and Lirung tongue, which both have a mean surface velocity below 3 m a −1 , it is therefore practically impossible to discriminate moving ice from quasi-stagnant ice.Following Scherler et al. (2011b), all glacier grid cells with a surface velocity of less than 2.5 m a −1 are therefore termed "stagnant" for simplicity.According to this definition, the tongue area classified as "stagnant" (Table 7) ranges from 20 % (Langshisha tongue) to 85 % (Ghanna tongue).
In our sample of five debris-covered glaciers, cliffs and lakes seem to appear more frequently on glaciers which are dynamically active.We identify a highly significant negative correlation (Pearson's linear correlation coefficient r = −0.99) between cliff area fraction per tongue and the percentage of stagnant tongue area."Lake area" and "% stagnant area" are also negatively correlated (r = −0.87).At the scale of individual tongues, a correlation between surface velocities and cliff appearance is evident at Shalbachum Glacier (Fig. 11c), where we identify a correlation of 0.85 (respectively 0.68) between the altitudinal velocity profile and cliff (respectively lake) areas per 50 m elevation band.Also on the two other large debris-covered tongues in the valley, the Langtang and Langshisha tongues, cliff appearance clearly decreases towards the termini where the glaciers are quasistagnant.However, on these two glaciers the surface velocities and cliff appearance are not linearly correlated since the highest cliff area densities are identified 200-300 m below www.the-cryosphere.net/10/2075/2016/The Cryosphere, 10, 2075-2097, 2016 Cliff and lake area corresponds to the percentage of 30 m pixels containing cliffs/lakes (median of six available cliff and lake maps from the period 2006-2015).Mean velocity is calculated on the basis of 2009-2010 surface velocities (Fig. 13).To discriminate moving ice from quasi-stagnant ice we use a threshold of 2.5 m a −1 (cf.Scherler et al., 2011b), which also corresponds to the approximate uncertainty of remote-sensing derived surface velocity.
the altitude ranges corresponding to maximum surface velocity.
To investigate a possible link between accelerated thinning and the presence of supraglacial lakes and cliffs we compare "cliff area" and "lake area" (as provided in Table 7) to changes in mean thinning rates per tongue ( h / t, difference between "1974-2006" and "ensemble mean 2006-2015" as provided in Table 4).Overall, the correlation coefficient between fractional cliff area per tongue and h/ t is −0.62 (and −0.50 between lake area and h / t).The likely reduced thinning rates on Ghanna tongue (Fig. 8e) indeed correspond to low cliff and lake area fractions (3.2 and 0.4 %, respectively).On Lirung, Shalbachum and Langshisha tongues thinning accelerated by 0.47-0.64m a −1 , whereas fractional cliff and lake areas are similar (cliff area: 8.0-10.5 %; lake area: 2.3-2.6 %).Also Langtang tongue is characterized by relatively high cliff and lake area fractions (10 and 3.3 %, respectively; Table 7) but the identified changes in thinning rates are only minor.The acceleration of mean thinning rates at Langtang tongue is significant at the 95 % confidence level (Fig. 8a), but the difference in mean thinning rates 1974-2006 and 2006-2015 is only −0.12 m a −1 (Table 4).
At locations where thinning rates did not increase significantly we mostly identify low cliff area fractions below 10 % (e.g., at Langtang tongue below 4750 m a.s.l., and above 5150 m a.s.l., at Shalbachum below 4500 m a.s.l. or on Ghanna tongue).Conversely, cliff area fractions are generally higher than 10 % where the 2006-2015 ensemble consistently indicates thinning acceleration (Fig. 11).Exceptions to this observation are the high cliff area fractions at Langtang Glacier (4750-4900 m a.s.l.), where thinning rates did not change significantly (Fig. 11a), and low cliff area fractions at Shalbachum Glacier (4750-4800 m a.s.l.), where thinning rates increased (Fig. 11c).Lirung tongue also shows a different behavior, except at the lowest elevation band.However, maximum thinning acceleration at 4300 m a.s.l.corresponds to a relatively high lake area fraction of 6 % (Fig. 11d).
Altitude bands with no significant increases in thinning rates on Langtang Glacier consistently coincides with relatively low surface velocities below 5 m a −1 .At Langshisha and Shalbachum tongues this is also the case (Fig. 11).Across all debris-covered glacier tongues, 77 % of all elevation bands where thinning accelerated ( h / t < −0.2 m a −1 ) are not stagnating (velocities above 2.5 m a −1 ); in 72 % of all elevation bands where thinning rates remained constant or declined ( h / t ≥ −0.2 m a −1 ) we observe velocities below 2.5 m a −1 .

Impacts of the April 2015 earthquake
We calculate a total volume of post-earthquake avalanche debris in May 2015 of 2.5 × 10 7 m 3 , which is equivalent to a cube length of 292 m. 40 % of the avalanche material re-   8).The two glaciers which were most affected by avalanches were Langtang Glacier (receiving 58 % of the total volume) and Lirung Glacier (29 %).The avalanche cone at Lirung Glacier piled up to a height of nearly 60 m, while the avalanche mate-rial at Langtang Glacier was more spread (Fig. 12).Consequently, more material remained until 6 October 2015 at Lirung Glacier (57 %), while at Langtang Glacier 31 % remained (Table 8).Field visits at the end of October 2015 to Lirung and Langtang glaciers revealed that a smooth dewww.the-cryosphere.net/10/2075/2016/The Cryosphere, 10, 2075-2097, 2016  bris layer melted out of the avalanche material and covered the surface uniformly with a thickness of a few centimeters (P.Buri and P. Egli, personal communication, 2016).
Volume loss from glacier area where avalanche material accumulated between May and October 2015 was 30 times higher than during an average ablation season (May 2015-October 2015: −1.5 × 10 7 m 3 ; October 2006-April 2014: −5 × 10 5 m 3 a −1 ).After this rapid initial downwasting the avalanche deposits diminished to a volume of 10 7 m 3 , equivalent to an average positive surface elevation change over all debris-covered glacier area of 0.52 ± 0.19 m (Table 8, or 0.09 m a −1 if divided over 6 years, which is the shortest time interval of any h / t map in the ensemble involving the October 2015 DEM).The avalanche impact on the Octo-ber 2015 DEM is thus within the uncertainty range associated with multiannual h / t values (± 0.12 m a −1 , Table 3) and justifies why the October 2015 DEM is considered for the 2006-2015 ensemble.

Elevation changes of debris-covered glaciers
Elevation changes in the debris-covered area are primarily independent of elevation (Fig. 9), as previously identified in the Langtang catchment (Pellicciotti et al., 2015) and elsewhere in high-mountain Asia (e.g., Bolch et al., 2011;Dobhal et al., 2013;Pieczonka et al., 2013;Pieczonka and Bolch, 2015;Ye et al., 2015).Such patterns have usually been explained by down-glacier increase of debris thickness and by ablation associated with supraglacial lakes and exposed ice cliffs.Our analysis shows that, with few exceptions, the highest thinning rates and the strongest increase in thinning rates can be associated with areas with a high concentration of ice cliffs and supraglacial ponds (Fig. 11 and Fig. S4 in the Supplement).While previous studies have pointed out that debriscovered areas with a large presence of supraglacial cliffs and lakes make a disproportionately large contribution to ablation (Reid and Brock, 2014;Buri et al., 2016;Miles et al., 2016a;Thompson et al., 2016), this is the first study which documents the relation between accelerations in thinning rates and the large presence of supraglacial cliffs and lakes.
Supraglacial cliffs seem to appear more frequently on slowly moving ice (5-10 m a −1 , Fig. 11) and not where the glacier is stagnant (Sakai et al., 2002;Bolch et al., 2008;Thompson et al., 2016).This can be explained by compres-sive stresses associated with flow deceleration that may initiate fracturing (Benn et al., 2009).Such stresses are usually not large enough to initiate open surface crevasses, but in combination with elevated water pressure due to local water inputs they lead to hydrologically driven fracture propagation (hydrofracturing) and englacial conduit formation (Benn et al., 2009).The collapse of large englacial voids destabilizes the debris layers and leads to the formation of new ice cliffs.Accordingly, accelerated thinning of debris-covered area in the upper Langtang catchment does not take place on stagnating parts of the tongues but rather at the transition between the active and the stagnant ice (Fig. 11).
The appearance of supraglacial lakes, however, is strongly related to the surface gradient (Sakai and Fujita, 2010;Miles et al., 2016b).Large supraglacial lakes can only form where the slope is less than 2 • (Reynolds, 2000) and where local water input is high.These conditions are not met on debriscovered glacier sections in the upper Langtang catchment, since local surface slope is consistently above 5 • (Pellicciotti  et al., 2015).It is interesting to note that the highest lake area fractions (lake area > 6 %) are found on avalanche deposition zones at Langtang Glacier (4750-4800 m a.s.l.; Figs.7a and  11a) and at Lirung Glacier (4300 m a.s.l.; Figs.7d and 11d).This is likely related to high local surface water inputs from melting of avalanche snow and ice.On Langtang Glacier frequent avalanche inputs may explain why thinning did not accelerate at the altitude range between 4750 and 4900 m a.s.l., in spite of the presence of exposed ice (cliff area > 13 %; Fig. 11a).
Several studies show that lakes and cliffs are important but cannot explain the mass loss alone (e.g., Sakai et al., 2002;Juen et al., 2014).The high thinning magnitudes on the upper sections of Shalbachum tongue (4750-4800 m a.s.l.) likely cannot be attributed to lakes and cliffs (cliff/lake area fractions are below 5 %; Fig. S4c in the Supplement), and thin layers of deposited debris in the upper sections of the glacier tongue could explain such patterns.
Reduced ice fluxes also contribute to thinning accelerations.To assess how much this factor contributes to the observed accelerations in thinning it would be necessary to quantify changes in ice flux over time (e.g., Nuimura et al., 2011;Berthier and Vincent, 2012;Nuth et al., 2012).However, information about the evolution of surface velocities over long time periods would be required, which our dataset cannot provide.Over the timescales considered in this study, however, high warming rates have been identified in this part of the Himalaya (Shrestha et al., 1999;Lau et al., 2010).The rise in air temperatures directly impacts glacier melt rates, and can explain rapid acceleration of thinning where ice is not insulated from warming by thick debris.Banerjee and Shankar (2013) numerically investigated the response of extensively debris-covered glaciers to rising airtemperatures and describe the dynamical response as follows: during an initial period the fronts remain almost stationary and in the ablation region a slow-flowing quasi-stagnant tongue develops.During this period, which may last more than 100 years, glaciers lose volume by thinning.After this initial period glaciers start to retreat with a higher rate, while annual volume loss decreases because of thickening debris layers.Since thinning rates near the fronts of the large debris-covered glaciers in the valley (Langtang, Langshisha and Shalbachum glaciers) have not yet started to significantly decrease (Fig. 11a-c) and the glacier tongues are still dynamically active (Fig. 13) it can be assumed that the quasistationary length period will persist for these glaciers in the near future.The model of Banerjee and Shankar (2013) does not account for supraglacial cliffs and lakes, which likely contribute to thinning acceleration (Fig. 11).However, we have shown that they primarily appear on parts of the glacier tongues which are still dynamically active (Table 7).Our results suggest that they become less abundant with decreasing flow.The presence of cliffs and lakes therefore does not interfere with the dynamical response of debris-covered glaciers as described by Banerjee and Shankar (2013).
Near the snout of Ghanna Glacier a deceleration in thinning rates by −80 % can be clearly identified (Figs. 9e and 11e,.Previous studies have provided evidence that ablation rates of debris-covered ice may decrease over time as a consequence of thickening debris cover, in spite of rising air temperatures (Banerjee and Shankar, 2013;Rowan et al., 2015).

Elevation changes of debris-free glaciers
The downwasting rates of 2006-2015 on Yala Glacier are 0.5-1.2m a −1 higher than on Kimoshung Glacier (Table 4).However, the two glaciers have a very different hypsometry (Fig. S5 in the Supplement).Kimoshung Glacier has a very steep tongue that reaches to similarly low elevations as the debris-covered glacier tongues (Table 1).The glacier is nearly in equilibrium with the climate (Table 4), which explains the low thinning rates at low elevations (Fig. S5 in the Supplement).Currently the estimated AAR of Yala Glacier is 40 % (Table 1), which is a common value in the HKH region (Kääb et al., 2012).The estimated AAR of 86 % at Kimoshung Glacier, in contrast, corresponds to an exceptionally high value for the HKH (Khan et al., 2015).The differences in thinning rates point to the role of glacier hypsometry for the response of debris-free glaciers to climatic changes (e.g., Jiskoot et al., 2009).Almost balanced mass budgets in recent years (Table 4) and only minor area changes (Table 6) are associated with Kimoshung Glacier.Thinning did not increase significantly with respect to the period 1974-2006 (Fig. 7g).Due to the steep tongue of this glacier the AAR is not sensitive to changes in the ELA (Table 5).Only a small fraction of area is therefore additionally exposed to temperatures above freezing level in case of atmospheric warming, which causes the glacier to be less sensitive to observed warming trends in the region (Shrestha et al., 1999;Lau et al., 2010).One possible explanation for the balanced condiwww.the-cryosphere.net/10/2075/2016/The Cryosphere, 10, 2075-2097, 2016 tions of Kimoshung Glacier could therefore be that precipitation in recent decades remained approximately stable, which agrees with the findings of studies on precipitation trends in this part of the Himalaya (Shrestha et al., 2000;Immerzeel, 2008;Singh et al., 2008).However, further analysis is required for justification.The mass balance of Yala Glacier, in contrast, is sensitive to fluctuations in temperature.A hypothetical rise of the ELA by 100 m at this glacier causes 30 % of its area to turn from accumulation into ablation area (Table 5).Accordingly, thinning of Yala Glacier is accelerating rapidly (Fig. 9f).Due to the common AAR of Yala Glacier and the extreme topography of Kimoshung Glacier it can be assumed that other debris-free glaciers in the region are also thinning and that balanced conditions such as observed on Kimoshung Glacier are exceptional.

Differences between debris-free and debris-covered glaciers
The response of debris-covered and debris-free glaciers to warming is substantially different, as described in the two sections above and exemplified by the altitudinal elevation change profiles in Fig. 10.Our observations do not support the findings of previous studies about similar present-day lowering rates of debris-covered and debris-free glacier areas at the same elevation (Kääb et al., 2012;Nuimura et al., 2012;Gardelle et al., 2013).Also for debris-covered elevation bands where up to 18 % of the area is covered by supraglacial cliffs and lakes (e.g., at Langtang tongue 5050 m a.s.l. or at Langshisha tongue 4750 m a.s.l.) thinning rates do not exceed 1.8 m a −1 , while for Yala Glacier the lowering rates are already above this value at 5250 m a.s.l. and further increase down glacier (Fig. 10).Within the same altitudinal range (5200-5300 m a.s.l.) thinning rates of debriscovered glaciers do not exceed 35-75 % of the thinning rates of Yala Glacier.
Regarding the mean surface elevation changes (Table 4), our observations reveal a heterogeneous response to climate of both the debris-free and the debris-covered glaciers.As discussed in the two sections above, there are examples for both types of glaciers where thinning has increased significantly or where thinning remained approximately constant.A significant difference in thinning trends between debrisfree and debris-covered glaciers in our sample cannot be identified.In our sample, the best predictor for thinning accelerations seems to be the altitude distributions of glaciers.Glaciers with a high AAR (Kimoshung) or which reach the highest elevations (Lirung) have the most balanced mass budgets and show no significant changes in volume loss over time (Fig. 7, Table 4).Glaciers which are most sensitive to ELA changes (more than ± 10 % AAR change in response to ± 100 m ELA uncertainty, Table 5) such as Yala, Langtang and Langshisha glaciers reveal the most significant thinning accelerations (Fig. 7, Table 4).However, debris-free Yala Glacier is currently downwasting at 60-100 % higher rates than the large debris-covered glaciers in the valley.Considering the common characteristics of Yala Glacier and given that this glacier has been denominated as a benchmark glacier for the Nepal Himalayas (Fujita and Nuimura, 2011), it seems important that future geodetic or field-based studies extend our analysis to larger glacier samples.

Post-earthquake avalanche impacts
Accumulation by debris-laden avalanches is one of the most important processes for debris-covered glacier formation (Scherler et al., 2011a).The tongue of Lirung Glacier would likely not exist without accumulation through avalanches (Ragettli et al., 2015).It is detached from the accumulation area (Fig. 1) and reaches 200-700 m lower elevations than all other debris-covered glaciers (Table 1).Our volume calculations of the post-earthquake avalanche impact allow quantifying the avalanche impact on mass balance and comparing it to mass loss during an average year.Given the avalanche deposits remaining on Lirung tongue by 6 October 2015 (divided by the area of the tongue: 3.87 ± 0.23 m, Table 8) and the average h / t rates between October 2006 and February 2015 of −1.64 ± 0.10 m a −1 (Fig. 8d), the avalanche after the earthquake compensated by 240 % the volume loss of 1 average year.At the scale of all debris-covered area in the valley this value amounts to 50 % (0.52 ± 0.19 m avalanche deposits and −1.02 ± 0.08 m a −1 average thinning).According to Scally and Gardner (1989) avalanche deposit density increases until the end of the ablation season to about 80 % of ice density.The mass deposits therefore compensate mass loss during a normal year by about 180 % at Lirung tongue (40 % at the catchment scale).Still, our analysis has revealed that the impacts are not significant in comparison to the 2006-2015 ensemble uncertainty (Sect.4.6, Fig. 8d and  f).

Comparison to other studies
The four largest debris-covered glaciers in the valley (Langtang, Langshisha, Shalbachum, Lirung) have been the focus of a recent geodetic mass balance study by Pellicciotti et al. (2015), who reconstructed elevation and mass changes using the 1974 Hexagon DEM which is also used in this study (spatial resolution 30 m) and the 2000 SRTM3 DEM (90 m).They found that all four glaciers lost mass over the study period but with different rates (on average −0.32 ± 0.18 m w.e. a −1 ).We find an overall glacier mass balance for the period 1974-2006 of the four glaciers which is slightly less negative (−0.22 ± 0.08 m w.e. a −1 ).However, the results match within the uncertainties.The lower uncertainty estimates by our study are justified by the high resolution and quality of the 2006 Cartosat-1 DEM (Table 3).Differences in the mass balance of Langtang, Lirung and Shalbachum glaciers are within uncertainty bounds and can be attributed to differences in used glacier masks, study The Cryosphere, 10, 2075-2097, 2016 period, outlier correction approaches, density assumptions and uncertainties regarding the penetration depth of the SRTM radar signal (Kääb et al., 2015).However, for Langshisha Glacier we calculate a mass balance which is substantially less negative than in Pellicciotti et al. (2015).While we identify almost balanced conditions for the period 1974-2006 (−0.10 ± 0.08 m w.e. a −1 , Table 4), the mass balance indicated by Pellicciotti et al. (2015) is very negative (−0.79 ± 0.18 m w.e. a −1 ).The discrepancy can be explained by the overestimated extent of the accumulation areas by Pellicciotti et al. (2015) (Fig. S3 in the Supplement) in combination with unrealistic lowering rates of up to −2 m a −1 at about 6000 m a.s.l.(Fig. 4d in Pellicciotti et al., 2015).The more realistic elevation change values obtained by the present study for the accumulation areas (ranging between −0.4 and +0.4 m a −1 , Fig. 10b) point to the need of restrictive outlier definitions and the advantage of having information from multiple datasets available for gap filling.
Yala Glacier has been frequently visited for field measurements in the last 25 years.Sugiyama et al. (2013) calculated mean thinning rates of Yala Glacier for the periods 1982-1996 (−0.69 ± 0.25 m a −1 ) and 1996-2009 (−0.75 ± 0.24 m a −1 ) on the basis of ground photogrammetry and GPS surveys.The values suggest a more moderate acceleration of thinning rates than in our study (−0.33 ± 0.06 m a −1 1974-2006 to −0.89 ± 0.23 m a −1 2006-2015, Table 4).However, similarly to our study Sugiyama et al. (2013) identified a rapid acceleration of thinning rates at the lowest elevations.At higher elevations the uncertainty of photogrammetric surveys increases because of low contrast due to homogeneous snow layers.
The acceleration in mass loss in recent periods identified by this study agrees with other studies from the Nepalese Himalaya which assess multitemporal elevation changes (Bolch et al., 2011;Nuimura et al., 2012).Bolch et al. (2011) identify an increase in mass loss rates by 0.47 m w.e. a −1 comparing the two periods 1970-2007 (−0.32 ± 0.08 m w.e. a −1 ) and 2002-2007 (−0.79 ± 0.52 m w.e. a −1 ).Nuimura et al. (2012) calculate increasing mass losses in the same study region between 1992-2008 (−0.26 ± 0.24 m w.e. a −1 ) and 2000-2008 (−0.45 ± 0.60 m w.e. a −1 ).However, the identified acceleration in glacier thinning is not significant given the largely overlapping error bounds.Moreover, similar mass loss rates were calculated by Bolch et al. (2011) for 1970-2007and by Gardelle et al. (2013) for the same 10 glaciers and for the period 2001-2011 (average of −0.41 ± 0.21 m w.e. a −1 ).The ensemble approach of this study can therefore substantially strengthen previous conclusions that mass loss of glaciers in the Central Himalaya is accelerating.The volume changes calculated over several multiyear periods between 2006 and 2015 consistently indicate that glacier thinning has indeed accelerated (Fig. 7h).

Conclusions
This study presents glacier volume changes of seven glaciers (five partially debris-covered, two debris-free) in the upper Langtang catchment in Nepal, using a DEM from 1974 stereo Hexagon satellite data and seven DEMs derived from 2006-2015 stereo or tri-stereo satellite imagery.We carefully selected elevation change maps which are least affected by uncertainty to obtain multiple independent DEM differences for the period 2006-2015.
Our results point to increasing thinning rates, from −0.24 ± 0.08 m a −1 in 1974-2006 to −0.45 ± 0.18 m a −1 in 2006-2015, where the estimated confidence level of accelerated thinning rates is higher than 99 %.This study therefore supports the findings of previous studies (Bolch et al., 2011;Nuimura et al., 2012) that glacier wastage in the Central Himalaya is accelerating.However, whereas a majority of glaciers in the study region are thinning rapidly, glaciers with a high accumulation area have almost balanced mass budgets and experience no or only insignificant accelerations in thinning.
Our observations also reveal that thinning has mostly accelerated in the upper reaches of the tongues (up to +150 %, comparing the periods 1974-2006 and 2006-2015), while the nearly stagnant areas near the terminus show constant or decreasing thinning rates (up to −80 %).The highest thinning rates and the strongest increase in thinning rates can be associated with areas with a high concentration of ice cliffs and supraglacial ponds.Constant or decelerating thinning rates can be associated with areas with relatively homogeneous debris layers near the termini of glaciers.We conclude that the response of extensively debris-covered glaciers to global warming is largely determined by feedback processes associated with different surface characteristics.
The behavior of glaciers in the study area is highly heterogeneous, and the presence or absence of debris is not a good predictor for mass balance trends.However, the spatial thinning patterns on debris-covered glaciers are fundamentally different than those on debris-free glaciers.Debrisfree glaciers in our sample present thinning rates that are linearly dependent on elevation, while debris-covered glaciers have highly nonlinear altitudinal elevation change profiles.Our observations do not provide evidence for the existence of a so-called debris-cover anomaly, where the insulating effect of thick supraglacial debris is compensated by enhanced melt from exposed ice cliffs or due to high energy absorption at supraglacial ponds.Within the same altitudinal range, lowering rates on debris-free Yala Glacier are 35-300 % higher than on debris-covered glacier area.On debris-free Kimoshung Glacier the thinning rates are similar to those of debris-covered area, but this result must be explained by compressive flows that compensate for surface lowering by ablation due to its exceptionally balanced conditions.
Geodetic mass balance studies such as this have been increasingly revealing heterogeneous patterns of changes and S. Ragettli et al.: Heterogeneous glacier thinning patterns over the last 40 years in Langtang Himal, Nepal a complex response of debris-covered glaciers that call for an enhanced understanding of processes over debris-covered glaciers.Their ablation, mass balance and response to climate is modulated by debris supply, transport, glacier flow, lake and cliff developments and a complex subglacial hydrology and hydraulics that all need to be understood in the future to be able to predict future changes of these glaciers over multiple timescales.
7 Data availability DEM data, surface velocity data and glacier/cliff/lake inventories are available upon request.Please contact Francesca Pellicciotti (francesca.pellicciotti@northumbria.ac.uk) for this purpose.Certain licence restrictions may apply to data derived from commercial satellite imagery.
The Supplement related to this article is available online at doi:10.5194/tc-10-2075-2016-supplement.

Figure 1 .
Figure1.Map of the upper Langtang catchment.The numbers on the map correspond to the glaciers listed in Table1.Monsoon snowcover frequency is based on Landsat 1999 to 2013 land cover classifications(Miles et al., 2016b).The 1974 glacier area (dotted lines) is shown for the seven studied glaciers only.

Figure 2 .Figure 3 .
Figure 2. (a) Uncertainty estimates of average elevation change rates ( h / t) per individual glacier and per debris-covered tongue.(b) Ensemble of stereo matching scores per individual glacier and debris-covered tongue and (c) per glacier area above and below the estimated ELA.The central marks correspond to the median of all h / t maps or DEMs in the ensemble.The edges of each box are the 25th and 75th percentiles.The whiskers extend to the most extreme data points.
Figure 5. Uncertainties in elevation change rates ( h / t) as a function of the time interval between DEMs ( t).Median results of all available 28 h / t maps.Error bars extend to the most extreme data points.

Figure 6 .
Figure 6.Supraglacial cliffs and lakes as identified from the October 2006 Cartosat-1 satellite image: (a) Langtang and Ghanna glaciers, (b) Shalbachum Glacier, (c) Lirung Glacier and (d) Langshisha Glacier.Cliff area shows the median fraction (%) of 30 m pixels per 50 m elevation band that contain cliffs, considering all six available cliff maps from 2006 to 2015.

Figure 7 .Figure 8 .
Figure 7. Mean elevation change rates ( h / t) per period and glacier.For better readability, only the maximum width of error bounds corresponding to individual periods 2006-2015 is shown.Note that the 1974-2006 timescale is not linear (dashed dark blue line).

Figure 9 .Figure 10 .
Figure 9. Altitudinal distribution of mean annual elevation change ( h / t) over 50 m elevation bands of debris-covered tongues (debriscovered area of each glacier excluding tributary branches).Uncertainty bounds correspond to uncertainty as a function of elevation derived for each h / t map individually (Fig. 3).

Figure 11 .
Figure 11.Altitudinal distribution of cliff and lake area fractions, glacier velocity and changes in thinning rates.Cliff and lake area is shown as % of 30 m pixels containing cliffs/lakes per 50 m elevation band, whereas the values represent the median of six available cliff and lake maps from the period 2006-2015.Glacier velocities (m a −1 ) represent the median per 50 m elevation band of data shown in Fig. 13 and error bars represent the standard deviation in pixel values per elevation band.Changes in thinning rates ( ( h / t) [m a −1 ]) are calculated comparing 1974-2006 and the 2006-2015 ensemble mean.Negative ( h / t) values represent thinning accelerations.Error bars represent the maximum variations in ( h / t) considering all individual periods within the 2006-2015 ensemble.

Figure 12 .
Figure 12.Avalanche affected sections of Lirung and Langtang glacier, before and after the earthquake on 25 April 2015, and corresponding surface elevation changes ( h).Imagery © Airbus DS 2014/2015.

Figure 13 .
Figure 13.Surface velocities 2009-2010 cropped to catchment boundaries.Values have units of meters per year and are derived by cross-correlation feature tracking.Off-glacier velocities are shown in transparent color.

Table 1 .
Characteristics of the studied glaciers in the upper Langtang catchment.The measures are based on the SRTM 1 arcsec Global DEM and glacier outlines of 2006.

Table 3 .
Mean * uncertainties associated with different sets of elevation change ( h / t) maps.Uncertainties associated with individual maps are shown in TableS1in the Supplement.
are printed in bold letters.
al., www.the-cryosphere.net/10/2075/2016/The Cryosphere, 10, 2075-2097, 2016 2007) and the available satellite imagery.The orthorectified Cartosat-1 Nov 2009 and ALOS-PRISM December 2010 images were used for this purpose.Other image pairs were not considered due to longer periods between acquisitions (leading to image decorrelation) or the presence of snow patches at lower elevations (SPOT6 April 2014, SPOT7 May 2015).

Table 4
ences can be explained by unrealistic patterns (strongly negative elevation changes above 5400 m a.s.l.) that are not identified as outliers with a 3σ threshold applied to areas below 5500 m a.s.l.Our analysis thus shows that elevation change www.the-cryosphere.net/10/2075/2016/The Cryosphere, 10, 2075-2097, 2016

Table 7 .
Characteristics of the debris-covered tongues (debris-covered glacier area excluding tributary branches).

Table 8 .
Elevation changes of debris-covered glacier tongues due to avalanches triggered by the Nepal earthquake on 25 April 2015.The first three data columns provide the volume changes of avalanche-affected area ( h > 5 m in May 2015) divided by the total debris-covered area (Table1).Only lower part (south of 28 • 19 N), upper part not on April 2014 scene * Estimation based on average annual thinning October 2006-April 2014 * *