Spatial patterns in glacier area and elevation changes from 1962 to 2006 in the monsoon-influenced eastern Himalaya

Introduction Conclusions References


Introduction
Himalayan glaciers have aroused a lot of concern in the last few years, particularly with respect to glacier area and elevation changes and their consequences for the regional water cycle (Immerzeel et al., 2010(Immerzeel et al., , 2012;;Kaser et al., 2010;Racoviteanu et al., 2013).Remote sensing techniques helped improve estimates of glacier area changes (Bajracharya et al., 2007;Bolch, 2007;Bolch et al., 2008a;Bhambri et al., 2010;Kamp et al., 2011), glacier lake changes (Wessels et al., 2002;Bajracharya et al., 2007;Bolch et al., 2008b;Gardelle et al., 2011) and region-wide glacier mass balance (Berthier et al., 2007;Bolch et al., 2011;Kääb et al., 2012;Gardelle et al., 2013), but significant gaps do remain.The new global Randolph inventory (Pfeffer et al., 2014) provides a global dataset of glacier outlines intended for large-scale studies; in some regions the quality varies and the outlines may not be suitable for detailed regional analysis of glacier parameters.Among recent glacier inventories, we cite those in the western part of the Himalaya (Bhambri et al., 2011;Kamp et al., 2011;Frey et al., 2012), and a few in the eastern Himalaya (Bahuguna, 2001;Krishna, 2005;Bajracharya and Shrestha, 2011;Basnett et al., 2013).Parts of the monsoon-influenced eastern Himalaya (the eastern extremity of Nepal, Sikkim and Bhutan) still lack comprehensive, multi-temporal glacier data needed for accurate change detection.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 satellite image acquisition.Furthermore, this area has very limited baseline topographic data needed for comparison with recent satellite-based data, as discussed in detail in Bhambri and Bolch (2009a).The earliest Indian glacier maps date from topographic surveys conducted by expeditions in the mid-nineteenth century (Mason, 1954), but these are limited to a few glaciers.The Geologic Survey of India (GSI) inventory based on Survey of India maps (Sangewar and Shukla, 2009) is not in the public domain.For eastern Nepal, 1970's topographic maps from Survey of India 1 : 63 000 scale are available, but their accuracy is not known with certainty.Given these limitations, de-Introduction

Conclusions References
Tables Figures

Back Close
Full classified Corona imagery from the 1960's and 1970's has increasingly been used to develop baseline glacier datasets, such as in the Tien Shan (Narama et al., 2007), Nepal Himalaya (Bolch et al., 2008a) and parts of Sikkim Himalaya (Raj et al., 2013).While some theoretical understanding of the climate-glacier relationship at a mountain-range scale is starting to emerge (Scherler et al., 2011;Bolch et al., 2012;Gardelle et al., 2013), updated, accurate glacier inventories are continuously needed to advance our understanding of the climatic, topographic and glaciological parameters that govern glacier fluctuations across the Himalaya.In this study, we evaluate spatial patterns in glacier area and elevation changes using multi-temporal satellite data from Corona KH4, Landsat 7 Enhanced Thematic Mapper Plus (ETM+), the Advanced Spaceborne Thermal Emission Radiometer (ASTER), QuickBird (QB) and WorldView-2 (WV2) sensors.We constructed a new "reference" geospatial glacier inventory based on 2000 Landsat ETM+ coupled with ASTER imagery, and extracted glacier parameters on a glacier-by-glacier basis.We evaluated spatial trends in glacier area and elevation changes over almost five decades based on high-resolution datasets (1962 Corona KH4 and2006 QB/WV data), and investigated them using statistical analysis.A particular emphasis was placed on the behavior of clean glaciers vs. debris-covered glaciers.For a subset of the debris-covered tongues, we computed glacier elevation changes (1960s to 2000) on the basis of topographic maps and recent DEMs and related them to ASTER-derived surface temperature and debris cover characteristics on a pixel-bypixel basis.In previous studies, Sakai et al. (1998Sakai et al. ( , 2002) ) and Fujita and Sakai (2014) illustrated the importance of supraglacial surface temperature as well as surface features (ice walls, ablation cones, supraglacial lakes) for inducing differential ablation on the surface of the debris-covered tongues.Here we complement these studies with remote sensing techniques to investigate the spatial variability of supraglacial debris, as a step towards inferring its thermal characteristics.Introduction

Conclusions References
Tables Figures

Back Close
Full

Study area
The study area encompasses glaciers in the eastern Himalaya (27 • 04 52 N to 28 • 08 26 N latitude and 88 • 00 57 E to 88 • 55 50 E longitude), located on either side of the border between Nepal and India in the Ganges and Brahmaputra basins (Fig. 1).
Relief in this area ranges from 300 m to 8598 m (Mt.Kanchenjunga).Long valley glaciers cover about 68 % of the glacierized area, mountain glaciers cover 28 %, and the remaining are 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, mapped in 1970 by Survey of India (Shanker, 2001;Sangewar and Shukla, 2009).The western part is located in eastern Nepal, and encompasses the Tamor and Arun basins.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 (Yanai et al., 1992;Benn and Owen, 1998).The Himalaya and Tibetan plateau (HTP) acts as a barrier to the monsoon winds, bringing about 77 % of precipitation on the south slopes of the Himalaya during the summer months (May to 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 per year, with annual averages of 3580 mm recorded at Gangtok station (1812 m) (1951 to 1980) (IMD, 1980), and 164 rainy days per year (Nandy et al., 2006).Mean minimum and maximum daily temperatures at this station were reported as 11.3 • C and 19.8 • C, with an average of 15.5 • C based on the same observation record (IMD, 1980).Introduction

Conclusions References
Tables Figures

Back Close
Full   -EROS 1996).The Corona KH4 system was equipped with two panoramic cameras (forward-looking and rear-looking 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.The nominal ground resolution reported for the KH4 mission is 7.62 m (Dashora et al., 2007); however, we calculated an actual pixel size of approximately 2 m using the scale of the photos and the scan resolution.Raw, unprocessed Corona images obtained from USGS are known to contain significant geometric distortions due cross-path panoramic scanning.Furthermore, the Frame Ephemeris Camera and Orbital Data (FECOD) 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 available.Orthorectification of Corona scenes is notoriously difficult due to the heavy distortions, so here we describe in detail the methodology used.
We defined a non-metric camera model in ERDAS Leica Photogrammetric Suite (LPS), with parameters (focal length, air photo scale, flight altitude) extracted from the declassified documentation of the KH4 mission (Dashora et al., 2007).The LPS bundle block adjustment procedure was used to estimate the orientation of all the CORONA Introduction

Conclusions References
Tables Figures

Back Close
Full stripes simultaneously, and model parameters were calculated on the basis of 117 ground control points (GCPs) extracted from the panchromatic band of the 2000 Landsat ETM+ image (15 m spatial resolution).GCPs and were identified on the Landsat image on non-glacierized terrain including moraines, river crossings, and outwash areas.Tie points (TPs) were digitized from one Corona strip to another as well as on the Landsat image.Elevations were extracted from the Shuttle Radar Topography Mission (SRTM) DEM v.4 (CGIAR-CSI 2004) and were used to correct effects of topographic terrain displacements.The Corona stripes were mosaicked in ERDAS LPS to produce the final orthorectified image.The horizontal accuracy (RMSEx, y) of the bundle block adjustment process was 10.5 m.An independent analysis of location accuracy using 30 check points chosen independently of the initial GCPs yielded an actual "ground" RM-SEx, y of the Corona images of ∼ 60 m.This is consistent with accuracies obtained on Corona images in Mexico, using a fitting software and the FECOD parameters (H.Snyder, personal communication, 2011).
(2) The orthorectified Landsat ETM+ scene from December 2000, obtained from the USGS Eros Data Center was used as "reference" dataset.The Landsat ETM+ scene has seven spectral channels at 30 m spatial resolution, a thermal channel at 60 m and a panchromatic channel at 15 m, a revisit time of 16 days and a large swath width (185 km).In addition, six orthorectified ASTER scenes from 2000 to 2005 were obtained at no cost through the Global Land Ice Monitoring from Space (GLIMS) project (Raup et al., 2007).The ASTER scenes have 3 channels in the visible wavelengths (15 m spatial resolution), 3 channels in the short-wave infrared at 30 m, and four thermal channels at 90 m, a swath width of 60 km and a revisit time of 16 days.Images were selected at the end of the ablation season for minimal snow, and had little or no clouds.
The ASTER scenes were used to complement the Landsat-based glacier delineation in challenging areas where shadows or clouds obstructed the view of the glaciers.In addition, the surface kinetic temperature product (AST08) from an ASTER scene from November 2001 was used in a debris cover delineation algorithm (Racoviteanu and Introduction

Conclusions References
Tables Figures

Back Close
Full Williams, 2012) and another ASTER scene from October 2002 was used to investigate surface temperature trends over debris covered tongues.
(3) Two QB scenes from January 2006 were obtained from Digital Globe as orthoready standard imagery (radiometrically calibrated and corrected for sensor and platform distortions) (Digital_Globe, 2007).These scenes cover an area of 1107 km 2 , and were well-contrasted and mostly snow-free outside glaciers.We orthorectified these scenes in ERDAS Imagine Leica Photogrammetry Suite (LPS) using Rational Polynomial Coefficients (RPCs) provided by Digital Globe, using an 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 to reduce disk space and processing time.One WorldView-2 (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 datasets were registered to UTM projection zone 45N, with elevations referenced to WGS84 datum.

Elevation datasets
Two elevation datasets were used in this study: (1) the hydrologically-sound CGIAR SRTM DEM (90 m spatial resolution) (CGIAR-CSI 2004) was used to extract glacier parameters for the 2000 decade.The SRTM accuracy and biases have been quantified in several studies (Berthier et al., 2006;Fujita et al., 2008;Gardelle et al., 2012).Its vertical accuracy in our study area, calculated as root mean square (RMSEz) 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.
(2) The Swiss topographic map (1 : 150 000 scale), compiled from Survey of India maps from the 1960s to 1970s, published by the Swiss Foundation for Alpine Research was used to extract 1960's glacier elevation contours.
The exact date of each quadrant is not known with certainty because the original largescale Indian topographic maps at this scale are restricted within 100 km of the Indian border (Srikantia, 2000;Survey_of_India, 2005).Introduction

Conclusions References
Tables Figures

Back Close
Full

Analysis extents
We defined two spatial domains for our study area (Fig. 1).Spatial domain 1 includes the Sikkim province of India, eastern Nepal (Tamor and Arun basins), as well as parts of Bhutan and China (Table 2).This spatial domain was split into four sub-regions 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) are used to characterize the sub-regions climatically.The dataset 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;Andermann et al., 2011;Palazzi et al., 2013), here we are not concerned with the absolute values of gridded precipitation, but only with characterizing the sub-regions in our study area using relative rainfall values.TRMM data integrated over 10 years (1998 to 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 mm yr −1 versus 805 mm yr −1 ).There is a more pronounced north-south gradient in precipitation, with the lowest amount of precipitation on the China side (146 mm yr −1 ) (Table 3).
Spatial domain 2 is a subset of spatial domain 1, and comprises 50 glaciers from Nepal (Tamor basin) and Sikkim (Zemu basin), clean as well as and debris-covered.These glaciers were used for a more thorough analysis of glacier area and elevation changes and their dependence on climatic and topographic variables using correlations and multiple regression analysis for two time steps: the 1960's decade (represented by Corona imagery), and 2010's decade (represented by QB and WV2).Introduction

Conclusions References
Tables Figures

Back Close
Full

Glacier delineation and analysis
For the 1960's decade, glacier outlines were extracted from the panchromatic Corona imagery by thresholding the digital numbers (DN > 200 = snow/ice) based on visual interpretation.A 5 × 5 median filter was used to partially remove remaining noise, as shown in other studies (Paul, 2007;Racoviteanu et al., 2008a).Manual corrections were applied subsequently on the basis of the Swiss topographic map using on-screen digitizing in areas of poor contrast or transient snow/clouds, which obstructed the view of glaciers.
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 1960's 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 combined with and topographic variables in a decision tree, described in Racoviteanu and Williams (2012).
For the 2006/09 QB/WV2 images, ice masses 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 dataset were delineated manually on the basis of supraglacial features (lakes and ice walls) visible on the high Introduction

Conclusions References
Tables Figures

Back Close
Full resolution images.We also mapped supraglacial lakes based on band ratios validated with texture analysis.For all inventories in spatial domain 1, 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 elevation, median elevation, average slope angle and average aspect were extracted on a glacier-by-glacier basis using zonal functions.Average glacier thickness were calculated from mass turnover principles and ice flow mechanics Huss and Farinotti (2011), based on the approach of Farinotti (2009).Their method uses glacier outlines and the SRTM DEM to derive thickness estimates iteratively based on Glenn'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 (Raup and Khalsa, 2007;Racoviteanu et al., 2009).
Glacier area changes (1962 to 2000) and their dependency on topographic and climatic variables were calculated on a glacier-by-glacier basis for the 50 glaciers in spatial domain 2. Elevation changes (1960 to 2000) were calculated for 21 debris-covered glacier tongues from 1960 to 2000 on the basis of the topographic map and the SRTM DEM.Elevation contours were digitized from the topographic map and were interpolated using the TopoGrid algorithm in ArcGIS, using a 90 m pixel size.We then computed elevation differences  over the debris covered tongues on a pixelby-pixel basis and examined their spatial patterns.ASTER-based surface temperatures from the 2002 ASTER image were extracted for each debris-cover tongue, and the relationship between elevation differences and surface temperature was evaluated along longitudinal transects on a pixel-by-pixel basis for the selected glaciers.
Uncertainties of glacier outlines derived using the semi-automated methods for Corona, Landsat/ASTER and Quickbird were ±3 %, ±6 and ±3 %, respectively.Uncertainties in elevation changes of the debris covered tongues were calculated as the Introduction

Conclusions References
Tables Figures

Back Close
Full 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 an area of 1463 km 2 ± 88 km 2 (Table 4a).Of the total glacierized area in spatial domain 1, a total of 160 km 2 ± 10 km 2 was covered by supraglacial debris (11 % of the glacierized area).Of the 487 glaciers in this spatial domain, 68 glaciers (13 %) had some percent of debris cover on their tongues.The debris cover distribution among the four regions in spatial domain 1 displays differences in the north and south directions.
In Sikkim, supraglacial debris covered an area of 78 km 2 ± 5 km 2 in 2000 (14 % of the glacierized area), in contrast with the northern side of the Himalaya (China, 2 % of the glacierized area).The prevalence of debris cover on the south side of the Himalaya can be explained in part by the geology and topographic patterns in the two regions.
The northern side of the divide is part of the Tibetan plateau, situated in a monsoon shadow, with gentler slope and lower rates of erosion because of the dry climate.The southern slopes of the Himalaya are steep, comprise of soft sedimentary rocks and Precambrian crystalline rocks (Mool et al., 2002).These slopes are prone to erosion and rock fall due to large amounts of moisture brought by the monsoon, which may explain the high amount of debris on the glacier surface on the south side of the divide.The behavior of debris covered vs. clean glaciers in relation to topographic/climatic patterns is explored in more detail in Sect.4.3 below.
In 2000, glacier size ranged from 0.05-105 km 2 , with an average size of 3 km 2 and a median size of 0.9 km 2 (Table 4b).These values for glacier area appear small, but are Introduction

Conclusions References
Tables Figures

Back Close
Full about three times larger than tropical glaciers for example, where the average glacier size was ∼ 1 km 2 (Racoviteanu et al., 2008a).The frequency histogram of glacier area (Fig. 4a) is skewed to the right (skewness = 8.4), showing that glaciers with area < 10 km 2 are predominant, and glacier size decreases non-linearly.The long right-tail extremes represent only a few glaciers with an area > 100 km 2 .The prevalence of small glaciers is a pattern faily common in inventories which include small glaciers, as we noted in the Cordillera Blanca of Peru in a previous study (Racoviteanu et al., 2008a).Glacier termini elevations in spatial domain 1 ranged from 3990 m to 5777 m, with an average of 4908 m; median glacier elevation ranged from 4515 m 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 ELA's on the northern slopes of the Himalaya (6000-6200 m) compared to ELA's the southern slopes (4600-5600 m), due to drier climate in the former.In our study area, we note a slight east-west gradient in glacier elevations, with higher glacier minimum and median elevations on the western side (Nepal) (+50 m) compared to the eastern side (Sikkim).This may be explained by the location of glaciers on the western side of the topographic divide, away from the monsoon.The north-south gradient in glacier elevations is stronger than the eastwest, with glacier termini and median elevations much higher on northern side of the divide (China) compared to the southern side (Nepal/Sikkim) (+700 m and +400 m respectively).These differences seem to be consistent with general air circulation patterns in the area, where the Asian summer monsoon brings large amounts of precipitation on the southern slopes of the Himalaya, favoring glacier growth at lower elevations.In contrast, the upper reaches of the valleys and the Tibetan plateau have a drier climate due to the monsoon being blocked by the topographic barrier (Clift and Plumb, 2008), causing glaciers to be found at higher elevations.
The average slope of glaciers in spatial domain 1 was 23 • , with a positive skew (skewness = 0.38) (Fig. 4b) and no significant differences among the four regions Introduction

Conclusions References
Tables Figures

Back Close
Full (p > 0.05) (Table 4b).Glacier length ranged from 0.08 km to 23 km (Zemu glacier), with an average of 2 km (Fig. 4c).Mool (2002) reported a length of 26 km for Zemu glacier in 1970's, and the Geologic Survey of India (Sangewar and Shukla, 2009) reported 28 km for the same glacier based on 1970 topographic maps.This suggests a decrease in length of ∼ 5 km in 30 years for this glacier.Glacier thickness ranged from 3 m to 144 m, with the highest frequency for thickness less than 30 m (Fig. 4d).The frequency distribution of both length and thickness were positively skewed with long tails, indicating the prevalence of short, shallow valley-type glaciers, respectively.Glacier aspect shows two predominant orientations of the glacier tongues: west-northwest (W-NW) and east-northeast (E-NE) (Fig. 5).Glaciers in Nepal had an average aspect of 237 • (SW), whereas glaciers of Sikkim for example had an average orientation of 131 • (SE).In a previous inventory, Mool et al. (2002) attributed the predominant orientation of glaciers (south, southwest, southeast and east) to the higher temperatures on the western slopes, which prevent glacier growth compared to the colder eastern slopes.

Glacier area changes 1962-2006 (spatial domain 2)
Out of the 487 glaciers in spatial domain 1, a sample of 50 glaciers from Tamor basin in Nepal and Zemu basin in Sikkim (spatial domain 2) were used for deriving area changes from 1962 (Corona imagery) to 2000 (Landsat/ASTER) and 2006 (Quickbird).
To minimize uncertainties due to differences in rock outcrops specific to each dataset, we used the same rock outcrops for each dataset.We obtained an area change of The area change results presented above (−0.36% yr −1 ± 0.17 % in the last four decaades) point to rates of retreat that are half the ones reported in other parts of the Himalaya.For example, an area change of −0.7 % yr −1 was reported by Kulkarni et al. (2007) for the western Himalaya and in other glacierized mountain ranges such as Alps (Kääb et al., 2002), the Tien Shan (Bolch, 2007) and the Peruvian Andes (Racoviteanu et al., 2008a).We speculate that eastern Himalaya glaciers may be less sensitive to climatic changes given their location in a monsoon-influenced area, which may provide enough precipitation to high altitudes to maintain accumulation.Basnett et al. (2013) examined annual climate records from the Gangtok station (1812 m), and found that mean annual temperatures increased at a rate of 0.04 • C yr −1 in the last four decades, with most warming observed during the winter, and a slight, non-significant decreasing trend in precipitation.However, given the location of the station well below the glacier termini, these data are not conclusive.
Correlations between topographic/climatic parameters and area change for the 50 glaciers in spatial domain 2 indicate that glacier maximum elevation, altitudinal range, median elevation, glacier area, percent debris cover, aspect and precipitation were correlated with the percent area change from 1962 to 2006 (correlation statistically significant at the 90 % confidence interval, p < 0.1) (Table 6).The strongest (negative) correlation with area change was found for maximum elevation and altitudinal range, indicating larger area changes for glaciers situated on lower summits, and with a small altitudinal range.In contrast, glacier location (latitude and longitude), slope, solar ra-Introduction

Conclusions References
Tables Figures

Back Close
Full diation and minimum elevations were not statistically significant, and they were only weakly correlated to area change (Table 6).On a glacier-by-glacier basis, glaciers lost −0.7 % to −64 % of their area, with a mean of 25 % (−0.56 % yr −1 ) from 1962 to 2006 (Fig. 6).Spatial patterns of these area changes (Fig. 6) show largest area changes (> 50 % area loss) for several glaciers in the northern and southern parts of the study area.A closer look at these glaciers showed that they are small (< 1 km 2 ), steep (average slope of 25 • ), have little debris cover (< 6 % of glacier area on average), and have termini elevations of 4423 to 5500 m.Correlation analysis suggests a greater glacier area loss for south-facing slope aspects.Furthermore, smaller glaciers situated at lower maximum elevations, with a smaller altitudinal range exhibited larger rates of area change.These results are consistent with observations from a different climatic area, the Cordillera Blanca of Peru (Racoviteanu et al., 2008a), indicating consistent patterns in glacier behavior.The relationship between the debris cover percent of glacier area and area loss was statistically significant at 95 % confidence interval (p < 0.05), i.e. glaciers with larger parts of their area covered by debris experienced less area loss.This is to be expected given the potential insulating effect of debris cover on glaciers above a certain "critical" debris thickness (Mihalcea et al., 2008;Zhang et al., 2011).There were differences in area changes between clean and debris covered glaciers.Figure 7a and b shows a larger spread and a higher percentage of area loss of clean glaciers compared to debriscovered glaciers.Between 1962 and 2006, clean glaciers lost 3.4 to 64 % of their area with a mean of 32.6 %, while debris-covered glaciers lost only 0.7 % to 61 % of their area, with a mean of 14.3 % on a glacier-by-glacier basis.The difference in mean rates of area change between clean and debris covered glaciers was statistically significant based on two-sample F test (p value < 0.05).
Topographic parameters of debris-covered vs. clean glaciers may explain some of these differences in area changes.Statistical analysis (two sample F test for variances) showed significant differences between the debris-covered and clean glaciers samples in terms of glacier area, area change, minimum elevation, altitudinal range and Introduction

Conclusions References
Tables Figures

Back Close
Full  7).Clean glaciers in this area are ∼ 10 times smaller (2 km 2 on average) than debris-covered glaciers (23 km 2 ), they have higher termini elevations (+535 m), and an elevation range about 3 times smaller than debris-covered glaciers (Table 7), which may all contribute to their higher rates of retreat.In Figs. 8 and 9 we illustrate a few of these changes in more detail, using high resolution Corona and QB imagery.Figure 8 a and b shows the rapid growth of the supra-glacier lake on S. Lonak glacier in Sikkim, also noted in Basnett et al. (2013).
Glaciers on the upper northern part of the image are most likely rock glaciers, where some of the ice underneath has collapsed, but the area changes are less visible, and more uncertain, as pointed above.A closer look at a subset area (Fig. 9) shows two adjacent glaciers (N.Lonak and S. Lonak) which underwent mass wasting.Their terminus retreated ∼ 650 m to 1.3 km respectively from 1962 to 2006, with visible growth of a supraglacial lake.Another branch of N. Lonak glacier has wasted significantly (∼ 1.5 km from 1962 to 2006), and a glacier outlet is now clearly visible.The terminus of the glacier in 2006 is uncertain, since some part of the glacier terminus may be inactive.In contrast, other debris-covered glaciers such as the near-by Jongsang glacier in the upper right part of the image does not show significant rates of terminus retreat (∼ 100 m in 44 years), in spite of increase in the area of supraglacial lakes.

