Spatial patterns in glacier characteristics and area changes from 1962 to 2006 in the Kanchenjunga–Sikkim area, eastern Himalaya

. This study investigates spatial patterns in glacier characteristics and area changes at decadal scales in the eastern Himalaya – Nepal (Arun and Tamor basins), India (Teesta basin in Sikkim) and parts of China and Bhutan – based on various satellite imagery: Corona KH4 imagery, Landsat 7 Enhanced Thematic Mapper Plus (ETM + ) and Advanced Spaceborne Thermal Emission Radiometer (ASTER), QuickBird (QB) and WorldView-2 (WV2). We compare and contrast glacier surface area changes over the period of 1962–2000/2006 and their dependency on glacier topography (elevation, slope, aspect, percent debris cover) and climate (solar radiation, precipitation) on the eastern side of the topographic barrier (Sikkim) versus the western side (Nepal).


Introduction
Himalayan glaciers have generated a lot of concern in the last few years, particularly with respect to potential consequences of glacier changes on the regional water cycle (Kaser et al., 2010;Immerzeel et al., 2010Immerzeel et al., , 2012Racoviteanu et al., 2013a). In the last decades, the availability of low-cost data from optical remote sensing platforms with global coverage provided opportunities for glacier mapping at regional scales. Remote sensing techniques have considerably helped improve estimates of glacier area changes (Bhambri et al., 2010;Bolch, 2007;Bajracharya et al., 2007;Kamp et al., 2011;Bolch et al., 2008a), glacier lake changes (Bajracharya et al., 2007;Bolch et al., 2008b;Gardelle et al., 2011;Wessels et al., 2002) and region-wide glacier mass balance Bolch et al., 2011;Gardelle et al., 2013;Kääb et al., 2012), but significant gaps do remain. The new global Randolph Glacier Inventory (RGI) v.4 (Pfeffer et al., 2014) provides a global data set of glacier outlines intended for largescale studies; however, in some regions the quality varies and the outlines may not be suitable for detailed regional analysis of glacier parameters. A new Landsat-based inventory Published by Copernicus Publications on behalf of the European Geosciences Union.
has been complied using imagery from 1999 to 2003 that, along with the current study, may help improve the accuracy in some areas of RGI (Nuimura et al., 2014). Some other regional glacier inventories have been constructed in the past, for example for the western part of the Himalaya (e.g., Kamp et al., 2011;Bhambri et al., 2011;Frey et al., 2012), but only a few are available for the eastern extremity of the Himalaya (e.g., Krishna, 2005;Basnett et al., 2013;Bajracharya et al., 2014;Kulkarni and Narain, 1990). The use of remote sensing for glacier mapping in this area is limited by frequent cloud cover and sensor saturation due to unsuitable gain settings and the persistence of seasonal snow, which hampers quality satellite image acquisition. Furthermore, this area has very limited reliable baseline topographic data needed for glacier change detection, as discussed in detail in Bhambri and Bolch (2009). The earliest Indian glacier maps date from topographic surveys conducted by expeditions in the mid-nineteenth century (i.e. Mason, 1954), but these are limited to a few glaciers. The Geologic Survey of India (GSI) inventory based on 1970s Survey of India maps (Sangewar and Shukla, 2009;Shanker, 2001) is not in the public domain. For eastern Nepal, 1970s topographic maps from Survey of India on a 1 : 63 000 scale are available, but their accuracy is not known with certainty. Given these limitations, declassified Corona imagery from the 1960s and 1970s has increasingly been used to develop baseline glacier data sets, for example in the Tien Shan (Narama et al., 2010), Nepal Himalaya (Bolch et al., 2008a) and parts of Sikkim Himalaya (Raj et al., 2013).
With a wide variety of satellite data becoming available, the topographic and climatic controls on glacier surface area have received increasing attention in recent studies, particularly with respect to debris-covered glaciers (Basnett et al., 2013;Thakuri et al., 2014;Bolch et al., 2008a;Salerno et al., 2008). Some studies have characterized the small-scale glacier surface topography of debris-covered glaciers using field-based surveys (Iwata et al., 2000;Sakai and Fujita, 2010;Zhang et al., 2011), while other studies focused on understanding patterns at the mountain-range scale (Scherler et al., 2011;Bolch et al., 2012;Gardelle et al., 2013;Racoviteanu et al., 2014). Glacier shrinkage and mass loss has been documented in the Himalaya concomitantly with an increase in debris cover Nuimura et al., 2012). However, the influence of debris cover on glacier mass balance remains debatable (Scherler et al., 2011;Kääb et al., 2012), and modeling of melt under the debris cover is subject to uncertainties due to limited field-based measurements of debris thickness needed for model parameterization (Zhang et al., 2011;Foster et al., 2012;Mihalcea et al., 2008a, b).
While significant progress has been made in recent years in remote sensing glacier mapping in the Himalaya, some of the subregions still need updated glacier area and surface characteristics including debris cover. The objective of this study is twofold: (1) to present the current glacier distribution and characteristics in a data-scarce area of the eastern Himalaya based on an updated 2000 Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and Advanced Spaceborne Thermal Emission Radiometer (ASTER) inventory, along with elevation data from the Shuttle Radar Topography Mission (SRTM); and (2) to investigate spatial patterns in glacier surface area changes from 1962 (Corona KH4) to 2000 (Landsat/ASTER) and 2006 QuickBird (QB)/2009 WorldView-2 (WV2) and their dependence on topographic and climatic factors, with a particular emphasis on debriscovered glacier tongues. These updated glacier data sets will help filling a gap in global glacier inventories such as the RGI The Cryosphere, 9, 505-523, 2015 www.the-cryosphere.net/9/505/2015/ A. E. Racoviteanu et al.: Patterns in glacier characteristics and area changes 507 (Pfeffer et al., 2014), as well as for subsequent future mass balance applications at regional scales.

