Moderate mass loss of Kanchenjunga Glacier in the eastern Nepal Himalaya since 1975 revealed by Hexagon KH-9 and ALOS satellite imageries

This study presents the geodetic mass balance of Kanchenjunga Glacier, a heavily debris-covered glacier in the easternmost Nepal Himalaya, between 1975 and 2010 using high-resolution (5 m) digital elevation models (DEMs) generated from Hexagon KH-9 and ALOS PRISM stereo images. The glacier velocities were calculated using a feature tracking method 10 with two ALOS orthoimages taken in 2010. The elevation difference between the two DEMs indicates considerable surface lowering across the debris-covered area, and slight thickening in the accumulation area between 1975 and 2010. The flow velocity throughout the debris-covered area is slow, in contrast with the faster flow velocity in the lower accumulation area. The rates of elevation change positively correlate with the elevation of the debris-free part, while they negatively correlate with the elevation of the debris-covered part, which might be due to the debris thickness distribution. The rate of elevation 15 change also positively correlates with the glacier flow velocity and longitudinal slope. Significant surface lowering was observed at supraglacial ponds, though the ponds should have short life spans. The geodetic mass balance of Kanchenjunga Glacier, which is the largest measured glacier in the neighbouring Himalayas, is estimated to be –0.18 ± 0.17 m w.e. a for the period of 1975–2010; it is less negative than that estimated for debris-covered glaciers.


