Surface lowering of the debris-covered area of Kanchenjunga Glacier 1 in the eastern Nepal Himalaya since 1975 , as revealed by Hexagon 2 KH-9 and ALOS satellite observations 3 4

This study presents the geodetic mass balance of Kanchenjunga Glacier, one of the largest debris-covered glaciers 9 in the easternmost Nepal Himalaya, which possesses a negative mass balance of –0.18 ± 0.17 m w.e. a for the 1975–2010 10 study period, estimated using digital elevation models (DEMs) generated from Hexagon KH-9 and ALOS PRISM stereo 11 images. Accurate DEMs, with a relative uncertainty of ±5.5 m, were generated from the intensive and manual editing of 12 triangulated irregular network (TIN) models on a stereo MirrorTM/3D Monitor. The glacier ice-flow velocity field was also 13 calculated using a feature-tracking method that was applied to two ALOS orthoimages taken in 2010. The elevation 14 differences between the two DEMs highlight considerable surface lowering across the debris-covered area, and a slight 15 thickening in the accumulation area of Kanchenjunga Glacier between 1975 and 2010. The magnitude and gradient of surface 16 lowering are similar among the six glacier tributaries, even though they are situated at different elevations, which may reflect 17 variations in the ice-flow velocity field. The pattern of surface lowering correlates well with the ice-flow velocity field over 18 the debris-covered portion of the main tributary, suggesting that the glacier dynamics significantly affect surface lowering by 19 altering the emergence velocity along the glacier, particularly in the compressive ablation area. Surface lowering patterns 20 partially correspond to the supraglacial pond area fraction of the glacier, with enhanced surface lowering observed in areas 21 that possess a larger pond area fraction. These findings support the hypothesis that supraglacial ponds may intensify ice 22 wastage and play a key role in the heterogeneous surface lowering of debris-covered glaciers. The estimated mass loss of 23 Kanchenjunga Glacier is moderate compared with other debris-covered glaciers in neighboring Himalayan regions, which 24 may be due to the lower pond area fraction of Kanchenjunga Glacier relative to other glaciers. 25


Introduction
Glaciers in the Himalayas have been undergoing extensive and ongoing mass loss in recent decades, with the extent of this loss exhibiting a high degree of spatial heterogeneity (Fujita and Nuimura, 2011;Bolch et al., 2012;Kääb et al., 2012;Gardelle et al., 2013).There has been a growing interest in studying these changing Himalayan glaciers, but notable spatial and temporal gaps exist in existing data sets.Glaciers in the Khumbu (e.g., Nakawo et al., 1999;Bolch et al., 2011;Nuimura et al., 2011Nuimura et al., , 2012;;Salerno et al., 2015) and Langtang (e.g., Pellicciotti et al., 2015;Ragettli et al., 2016;Nuimura et al., 2017) regions of Nepal have been frequently studied compared with other regions in the Himalayas.Khumbu Glacier is one of the most extensively studied glaciers in the Himalayas, largely due to its proximity to Mt. Everest, which means that better logistical facilities provide easier access to the glacier.While these repeated investigations in a particular region or on a particular glacier help to strengthen our understanding of cryospheric processes, records from data-scarce regions are equally important.Glaciers in the Kanchenjunga region have received little attention and are thus poorly understood, although Basnett et al. (2013) and Racoviteanu et al. (2015) provided some constraints on glacier surface area changes.Declassified US spy satellite data (e.g., Corona KH-4 and Hexagon KH-9 stereo images) are now available, spanning the 1960s to mid-1980s time frame for many glacierized areas of the globe, which can provide valuable information about remote regions.Digital elevation models (DEMs) can thus be generated from these newly accessible data sources, allowing us to investigate the multidecadal mass balance of glaciers (e.g., Pieczonka et al., 2013;Pieczonka and Bolch, 2015;Maurer et al., 2016;Ragettli et al., 2016).
Many Himalayan glaciers are characterized by supraglacial debris cover in their ablation areas.These debris-covered glaciers possess more complex and variable responses to climate change than debris-free glaciers because the debris mantle can either insulate the ice or accelerate ice melting, depending on its thickness (e.g., Mattson et al., 1993;Mihalcea et al., 2008;Foster et al., 2012;Juen et al., 2014).This inhomogeneous distribution of debris across debris-covered glaciers can promote different rates of ice melting (e.g., Han et al., 2010;Zhang et al., 2011;Fujita and Sakai, 2014).Furthermore, the presence of supraglacial ponds and ice cliffs can lead to significant and intense ice wastage (e.g., Sakai et al., 2000Sakai et al., , 2002;;Steiner et al., 2015;Miles et al., 2016).However, detailed observations of ice cliffs and supraglacial ponds on debris-covered glaciers are often limited to small areas of the glacier surface to elucidate the importance of these small-scale features, since they require high-resolution observations.This is further complicated by the fact that many ice cliffs and supraglacial ponds on large and dynamic debris-covered glaciers are either physically inaccessible or too hazardous for conducting direct measurements.This study aims to provide new details on the evolution of Himalayan debris-covered glaciers by employing satellite imagery to estimate the recent geodetic mass change of Kanchenjunga  • N, 88.0-88.2• E; Fig. 1), a large and heavily debris-covered glacier in a data-scarce region in the easternmost Nepal Himalaya.This study also aims to provide details of the effects of supraglacial ponds and ice-flow velocity on the ongoing change of debris-covered glaciers.
2 Data and methods

