Mapping ice cliffs on debris-covered glaciers using multispectral satellite images Remote Sensing of Environment

Ice cliffs play a key role in the mass balance of debris-covered glaciers, but assessing their importance is limited by a lack of datasets on their distribution and evolution at scales larger than an individual glacier. These datasets are often derived using operator-biased and time-consuming manual delineation approaches, despite the recent emergence of semi-automatic mapping methods. These methods have used elevation or multispectral data, but the varying slope and mixed spectral signal of these dynamic features makes the transferability of these ap- proaches particularly challenging. We develop three semi-automated and objective new approaches, based on the Spectral Curvature and Linear Spectral Unmixing of multispectral images, to map these features at a glacier to regional scale. The transferability of each method is assessed by applying it to three sites in the Himalaya, where debris-covered glaciers are widespread, with varying lithologic, glaciological and climatic settings, and encompassing different periods of the melt season. We develop the new methods keeping in mind the wide range of remote sensing platforms currently in use, and focus in particular on two products: we apply the three ap- proaches at each site to near-contemporaneous atmospherically-corrected Pl ´ eiades (2 m resolution) and Sentinel-2 (10 m resolution) images and assess the effects of spatial and spectral resolution on the results. We find that the Spectral Curvature method works best for the high spatial resolution, four band Pl ´ eaides images, while a modification of the Linear Spectral Unmixing using the scaling factor of the unmixing is best for the coarser spatial resolution, but additional spectral information of Sentinel-2 products. In both cases ice cliffs are mapped with a Dice coefficient higher than 0.48. Comparison of the Pl ´ eiades results with other existing methods shows that the Spectral Curvature approach performs better and is more robust than any other existing automated or semi-automated approaches. Both methods outline a high number of small, sometimes shallow-sloping and thinly debris-covered ice patches that differ from our traditional understanding of cliffs but may have non-negligible impact on the mass balance of debris-covered glaciers. Overall these results pave the way for large scale ef- forts of ice cliff mapping that can enable inclusion of these features in debris-covered glacier melt models, as well as allow the generation of multiple datasets to study processes of cliff formation, evolution and decline.


Introduction
Debris-covered glaciers are found in all mountain ranges of the world. Despite the insulating effect of a thick debris layer (Nicholson and Benn, 2006;Ostrem, 1959) the mass losses of debris-covered glaciers have been shown to be higher than predictions from melt models (Anderson et al., 2019b;Brun et al., 2019;Gardelle et al., 2013;Kääb et al., 2012;Nuimura et al., 2012;Ragettli et al., 2016). This enhanced mass loss has been partially attributed to the supraglacial features specific to these glaciers, namely supraglacial ice cliffs and ponds Salerno et al., 2017). These features act as enhanced energy exchange zones, as they are directly exposed to incoming shortwave radiations from the sky and longwave radiations from their surroundings (Buri et al., 2016a;Miles et al., 2016), unlike the surrounding debris-covered ice. Supraglacial ice cliffs stand out at the surface of debris-covered glaciers, appearing in the rocky landscape as steep ice slopes sometimes covered by a thin layer of wet dust. Both empirical and numerical modeling studies have shown that ice cliffs melt at an enhanced rate (Anderson et al., 2019a;Brun et al., 2018;Buri et al., 2016aBuri et al., , 2016bHan et al., 2010;Juen et al., 2014;Reid and Brock, 2014;Sakai et al., 1998;Steiner et al., 2015;Watson et al., 2017b). Cliffs have been observed to melt between 2.5  and 35 (Sakai et al., 1998) times more than the surrounding debris-covered ice. As a result, they are important contributors to the mass loss of debriscovered glaciers, and can account for 20-40% of melt in the debriscovered area despite typically covering less than 10% of this area (Anderson et al., 2019b;Brun et al., 2019;Immerzeel et al., 2014;Pellicciotti et al., 2015;Sakai et al., 2002;Thompson et al., 2016). Ice cliffs are active features and can evolve, appear or disappear rapidly (Buri and Pellicciotti, 2018;Steiner et al., 2019), and past studies to assess their distribution have been challenged by the difficulty of mapping these features (Steiner et al., 2019;Watson et al., 2017a). Improved understanding of the spatial distribution and temporal evolution of cliffs is required to assess their role and melt contribution at the glacier and regional scales, but requires a mapping scheme that is consistent, automated and transferable.
The availability of a growing number of satellite images of different spatial resolution and spectral properties now allows the development of new methods based on multispectral analysis. Here we test three new methods based on multispectral images with the aim to map ice cliffs at the surface of debris-covered glaciers in an automated, transferable and objective way. Transferability between sites and images is assessed by applying all three methods to three debris-covered glaciers of the Central and Eastern Himalaya with varying lithology, altitude range and climatic settings. The images used correspond to distinct periods of the monsoon-affected melt season  to ensure transferability of the methods between different climatic conditions. We also test transferability between sensors by comparing results from nearcontemporaneous Pléiades and Sentinel-2 scenes for each site. Finally, we compare our three approaches with existing semi-automated and automated mapping methods.