20
In recent decades, glaciers in the Himalayas, once one of the mountainous areas in the world with the most glaciers, have been widely losing mass, though they exhibited a high degree of spatial heterogeneity Bolch et al., 2012;Kääb et al., 2012;Gardelle et al., 2013). Studies on the change of Himalayan glaciers changes have increased significantly in recent years; however, there are still notable gaps, spatially and temporally. For instance, glaciers in the Khumbu and Langtang regions of Nepal have been frequently studied (e.g., Nakawo et al., 1999;Bolch et al., 2011;Nuimura et al., 2011Nuimura et al., , 25 2012Pellicciotti et al., 2015;Salerno et al., 2015;Ragettli et al., 2016). In particular, Khumbu Glacier is one of the most extensively studied glaciers in the Himalayas, which might be due to easier access to the glacier, better logistic facilities, and the existence of Mt. Everest. While repeated investigations in a particular region or of a particular glacier help to deepen our understanding of cryospheric processes, records from data-scarce regions are also important. Glaciers in the Kanchenjunga region (present study site) have been studied little and are thus poorly understood, though Racoviteanu et al. (2015) 30 investigated glacier surface area changes between the 1960s and 2000s. In-situ monitoring programs often have shortcomings pertaining to their limited temporal and spatial coverage; modern satellite observations are mostly limited to the year after 2000, despite their large spatial coverage. Declassified US spy satellite data (e.g., Corona KH-4 and Hexagon KH-9 stereo images) are available between the 1960s and the mid-1980s for many glacierized areas of the globe, which provide valuable information about remote regions. Digital elevation models (DEMs) can be generated from both of these data sources and allow us to investigate the multidecadal mass change of glaciers (e.g., Pieczonka et al., 2013;Maurer et al., 2016;Ragettli et al., 2016).
Many glaciers in the Himalayas are characterized by the presence of supraglacial debris cover in their ablation areas (i.e., 5 debris-covered glaciers); their responses to climate change are more complex than those of debris-free glaciers because of the insulation effect of the debris mantle (e.g., Mihalcea et al., 2008;Foster et al., 2012), inhomogeneous distribution of debris thickness (Han et al., 2010;Zhang et al., 2012), and presence of supraglacial ponds and ice cliffs (Sakai et al., 2000(Sakai et al., , 2002. Earlier studies have reported that the variability in the debris thickness resulted in different rates of ice melting. In particular, ice cliffs and supraglacial ponds are primarily responsible for significant contributions to intense ice wastage (Sakai et al., 10 2000;Röhl, 2008;Han et al., 2010;Steiner et al., 2015). However, those observations were not performed on the spatial scale of an entire glacier but were limited to several ice cliffs and supraglacial ponds, which are small-scale landform features of the glacier surface, because they require high-resolution (e.g., in-situ) observations. Many ice cliffs and supraglacial ponds of large and dynamic debris-covered glaciers, however, are either physically inaccessible or too hazardous to conduct direct measurements.

15
To fill the spatial gap of the decadal change of debris-covered glaciers of the Himalayan, this study aims to estimate the recent geodetic mass change of Kanchenjunga ), a large and heavily debris-covered glacier in a data-scarce region in the easternmost Nepal Himalaya (Fig. 1), using multitemporal satellite imagery. This glacier is covered by supraglacial debris mantle from the terminus to 19 km up-glacier, covering ~15 km 2 of the surface area, and has six major tributaries (defined clockwise as T1 to T6, Fig. 1). This study also aims to better understand the effects of 20 supraglacial ponds, flow velocity, and topographic factors on ongoing changes in debris-covered glaciers under the recent climatic change.

Data
Two sets of optical stereoscopic images (hereafter stereo images) obtained using the Hexagon KH-9 and Advanced Land 25 Observing Satellite -Panchromatic Remote-sensing Instrument for Stereo Mapping (ALOS PRISM), are used in this study.
Hexagon KH-9 has a spatial resolution (6 to 9 m) and wide geographic coverage (125 km × 250 km). The mapping camera of Hexagon KH-9 took consecutive nadir images of the ground with approximately 70% overlap (Surazakov and Aizen, 2010;Pieczonka et al., 2013). Three 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 for DEM creation and mapping of the glacier and 30 supraglacial ponds. Similarly, three pairs of ALOS PRISM images (spatial resolution of 2.5 m) and rational polynomial coefficients (RPCs) data acquired from the Remote Sensing Technology Center of Japan were utilized (Table 1)

Glacier delineation and hypsometry
The glacier outlines for 1975 and 2010 were manually delineated using 3-D stereoviews of orthorectified images, the same stereo models that were utilized for generating/editing DEMs. With the help of the high spatial resolution, cloud-free and minimal seasonal snow cover in the chosen images, and 3-D viewing, the identification and delineation of the glacier boundary were feasible. In areas with poor contrast, shade, and steep snow-covered slopes, the topographic (i.e., slopes and 5 contour lines) and geomorphic features (i.e., surface roughness and crevasses) were carefully checked for the mapping based on our experiences in high mountain Asia including the Himalayas (Nuimura et al. 2015;Nagai et al. 2016;Ojha et al., 2016).
The uncertainty associated with the glacier surface delineation was estimated as an error of ±1 pixel along its perimeter (Fujita et al., 2009;Ojha et al., 2016). The glacier was divided into a set of elevation bands with intervals of 50 m to calculate the area-weighted average of the elevation and thus the volume change. Because the DEMs created in this study do not cover 10 the entire glacier, ASTER GDEM2 (Tachikawa et al., 2011) was used to obtain the glacier hypsometry.

DEM generation from ALOS PRISM imagery
The ALOS PRISM images were processed with RPCs data, which contain information about the interior (specifics of the internal geometry of cameras or sensors such as focal length and principle point) and exterior orientations (e.g., position and tilt of camera/sensor) of the acquired images. The use of stereo images and RPCs makes geometric modelling feasible, 15 thereby allowing the creation of DEMs and orthoimages, even without supplying ground control points (GCPs). The ALOS stereo models and triangulated irregular network (TIN) model were generated using the RPCs data and the Leica Photogrammetric Suite (LPS) Workstation. However, automatically extracted TIN models often contain many errors, especially in areas with highly irregular and abruptly changing topography, such as high relief, shadows, and low contrast in the images, leading to an inaccurate terrain representation (Lamsal et al., 2011;Sawagaki et al., 2012). Editing of mass points, 20 which are sets of vertices with XYZ coordinate values defining the vector-based terrain surface, is therefore an important post-process to obtain a precise terrain representation. The stereo MirrorTM/3D Monitor and Leica 3D TopoMouse were used for the GCP collection and terrain editing. The number of mass points that are required to represent a terrain or a particular feature on a terrain surface depends on several factors including the regularity or uniformity and size of the feature, and desired accuracy. The TIN model was extensively edited upon viewing the 3D images until the terrain representation 25 using mass points was satisfactorily achieved. The major editing tasks included the removal of false spikes (mass points in the air) and depressions (mass points below the actual surface). Thus, adequate and representative mass points were placed exactly on the terrain and glacier surfaces including supraglacial ponds and moraine ridges. The LPS Terrain Editor was used to correct and minimize errors of the automatically generated DEMs. Finally, the extensively edited TIN model was gridded into a DEM with a spatial resolution of 5 m (hereafter ALOS-DEM). The detailed procedures are described in Lamsal et al. 30 (2011) and Sawagaki et al. (2012).

DEM generation from Hexagon KH-9 imagery
Hexagon KH-9 images contain distortions due to the development and duplication of the film and long-term storage. The images become suitable for DEM extraction by correcting the distortions with the aid of crosshairs in the images, which allows the accurate reconstruction of the image geometry (Surazakov and Aizen, 2010). Because RPCs are unavailable for 35 4 Hexagon KH-9 images, GCPs were collected from the edited ALOS stereo-model. Because both images used here have a high spatial resolution, terrain features, such as boulders, trail intersections, and sharp notches on moraines, were identifiable.
In total, 21 GCPs were extracted from the off-the-glacier area for which an unchanged surface elevation is expected (Fig. 1).
Five of the 21 GCPs were randomly selected as check points and then used to independently verify the quality of the aerial triangulation (GCPv in Fig. 1). Unlike in the case of the GCPs, the XYZ coordinates of the five GCPv were not used for the 5 aerial triangulation; the image coordinates were used to compute new XYZ coordinates. The root mean square errors of the new and original coordinates of the five GCPv, which are residuals of the aerial triangulation, are 7.7 m horizontally and 4.4 m vertically, respectively. The same TIN editing was performed as that used for the ALOS DEM (Sect. 2.3). The generated DEM was resampled at the same spatial resolution of 5 m (hereafter Hex-DEM) as ALOS-DEM.

10
The upper accumulation area of Kanchenjunga Glacier is extensively covered by snow; the high brightness and poor contrast of both ALOS PRISM and Hexagon KH-9 images preclude the creation of DEMs for the entire glacier. Therefore, DEMs are available only for the lower 22.6 km 2 section (38% of the entire glacier); the debris-covered and -free area is 15.0 km 2 and 7.6 km 2 in size, respectively. Patchy DEM creation is attainable at nine locations in the upper area due to a better local image contrast (hereafter point sites, Fig. 2a). Two assumptions were used to estimate the rate of elevation change in the upper 15 accumulation area higher than 6100 m a.s.l. when DEMs were unavailable: Case 1) the average of rate of elevation change derived from the nine point sites (+0.01 m a -1 ) and Case 2) a fitting curve for the debris-free part above 5800 m a.s.l., where the elevation change is positive (Fig. 3a). Case 1 corresponds to the "no change assumption" in which voids are replaced with zero. (Pieczonka et al., 2013;Maurer et al., 2016), while Case 2 is an alternative of the interpolation approach, in which voids are replaced by the regional mean of the corresponding elevation band (Gardelle et al., 2013). However, the elevation change 20 above 6100 m a.s.l. is unavailable in this analysis so that the upper bound of the exponential fitting curve is assumed to be 0.3 m a -1 , which is based on similar studies showing profiles of the elevation change for glaciers in the eastern Himalayas (Nuimura et al., 2012;King et al., 2016;Maurer et al., 2016). The elevation change of the unmeasured area below 6100 m a.s.l. is assumed to be same as the measured elevation change on the debris-free part. Thus, the two cases were applied to the upper 14.3 km 2 (19% of the entire area). The surface elevation change calculated from the two DEMs (dh/dt, m a -1 ) was 25 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.), r is the density of ice or firn, A z and A T are the areas at a given elevation (50-m elevation 30 band) and of the entire glacier (km 2 ), and r w is the density of water (1000 kg m -3 ). The accumulation area of glaciers often consists of snow and ice, whereas the ablation area mostly contains ice, which means that the elevation change in the ablation area results in denser and thus more mass change than that in the accumulation area. Many previous studies used 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), while two different density assumptions for ablation and accumulation areas were applied in recent studies (e.g., Kääb et al., 2012; 35 5 Pellicciotti et al., 2015). In this study, two density scenarios are considered: Scenario 1) 900 kg m -3 for the entire glacier area, and Scenario 2) 900 and 600 kg m -3 for the ablation and accumulation areas, respectively. For Scenario 2, an equilibrium line altitude (ELA) is required to divide the accumulation and ablation areas, though it is unmeasured. The ELA for the studied period, above which the elevation change of debris-free surface is positive, is assumed to be 5850 m a.s.l. (Fig 3a). The mean of the four combinations, which are based on two cases of elevation change in the unmeasured accumulation area and two 5 density scenarios, is assumed to be the most plausible geodetic mass balance .
The accuracy of the geodetic mass balance (σ g ) is quantified by considering two sources of uncertainty as: (2) 10 where r 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 of 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 of the off-glacier area (red polygon in Fig. 2a) and then uniformly applied to the area below 6100 m a.s.l. for which measured data are available (Fig. 3a). We believe that the TIN editing performed on the Hexagon KH-9 and ALOS PRISM images guarantees 15 the same degree of uncertainty over the measured area and 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 the use of standard deviations of the glacier elevation change within a given elevation band to be an extrapolation error; this approach is often used in studies (e.g., Maurer et al., 2016). In this study, however, the estimation of this extrapolation error is hampered by the bright and poorly contrasted snowfields.
Therefore, as an alternative for the extrapolation error, the error based on the two assumptions of the elevation change of the 20 upper accumulation area (σ z ) is defined as the difference of the elevation change between Case 1 (+0.01 m a -1 ) and Case 2 (fitting curve in Fig. 3a).

