Explorative Study on Mapping Surface Facies of Selected Glaciers from Chandra Basin, Himalaya Using WorldView-2 Data

: Mapping of surface glacier facies has been a part of several glaciological applications. The study of glacier facies in the Himalayas has gained momentum in the last decade owing to the implications imposed by these facies on the melt characteristics of the glaciers. Some of the most commonly reported surface facies in the Himalayas are snow, ice, ice mixed debris, and debris. The precision of the techniques used to extract glacier facies is of high importance, as the result of many cryospheric studies and economic reforms rely on it. An assessment of a customized semi-automated protocol against conventional and advanced mapping algorithms for mapping glacier surface facies is presented in this study. Customized spectral index ratios (SIRs) are developed for e ﬀ ective extraction of surface facies using thresholding in an object-based environment. This method was then tested on conventional and advanced classiﬁcation algorithms for an evaluation of the mapping accuracy for ﬁve glaciers located in the Himalayas, using very high-resolution WorldView-2 imagery. The results indicate that the object-based image analysis (OBIA) based semi-automated SIR approach achieved a higher average overall accuracy of 87.33% ( κ = 0.85) than the pixel-based image analysis (PBIA) approach. Among the conventional methods, the Maximum Likelihood performed the best, with an overall accuracy of 78.71% ( κ = 0.75). The Constrained Energy Minimization, with an overall accuracy of 68.76% ( κ = 0.63), was the best performer of the advanced algorithms. The advanced methods greatly underperformed in this study. The proposed SIRs show a promise in the mapping of minor features such as crevasses and in the discrimination between ice-mixed debris and debris. We have e ﬃ ciently mapped surface glacier facies independently of short-wave infrared bands (SWIR). There is a scope for the transferability of the proposed SIRs and their performance in varying scenarios. m) of this sensor o ﬀ ers a spatial resolution of 0. m, making it an almost indispensable component for visual interpretative applications and pansharpening of the multispectral product. The WV-2 products used in this study were acquired on 16 October 2014. The spectral resolution of this sensor permits this study to explore new spectral index ratios (SIRs) and thresholds for improved classiﬁcation techniques. In conjunction to the WV-2 data products, an ASTER Global Digital Elevation Model (GDEM) v2, with a spatial resolution of 30 m (downloaded from gdex.cr.usgs.gov / gdex / ), was utilized to generate a 3-dimensional surface by draping a pansharpened WV-2 image. This surface was used as the central reference for precise delineation of glacial boundaries.


Introduction
Snow accumulated on a glacier will metamorphose into ice on exposure to meteorological and geological alterations. This transformation leads to superficial expressions on the surface of the glacier and these expressions are as diverse in attributes as the processes that mold this transformation are intricate. These expressions are known as facies. The structural differences in these facies permit their individual detection and identification. Several facies on annexation categorically fall under two general zonations of a glacier, viz., the zone of accumulation and the zone of ablation. The equilibrium

Glacier Facies via Remote Sensing
Remote identification and extraction of surface glacier facies necessitates an emphasis on quality spectral as well as spatial resolution. Multispectral satellite imagery provides a wide range of techniques and alternatives for demonstrating the extraction of available glacier facies [18][19][20][21][22][23][24]. Radar and its associated approaches have been widely employed in this domain [25][26][27][28][29][30]. Some of the early developmental techniques for the extraction of glacier facies relied on singular datasets, such as the European remote sensing (ERS) satellites, i.e. ERS 1 and ERS 2 imagery [30,31] and Landsat Thematic Mapper (TM)/Multispectral Scanner (MSS) [22,24,32]. Singular use of the visual to near infrared (VNIR) and short wave infrared (SWIR) component of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) sensor has also been demonstrated [33]. However, a large number of studies have now begun to integrate varying sensors, data products, field spectra, and freely available digital elevation models (DEMs) into pixel-and object-based environments to extract surface glacier facies [34][35][36][37][38][39][40]. These studies provide a vast working knowledge base to design and implement an effective facies mapping technique. Although the accurate mapping of glacier facies is incorporated into a large number of glaciological studies [3], a relatively fewer number of studies have focused solely on the mapping of glacier facies. This study aims to singularly outline a robust and promising technique for mapping of glacier facies by first examining classification techniques.
Image classification is easily the most widely adopted approach for information extraction from satellite imagery [41]. The pixel, being the basic unit of digital remote sensing, is considered the primary unit of classification. Any classification scheme is designed to allocate either the complete pixel, or its deconstructed forms into groups or classes. This allocation is performed through the means of classifiers. Traditionally, two broad categories of pixel-based image analysis (PBIA) methods are utilized, i.e., supervised and unsupervised. Conventional PBIA performs supervised classification through the training of the classifier by the user, whereas unsupervised classification does not require user delegation [42]. Although PBIA has been used widely over the years, the concerns associated with these approaches range from the mixed pixel problem of coarse and moderate resolution imagery [43] to the problem of data redundancy, the salt and pepper effect, and inaccuracy when applied to high-and very high-resolution imagery [44]. The classification of a pixel does not always guarantee the accurate classification of the feature(s) within it. A pixel is classified purely based on pattern recognition of the radiative response of the object(s) represented by it [45]. Nonetheless, the spectrum of coarse to medium spatial resolution satellites is an area-weighted average of the spectra of the several possible objects contained within it [46]. With finer resolution and smaller pixel size, each pixel is more capable of providing the spectra of individual objects. However, pixels are not true representatives of geographic objects, as naturally occurring phenomena are rarely of the same geometric proportions as digital pixels. Therefore, the technique of geographic object-based image analysis (GEOBIA) has developed over the years to isolate the heterogeneity of earth objects to enable highly accurate classification [47].
Object-based image analysis (OBIA) occurs primarily through the implementation of segmentation algorithms. Segmentation categorically permits the creation of image objects (classification units) based on homogenous criteria. Unlike PBIA, the criteria incorporated into segmentation consist of spectral as well as spatial and contextual information [48]. The objects created by segmentation can replicate regional as well as local phenomena depending upon the scale factors [49]. Liu and Xia [50] state that the maximum accuracy achievable by OBIA is equivalent to the maximum accuracy of the segmentation, making it one of the most crucial factors in the classification process. Comparative assessments between OBIA and PBIA have indicated that the high accuracies achieved by OBIA hold steadfast when highand very high-resolution imageries are under consideration [51,52]. With the rise in cryosphere monitoring using very high resolution (VHR) imagery, the efficacy of OBIA against PBIA for glaciological applications has been tested with impressive results [53][54][55][56]. A large number of glacier-related information extraction from images is carried out through techniques that involve spectral band manipulations, such as band ratioing and algorithms [57]. As these techniques depend on the accuracy of the derivable spectra, the quality of the spectral signature is of high importance. Therefore, removing the atmospheric effects from satellite imagery (radiometric correction) is necessary. The primary objective for performing radiometric correction is to transform the digital brightness values per pixel into surface reflectance [58,59]. Performed as a two-part process, radiometric correction involves the conversion of image digital number values to sensor radiance and then to surface reflectance [60]. Several physics-based, radiative transfer atmospheric models have been developed to identify and retrieve surface reflectance [61,62]. One of these models, the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model [63] has been tested extensively for its capability to retrieve surface reflectance [64][65][66][67][68].