Background
Despite their implications for ablation, few studies have assessed ice cliff characteristics, spatial distribution and temporal evolution at the scale of an entire glacier (Steiner et al., 2019;Thompson et al., 2016;Watson et al., 2017a). Most previous studies have identified cliffs using manual delineation, which is the most common cliff mapping approach Buri and Pellicciotti, 2018;Han et al., 2010;Juen et al., 2014;Steiner et al., 2019;Thompson et al., 2016;Watson et al., 2017a), and is also needed to train and evaluate semi-automated and automated mapping approaches (Anderson et al., 2019b;Herreid and Pellicciotti, 2018;Kraaijenbrink et al., 2016). Manual delineation of cliffs is based on the expert knowledge of the operator who can identify cliffs with varying shapes and surface characteristics reasonably accurately. However, it is a repetitive and time-consuming task and one that can be biased by an erroneous representation of a cliff leading to misinterpretation of the information contained in the remote sensing product, and can therefore be qualified as subjective. Furthermore, the quality of the manual cliff outlines will vary across a scene due to operator fatigue, surface properties and illumination. It sometimes requires several operators to validate the outlines and reduce the delineation bias (e.g. Anderson et al., 2019aAnderson et al., , 2019bSteiner et al., 2019). The manual delineation suffers from variability in cliff illumination, which partly conceals the cliffs in their own shadows. This can be dealt with by delineating only the top of the cliffs (Thompson et al., 2016;Watson et al., 2017a), but results in reduced information regarding their geometries.
More objective and less time consuming methods have been sought using classification criteria based on the slope (Herreid and Pellicciotti, 2018;Reid and Brock, 2014) or on the spectral signatures (Anderson et al., 2019b;Kraaijenbrink et al., 2016) of ice cliffs. The slope-based approaches rely on elevation data from a Digital Elevation Model (DEM) and are based on the assumption that debris cannot be sustained above a certain slope value. This threshold can be taken as a fixed value (Reid and Brock, 2014), but the stability of supraglacial debris depends on many other factors including water content, shape of the slope or status of the base of the slope (Moore, 2018). Surface slope thresholds thus depend on the glacier, climatic and geomorphological settings, and the transferability of this approach from site to site can be improved by automatically selecting threshold values (Herreid and Pellicciotti, 2018). Due to the steep, sometimes overhanging cliff faces and their narrow shapes, this approach requires fine resolution DEM products (≤ 5 m) and will be affected by DEM spatial resolution, even though efforts have been made to take into account areas at the ends of the cliffs where the slope signal is saturated (Herreid and Pellicciotti, 2018).
Alternatively, studies have tried to reproduce the classification of cliffs as made by the expert's eye by disentangling the spectral signature of cliffs from the surrounding debris and supraglacial ponds (Anderson et al., 2019b;Kraaijenbrink et al., 2016). These approaches are challenged by the varying and often heterogeneous nature of the cliff surface, from thinly debris-covered to bright bare ice and penitents (Fig. 1). An attempt was made using object-based image analysis (OBIA) applied to high spatial resolution UAV data (Kraaijenbrink et al., 2016). Although the OBIA results are promising, they were limited to a small portion of Langtang glacier with homogeneous cliff surface characteristics ( Fig. 1). Transferability of the complex training and classification scheme to other sites or coarser data is thus questionable, and probably requires new training datasets for each new scene (Kraaijenbrink et al., 2016). The OBIA approach is initiated by the segmentation of the image into near-homogeneous groups of pixels, and Mölg et al. (2020Mölg et al. ( , 2019 exploited this initial segmentation to reduce the subjectivity and workload of the manual delineation approach, even though considerable work remains to identify the cliff objects and adjust the outlines. A simpler approach successfully used the changes in brightness to pick out cliffs from the surrounding debris on one Alaskan glacier (Anderson et al., 2019b), but this requires scenes with little shading and new training data for each new scene due to the varying nature of the cliffs and debris, even at the scale of one glacier (Fig. 1). Furthermore, this method is based on the assumption that cliff brightness is lower than the brightness of the surrounding debris, which is not the case when there is bare ice at the cliff surface, or when the debris is wetted (Fig. 1, Anderson et al., 2019b). All these approaches can be qualified as 'semiautomated' since they rely on empirical thresholds that are optimized for a scene and may be more or less transferable to other scenes, as compared with the adaptive slope threshold approach which determines threshold values internally and is therefore fully automated (Herreid and Pellicciotti, 2018).
Here, we start from the hypothesis that multispectral information can be exploited in alternative ways to map cliffs in a robust manner that is transferable, objective and semi-automated. Spectral indices based on simple band ratios have a rich history in producing binary maps of target features, and have proven successful for snow or ice (e.g. Girona-Mata et al., 2019;Paul et al., 2016), water bodies (e.g. Huggel et al., 2002;Watson et al., 2018) and debris (e.g. Casey et al., 2012;Pope and Rees, 2014) in similar environments. Linear spectral unmixing (LSU), an approach that decomposes the measured spectrum of each pixel of an image as a linear combination of the spectral end-members of its constituent surfaces, can make use of the increasing amount of information contained in modern satellite products for improved mapping (e.g. Keshava and Mustard, 2002). This method has been especially successful when applied to features with a highly mixed signal (e.g. Kopačková and Hladíková, 2014;Vikhamar and Solberg, 2003). Cliff mapping is needed at different spatial and temporal scales. At the glacier scale, fine spatial resolution imagery is preferable to accurately assess cliff distribution and evolution, in order to understand the cliff dynamics in detail. At the larger scale, such as as required to inform regional melt modeling, an estimate of cliff coverage from coarser imagery may be adequate, if its bias is known. This dataset needs to rely on images that are freely available over large areas. Here we develop and assess the potential of methods that can be used for both research questions. Two distinct types of multispectral satellite imagers are suited to these distinct objectives: commercial platforms with a limited number of spectral bands but an extremely fine spatial resolution such as WorldView (Anderson et al., 2019b) or Pléiades (e.g. Berthier et al., 2014), and freely available products with coarser spatial resolution but more spectral information, like Landsat or Sentinel-2. To address these objectives in a manner that allows comparison between methods we use near-contemporaneous scenes from Pléiades and Sentinel-2 and consider the tradeoffs of spatial and spectral resolution between these data sources.

Multispectral data
We first apply the cliff delineation methods to three Pléaides acquisitions covering three different areas of the Himalaya encompassing one or more debris-covered glaciers (Fig. 2). One scene is from the premonsoon (May), one is from the beginning of the monsoon (June) and the third one is from the post-monsoon season (September, Fig. 2). We also apply the cliff delineation to three Sentinel-2 images taken within a week of the Pléiades scenes to avoid major changes in conditions and in the surface topography of the glacier, and using the same extents (Fig. 2, Table S1). Snow-and cloud-covered areas are discarded manually from the analysis (Fig. 2). The Pléiades and Sentinel-2 sensors sample the visible and near infrared spectrum, with a few observations at longer wavelengths in the case of Sentinel-2. Pléiades has a ground resolution of ~2 m for the multispectral orthoimages. For Sentinel-2 we use the four 10 m resolution bands that coincide with the four Pléiades bands in terms of central wavelength, as well as the six 20 m resolution bands (Table S2).
The Pléiades images were stereo-processed to generate 2 m resolution DEMs and orthorectified multispectral images using Rational Polynomial Coefficients (RPCs) with the AMES stereo pipeline (Beyer et al., 2018;Shean et al., 2016) for the Langtang scene, and with the Semi Global Matching (SGM) correlation algorithm (Hirschmüller, 2007) for the Satopanth scene. The Khumbu Pléiades image (16/05/ 2016) was processed to a 2 m resolution orthorectified multispectral image using the photogrammetry module of ERDAS Imagine 2015 software. The acquisition was not stereo so it was orthorectified using an earlier DEM generated from stereo Pléiades imagery (07/10/2015). The image RPCs and eight ground control points were used for georeferencing. The multispectral images were then processed to surface reflectance using the GRASS GIS (Neteler et al., 2012) i.atcorr tool, which uses the 6S radiative transfer algorithm (Kotchenova et al., 2006;Vermote et al., 1997). For the Sentinel-2 data, the images were corrected to surface reflectance using the MAJA atmospheric correction processor (Hagolle et al., 2015).