Data
Two sets of optical stereo images, obtained by Hexagon KH-9 and Advanced Land Observing Satellite-Panchromatic Remote-sensing Instrument for Stereo Mapping (ALOS PRISM), are used in this study.The Hexagon KH-9 images have a spatial resolution of 6-9 m and a wide geographic coverage (125 km × 250 km), with consecutive ground nadir images possessing ∼ 70 % overlap (Surazakov and Aizen, 2010;Pieczonka et al., 2013).Three Hexagon KH-9 images, taken in December 1975 (Table 1), were obtained from the Center for Earth Resources Observation and Science of the U.S. Geological Survey, to create a DEM of Kanchenjunga Glacier and map its surface features.Three pairs of ALOS PRISM images (2.5 m spatial resolution) and rational polynomial coefficient (RPC) data were also acquired from the Remote Sensing Technology Center of Japan (Table 1), with the 2010 images analyzed to generate a DEM of Kanchenjunga Glacier, obtain surface ice-flow velocity measurements, and map its surface features.

Glacier delineation and hypsometry
The 1975 and 2010 glacier outlines were manually delineated from the orthorectified 3-D stereoviews, which were then used to generate and edit the glacier DEMs.The identification and delineation of the glacier boundary were feasible where the images were cloud-free and possessed minimal snow cover.However, areas with poor contrast, shade, and steep snow-covered slopes, which were generally associated with topographic features (i.e., slopes and contour lines) and geomorphic features (i.e., surface roughness and crevasses), were carefully checked for proper glacier surface delineation based on our experience with remote sensing analysis in high-mountain Asia (Nuimura et al., 2015  a For pond delineation for Khumbu Glacier.b For hypsometry.Nagai et al., 2016;Ojha et al., 2016Ojha et al., , 2017)).The uncertainty associated with the glacier surface delineation was estimated to be ±1 pixel along its perimeter (Fujita et al., 2009;Ojha et al., 2016).The glacier was divided into 50 m elevation bands to calculate the area-weighted average and thus volume change between the two DEMs.Because the 1975 and 2010 DEMs do not cover the entire glacier, ASTER GDEM2 (Tachikawa et al., 2011) was used to calculate the glacier hypsometry, as well as the elevation of boundary between the debris-covered and debris-free glacier surfaces for the six tributaries of Kanchenjunga Glacier (labeled T1 to T6 in Fig. 1).

DEM generation from ALOS PRISM imagery
The ALOS PRISM images were processed with their RPC data, which contain information about the interior (e.g., focal length and principle point of the camera/sensor) and exterior orientations (e.g., position and tilt of the camera/sensor) of the acquired images.The joint use of stereo images and RPCs makes geometric modeling feasible, thus removing the need for ground control points (GCPs) when generating DEMs and orthoimages.The ALOS stereo models and triangulated irregular network (TIN) model were produced from the RPC data and the Leica Photogrammetry Suite (LPS) workstation.However, automatic terrain extraction of TIN models often contains many errors, especially in areas of highly irregular and abrupt changes in topography, which typically consist of high relief, shadowed, and low-contrast regions in the images, leading to an inaccurate terrain rep-resentation (Lamsal et al., 2011;Sawagaki et al., 2012).It is thus necessary to edit the mass points, which are the sets of vertices in XYZ space that define the vector-based terrain surface, to obtain an accurate terrain representation.A StereoMirrorTM/3D monitor and Leica 3D TopoMouse were used for the GCP collection and terrain editing.The number of mass points defining the terrain surface depends on several factors, including the irregularity or uniformity and size of the feature, and the desired accuracy.Here we employed an average mass point density of 0.31 points per 100 m 2 (∼ 70 000 points over 22.6 km 2 ), ranging from ∼ 0.05 points per 100 m 2 on the gentle slopes around the upper boundary between debris-covered and debris-free surfaces of the glacier to ∼ 0.8 points per 100 m 2 on the bumpy debriscovered surface of the glacier.The LPS Terrain Editor was used to edit the TIN model until the terrain representation using mass points was satisfactorily achieved, thus minimizing the errors in the DEMs.The major editing tasks included the removal of false spikes (mass points above the actual surface) and depressions (mass points below the actual surface).Then, adequate and representative mass points were placed exactly on the terrain and glacier surfaces, including supraglacial ponds and moraine ridges.The edited TIN model was then gridded into a DEM with a spatial resolution of 15 m to reduce the effect of the different levels of resolution in the TIN model (hereafter ALOS-DEM).Further details on the TIN editing and DEM creation procedures are described by Lamsal et al. (2011) and Sawagaki et al. (2012).The Hexagon KH-9 images contain distortions from both the development and duplication of the films, as well as their long-term storage.These image distortions were corrected with the aid of the crosshairs in the images to make the images suitable for DEM extraction (Surazakov and Aizen, 2010).Since RPCs are unavailable for the Hexagon KH-9 images, GCPs were collected from the edited ALOS stereomodel, with a focus on boulders, trail intersections, and sharp notches on moraines that were clearly visible in both the ALOS PRISM and Hexagon KH-9 images.A total of 21 GCPs were extracted from the ice-free terrain surrounding Kanchenjunga Glacier (Fig. 1), with 5 GCPs (GCPv) randomly selected to independently verify the quality of the aerial triangulation (GCPv in Fig. 1).Comparison of the TIN models created with the inclusion of all 21 GCPs and with the exclusion of the 5 GCPv yielded vertical and horizontal root-mean-square errors of 4.4 and 7.7 m, respectively, thus validating the use of the GCPs in creating the Hexagon KH-9 TIN models.The Hexagon KH-9 TIN editing and DEM generation process followed that employed to create the ALOS-DEM (Sect.2.3).The generated Hexagon KH-9 DEM (hereafter Hex-DEM) was resampled at a spatial resolution of 15 m, following the ALOS-DEM, for consistency in comparing the DEMs.