Study area
The study area encompasses glaciers in the eastern Himalaya (27 • 04 52 -28 • 08 26 N latitude and 88 • 00 57 -88 • 55 50 E longitude), located on either side of the border between Nepal and India in the Kanchenjunga-Sikkim area (Fig. 1). Based on SRTM data, relief in this area ranges from 300 m at the bottom of the valleys to 8598 m (Mt. Kanchenjunga). Valley glaciers cover about 68 % of the glacierized area, mountain glaciers cover 28 % and the remaining area is made up of cirque glaciers and aprons (Mool et al., 2002). The glacier ablation area is typically covered by heavy debris cover originating from rockfall on the steep slopes (Mool et al., 2002), reaching up to a thickness of several meters at the glacier termini (Kayastha et al., 2000). The eastern part of this area constitutes the Sikkim province of India, and the western part is located in eastern Nepal and encompasses the Tamor basin and parts of the Arun basin. Climatically, this area of the Himalaya is dominated by the South Asian summer monsoon circulation system (Bhatt and Nakamura, 2005) caused by the inflow of moist air from the Bay of Bengal to the Indian subcontinent during the summer (Benn and Owen, 1998;Yanai et al., 1992). The Himalaya and Tibetan plateau act as a barrier to the monsoon winds, bringing about 77 % of precipitation on the south slopes of the Himalaya during the summer months (May-September) (Fig. 2). This climatic particularity causes a "summer-accumulation" glacier regime type, with accumulation and ablation occurring simultaneously in the summer (Ageta and Higuchi, 1984). In Sikkim, rainfall amounts range from 500 to 5000 mm yr −1 , with annual averages of 3580 mm recorded at Gangtok station (1812 m;1951-1980) (IMD, 1980, and 164 rainy days per year (Nandy et al., 2006). Mean minimum and maximum daily temperatures at Gangtok station were reported as 11.3 and 19.8 • C, with an average of 15.5 • C based on the same observation record (IMD, 1980).

Satellite imagery
Remote sensing data sets used in this study are summarized in Table 1 and included (1) baseline remote sensing data from Corona declassified imagery (year 1962),  The graph shows the monsoon period from June to September, with a peak precipitation in July, and the influence of the northeastern monsoon during the winter/early spring (January-March).
1. Corona KH4 scenes (1962) were obtained from the US Geological Survey EROS Data Center (USGS-EROS, 1962). The Corona KH4 system was equipped with two panoramic cameras (forward-looking and rearlooking with 30 • separation angle) and acquired imagery from February 1962 to December 1963 (Dashora et al., 2007). We chose images from the end of the ablation season (October/November in this part of the Himalaya), suitable for glaciologic purposes. Six Corona stripes were scanned at 7 microns by USGS from the original film strips, with a reported nominal ground resolution of 7.62 m (Dashora et al., 2007). Corona images are known to contain significant geometric distortions due to cross-path panoramic scanning. The Frame Ephemeris Camera and Orbital Data camera/spacecraft parameters (roll, pitch, yaw, speed, altitude, azimuth, sun angle and film scanning rate) for Corona missions, needed to construct a camera model and to correct these distortions, are not easily available. To orthorectify the scenes, we defined a non-metric camera model in ERDAS Leica Photogrammetric Suite (LPS), with focal length, air photo scale and flight altitude extracted from the declassified documentation of the KH4 mission (Dashora et al., 2007). We used the bundle block adjustment procedure in LPS to simultaneously estimate the orientation of all the CORONA stripes on the basis of 117 ground control points (GCPs). Latitude and longitude (x, y) information (of the GCPs) were extracted from the panchromatic band of the 2000 Landsat ETM+ image (15 m spatial resolution) on nonglacierized terrain including moraines, river crossings and outwash areas, whereas elevation information (z) were extracted from the SRTM DEM v.4 (CGIAR-CSI 2004).
www.the-cryosphere.net/9/505/2015/ The Cryosphere, 9, 505-523, 2015 Tie points were automatically extracted in LPS on overlapping Corona strips, and visually checked on the Landsat image. The Corona stripes were mosaicked in ERDAS LPS to produce the final orthorectified image, with a horizontal accuracy (RMSE x, y) of the bundle block adjustment of 10.5 m. The orthorectification process of the 1962 Corona yielded a RMSE x, y error of ±10 m and the actual "ground" RMSE x, y of the Corona block of ∼ 60 m. A trend analysis on the horizontal shifts between Corona and the reference Landsat scene showed that the largest errors occurred towards the edges of the images, mostly outside the glaciers, and did not impact the area change analysis.
2. The orthorectified Landsat ETM+ scene from December 2000, obtained from the USGS Eros Data Center, was the main data set for the updated glacier inventory. In addition, six orthorectified ASTER products (2000)(2001)(2002) were obtained at no cost through the Global Land Ice Monitoring from Space (GLIMS) project . Images were selected at the end of the ablation season for minimal snow and had little or no clouds. Five of these scenes were used for on-screen manual corrections of the Landsat-based glacier outlines in challenging areas where shadows or clouds obstructed the view of the glaciers. In addition, the surface kinetic temperature product (AST08) from the 27 November 2001 ASTER scene was used for clean ice delineation of debris cover along with topographic information using a decision-tree algorithm (Racoviteanu and Williams, 2012). The 29 October 2002 scene, covering the Kanchenjunga-Sikkim area east and west of the topographic divide, was used to investigate the spatial distribution of surface temperature over selected debriscovered tongues.
3. Two QB scenes from January 2006 were obtained from Digital Globe as ortho-ready standard imagery (radiometrically calibrated and corrected for sensor and platform distortions) (Digital Globe, 2007). These scenes, covering an area of 1107 km 2 were well contrasted and mostly snow-free outside the glaciers. We orthorectified these scenes in ERDAS Imagine Leica Photogrammetry Suite (LPS) using rational polynomial coefficients (RPCs) provided by Digital Globe and the SRTM DEM and mosaicked them in ERDAS Imagine. The scenes were resampled to 3 m pixel size during the orthorectification process using the cubic convolution method suitable for continuous raster data in order to reduce disk space and processing time. One WV2 panchromatic, ortho-ready scene at 50 cm spatial resolution from 2 December 2010 was also obtained to cover the terminus of Zemu glacier, which was missing from the QB extent. All data sets were registered to UTM projection zone 45N, with elevations referenced to the WGS84 datum.

Elevation data sets
Two elevation data sets were used in this study: 1. The hydrologically sound, void-filled CGIAR SRTM DEM (90 m spatial resolution) (CGIAR-CSI, 2004) was used to extract glacier parameters for 2000. The SRTM data set is known to have biases on steep slopes and at higher elevations (Fujita et al., 2008;Berthier et al., 2006;Nuth and Kääb, 2011), as well as due to radar penetration on snow (Gardelle et al., 2012b). For this area, the vertical accuracy of the SRTM DEM, calculated as root mean square (RMSE z) with respect to 25 field-based GCPs, was 31 m ± 10 m. The GCPs were obtained in the field on non-glacierized terrain including roads and bare land outside the glaciers using a Trimble Geoexplorer XE series.
The Cryosphere, 9, 505-523, 2015 www.the-cryosphere.net/9/505/2015/  2. The Swiss topographic map (1 : 150 000 scale), compiled from Survey of India maps from the 1960s and published by the Swiss Foundation for Alpine Research, was used for manual corrections of the 1962 Corona glacier outlines to discard any seasonal snow and to correct shadow areas or bright water bodies that could be misclassified as ice. The exact month or year of each quadrant or of the original air photos is not known with certainty because the original large-scale Indian topographic maps at this scale are restricted to within 100 km of the Indian border and are therefore inaccessible (Srikantia, 2000;Survey of India, 2005); however this map was useful for manual corrections of the Corona outlines.

Analysis extents
We defined three analysis extents for our study area ( Fig. 1 and Table 2): 1. The Landsat/ASTER domain includes the Sikkim province of India, parts of eastern Nepal (Tamor and Arun basins) and parts of Bhutan and China (Table 2). This domain was used to construct the updated 2000 glacier inventory.
The Landsat/ASTER domain was split into four subregions on the basis of east-west and north-south climate/topographic/political barriers, as shown in Fig. 3. Rainfall averages from the Tropical Rainfall Measuring Mission (TRMM) data 2B31 product (Bhatt and Nakamura, 2005;Bookhagen and Burbank, 2006) were used to characterize the subregions climatically. The data set contains rainfall estimates calibrated with groundcontrol stations derived from local and global gauge stations (Bookhagen and Burbank, 2006) with a spatial resolution of 0.4 • , or ∼ 5 km. Given the well-known biases in the TRMM data (Bookhagen and Burbank, 2006;Palazzi et al., 2013;Andermann et al., 2011), here we are not concerned with the absolute values of gridded precipitation but only with characterizing the subregions in our study area using relative rainfall values. TRMM data integrated over 10 years (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) show differences in precipitation patterns among the four regions and justifies our choice of spatial domains (Table 3). The eastern side of the study area (Sikkim) receives higher precipitation amounts than the western side (Nepal) (977 versus 805 mm yr −1 ). There is a pronounced north-south gradient in precipitation, with the lowest amount of precipitation noticeable on the Chinese side (146 mm yr −1 ) (

Glacier delineation and analysis
For the 1960s, clean glacier outlines were extracted from the panchromatic Corona imagery by thresholding the digital numbers (DN > 200 = snow/ice), chosen based on visual interpretation. Debris-covered glacier tongues were delineated manually on the basis of lateral moraines and other visual clues such as supraglacial lakes. A 5 × 5 median filter was used to remove noise (isolated pixels from snowfields or internal rocks), as recommended in other studies (Racoviteanu et al., 2009;Andreassen et al., 2008). Ice polygons with area < 0.02 km 2 were not considered valid glaciers and were excluded from the analysis. Manual corrections were applied subsequently on the basis of the topographic map using on-screen digitizing in areas of poor contrast or transient snow/clouds, which obstructed the view of glaciers. For the 2000s, glaciers were delineated from the Landsat ETM+ scene using the normalized difference snow index (NDSI) (Hall et al., 1995), with a threshold of 0.7 (NDSI > 0.7 = snow/ice). The NDSI algorithm relies on the high reflectivity of snow and ice in the visible to near infrared wavelengths (0.4-1.2 µm), compared to their low reflectivity in the shortwave infrared (1.4-2.5 µm) (Dozier, 1989;Rees, 2003). Compared to other band ratios (Landsat 3/4 and 3/5), the NDSI glacier map was cleaner and less noisy and was therefore preferred (Racoviteanu et al., 2008b). A 5 × 5 median filter was used here as well to remove remaining noise, and a few areas were adjusted manually on the basis of the ASTER images, notably frozen lakes misclassified as snow/ice and some glaciers underneath low clouds in the southern part of the image. Some transient snow persisting in the deep shadowed valleys was manually removed from the glacier outlines on the basis of the topographic map. Debris-covered glacier tongues were delineated using multispectral data (band ratios, surface kinetic temperature and texture) from the 27 November 2001 ASTER scene com-bined with topographic variables in a decision tree, as described in Racoviteanu and Williams (2012).
For the QB (2006) image, clean ice surfaces were delineated using band ratios 3/4, then ISODATA clustering with a threshold of 1.07 (snow/ice > 1.07) and a majority filter of 7 × 7 to remove noise. Debris-covered tongues for this data set were delineated manually on the basis of supraglacial features (lakes and ice walls), along with lateral and frontal moraines visible on the high-resolution images. We also mapped supraglacial lakes from this high-resolution data based on band ratios along with texture analysis.
For all inventories in the Landsat/ASTER domain, ice masses were separated into glaciers on the basis of the SRTM DEM, using hydrologic functions in an algorithm developed by Manley (2008), described in Racoviteanu et al. (2009). Glacier area, terminus elevation, maximum and median elevation, average slope angle and aspect were extracted on a glacier-by-glacier basis using zonal functions on the SRTM DEM. Average glacier thickness and length were calculated from mass turnover principles and ice flow mechanics by Huss and Farinotti (2012), based on the approach of Farinotti (2009). Their method used our glacier outlines and the SRTM DEM to derive thickness estimates iteratively based on Glen's flow law and a shape factor (Paterson, 1994). For simplicity and consistency for change analysis, we assumed no shift in the ice divides over the period of analysis and excluded all nunataks and snow-free steep rock walls from the glacier area calculations. Bodies of ice above the bergschrund were considered part of the glacier (Racoviteanu et al., 2009;Raup and Khalsa, 2007). Glacier area changes  and their dependency on topographic and climatic variables were calculated on a glacier-by-glacier basis for the 232 glaciers in spatial domain 2 using linear regression.

Uncertainty estimates
Glacier outlines derived from remote sensing data at various spatial and temporal resolution are subject to various degrees of uncertainty, as discussed in recent studies (Paul et al., 2013;Racoviteanu et al., 2009). This becomes an important issue in glacier change analysis, where errors from various data sources accumulate at each processing step. The main sources of uncertainty considered here are (1) image classification errors (positional errors and/or errors due to the semi-automated glacier mapping method) and (2) conceptual errors associated with the definition of a glacier, including mapping of ice divides, mixed pixels of snow and clouds and internal rock differences, which propagate to the glacier change analysis, all described in detail in Racoviteanu et al. (2009).
1. The errors in remote sensing glacier surface areas due to classification (E classif ) were estimated using Perkal's epsilon band around each glacier outline data The Cryosphere, 9, 505-523, 2015 www.the-cryosphere.net/9/505/2015/ set (Racoviteanu et al., 2009;Bolch et al., 2010), with a ∼ 1-pixel variability (Congalton, 1991). Using ±30 m for Landsat/ASTER, ±6 m for Corona and ±3 m for QB outlines, the area uncertainty was ±3, ±6 and ±2 % of the glacierized area for Corona, Landsat/ASTER and QuickBird, respectively. The Perkal method is known to slightly overestimate the errors, as described in Burrough and McDonnel (1998). Recent glacier analysis comparison experiments reported a range of uncertainty of < 5 % for remote sensing glacier outlines compared to high-resolution imagery Paul et al., 2013). For manually adjusted glacier outlines, particularly debris-covered tongues, we used screen digitizing in streaming mode with a high density of vertices to minimize area errors (B. Raup, National Snow and Ice Data Center, personal communication, 2014).
2. Uncertainties due to different digitization of internal rocks (E rocks ) were derived by comparing area changes computed with internal rocks specific to each data set versus "merged" internal rocks from all data sets. The differences in glacier data sets due to rock inconsistencies amounted to ∼ 2 % of the glacier area. To minimize uncertainties in the glacier area change, we merged rock outcrops from each data set and removed them from all the area calculations. The "inactive" bodies of ice above the bergschrund were included as part of the glacier (Racoviteanu et al., 2009). For simplicity, we neglected the area change that might be due to exposure of new internal rock due to glacier ice thinning.
Total errors in glacier area estimate for each data set (E) were calculated as RMSE of the classification (E classif ) and the internal rocks (E rock ): Errors in glacier surface area change (E change ) from 1962 to 2000 were computed as RMSE of the total error for each time step calculated in Eq. 1: 4 Results

The 2000 Landsat/ASTER glacier characteristics
The 2000 glacier inventory based on Landsat and ASTER yielded 487 glaciers (of which 162 were situated in Nepal, 186 in Sikkim, 30 in Bhutan and 109 in China), covering a total surface area of 1463 ± 88 km 2 (Table 4a). Of the 487 glaciers in this spatial domain, 68 glaciers (13 %) had debris cover on their ablation areas. Supraglacial debris covered 160 ± 10 km 2 (11 % of the glacierized area in spatial domain 1), with some differences between north and south slopes of the study area (discussed in Sect. 5.1). In Sikkim, supraglacial debris covered an area of 78 ± 5 km 2 in 2000 (14 % of the glacierized area).
In 2000, glacier size ranged from 0.02 to 105 km 2 , with an average size of 3 km 2 and a median size of 0.9 km 2 (Table 4b). The histogram of glacier area (Fig. 4a)  the right (skewness = 8.4), showing that glaciers with area < 10 km 2 are predominant in this region, and glacier size decreases non-linearly. The long right-tail extremes represent only a few glaciers, such as Zemu with an area > 100 km 2 .
The average slope of all glaciers in the inventory was 23 • with a positive skew (skewness = 0.38) (Table 4b and Fig. 4b) and no significant differences among the four regions (p > 0.05) (Table 4b). Glacier length ranged from 0.08 to 23 km (Zemu glacier), with an average of 2 km (Table 4b and Fig. 4c). Glacier thickness ranged from 3 to 144 m, with the highest frequency for thicknesses less than 30 m (Fig. 4d). The frequency distribution of both glacier length and thickness were positively skewed with long tails, indicating the prevalence of short, shallow valley-type glaciers. Glacier aspect shows two predominant orientations: westnorthwest (W-NW) and east-northeast (E-NE), following the topographic divide (Fig. 5). On average, glaciers on the Nepal side had an average aspect of 237 • (SW), whereas glaciers on the Sikkim side had an average aspect of 131 • (SE), consistent with local topography (Table 4b).
Glacier terminus elevations in the Landsat/ASTER domain ranged from 3990 to 5777 m, with a mean of 4908 m; median glacier elevation ranged from 4515 to 6388 m, with a mean of 5702 m (Table 4b). Considering glacier median elevation as a coarse approximation of glacier equilibrium line altitude (ELA), our results are in agreement with Benn and Owen (2005), who documented higher ELAs on the northern slopes of the Himalaya (6000-6200 m) compared to ELAs on the southern slopes (4600-5600 m).
The Cryosphere, 9, 505-523, 2015 www.the-cryosphere.net/9/505/2015/ On a glacier-by-glacier basis, glaciers in the Corona domain lost 2-95 % of their area, with a mean of 32 % from 1962 to 2000 (Fig. 6). The spatial distribution of these area changes, illustrated in Fig. 6, shows that the largest area changes (> 70 % area loss) occurred for only a few isolated glaciers in the northern and southern extremities of the study area (17 glaciers). A closer examination of these glaciers revealed that these were small clean glaciers (< 0.1 km 2 ) with steep slopes (mean of 26 • ), indicating a need to investigate the topographic controls on area change and clean versus debris-covered glaciers separately.
Clean glaciers lost more of their area from 1962 to 2000 (34 %) compared to debris-covered glaciers (22 %) across the region, with few differences from east to west (Nepal and Sikkim) ( Table 6). The difference in mean rates of area change between clean and debris-covered glaciers was statistically significant based on a two-sample F test (p < 0.05). Figure 7a and b show a larger spread and a higher percentage of surface area loss of clean glaciers compared to debriscovered glaciers. For both glacier types, however, there is a high variability in percent area change, perhaps due to other factors such as local topography.
Linear regression analysis showed that percent area change per glacier was negatively correlated to glacier area, altitudinal range, and glacier median and maximum elevation (significant correlations at 99 % confidence interval, p < 0.01) ( Table 7). Glacier minimum elevation and slope were significant positive controls on glacier area change at 95 % confidence interval (p < 0.05). Solar radiation, precipitation and percent debris were not statistically significant controls on glacier area change (p > 0.1, confidence interval 90 %) ( Table 7). These are discussed in Sect. 5.3. A further analysis of the clean and debris-covered glaciers showed significant differences in terms of glacier area, area change, minimum elevation, altitudinal range and length based on a two-sample F test for variances) (p < 0.05) (Table 8). Clean glaciers in this area are ∼ 12 times smaller (1 km 2 on average) than debris-covered glaciers (15 km 2 ) and they have higher termini elevations (+391 m) and an overall altitudinal range about 3 times smaller than debriscovered glaciers (Table 7). On a glacier-by-glacier basis, clean glaciers lost more area (34 %) than debris-covered glaciers (22 %) from 1962 to 2000. Clean glaciers with smaller altitudinal range tend to display more area loss compared to debris-covered glaciers.

Spatial distribution of glacier characteristics across the study area
One of the important steps in utilizing our glacier inventory data is to understand spatial patterns in glacier characteristics across the region. Our study area displays region-wide consistency in glacier characteristics, notably glacier area, elevation and topography across the four subregions based on the 2000 glacier data (Table 4). For example, the prevalence of small glaciers noted in this area is consistent with worldwide patterns also observed for the Cordillera Blanca of Peru in a previous study (Racoviteanu et al., 2008a). There is however variability in glacier size within eastern Himalaya: for example, the mean glacier size reported in this study area (3 km 2 ) is double compared to the Khumbu region west of our study area (1.4 km 2 ) (Bajracharya and Shrestha, 2011). The glacier slope across our study area (23 • ) is consistent with average glacier slopes reported for the Khumbu region in Nepal (22 • ) (Bajracharya and Shrestha, 2011;Salerno et al., 2008), indicating a general tendency for steep glaciers across the region. There are only a few large, long glaciers in the area, such as Zemu glacier (103 km 2 , 23 km in 2000). With respect to glacier aspect, we also note similar predominant orientations of glaciers southwards, in the direction of the prevailing monsoon circulation consistent with other studies such as the Khumbu region (average aspect 181 • ) (Mool et al., 2002;Salerno et al., 2008). The comparison of glacier characteristics across subregions points to a pronounced gradient from north to south (Bhutan/China subregions compared to Sikkim/Nepal), particularly with respect to glacier elevations and debris cover. Glaciers on the northern side of the divide (China) have higher glacier termini and median elevations compared to the southern side (Nepal and Sikkim) (+700 and +400 m, respectively) (Table 4). These differences seem to be consistent with general air circulation patterns in the area. The Asian summer monsoon brings large amounts of precipitation on the southern slopes of the Himalaya, favoring glacier growth at lower elevations and a lower ELA. In contrast, in the upper reaches of the valleys and on the Tibetan plateau, the monsoon is blocked by the topographic barrier (Clift and Plumb, 2008), causing a drier climate and higher glacier ELAs. There is a much less pronounced east-west gradient in glacier elevations, with higher glacier minimum and me- dian elevations on the western side (Nepal) (+50 m) compared to the eastern side (Sikkim). This may be explained by the location of Nepalese glaciers on the western side of the topographic divide, away from the prevailing monsoon. Debris coverage also shows a pronounced variability north to south of the topographic divide. Himalayan glaciers are often referred to as "heavily" debris-covered, but the percent glacierized area covered by supraglacial debris varies across the mountain range. In our study area, debris cover is more prevalent on the southern side of the divide (Sikkim, 14 % of glacierized area) compared to the northern one (China, 2 % of the glacierized area), perhaps due to different geologic and topographic patterns. The northern side of the divide, which is part of the Tibetan plateau, is situated in a monsoon shadow and is therefore dry; the gentler slopes induce lower rates of erosion. In contrast, the southern slopes of the Himalaya tend to be heavily covered with debris due to the abundance of rock material from the steep slopes. The steep slopes made of soft sedimentary rocks and Precambrian crystalline rocks (Mool et al., 2002) and are prone to high rates of erosion, particularly with large amounts of monsoon moisture. This north-south difference in debris cover amount was also noted in other studies (Scherler et al., 2011). In our study, we found a lower percent of debris coverage (21 %) than in the entire central/eastern Himalaya reported in Scherler et al. (2011) (36 % debris cover), or in the Khumbu region, west of our study area, reported by Fujii and Higuchi (1977), Nuimura et al. (2012) (34.8 %), Racoviteanu et al. (2013a) (27 %) and Thakuri et al. (2014) (32 %).

Regional glacier area changes
The overall rate of surface area loss of 0.5 ± 0.2 % yr −1 from 1962 to 2000 for Sikkim and eastern Nepal obtained here is in agreement with other studies from the southern slopes of the Himalaya. Similar rates of area loss (0.1-0.3 % yr −1 ) were reported from the Khumbu and Garhwal regions, west of our study area, for approximately the same time period (Thakuri et al., 2014;Basnett et al., 2013;Nuimura et al., 2012;Bolch et al., 2008a;Bhambri et al., 2011). Similarly, for glaciers of Bhutan, east of our study area, Karma et al. (2003) found an average surface glacier area loss of 0.3 % yr −1 from 1963 to 1993. It is worth mentioning that these rates of area change are lower than those previously reported for the drier monsoon-transition zone in the western Himalaya (Himachal Pradesh) (0.7 % yr −1 ) by Kulkarni et al. (2007). In a more recent study, Bahuguna et al. (2014) found lower rates of glacier area loss (0.4 % yr −1 ) for the same area, which is in agreement with rates of area loss we report here for eastern part. Updated glacier area changes from other recent studies (Racoviteanu et al., 2014;Bolch et al., 2012) also point at lower rates of area loss than previously reported, particularly for the Indian Himalaya. The similar overall rate of glacier area change in the eastern compared to the western part for both debris-covered glaciers and clean glacier types suggests consistent patterns across the region ( Table 5).
The smaller glacier area loss for debris-covered glaciers noted in our study is in agreement with studies from Khumbu (Nuimura et al., 2012;Thakuri et al., 2014) or other studies in the central/eastern Himalaya (Bolch et al., 2008a;Thakuri et al., 2014;Bhambri et al., 2011). These studies also reported lower rates of glacier surface area loss and even stable or less retreating glacier termini for debris-covered tongues compared to clean glaciers (Scherler et al., 2011). Area changes for debris-covered glaciers need to be interpreted with caution due to the wide variability in debris cover characteristics such as thickness. Furthermore, these stagnating or less changing tongues may not reflect the true state of the glaciers, for example patterns of glacier thinning which may occur at similar rates to clean glaciers (Gardelle et al., 2012a;Kääb et al., 2012).

Topographic and climatic controls on area changes
While the consistent area change patterns across the subregions (east to west) are useful for comparison with larger areas, these patterns cannot be used to understand glacierby-glacier variability in area changes, which may be controlled by local topography and climate. In this study, we found that topographic factors, notably glacier size, glacier altitude (maximum, median, altitudinal range) and aspect, were most important in determining rates of glacier area loss in spatial domain 1. Glacier size plays a significant role in determining area change, i.e., smaller glaciers experienced higher rates of area loss ( Table 7). The tendency of larger glaciers to lose less area (> 20 km 2 ) was observed in various studies (Racoviteanu et al., 2008a;Salerno et al., 2008;Loibl et al., 2014), though in the case of Salerno et al. (2008) the correlation was not statistically significant for the Khumbu.
Higher glacier elevations and larger altitudinal ranges significantly reduce the rates of area loss, as was also noted in the Khumbu region in Nepal and elsewhere. The dependency of area change on glacier size and elevation is also consistent with observations from the Cordillera Blanca of Peru (Racoviteanu et al., 2008a) in the outer tropics, indicating consistent patterns in glacier area changes worldwide. Glacier slope also plays a significant role in determining glacier area change, i.e., the steeper the glacier, the larger the area loss observed in our study. The same tendency was observed in the Khumbu area (Salerno et al., 2008), but the correlation is less significant than the glacier altitude (p < 0.05). The presence of gentle slopes covered with supraglacial debris in the ablation areas of glaciers, fairly common in this area, may have reduced the strength of the correlation. Glacier aspect was also found to be a significant control on area change, with more area loss for glaciers oriented southwards and southwest (p < 0.01). This is in agreement with findings from Salerno et al. (2008) for the Khumbu region but is in contrast with results from Loibl et al. (2014) for the Nyainqêntanglha Range in southeastern Tibet, about 600 km east of our study area, who found that south-facing glaciers experienced lower rates of terminus retreat. Percent debris cover was a negative control on area change, i.e., glaciers with more extensive debris cover on their areas tend to lose less area overall, but this was not statistically significant (p > 0.05). Debris-covered glaciers may benefit from the insulating effect of debris cover above a certain "critical" debris thickness (Mihalcea et al., 2008a;Zhang et al., 2011), which needs to be further investigated.
Geographic location (latitude and longitude) were negative controls on glacier area change, suggesting that glaciers located north and eastwards of the study area tend to lose more area, but only latitude was statistically significant (p < 0.05). Climate indices (precipitation and solar radiation) were not significant factors controlling glacier area loss. In contrast, Loibl et al. (2014) showed that glaciers located in a monsoon-influenced area were more sensitive to cli- mate change. This is in agreement with larger-scale studies (Gardelle et al., 2013) that indicated a tendency for enhanced glacier wastage in the eastern, monsoon-influenced parts of the Himalaya. With respect to climatic factors in this area, Basnett et al. (2013) reported an increase in mean annual temperature, more significantly in the winter (+2 • C yr −1 in the last four decades). Increasing temperatures on the south slopes of the Himalayas were also noted in other studies (Shrestha et al., 2000;Thakuri et al., 2014) based on instrumental data but were estimated to have less effect on glacier area than changes in precipitation because of the orientation of these glaciers towards the prevailing monsoon circulation. In our study, the climatic control on glacier area is not conclusive, and finer-resolution, more accurate temperature and precipitation data sets would be needed. Furthermore, similarly to areas further east (Loibl et al., 2014), average annual solar radiation and latitude were not found to be significant controls on glacier area change in our study.

Surface temperature distribution on debris-covered tongues
Lower rates of area change for the debris-covered glaciers may be further explained by surface characteristics of debris cover (thickness, thermal conductivity and resistance). However, these are not easily available in this area due to the lack of field-based measurements and the difficulty of conducting surveys on the debris-covered tongues. Therefore, in this study, we are only qualitatively showing the distribution of surface temperature on selected debris-covered tongues in spatial domain 3 based on the 2002 ASTER scene. Figure 8 shows a high variability in supraglacial surface temperature at 90 m spatial resolution, but there is no clear general temperature trend for the eastern slopes (Sikkim side) versus the western slopes (Nepal) side. The fluctuations in surface tem- The Cryosphere, 9, 505-523, 2015 www.the-cryosphere.net/9/505/2015/ peratures along transects are clearly visible in Fig. 9, with some sharp spikes of high and low temperatures, particularly for Kanchenjunga and Yalung glaciers (labeled "A" and "B" in Fig. 8). This strong variability in supraglacial debris temperatures may be due to the presence of surface features such as debris thickness, size of the debris particles and thermal resistance and conductivity of the debris. For the debriscovered tongues investigated here, the supraglacial temperatures range from 0 to 30 • C, suggesting that the supraglacial debris heats up considerably during the day. At the glacier scale, temperature drops over supraglacial features such as ice walls and supraglacial lakes, which tend to be colder than the surrounding debris, and this is visible even at the coarse spatial resolution of the temperature data (90 m). In Fig. 9 we note the slight upward trend of the supraglacial temperature towards the glacier termini, particularly for Zemu glacier ("C" in Fig. 8). For this glacier, the upper-middle part of the debris surface is colder (−3 to 5 • C) than the last 10 km towards the glacier terminus (5 to 14 • C) (Fig. 9). In a different paper (Racoviteanu and Williams, 2012), we found similar patterns of surface temperature increasing towards the glacier terminus for the same glacier but based on a different ASTER scene (November 2001), indicating consistent patterns for this glacier. The higher surface temperatures towards the glacier terminus may indicate a thicker debris cover that insulates the ice underneath, noted in other studies (Mihalcea et al., 2008a). The daytime debris temperature ranges and strong spatial variability noted here are similar to the ones we found for Khumbu, west of this study area (−3 to 17 • C) (Racoviteanu et al., 2013b). In Khumbu, we found that supraglacial debris had a distinct temperature signal compared to other surfaces such as non-ice moraine, clean ice and supraglacial/proglacial lakes, with more pronounced differences among these three during the daytime. The suitability of ASTER-based surface for inferring debris characteristics, most notably thickness, has been demonstrated in other studies (Mihalcea et al., 2008a;Zhang et al., 2011). For this study area, there were no field measurements available to test the validity of ASTER temperatures for quantifying supraglacial debris characteristics. However, in a different study (Racoviteanu et al., 2013b), we validated ASTER-based surface temperatures extracted from nine night scenes from 2010 to 2011 for the Khumbu by inverting field-based long-wave radiation (L out ) using the Stefan-Boltzmann law (L out = εT 4 ). The measurements were from the automatic weather station installed on Changri Nup glacier (Wagnon et al., 2013). We found a good agreement between ASTER temperatures and field-based measurements (R 2 = 0.92) using a sensitivity analysis (ε = 0.97 ± 0.1) to account for small-scale variability in emissivity. Given that the Kanchenjunga-Sikkim area has similar characteristics to Khumbu in terms of debris cover and geographic location and that the images were acquired around the same time of the year as the Khumbu (November-January), we consider Figure 9. Surface temperature distribution along longitudinal transects from selected glaciers (shown in Fig. 8): A is Kanchenjunga glacier, B is Yalung glacier and C is Zemu glacier. Distance on the x axis is measured from the upper part of the debris-covered area down glacier to the terminus. that this validation may be applicable to the present study area.

The role of glacier lakes
The role of supraglacial/proglacial lakes for glacier area change in this area of the Himalaya was addressed in detail in recent studies (Basnett et al., 2013;Bajracharya et al., 2014). Gardelle et al. (2011) also pointed out the increased formation of supraglacial lakes particularly in the eastern part of the Himalaya. A quantitative assessment of lake formation is beyond the scope of this paper; here we only illustrate qualitatively some of the changes occurring on glaciers with supraglacial or proglacial lakes using high-resolution Corona and QuickBird imagery. and Bajracharya (2003) inventoried 266 glacier lakes covering a total area of 20 km 2 (3.5 % of the glacierized area) based on 2000 Landsat ETM+ imagery. For spatial domain 3, we estimated that glacier lakes covered 1.3 % of the total debris-covered glacier area, or 5.8 % of the area if we consider only the debris-cover (ablation) part, based on the QB/WV2 imagery. Salerno et al. (2012) reported similar percentage for the area of supraglacial lakes, i.e., 0.3-2 % of the overall glacierized area for the Khumbu region. While supraglacial lakes do not cover extensive areas of the glacierized surface, they were shown to increase surface ablation rates in this part of the Himalaya ;  Sakai et al., 2002). It was also shown that supraglacial lakes located at the glacier terminus tend to merge to create large, fast-growing proglacial lakes that accelerate glacier area loss (Basnett et al., 2013;Bajracharya et al., 2014). Some of the proglacial lakes in our study area are visible in Figs. 10 and 11 for the northern part of spatial domain 3 (Changsang, east Langpo, Jongsang, middle Lhonak, south Lhonak). Most of these lakes are moraine-dammed lakes considered dangerous for potentially inducing glacier lake outburst flood events and were shown to accelerate the glacier area loss in recent decades . Figure 10a-b show the evolution of the proglacial lake on North and South Lhonak glaciers in Sikkim, also noted in Basnett et al. (2013). A closer look at a subset area (Fig. 11) shows the visible growth of a proglacial lake for the adjacent North Lhonak and South Lhonak glaciers. We estimate that these two glaciers retreated ∼ 650 m and 1.3 km from 1962 to 2006, respectively. Another branch of the North Lhonak glacier has wasted significantly by ∼ 1.5 km from 1962 to 2006, and a glacier outlet is now clearly visible. The northern branch of Jongsang glacier was entirely covered by a supraglacial lake in 2006, while another part shows less significant rates of terminus retreat (∼ 100 m in 44 years). A part of the Jongsang glacier shows a slight "false" glacier tongue advance, most likely due to uncertainties in the mapping of Corona imagery. While our purpose here is not to present glacier length changes, we note that these estimates are in agreement with trends of glacier thinning and increased glacier lake formation reported in this area of the Himalaya previously (Gardelle et al., 2011;Basnett et al., 2013;Kääb et al., 2012).