Goals of This Study
Pope et al. [69] analyzed the effect of atmospheric correction algorithms on the derivable glacial surface albedo for multi-spatial resolution datasets. The present study performs radiometric correction on WorldView-2 (WV-2) very high-resolution satellite imagery, which was acquired during early winter (16 October 2014), thereby presenting unique conditions for the development of classification methods. The imagery is subjected to PBIA using conventional as well as advanced classifiers to test the suitability of such atmospherically corrected datasets with conventional hard pixel supervised classification for extraction of glacier facies. The domain of OBIA is utilized to test the robustness of OBIA in its capacity to formulate new customized spectral index ratios with the atmospherically corrected multispectral range of WV2. These new indices are thresholded to outline distinctions between glacier facies for their semi-automated extraction. The accuracy of these classification techniques are determined by error matrices based on sample points assigned through the scrutiny of image-derived spectral plots. The study tests the efficacy of radiometrically calibrated very high-resolution imagery for comparative classification using PBIA and OBIA techniques. The end goal of this study is the accurate identification and extraction of glacier surface facies.

Study Area and Data
This study was carried out on five relatively unmonitored glaciers of the upper Chandra basin in the Himalayas, which lie within the district of Lahaul and Spiti, Himachal Pradesh, India ( Figure 1). This basin houses India's Himalayan research base, Himansh, at an altitude of 4080 m above mean sea level. This base is not visible in the present imagery as it is situated on the Sutri Dhaka glacier, which is beyond the current image extent. Unmonitored glaciers were chosen to test the applicability of the proposed methods on glaciers of several shapes and sizes. The G1 (GLIMS Ref. The central dataset utilized for classification is the WorldView-2 (WV-2) multispectral data product, whose potential for cryospheric applications has been previously investigated with promising results [70][71][72][73][74]. In addition to high spatial resolution, the success of remote detection of snow and its metamorphosed forms require sophisticated spectral modulations. The multispectral range of WV-2 consists of specific bands that are relatively unexploited to their maximum potential for glaciological applications. As stated by Aguilar et al. [75], the most relevant innovation offered by WV-2 is its potential spectral performance. With a multispectral spatial resolution of 2 m, it can enable identification of even the most obscure details. These bands are coastal (0.40-0.45 µm), blue (0.45-0.51 µm), green (0.51-0.58 µm), yellow (0.565-0.625 µm), red (0.63-0.69 µm), red edge (0.705-0.745 µm), near infrared 1 (NIR1) (0.770-0.895 µm) and near infrared 2 (NIR2) (0.86-1.04 µm). The panchromatic product (0.45-0.80 µm) of this sensor offers a spatial resolution of 0. m, making it an almost indispensable component for visual interpretative applications and pansharpening of the multispectral product. The WV-2 products used in this study were acquired on 16 October 2014. The spectral resolution of this sensor permits this study to explore new spectral index ratios (SIRs) and thresholds for improved classification techniques. In conjunction to the WV-2 data products, an ASTER Global Digital Elevation Model (GDEM) v2, with a spatial resolution of 30 m (downloaded from gdex.cr.usgs.gov/gdex/), was utilized to generate a 3-dimensional surface by draping a pansharpened WV-2 image. This surface was used as the central reference for precise delineation of glacial boundaries.

Preprocessing
The WV-2 imagery utilized for this study (courtesy of DigitalGlobe) was obtained as ortho-rectified individual tiles. The entire study region comprised of 21 such tiles. The end goal of this study is the accurate identification and extraction of glacier surface facies. Successful classification of multispectral imagery centers around the accuracy of the separable spectra [76] associated with the features under consideration. This requires thorough correction for atmospheric attenuation. The individual tiles were subjected to an extensive preprocessing protocol to attain optimal image quality and derive the final glacier extents. The protocol devised for preprocessing is illustrated in Figure 2.
Engaging the protocol outlined and tested elsewhere [53][54][55]60], the individual ortho-rectified corrected WV-2 raw tiles were subjected to a two-step radiometric corrective procedure. The steps were as follows: (1) Transformation of the raw digital number (DN) values to at sensor spectral radiance; (2) transformation of the at-sensor spectral radiance to apparent surface reflectance (atmospheric correction). The importance of atmospheric correction is significant when the methods of mapping emphasize the utilization of bands within the visible spectrum [77]. This study not only incorporates atmospherically corrected bands of the visible spectrum but also of the near-infrared region. Therefore, after investigations of relevant literature [78][79][80], this study employed the FLAASH (Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes) atmospheric correction module [40,66,69,81] to perform this task, which is based on the moderate resolution atmospheric transmission 4 (MODTRAN4) radiative transfer code [82]. This consisted of two steps; the first was the retrieval of atmospheric parameters, including an aerosol description (most importantly, the visibility or optical depth) and the column water vapor amount (Supplementary Material: Equation (S1)). The model atmosphere [63] and the aerosol/haze [83] were then used for the solution of the radiative transfer (RT) equation (Supplementary Material: Equation (S2)) to convert radiance to reflectance [84]. The FLAASH module in Environment for Visualizing Images (ENVI) v5.3 computed the parameters necessary for atmospheric correction, such as sensor altitude, view zenith angle, etc., through the selection of the sensor name and sensor type. The atmospheric and aerosol models were assigned, as prescribed by Abreu and Anderson [85], while other associated parameters were user-defined (Supplementary Table S1). This process delivered 21 atmospherically corrected multispectral and panchromatic WV-2 tiles, which were mosaicked using ENVI 5.3 to obtain the final image.
The simple seamless mosaic was suitable as radiometric calibration removed any unwanted distortions, rendering only the union of the tiles to be implemented. Efficient extraction of glacial boundaries via manual digitization would ideally require high spatial resolution to avoid ambiguity. Increasing the spatial resolution through pansharpening [86] would therefore ideally increase the efficacy of digitization. The Gram Schmidt (GS) pansharpening technique [87], touted to be one of the most widely used methods [88], was used in ENVI v5.3.
The precise identification of ice divides is crucial to separate glacial boundaries. A planar image is often incapable of providing sufficient detail for this task ( Figure 3a); in such scenarios, the use of elevation data is necessary to facilitate flow direction [89]. An ASTER GDEM v2 was elevated in ArcGIS 10.3. The pan-sharpened image was draped on this elevated GDEM with a vertical exaggeration of 0.7. This resulted in a 3-dimensional surface, essentially providing a topographic reference for manual digitization of the glacial boundaries. Figure 3 illustrates the efficacy of the topographic reference in enabling manual digitization of ice flow divides. The delineated boundaries were then used as geo-processing extents to extract the selected glaciers in ArcGIS 10.3.  Table S1). This process delivered 21 atmospherically corrected multispectral and panchromatic WV-2 tiles, which were mosaicked using ENVI 5.3 to obtain the final image. The simple seamless mosaic was suitable as radiometric calibration removed any unwanted distortions, rendering only the union of the tiles to be implemented. Efficient extraction of glacial boundaries via manual digitization would ideally require high spatial resolution to avoid ambiguity. Increasing the spatial resolution through pansharpening [86] would therefore ideally increase the efficacy of digitization. The Gram Schmidt (GS) pansharpening technique [87], touted to be one of the most widely used methods [88], was used in ENVI v5.3. elevation data is necessary to facilitate flow direction [89]. An ASTER GDEM v2 was elevated in ArcGIS 10.3. The pan-sharpened image was draped on this elevated GDEM with a vertical exaggeration of 0.7. This resulted in a 3-dimensional surface, essentially providing a topographic reference for manual digitization of the glacial boundaries. Figure 3 illustrates the efficacy of the topographic reference in enabling manual digitization of ice flow divides. The delineated boundaries were then used as geo-processing extents to extract the selected glaciers in ArcGIS 10.3.