Glacier ice-flow velocity
We calculated the surface ice-flow velocity field for Kanchenjunga Glacier using a feature-tracking method (Heid and Kääb, 2012) and then analyzed this ice-flow velocity field in combination with our DEMs to investigate the effects of topography and ice dynamics on glacier elevation changes.Orthorectified pairs of the ALOS PRISM images acquired in March and December 2010 (9-month gap) were processed using the Co-Registration of Optically Sensed Images and Correlation (COSI-Corr) algorithm, which was chosen because of its proven application in deriving terrain displacements, including glacier ice-flow velocities, in mountainous regions (e.g., Leprince et al., 2007;Scherler et al., 2008).The COSI-Corr estimates the phase difference in the Fourier domain, and then computes the relative surface displacement between the initial (reference image) and final image (search image).The stereo images were first orthorectified and then co-registered to ensure that the corresponding pixels of the ice-free terrain in each image overlap exactly, which is required to initialize the matching process.Ice-flow velocities were computed using a correlation window of 64×64 pixels, which corresponds to a 160 m × 160 m area, a robustness iteration of 4, and a mask threshold of 0.9 for noise reduction, following Leprince et al. (2007).To ensure the quality of the velocity map, poor matching in the surface displacement field was removed by applying a correlation threshold of 0.6, which resulted in some voids in the shaded area but also ensured that a robust ice-flow velocity map was produced (Fig. 2b).Further details on the COSI-Corr algorithm and its application to ice-flow velocity computations are described by Leprince et al. (2007) and Scherler et al. (2008).Uncertainties in the ice-flow velocity measurements were evaluated through analysis of the displacement field for the surrounding ice-free terrain possessing surface slopes gentler than 25 • .

Mapping supraglacial ponds
Supraglacial ponds and ice cliffs can enhance melt on debriscovered glaciers by absorbing radiative heat and essentially producing hot spots along the glacier surface (e.g., Sakai et al., 2000Sakai et al., , 2002;;Steiner et al., 2015;Miles et al., 2016).Constraints on their spatial distribution and density are thus important in determining the role of these surface features in modifying the glacier surface.We delineated supraglacial ponds along the glacier surface during the DEM creation process (Sect.2.3 and 2.4).Supraglacial ponds appear as distinct flat terrain features in the stereo images, whereas it is difficult to distinguish ice cliff and debris-covered steep slopes in the panchromatic images; we thus focused our analysis solely on the distribution and evolution of supraglacial pond coverage from the satellite images.We further limited our analysis to pond sizes of > 0.001 km 2 (12 × 12 pixels in the ALOS images) to avoid misinterpretations due to topographic features.We also performed a pond delineation analysis for Khumbu Glacier using the ALOS PRISM images taken in October 2008 (Table 1) for comparison.We calculated the fraction of pixels containing ponds per 50 m elevation band, following Ragettli et al. (2016).We also computed the rate of elevation change at each of the mapped ponds on a pixel-by-pixel basis between the two DEMs, and constructed a polygon-based map of supraglacial ponds in 2010 to derive the average area of each pond.