Glacier elevation changes for debris-covered tongues 1960-2000
Given that some of the glaciers investigated here show little area change, the question remains whether some of these glaciers are experiencing elevation changes (thinning/thickening), and what may govern these changes (for example, debris cover or supra glacial lakes).In this area, supraglacial lakes covered only 5.8 % of the area of the debris-covered tongues in 2006/09, as delineated from QB/WV2 imagery.Some of the glaciers, for example Zemu, the largest glacier in this area, had 27 % of its area covered by supraglacial debris in 2000.Out of the 50 glaciers analyzed in spatial domain 2, 21 glaciers had debris cover in their ablation zones (an average of 21 % of the Introduction

Conclusions References
Tables Figures

Back Close
Full glacierized area on a glacier by glacier basis).We chose these 21 glacier tongues to evaluate elevation differences from 1960 to 2000 on a pixel-by-pixel basis (Fig. 10).We found that debris-covered tongues experienced an average lowering of −30.8 ± 39 m from 1960's to 2000 (−0.8 m ± 0.9 m yr −1 ), which are similar to annual rates reported in Gardelle et al. (2013) in the same climatic zone (Bhutan, −0.5 ± 0.13 yr −1 and Hengdun Shan, −0.81 ± 0.13 yr −1 ) for the period 2000-2011.Our rates are lower than surface lowering in the ablation zones of glacier further west (Karakoram, −0.33±0.16yr −1 ) (Gardelle et al., 2013).We obtained higher thinning rates for glaciers in contact with a growing pro-glacial lake, such as S. Lonak glacier discussed above ("A" in Fig. 10).Other glaciers such as Zemu ("B" on Fig. 10) display heterogeneity in the elevation changes, with less distinct spatial patterns.In Fig. 10 we note larger elevation differences in the mid-upper zone of the ablation area of debris-covered tongues, and generally thickening towards the terminus, for example Yalung glacier ("C") and Kanchenjunga glacier ("D").The latter shows an interesting pattern of "bulging" of debris in the mid-lower part of the debris-covered area, which we speculate might be due to a thickening wave that is reaching the glacier terminus.Due to the lack of velocity or mass balance measurements in this area, we cannot assert this.Here we consider that high rates of thinning of > 150 m, which are observed towards the rock walls in the upper, steep parts of the debris-covered area, are most likely due to errors in the topographic map in these areas.

