Glacier facies characterisation in transboundary West Sikkim Himalaya from TerraSAR-X; GLCM based classification approach

ABSTRACT The proposed facies classification approach uses three features of the GLCM with multitemporal X-band SAR data of TerraSAR-X combined with k-means clustering. This approach was tested on 15 transboundary Himalayan glaciers. They were classified into six facies namely debris, dry snow, wet snow, bare ice, percolation zone 1 and percolation zone 2 and the covered areas for these classes are 20.56%, 12.50%, 20.26%, 13.16%, 25.49% and 8.03% of total area respectively. After classification, the signature for classes are validated using the benchmark outline of East Rathong glacier, India. The accuracy achieved is 89.56% with negligible variance and kappa of 0.87.


Introduction
Glacier facies are the zones containing important information related to any glacier such as wet snow, dry snow, superimposed ice and many more.The division of glaciers into different zones was first popularised by Benson (1960).The pits and Ramsonde profiles were considered in the field for characterisation.The proposed approach is able to define different zones of the Greenland Ice Sheet.These zones are segregated into four facies as 'the ablation facies', 'the soaked facies', 'the percolation facies' and 'the dry snow facies' using measurements from stratigraphic studies.In the mentioned study ablation facies correspond to the zone of ablation in a glacier, soaked facies are snow with a moisture content or indirectly wet snow, and percolation facies are the zones where meltwater percolates into porous snow, followed by dry snow facies corresponding to the zone with dry snow all over.In a study, Müller (1962) estimated five different zones in Axel Hieiberg Island, Canadian Arctic Archipelago (Supplemental Figure S1), using different techniques.The author classified different zones through mass budget-related studies in the study area.
Facies characterisation is very important for mass balance-related studies as it not only delineates the snow-line altitude (SLA) but provides important health parameters like CONTACT Arpan Sharma arpansharma@outlook.com Supplemental data for this article can be accessed online at https://doi.org/10.1080/14498596.2022.2164085.currently available areas under different zones.With the information on SLA or equilibrium-line altitude (ELA) a proxy method of mass balance study can be applied.For example, Kulkarni (1992) used ELA to develop a method for estimation of mass balance in Himalayan glaciers.The author used the accumulation area ratio (AAR) to establish its empirical relation with the ELA.The AAR method uses the ratio of the accumulation area to the total area of a glacier.The author found that the average value of 0.44 had zero mass balance in the glaciers of the western Himalayas, though the studies mentioned require a field survey.
Following the Landsat mission by NASA, glacier facies characterisation using remote sensing and GIS is considered more efficient because in-field survey of all glaciers is not possible, especially in the Himalayan region.In this regard, Williams et al. (1991) used the Landsat Multispectral Scanner (MSS) and the Thematic Mapper (TM) with visible and infrared bands to characterise glacier facies in Bruarjokull, a surge-type outlet glacier in the Vatnajokull ice cap, Iceland.They identified morainal deposits, ice facies, transient snow line, slush zone, slush limit (diffuse) and snow facies (wet snow and frozen snow) through satellite imagery.Adam et al. (1997) used SAR imagery to map the glacier snow line in Place glacier basin, British Columbia.The data were acquired from the European Remote Sensing (ERS-1) Satellite in 1992.A supervised classification using a maximum-likelihood classifier was used to map the study area.Wet snow, glacier ice, and bed rock zones were successfully visualised using the classification approach.However, the supervised classification was successful because of the fieldwork done over different periods in the study area.Partington (1998) used ERS-1 imagery to characterise glacier facies in the Greenland ice sheet and the Wrangell St. Elias mountains, Alaska.They used three multitemporal SAR data from 1993 to create a visual channel in RGB.Facies were visible in the RGB channel.The author observed that distinct zones and facies are visible in the multi-temporal SAR data.The characterised facies are as follows; dry snow facies, ice facies, melt areas and moraine, followed by combined percolation and wet snow facies.Their study is limited to the traditional approach and they didn't develop an automated approach for direct characterisation.
Since the study by Partington (1998) the use of SAR and radar data has been more common as it has all-weather capability and can penetrate the cloud.In similar regard, Rau et al. (2000) used backscatter values of different ice parameters like stratigraphy, density, grain size, liquid water content, surface roughness, frequency, polarisation and incidence angle as suggested by Friedrich (1996) to develop a backscatter-based model for characterisation of facies.Ground-penetrating radar (GPR) based surveys are also preferred sometimes to characterise glacier facies.The use of GPR would give significant results but there is a requirement for a field-based estimation, which, as mentioned, is not possible in maximum glaciers (Langley et al. 2008).
In a more advanced manner for better accuracy, machine learning has been in use for the characterisation of facies.Machine learning-based methods are reliable as they learn the pattern of appearance of different characteristics within the data.Classification of facies through SAR backscatter, while utilising machine learning, could be useful in getting precision.Li et al. (2012) developed a model based on support vector machines (SVM) to classify wet snow, ice and others.SVM uses different hyperplanes to separate the classes (Khan et al. 2020).Further, Li et al. (2012) delineated the boundaries of wet snow and ice using the same technique.The study was carried out in the Dongkemadi glacier.
However, it has been observed that there are very few studies in the Himalayan region in regard to glacier facies, related to either remote sensing methods or field-based methods (Patel et al. 2021).This is because of the inaccessibility due to rugged terrain and uncertain weather patterns for fieldwork followed by the high computational power and specialisation required.On a general consideration that climate change is real and the role of academic research is important for policy formation, especially at the regional level, a better objective can be created that both satisfies the accuracy requirement and removes major complexities for computation.
The year 2014 brought up a new way of earth observation when the European Space Agency (ESA) launched a satellite based on C-band SAR with a temporal resolution of 6 days.Sentinel satellites have better temporal, ground and azimuth resolution as compared to previous ERS satellites.The data are provided free under the Copernicus programme to make earth observation accessible to researchers.In a similar vein, Zhou and Zheng (2017) used Sentinel 1 data to map radar glacier zones of the Antarctic peninsula.The authors used the backscattering coefficient value to visualise the zones.Threshold parameters were set for different zones to classify the facies.in a more precise manner used machine learning and polarimetric SAR (PolSAR) and classified glacier facies in the Antarctic Peninsula.The authors used a supervised SVM technique and a decisiontree classification method.PolSAR was considered because it provides enhanced backscatter values of the target.However, supervised classification becomes difficult in Himalayan glaciers due to variations in slope parameters where obtaining supervised training data is problematic; therefore, this cannot be tested on a large scale in these glaciers.To date, in supervised classification; the random forest classifier is in use for the classification of debris cover and non-debris cover only but not overall facies classification (Lu et al. 2020(Lu et al. , 2022)).In Himalayan glaciers there are fewer studies related to glacier facies characterisation.In the present context field-based facies classification in the majority of glaciers is not possible because of inaccessibility due to rugged terrain and uncertain weather patterns.This study uses the novel approach of grey-level co-occurrence matrixbased unsupervised machine learning for glacier facies characterisation in the study area to overcome field-based uncertainties.The objective of this study is to create baseline data for facies in transboundary eastern Himalaya.So for better utilisation of research for policy formation, an automated classification approach is considered while adding texture as suggested by Adam et al. (1997).

