Glacier Changes in the Chhombo Chhu Watershed (CCW) of Tista basin between 1975 and 2018, Sikkim Himalaya, India

. Glaciers of the Tista basin represent an important water source for mountain communities and a large population 10 downstream. The present article attempts to assesspresents the observable changes in the Chhombo Chhu Watershed glacier area (CCW) of the Tista basin, Sikkim Himalaya. The CCW contains 74 glaciers (>0.02 km 2 ) with a mean glacier size of 0.61 km 2 . We have obtained the change of such glacier outlinesdetermined changes in glaciers from the declassified hexagon KH-9 (1975), Landsat 5 TM (1989), Landsat 7 ETM+ (2000), Landsat 5 TM (2010), and Sentinel 2A (2018) images. The total glacier area in 1975 was 62.6 ± 0.7 km 2 ; and by 2018, the area had decreased to 44.8 ±1.5 km 2 , an area loss of 17.9 ± 1.7 km 2 (0.42 ± 0.04 km 2 a –1 ). Clean 15 glaciers exhibited more area loss of by 11.8 ± 1.2 km 2 (0.27 ± 0.03 km 2 a –1 ) followed bthan y partially debris-covered and fully debris-covered glaciers. The area loss is (5.0 ± 0.4 km 2 or ( 0.12 ± 0.01 km 2 a –1 ) for partially covered glaciers and maximum debris-covered (1.0 ± 0.1 km 2 (or –0.02 ± 0.002 km 2 a –1 ) fully debris covered glaciers. glaciers. The quantum of glacier area loss in the CCW of Sikkim Himalaya is 0.62 ± 0.5 km 2 a –1 0 took its pace during 2000–2010 and it is 0.77 ± 0.6 km 2 a –1 (0.62 ± 0.5 km 2 a –1 ) and during 2010–2018. (0.77 ± 0.6 km 2 a –1 ) timeframes. Field investigations of selected glaciers and climatic records also support the glacier recession in the CCW due to a significant increase in temperature (0.25 °C a − 1 ) and more or less static precipitation since 20 1995. The dataset is now available from the Zenodo web portal: http://doi.org/10.5281/zenodo.4457183 (Chowdhury et al., 2021). (1.4 km 2 ), Baspa basin (1.7 km 2 ) Saraswati/ Alaknanda basin (3.7 km 2 ) Bhambri 2012; 2015; Das and Sharma, 2018; Debnath The present study shows a significant total glacier area loss of 17.9 ± 1.7 km 2 (28.5 ± 3.6%) with a recession rate of 0.42 ± 0.04 km 2 a –1 (0.66 ± 0.1% a –1 ) from 1975 to 2018 research presents an integrated watershed-based study of glacier change across the CCW in Sikkim Himalaya from 1975 to 2018. This glacier analysis comprises 74 glaciers with a total area of 44.8 ± 1.5 km 2 , including 64 clean with an area of 28.4 ± 1.1 km 2 (63.4% of total glacier in 2018. The mean glacier area of the watershed stands at 0.61 , with the dominance of small-sized glaciers. Our mapping reveals that there has beeis n a glacier area ± 1.7 2 in the analysis period, an equivalent to 28.5 ± 3.6% shrinkage. The overall rate of glacier area loss between 1975–2018 is 0.42 2 a − ± 0.1% a − 1 consistent other present and have witnessed a higher ± km 2 − −


Introduction
25 The best manifestation of climate change, either positive or negative, is the advance and retreat of the glaciers. Both these processes are time-dependent, for they will be determined by the size of a glacier and the related influencing factors.
Response of negative mass balance is shown in the frontal area as recession supported by identifiable landforms. Glaciers play an important role as response indicators of climate change (Bahuguna et al., 2014).
Most of the Himalayan glaciers are losing mass at different rates, such as Khatling glacier in the Garhwal Himalayas 30 and Zemu glacier in the Sikkim Himalayas are retreating at the rate of 88 ± 0.3 m a −1 and 9.1 m a −1 respectively (Raj et al., 2017;Rashid and Majeed, 2020). Exceptional trends in glaciers in the central and western part of Karakoram Himalaya also prevail (Bajracharya and Shrestha, 2011;Hewitt, 2011;Mehta, 2011. However, the future stability of such a trend on glacier status is debatable under different climate scenarios (Ruosteenoja et al., 2003;Immerzeel, 2008). Since 1901, a significant rise in daily maximum temperature of 1°C in North-East India and basin-wide warming trend of 0.6°C in the 35 Brahmaputra basin are has been documented (Dash et al., 2007;Immerzeel, 2008). Moreover, the station-based meteorological Eastern Himalayas (CISMHE, 2007;Basnett et al., 2013;Debnath et al., 2018). Some of the important glaciers in this watershed are Tista Khangse (TK), Gurudongmar, Kangchengyao, Chhuma Khangse, Tasha Khangse, and Yulhe Khangse.

Figure 1
It is now known that the climate of the Sikkim Himalaya influenced by the triple system, primarily dominated by the Indian Summer Monsoon (ISM) that is the main source of precipitation (Racoviteanu et al., 2015), and additionally receive 80 occasional precipitation through North-east (winter) monsoon and hardly from mid-latitude westerlies (Murari et al. 2014;Kumar et al., 2020). At the higher elevations, most of the precipitation is in the form of snow at any time of the year. Two types of climatic conditions govern the study region (CCW), i.e., the higher reaches of the mountain to the north of Thangu valley and the northern slope of Kangchengyao Massif is dominated by the cold semi-arid deserted climate; whereas the lower part of the watershed from Thangu is characterized by Temperate climate. The higher reaches of the mountain system 85 to the north is an extension of Trans-Himalaya of the Tibet region. It is a comparatively cold semi-arid deserted type climate, similar to the Ladakh region in the west (Brazel and Marcus, 1991) and is nearly devoid of vegetation, except for the trim shrubs close to Gurudongmar and Tso Lhamo lakes. Glaciers on the northern slope of Kangchengyao massif are in a rain-shadow area, having scanty precipitation, whereas the glaciers on the southern slopes of the massif receive abundant precipitation from the ISM during June-September. Therefore, it is assumed that the glaciers of the Sikkim 90 Himalaya are generally summer accumulation irrespective of the western Himalaya and Karakoram region (Ageta and Higuchi, 1984). The climograph (Fig. 2b) shows a long-term average temperature that does not exceed 10°C, whereas mean monthly precipitation ranges from 1.3 mm to 104 mm based on the nearest meteorological station data of Pagri (54 km from the centroid of the basin) (Fig. 2a).   A variety of remotely sensed satellite images, with different temporal, multi-spectral and medium to high spatial resolution have beenis used to delineate glacier outlines, and change detection analysis in the study region (Table 1). For this study, these images have beenare selectively chosen only for the end of the ablation period, with minimal seasonal snow and cloud coverage (November to December) and a vertical solar position to avoid shadows (Chand and Sharma, 2015). These images include; panchromatic declassified hexagon (KH9-11;1975), Landsat 5 Thematic Mapper (TM;1988 and2010), Landsat 7 Enhanced Thematic Mapper Plus (ETM+;, and Sentinel 2A (2018), and further crossverified with the high-resolution images on the Google Earth (GE) platform. Sentinel 2A has a much better spatial resolution than any other freely available data sources for glacier analysis (Paul et al., 2020). In addition, the Survey of India (SoI) topographical sheets at 1:50,000 scale with 40 m contour interval have beenare used for obtaining topographic information; and Shuttle Radar Topography Mission (SRTM; 2000) Digital Elevation Model (DEM) as a reference DEM for the delineation of drainage basins and extraction of glacier topographic parameters, respectively (Table 1). In our study area, it is statistically proven that the SRTM (30 m) has an overall edge over ASTER and CARTOSAT-1 DEMs. All the satellite images and SRTM DEM are freely obtained from the U.S. Geological Survey (USGS; http://earthexplorer.usgs.gov/) and the monthly average temperature (°C) and precipitation (mm) for Pagri Meteorological Station (PMS) from 1960-2013 from Climate Research Unit (CRU), University of East Anglia (http://www.cru.uea.ac.uk/data). By the end of the 1960s, when Corona reached its technical limits, a series of photographic reconnaissance satellites (KH9) with higher spatial resolution and broader coverage were launched by the U. S. between 1971and 1984(Surazakov and Aizen, 2010. In this study, the 10 small subsets of declassified hexagon (KH9-11) image (original frame ~ about 250×125 km on the ground) acquired during KH9 mission (1211) in 1975 has been spatially georeferenced and coregistered to Sentinel 2A (S2A) (base image) using Spline Transformation Method (STM), along with adjustment rectification algorithm, to reconstruct the correct image geometry that existed at the moment of film exposure. All of the semi-automated geometric calibration procedures have been are processed in ESRI ArcGIS 10.1 software. We have acquired 286 GCPs all over the subset image that surrounds the glacier, snow-covered mountain range, consisting of prominent peaks, stable river junctions and roads, unchanged lineament structures, and lakeshores and rock outcrops for image orientation (Bhambri et al., 2011). The quality of the KH9-11 image was excellent as well as nearly scratch and noise-free. Fig. 4a shows the distortion fields of KH9-11 that are have been noticed in the corner areas near the mountain ridges, similar to previous studies (Surazakov and Aizen, 2010) and later geometrically rectified (Fig. 4b). The horizontal shift or positional accuracy between KH-9 and S2A images has been estimated at ~3.8m (∼0.94 pixels). As a fundamental step in radiometric preprocessing, atmospheric corrections, and topographically-induced illumination variations, all the different date datasets of Multi-spectral (MS) bands of Landsat series have beenare converted from DN values into a standard meaningful physical unit such as spectral radiance (Lλ) and top of atmospheric (TOA) reflectance (ρp) (Chuvieco and Huete, 2009). Finally, the dark object subtraction (DOS) method has been adopted for the radiometric correction of each image (Chavez, 1988) using ENVI 5.1 software. All the stacked Landsat images are have been co-registered with the orthorectified S2A image. The horizontal shift of Landsat TM and pan-sharpened ETM+ has beenis calculated at ∼11.7 (∼0.39 pixels) and ∼10.2 m (∼0.68 pixels), respectively, to the base image. Finally, all the visible, NIR, and SWIR bands of the S2A image have been resampled to 10 m resolution and layer stacked in ESA SNAP 5.0 sare to convert into a Geo-tiff image file for the extraction of updated glacier inventory (2018) in the CCW.
Algebraic algorithms for image enhancement (viz. NDVI, NDWI, and NDSI) are have been applied on different spectral bands of the Landsat and Sentinel 2A images to semi-automatically map of clean glacier ice (Bolch and Kamp, 2006;Racoviteanu et al., 2008Racoviteanu et al., , 2009Bhambri et al., 2011). However, these algorithms couldn't differentiate the debris-covered glaciers correctly from the surrounding moraines or rock outcrops in the region, and therefore, they have Earth (≤5 m), Sentinel 2A (10 m), PAN-sharpened multispectral images of Landsat ETM+ (15 m), and the outcome of topographic parameters such as slope, aspect, and hill shaded relief maps, calculated from SRTM DEM, are have been used for visual rectification and check to map the glaciers (Schmidt and Nüsser, 2012). The seasonal snow cover and shadow areas have also beenare eliminated. The determination of glacier termini and glacial lakes of some important large valley glaciers (e.g., Tista Khangse, Gurudongmar, and Kangchengayao-2) have beenare mapped during multiple field expeditions, and repeat photographs taken between 2017 and 2018 (Fig. 3). Finally, the glacier vector outlines (≥ 0.02 km 2 ) has been prepared in our inventory (Chand and Sharma, 2015). For the minimum size threshold, the RGI considers 0.01 km 2 as the minimum size, but recent inventories have considered 0.02 km 2 (Frey et al., 2012;Chand and Sharma, 2015).

160
The impact of topographic control on glacier fluctuations and distribution has been is statistically evaluated. We have used the morphological classification of glaciers in the watershed using the Global Land Ice Measurements from Space (GLIMS) guidelines and Glacier Atlas of India (Rau et al., 2005;Raina and Srivastava, 2008;Chand and Sharma, 2015). We have considered the glacier with more than 50% debris-cover as termed as maximum debris-covered, less than 50% is called partially debris-covered for a qualitative differentiation and debris-free glacier has been termed as clean glacier.
The classified morphological glaciers are categorized as Simple-basin valley, Compound-basin valley, Simple mountain-basin, Cirque, Niche and Glacieret (Ice aprons & Snowfields). In this study, Rock glaciers (RG) are excluded.  : 1975-1988, 1988-2000, 2000-2010, and 2010-2018. Glacier area changes have beenare directly calculated by subtracting the total area of the recent year (2018) from the initial year (1975). Absolute and relative changes have also beenare also calculated for these periods. Fig. 12 is a visual representation of different glacier's area change in the study region.

Glacier change and uncertainty estimation
Mapping uncertainty estimation is necessary to assess the significance of the results to avoid misinterpretation of mapping of the glacier area. For each glacier, the error has been is calculated based on a buffer drawn around the outlines of the glacier using ArcGIS 10.1 software, as suggested by Granshaw and Fountain (2006) and Bolch et al. (2010a,b).
For example, 5 m of buffer size (i.e. ½ of a pixel) has beenis drawn around the original glacier outlines of Sentinel 2A image. Similarly, in the case of Landsat 5 (TM), 15 m buffer size have beenis used for the glacier outlines. For pansharpened Landsat 7 ETM+ and KH9-11 images, the buffer size were 7.5 m and 2 m, respectively. It is also observed that larger glacier outlines had relatively very small errors than the small glacier patches (Bolch et al., 2010a). In this study, the mapping uncertainty of the total glacier area were are calculated as ±0.7 km 2 (~1.2 %), ±5.3 km 2 (~8.9%), ±2.6 km 2 (~4.5%), ±4.8 km 2 (~9.5%), and ±1.5 km 2 (~3.4%) mapping uncertainty for the images of KH9-11 Hexagon (1975), Landsat TM 5 (1988), pan-sharpened Landsat ETM+ (2000, Landsat TM 5 (2010), and Sentinel-2A (2018) respectively. Glacier area change uncertainty ( ac) was also estimated following Eq. (1) (Hall et al., 2003  The long-term (1960-2013) CRU monthly mean temperature and precipitation data of the Pagri Meteorological Station (PMS) (4330 m), located approximately 54 km east from the centroid of the basin, have beenare incorporated for the climate trend analysis (Fig. 2a). The Mann-Kendall (MK) statistical method has been employed to assess the climatic trend of mean monthly temperature and precipitation, and its impact on glacier fluctuation in the study region (Kendall, 1975;Mann, 1945). A positive value of Mann-Kendall Z statistic (ZMK) indicates an increasing trend and vice versa (Cengiz et al., 2020). The magnitude of a trend (i.e., rate of change per year) for time series analysis has been carried out using Sen's slope (Q) method (Sen, 1968).

Field measurement
Field surveys have beenare carried out during the pre-monsoon (May-June) and post-monsoon (October-December) seasons between 2017 and 2018 for the validation of maps of selected glaciers (e.g., Tista Khangse, Gurudongmar, Kangchengyao-2, and many unknown cirques and niche glaciers) of the study region. It is observed that the northern slope of Kangchengyao massif has a large valley glacier mostly without much debris-cover except Kangchengyao-2 glacier ( Fig. 3a-2), with enormous proglacial lakes, crevasses, and braided glacial streams. We did not measure the termini positions of the glaciers mentioned above, as these are connected with large proglacial lakes (see  (Fig. 3). The benchmark glacial lakes (i.e., Gurudongmar lakes, Khangchung, and Tso Lhamo), glaciers and associated morphology have been verified using a handheld Garmin eTrex 30x GPS with positional accuracy of ± 3 m (WAAS-enabled). The land surface temperature on the northern and southern part of the watershed for bare rocks, grasses, moraine, water etc. has beenis measured using the Fluke infrared thermometer. The Digi-Schmidt hammer-2000 has beenis used to measure the relative hardness of boulders that infers about the relative weathering and relative time of glacier recession from the places.

Glacier inventory and its distribution in 2018
In this study, we have identified and mapped 74 glaciers larger than 0.02 km 2 (minimum size threshold) using Sentinel 2A MSI (2018) and corresponding Google Earth platform images that cover an area of 44.8 ± 1.5 km 2 ( Table 2). The current study is the latest glacier inventory map (2018) of CCW in the Sikkim Himalayas (Fig. 1b). Glaciers in the CCW range from small glacieret to large valley types, with a size range from 0.02 to 6.7 km 2 . The histogram and normal Q-Q plot suggest that the glacier area in the CCW is not normal distribution ( Fig. 5a-b). The higher values of skewness (3.96) and kurtosis (19.32) and Shapiro-Wilk normality test (0.53) at 0.00 significance level also confirm the previous assumption.
The median size of the glacier area is 0.31 km 2 in the watershed (Fig. 5c).

220
Tista Khangse is the largest glacier with an area of 6.7 ± 0.1 km 2 , which contributes as a prime source and origin of Tista River in the northeastern most corner of the watershed in Sikkim. The glacier area are divided into five different size classes such as <0.5, 0.5-1, 1-1.5, 1.5-2 and >2 km 2 (Table 2; Fig. 6a). Out of these, <0.5 km 2 glacier size class has the maximum number of glacier count (51) with an area of 10.55 ± 0.6 km 2 followed by 13, 4, 1, and 5 glaciers in the size class of 0.5-1, 1-1.5, 1.5-2 , and >2 km 2 covering an area of 8.84 ± 0.3, 4.86 ± 0.2, 1.86 ± 0.04, and 18.7 ± 0.4 km 2 respectively. The glacier size class >2 km 2 is mostly dominated by the valley (compound basin and simple basin) and mountain glaciers (simple basin). All other topographically controlled parameters according to glacier size classes are mentioned in Table 2.
The morphological classification pattern of glaciers in Sikkim (Eastern Himalaya) is different from that in the central and western Himalayan regions (Raina and Srivastava, 2008). Most of the valley (simple basin) glaciers are clean, with an area of 13.5 km 2 (75%), only 16%, and 9% are partially debris-covered (2.9 km 2 ) and maximum debris-covered (1.6 km 2 ) glaciers, respectively. Mountain (simple basin) glaciers are clean with an area of 7.8 km 2 (59%), but 41% are partially debris-covered (5.5 km 2 ). Small glaciers (viz., cirque, niche, and glacieret) are completely clean in the watershed.
Distribution of glaciers according to morphological types is also shown in Fig. 9.

250 255
Topographic parameters such as elevation and slope also play crucial roles in the variations of regional characteristics of glacier distribution (Bhambri et al., 2011). The mean elevation of the glaciers ranges from 4846 to 6691 m, with an average of 5598 m (Fig. 6a). The mean elevation is higher than that of the other basins of the central and western Himalayas, such as Alaknanda (5380 m (Bhambri et al., 2011;Frey et al., 2012;Schmidt and Nüsser, 2012;Chand and Sharma, 2015;Das and Sharma, 2018), which certainly reflects the fact that this study region located at an extreme northeast part of Sikkim Himalaya has a cold semi-arid climatic condition. The studied glaciers have almost symmetrical average median elevation (5601 m) to mean elevation and ranges from 4846 to 6666 m, which is also a good proxy for the long-term equilibrium-line altitude (ELA) estimations based on topographic parameters (Braithwaite and Raper, 2009) (Fig. 7a). The glacier termini are located around an average minimum elevation of 5348 m with varying ranges from 4688 to 6529 m. Hence, large valley glaciers (simple basins) extend to the lower elevations, while smaller glaciers have higher termini. All other topographic parameters (minimum and maximum elevation, average mean and range elevation) vary according to the size class and morphological types (Table 2 & 3). Fig. 6b shows the distribution of glaciers according to the elevation size classes and morphological types.
Glaciers in the watershed are distributed with a mean slope of 26°, ranging from 12° to 45°. It is evident that the larger glaciers have gentle slopes than the smaller ones, and morphologically, the Mountain (simple basin), Cirque, Niche, and glacieret (Ice aprons & Snowfields) glaciers have steeper slope than valley (compound basin and simple basin) glaciers (  Fig. 6 & 7). Similar average slopes were observed in the other basins of Sikkim in Eastern Himalaya (Racoviteanu et al., 2015;Debnath et al., 2019) and Western Himalayas (Frey et al., 2012).
The glacier area changes according to morphological types are tabulated in Table 6. Valley (simple basin) glaciers have lost maximum percentage area of -8.4 ± 0.6 km 2 (-0.19 ± 0.01 km 2 a -1 ) in the watershed than the other glacier morphological types. It is found that the maximum glacier area loss is taking place at the lower terminus elevation and comparatively lesser slope of large valley (simple basin) glacier than the other morphological types (Table 6; Fig. 11b-d). This study records the maximum percentage of glacier area recession and its annual change rates from the northwest (-39.8 ± 6.6% or -0.93 ± 0.2% a -1 ) aspect, followed by the west (-34.2 ± 3.7% or -0.80 ± 0.1% a -1 ) and east (-33.6 ± 5.2% or -0.78 ± 0.1% a -1 ) (Fig. 11f). But in terms of absolute area change, glaciers with northern aspect lost a maximum area of -4.6 ± 0.5 km 2 , followed by glaciers with the south (-3.1 ± 0.3 km 2 ), west (-2.8 ± 0.2 km 2 ), southeast (-2.2 ± 0.2 km 2 ) and northeast (-1.8 ± 0.1 km 2 ) aspect. Since 1975, the glaciers with a mean slope between 15° and 30° exhibit a relatively higher retreat rate (Fig. 11b). All these above factors suggest that the glacier's response to climate change is largely controlled by individual glacier morphology, clean ice and debris-cover characteristics, and its associated topographical parameters (e.g., size, length, slope, minimum, mean and range elevation, and orientation) on their surface within the study region.    (Raina and Srivastava, 2008), primarily based on the SOI topographical sheets and aerial photographs. We also compared our results with the glacier outlines compiled by ICIMOD (Bajracharya and Shrestha, 2011) and RGIv6.0 (Mool and Bajracharya, 2003;Nuimura et al., 2015). Table 7 shows that the results of observed total glacier count and their area in this current study (2018) were much different from the ICIMOD (2005 ± 3) database (79 glaciers covering a total area of 45.8 km 2 ) almost after a decade later. In comparison with RGI v6.0 (2000 ± 3) glacier inventory data (90 glaciers covering a total area of 51.1 km 2 ), glacier count was overestimated (7 glaciers), while the total area was underestimated (-6 km 2 ) in the CCW. The major reasons for these different results in the RGI v6.0 database were due to misinterpretation of single glacier domain (based on slope, aspect, and glacier divide) into multiple outlines derived from automated mapping. We did not delineate some glacier boundaries, mainly in the western part of the watershed (e.g. Lachung Khangse and other unnamed glaciers) (Fig. 13a-b).
The GSI glacier inventory presented an overestimated total glacierized area (84 glaciers covering a total area of 80.7 km 2 ) compared to our estimated area in 1975 (a decade later) using the KH-9 image. The possible reasons for this overestimation were due to topographic map scale limitations (1:50,000) and uncertainties of glacier outlines due to seasonal snow and debris cover (Raina and Srivastava, 2008;Chand and Sharma, 2015). A previous study by Bhambri and Bolch (2009) also found certain inaccuracies with these topographic maps, originally derived from aerial photographs acquired during March-June (1964 ± 2), when seasonal snow can create a problem for the correct delineation of glaciers.
Thus, many researchers earlier reported similar disadvantages, such as (i) misinterpretation of clean and debris-covered glaciers derived from automated mapping; (ii) misclassification of seasonal snow as glaciers; (iii) temporal and spatial differences of acquired images and mapping period; (iv) division of single ice mass into multiple glaciers domain and viceversa without checking the local topographical distribution (i.e., slope, orientation, and glacier divide) for delineation process (Bhambri et al., 2011;Chand and Sharma, 2015;Das and Sharma, 2018).
The number of glaciers also reduced from 83 (1975)    The Spatio-temporal fluctuations of the glacier are solely dependent on the local topographic parameters such as glacier size, elevation (minimum, maximum, mean, and range), slope and orientation, and debris cover patterns under existing similar climatic conditions. The high mountain divide (i.e., Kanchengyao-Pauhunri massif) plays a significant role in imposing strong topographic influences between the northbound and southbound glacier fluctuations in the region.
Glacier size is also a determining factor for area change. The area of size class (<0.5 km 2 ) has gained mainly due to the fragmentation and area reduction of existing glaciers in the watershed ( Table 5). The overall tendency of glacier area loss <2 km 2 reveals that small-size glaciers are more sensitive and good indicators of climate change. It is because of their faster response to relatively short-term climate variations (Paul et al., 2004b). We also found a positive correlation (R 2 = 0.47) between glacier area and absolute glacier area change (Fig. 11a). About 12.1% of the total glaciers with an area of 2.1 km 2 have entirely disappeared, and 1.2% of individual glaciers have formed out of fragmentation from the main glacier mass in the study region since 1975. The debris-cover characteristics play a significant role in the glacier area change in this study as well. This study observed that the number of clean glaciers is dominant in the CCW of Sikkim Himalaya. The clean glacier area changed at a higher rate of -0.27 ± 0.03 km 2 a -1 than the partially debris-covered (-0.12 ± 0.01 km 2 a -1 ) and maximum debriscovered (-0.02 ± 0.002 km 2 a -1 ) glaciers. In this study, partially debris-covered glaciers exist on the southern (~25%) slope, whereas such debris-coverage is insignificant for the north-facing glaciers, as most of these (~26%) are clean glaciers, and such kind of observations are also found in Bhutan Himalaya by Kääb (2005). These intensive debris supply probably originated from the surrounding steep rock faces of the glaciers in the southern slopes (Frey et al., 2012), which in turn insulate the debris-covered glacier ice to form several dynamic thermokarst features such as depressions and supraglacial ponds in the lower part of the glacier tongues (Kääb, 2005). This study observed that a relatively gentle slope with lower elevations in the ablation areas could be favourable for the abundant supply of supraglacial debris.
The difference of glacier area loss on the western (2.8 ± 0.2 km 2 or 0.06 ± 0.005 km 2 a -1 ) and eastern (1.1 ± 0.1 km 2 or 0.06 ± 0.003 km 2 a -1 ) aspects during 1975-2018 can be described as more effective melting on the western slopes taking place during the afternoon due to combination of more incident solar radiation and warmer air temperature (Evans, 2006).
Besides, the north-facing glaciers (including northwest, north, and northeast) in this region are more susceptible to area loss (7.1 ± 0.7 km 2 or 0.17 ± 0.02 km 2 a -1 ) than the south-facing glaciers (including southeast, south, and southwest), which had a total area recession of 6.9 ± or 0.16 ± 0.01 km 2 a -1 . North facing glaciers are mostly clean glaciers that are highly variable with short-term temperature and precipitation changes. Moreover, direct insolation on the northern aspects glaciers lying at the fringe of the Tibetan plateau might have retreated faster under the rising temperature. Here, the "plateau flat surface heating" effect over the higher elevated semi-arid terrain reported by Brazel and Marcus (1991) can be an effective controlling parameter for the larger amount of glacier area losses from the northern aspects. It is because of the glaciers located on the higher elevations above 5300 m.a.s.l on the north face of Kanchengyao-Pauhunri massif (i.e., Gurudongmar and Tista Khangse) is an extension of Trans-Himalaya of Tibetan plateau, an undulating flat surface mostly devoid of vegetation as compared to the southern part in the Thangu valley. Semi-arid higher elevated plateau type topography receives higher insolation and as a result of warmer land surface temperatures (Brazel and Marcus, 1991).

400
Our field measurements during 1-4 December 2018 (10 am to 1 pm) confirms that mean infrared temperatures of rock (14°C), water (2.7 °C), and grass (11.5°C) on the flat surfaces in the proximity of Gurudongmar lake region (>5150 m) are higher due to more incoming solar radiation than the sloping mountainous terrain near Thangu valley (3900 m) in the south.
For example, the mean infrared temperature of rocks, water, and grasses near Thangu valley was measured as 4.2°C, 5.1°C, and 11.3°C, respectively. These recorded data reveal the concept of the "plateau flat surface heating" effect and a direct impact of the incident solar radiation and effects of shadows on the glacier area change on the northern aspects (Schmidt and Nüsser, 2012).
On the other hand, the glaciers of the Karakoram region behave distinctively different from the Hindu Kush Himalaya (HKH) (Hewitt, 1969(Hewitt, , 2005Bahuguna et al., 2014). In another study, Bhambri et al. (2017) reported that 221 glaciers in the Karakoram region have surge-type behaviour, covering an area of 7734 ± 271 km 2 (~43% of the total Karakoram glacierized area). Hewitt (2005) mentioned the probable causes of this glacier mass gain in very few areas of central Karakoram Himalaya are due to the localized impact of higher elevation and relief and a different climatic anomaly involved.
The quantum of glacier area recession rate in the CCW (0.42 ± 0.04 km 2 a -1 ) is similar to that of Changme Khangpu basin of Sikkim Himalaya (0.45 ± 0.001 km 2 a -1 ). We conclude that glaciers in Sikkim (Eastern Himalaya) are retreating at a higher magnitude as compared to the counterparts (i.e. other Himalayan regions) (Bahuguna et al., 2014) and this major shift in the glacier behaviour of Sikkim Himalaya is also noticeable since 2000 in this study (Table 4). The glaciers of the Eastern Himalaya, including Sikkim, are sensitive to climate warming due to increased precipitation from ISM and the Northeast winter monsoon (Ageta and Higuchi, 1984;Benn and Owen, 1998;Murari et al. 2014). The climatic condition of the Karakoram and western Himalaya are different from that of the counterpart (Bhambri et al., 2011). The rising temperature contributes to glacier shrinkage and thinning over the entire Tibetan Plateau (Fujita and Ageta, 2000;Yao et al., 2012) as well as in the different Himalayan regions (Ageta and Higuchi, 1984;Bhutiyani et al., 2010;Pratap et al, 2016;Das and Sharma, 2018). Table 8 shows the annual and season-wise analysis of MK, Sen's slope (Q), and linear regression tests. The mean annual temperature experienced a significant positive trend (↑) at a rate of 0.249 °C a −1 (0.05 significance level) between 1960-2013. It is important to note that the winter season experienced the maximum rising trend at the rate of temperature change (0.081 °C a −1 ), followed by pre-monsoon (0.063 °C a −1 ), monsoon (0.057 °C a −1 ) and post-monsoon (0.050 °C a −1 ) seasons. Moreover, the Himalayan glaciers are pretty sensitive to precipitation, directly or indirectly, through the albedo feedback mechanism on the short-wave radiation balance (Azam et al., 2018). Thus, an increasing trend (↑) in the annual precipitation (0.639 mm a -1 ) in the region confirms the fact of glacier shrinkage to some extent (Treichler et al., 2019). A significant rising trend (↑) of precipitation change rate was observed during pre-monsoon (0.270 mm a −1 ), followed by winter (0.228 mm a −1 ; 0.05 significance level) and post-monsoon (0.104 mm a −1 ) seasons.

475 480
But on the other hand, the precipitation change rate during the monsoon season experienced a decreasing trend of -0.197 mm a −1 . Such significant temperature and precipitation changes have immense impacts on glacier deglaciation in the investigated region since 1975. Fig. 14 confirms that there is a substantial increase in the temperature trend and more or less a static precipitation scenario since 1995, and the rate of glacier area loss also took its pace during 2000-2010 (0.62 ± 0.5 km 2 a −1 ) and 2010-2018 (0.77 ± 0.6 km 2 a −1 ) timeframes (Table 4). Based on dendrochronology, Yadava et al, (2015) have reported that 1996-2005 was the warmest period for Lachen and Lachung valley (North Sikkim) derived from mean late summer (July-September) temperature reconstruction (AD 1852(AD -2005. Table 8 Moreover, this significant trend of decreasing precipitation during summer monsoons and warming temperature during all seasons (maximum in winter season) restricts the refreezing of precipitation to form solid ice, thus reduce the degree of accumulation rates than the ablation on the existing glacier surfaces (Benn et al., 2001;Fujita and Ageta, 2000).
Similarly, Basnett and Kulkarni (2019) also reported a reduction in snowfall with a declining rate of 2.81 ± 2.02% (−0.3 ± 0.18% a −1 ) over the Sikkim Himalaya during 2002-2011 caused due to the influence of overall warming climate, especially the rise in winter minimum temperature and this correlate the finding of the present study. Thus, we can infer that glacier area shrinkage rates are also dependent on climatic factors.
Since the glacier response time is of many decades and centuries, the contemporary data may not have any perceptible reflections immediately but over the decades, it is also possible that the present trend of loss in the glaciers mass may be a response to climatic conditions that existed a few centuries ago, given the mass transfer with respect to velocity, ranging from ~5-25 m a -1 in most the Himalayan glaciers. This research presents an integrated watershed-based study of glacier change across the CCW in Sikkim Himalaya from 1975 to 2018. This glacier analysis comprises 74 glaciers with a total area of 44.8 ± 1.5 km 2 , including 64 clean glaciers with an area of 28.4 ± 1.1 km 2 (63.4% of total glacier area) in 2018. The mean glacier area of the watershed stands at 0.61 km 2 , with the dominance of small-sized glaciers. Our mapping reveals that there has beeis n a glacier area recession of 17.9 ± 1.7 km 2 in the analysis period, an equivalent to 28.5 ± 3.6% shrinkage. The overall rate of glacier area loss between 1975-2018 is 0.42 km 2 a −1 (0.66 ± 0.1% a −1 ) in the CCW, consistent with the other basins in Sikkim Himalaya. The present century, i.e., 2000The present century, i.e., -2010The present century, i.e., and 2010The present century, i.e., -2018, have witnessed a higher rate of area shrinkage of 0.62 ± 0.5 km 2 a −1 and 0.77 ± 0.6 km 2 a −1 respectively, as compared to previous timeframes between 1975-1988 (0.24 ± 0.4 km 2 a −1 ) and 1988-2000 (0.20 ± 0.5 km 2 a −1 ). Hence, the glaciers in the CCW of Sikkim (Eastern Himalaya) are retreating at a higher rate than the other parts (i.e., western and central Himalayas). Larger glaciers (>2 km 2 ) have lost greater area (7.7 ± 0.4 km 2 or -0.18 ± 0.01 km 2 a −1 ) than the small size classes since 1975. We suspect that such an anomaly might have been produced by a combined effect of higher solar incidence energy, lower terminus elevation, gentle slopes, and associated warmer air temperature. Morphologically, the valley glaciers (simple basins) have lost the maximum ice cover in the CCW. As a whole, the north-facing (including northwest, north and northeast) glaciers shrank at a higher rate of 0.17 ± 0.02 km 2 a -1 than in the other aspects. The number of clean glaciers shows a decreasing trend in loss, with a maximum area loss of 11.8 ± 1.2 km 2 (~0.27 ± 0.03 km 2 a −1 ) than the partial and maximum debris-covered glaciers in the CCW.
These factors suggest that the glacier's response to climate change is primarily controlled by glacier morphology, surface cover characteristics, and associated local topographical parameters (i.e., size, length, elevation, slope, and aspect) within the CCW. Mean annual and seasonal air temperature shows a significant positive trend since the 1960s. It appears that the rising trend of temperature during the winter season and the declining trend of precipitation in the summer season has been is onf of the the key driving factors for glacier changes in the CCW. However, the timeframe of our analysis is much smaller than the response time in the glaciers. In case of the continuation of this trend into the future, similar variability in the region will undoubtedly affect the hydroelectric power projects due to the mobilization of the vast Quaternary deposits stored all along on account of fluctuations in the discharge in this basin.

Data availability
The temporal datasets developed in this current study for glacier changes in the Chhombo Chhu Watershed of Tista basin, Sikkim Himalaya, India between 1975 and 2018 can now be accessed at Zenodo web portal: http://doi.org/10.5281/zenodo.4457183 (Chowdhury et al., 2021).

Author contributions
AC has initially designed the idea, collected all primary and secondary data through different sources, gone through the field for detailed investigation and ground truth verification, prepared the final datasets and figures and wrote of the draft manuscript. MD has accompanied during the field survey and laboratory activities. MCS, as a topical expert and joint supervisor has helped in conceptualizing, writing-review and editing. SKD, as an academic supervisor and mentor, has contributed in the entire project administration, supervision, finalizing the structure of the research article and final editing. All the authors have significantly contributed to the final form of this manuscript.
We would like to thank the Forest Department, Govt. of Sikkim for issuing the permit for field visits. The authors are also obliged to North-Eastern Hill University (Shillong) and Jawaharlal Nehru University (New Delhi) for providing different laboratory facilities and the ArcGIS 10.1 software. We thank the USGS for providing declassified Hexagon (KH-9), Landsat TM/ETM+, Sentinel 2A MSI and SRTM DEM data free of cost. We are also thankful to Mr. Tejashi Roy (Jawaharlal Nehru University, New Delhi) for helping in the field measurements during November 2018. Sincere thanks to Dr. Ian Harris, University of East Anglia (UK) for providing the CRU station data and valuable suggestions. We greatly appreciate the efforts of anonymous reviewers for constructive suggestions to improve the content and quality of our paper.    Table 5. Glacier area changes as per size class  in the CCW