New cliff mapping approaches
We develop three new approaches (Fig. 3) to map cliffs on debriscovered glaciers. As cliffs are often associated with supraglacial ponds (Steiner et al., 2019;Watson et al., 2017a), mapping them successfully also requires the separation of partially frozen ponds from cliffs, and careful delineation of the cliff-pond boundary in joint systems. Therefore, supraglacial ponds are also delineated as part of each approach. As the shape and size of cliffs can vary widely (Fig. 1), we avoid excessive morphological filtering that could bias the results and use the same basic morphological filters for all approaches to make the results directly comparable. We discard cliffs and ponds of size less than or equal to 5 pixels (20 m 2 ) for the Pléiades and 1 pixel (100 m 2 ) for the Sentinel-2 imagery, reflecting likely spectral mapping errors (Salerno et al., 2012;Watson et al., 2016). We also morphologically fill the ponds, to account for partially frozen ponds with ice in the center.

Spectral curvature (SC)
Ice cliffs can be highly variable in their surface condition (Herreid and Pellicciotti, 2018). Although literature tends to characterize them as 'exposed' glacier ice, they are often covered by a thin veneer of debris. They also experience enhanced melt when not covered by snow, leading to a wetted, dirty ice appearance in the melt season ( Fig. 1). We therefore expect the spectra of the cliffs' pixels to be a mixture of 'pure' spectra of ice, water and debris in varying proportions. The rock debris component can also have varying reflectance values depending on factors such as its composition or water content. A study of the mineral content of the debris on the Ngozumpa and Khumbu glaciers showed varying values of reflectance depending on the lithology (Casey et al., 2012, Fig. 4), but the spectra of rock debris in the visible and near infrared are generally relatively flat with a small positive slope (Fig. 4). This contrasts strongly with the spectra of ice and water, which despite different reflectance values, both have a concave shape in the visible range (Fig. 4). The values for ice are typically higher than for water, although varying turbidity of supraglacial ponds can increase brightness (Kraaijenbrink et al., 2016;Wessels et al., 2002), confounding simple classifiers. The shape and values of the cliff spectra may vary considerably depending on its surface characteristic, water content, and mineral composition of the debris veneer ( Fig. 1). A cliff with a bare ice surface should have high reflectance values for visible bands but lower reflectance in the near infra-red wavelengths, as for glacier ice (Paul et al., 2016), while this contrast can be attenuated if the cliff is covered by a thin debris layer (Fig. 4).
We first take advantage of the 0.6 μm reflectance peak that is evident in the spectral signature of water in the visible range ( Fig. 4) to map the supraglacial ponds using the Normalized Difference Water Index (NDWI) (Fig. 5), a metric used to map water bodies (McFeeters, 1996) and that was previously applied successfully to map lakes and ponds on debris-covered glaciers (Bolch et al., 2008;Gardelle et al., 2011;Huggel et al., 2002;Liu et al., 2015;Miles et al., 2018;Watson et al., 2018;Wessels et al., 2002): where R(λ i ) is the reflectance value at the central wavelength λ i . λ 1 , λ 2 , λ 3 , λ 4 correspond to the four Pléiades bands in an ascending order and to the corresponding Sentinel-2 Bands 2, 3, 4 and 8 (Fig. 4, Table S2). Water bodies have higher NDWI values than cliffs because the reflectance in the near infrared is lower than for ice. We define the supraglacial ponds as the pixels with a NDWI higher than a certain threshold T NDWI , which is optimized for each scene (see Section 3.3): We then assess if the remaining pixels correspond to debris or cliffs by using the curvature C of the ice cliff spectra (Jia et al., 2019;Lee and Carder, 2000;Fig. 5), which is less pronounced than for pure ice or water but highlights the curved cliff spectra relative to the flat debris spectra ( Fig. 4): Since the absolute value of the curvature varies across the glacier depending on the cliff characteristics, we first filter the curvature by calculating the difference to the median over a moving window of 100 by 100 m. The size of the area is chosen so that a majority of pixels within this zone are debris-covered. We then identify the cliffs as the pixels with a filtered curvature C filt lower (more negative) than a certain threshold T curv , which is optimized for each scene (see Section 3.3):

Linear spectral unmixing
The LSU approach aims to provide the composition of each pixel based on its spectral signature, and can make use of the information contained in a large number of bands, while an approach using spectral indices is usually limited to three or four bands. In the case of Sentinel-2, additional information can be obtained from the six 20 m resolution bands, Bands 5, 6, 7, 8A, 11 and 12 (Fig. 4, Table S2). For this purpose, we re-sample these to 10 m using a nearest-neighbor approach, and apply the LSU to the ten resulting bands.
The LSU decomposes the spectrum of each pixel in the image as a linear combination of a set of pure spectral end-members. The coefficients of the linear combination can be interpreted as the abundance values of each of the spectral end-members in the image (Keshava and Mustard, 2002;Kopačková and Hladíková, 2014;Vikhamar and Solberg, 2003). These spectral end-members can be hand-picked directly from within a scene, at locations where it is assumed that the pixel is only composed of one of these pure elements. These surface reflectance values should be identical from one image to the other, and differences should reflect differences in surface composition. In all the following LSU and LSU with scale (LSU-s) delineations applied to the Pléiades  images we use the same spectral end-members extracted from the May 2016 Khumbu image. Similarly, for all Sentinel-2 scenes we use the same spectral end-members extracted in the Khumbu May 2016 Sentinel-2 image, at the same locations as the Pléiades spectral end-members. To extract these spectral end-members we manually select five pixels for each of the pure elements as our spectral end-points and take the average after confirming that the reflectance values are consistent between samples (Vikhamar and Solberg, 2003). For water, we sample non-frozen lakes outside of the glacier extents to avoid variations due to turbidity. Shallower lakes and presence of ice at the bottom can lead to a 60% increase in surface reflectance for on-glacier lakes (Fig. S1A), but the shape of the spectrum is preserved, so the impact on the LSU is assumed to be limited. For ice, we sample bare ice patches in the ablation area with no identifiable surface dust. We use two spectral endmembers for the debris to represent distinct lithologies, one for the light brown granitic debris that covers much of Khumbu Glacier and across the region (Hambrey et al., 2009) and one for the darker schistic debris which is predominant in the upper central area of Khumbu Glacier (Casey et al., 2012) and in some locations of the Langtang and Satopanth glaciers (Fig. 2, Fig. S1C). For the light debris, we sample the side moraines and dry debris deposits while for the dark debris we sample the debris in the upper central area of Khumbu.
The spectrum of each pixel in a multispectral image can be expressed as a linear combination of these spectral end-members (Keshava and Mustard, 2002;Kopačková and Hladíková, 2014;Vikhamar and Solberg, 2003). For each pixel, the reflectance of each band can be expressed as the weighted sum of the reflectance values of the spectral end-members for this wavelength (Kopačková and Hladíková, 2014): where R(λ i ) is the reflectance at central band wavelength λ i , n is the number of spectral end-members, r j (λ i ) are the reflectance values of endmembers j at band λ i and the coefficients α j are calculated using a linear least squares approach to solve the system. The choice of the best coefficients α j is constrained by imposing positive values: The result is characterized by an error equal to the difference of the initial and calculated spectra. The coefficients α j of the linear spectral unmixing are scaled so that the sum of the coefficients is equal to one: where: The final coefficients α j ′ can then be interpreted as the abundance values of each pure element composing the pixel spectrum. The scaling factor s = ∑ n j=1 α j gives information on the absolute reflectance values in the pixel relative to the reflectance values of the spectral end-members. A high scaling factor means that the reflectance values of the initial spectrum are high compared to the reflectance values of the spectral end-members that fit the shape best. This would be the case for example for a turbid pond (Fig. 4), which would be mapped as water but with a high scaling factor to account for its higher reflectance. Thus, the results of the linear spectral unmixing are maps of the proportion of each spectral end-member within each pixel (Fig. 6). We use the water map to delineate supraglacial ponds using a threshold value of water T water over which the pixels are considered to be ponds: Once the ponds are delineated, we use a threshold value of ice T ice to map ice cliffs: Both thresholds are optimized for each scene (Section 3.3). The ponds often have a non-negligible signature of ice in their spectrum so it is important to map them first to avoid them being mapped as cliffs.