Glacier velocity
To investigate the effects of the topography and dynamics on glacier elevation changes, the distribution of the surface velocity of Kanchenjunga Glacier was calculated using a feature tracking method (Heid and Kääb, 2012). Orthorectified pairs 25 of ALOS PRISM images acquired in March and December of 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 applications in deriving the terrain displacement including the glacier velocity 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 allows the computation of the relative surface displacement between the initial (reference image) and final image (search image). The stereo images were firstly 30 orthorectified and then co-registered to ensure that the corresponding pixels in each image overlap exactly, which is required to initialize the matching process. The glacier velocity was computed using a correlation window of 64 × 64 pixels, corresponding to 160 × 160 m, a robustness iteration of 4, and a mask threshold for the noise reduction of 0.9 (Leprince et al., 2007). To ensure the quality of the velocity map, poor matching in the surface displacements was removed by applying a correlation threshold of 0.6, which results in some voids in the shaded area (Fig. 2b). The details of the process and 6 methodological background of COSI-Corr and its application to glacier-velocity computations are described in Leprince et al. (2007) and Scherler et al. (2008). The uncertainty of the flow velocity was evaluated using the displacement of the off-glacier area by constraining a surface slope gentler than 25°.

Supraglacial ponds
Supraglacial ponds and ice cliffs can enhance the melt of debris-covered ice by absorbing radiative heat as hot spots (Sakai et 5 al., 2000(Sakai et 5 al., , 2002Steiner et al., 2015). Therefore, their distribution density is important to better understand the mechanisms of glacier surface lowering. Based on the TIN editing used to create the above-mentioned DEMs, we simultaneously delineated supraglacial ponds on the glacier, which are easily identifiable as flat terrain features in the stereo images, whereas it is difficult to distinguish the ice cliff and debris-covered steep slopes in panchromatic images. To avoid the misinterpretation of the topography, we delineated the ponds greater than 0.001 km 2 (12 × 12 pixels of ALOS image). For comparison, the pond 10 delineation was also conducted for Khumbu Glacier using the ALOS PRISM images taken in October 2008 (Table 1). We The distribution map of the rate of elevation change was derived from Hex-DEM and ALOS-DEM (Fig.2a). The most pronounced surface lowering was observed between 4700 and 5500 m a.s.l. (Fig. 3a). Overall, the glacier surface slightly lowered in the lowermost ablation area, remarkably lowered in the middle ablation area, and slightly thickened in the upper 25 area, though spatial variability is notable. Despite the debris cover in the ablation zone, which causes an insulation effect, the glacier surface lowering is remarkable, except for the lowermost reach (Figs. 2a and 3a). The debris-covered section of the glacier extends from 4500 to 5850 m a.s.l., covering 14.9 km 2 (25% of the entire glacier); the area-weighted average of the rate of elevation change is -0.55 m a -1 . Furthermore, the debris-covered area shows a more negative trend (-0.51 m a -1 ) than that of debris-free ice at ~5100-5850 m a.s.l. (-0.30 m a -1 ), where the two rates coexist (Fig. 3a). Although the lowering of 30 the glacier surface is affected not only by surface melting but also by glacier dynamics, this suggests that the debris-covered surface has also experienced thinning comparable to that of the debris-free ice surface, which is consistent with the results for regions in the Himalayas Nuimura et al., 2012;Gardelle et al., 2013).
To estimate the area-weighted geodetic mass balance of the entire glacier, the measured rate of elevation change at respective elevation bands was used for the lower ablation area (Fig. 3a). The two assumptions of the elevation change (constant value 7 of +0.01 m a -1 as Case 1 and a fitting curve in Fig. 3a as Case 2) and the two density scenarios were adopted for the unmeasured upper accumulation area above 6100 m a.s.l. (19% of the entire area). The geodetic mass balance of Kanchenjunga Glacier was finally estimated to be -0.18 ± 0.17 m w.e. a -1 for the period 1975-2010; it ranges from -0.15 to -0.20 m w.e. a -1 (Table 2).
We tested the sensitivity of the geodetic mass balance to the assumptions. The density scenario does not significantly affect 5 the resulted mass balance in both cases (< 0.02 m w.e. a -1 , Table 2). 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 greatest influence is derived from how the elevation change in the unmeasured area is assumed (Cases 1 or 2), the difference of the final mass balance is ~0.04 m w.e. a -1 (Table 2), which is one fourth of the uncertainty (±0.17 m w.e. a -1 ). where glacier surface changes from debris-covered to debris-free ones (Fig. 4b). This is also true for the longitudinal slope, which is gentler than 15° in the entire area (Fig. 4c). The flow velocity of the main tributary (T3) up to 18 km from the terminus is slow (< 20 m a -1 ), whereas the flow velocity of the tributaries significantly increases toward the up-glacier, except for T2 (Fig. 4d). The longitudinal pattern of the flow velocity is similar to that of the slope gradient (Figs. 4c and 4d). The gradient of the flow velocity, where a positive value implies divergence and vice versa, is noisy without a significant trend, 35 though a compressive (negative) trend is expected in the ablation area (Fig. 4e).
Here, we compare the rate of elevation change with the elevation, slope and flow velocity, which are expected to relate to the melting of glacier ice and the emergence velocity determining the elevation change of the glacier surface. In general, the mass balance of the debris-free glacier is more negative at lower elevations. On the other hand, the debris thickness increases down-glacier with considerable spatial variability (Kirkbride and Warren, 1999;Nakawo et al., 1999). The melt rate is high if the debris cover is thin because of the effect of the low albedo, while it greatly decreases under thicker debris because of the 5 insulating effect (Mattson et al., 1993). Consequently, the high melting rate of the debris-covered glaciers is typically found in the mid-ablation areas characterized by a thin debris cover, while the melting rate substantially decreases at the lowermost terminus because of thicker debris and toward higher elevations because of cooler air temperatures. Both elevation and longitudinal profiles of the rate of elevation change follow this trend (Figs. 3a and 4a). Although the rates of elevation change along the tributaries show different profiles, the gradients against the elevation are similar (1.86 ± 0.56 × 10 -3 m a -1 m -1 ), 10 suggesting that the glacier thickness decreases towards the glacier confluences (thick solid lines in Fig. 5a). On the other hand, the rate of elevation change increases toward zero along the main tributary (T3) below 5400 m a.s.l. (thick black broken line in Fig. 5a). Despite this large variability, the rate of elevation change and elevation show a significant correlation in the entire domain, where data are available (n = 197, r = 0.35, p < 0.001; thin line in Fig. 5a).
In the rigid definition, the difference among the ice fluxes at a given section of glacier, for which the flow velocity and ice 15 thickness are required generates the emergence velocity. In this regard, the gradient of the flow velocity is expected to correlate with the rate of elevation change (less gradient, more lowering). However, the gradient of the flow velocity is very noisy (Fig. 4e) and shows no relationship with the rate of elevation change (not shown), while the flow velocity exhibits a statistically significant correlation with the rate of elevation change (n = 195, r = 0.48, p < 0.001; Fig. 5b). On the other hand, the large variability of the rate of elevation change at slow speed surfaces should be affected by the large variability of the 20 surface mass balance because the emergence velocity is expected to be negligible. The area with a slow flow velocity (< 10 m a -1 ) extends up to 19 km from the terminus (Fig. 4d); the surface lowering decreases toward the terminus (Figs. 4a and 5a).
This suggests that the distribution of the debris thickness more likely results in the observed distribution of the surface lowering of the debris-covered area.
The longitudinal slope statistically correlates with the rate of elevation change (n = 186, r = 0.40, p < 0.001; Fig. 5c). This is 25 due to a significant relationship between the slope and flow velocity (n = 186, r = 0.65, p < 0.001; Fig. 5d), which is generally expected based on glacier dynamics (Cuffey and Paterson, 2010).