Supraglacial surface temperature patterns
To explain some of these difference in thinning/thickening patterns, we extracted surface temperature along longitudinal transect for a few glaciers, as shown in Fig. 10.
Here we are presenting one of the glaciers, Zemu glacier ("B" in Fig. 10).Surface temperature patterns derived from 2002 ASTER imagery for Zemu glacier (Fig. 11) display a large variability, particularly for the 10 km towards the glacier terminus.As a general trends, we note higher surface temperatures from the terminus to ∼ 12 km upwards, indicating a thicker debris which insulates the ice underneath, then temper-Introduction

Conclusions References
Tables Figures

Back Close
Full atures start decreasing from about 12 km from the glacier terminus to the limit with clean ice (middle of the ablation zone), indicating a thinner debris cover.Elevation differences display a high variability, with spikes of large changes in areas around supraglacial lakes and ice cliffs, which in general accelerate ablation (Sakai, 2002;Fujita and Sakai, 2014).In a different paper (Racoviteanu and Williams, 2012), we found similar patterns of surface temperature increasing towards the glacier terminus, indicating the presence of thicker supraglacial debris.On Fig. 11, the relationship between elevation changes and surface temperature is more clear in the middle-upper debris area mentioned above, where we also note larger elevation differences (thinning).There is less variability of surface temperature than the lower part, probably associated with thin supra-glacial debris in this area.Regression analysis using surface temperature as explanatory variable for Zemu showed a non-significant dependency of elevation changes on surface temperature (p > 0.05).An ordinary least-squared regression using all 21 debris-covered tongues showed a weak dependency of elevation changes on surface temperature (R 2 = 0.01).Standard residuals are close to normally distributed, with larger residuals in the upper part of the debris covered tongues (close to steep walls), discussed above.