Uncertainty and limitations
Inconsistencies in glacier area change estimates have been pointed out in other studies, for the Himalaya and elsewhere (Racoviteanu et al., 2008a, b), and are also noted in the current study. Glacier area changes in the Himalaya are heterogeneous and depend on a variety of factors including local topography and climate, so some caution should be applied when comparing rates of area changes from one area to other areas, even in the same climatic zone. For example, for Sikkim we estimated a surface area loss of 88.9 ± 5 km 2 (13.5 % from 1962 to 2006, or 0.36 ± 0.17 % yr −1 ). Other studies in this area point to contrasting results. For the same geographic area, Basnett et al. (2013) reported an area loss of 0.16 ± 0.10 % yr −1 from 1989/1990 to 2009/2010), which is about half of the area change in our findings. In contrast, a recent study  reported the highest rates of area loss (about 0.8 % yr −1 ) for the last decade, even higher than rates reported previously for the western Himalaya by Kulkarni et al. (2007). We speculate that such large differences might be due to errors inherent in the baseline data sets coupled with misclassification due to snow cover or debris-covered areas. Glacier area changes reported for Sikkim in different studies, using a variety of data including topographic maps (Table 9), illustrate this point. For example, for Sikkim, our study estimated 569 ± 70 km 2 of glacierized area in 2000 based on Landsat/ASTER data. For the same time period, Mool et al. (2002) reported an area of 577 km 2 based on the same source imagery (Landsat ETM+) (Table 9). These two area estimates differ only by 8.2 km 2 (1.4 %) of our estimated area and only the number of glacier differs substantially (186 glaciers in our study compared to 285 glaciers in ICIMOD study) most likely due to the way in which ice masses were split and how glaciers were counted. Methodology differences and inconsistencies in glacier estimates are quite common in multi-temporal image analysis performed by different analysts and were previously noted in other areas of the world (Racoviteanu et al., 2009). Similarly, for the 1962 decade, our analysis of Corona 1962 imagery for Sikkim yielded 178 glaciers with an area of 658 ± 20 km 2 . In a recent publication (Racoviteanu et al., 2014), we reported 158 glaciers with an area of 742 km 2 for the 1960s based on the Swiss topographic map. The GSI (Sangewar and Shukla, 2009) reported 449 glaciers with an area of 706 km 2 for the 1970s based on topographic maps. Our 1962 Corona glacier inventory yields a smaller total glacier area than the one based on the topographic map (84 km 2 , or 11 %) (Racoviteanu et al., 2014) or the GSI inventory based on topographic maps (48.3 km 2 , or 7 % area) (Sangewar and Shukla, 2009). We consider that both of these latter mentioned studies overestimated the glacier area in the 1960s, perhaps due to the presence of persistent snow in the source aerial imagery. We consider the Corona 1962 data set to be more reliable than the inventories based on topographic maps, and hence we used this data set as baseline for comparison with the recent imagery.
Glacier inventories in Sikkim for the recent decades also point to contradictory patterns. For the 1980s, another study (Kulkarni, 1992) reported a glacierized area of 431 km 2 for 1987/1989 based on Indian IRS-1A and Landsat data. Considering our 1962 Corona inventory, this would imply an area loss of 42 % since 1962 (2.1 % yr −1 ) followed by a strong increase in glacier area (+33.5 %, or +3 % yr −1 ) from 1987/1988 to 2000 (based on our Landsat analysis), which is undocumented and unlikely in this area. We consider the 1987/1989 estimates to be highly unreliable, given that there are no glacier surges that might induce an apparent "glacier growth". In some areas, we noted omissions of some debriscovered tongues from the glacier maps, which might explain some of the differences.