Role of supraglacial ponds on the surface lowering
During the study period from 1975 to 2010, the supraglacial ponds of Kanchenjunga Glacier increased in number by a factor > 4 from 8 to 35 and in average size from 1710 to 4580 m 2 ( Table 3). Most of the ponds formed in the mid-ablation area 30 between 4700 and 5400 m a.s.l. or between 5 and 18 km from the terminus. The lifetime of a supraglacial pond is expected to be a few to several years (Miles et al., 2017), though it is not fully investigated on the decadal scale. Therefore, the 1975 ponds disappeared either by deformation due to glacier dynamics or by drainage through englacial channels (Sakai et al., 2000). We observed a significant negative correlation between the size of supraglacial ponds in 2010 and rate of elevation change, indicating that at sites, where larger supraglacial ponds exist, the glacier surface is more down-wasted than at sites 35 with smaller ponds (n = 35, r = 0.45, p < 0.01; Fig. 6a). In addition, many ponds exhibit a greater lowering than the lower bound of the debris-covered surface (Fig. 6b). The illustrated box plots in Fig. 6b show that the rate of elevation change for supraglacial ponds (-1.25 ± 0.34 m a -1 ) is characterized by more lowering than that of the other debris-covered surfaces up to 20 km from the terminus (-0.74 ± 0.42 m a -1 ).
Supraglacial ponds and ice cliffs, which often coexist, are the principal spots of extensive heat absorption (Sakai et al., 2000(Sakai et al., , 2002Röhl, 2008;Steiner et al., 2015). The significant lowering of the glacier surface, where supraglacial ponds exist, which 5 do not persist at the same location for decades but for a few years (Miles et al., 2017), suggests that these hot spots may repeat a development and decay cycle over a few to several years by changing their location; these ponds enhance melting, and then leave their traces in form of significant lowering. This mechanism might be one of the causes of the surface lowering of debris-covered areas comparable to that of a debris-free glacier.
The substantial lowering of the mid debris-covered surface (Fig. 4a) could initiate the formation of additional ponds or the 10 expansion of existing ponds because the slow flow velocity allows the ponds to stay at the same locations (Fig. 4d) (Miles et al., 2017). In the Langtang region, many ponds on debris-covered glaciers were found in the stagnant area with a small glacier velocity (Ragettli et al., 2016); they persisted for several years (Miles et al., 2017). Watson et al. (2016) investigated nine debris-covered glaciers in the Everest region and reported a net increase of the pond area with large inter-and intraannual variability. Recently, the total number, surface area, and average size of the supraglacial ponds significantly increased 15 in the mid-ablation area (Table 3). This further enhanced the surface lowering by increased melting. Considering the high rate of surface lowering, negligible flow velocity, gentle surface slope, and several existing supraglacial ponds of notable size, a large supraglacial lake is expected to form in the mid debris-covered area when the calving of the glacier terminus begins (Röhl, 2008;Sakai et al., 2009;Sakai and Fujita, 2010).

20
In this section, we compare the geodetic mass balance of Kanchenjunga Glacier (-0.18 ± 0.17 m w.e. a -1 ) with that of other studies analyzing Himalayan glaciers with similar remotely sensed data (e.g., Corona or Hexagon spy satellites and recent satellites). Khumbu Glacier is the most investigated "iconic" debris-covered glacier in the Himalayas, which exhibits a more negative mass balance for the period 1970-2007 (-0.27 ± 0.08 m w.e. a -1 ; Bolch et al., 2011) than that of Kanchenjunga Glacier. We discuss two possible causes of the differences: 1) contributions of accumulation and debris-covered areas and 2) 25 the density of supraglacial ponds. We use the hypsometry of Khumbu Glacier delineated in Nuimura et al. (2012). If we assume a threshold elevation of 5800 m a.s.l., above which the glacier surface was lifted in the study period, the area above the threshold elevation (54% of the entire glacier) contributes to suppress the negative mass balance of Kanchenjunga Glacier, which is larger than that of Khumbu Glacier (41%). On the other hand, the debris-covered area ratio of Kanchenjunga Glacier (25%, 15.0 km 2 ) is much smaller than that of Khumbu Glacier (40%, 8.2 km 2 ). A greater accumulation area and smaller 30 debris-covered area in the entire region of Kanchenjunga Glacier would cause a less negative geodetic mass balance compared with that of Khumbu Glacier. In addition, supraglacial ponds might play a key role in the surface lowering of the debris-covered area. Supraglacial ponds on Khumbu Glacier delineated from ALOS PRISM images (October 2008) were therefore compared with those of Kanchenjunga Glacier (Table 3). The major variables of supraglacial ponds, such as number, total area, pond area ratio to debris-covered area, and pond density, are roughly two to four times greater for 35 Khumbu Glacier than that of Kanchenjunga Glacier, though the average pond size is similar (Table 3). This high pond density likely contributes to the more negative geodetic mass balance of Khumbu Glacier. The pond density and inter-annual variability in size and location have been recently investigated in the Langtang (Miles et al., 2016) and Everest regions (Watson et al., 2016) using high-resolution remote sensing images. This type of analysis in wider regions might be an important factor for regionally contrasting changes in debris-covered glaciers of the Himalayas-Karakorum range Gardelle et al., 2013).