Uncertainty estimates
Glacier outlines derived from remote sensing data, even when using well-established semi-automated methods, are known to be subject to various degrees of uncertainty, as discussed in recent studies (Racoviteanu et al., 2009;Paul et al., 2013).This becomes an important issue in glacier change analysis, where errors from various data sources accumulate at each processing step.In this study, we combined remote sensing data at various spatial and temporal resolutions with data from topographic maps and DEMs to derive glacier parameters and estimate of area/elevation changes.The main sources of uncertainty here arise from: (1) errors inherent in the source data (geolocation errors, sensor limitations and ambiguity in the satellite signal, and/or accuracy of the topographic maps), (2) image classification errors (positional errors and/or errors Introduction

Conclusions References
Tables Figures

Back Close
Full due to the algorithms used for glacier mapping); (3) conceptual errors associated with the definition of a glacier, including mapping of ice divides, mixed pixels of snow and clouds, and internal rock all described in detail in Racoviteanu et al. (2009); (4) elevation errors due to inherent DEM uncertainty.Other types of errors, not discussed here, include errors due to atmospheric and topographic effects, resampling, image co-registration, spectral mixing or raster to vector conversion.Instead, we focus on the main types of errors, with an emphasis on classification errors and conceptual errors, which we consider key for glacier area change estimates.
1.The orthorectification process of the 1962 Corona, using GCPs from a Landsat scene, yielded an error of ±10 m.Geolocation errors on the ground were estimated using a trend analysis on the horizontal shifts between Corona and the reference Landsat scene, and resulted in ∼ 60 m, with the largest errors were concentrated towards the edges of the images, mostly outside the glaciers.We consider that these errors minimally affect our calculations of area change, since they do not impact the differences in area.The geolocation accuracy of the standard QB scenes is estimated as RMSE x , y of 14 m globally (Digital_Globe, 2007).
2. The accuracy of remote sensing glacier outlines is often estimated using the "Perkal epsilon band" around glacier outlines (Racoviteanu et al., 2009;Bolch et al., 2010), using a ∼ 1-pixel variability (Congalton, 1991).Recent glacier analysis comparison experiments estimated the accuracy of automated glacier outlines using manually-derived outlines at multiple times and various analysts (Raup et al., 2007;Paul et al., 2013), and obtained a range of uncertainty of < 5 % for remote sensing glacier outlines compared to high resolution imagery.In this study, we could not perform a validation using high-resolution imagery, since this was used as source imagery for area change.Instead, using the 1-pixel variability (±30 m for Landsat/ASTER, ±6 m for Corona and ±3 m for QB outlines), we obtained a total estimate of image classification errors of 3-6 % for our datasets.Highest uncertainties were associated with the larger pixel size (Landsat).The