Linear spectral unmixing with scale
Mapping cliffs with the LSU approach can be challenging, as the 'ice' component of the cliff area can be small due to their mixed surface composition (Fig. 5), and the 'water' and 'ice' spectral end-members are not independent (Fig. 4). Furthermore, the reflectance values obtained for bare ice can change by more than 100% from site to site depending on the amount of englacial debris, and differences in foliation, orientation or density (Fig. S1B). Thus, in this alternative implementation of LSU, we consider ice cliffs as falling into two categories: one in which cliffs consist of bare ice and snow, exhibiting high reflectance values, and one in which they are lightly debris covered and wet, and thus exhibit very low reflectance values. In both cases, the resulting cliff spectrum can be described as a combination of water and debris spectral end-members, but only by removing the constraint that subpixel components must sum to one (Eq. (7)). In the case of bright cliffs, the required combination would sum to much greater than one, corresponding to a high scale factor s while dark, thinly debris-covered cliffs scale factor s would be very low.
Therefore, in this linear spectral unmixing with scale (LSU-s) approach, we only use the pure spectra of water, light debris and dark debris for the unmixing, and use the scale map s as an additional criterion to map cliffs (Fig. 7). As s varies proportionally to one, we take the logarithm of its value and filter it similarly to the Spectral Curvature by subtracting the value of the median calculated over a 100 by 100 m moving window. The cliffs are then delineated first, using two scale thresholds, T bright for the bright cliffs and T dark for the dark cliffs: In a second step, the ponds are delineated using a NDWI threshold. A simple water threshold based on the ratio of the near infrared and the green bands was also tested for the ponds but the results were never as good as those obtained with the NDWI.

Optimization
We evaluate the results of the different mapping methods against a set of manually delineated cliffs and ponds. In each of the Pléiades scenes, we manually outline a subset of 20-40 cliffs and 25-30 ponds of different surface characteristics to serve as a validation dataset. Illumination conditions or high-density cliff zones make it difficult to pick out the individual features, so we only outline features for which our confidence in the manual delineation is high (Fig. 8B). We then define a test area with a fifty-meter buffer around all mapped features and manually map the remaining features within this buffer zone (Fig. 8C). These buffers maintain a balance between the proportion of cliff or pond and debris in the area tested. The outlines were checked independently by a second operator who modified them when deemed necessary. The test area thus defined corresponds to 5% to 13% of the total area where we applied the cliff delineation methods and is distributed across the entirety of the study glaciers to exclude bias due to particular conditions ( Fig. S2, S3, S4).
The results of each cliff mapping method are compared to the manual results for the full test area, for both cliff and pond outputs. Each pixel within the buffer zones is then categorized as true positive, true negative, false positive or false negative as part of the confusion matrix for both pond and cliff coverages, in order to quantify our results (Fig. 8D). The best threshold parameters are determined by optimizing the Dice coefficient defined as: where TP is the number of true positive pixels, FP the number of false positive pixels and FN the number of false negative pixels. As such, the Dice coefficient penalizes both missing cliff pixels and debris pixels falsely detected as cliff, and is not biased by TN, the number of true negative pixels (Buri et al., 2016a;Dice, 1945;Rittger et al., 2013). This contrasts with metrics such as the accuracy (Eq. (13)) which produces a value with direct meaning (the portion of correctly mapped pixels), but which can lose meaning for mapping low-density features. As the cliff or pond area is generally small compared to the debris area in the buffer zones (less than 10% of area), the number of true negative pixels is usually high, and results in high accuracy values regardless of the false negative and positive values: Manually mapping cliffs in the Sentinel-2 images is difficult and highly uncertain due to the coarse spatial resolution of the sensor, so we instead degrade the validation datasets derived from the corresponding Pléiades scenes. The Pléiades and Sentinel-2 images are first coregistered using a normalized cross-correlation approach with ImGRAFT (Messerli and Grinsted, 2015) by taking the average displacement value to translate the Sentinel-2 image. We then degrade the Pléiades validation dataset to the Sentinel-2 resolution of 10 m, and reclassify as cliffs or ponds the pixels that contain more than 50% of the original cliff or pond coverage.