5
We further compared the geodetic mass balance of Kanchenjunga Glacier with other regionally analyzed debris-covered glaciers. By excluding debris-free and lake-terminating glaciers, five glaciers in Langtang (1974-2006, Ragettli et al., 2016), eight glaciers in Khumbu (1970-2007, Bolch et al., 2011, and five glaciers in Bhutan (1974-2006, Maurer et al., 2016 were selected (Fig. 7). The geodetic mass balance of rather small glaciers (< 20 km 2 ) shows a large variability, while that of the four large glaciers (two Bhutanese glaciers and Langtang and Kanchenjuga Glaciers) exhibits a moderate mass loss (-0.15 to 10 -0.24 m w.e. a -1 ). Because of their large size, the area-weighted regional means of the geodetic mass balance have values similar to those of the large glaciers (open squares in Fig. 7). Therefore, the mass balance of Kanchenjunga Glacier could be representative for the regional mean in the easternmost Nepal Himalaya, though this study deals with a single glacier due to the time-consuming TIN editing method.

15
This study provides the geodetic mass balance of Kanchenjunga Glacier, a large debris-covered glacier situated in the easternmost Nepal Himalaya, using two high resolution DEMs (5 m) spanning 35 years, created from Hexagon KH-9 (1975) and ALOS PRISM (2010) stereo images. The rate of elevation change of the debris-covered surface is similar to that of the debris-free surface. Considering the lowering in the mid-ablation area, the debris-covered surface also experienced thinning comparable to that of the debris-free glacier surface. The rate of elevation change shows a significant correlation with the 20 flow velocity. The large variability of the rate of elevation change at slow speeds suggests that the surface mass balance, which may be affected by the debris distribution, is more significant for the observed distribution of surface lowering of Kangchenjunga Glacier. Along the longitudinal profile of the rate of elevation change, many areas, including ponds, show a greater lowering than in the debris-covered area. This supports the assumption that the supraglacial ponds and associated ice cliffs might be a major constituent of the intense ice wastage and play a key role in the heterogeneous surface lowering of 25 debris-covered glaciers. The lower pond density might also be important for the less negative geodetic mass balance of Kanchenjunga Glacier compared with Khumbu Glacier.
Although the generated DEM covers only 40% of the entire glacier and thus the uncertainty of the estimated mass balance is large, we think that the TIN editing method results in a better accuracy for the generated DEM (   . Geodetic mass balance of the debris-covered glaciers (solid circles) in the Himalayas obtained from the difference of two DEMs generated from declassified spy satellites and recent ones. The red, blue, green, and black circles denote eight glaciers in Khumbu (Bolch et al., 2011), five glaciers in Langtang (Ragettli et al., 2016), five glaciers from Bhutan (Maurer et al., 2016), and Kanchenjunga Glacier (this study), respectively. The open squares with thick error bars are the area-weighted 5 regional means and standard deviations. Note that the error bar of Kanchenjunga Glacier is the standard deviation, while the others represent the standard error.
Half width