Conclusions References
Tables Figures

Back Close
Full method is known to slightly over-estimate the errors described in Burrough and McDonnel (1998), so we consider our overlal area accuracy estimates to be rather conservative, and accommodate other errors not considered here (higher uncertainty of debris-covered tongues, GIS operations etc.).For manually-adjusted glacier outlines, including some of the debris-covered tongues, we minimized errors by using screen digitizing in streaming mode, with a high density of vertices.
3. Conceptual errors are relevant here in the context of area change and relate to how the glaciers were defined.The glaciologic community has come to some consensus on how these rules should be applied (Racoviteanu et al., 2009), and here we comply to these guidelines.For consistency, we removed internal rocks from all the area calculations; we manually removed perennial snowfields from the glaciers; we included "inactive" bodies of ice above the bergschrund as part of a glacier.Uncertainties in area change were computed as the RMSE of the uncertainties embedded in each dataset, from classification errors above, and amounted to 3-6 % of the glacierized area.Uncertainties related to differences in internal rocks in each dataset were derived by comparing area changes computed with internal rocks specific to each dataset, vs. "merged" internal rocks from each datasets.The differences in glacier datasets due to rock inconsistencies amounted to ∼ 2 % glacier area, inducing differences in the rates of change in area of ±0.05 % yr −1 .For simplicity, here we neglected the area change that might be due to exposure of new internal rock due to glacier ice thinning.The glacier area changes computed from 1962 to 2006 suggest a higher rate of retreat in the last decade (−0.43 % yr −1 ± 0.9 % from 2000 to 2006) compared to the previous period (−0.20 % yr −1 ± 0.16 % from 1962 to 2000).We speculate that some of this apparent "accelerated" rate of glacier retreat in the last decade might be due to the difference in spatial resolution of the imagery.Due to the short time period (2000)(2001)(2002)(2003)(2004)(2005)(2006), area uncertainties for this time step might be larger than the area change.