Geodetic mass balance and uncertainty estimates
The upper accumulation area of Kanchenjunga Glacier possesses extensive snow cover, which corresponds to high brightness and poor contrast regions in both the ALOS PRISM and Hexagon KH-9 images.These features precluded the creation of comprehensive DEMs for the entire glacier.However, patchy DEMs were generated in the upper accumulation area from nine locations that possessed better local image contrasts (hereafter named point sites, Fig. 2a).Two cases are employed to estimate the rate of elevation change in the upper accumulation area above 6100 m a.s.l.
where DEMs were unavailable: Case 1 employs the average of the rate of elevation change derived from the nine point sites (+0.01 m a −1 ) and Case 2 uses the bestfit curve for the debris-free part above 5800 m a.s.l., where positive elevation changes are observed between 5800 and 6100 m a.s.l.(Fig. 3a).Case 1 corresponds to the "no change assumption", with voids in the DEMs given values of zero  ( Pieczonka et al., 2013;Maurer et al., 2016), whereas Case 2 provides an alternative interpolation approach, with voids replaced by the regional mean of the corresponding elevation band (Gardelle et al., 2013).However, since there are no constraints on the rate of elevation change above 6100 m a.s.l., we placed an upper bound of 0.3 m a −1 to the curve fitting in Case 2, which is based on the results of similar studies that have analyzed elevation changes for glaciers in the eastern Himalayas (Nuimura et al., 2012;Maurer et al., 2016;King et al., 2017).The elevation change of the unmeasured area below 6100 m a.s.l. is assumed to be the same as the measured elevation change on the debris-free part.These two cases were thus applied to the upper 14.3 km 2 (19 % of the entire area) of Kanchenjunga Glacier to provide estimates on the surface elevation changes in the accumulation area.The surface elevation change calculated from the two DEMs (dh/dt, m a −1 ) was then converted into the geodetic mass balance of the entire glacier (b g , m w.e. a −1 ) using the following equation: where z is the elevation (m a.s.l.), ρ is the density of ice or firn, A z and A T are the areas at a given elevation (50 m elevation band) and of the entire glacier (km 2 ), respectively, and ρ w is the density of water (1000 kg m −3 ).The accumulation area of mountain glaciers often consists of snow and ice, whereas the ablation area is largely ice, which means that the elevation change in the ablation area relates more to mass change than that in the accumulation area.Numerous mass balance studies have employed a uniform ice density of ∼ 900 kg m −3 for entire glaciers (e.g., Bolch et al., 2011;Nuimura et al., 2012;Gardelle et al., 2013), whereas different density assumptions for the ablation and accumulation areas have been applied in recent studies (e.g., Kääb et al., 2012;Pellicciotti et al., 2015).We considered two density scenarios in this study: Scenario 1 assumes 850 ± 60 kg m −3 (Huss, 2013) for the entire glacier area, whereas Scenario 2 assumes 900 and 600 kg m −3 for the ablation and accumulation areas, respectively.An equilibrium line altitude (ELA) is required for analysis of the accumulation and ablation areas in Scenario 2. We assumed an ELA of 5850 m a.s.l. for the study period since the debris-free surface of Kanchenjunga Glacier possesses a positive elevation change (Fig. 3a).The mean of the four potential elevation changes (Cases 1 and 2) and density (Scenarios 1 and 2) combinations is assumed to represent the most plausible geodetic mass balance of the glacier (Kääb et al., 2012).The geodetic mass balance of the entire glacier was estimated as follows: the measured rate of elevation change at each 50 m elevation band was used for the lower ablation area (Fig. 3a), with the two assumptions of elevation change (constant value of +0.01 m a −1 for Case 1 and the best-fit curve in Fig. 3a for Case 2) and the two density scenarios adopted for the upper accumulation area above 6100 m a.s.l.(24 % of the entire area) where there were no available measurements.The averaged profile of elevation change for the debris-free glacier surface was then applied to rest of the unmeasured area between 5250 and 6100 m a.s.l.
The accuracy of the geodetic mass balance (σ g ) is a function of the two data sets that shape our mass balance calculation and their respective uncertainties, as follows: where ρ i is the ice density of 900 kg m −3 , σ a is the relative vertical accuracy between the two DEMs, and σ z is the difference in the assumed elevation changes of the unmeasured upper accumulation area (Cases 1 and 2).The relative vertical accuracy is evaluated as the standard deviation of the elevation difference (5.5 m and thus 0.16 m a −1 ) of the DEMs for the ice-free terrain (red polygon in Fig. 2a) and then uniformly applied to the glacier area below 6100 m a.s.l., which possesses approximately continuous data coverage (Fig. 3a).
We believe that the TIN editing performed on the Hexagon KH-9 and ALOS PRISM images guarantees the same degree of uncertainty over the measured area, and thus assume that it can be applied to the unmeasured area at the same elevation below 6100 m a.s.l.Gardelle et al. (2013) proposed employing the errors in the glacier elevation change within a given elevation band as the extrapolation error.However, the bright and poorly contrasted snowfields observed in our images hamper the estimation of this extrapolation error.We thus defined the error in the elevation change of the upper accumulation area (σ z ) as the difference in the elevation change between Case 1 (+0.01 m a −1 ) and Case 2 (best-fit curve in Fig. 3a).