The Identification of Glacier Surface Facies
Mapping glacier surface facies, through either pixel-or object-based methods ( Figure 2) required precise identification of these facies, both visually and spectrally, over the entire study area. This step was crucial to provide sufficient and accurate inputs for classification.

Visual Identification of Surface Facies
The WV-2 imagery displayed in Figure 4 used an RGB combination of band 5, band 3, and band 2 of the atmospherically corrected, pan-sharpened imagery.

The Identification of Glacier Surface Facies
Mapping glacier surface facies, through either pixel-or object-based methods ( Figure 2) required precise identification of these facies, both visually and spectrally, over the entire study area. This step was crucial to provide sufficient and accurate inputs for classification.

Visual Identification of Surface Facies
The WV-2 imagery displayed in Figure 4 used an RGB combination of band 5, band 3, and band 2 of the atmospherically corrected, pan-sharpened imagery.
Snow is characterized in this image by a smooth and continuous texture. It is much brighter than the other observable facies and contains little to no variation in terms of smoothness. It is observed to be brightest at the highest elevation with decreasing brightness downward. However, the visible smoothness of snow remains constant. Crevasses are easily identifiable through their characteristic troughs and disheveled appearance. Ice, on the other hand, is of a rougher texture than snow and consists of visible striations due to its flow. Its brightness is slightly lower than snow. Ice-mixed debris (IMD) is of very low brightness, is much lower in elevation, and consists of visible sporadic surface distributions of ice. Debris consists of no visible ice and therefore is of lower brightness than IMD. The shadowed area is of the lowest brightness and, while the example site shown in Figure 4 is of a constant dark tone, the shadows over the entire imagery consist of several tonal Snow is characterized in this image by a smooth and continuous texture. It is much brighter than the other observable facies and contains little to no variation in terms of smoothness. It is observed to be brightest at the highest elevation with decreasing brightness downward. However, the visible smoothness of snow remains constant. Crevasses are easily identifiable through their characteristic troughs and disheveled appearance. Ice, on the other hand, is of a rougher texture than snow and consists of visible striations due to its flow. Its brightness is slightly lower than snow. Ice-mixed debris (IMD) is of very low brightness, is much lower in elevation, and consists of visible sporadic surface distributions of ice. Debris consists of no visible ice and therefore is of lower brightness than IMD. The shadowed area is of the lowest brightness and, while the example site shown in Figure 4 is of a constant dark tone, the shadows over the entire imagery consist of several tonal variations. Supplementary Figure S1 illustrates the surface facies depicted in Figure 4 for different band combinations to highlight visual differences of the facies in each band combination.

Spectral Characteristics of Surface Facies
It is possible to perform mapping procedures over raw digital numbers (DNs), however, it is not necessary that objects having the same DNs should have the same spectral reflectance [90]. Separability of surface facies is a function of spectral characteristics as well as textural variations. Therefore, identification of the spectral differences between facies is vital to their mapping.

Spectral Characteristics of Surface Facies
It is possible to perform mapping procedures over raw digital numbers (DNs), however, it is not necessary that objects having the same DNs should have the same spectral reflectance [90]. Separability of surface facies is a function of spectral characteristics as well as textural variations. Therefore, identification of the spectral differences between facies is vital to their mapping. Figure 5 presents the average scaled spectral reflectance graph of the observed facies in this study. This graph uses the spectra collected from example sites, such as those illustrated in Figure 4, from each of the five glaciers in this study to obtain the mean spectral reflectance for the whole region. Snow is identified by its extremely high reflectance, whereas ice is characterized by a lower reflectance than snow, but significantly higher than the low reflectance of other facies. Crevasses, although much lower in reflectance than ice, maintain a similar trend as ice. The distinction between ice-mixed debris and debris is only based on the slightly higher reflectance of IMD. The shadowed area maintains low reflectance throughout the multispectral range. The identification of each of these facies and their distinctions from snow and ice are based on their variations in reduced reflectance [40]. A detailed graph with enlarged scales for each individual facies can be found in Supplementary Figure S2. The effectivity of the atmospheric correction is observable when the pre-corrected and post-corrected reflectance graphs are compared. Therefore, while Supplementary Figure Figure 5 presents the average scaled spectral reflectance graph of the observed facies in this study. This graph uses the spectra collected from example sites, such as those illustrated in Figure 4, from each of the five glaciers in this study to obtain the mean spectral reflectance for the whole region. Snow is identified by its extremely high reflectance, whereas ice is characterized by a lower reflectance than snow, but significantly higher than the low reflectance of other facies. Crevasses, although much lower in reflectance than ice, maintain a similar trend as ice. The distinction between ice-mixed debris and debris is only based on the slightly higher reflectance of IMD. The shadowed area maintains low reflectance throughout the multispectral range. The identification of each of these facies and their distinctions from snow and ice are based on their variations in reduced reflectance [40]. A detailed graph with enlarged scales for each individual facies can be found in Supplementary Figure S2. The effectivity of the atmospheric correction is observable when the pre-corrected and post-corrected reflectance graphs are compared. Therefore, while Supplementary Figure S2 displays the range of atmospherically corrected reflectance, Supplementary Figure S3 displays the range of at sensor radiance for each facies. An evaluation of the derived spectral reflectance graph is provided in Section 3.1.

Object-Based Image Analysis
Classification through the OBIA technique is performed through the succession of several steps. These steps include (a) multiresolution segmentation, (b) creation of customized spectral index ratios (SIRs), and (c) development of rule sets and threshold adjustments. Figure 2 encapsulates the progression used to perform this classification.

Object-Based Image Analysis
Classification through the OBIA technique is performed through the succession of several steps. These steps include (a) multiresolution segmentation, (b) creation of customized spectral index ratios (SIRs), and (c) development of rule sets and threshold adjustments. Figure 2 encapsulates the progression used to perform this classification.