TCD Introduction Conclusions References
Tables Figures

Back Close
Full 4. Elevation errors in the glacier vertical changes were calculated using the RMSE z , assuming an error of 25 m for the 1960's topographic map (1/4 of the contour interval) and the calculated error for SRTM of 31 m, yielding an elevation change of 30.5 m ± 39 m.The errors associated with digitizing of the topographic map are considered smaller than errors associated with the original map (Paul, 2003).

Comparison with other areas
Comparison of our results to other studies in this area of Himalaya is subject to inconsistencies and inaccuracies in the classification methods used.For example, in Sikkim, our study estimated 569 km 2 ± 70 km 2 of glacierized area for Sikkim in 2000 period based on Landsat/ASTER data (Sect.4.1).For the same time period, Mool et al. (2002) reported 285 glaciers with an area of 577 km 2 based on the same source imagery (Landsat ETM+) (Table 8).These two area estimates differ only by 8.2 km 2 (1.4 %) of our estimated area.The significant difference in glacier numbers (186 glaciers in our study compared to 285 glaciers in ICIMOD study) are 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 km 2 ±20 km 2 .
In a recent publication (Racoviteanu et al., 2014), we reported 158 glaciers with an area of 742 km 2 in 1960s based on the Swiss topographic map, which is 13 % higher than our new estimates based on Corona imagery from 1962.We consider the Corona estimates more reliable than the topographic map, and we use this dataset as the baseline for comparison with more recent imagery.Our 1962 Corona glacier inventory differs by 48.3 km 2 , or ∼ 7 % from the inventory published by the Geological Survey of India (GSI) (Sangewar and Shukla, 2009), which reported 449 glaciers with an area of 706 km 2 for approximately the same time period.We consider both the GSI and the Swiss topographic map to overestimate the glacier area, potentially due to classifying 3970 Introduction