Study area
Sikkim is situated in the landlocked transboundary border state between Nepal, China and Bhutan.The western part of Sikkim has Nepal and Tibet (China) as its border.It's boundareies are 28° 06′41.25′′N and 26° 58′06.66′′N and 88° 53′41.57′′E and 88° 03′18.17′′E. The altitude ranges from 300 m near the river to 8586 m at the Kanchenjunga, which is the third-highest mountain peak in the world.Sikkim's climate ranges from subtropical in the south to tundra in the north.Landscapes above 1900-2000 m altitude experience snowfall during the glacial accumulation season, which is from October to May.As per the meteorological department of India, Sikkim receives annual rainfall ranging from 2000 to 2800 mm per year (observed rainfall variability of the Sikkim state, IMD, Pune, retrieved 18 November 2021, https://imdpune.gov.in/hydrology/rainfall%20variability%20page/sikkim_final.pdf).The state is the smallest in terms of population and the second smallest in terms of area in India.This region is very important in terms of strategic inter-boundary relations.All the major glaciers in the western part of Sikkim having transboundary significance were selected for the study (Table 1).
These glaciers are: Changshang, Rathong, South Lhonak, South Simvo, Talung, Tongshiong, Yalung, Zemu, Glacier 1, Glacier 2, Kaer,Ktr Gr 193 (Ktr 193), Middle Lhonak, North Lhonak and Lhonak Nepal (Ktr Gr 171).In the whole study area, Ktr Gr 193 is the largest glacier while Glacier 2 is the smallest.East Rathong glacier is the benchmark glacier for the field studies while Zemu glacier is the origin of the Teesta River that flows to Bangladesh via Sikkim and West Bengal.Three of the glaciers, namely Yalung, Zemu and Tongshiong, are situated on the third-highest mountain peak in the world, Kanchendzonga.The transboundary significance is given in Table 1.
The proposed work uses the GLIMS outline for the boundary delimitation.This is an online platform to view and download worldwide glacier data.However, due to the swathe limitation on the data acquired, Zemu and Talung cover only 98% of the glacier area.Hence, the terminus portion is ignored for these glaciers (Figure 1).

Datasets
TerraSAR-X Stripmap products in the range and azimuth resolution of 1.36 × 2.23 m were acquired by the European Space Agency from the German Aerospace Center (DLR) in different months of 2020-2021.TerraSAR-X satellite is an X-band SAR-based satellite operating on a 31 mm wavelength and 9.6 GHz frequency.The multitemporal data from different seasons are considered for better classification accuracy.Each dataset is acquired in the ascending pass with an incidence angle of 43.36-45.47degrees.As mentioned earlier, Partington (1998) observed that distinct zones and facies are visible in the multi-temporal SAR data.So, data from different seasons of the year 2020-2021 as mentioned in Table 2 are chosen for the study area.For validation the latest Cartosat 1 imagery provided by the Indian Space Research Organization is used.

Pre-processing
Since the data were acquired in the slant range, classification required conversion of the slant range product into a ground range product.This is to analyse and understand the textural features in the data.The first step is to convert the slant range into the ground range.
It is known that ground penetration of SAR signals induces vertical displacements which create a loss of coherence in multi-temporal slant range SAR images (Zhao et al. 2013).The acquired datasets are in a complex format having slant range information that is not important for our study.For changing a slant range to a ground range, a warp polynomial of order 3 is assigned to map the ground range pixels to slant range pixels.Each ground range pixel computes the corresponding pixel location in the slant range image using a warp polynomial.Further, a pixel value is assigned using linear interpolation in the new image.
Since the ground range products need backscatter conversion to clarify the characteristics of backscatter coefficients in each facies class, radiometric calibration is applied.Sentinels Application Platform (SNAP) uses the following equation for conversion.
Here DN i is digital number and A 2 i is the square of the amplitude value.While converting the given product; calibration constant correction, incidence angle correction and noise equivalent beta naught (NEBN) correction are also applied automatically using SNAP software.A radar backscatter has a high sensitivity to fine variations for clarifying the  characteristics of ground objects.So, its use for studying minute changes in glacier facies is ideal (Forster et al. 1996).Also, observed variations of backscatter signals in any type of snow class and crystallisation type help to classify the glacier facies more appropriately (Cuffey and Paterson 2010).
After processing the data into backscatter values, all products are saved in the linear scale.However, for visual and statistical interpretation all are converted into a decibel (dB) scale.
It is well known that SAR data have a salt-and-pepper effect that damages the interpretation of the images, making statistics unsuitable.This speckle usually occurs because of random constructive and destructive interference of the de-phased signals.It is important to despeckle the data properly before statistical analysis.
A Lee 3 × 3 speckle filter is applied for the removal of noise, although to ensure no damage to image texture a speckle filter is not applied to a copy of the SAR data to be used for GLCM estimation.A speckle filter is applied only after the GLCM estimation for this copy of the SAR data.After the filtering of multitemporal data and GLCM textures, co-registration is done for further classification.

Grey level co-occurrence matrix (GLCM)
GLCM was first used by Haralick et al. (1973), where a test was done on geological, natural and man-made features.GLCM is widely considered to estimate the secondorder texture features.Various features as suggested by Haralick et al. (1973) are as follows: angular second moment, contrast, correlation, entropy, variance, inverse difference moment, difference average, difference variance, difference entropy, sum average, sum variance and sum entropy.Every feature computes the statistical relation between the co-occurrence of different pixels in a moving window.
As per the SNAP software, the features are categorised into three groups.In the first, a Contrast Group creates three bands namely, Contrast, Dissimilarity, and Homogeneity (HOM), in terms of orderliness the bands formed are Angular Second Moment, Energy, Maximum Probability, and Entropy and finally, in the statistical group, GLCM Mean, GLCM Variance and GLCM Correlation are formed.Homogeneity, Entropy, and GLCM mean are selected for further processing to avoid redundancy.
The equations shown in Table 3 are used for the estimation of GLCM textures (Hall-Beyer 2017).

Classification and post-classification analysis
There are two classification approaches based on machine learning in remote sensing and GIS.They are supervised and unsupervised classification.Supervised classification requires a priori knowledge to derive the training classes and unsupervised classification uses clustering to differentiate the classes.However, in our case, a priori knowledge does not apply because a reconnaissance survey is almost impossible due to safety concerns over rugged terrain, uncertain weather, avalanches over high altitudes along with transboundary sensitivity.Unsupervised classification is considered widely in the above-mentioned cases.Barcaza et al. (2009), based upon normalised differential indices (NDI), used k-means unsupervised classification to group five different classes, namely non-ice, shadow, bare ice, slush and wet snow in the Northern Patagonia Icefield.There are successful attempts in classifying glacier facies using unsupervised k-means clustering in Vestfonna (Svalbard) and the Peruvian Andes and many other glaciers (Ayma et al. 2019, Barzycka et al. 2019).However, as a novel contribution high-resolution SAR data are tested on Himalayan glaciers for characterisation of glacier facies.
K-means clustering is chosen where k pixels are automatically selected to define the initial cluster centres.Further, each pixel to the nearest cluster centre defined by the Euclidean distance is assigned.Cluster centres are calculated as the means of all samples from the pixels in a cluster.When convergence is achieved after a certain iteration limit the classification is complete.The classification is done in SNAP software (Figure 2).
A post-classification signature acquisition is performed to get the possible average dB values in each band and each class.Average values with horizontal histograms are estimated using the statistics tool in SNAP for each signature.Further, majority, minority and median analyses are performed to understand the characteristics of each class in all the glaciers.The presence of the majority pixels of a single class denotes that the class is currently present in a higher number than the rest of the classes.For example, debris in the majority of pixels in glacier A is indicative of higher debris cover.The presence of median pixels for a class in glacier B is indicative of the presence of the class being neither majority nor minority.Further, the presence of a minority class is indicative of lesser occurrence in a particular glacier.
Table 3. Equations used for estimation of GLCM textures where i and j are the row and column numbers; P i;j is probability for cell i,j; n is number of rows and columns; μ i and μ j are the mean values of i and j; σ 2 i and σ 2 j are variances of i and j.

Validation
Sample facies are selected through digitised polygons from the latest available Cartosat 1 imagery acquired by the Indian Space Research Organization (ISRO) for our benchmark location, East Rathong glacier (Figure 4).Selected facies are also as per the field observation from the ablation zone of the benchmark outline as the location is accessible in the ablation zone only.Visual comparison of the acquired imagery with field observations is also made to improve the quality of validation.
A random pixel value of 0 is assigned to the manually digitised polygon while converting it to raster so that the value of 0 matches the assignment in the classified raster.The converted raster and classified results are re-sampled to the same value.r. kappa in GRASS is run to assess the classification accuracy.It estimates the error matrix of the classified results.

GLCM, classification and post classification
The GLCM textures classified are Contrast, Dissimilarity, ASM, Homogeneity, Energy, MAX, Entropy, GLCM mean, GLCM Variance and GLCM Correlation.All the textures contain lighter and darker pixels for visual characterisation.It has been observed that GLCM textures follow a similar pattern to the bands.For example, Energy, Max, ASM and Homogeneity bands show similar kinds of dark pixels.Between Contrast and Dissimilarity, the latter shows better textures than the former.The entropy band contains more light pixels while the statistical group GLCM mean has a better combination of light and dark pixels (Supplemental Figure S2).Dissimilarity, Homogeneity and GLCM mean are used for the classification, having all the necessary textures combining both the dark and light pixels.For example, the presence of debris can be visualised as a combination of lighter and darker pixels in the Dissimilarity band where the former is in a higher number, while in the Homogeneity band darker pixels are higher in number for the Debris class.In the GLCM mean band, a lighter shade shows the presence of the debris class.These three textural features along with the multi-temporal SAR data are chosen to classify all other glacier facies along with the debris.

Debris cover
Total classified pixels to the feature value are shown in Supplemental Figure S3.The mean values in dB scale for 3 October 2020, 10 January 2021, 17 March 2021, 13 June 2021and 20 September 2021 are −9.73, −6.33, −7.11, −12.91 and −10.49 with maximum error of 0.06, 0.04, 0.05, 0.06 and 0.06 respectively.The average feature value for Dissimilarity, Homogeneity and GLCM mean for debris cover are 13.86, 0.31 and 43.99 with a maximum error of 0.03, 0.001 and 0.03 respectively It is estimated that in Changshang, Glacier 1, Kaer, Lhonak (Nepal), Rathong, South Simvo, Zemu and Ktr Gr 193, debris cover is in the majority.The total debris cover in all the glaciers is 20.56%.A study by Alifu et al. (2020) classified debris cover in the northwestern Karakoram region of Gilgit-Baltistan, Pakistan, using machine learning.Though the authors use a supervised classification approach, the prediction is better in the random forest approach.Their main focus was to test only the debris cover while all the facies are considered in our case, as stated in the introduction, and unsupervised learning fits our purpose.Also, their method is a hybrid method based on both optical and microwave data.The advantage of our approach is the use of high-resolution SAR data, which is better in terms of texture.Generally, it is the major requirement for better interpretation that is fulfilled by GLCM.A study by Sharma et al. (2013) had noted that in the Indus basin the debris cover in the ablation region is 20.6%, which is higher than the 26.3% and 25.6% of the Ganga and Brahmaputra basins respectively.Their study suggested that an average of 24.16% of the ablation area is covered with debris in these three basins.This study was also supported by the study by Pandey et al. (2021), which suggested that around onefourth of the Chandra basin is debris-covered.It is observed that around one-fifth of glaciers in the transboundary Sikkim Himalaya are debris-covered.

Bare ice
Total classified pixels to the feature value are shown in supplemental Figure S4.The mean values in dB scale for 3 October 2020, 10 January 2021, 17 March 2021, 13 June 2021and 20 September 2021 are −15.73, −12.52, −11.38, −17.75 and −16.20 with maximum error of 0.05, 0.04, 0.08, 0.06 and 0.05 respectively.The average feature value for Dissimilarity, Homogeneity and GLCM mean for bare ice are 13.88, 0.27 and 18.73 with a maximum error of 0.02, 0.0007 and 0.02 respectively.
It is observed that in Glacier 2, Middle Lhonak, North Lhonak and Tongshiong bare ice is in the majority.In Yalung glacier, it is in the minority.Total bare ice is around 13.16% of the study area.A study by Kundu and Chakraborty (2015) noted that in optical data, superimposed ice is not detectable due to its having a similar spectral signature to bare ice.From this and other studies, it is noted that bare ice is mostly visible post the melting season (Shimada et al. 2016).Earlier, König et al. (2002) noted that differences in the roughness of the surface of bare and superimposed ice could be visible only in late summer in synthetic aperture radar (SAR) images.In our case, bare ice is detected without any complication in multi-temporal data.Also, the classification of superimposed ice is rare in multi-temporal data.

Percolation zone 1
Total classified pixels to the feature value are shown in Supplemental Figure S5.The mean values in dB scale for 3 October 2020, 10 January 2021, 17 March 2021, 13 June 2021and 20 September 2021 are −12.80, −9.49, −8.70, −15.27 and −13.21 with maximum error of 0.06, 0.04, 0.05, 0.06 and 0.06 respectively.The average feature value for Dissimilarity, Homogeneity and GLCM mean for the percolation zone 1 are 16.54, 0.22, and 29.31 with a maximum error of 0.02, 0.0004 and 0.03 respectively.
In Lhonak (Nepal) and Yalung, percolation zone 1 lies in the median.As observed by Pandey et al. (2020), the lower percolation zone is related to the wet snow.It is in neither the majority nor the minority in any glacier in our results, though 25.49% of the area comes under percolation zone 1.
It is observed that in Talung glacier and South Simvo, percolation zone 2 is in the minority.The upper percolation zone, which is named percolation zone 2, is the area where dry snow melts.Percolation zone 2 is found in the places where there is the presence of the class dry snow.It is 8.03% of the study area.A similar kind of result was reported by Pandey et al. (2020) in surge-type glaciers of Karakoram Himalaya.Also, Scher et al. (2021) reported that SAR data are useful in the study of melt characteristics.Here we agree with the authors that SAR data are useful in understanding the melt characteristics in Himalayan glaciers.

Dry snow
Total classified pixels to the feature value are shown in Supplemental Figure S7.The mean values in dB scale for 3 October 2020, 10 January 2021, 17 March 2021, 13 June 2021and 20 September 2021 are −17.37, −14.08, −12.94, −18.85 and −17.66 with maximum error of 0.05, 0.03, 0.04, 0.05 and 0.05 respectively.The average feature value for Dissimilarity, Homogeneity and GLCM mean for dry snow are 13.15, 0.29 and 16.75 with a maximum error of 0.02, 0.0007 and 0.03 respectively.
It is observed that in Changshang, Glacier 1, Glacier 2, Kaer, Lhonak (Nepal), Rathong, South Simvo, Tongshiong and Zemu dry snow is in the minority, and in Ktr Gr 193 it is in the median, while in Talung glacier it is in the majority.The presence of dry snow varies in different glaciers.Thakur et al. (2018) reported that between the Ganga Basin and Beas Basin there is a difference in the presence of dry snow in different glaciers.For example, Ganga Basin has dry snow of 23.61% and at Beas Basin it is around 29.38%; in our case, it is 12.50% for all the 15 glaciers in the study area.However, the study by the authors mentions the percentage of dry snow and wet snow without mentioning the percentage cover of percolation facies, which is indicative of overall snow cover classification.Our results are not based on overall snow cover characterisation.The upper percolation zone is dry in our case and the lower zone is slightly wet.Kundu and Chakraborty (2015) reported that dry snow develops in areas where the snowpack never melts; we observe it in similar areas near the peak altitude.

Wet snow
Total classified pixels to the feature value are shown in Supplemental Figure S8.The mean values in dB scale for 3 October 2020, 10 January 2021, 17 March 2021, 13 June 2021and 20 September 2021 are −13.57, −10.49, −10.22, −15.83 and −13.86 with maximum error of 0.05, 0.04, 0.05, 0.05 and 0.05 respectively.The average feature value for Dissimilarity, Homogeneity and GLCM mean for wet snow are 16.84, 0.22 and 13.78 with a maximum error of 0.02, 0.0006 and 0.04 respectively.
It is observed that in Changshang, Glacier 1, Glacier 2, MIddle Lhonak, North Lhonak, Rathong, South Lhonak, Tongshiong and Zemu wet snow is in the minority, while in Yalung it is in the majority.Total wet snow is around 20.26% in all the study areas.The signature of wet snow as per the study of Kundu and Chakraborty (2015) is similar to our study.

Validation and accuracy assessment
For accuracy assessment and validation, the Rathong glacier benchmark outline of the Science and Technology Department, Government of Sikkim (dstsikkim.gov.in), is used.The reconnaissance survey reach has been up to the wet snow zone only as upper zones are inaccessible due to security concerns.Two zones, namely debris and wet snow, are observed directly and other zones are observed remotely in the field followed by digitisation of polygons in Cartosat 1 imagery.In the Cartosat 1 imagery, different textures are identified through direct and remote observations from the field.A random pixel value of 0 is assigned to the manually digitised polygon while converting it to raster so that the value of 0 matches the assignment in the classified raster.The converted raster and classified results are re-sampled to the same value.r.kappa in GRASS is run to assess the classification accuracy.It estimates the error matrix of the classified results.An overall kappa of 0.87 is achieved with negligible variance.The overall accuracy for the classification is 89.56% (Table 4 and 5).

Limitations and further implications
The classification method used in our study is k-means clustering, which is unsupervised classification.Random forest classification in the supervised classification is tested with the digitised polygons of the benchmark outline, which were treated as training data, although the results are not significant due to the lower ratio of training and testing data (Zhu 2013).As mentioned earlier, field-based data are rare in all other glaciers in the study area except the Rathong glacier benchmark outline, which is accessible in the Ablation Zone only, so the training data are fewer for a priori knowledge-based classification.Another unsupervised classification based on clustering, i.e. expectation-maximisation (EM), is also tested.EM clustering is not efficient in our case for two reasons: firstly computational time is comparatively higher than k-means clustering.Secondly, the overall accuracy is lower, at 79%.EM clustering does not perform well in terms of computational efficiency and accuracy.A computational study by Jung et al. (2014) compared the performance of k-means clustering and EM clustering on non-raster datasets.It was found that k-means performs better than EM clustering in such datasets also.However, the computational time is higher in this case.The traditional field-based approaches used by Benson (1960), Müller (1962) and Kulkarni (1992) also have limitations in the study area as field-work is not suitable because of the reasons mentioned in the introduction section.Hence, the approach developed here could help remove field-based complexities for the creation of baseline data in Himalayan regions.
Glacier facies characterisation is useful for mass balance estimation using remote sensing and GIS or related geodetic methods.However, it is noted that the spatial heterogeneity of the glacier surface still hampers identifying the glaciers (Snehmani et al. 2015).This makes it more difficult to observe and understand glacier changes (Zhang et al. 2019).Also, as mentioned earlier, it is more challenging to find cloud-free data, especially in the larger regions of the Himalayas.A study noted that in the highmountain Asian glaciers differentiation of debris cover is challenging (Mölg et al. 2018).Our study addresses these challenges by identifying debris-covered glaciers by utilising high-resolution SAR data.The glaciers Changshang, Glacier 1, Kaer, Lhonak (Nepal), Rathong, South Simvo, Zemu and Ktr Gr 193 have high debris cover.Out of these glaciers, Zemu glacier is noted in other studies as also having high debris cover (Rashid and Majeed 2020).

Conclusion
The grey-level co-occurrence matrix is used to extract different textural features of the ice and snow in 15 glaciers of the transboundary Sikkim Himalayas.To avoid redundancy three features, namely Dissimilarity, Homogeneity and GLCM mean, were used along with the multitemporal TerraSAR-X data to classify the glacier facies (RGZs) in the study area.The six glacier facies are debris, dry snow, wet snow, bare ice, and percolation zone 1 and percolation zone 2. The 15 glaciers in the study are Changshang, Rathong, South Lhonak, South Simvo, Talung, Tongshiong, Yalung, Zemu, Glacier 1, Glacier 2, Kaer, Ktr Gr 193,Middle Lhonak,North Lhonak and Lhonak (Nepal).Areas covered for these classes in the 15 glaciers are 20.56%,12.50%, 20.26%, 13.16%, 25.49% and 8.03% of total area respectively.After the classification signature acquisition is performed.majority, minority and median analyses are also carried out to understand the characteristics of each class in different glaciers.A benchmark glacier outline of East Rathong glacier is used to validate the results.An overall accuracy of 89.56% with negligible variance and a kappa of 0.87 is achieved.

Table 1 .
Details of the glacier in the study area as per the GLIMS glacier database along with the observed transboundary significance.The largest glacier in the study area.It is in Nepal but borders both India and China.

Table 4 .
Confusion matrix for classes validated using East Rathong glacier outline.

Table 5 .
Accuracy results for classes validated using East Rathong glacier outline.