Application of other existing methods
We compare our SC, LSU and LSU-s approaches with other existing methods by applying these to the Langtang and Satopanth Pléiades 2 m multispectral images or DEMs when required, using the same parameters as in the original studies as much as possible. No DEM is available for the Khumbu scene so we do not test the slope-based methods for this scene. We test: 1) a fully automated adaptive slope threshold (AST) approach that also extends the narrowing ends of cliffs (Herreid and Pellicciotti, 2018), 2) a semi-automated simple slope threshold (SST) approach (Reid and Brock, 2014), 3) a semi-automated adaptive binary threshold (ABT) approach (Anderson et al., 2019b), and 4) an OBIA segmentation combined with a manual delineation (OBIA-m) (Kraaijenbrink et al., 2016;Mölg et al., 2020Mölg et al., , 2019. For the automated AST, we use the same set of initial parameters as applied on Canwell and Ngozumpa glaciers (Herreid and Pellicciotti, 2018). Due to computational limitations we run this approach using DEMs bilinearly resampled to 5 m resolution, as in the original study. For the semiautomated SST and ABT approaches (Anderson et al., 2019b;Reid and Brock, 2014), we recalibrate the parameters for each image as for our own methods. For the ABT, we emulate the original implementation by defining the brightness as the mean value of the four bands, and we consider as cliff the pixels for which the difference between brightness and mean brightness across a certain window is lower than a certain threshold (Anderson et al., 2019b): We optimize the window size along with the brightness threshold. As a result, the ABT approach we test here also has less parameters than in the original study (Anderson et al., 2019b). For the OBIA-m approach (Kraaijenbrink et al., 2016) we use the open-source software SAGA GIS to segment the multispectral image (Conrad et al., 2015). We do not reproduce the classification step detailed in the original study (Kraaijenbrink et al., 2016), but rather use a manual delineation applied to the results of the segmentation step (Mölg et al., 2020(Mölg et al., , 2019. Results of each method are assessed against the same validation datasets that we use for all three scenes. We also take into account the results provided in the different studies and express them in terms of Dice coefficient to make them comparable.

Spectral curvature
Using the NDWI and Curvature, we obtain maximum values of the Dice coefficient between 0.85 and 0.9 for ponds for the three scenes (Fig. 9A) and between 0.50 and 0.69 for cliffs (Fig. 9B) for the Pléiades image. For Sentinel-2, the maximum Dice coefficient is also high for ponds, between 0.78 and 0.89, but lower for cliffs, between 0.44 and 0.45 (Fig. 9D-E). The peaks in Dice coefficients are obtained for very similar threshold values in each image, and more broadly speaking, there are consistent ranges of values for which the Dice coefficient is close to its maximum value in all three images for both cliffs and ponds of both the Pléiades and the Sentinel-2 images. For cliffs in the Pléiades images, this happens for filtered Curvature between − 0.04 and − 0.02. There is a slight shift in NDWI for Khumbu (Fig. 9A, D), where there is a lot more dark, wetted debris, which tends to also have higher NDWI values (Fig. 5). For this image, the range of values maximizing the Dice for cliffs is also shifted to slightly higher values and the maximum Dice is lower than in the other images (Fig. 9B). Overall, the areas mapped based on Sentinel-2 data are in good agreement with those based on Pléiades but the cliffs do not stand out as clearly (Fig. 9C, F).

LSU
Applying LSU to the Pléiades images, we obtain maximum Dice coefficients ranging between 0.60 and 0.88 for ponds and between 0.12 and 0.49 for cliffs ( Fig. 9G-H). For Sentinel-2, the maximum Dice ranges between 0.55 and 0.62 for ponds, and between 0.16 and 0.49 for cliffs. Furthermore, the range of optimal NDWI and curvature values changes from one image to the other. We obtain low Dice values for Khumbu ponds, while for Langtang and Satopanth, very low optimal threshold values show that all pixels with any water content are categorized as ponds (Fig. 9G). Thus, for Satopanth and Langtang all the cliffs are mapped as ponds and cannot be characterized properly by an ice threshold (Fig. 9H). For the Khumbu image, the thinly debris-covered cliffs are mapped as debris and only the bare ice patches are categorized as cliffs (Fig. 9I). Similar results are observed for the Sentinel-2 images, although with some improvement in the Dice values for the cliffs (Fig. 9J-L).

LSU-s
For the LSU-s approach we obtain Dice coefficients ranging from 0.61 to 0.83 for ponds (Fig. 9M) and from 0.38 to 0.53 for cliffs (Fig. 9N) using the Pléiades images. For Sentinel-2, the maximum Dice coefficient is between 0.58 and 0.75 for ponds, and between 0.51 and 0.52 for cliffs. As with the SC results, there are ranges of values that maximize the Dice for each of the images, especially for Sentinel-2. For the cliffs, the main control on the results comes from the lower scale threshold. As long as the upper scale threshold is around 0.20, the change in Dice values is relatively small except for Khumbu where bright cliffs are prevalent (Fig. 9N, Q). The Khumbu Dice curve for cliffs is flatter for low values of lower scale threshold, highlighting that 'dark' cliffs have very low scale values in this image. The ponds are mapped after the cliffs here, and NDWI mapping results are slightly different from the ones obtained with the SC approach with a decrease of the maximum dice coefficients by 0.1 to 0.3 (Fig. 9A, D, M, P).

Comparison of methods
The highest Dice coefficients are obtained for cliffs and ponds mapped with the SC approach applied to the Pléiades scene (Figs. 9, 10,  Table 1). These values are high relative to other published approaches, but some noise remains and there are non-negligible false positive patches in the debris area compared to the manual delineation (Figs. 9C, 10B). The LSU method performs very poorly for both Pléiades and Sentinel-2 (Figs. 9I, 10C, G), but the LSU-s method has the best results for Sentinel-2 (Figs. 9P-R, 10H) with a relatively low number of false positives.
We apply the other existing methods to Langtang and Satopanth Pléiades scenes, and to the Khumbu Pléiades scene for the ABT approach (Table 1). For the images tested, the OBIA-m gives the same results as the manual delineation, but is faster due to the initial automated segmentation of the orthoimage in objects with similar spectral characteristics. The optimized Dice coefficients of the SST are similar from one scene to the other (Table 1), but are not obtained for the same slope threshold values (Fig. S5). The optimized Dice coefficients of the ABT vary between 0.24 and 0.64 depending on the scenes (Table 1, Fig. S6). These results can be improved when delineating the ponds first with the NDWI (Table 1). The ponds are also not mapped by the slope-based approaches but this has a limited effect on the classification (Fig. 11,  S3). The AST approach has a dice coefficient ranging from 0.40 to 0.55 obtained with the same initial parameters as in the original study (Herreid and Pellicciotti, 2018).

Cliff and pond coverage at each site
We report the count, total projected area and density (total feature projected area as a percentage of area) of cliffs and ponds obtained by applying the SC approach to the Pléiades and the LSU-s to the Sentinel-2 scenes as these are the best performing methods for the corresponding sensors ( Table 2). The cliff density observed in the Pléiades scenes is the  Table 1 Evaluation of different mapping methods for cliffs. In color are the methods tested in this study, including those requiring manual delineation (red), and the automated and semi-automated ones (purple). For the three methods developed in this study, we chose the same threshold values for all three Pléiades scenes. Transferability is a qualitative assessment characterizing the ability of the methods to produce consistently good performance from scene to scene. same for Khumbu and Satopanth (9.2%), but lower for Langtang (3.3%). The cliff density is higher than the pond density in all scenes, respectively 1.9, 3.7 and 13.1 times higher than the pond density on Langtang, Khumbu and Satopanth. The pond density is between 1.7 and 2.5% for Khumbu and Langtang, but only 0.7% for Satopanth.
There are 1.2 to 7.2 times fewer features detected with the LSU-s applied to the Sentinel-2 images compared to the SC applied to the Pléiades images, with this difference being higher for cliffs ( Table 2). The pond density is higher for Sentinel-2 than for Pléiades, but for cliffs, this is only the case on Langtang glacier (Table 2). There are a high number of cliffs smaller than 100 m 2 outlined using the Pléiades imagery, accounting for 11 to 22% of the total area of the outlined cliffs in one image (Fig. 12), so the area per feature is higher for the Sentinel-2 outlines.

Evaluation of methods
Here we assess our new methods along with other published approaches in terms of their performance, transferability and ease of implementation.

Comparison with past studies
Our results for both Khumbu and Langtang using the SC method with Pléiades data are consistent with previous reported cliff and pond density (planimetric area of cliffs or ponds divided by the total area of the debris-covered zone where the delineation was applied) values using manually delineated outlines. Pond density on Khumbu in May 2009 was 2.9%, which is consistent with the 2.5% value that we find in our Khumbu pre-monsoon scene, in the upper portion of the glacier (Watson et al., 2017a). For Langtang, the pond coverage peaks around 2% in June (Miles et al., 2017;Steiner et al., 2019), also consistent with our 1.7% density. We find a cliff density of 3.3% on Langtang, which agrees with previously reported values of 3.4% (+/− 0.9%) (Steiner et al., 2019). Previous studies have reported a 4% cliff density on Khumbu (Watson et al., 2017a), but it was calculated based on a digitization of cliffs as lines with no consideration of cliff area, this number is therefore not comparable with the 9.2% cliff density that we find. The pond density values obtained with the Sentinel-2 images are within the range of previously reported values for Khumbu and Langtang, and slightly higher than previously reported for the Langtang cliffs.

Performance of methods
Our three new methods map ice cliffs simultaneously with ponds, which are detected with high accuracy across scenes and methods (Fig. 9), and which we briefly discuss here. Ponds are best mapped on both Pléiades and Sentinel-2 images (Dice coefficients higher than 0.85 and 0.75 respectively) using a NDWI threshold. As such, using a 10 m instead of a 2 m resolution image has relatively little impact on the pond classification accuracy, although the pond density is 25 to 30% higher with Sentinel-2 than with Pléiades in our results, highlighting a possible overestimation of the pond density with Sentinel-2 likely due to mixed marginal pond pixels (Table 2, Fig. 9A, D, M, P). The suitability of Sentinel-2 images for mapping ponds was identified by Watson et al., 2018 and is promising for extensive and repeated mapping of ponds. The range of threshold values maximizing the Dice (0.05-0.2) is consistent from image to image and sensor to sensor, and a value of 0.1 is the best Table 2 Cliffs and ponds metrics obtained from the SC approach applied to the Pléiades images, and from the LSU-s approach applied to the Sentinel-2 images. The density corresponds to the total cliff projected area divided by the total debris-covered area.

Glacier
Khumbu 05  compromise to delineate ponds automatically in all images. It is worth noting that images with large extents of dark, wetted debris exhibit higher NDWI and Curvature values (Fig. 5), as in the Khumbu scene. This results in a shift of the optimal threshold value from 0.05 to 0.2 to avoid mapping parts of the debris as ponds (Fig. 9A, D). Considering our primary objective of mapping cliffs, we obtain satisfactory results for cliffs using the SC, with a Dice higher than 0.45 for the Sentinel-2 and higher than 0.5 for the Pléiades images. Performance is higher (Dice >0.6) for the Pléiades images when the extent of dark wetted debris is lower (Fig. 9B). Wetted debris patches can indeed be erroneously identified as cliffs due to their increased spectral curvature values. We obtain poor results for mapping cliffs using the LSU approach for both Pléiades and Sentinel-2. In this case, unless there is bright bare ice, the cliffs are identified as water or debris, which is the result of the mixed signals of wetted debris on an ice face. Cliffs with bright bare ice are more frequent on Khumbu, where we obtain relatively better results (Fig. 9H, I, K, L), while in the Langtang and Satopanth scenes almost all cliffs are categorized as ponds in the first step (Fig. 9G, H, J, K). We also found higher reflectance values for bare ice in the Khumbu scene than in the Langtang scene (Fig. S1B), which indicates that the ice end-member we used is not representative enough to map the Langtang and Satopanth cliffs. This could come from changes in englacial debris concentration, and differences in orientation, density or foliation and these variations in the bare ice spectra from glacier to glacier and from scene to scene make the LSU approach difficult to apply in general. The LSU-s does not depend on the bare ice spectrum and uses the additional information from the Sentinel-2 re-sampled 20 m bands to maximize the Dice coefficient around 0.5 for each of the three scenes, making it the best method to map cliffs with this sensor (Fig. 9Q). In some cases of extensive darker debris patches (such as on Khumbu and locally on Satopanth), the scale filtering may not be sufficient to outline all the cliffs from the Sentinel-2 image, resulting in an underestimation of the total cliff area and density ( Table 2). The LSU-s method does not perform as well for Pléiades as for Sentinel-2, especially in the Langtang and Satopanth scenes (Fig. 9N), for which there are false positive results for locally darker debris patches that are spectrally similar to thinly debris-covered cliffs (Fig. 9O).
We obtain varying results from the other cliff mapping approaches. A SST approach can produce a reasonably accurate cliff distribution with Dice coefficients between 0.43 and 0.51 for Langtang and Satopanth, but requires recalibration for each scene since the optimal slope threshold varies from 29 • for Langtang to 40 • for Miage (Fig. S5, Reid and Brock, 2014). This recalibration can be avoided by using an automated method such as AST with its original parameters (Herreid and Pellicciotti, 2018), which also produces slightly higher Dice coefficients between 0.40 and 0.55 (Table 1, Fig. 11, S7). The dependence of this method on the DEM spatial resolution is partially compensated by extending the cliff ends, where the cliff is narrowest and the slope susceptible to be saturated (Herreid and Pellicciotti, 2018). The ABT (Anderson et al., 2019b) requires particular conditions to produce reasonable results (Table 1, Fig. 11F). We optimized the brightness threshold and the window size of this method for the three Pléiades scenes and get a Dice of 0.64 for Langtang and 0.49 for Satopanth, but for Khumbu the Dice is lower than 0.25 (Table 1, Fig. 11F, Fig. S6). This approach works relatively well for Langtang since the cliffs are generally much darker than the surrounding debris, but performance suffers in the Khumbu scene and parts of the Satopanth scene because 1) the debris can be as dark as the wetted cliff surfaces and 2) many cliffs have very bright bare ice patches (Fig. 11, (Anderson et al., 2019b). A second brightness threshold would be necessary to identify these cliffs. This method also only maps cliffs without taking the ponds into account, which results in most of the ponds being also mapped as cliffs due to their similar brightness values (Fig. 11, S7). The ABT therefore benefits from pond mapping with NDWI prior to its application (Table 1), producing Dice coefficient values nearly as high as our SC approach for the Langtang scene.

Transferability of methods
We define transferability as the capability of a method to produce consistently good performance for different sites or input data. Transferability is therefore important to enable widespread use of a method. We assess each method's transferability by examining its performance across three scenes with different climatic and geomorphological settings.
We consider the SC approach to be transferable for Pleiades data as the Dice is maximized for the same range of curvature thresholds (− 0.04 to − 0.02). Based on the optimized Dice coefficients, a curvature value of − 0.03 is a good compromise to map cliffs in all scenes (Fig. 9B). Similarly the LSU-s appears to be transferable for Sentinel-2 data, as an upper filtered scale threshold of 0.2 (log scale) and a lower filtered scale threshold of − 0.2 (log scale) maximize the Dice coefficient of mapped cliffs around 0.5 in all three Sentinel-2 scenes while a NDWI threshold of 0.1 maximizes the Dice coefficient of ponds. Both slope-based approaches are also transferable, with Dice coefficients between 0.40 and 0.55 for all the scenes where they were applied, including on an Alaskan glacier for the AST (Herreid and Pellicciotti, 2018). This is not the case for the ABT (Anderson et al., 2019a(Anderson et al., , 2019b, which does not work when the debris is as dark or darker than the cliffs. Our assessment of transferability encompasses climatic and geomorphological settings of Nepal and Indian Himalaya, including images taken in the pre-monsoon, monsoon and post-monsoon season. We expect that the approaches classified as transferable in this study should also be transferable to other glaciers in other mountain ranges, but accounting for variable lithology may require some adjustment to the LSU-s as well as to the SC approach, and the transferability beyond the Himalaya needs to be evaluated. Furthermore, application of any cliff mapping method is limited to cloud-and snow-free debris-covered areas, which need to be accurately outlined prior to the mapping. This was done manually in this study since we focused on a limited number of glaciers, but is an important prerequisite to be considered when trying to map cliffs at the regional scale in an automated way. For this purpose, the debris-covered area of a glacier can be mapped automatically using glacier outlines and a band ratio Kraaijenbrink et al., 2017), which has led to the recent release of global maps of supraglacial debris-cover extents (Herreid and Pellicciotti, 2020;Scherler et al., 2018). Snow-covered areas can be mapped automatically using an adaptive normalized difference snow index (NDSI) (e.g. Girona-Mata et al., 2019;Rastner et al., 2019) and automated approaches exist to mask cloud-covered areas along with zones with deep shadows (Chen et al., 2013;Miles et al., 2017;Zhu et al., 2015).

Efficiency of methods
Manual delineation of cliffs over the entirety of a glacier is very timeconsuming. For one or a few glaciers, applying the OBIA-m approach is a worthwhile improvement as it gives results equivalent to a full manual delineation while being several times faster. Indeed, once the segmentation has been applied, the operator only needs to select the clusters comprising cliffs and ponds and adjust the outlines (Table 1, Mölg et al., 2020Mölg et al., , 2019. For a larger domain, using a (semi-)automated and transferable approach can be advantageous as in theory no recalibration needs to be conducted. The AST approach detects cliffs iteratively using a set of precalibrated parameters and is as such entirely automated. However it is computationally expensive even for one glacier regardless of its setting (Herreid and Pellicciotti, 2018); we could only run it for 5 m resolution degraded DEMs (Table S3). This method uses the abundance of cliffs to optimize the slope threshold so it also requires, contrary to the other approaches, very accurate outlines of the debris-covered area, otherwise it identifies steep surrounding topography as large cliffs, which can impact the slope threshold identifying ice cliffs (Herreid and Pellicciotti, 2018). The SC and LSU-s approaches are both transferable at the scale of the Himalaya, semi-automated and computationally much faster than the AST. The SC approach is the faster of the two, as it does not require the spectra of each pixel to be decomposed (Table S3). The SST and ABT are also computationally efficient approaches (Table S3) but less transferable than the SC or the LSU-s. When applying the SST and ABT approaches to glaciers with different surface characteristics, or the SC and LSU-s outside of the Himalaya, recalibrating the parameters of the method is necessary. This implies that for each scene a dataset of 20-30 cliffs and ponds needs to be manually outlined. This is feasible when dealing with a few satellite images, but becomes difficult when dealing with many glaciers with very different surface characteristics. Furthermore, this requires fine spatial resolution and high-quality spectral data for which manual delineation is possible.

Definition of ice cliffs
The different mapping methods have each been developed to target a specific definition of ice cliffs. While the slope-based approaches assume that the main characteristic of a cliff is its steepness (Herreid and Pellicciotti, 2018;Reid and Brock, 2014), the multispectral approaches, including the new approaches developed in this study and Anderson et al. (2019b), assume that surface reflectance characterizes cliffs and map them based on their spectral signature, without taking into account the slope. Other approaches such as manual delineation may use a combination of slope and spectral signal criteria, for example mapping first with the spectral information and refining the outlines by discriminating ice cliffs and shallow-sloping ice patches (Kraaijenbrink et al., 2016;Steiner et al., 2019;Watson et al., 2017a). These variations in the definition used may lead to mismatch between the resulting datasets, even though they should broadly agree, as the surface exposure of ice within the debris-covered area implies that the debris has been displaced from this zone, which requires the presence of a slope, sometimes combined with the action of water (Moore, 2018).
Both spectral and slope-based delineation methods outline zones that do not look like cliffs and that an operator would classify as false positive objects, but that contain the spectral signature or high slope exhibited by cliffs, and sometimes both (Fig. 11). The semi-automated approaches applied to the Pléiades and, to a lesser extent, to the Sentinel-2 images, identify many small ice patches that would be impossible for an operator to map in a consistent way. These features considerably increase the number of mapped features, and slightly increase the total cliff density (Figs. 11, 12, Table 2). They are small enough that they are often not represented by the slopes calculated from the DEMs, even if they are steeper than the debris. Depending on the spatial resolution, cleaning the results of the semi-automated approaches with a size threshold may be tempting to produce cliff outlines closer to what an operator sees, but their spectral signature does suggest that some ice is exposed or thinly debris-covered at these locations. Thus, these findings highlight a potential bias of the manual delineation, which only works for large features that have an identifiable shape and can cast a shadow, thus having a 'cliff' appearance. This is also true for the OBIA-m approach, for which it is difficult to apply an initial multi-scale segmentation to get the outlines of both big cliffs and small ice patches directly. Moreover, these small areas of shallow-sloping ice have seldom been counted as ice cliffs in previous studies (Anderson et al., 2019b), nor have they been accounted for when modeling the melt of debris-covered glaciers, but may account for substantial mass loss, as if they are partially or thinly debris-covered, they are susceptible to act as areas of enhanced melt (Evatt et al., 2015;Fyffe et al., 2020;Ostrem, 1959;Reid and Brock, 2010). Thus, it seems important from a melt-modeling perspective to investigate the distribution and density of these small areas of shallowsloping ice. However, their limited size and slope most probably prevents them from surviving due to debris redistribution, unless they develop into an actual cliff.
It is apparent that there is no broadly agreed definition of ice cliffs that can be used to tailor a method to their detection, and this seems an important step towards large-scale classification efforts. A plausible definition of an ice cliff could be a feature within an otherwise debris-covered area containing ice in its spectral and/or thermal signature (Herreid and Pellicciotti, 2018), but large and steep enough to survive over at least one melt season. However, the presence of the smaller areas of exposed ice detected in this study opens a new perspective. These patches should be investigated to understand their nature and distribution, their evolution (can they act as seeds for larger cliffs, or are they ephemeral features), and their contribution to glacier melt.

Future applications
A point to bear in mind when mapping steep features like ice cliffs from satellite imagery is that the final planimetric area will be impacted by possible shadowing and steep cliffs being concealed from the satellite view angle. Shadowing leads to an overestimation of the cliff planimetric area for all approaches based on the multispectral images. The pixels in the shadow are indeed categorized as ice cliff due to the low brightness and scaling factor values for the ABT and the LSU-s. The SC also maps shadows as cliffs due to its high sensitivity to small shifts in reflectance between the different band values. However, shadowing only affects steep north-facing cliffs (in the Northern Hemisphere), resulting in a limited impact on the cliff distribution at the glacier scale. All approaches are also equally affected by 1) the satellite viewing angle that results in cliffs with overhangs or steeper than the satellite viewing angle being concealed, 2) the spatial resolution of the image that prevents small steep features from being mapped. It is therefore key to use images with high satellite and sun elevation angles. In this study these angles where higher than 70 • and 55 • respectively for all images.
This being said, the transferability and efficiency of the semiautomated SC and LSU-s approaches opens up the possibility to systematically map ice cliffs on debris-covered glaciers at other sites, including in other mountain ranges, to study the distribution and evolution of ice cliffs across scales larger than an individual glacier or catchment. Although our confidence in the results of these approaches is already high, they could still be improved by combining spectral and slope-based approaches. This could be particularly important to discriminate between ice cliffs and seasonal snow patches, which were not a problem for our three scenes. Combining the spectral approaches with the OBIA segmentation could also be something to explore, but would require the development of a multi-scale segmentation to account for cliffs of all shapes and sizes. Finally, adding a second brightness threshold to the ABT to account for bare ice portions of ice cliffs would be a useful improvement for this method. LSU-s enables the mapping of cliffs at the large scale using freely available Sentinel-2 data to study the variations of ice cliff coverage from glacier to glacier, region to region and potentially over time. Due to the coarser spatial resolution of the images, the Sentinel-2 cliff outlines from the LSU-s approach can be interpreted as a first-order map of the larger cliffs. Indeed, the Sentinel-2 LSU-s approach only maps the larger features and the resulting outlines are relatively crude (Table 2, Figs. 9R,10,12). In the case of extensive dark debris patches like on Khumbu, and to a lesser extent on Satopanth, the filtering of the scale may in some places not be strong enough for cliffs to contrast with their background, resulting in a high number of false negatives and an underestimation of the cliff density with regards to the outlines from SC applied to the Pléiades (Table 2). However, the larger cliffs are mapped and provide a close estimate of the total planimetric cliff area for the glacier. This is therefore a promising approach to assess ice cliff distribution across Himalayan debris-covered glaciers, given the availability of Sentinel-2 imagery covering large areas with high repeat frequency.
On the other hand, using fine spatial resolution imagery with the SC approach gives a highly detailed map of the ponds, ice cliffs, and smaller ice patches across the debris-covered area of a glacier. Supraglacial cliffs and ponds have significant implications for the glacier mass balance, and the effect of smaller bare ice patches is probably also important. There are already existing independent melt models of supraglacial ponds and cliffs (Buri et al., 2016a(Buri et al., , 2016bMiles et al., 2018Miles et al., , 2016 and the accurate mapping of all these features is a first and decisive step to quantify their impact on melt at the scale of a glacier or a catchment using a glacier melt model integrating these features.

Conclusion
This study has explored the potential of three new semi-automated, objective approaches to enable accurate mapping of cliffs and ponds from multispectral images within one glacier and from glacier to glacier. We also assessed their performance when applied to both fine spatial resolution, commercial (Pléiades, 2 m resolution, 4 bands) and freelyavailable, coarser spatial resolution images with more spectral information (Sentinel-2, re-sampled to 10 m resolution, 10 bands). Our main goal was to develop a method that can be applied across scales, and evaluate its potential in relation to existing approaches. The transferability of the new methods is assessed using three scenes with different climatic settings, elevations, debris types, cliff and pond characteristics, in Central and Eastern Himalaya.
We find that the Spectral Curvature (SC) approach, which uses Spectral Curvature to map cliffs after having removed ponds using NDWI, outperforms all other methods for the fine spatial resolution Pléiades data. This approach is transferable, semi-automated, computationally efficient, and is the best performing for each of the three scenes. The Linear Spectral Unmixing with scale (LSU-s) approach, an adapted linear spectral unmixing method, consistently produces the best performance for Sentinel-2 images by taking advantage of their additional spectral information. It is transferable and computationally efficient, making it the most appropriate method to map cliffs with this type of sensor. The two methods classify cliffs that occupy between 3 and 9% of total planimetric area depending on the site and scene. We thus propose both as a way forward for either mapping repeat cliffs inventories at the scale of select glaciers to study processes of cliffs formation and evolution (SC approach); or to establish large scale inventories of cliffs features across distinct regions (LSU-s), to enable novel analyses of much larger cliffs populations than tested to date, and quantify their characteristics, distribution and patterns over large scales.
Other existing methods all have potential for distinct applications, but seem less suitable for the two goals above. The only other automated method that is also transferable despite lower performances than the SC approach, the AST, is very slow from a computational standpoint. When seeking to identify cliffs for glaciers with consistent surface characteristics, the semi-automated SST and ABT approaches may be considered, even though they do not score as high or as not as transferable as the SC approach for the Pléiades scenes. The ABT may fail for glaciers with cliffs as bright as or brighter than the surrounding debris, but otherwise these methods generally perform well and are computationally very efficient. All of these methods are faster and more objective than a timeconsuming manual delineation, even if an initial OBIA segmentation is used. None of the other existing approaches can be applied to the coarser Sentinel-2 data: the slope-based approaches require elevation data and are therefore not applicable to the Sentinel-2 images, while the coarse resolution of this sensor prevents the application of the other semiautomated or manual approaches to map cliffs.
Our results interestingly highlight a high number of smaller, sometimes shallow-sloping ice patches, which had not been previously accounted for in the assessment of ice cliffs and their contribution to melt. These smaller features reveal a bias in the manual delineation of cliffs towards identifying large features, and reflect the need for a clear definition of 'ice cliffs'. Although representing a small fraction of the total planimetric cliff area, these small ice patches could have nonnegligible impacts on the melt of debris-covered glaciers.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.