Conclusions References
Tables Figures

Back Close
Full persistent snow as glacial ice in the 1960's-1970's source aerial imagery.In contrast, another study (Kulkarni, 1992b), reported a glacierized area of 431 km 2 for Sikkim in 1987/88 based on Indian IRS-1A and Landsat data.Compared to our Corona inventory, this would suggest an area loss of 42 % since the 1960's-1970's (2.1 % yr −1 ) followed by a subsequent glacier growth in the 2000s decade (based on our Landsat analysis).
We consider the 1987/89 estimates to be highly unreliable, given that there are no glacier surges that might induce an apparent "glacier growth".Besides, the 1987/89 area estimate is smaller than the 2000 area estimated both by our study and by Mool et al. (2002).We speculate that these differences may be due to omissions of some debris-covered tongues from the glacier maps.

Summary and conclusions
In this study we combined remote sensing data from various sensors to construct a new glacier inventory for the eastern Himalaya, and to quantify glacier area and elevation changes in the last four decades.We have constructed two updated glacier inventories for this part of the Himalaya, based on 1962 Corona and 2000 Landsat/ASTER imagery, respectively.Spatial trends of glacier area and elevation changes in the last decades include: -Larger percent of debris cover on glaciers on the southern slopes of the Himalaya (14 % of the glacierized area) compared to northern slopes (2 %), with supraglacial lakes constituting about 6 % of the debris covered area; a glacier's headwater elevation, glacier area, debris cover, aspect and precipitation; -Higher rates of retreat for clean glaciers (−34.6 %, or −0.7 % yr −1 ) on a glacier-byglacier basis, compared to debris-covered glaciers (−14.3 % or −0.3 yr −1 ) in the last decades, as noted also in other studies elsewhere (Racoviteanu et al., 2008a;Basnett et al., 2013); -General trends of thinning of debris-covered tongues (−30.8 m ± 39 m) on average over the last four decades), with thickening towards the terminus for some glaciers, and rapid growth of pro-glacier lakes for others.
In this study we showed that in spite of intensive, time-consuming pre-processing steps of the Corona imagery for high altitude rugged terrain, declassified imagery from the 1960's has an important potential for glacier change detection in data-sparse areas of the Himalaya, as pointed in other studies (Narama et al., 2007;Bolch et al., 2008a).Further work would be needed on order to use Corona stereo imagery to investigate elevation changes.In this study, we did not use a DEM constructed from Corona imagery to investigate glacier elevation changes due to the time constraints and the lack of availability of a high-resolution DEM from the 2000's.By contrast, topographic maps proved to be an important data source for computing elevation changes when georeferenced and checked carefully.We also showed the potential of ASTER day surface temperature for investigating supraglacial features and as a proxy for glacier thickness.
There is a general tendency of thicker supra-glacier debris for most debris-covered tongues, as suggested by surface temperature trends, but the link between surface elevation changes and surface temperature is not yet conclusive and requires additional investigation.The geospatial datasets and the topographic/climatic links developed in this study can be further to understand the behavior of Himalayan glaciers, such as spatial patterns of glacier melt, and the contribution of glaciers to water resources.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full    Full  Full      .Glacier elevation changes for the selected debris-covered tongues from 1960s to 2000, based on the topographic map and SRTM DEM.We note higher rates of thinning in the upper part of the debris-covered area (red areas), with a general tendency of thickening at the glacier termini and for some glaciers in the low-middle part of the debris (blue areas).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Remote sensing datasets used in this study are summarized in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |root mean square error (RMSE z ), using the accuracies of the SRTM DEM (±31 m) and the error in the digitization of the topographic map (±25 m).Sources of uncertainty are discussed in detail in Sect.5.4.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