Summary and outlook
In this study we combined remote sensing data from various sensors to construct a new glacier inventory for the Kanchenjunga-Sikkim region in the eastern Himalaya. Based on 1962 Corona and 2006 QuickBird imagery, we found an overall negative glacier surface area loss of 0.5 ± 0.2 % yr −1 since 1962, in agreement with those noted in other studies in the Himalaya. The area loss rates reported here are lower than the average rate of 0.7 % yr −1 reported in other glacierized areas of the world such as the Alps (Kääb et al., 2002), the Tien Shan (Bolch, 2007) and the Peruvian Andes (Racoviteanu et al., 2008a). Glaciers exhibit heterogeneous patterns of area change depending on topographic and climatic factors, most notably glacier altitude (maximum, median, altitudinal range), glacier size, slope and aspect. Glacier area changes depend strongly on glacier size and elevation, which is consistent with other areas in the central/eastern Himalaya (Thakuri et al., 2014) or elsewhere, for example in the outer tropics (Racoviteanu et al., 2008a). The conclusions drawn with respect to spatial patterns in glacier characteristics, glacier area loss and their topographic and climatic dependency are as follows: -strong north-south gradient in terms of glacier elevations and debris cover, with larger percent of area covered by debris and higher glacier elevations on the northern side for the divide but less east-west gradient in these characteristics; -supraglacial debris cover is prevalent on the southern slopes of the Himalaya (14 % of the glacierized area) compared to northern slopes (2 %); -supraglacial lakes constitute about 6 % of the debriscovered area, and some of these supraglacial lakes have merged to form proglacial lakes; www.the-cryosphere.net/9/505/2015/ The Cryosphere, 9, 505-523, 2015 -higher rates of area loss for clean glaciers (34 %, or 0.7 % yr −1 ) compared to debris-covered glaciers (14.3 % or 0.3 % yr −1 ) across the subregions on a glacier-by-glacier basis; -the amount of glacier area loss is partly controlled by a glacier's headwater elevation, altitudinal range, glacier area, slope and aspect, with the largest area loss observed for small, steep glaciers with a smaller altitudinal range and less debris cover; -while Himalayan glaciers are undoubtedly undergoing negative area change, the rates of area loss noted in this study (0.5 % yr −1 ) as well as other recent studies in the area (0.2-0.4 % yr −1 since the 1960s) are lower than other glacierized areas worldwide (0.7 % yr −1 ).
The glacier area change estimates reported here are subject to uncertainties most notably with respect to early topographic maps and declassified Corona imagery; therefore a considerable effort was given to minimizing errors by multiple re-iterations of the glacier outlines. The understanding of the spatial patterns of glacier changes in the current study is limited by (1) a lack of a baseline elevation data set for the 1960 to compute glacier elevation changes from 1960s to 2000 and (2) lack of field-based measurements to validate debris-cover mapping and surface temperature distribution. With respect to the latter, while surface temperature trends show a slight increase towards the terminus, suggesting a thicker debris cover, the supraglacial surface temperatures are highly heterogonous and require additional investigation. A further improvement in the current study will be to include the supraglacial and proglacial lakes and surface temperature as determinant factors for the glacier area change, per-haps in a more sophisticated multivariate regression model. The glacier data sets constructed in this study can be further utilized to understand the behavior of glaciers in this little-investigated area of the Himalaya, particularly with respect to spatial patterns of glacier melt, and the contribution of glaciers to water resources.