Multiresolution Segmentation
Improving resolution results in increasing complexity of per pixel characteristics. Conventional thematic classification often associated with such complexities is insufficient for backtracking individual spectra of multiple features. This causes a pixel to be thematically classified as an average of the spectra of the features contained within it [91]. Utilization of the spatial and contextual characteristics in conjunction with the spectral signatures of pixels permits the identification of segments that are homogeneous [45]. Several segmentation algorithms have been developed over the years to segment images into objects. This study utilizes the multiresolution segmentation algorithm present in eCognition Developer 64 to segment the imagery into homogeneous objects. This algorithm provides the opportunity to place individual weightage on parameters such as layer weights, scale, shape, and compactness ( Figure 2) relative to the spectral and/or spatial characteristics based on the importance of each for creation of the homogeneous segments [92].
As illustrated in Figure 2, NIR1 and NIR2 were given the highest importance, i.e., 4 and 3, respectively, because of their absorptive capacities by water molecules [93]. The NIR bands are therefore placed in a position to enhance the discriminative prowess for isolation of homogeneous objects based on surficial moisture. The optical layers were allotted a weight of 2 for red edge, coastal, and green bands, while the yellow, blue, and red bands were weighed at 1 respectively. These weights were allotted considering that the impact of impurities present in snow prominently affects reflectance in the optical range [94]. The scale parameter, which was set to 200, was defined based on its ability to generate objects of the biggest possible scale to discriminate different image regions [95] (Figure 2). This study utilizes not only the spectral homogeneity, but also the spatial heterogeneity to isolate finer features, such as crevasses. Taking advantage of the very high spatial resolution offered by the WV-2 imagery segmentation was performed using the shape criterion of 0.4. Identification of homogeneous object (Supplementary Figure S4) characteristics necessitated a prominent discrimination through compactness, which was assigned as 0.8.

Customized Spectral Index Ratios and Thresholding
Multispectral satellite datasets provide band ratioing and index formulation capabilities. Supplementary Table S2 presents the mathematical expressions of indices conventionally used to characterize snow and ice and to separate valley rock from glacial bodies. These indices are mainly the traditional TM band ratios [24,96,97], the NDSI [98], the NDSII1 [99] and NDSII2 [33], the S3 [100] and the NDGI [33]. While these indices have been extensively used and manipulated through different sensors and band characteristics, customized indices from VHR sensors are yet to be exploited to their maximum potential. Jawak et al. [72] customized the standard Normalized Difference Snow Ice (NDSI) ( Supplementary Table S2) model into 14 different combinations using WV-2 sensor data. Robson et al. [34] combined optical, topographic, and Synthetic Aperture Radar (SAR) data to differentiate between clean ice and debris covered ice, using the customization of existing indices in an OBIA environment. The three customized spectral index ratios developed in this study ( Figure 2) are not based on manipulations of any pre-existing index but are purely developed in conjunction with the OBIA protocol to facilitate accurate classification.
The objects resulting from multiresolution segmentation were critically evaluated for understanding the mean reflectance per object. This was a time-consuming process as the total number of objects defined was 3166. Identification of the appropriate bands to formulate indices through this process would be exhaustive. Therefore, to expedite this process, the average scaled spectral reflectance plot ( Figure 5) was used to identify the mean spectral response functions for different areas of the image that visually corresponded to the identifiable surface facies types ( Figure 4). Three indices, namely SIR1, SIR2, and SIR3 were computed for performing classification. Figure 6 shows the minimum to maximum threshold range of the indices (represented as gray scale) and the mathematical expression of each index. Classification was performed using the rule set method through threshold adjustments of the devised SIRs. Each index was thresholded to suit its applicability for each facies. The established thresholds are listed in Table 1.   65 Area =300 pixel Shadowed areas in images of alpine glaciers occur due to the sensor angle and the solar azimuth. Shadowed areas in this image could not be extracted specifically, as the spectral response pattern of the various shadows across the glacier varies both within themselves as well as among each other. Therefore, a certain threshold was not definable within the customized SIRs for their extraction. Figure 7 illustrates the differential object properties of the shadowed regions, as highlighted by each of the customized SIRs. To overcome this problem, shadowed regions were extracted through manual delineation in ArcGIS 10.3. A similar problem was observed with a very minute portion of exposed valley rock within the imagery. Therefore, it was subjected to manual digitization. An aggregate of seven classes was identified and extracted by this entire procedure, namely, snow, ice, IMD, debris, crevasses, shadowed area, and valley rock. SIR1 <0.67 and SIR2 <0.8 Shadowed areas in images of alpine glaciers occur due to the sensor angle and the solar azimuth. Shadowed areas in this image could not be extracted specifically, as the spectral response pattern of the various shadows across the glacier varies both within themselves as well as among each other. Therefore, a certain threshold was not definable within the customized SIRs for their extraction. Figure 7 illustrates the differential object properties of the shadowed regions, as highlighted by each of the customized SIRs. To overcome this problem, shadowed regions were extracted through manual delineation in ArcGIS 10.3. A similar problem was observed with a very minute portion of exposed valley rock within the imagery. Therefore, it was subjected to manual digitization. An aggregate of seven classes was identified and extracted by this entire procedure, namely, snow, ice, IMD, debris, crevasses, shadowed area, and valley rock.

Pixel-Based Classification
To provide a comparative assessment against the OBIA, this study pursued pixel-based classification primarily through supervised classification techniques in ENVI 5.3. The methodology for supervised classification is depicted in Figure 2

Pixel-Based Classification
To provide a comparative assessment against the OBIA, this study pursued pixel-based classification primarily through supervised classification techniques in ENVI 5.3. The methodology for supervised classification is depicted in Figure 2.

Selection of Training Data
Supervised classification permits a user to select polygonal regions of interest (ROIs) to train the classification algorithm [45]. This study, in its effort to obtain accurate classification, used the visual differences (highlighted in Figure 4) and the spectral variations ( Figure 5 and Supplementary Figure S2) to identify training sites. Valley rock was excluded from the combined training process as it was solely exposed in the subset of the G5 glacier. The training sites outlined for each facies are illustrated in Figure 8. Multi-polygonal ROIs for each facies needed to be assigned to accommodate the spectral variations in each facies (Supplementary Figure S2) over the entire surface of the glaciers.