Results
The surface area of Kanchenjunga Glacier decreased from 60.5 ± 1.6 km 2 in 1975 to 59.1 ± 0.5 km 2 in 2010, revealing a 1.4 ± 0.1 km 2 (0.070 ± 0.006 % a −1 ) area loss over the 35year study period.The average of the 1975 and 2010 areas (59.8 ± 1.1 km 2 ) was used to estimate the mass change of the glacier.No terminus retreat was noticeable; however, two minor tributaries in the upper catchments, which were connected to the major tributaries in 1975, retreated and were disconnected by 2010 (T1 and T6; a close-up of the T6 change is shown in Fig. 1), leading to very small deceases in glacier size (0.15 and 0.33 km 2 ).Considering the uncertainty in the area delineation (±0.5-1.5 km 2 ) and the 35year measurement interval, our estimated surface area loss of Kanchenjunga Glacier is negligible (2.3 ± 0.2 %).
The spatial distribution of the elevation change derived from Hex-DEM and ALOS-DEM is shown in Fig. 2a.DEMs are only available for the lower 22.6 km 2 section (38 % of the entire glacier); the debris-covered and debris-free areas are 15.0 and 7.6 km 2 in size, respectively.The most pronounced surface lowering was observed between 4700 and 5500 m a.s.l.(Fig. 3a).The general observations are a slight lowering of the glacier near the terminus (−0.4 to 0.0 m a −1 ), significant lowering across the main ablation area (−0.7 to −1.2 m a −1 ), and slight thickening in the uppermost debrisfree areas (0.0 to +0.4 m a −1 ), although considerable spatial variability is present across the glacier.
Figure 2b shows the spatial distribution of the ice-flow velocity of Kanchenjunga Glacier derived from the two ALOS PRISM orthoimages acquired in March and December 2010.The spatial distribution of displacement on the off-glacier area, which we assume to be the uncertainty in the ice-flow velocity, and its histogram, suggests that the uncertainty in the ice-flow velocity field is ∼ 2.7 m a −1 .Ice-flow velocities are almost negligible in the lowermost areas, moderate in the ablation areas, and fastest in the mid-glacier areas, varying from 0 to 72 m a −1 .However, distinct ice-flow velocities are observed between the six tributaries.Furthermore, the iceflow velocities for the main tributary (T3) show a variable pattern of speed-up within the central debris-covered area (> 10 m a −1 ) that is surrounded by significantly slower iceflow velocities on either side (< 10 m a −1 ).
Both the number and average size of supraglacial ponds across Kanchenjunga Glacier increased from 1975 to 2010, resulting in a marked increase in pond area fraction (ratio of the pond to debris-covered areas) from 0.1 to 1.1 % (Table 2).Most of the ponds formed in the ablation areas between 4700 and 5400 m a.s.l., coincident with the most significant lowering of the debris-covered surface (Fig. 3a).We observed a significant negative correlation between the size of the 2010 supraglacial ponds and the rate of elevation change, with the larger supraglacial ponds being associated with greater surface lowering (n = 35, r = 0.45, p < 0.01; Fig. 4).Geodetic mass balance of Kanchenjunga Glacier ranged from −0.15 to −0.20 m w.e. a −1 for the 1975-2010 study period (−0.18 ± 0.17 m w.e. a −1 , Table 3).We tested the sensitivity of the geodetic mass balance based on our assumptions.The density scenario does not significantly affect the resulting mass balance in both cases (< 0.02 m w.e. a −1 , Table 3).The assumption of the upper limit in Case 2 (0.3 m a −1 ) alters the mass balance only by 0.01 m w.e. a −1 , even if the upper limit is changed by ±0.1 m a −1 (gray shading in Fig. 3a).Although the mass balance estimates are largely influenced by the scenario used (Cases 1 or 2), the difference in the final mass balance is only ∼ 0.04 m w.e. a −1 (Table 3), which is one-fourth of the mass balance uncertainty (±0.17 m w.e. a −1 ).Unmeasured accumulation area due to poorly contrasted bright snow surface is a common issue in recent, similar studies (e.g., Maurer et al., 2016;King et al., 2017).

Contrasting surface lowering among tributaries
Profiles of elevation change over the debris-covered and debris-free surfaces of Kanchenjunga Glacier suggest that the estimated surface lowering is largely independent of elevation (Fig. 3a), in agreement with other studies across highmountain Asia (e.g., Bolch et al., 2011;Gardelle et al., 2013;Pieczonka et al., 2013;Ragettli et al., 2016).However, distinct relationships between surface lowering and elevation emerge when analyzing the surface elevation change data for each tributary (Fig. 5a).Elevation gradients of surface low- ering range from 1.58 to 2.36 m a −1 km −1 , with these observations confined to > 5500 m a.s.l.along T3 (upper T3), as opposed to the bulk of the elevation range for each of the other tributaries.The boundary separating the debris-covered and debris-free areas is found at variable elevations among the tributaries, suggesting that T1, upper T3, T5, and T6 largely possess debris-free surfaces, whereas T2 and T4 possess debris-covered surfaces for approximately half of their respective measured sections (Fig. 5a).However, it does not appear that surface lowering is affected by debris cover.Even though T5 and T6 are at lower elevations than T1, T2, and upper T3, the magnitude and gradient of surface lowering is similar for all six tributaries.Since greater surface melt is expected along T5 and T6 due to their lower elevations, we suggest that the higher ice-flow velocities along T5 and T6 are compensating for the anticipated higher melt rates at lower elevations by increasing the emergence velocity (Fig. 5b).The convergence of ice fluxes along lower T3 also affects the increased ice-flow velocities, thus suppressing the surface lowering at ∼ 5000 m a.s.l.These observations and inferences are in agreement with recent studies in the Himalayas that reported comparable surface-lowering trends along both debris-covered and debris-free surfaces (e.g., Kääb et al., 2012;Nuimura et al., 2012;Gardelle et al., 2013).Furthermore, we observe that the debris-covered area exhibits a more negative trend (−0.51 m a −1 ) than the debrisfree ice at ∼ 5100-5850 m a.s.l.(−0.30m a −1 ), where both debris-covered and debris-free areas exist along Kanchenjunga Glacier (Fig. 3a).