−
10.3 % (−0.24 % yr −1 ± 0.08 % yr −1 ) from 1962 to 2006 for the 50 glaciers (Table5), with double the rates from 2000 to 2006 (−0.43 % yr −1 ± 0.9 % yr −1 ) compared to the previous period 1962 to 2000 (−0.20 ± 0.16 % yr −1 ).Our rates of change are within the range of those reported in other studies, if we consider the uncertainties in the change estimates.For example, for Sikkim, our study yielded an area loss of −88.9 ± 5 km 2 , or −13.5 % of the glacierized area (−0.36 % yr −1 ±0.17 %) for the last 38 years.In a recent study, Basnett et al. (2013) reported a glacier area change of 0.16 % yr −1 from 1989/90 to 2010 based on analysis of a few selected glaciers in Sikkim, which is about 20 % Discussion Paper | Discussion Paper | Discussion Paper | lower than our estimate.The Basnett et al. rates are also lower than the ones reported elsewhere in the eastern Himalaya by Bolch et al. (2008a) (−0.25 % yr −1 in the Khumbu region of Nepal from 1962 to 2001), using similar methodology and data sources.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | length (p < 0.05) (Table Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

-
Glacier area change amounts to −0.24 % ± 0.08 % yr −1 from the 1960's to the 2006's, with a higher rate of retreat in the last decade (−0.43 % yr −1 ± 0.1 % from 2000 to 2006) compared to the previous period (−0.2 % yr −1 ± 0.06 % from 1962 to 2000); -Greater glacier area changes for small, steep glaciers with a smaller altitudinal range and less debris cover; the amount of glacier retreat is partly influenced by Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 4 .
Figure 4. Frequency distribution of glacier parameters for the 487 glaciers in spatial domain 1 based on Landsat/ASTER analysis: (a) area; (b) slope; (c) length and (d) thickness.Glaciers smaller than 10 km 2 in area, < 2 km in length and < 30 m thickness are prevalent, with an average slope of 23• .

Figure 5 .
Figure 5. Aspect frequency distribution of the 487 glaciers in spatial domain 1 based on Landsat/ASTER analysis.On average, glaciers in this area are preferrentially oriented towards two directions: NW (300 • ) and NE (60 • ) corresponding to topographic/climatic barriers.

Figure 6 .
Figure 6.Spatial patterns in glacier area change derived from 1962 Corona and 2006 QB/WV2 data, shown on a glacier-by-glacier basis.The largest area loss is observed for a few glaciers in the southern and northern parts of the image, most likely due to uncertainties in the baseline data.

Figure 10
Figure10.Glacier elevation changes for the selected debris-covered tongues from 1960s to 2000, based on the topographic map and SRTM DEM.We note higher rates of thinning in the upper part of the debris-covered area (red areas), with a general tendency of thickening at the glacier termini and for some glaciers in the low-middle part of the debris (blue areas).

Table 1 .
Summary of satellite imagery and topographic maps used in this study.

Table 3 .
Precipitation patterns averaged for the period 1998-2010 in the four climatic/topographic zones in spatial domain 1, derived from TRMM 2B31 data.

Table 4 .
Topographic parameters for glaciers in spatial domain 1 and sub-regions based on ∼ 2000 Landsat/ASTER analysis.All parameters are presented on (a) region-by-region and (b) glacier-by-glacier basis from the SRTM DEM.Debris cover fraction is calculated as % glacier area of the debris covered tongues only.

Table 7 .
Comparison of glacier parameters for clean glaciers vs. debris-covered glaciers.

Table 8 .
Glacier area change in Sikkim based on previous studies.The percent area change is given with respect to the 1962 Corona glacier inventory from this study.