Application of Supervised Classifiers
The PBIA methods selected for this study consist of commonly used/conventional classifiers as well as advanced classification algorithms ( However, the MF method necessitated the application of a minimum noise fraction (MNF) transformation before classification. The MNF transform is oversimplified as a two-step principal component analysis, which enhances the quality of multispectral data by reducing computational requirements. This is performed by reducing data dimensionality and segregating the noise to yield higher order components. The higher order components consist of noise-free coherent eigen images. When used separately for classification these enhance the spectral classification results [100][101][102].
The conventional classifiers were processed using the default parameters available within the TERCAT workflow in ENVI. The stretch and rule thresholds applied to the advanced classifiers were set to default within the target detection workflow. To avoid analyst bias, the classified outputs (i.e., each facies) from a single classifier were all set to the same threshold. Thus, each classifier was assigned its own unique (default) rule threshold. The rule thresholds are listed in Table 2.

Application of Supervised Classifiers
The PBIA methods selected for this study consist of commonly used/conventional classifiers as well as advanced classification algorithms (   Post classification processing was avoided as each classifier would require different parameters and, as these parameters are mostly user-defined, it would have led to unintentional bias. The classification achieved by both the pixel-and object-based methods was subsequently tested for their thematic accuracy.

Measures of Accuracy
Assessing the accuracy of any classification via statistical analysis necessitates a distribution of reference samples. These samples first require selection of an appropriate sampling unit (pixel/block of pixels/polygon) [41]. This study assigns pixel-sampling units based on equalized random sampling. Ten sample points were distributed per facies on each glacier based on visual and spectral distinction, thus assigning each stratum an equal number of points [42,103]. Ten points were selected, as the area covered by the glaciers in the study region is small and certain facies, such as ice-mixed debris and debris, were much lower in spatial coverage than snow or glacier ice. Therefore, 300 points were used in total determine the accuracy. Foody [104] described the necessity of obtaining ground reference information as close to the date of acquisition as possible. However, this imagery was acquired at the onset of early winter. Given the inaccessibility and near-hostile condition of the Chandra-Bhaga basin during this period, field data could not be collected for validation. To overcome this setback, an equalized random sampling approach was adopted and validation was performed based on the interpretation of the spectral plots and visual analysis [33]. Avoiding any bias that may arise due to the use of spectral plots for both development of training sites and assignation of reference information, this study defined training sites as polygons and used individual pixels for ground reference information. In addition to the difference in sampling unit between reference points and training sites, both sets of data were developed mutually exclusive to each other. The overall distribution of sample points is depicted in Figure 9. The commonly used products of an error matrix, namely, error of commission (EC) (overestimation), error of omission (underestimation) (EO), overall accuracy (OA), and kappa statistic (KS) [105,106], were used to gauge the classifiers.

Results
This study aims to utilize the VHR and multispectral capabilities of the WV-2 sensor for the extraction of surface facies of selected Himalayan glaciers. As the spectra of facies are the backbone of this study, this section begins by first examining the scaled spectral reflectance graph of the surface facies. This is followed by a qualitative and quantitative evaluation of the mapped facies by the OBIA and PBIA, based on statistical results of the measures calculated from the error matrices. The measures (described in Section 3.5) were computed for every facies over each glacier individually. These were subsequently averaged to obtain the values listed in Supplementary Table S3.

Results
This study aims to utilize the VHR and multispectral capabilities of the WV-2 sensor for the extraction of surface facies of selected Himalayan glaciers. As the spectra of facies are the backbone of this study, this section begins by first examining the scaled spectral reflectance graph of the surface facies. This is followed by a qualitative and quantitative evaluation of the mapped facies by the OBIA and PBIA, based on statistical results of the measures calculated from the error matrices. The measures (described in Section 3.5) were computed for every facies over each glacier individually. These were subsequently averaged to obtain the values listed in Supplementary Table S3. Figure 5 and Supplementary Figure S2 portray the range of scaled reflectance as derived over each of the glaciers in the WV-2 imagery. The lack of short-wave infrared (SWIR) bands prevents comparison to in situ derived reflectance measurements, such as that performed by Gore et al. [40] and Negi et al. [107]. The graph of scaled reflectance ( Figure 5) obtained here holds a striking resemblance to that derived by Shukla et al. [90], in terms of the reflectance of snow. Although their study was based on Advanced Wide Field Sensor (AWiFS) bands, the trend of high reflectance of snow from the green to the NIR bands resembles the reflectance of snow derived here. Anisotropic reflectance corrections for this image in terms of the variations in bidirectional reflectance function (BDRF) were not performed for this imagery, as every glacier except G5 is of a relatively smaller size and the study focuses on five individual glaciers rather than an entire region. Such rectifications are much more important when mapping is performed over large regions [33]. The basis of differentiation between the reflectance characteristics of snow, ice, ice-mixed debris, crevasses, debris, and shadowed area (elaborated in Section 3.2.2) is their reduced reflectance properties [40]. The NIR-2 band however displays a higher average reflectance than is usually observed. This trend appears to be consistent for each facies other than snow. This could be due to the acquisition of the image at the onset of winter, as previous fresh snowfall and its subsequent freezing does affect the resultant reflectance [108]. An alternate reasoning for this could be the lack of BDRF corrections. The effect of this band is observed over each glacier. Therefore, as the compounding effect of this peculiarity will be uniform throughout the classifiers, and as the comparison of classification is with sample points derived from the same imagery, this peculiarity should not hamper the results of this study. Furthermore, the overlap of ice-mixed debris, debris, and shadowed area between bands 1 and 3 ( Figure 5) causes the dark and nearly similar appearance of these three facies (illustrated in Figure 4). Mapping of these facies would therefore depend on the other bands, which allow a degree of spectral separability between them.

Examining the Spectra
Yellow (B4) and red edge (B6) deliver higher average reflectance ( Figure 5 and Supplementary Figure S2), whereas, NIR-1 and NIR-2 show a drop and subsequent rise in average reflectance for most of the observable facies (excluding snow and shadow). This was incorporated into the development of two customized SIRs by individually dividing the B4 and B6 bands by the mean of the NIR bands ( Figure 2). A third index, generated by dividing the blue (B2) and the mean of the NIR-1/2 bands, engages the difference in range and amount of reflectance between the blue and the NIR-1/2 bands. Training sites for the supervised classifiers were allotted based on interpretations of the visual descriptors used for facies identification (Figure 4) and the spectra ( Figure 5 and Supplementary Figure S2). Mapping glacier surface facies, through either the pixelor the object-based methods (Figure 2) requires precise identification of these facies, both visually and spectrally, over the entire study area. This step is crucial to provide sufficient and accurate inputs for classification.

Evaluating the Mapped Surface Facies
It is common for studies to assess their classification using manually digitized boundaries and error-matrix-based overall producer's and user's accuracy [3,33,36,39,109]. However, such assessment does not provide a broader perspective of the errors within each class. While the kappa statistic is known to utilize the errors as well as the accuracies to provide a generally reasonable estimate of classification accuracy [110], it would be much more prudent to first evaluate the errors of omission and commission of the classification and then use the overall accuracy and kappa statistic to ascertain the final accuracy. Therefore, this study presents its results in a similar fashion.

Snow and Ice
Snow and ice in the present imagery is significantly distinct both spectrally ( Figure 5) as well as visually (Figure 4). Therefore, it should be self-evident that these two must be the most accurately mapped facies. In the object-based domain, snow was mapped solely by thresholding SIR2 (Table 1). The glaciers G1, G3, G4, and G5 required two thresholds, i.e., "SIR2 < 0.74" and "SIR2 < 0.75", with the latter accompanied by the contextual constraint "relative border to snow >0.3". A single threshold of "SIR2 < 0.70" was used to map snow for glacier G2. This parameter was assigned to prohibit close objects pertaining to the ice class from being misclassified as snow. Thresholding for the ice facies was performed using three thresholds combining SIR1 and SIR2 (Table 1). No contextual parameter was necessary for ice, as its spectrum did not overlap with any other facies.
The OBIA protocol resulted in no overestimation (EC) or underestimation (EO) for snow, i.e., 0.00% (Supplementary Table S3). The same holds true for the conventional MHD classifier. However, deliverances of the MXL and MD include overestimation of snow by 4.00% by the former, with the latter delivering an underestimation of 3.33%. The advanced classifiers greatly overestimated snow with a maximum EC of 82.00% and EO of 20.00% by the OSP. The SAM, MF, and CEM classifiers produced ECs of 22.00% and lower, whereas the resultant EO was nil by the MF, 20.00% for SAM, and 4.00% for CEM. The worst performer for snow was the OSP and the best was the OBIA protocol. The trend for the classifier performance for snow was OBIA > MHD > MD > MXL > CEM > MF > SAM > OSP. Ice was not overestimated by the OBIA protocol (EC = 0.00%), however, it was underestimated by an average EO of 14.38%. The MXL delivered no EC, but it underestimated ice by 4.62%. The MD provided no underestimations but gave an EC of 6.00%. The MHD however, gave an EO of 10.79% and an EC of 6.00%. The SAM did not underestimate ice, but rather completely overestimated it (EC = 100%). The CEM, MF, and OSP delivered ECs of 14.00%, 12.00%, and 92.00%, respectively, whereas their ECs were 27.42%, 4.04%, and 4.00%, respectively. The classifiers performed in the order of MXL > MD > OBIA > MHD > MF > CEM > OSP > SAM.

Ice-Mixed Debris and Debris
The margin of difference of the scaled spectral reflectance between the IMD facies and the debris facies was <0.1. The margin holds true for all the spectral bands, is clearly visible in Figure 5, and is highlighted in Supplementary Figure S2. The trend of reflectance for IMD and debris was similar, with the IMD reflecting slightly greater than debris ( Figure 5). The main reason for this appears to be the distribution of ice in the debris (Figure 4). Such a low margin of differentiation raises the highest probability of misclassification. Different thresholds of SIR3 made this distinction in the OBIA process (Table 1). With debris being the class providing the lowest reflectance, it was classified first using two thresholds, i.e., "SIR3 < 0" and "SIR3 < 0.1". IMD was extracted via two thresholds as well ("SIR3 < 0.35" and "SIR3 ≤ 0.455").
The OBIA method overestimated IMD by 10.00% but yielded an average underestimation of 30.44%. The MD performed better than the MHD and MXL for extracting IMD by delivering an EC of 38.00% and an EO of 16.67% (Supplementary Table S3). The MHD produced an EC of 30.00% but underestimated IMD by 38.86%, whereas the MXL delivered an EC of 38.00% and an EO of 47.40%. The MF and SAM completely overestimated IMD by 100.00%, leaving no underestimation (0.00%). The CEM and OSP provided ECs of 80.00% and 78.00%, respectively, while yielding EOs of 70.91% and 31.67%, respectively. The performance of the classifiers for IMD is thus given as OBIA > MD > MXL > MHD > OSP > CEM > MF = SAM. The OBIA protocol overestimated debris by 36.00% and underestimated it by 8.50%. The MD achieved the lowest EC of 14.00% for debris while underestimating it by 26.12%. The MHD and MXL recorded ECs of 40.00% and 65.33%, respectively, and EOs of 38.43% and 25.47%, respectively. The SAM entirely overestimated debris (100.00%) while underestimating it by 20.00%. The OSP followed by delivering an EC of 80.00% and EO of 55.00%. The CEM and MF provided ECs of 34.00% and 56.00%, respectively, with EOs of 52.06% and 30.96%, respectively. Therefore, the trend of classifier performance for the extraction of debris rests at MD > OBIA > MHD > CEM > MF > MXL > OSP > SAM.

Crevasses
In the visible spectrum, crevasses were differentiated from IMD and debris by a greater margin (>0.1) than that between IMD and debris. The NIR spectrum revealed an overlap of spectra for crevasses and IMD ( Figure 5). Visually, crevasses were easily identified by their distinct shape and disheveled appearance (refer to Section 3.2.1. and Figure 4). In such a scenario, while there may have been some possibility of misclassification, the use of the spatial characteristics of crevasses aided in their classification. An assignation of combined thresholds of SIR1, SIR2, and area-based parameters were necessary to map crevasses in the OBIA process (Table 1). Four threshold variations, namely, "SIR > 0.7"; "SIR1 > 0.65 and Area < 400 pixel"; "SIR1 > 0.65 and Area = 300 pixel"; and "SIR1 < 0.67 and SIR2 < 0.8", were established.
The OBIA overestimated crevasses by 20.00% but achieved no underestimation (0.00%). The MXL classified crevasses with an EC of 4.00% and an EO of 21.11%. The MHD and MD, on the other hand, overestimated crevasses by 20.00% and 34.00%, respectively, and underestimated it by 6.67% and 56.67% (Supplementary Table S3), respectively. Among the advanced classifiers, the SAM performed the worst by delivering an EC of 94.00% and an EO of 51.00%. The MF, OSP, and CEM yielded overestimations of 78.00%, 62.00%, and 40.00%, respectively, whereas they delivered EOs of 50.63%, 33.57%, and 34.23%, respectively. The trend of classifier performance for the extraction of crevasses is therefore given

Shadowed Area
Shadowed regions appeared visually similar to debris ( Figure 4) and they overlapped the spectral reflectance of IMD and debris at bands B1-B3 ( Figure 5). The reflectance plot for the shadowed area remained low throughout the multispectral range. The only point of differentiation between debris facies and the shadowed region was the difference in reflectance at the NIR-1/2 bands ( Figure 5). It would appear that as the developed indices incorporated a mean of the NIR bands, classification through the indices within the OBIA process should have been relatively easy. However, that was not the case, as described in Section 3.3.2 and illustrated in Figure 7. Hence, manual digitization was used to mask shadowed areas. However, to check the effectivity of the pixel based classifiers in such a scenario, training sites were assigned for shadowed areas as well (Figure 8).
The OBIA process therefore invariably achieved zero overestimation and underestimation for shadowed areas because of digitization. The MF and SAM achieved 0.00% of overestimation, but underestimated shadowed areas by 74.62% and 79.01%, respectively. The OSP yielded an average EC of 22.00% and an EO of 74.61%. The MXL and CEM registered ECs of 8.00% and 6.00%, respectively, while delivering EOs of 1.82% and 25.83%, respectively. The MHD and MD achieved a common EC of 4.0% and EOs of 16.17% and 34.02%, respectively. The trend of classifier performance for the extraction of shadowed areas is given as OBIA > MXL > MHD > CEM > MD > MF > OSP > SAM.

Overall Accuracy and Kappa Statistic
The OBIA method achieved the highest overall accuracy (Figure 10) for all glaciers, except G2. The average overall accuracy for all the five glaciers was found to be 87.33%. The average overall accuracy of the MHD yielded 78.67%. The MXL classifier outperformed the OBIA specifically for the G1 glacier (overall accuracy = 88.33%). However, when averaged over all the glaciers, it delivered an average of 78.71%. The MD classifier yielded an average overall accuracy of 72.81%. The CEM classifier delivered an average overall accuracy of 68.76%. The MF classifier provided an average overall accuracy of 58.48%. The OSP classifier yielded an overall accuracy at 32.00%. The SAM classifier was the poorest classifier among all the supervised classifiers, with an average overall accuracy of 30.48%. The average kappa statistic for each classifier for all the glaciers were also calculated along with the statistical error margin based on standard deviation to better ascertain the reliability of each classifier for such complex cryospheric applications ( Figure 11).     The overall accuracy and the kappa statistic, coupled with the interpretations of the errors of commission and omission and the producer's and user's accuracy, were used to analyze each classifier with respect to the glacier in particular and the overall study area in general. Our results indicate that the OBIA method yielded the highest average kappa value (κ = 0.86), which testifies that it is the most reliable classification method in this study. The lowest deviation of the OBIA also indicates its potential for transferability to other regions. The MHD and MXL obtained a common average kappa value of 0.75. However, the MXL encompasses a higher margin of error than the MHD. The MD (κ = 0.72), CEM (κ = 0.63), MF (κ = 0.48), OSP (κ = 0.19), and SAM (κ = 0.17) have an orderly detrimental average kappa and high instability. Therefore, the trend of classifier reliability and accuracy is found to be OBIA > MHD > MXL > MD > CEM > MF > OSP > SAM ( Figure 11). The results of each classification algorithm for every glacier are displayed in Figures 12-16.

Errors before Accuracies
It is common for studies reporting glacier classification results to present their evaluation based on either the overall accuracy, producer's accuracy or the kappa statistic, such as that given by Shukla and Ali. [109], Ali et al. [3], and Bhardwaj et al. [36]. These studies also prefer using a "globally" accepted overall accuracy value of 85% for appreciable classification results. However, Congalton and Green [106] describe that this minimum acceptable standard is not true for every classification scheme across various applications. The "minimum" acceptable standard depends on the scene and the methods employed, in conjunction with the end-application. Foody [111] elaborates on the details concerning the origin and development of this "standard" and how it does not apply in every scenario. Nevertheless, both Congalton and Green [106] and Foody [111] maintain the robustness of the error matrix and its derivatives for evaluation of accuracy. Therefore, the results presented in Section 3 are broken down facies-wise with the assessment of over/underestimation as the primary evaluator. The overall accuracy and the kappa statistic are then utilized to gain a cumulative perspective of the results. It is furthermore common for studies to state the accuracy by comparing derived boundaries of glaciers against manually digitized or globally available digital boundaries [3,39,109]. However, we have used the manually digitized boundary (refer to Section 3.1) to extract the glacier from the surrounding terrain; this negates the probability of any form of confusion between debris classes and the surrounding valley.
The resultant classifier trends are in accordance with their applicability for each facies, rather than just on the overall accuracy. It is beyond the scope of this study to determine the algorithmic cause for the poor performance of the pixel classifiers, as this would require elaborations on the mathematical workings of each classifier. The evaluation presented here is to determine which algorithm/classifier/method performs better in the given scenario for each facies.

Comparative Analysis
As per our knowledge, the facies of the five selected glaciers of this study have not been mapped before. This makes close-quarter comparison a daunting task. However, Keshri et al. [33], Shukla et al. [90], Bhardwaj et al. [36], and Gupta et al. [112] have mapped glaciers within the same area and with approximately similar classes, hence a relatable evaluation of the results of this study is possible. Studies indicate that the most commonly used classifier for glacier classification is the MXL classifier [90,109,113,114]. Although it is traditionally employed for this purpose, hosts of other classifiers are available as of today. The efficacy of such a wide variety of classifiers has not been tested for mapping glacier surface facies. Given the importance of glacier facies, it is imperative to identify an extraction method which can be widely applied across glaciered regions. While most of the available research has reported facies like snow, ice, IMD, and debris [33,90,115], none have utilized the VNIR spectrum to extract crevasses. This study is the first of its kind to use OBIA on VHR VNIR imagery to isolate crevasses. In comparison to the PBIA approach, the OBIA successfully derived crevasses with 0% underestimation (Supplementary Table S3). SIR1 and SIR2 (Figure 2), combined with the areal constraints (Table 1) proved to be efficient at this task. While the SIRs may have identified crevasses through thresholding, precise segmentation permitted the successful mapping (Supplementary Figure S4). The MXL was the only supervised classifier to achieve comparable results. However, it underestimated crevasses by 21.11%, thus reducing its suitability.
Snow and ice are two of the most commonly classified facies/zonations on glaciers. Several studies have mapped snow and ice using the NDSI, NDGI, customized ratios, sensor-specific methodologies, and in situ spectra-based differentiation [4,18,24,32,33,108,112,116]. To that extent, even the OBIA approach has been tested [39,56]. Several degrees of accuracy have been achieved using a variety of methods for the mapping of snow and ice; however, we intend to explore the ability of our customized indices against the PBIA approach. In this study, the non-overlap of snow and ice spectra should provide suitable mapping for all the classifiers. However, it appears that while the OBIA and MHD approach achieved 100% accurate classification of snow, ice was more successfully extracted by the MXL. This could be in part due to the misclassification of ice pixels into crevasses by the OBIA, though this was limited to 20%. The advanced classifiers failed to map snow and ice with sufficient accuracy irrespective of the wide spectral differentiation. The reliability of the MXL in this regard supports even some of the earlier assertions of its effectivity [32]. Most Himalayan glaciers are said to be debris covered. While the mapping of debris-covered glaciers is a problem that has been approached from various angles, thermal data and/or elevation data seem to be the most widely adopted method [35,36,117]. Elevation is found to be useful even when coupled with Synthetic Aperture Radar SAR in an object-based domain [39]. The present imagery, being of early winter, consists of widespread ice and snow. This prohibits the strict definition of supraglacial and periglacial debris. Therefore, to solve this problem and to check the capacity for debris characterization in such a scenario, this study maps ice-mixed debris (IMD) and debris. IMD is reported in several glacier classification schemes [3,33,109,118,119] and is consistent with its separability from debris due to the mixing of snow (refer to Sections 3.2.1 and 3.2.2). However, as the margin of differentiation between IMD and debris is low, their characterization is bound to consist of misclassifications. This stands true for all of the classifiers, including OBIA. The spectral overlap between IMD, debris, and shadowed areas in bands B1-B3 ( Figure 5) led the OBIA protocol to generate complex objects within the shadowed region. It was initially assumed that this overlap would hinder the characterization of IMD and debris. However, SIR3 (Figure 2), which uses B2 over the mean of B7 and B8, was used to differentiate between IMD and debris. This was interesting, as the blue band (B2) falls in the middle of the overlap between IMD, debris, and shadowed area. This distinction appears to be due to the presence of ice, causing the spectral contrast between IMD and debris. While the OBIA performed better than the PBIA for extracting IMD, the MD was the best classifier of debris. This was a surprising result as it has underperformed for every other facies. This may point toward a prospective combination of SIR3 and MD for future characterization of debris classes. Due to the complex properties of objects generated within the shadowed areas, they were not effectively mapped using the SIRs and therefore the OBIA approach required manual digitization (refer to Section 3.3.2 and Figure 7). Hence, while the MXL proved to be the most reliable identifier of snow without any manual corrections, the OBIA achieved statistically better results due to the manual intervention. None of the advanced classifiers could map shadowed areas with significant accuracy.
Although the OBIA method achieved the maximum observable accuracy, it required ≈23 minutes for processing (Supplementary Table S4). The MXL, which is the best PBIA classifier of glacier facies in this study needed ≈14 minutes to perform classification, whereas the MHD required about 13.55 minutes to complete the classification. The MF method took the longest time to process the imagery, as it required an additional step of generating a MNF transform. Therefore, as far as the temporal constraint of each method is concerned, the MHD would be a better bargain if the accuracy need not be very high. This could be for cases when regional level mapping of facies is concerned. However, when the overall accuracy of the derived facies is the most important parameter, the OBIA should invariably be the method of choice.

Significances
Precipitation over the surface of a glacier is not limited to the accumulation zone; it is spread over the entire glacier. Although this imagery is of early winter, a large distribution of snow and ice is observed (Figure 4). This is indicative of a precipitation event and subsequent refreezing of snow over the surface of the glacier. However, even in such scenarios, the sporadic distribution of snow over the glacier tongue is rarely mapped by coarser resolution and conventional methods. The OBIA approach outlined in this study has successfully managed to map this distribution on the selected glaciers (Figures 12-16). The FLAASH atmospheric correction played a key role not only in deriving spectral reflectance, but also in improving spectral separability ( Supplementary Figures S2  and S3). Comparative studies of the FLAASH model [64,66] have found it to be comparable with the result of other models such as the Atmospheric Removal Program (ATREM), the Atmospheric Correction Now (ACORN), and the Quick Atmospheric Correction (QUAC). It offers an important feature for correcting light scattered from adjacent pixels into the field-of-view [120]. It is noted that when the input parameters of atmospheric correction are reliable and accurate, the FLAASH model provides more accurate results than faster methods such as the QUAC model [64]. In this study, all relevant input parameters are available (via imagery metadata) and as the atmospheric parameters could be assigned based on prior knowledge, the FLAASH model was selected for correcting atmospheric effects. Following the similarity of the extracted scaled spectral reflectance curve to Shukla et al. [90] and the criteria for differentiation between facies based on their reduced reflectance properties [40], this study identified the spectral properties of the observable facies across all five glaciers with relative ease. The resultant classification of the OBIA method is comparable to that of similar studies such as [33,90], who have identified similar classes in the same region. While the OBIA successfully extracted snow, crevasses, IMD, and shadowed areas, the PBIA-based MXL and MD performed better for ice and debris characterization. While the OBIA did not perform well for the ice and debris facies, it cannot be dismissed due to its effectivity in the separability of snow, crevasses, and IMD. Most studies depend on the usage of SWIR bands for facies and thermal data with elevation information for mapping debris classes [3,35,36,90,109,117]. However, we have mapped facies using only the VNIR spectrum. In this context, our results match the observation of Pope and Rees [4], who emphasize the use of VNIR wavelengths for mapping facies. Mapping of facies via any OBIA method ultimately depends upon the segmentation and the thresholds established. The thresholds in Table 1 are common for all five glaciers, except for the extraction of snow in glacier G2. While the thresholds may be common in this image due to similar sensing conditions for the entire imagery, the indices could potentially be transferable to other glaciers with scene-dependent threshold adjustments, a common observation over glaciers [118,121].
If we only compare the OBIA and PBIA approaches based on the errors and accuracies, it would appear that the OBIA approach is significantly better than the PBIA (Sections 3 and 3.3). However, as it is elaborated in this study (Section 3), the complete picture is slightly different. Given the effectivity of the OBIA for snow, crevasses, and IMD, including manual digitization of boundaries and shadowed regions as well as the performance of the MXL and MD for ice and debris characterization, a synergistic combination could be an important future direction of research. Sahu and Gupta [122] propose such a merger of OBIA and PBIA methods for cryospheric applications. An effective methodology that combines the versatility of OBIA and the robustness of the PBIA could reduce the dependability of multiple kinds of ancillary inputs. Shukla and Yousuf [123] explain the necessity of ancillary inputs for the reduction of misclassification errors. However, this study has successfully mapped facies not only without multiple ancillary inputs, but also without dependence on SWIR bands. Therefore, rather than focusing on additional datasets, we propose the combination of our SIRs in an OBIA domain coupled with PBIA on VNIR spectrum datasets for the creation of a robust and transferable mechanism for mapping surface glacier facies.

Limitations
Saturation of snow in the VNIR spectrum is a recurring problem observed in similar studies. Pope and Rees [37] observed that even 16-bit Moderate Resolution Imaging Spectroradiometer (MODIS) data contained saturated snowy pixels. The WV-2 product used in this study is of 16-bit quantization; however, we have also observed saturated pixels in some parts of the snow facies. This is contributive to the characteristically high reflectance of snow observed in this study as well as by others [90]. Shadowed regions did not respond well to the indices and presented complex attributes (Figure 7) to both ice (when the shadows were lighter) and debris (when the shadows were much darker). Bhardwaj et al. [36] and Gore et al. [40] have observed similar problems in the shadowed regions. Gore et al. [40] points towards the use of higher radiometric resolution to solve this problem, however, our data is already 16-bit quantized. The non-availability of a higher radiometric product led to the use of manual digitization to identify and remove the shadowed area errors [36]. Although cumulatively the OBIA method was superior to the PBIA in terms of overall accuracy and the kappa statistic (refer to Section 3.3), it was more time-consuming and required consistent diligence. The PBIA, although less accurate, was far simpler in execution. Rastner et al. [56], who compared OBIA and PBIA for glacier mapping, have also found the limitation of time to be the most significant inhibitor in the usage of OBIA, a similar conclusion found in this study. Furthermore, the lack of any higher resolution imagery required that the same VNIR imagery be used for reference sample point generation. However, we followed the equalized random sampling approach based on Keshri et al. [33] and, as elaborated in Section 3.5, it was developed independently of the training sites used for the PBIA. The indices formed here are based on the spectral properties of WV-2. While it will be possible to test the indices on WV-2/3/4 imagery, inter-sensor transferability will not be easy. However, given the higher feature distinction capability and promising results with the WV-2 imagery, this study reiterates the importance and usefulness of sensor-specific indices for extracting maximum information from VHR multispectral data [124].

Conclusions
The present study tests the efficiency of OBIA and PBIA approaches for mapping surface glacier facies using VHR WorldView-2 imagery. Using a thorough preprocessing protocol to calibrate and correct the imagery for image analysis, subsets of five selected glaciers were extracted using manually digitized boundaries on a topographic reference. The OBIA approach permitted the coupling of segmentation and band ratioing to deliver three customized spectral index ratios that have the potential to be transferable across glaciered regions. While the OBIA approach achieved superior results for the identification of snow, crevasses, and IMD, the PBIA mapped ice and debris with the greatest accuracy. Shadowed regions were extracted manually via the OBIA, whereas the PBIA needed no such rectification. The differential capacities of the OBIA and PBIA techniques in mapping various facies have not been previously tested and no study has managed to map crevasses using the VNIR spectrum. One of the most important features of this study is its independence of SWIR bands and ancillary inputs. Among the seven classifiers tested of the PBIA approach, only MXL and MD managed to map ice and debris with high accuracies; no other classifier produced reasonable results. This study acknowledges its lack of field data for validation. However, the unbiased approach to sampling reference points highlights the strict validation performed. The focus on errors instead of accuracies scrutinizes every classifier for not only its overall performance, but also its performance for every thematic class. This heavily weighs down on the reliability of each method. The result concludes with the trend for overall classification accuracy and reliability as OBIA > MHD > MXL > MD > CEM > MF > OSP > SAM. However, as far as reliability is concerned facies-wise, the OBIA is successful for snow, crevasses, IMD, and even shadowed areas, whereas the PBIA is successful for ice and debris. While the OBIA approach is time-consuming and requires greater user intervention, the PBIA method is far less suitable for detailed feature mapping. Therefore, this study proposes a synergistic combination of OBIA and PBIA approaches to map surface glacier facies. Future research can focus on testing such a combined approach to provide detailed mapping of surface glacier facies, hence its transferability across varying scenarios and regions to fully realize their wide-ranging potential.