Influence of ice-flow velocity and supraglacial ponds on surface lowering
We compare the surface lowering with ice-flow velocities and supraglacial pond area fraction at 50 m elevation bands, following Ragettli et al. (2016).Figure 6 highlights the strong correlation between the surface lowering of the debriscovered area of Kanchenjunga Glacier and ice-flow velocity.
A strong positive correlation is observed between 4750 and 5600 m a.s.l.(n = 17, r = 0.74, p < 0.001), indicating that slower ice-flow velocities lead to increased surface lowering.Since the glacier surface is covered with debris mantle, the degree of surface melt is not a function of elevation, as is the case for debris-free ice.This high correlation between surface lowering and ice-flow velocity thus suggests that the observed pattern of surface lowering is affected primarily by the glacier dynamics, which experiences reduced surface lowering by increasing the emergence velocity in this area of the ablation zone.However, from the glacier terminus to ∼ 4750 m a.s.l., the opposite trend is observed (n = 5, r = 0.88, p < 0.05), likely due to a thicker debris mantle preventing surface melt, since the slower ice-flow velocities do not supply the necessary ice flux to compensate for the degree of surface melt anticipated for debris-free ice.The supraglacial pond area fraction at 50 m elevation bands shows a discontinuous, but partial, effect on surface lowering (Fig. 6).Four of the five elevation bands that possess a pond area fraction that is notably greater than the surrounding elevation bands align with a similar marked increase in surface lowering (light blue bars in Fig. 6).
It should be noted that both the ice-flow velocity field and supraglacial pond coverage for Kanchenjunga Glacier were obtained from the 2010 ALOS PRISM images, while the surface lowering was estimated for the 1975-2010 study period.The observed surface lowering in the ablation area likely results from the imbalance between the surface negative mass balance and emergence velocity, which produces a vertical upward motion due to the compressive flow regime in this region (Cuffey and Paterson, 2010), such that there is a direct correlation between the ice flux gradient and surface elevation changes (Nuimura et al., 2011(Nuimura et al., , 2017;;Vincent et al., 2016).However, since the middle section of the debris-covered area, at ∼ 5000 m a.s.l., exhibits an increase in the ice-flow velocity field where three tributaries coalesce with the main tributary T3 (T4, T5, and T6, Fig. 2b), the surface lowering here is suppressed by the increased ice flux and thus increased emergence velocity (Sect.4.1).These spatial patterns of the ice-flow velocity field would not change drastically though glacier thinning should yield a localized slowdown in the ice-flow velocity field.Even though we do not fully understand the evolution of supraglacial ponds on Kanchenjunga Glacier, it is unlikely that these ponds have persisted in the same locations for more than 3 decades, which is supported by the observed increase in pond-related variables between 1975 and 2010 (Table 2).Miles et al. (2017) revealed that supraglacial ponds on Langtang Glacier tended to form in the same locations along the glacier, and concluded that the geometry of the glacier influenced this tendency.Given the distinct increase in supraglacial pond coverage along Kanchenjunga Glacier between 1975 and 2010, the increased surface lowering where there is a higher pond fraction (Fig. 4), the persistent formation of supraglacial ponds (Miles et al., 2017), and the enhanced melting effects of these ponds on debriscovered glaciers (Sakai et al., 2000;Miles et al., 2016), the observed surface lowering may have been accelerated by recent supraglacial pond formation, given the drastic increase in pond area fraction between 1975 and 2010.Continued observations of the surficial evolution of Kanchenjunga Glacier could provide an indication of whether our observations are shaped by such a recent acceleration in supraglacial pond formation.
While there are a growing number of studies that focus on surface elevation changes along Himalayan glaciers (e.g., Maurer et al., 2016;Ragettli et al., 2016;Bolch et al., 2017;Brun et al., 2017;King et al., 2017;Lin et al., 2017), few studies have attempted to infer a direct link between surface lowering, ice-flow velocity, and supraglacial ponds and cliffs (Ragettli et al., 2016;Watson et al., 2017).We thus compare our results with those of Ragettli et al. (2016) for the glaciers in the Langtang region (Fig. 7).While the relationships between surface lowering, ice-flow velocity, and supraglacial pond fraction for the glaciers in the Langtang region are not as distinct as those observed for Kanchenjunga Glacier, weak correlations between surface lowering and ice-flow velocity are observed.Langshisha and Shalbachum glaciers show weak negative correlations between surface lowering and iceflow velocity (Fig. 7b and c).There is also some indication that greater supraglacial pond fractions may enhance surface lowering (5050 m a.s.l. at Langtang Glacier, 4700 m a.s.l. at Langshisha Glacier, 4500 m a.s.l. at Shalbachumu Glacier, and 4150 m a.s.l. at Lirung Glacier; light blue bars in Fig. 7), as found along Kanchenjunga Glacier, but these observations are not robust enough to draw any conclusive links between supraglacial ponds and surface lowering in the Langtang region.Watson et al. (2017) demonstrated a strong positive correlation between cliff density and recent surface lowering at nine debris-covered glaciers in the Everest region, with a > 50 % likelihood of supraglacial pond and cliffs coexisting along the glacier surface.Furthermore, Salerno et al. (2017) revealed that the surface lowering of 28 glaciers in the Khumbu region was due primarily to the surface slope gradient, which exhibited strong correlations with ice-flow velocity and secondary to proglacial ponds.

Comparison of geodetic mass balance with other regions
Here we compare the geodetic mass balance of Kanchenjunga Glacier (−0.18 ± 0.17 m w.e. a −1 ) with the results of www.the-cryosphere.net/11/2815/2017/The Cryosphere, 11, 2815-2827, 2017 other Himalayan glaciers that employed similar data sets (e.g., declassified Corona or Hexagon images and recent satellite images).Khumbu Glacier is the most intensively studied debris-covered glacier in the Himalayas, and is experiencing a greater degree of mass loss than Kanchenjunga Glacier, with a mass balance of −0.27 ± 0.08 m w.e. a −1 for the 1970-2007 period (Bolch et al., 2011), though these studies overlap within the uncertainty ranges.We suggest two possible causes for the mass balance differences between these two glaciers: (1) varying contributions of the accumulation and debris-covered areas and (2) varying density of supraglacial ponds.Here we use the hypsometry of Khumbu Glacier from Nuimura et al. (2012) for comparison.If we place the ELA for both glaciers at 5850 m a.s.l., as done for Kanchenjunga Glacier in our study, then a larger portion of Kanchenjunga Glacier (54 % of the entire glacier) lies above the ELA than observed on Khumbu Glacier (41 %), thereby highlighting that Kanchenjunga Glacier possesses a larger accumulation area that will suppress the negative mass balance of Kanchenjunga Glacier relative to Khumbu Glacier.The debris-covered area ratio of Kanchenjunga Glacier (25 %, 15.0 km 2 ) is also much smaller than that of Khumbu Glacier (40 %, 8.2 km 2 ).If the hypsometry of Kanchenjunga Glacier (Fig. 3b) is replaced by that of Khumbu Glacier (Nuimura et al., 2012), we will obtain more negative geodetic mass balance (−0.30± 0.17 m w.e. a −1 ), which is a similar value to the one observed for Khumbu Glacier (−0.27 ± 0.08 m w.e. a −1 ; Bolch et al., 2011).This combination of a larger accumulation area and smaller debris-covered fraction for Kanchenjunga Glacier would thus yield a less negative geodetic mass balance compared with Khumbu Glacier.Supraglacial ponds may also play a key role in the surface lowering of the debris-covered area.We delineated supraglacial ponds on Khumbu Glacier from ALOS PRISM images acquired in October 2008 for a direct comparison with our Kanchenjunga Glacier results (Table 2).While the average pond size was similar on the two glaciers, all other pond-related parameters are 2-4 times higher for Khumbu Glacier (Table 2).This high pond area fraction (3.7 %) likely contributes to the more negative geodetic mass balance of Khumbu Glacier.Recent investigations of the pond area fraction for glaciers in the Langtang (Ragettli et al., 2016) and Everest regions (Watson et al., 2017) yielded pond area fractions in the 2.3-3.3 % range for 4 debris-covered glaciers in the Langtang region (Ragettli et al., 2016) and in the 1-7 % range for 14 glaciers in the Khumbu region (Watson et al., 2017).We thus conclude that the moderately negative mass balance of Kanchenjunga Glacier is partially due to its smaller pond area fraction relative to other Himalayan glaciers.
We also compared the geodetic mass balances of debriscovered glaciers across the Himalayas (Fig. 8).The geodetic mass balances of smaller glaciers (< 20 km 2 ) possess a large degree of variability, while the four largest glaciers analyzed Figure 8. Geodetic mass balance of the debris-covered glaciers (solid circles) in the Himalayas, obtained from the difference between two DEMs generated from declassified satellite images and more recent satellite data.The red, blue, green, and black circles denote eight glaciers in the Khumbu region (Bolch et al., 2011), five glaciers in the Langtang region (Ragettli et al., 2016), five glaciers in Bhutan (Maurer et al., 2016), and Kanchenjunga Glacier (this study), respectively.Open squares with thick error bars are the areaweighted means and 1 standard deviation for the regional data sets.
(> 20 km 2 : two glaciers in Bhutan, one glacier in Langtang, and Kanchenjunga Glacier) exhibit a moderate mass loss (−0.15 to −0.24 m w.e. a −1 ; Fig. 8).Because of their large size, the area-weighted regional means of the geodetic mass balances yield values similar to those of the large glaciers (open squares in Fig. 8).Although we believe that the mass balance of Kanchenjunga Glacier could thus be viewed as representative of the region, more measurements should be accumulated for the regional mass balance in the easternmost Nepal Himalaya.

Conclusions
The geodetic mass balance of Kanchenjunga Glacier, one of the largest debris-covered glaciers in the easternmost Nepal Himalaya, is −0.18 ± 0.17 m w.e. a −1 for the 1975-2010 study period, as estimated from DEMs generated from Hexagon KH-9 (1975) and ALOS PRISM (2010) stereo images.While the TIN editing method employed in this study greatly improves the relative accuracy of the generated DEMs (5.5 m or ∼ 0.16 m a −1 ) and thus provides robust DEMs, the time-consuming manual editing process limited us to generating DEMs for only a single glacier.However, the DEMs in this study are valuable for validating the quality of regional DEMs generated from an automated method.Another shortcoming of this study is that the generated DEMs cover only 40 % of the entire glacier, leading to a large uncertainty in the estimated geodetic mass balance due the high The Cryosphere, 11, 2815-2827, 2017 brightness and poor contrast observed across the upper accumulation area in both the ALOS PRISM and Hexagon KH-9 images.This issue of unmeasured accumulation area is pointed out in recent, relevant studies for Himalayan glaciers.Even with these increased uncertainties regarding the geodetic mass balance of Kanchenjunga glacier, our estimate appears to be representative of the mass balance for Himalayan glaciers.
We observed that both the ice-flow velocity field and the presence of supraglacial ponds influenced surface lowering across the debris-covered area of Kanchenjunga Glacier.Although the six tributaries are situated at different elevations, the magnitude and gradient of surface lowering are similar for each, likely due to the varying ice-flow velocity field across the tributaries.Furthermore, surface lowering along the main tributary (T3) is highly correlated with the ice-flow velocity field, suggesting that the ice flux along T3 influences the emergence velocity and thus the degree of surface lowering, with vertical uplift in the ablation area countering some of the expected surface lowering along the tributary.Surface lowering generally increases along sections of the glacier that possess a larger supraglacial pond area fraction, supporting the hypothesis that supraglacial ponds contribute to localized intense ice wastage and play a key role in the heterogeneous surface lowering of debris-covered glaciers.While supraglacial ponds seem to accelerate the surface lowering of Kanchenjunga Glacier, the entire mass loss of Kanchenjunga Glacier is moderate compared with the other debris-covered Himalayan glaciers, which may be due to the lower pond area fraction and larger accumulation area ratio of Kanchenjunga Glacier relative to other glaciers.Since similar observations from the Langtang region in Nepal have shown rather equivocal relationships between surface lowering, ice-flow velocity, and supraglacial ponds, further analyses are needed to better understand the mechanisms that influence Himalayan debriscovered glaciers.

Figure 1 .
Figure 1. Outline of Kanchenjunga Glacier (blue), with 21 ground control points (GCPs), shown on a Hexagon KH-9 image taken in 1975.Five of the 21 GCPs were used to validate the accuracy of the transformation (GCPv).The major tributaries of Kanchenjunga Glacier are labeled T1 to T6.The inset map shows the location of Kanchenjunga Glacier (KJ, black box), as well as other glaciers (open boxes) in the Langtang region (LT), Khumbu region (KB), and Bhutan (BT), which were compared to this study (Fig. 7).The inset panel provides an example of the change in the glacier outline between 1975 (blue) and 2010 (orange).

D.
Lamsal et al.: Surface lowering of Kanchenjunga Glacier 2.4 DEM generation from Hexagon KH-9 imagery

Figure 2 .
Figure 2. (a) Rate of elevation change for the 1975-2010 study period and (b) ice-flow velocity between March 2010 and December 2010 for Kanchenjunga Glacier.The elevation difference (a) and displacement (b) in the ice-free terrain were used to evaluate the uncertainties of the two DEMs and ice-flow velocity (inset graphs).The black crosses in (a) denote the point measurements of elevation change in the upper accumulation area.The vectors and points in (b) are depicted at a 200 m spatial interval for better visibility.The box, thick line, circle, and whiskers in the inset graphs denote the interquartile, median, mean, and 1 standard deviation, respectively.

Figure 3 .
Figure 3. Elevation profiles of the (a) rate of elevation change (dh/dt) and (b) hypsometry of Kanchenjunga Glacier at 50 m elevation bands.The bars represent 1 standard deviation within each respective band.The green crosses and line denote the point measurements in the upper accumulation area (Fig. 2a) and the area distribution of the unmeasured part, respectively.The black line with gray shading is the best-fit curve for estimation of the elevation change in the unmeasured area (Case 2).

Figure 4 .
Figure 4. Rate of elevation change (dh/dt) of the supraglacial ponds observed on Kanchenjunga Glacier in 2010.

Figure 5 .
Figure 5. Elevation profiles of (a) the rate of elevation change (dh/dt) and (b) ice-flow velocity along the six tributaries of Kanchenjunga Glacier versus elevation.The tributaries are defined in Fig. 1.The elevation of the debris boundary of each tributary is indicated in (a).

Figure 6 .
Figure 6.Elevation profiles of the change in the rate of elevation change (dh/dt, black), ice-flow velocity (orange), and supraglacial pond area fraction (blue) at 50 m elevation bands over the debriscovered area of Kanchenjunga Glacier.The light blue bars denote the peaks in supraglacial pond area fraction mentioned in the main text.

Figure 7 .
Figure 7. Elevation profiles of the rate of elevation change (dh/dt, black), ice-flow velocity (orange), and supraglacial pond area (blue) versus elevation (in 50 m elevation bands) for the debris-covered area of four glaciers in the Langtang region, Nepal Himalaya, modified from Ragettli et al. (2016).The light blue bars in each panel denote peaks in the supraglacial pond area fraction that are mentioned in the main text.

Table 1 .
Details of remote sensing data used in this study to generate the digital elevation models (DEMs), map debris-covered and debris-free glacier surfaces and supraglacial ponds, and compute the surface ice-flow velocity field and hypsometry.

Table 2 .
Statistics of supraglacial ponds (≥ 1000 m 2 ) in the debris-covered areas of the Kanchenjunga and Khumbu glaciers.
n: number of ponds.A p : total area of ponds.a p : average pond area.A d : debris-covered area.R p : pond area fraction.D p : pond density.

Table 3 .
Area-weighted geodetic mass balance of Kanchenjunga Glacier for the 1975-2010 study period (m w.e. a −1 ).Two cases for the rate of elevation change, where DEMs are unavailable, and two density scenarios are assumed in our study.Scenario 1 assumes a constant density of 850 ± 60 kg m −3 for the entire glacier, whereas Scenario 2 assumes densities of 900 and 600 kg m −3 for the ablation and accumulation areas of the glacier, respectively.