An Automated Approach for Estimating Snowline Altitudes in the Karakoram and Eastern Himalaya From Remote Sensing

The separation of fresh snow, exposed glacier ice and debris-covered ice on glacier surfaces is needed for hydrologic applications and for understanding the response of glaciers to climate variability. The end-of-season snowline altitude is an indicator of the equilibrium line altitude of a glacier and is often used to infer the mass balance of a glacier. Regional snowline estimates are generally missing from glacier inventories for remote, high-altitude glacierized areas such as High Mountain Asia. In this study, we present an automated, decision-based image classification algorithm implemented in Python to separate snow, ice and debris surfaces on glaciers and to extract glacier snowlines at monthly and annual time steps and regional scales. The method was applied in the Hunza basin in the Karakoram and the Trishuli basin in eastern Himalaya. We automatically partitioned the various types of surfaces on glaciers at each time step using band ratios combined with topographic criteria based on two versions of the Shuttle Radar Topography Mission elevation dataset. Snowline altitudes were extracted on a pixel-by-pixel basis using a “buffer” method adapted for each elevation dataset. Over the period studied (2000 to 2016), end-of-the-ablation season annual ELAs fluctuated from 4,917 m a.s.l. to 5,336 m a.s.l. for the Hunza, with a 16-year average of 5,177 ±108 m, and 5,395 m a.s.l. to 5,565 m a.s.l. for the Trishuli, with an average of 5,444 ± 63 m a.s.l. Snowlines were sensitive to the manual corrections of the partition, the topographic slope, the elevation dataset and the band ratio thresholds particularly during the spring and winter months, and were not sensitive to the size of the buffer used to extract the snowlines. With further refinement and calibration with field measurements, this method can be easily applied to Sentinel-2 data (5 days temporal resolution) as well as daily PlanetScope to derive sub-monthly snowlines.


INTRODUCTION
Identifying various surfaces on glaciers (fresh snow, clean glacier ice, and supra-glacial debris cover) and extracting glacier snowlines are needed for glacier mass balance calibration and validation, as demonstrated in a growing body of literature (Rabatel et al., 2005(Rabatel et al., , 2008(Rabatel et al., , 2017Gardelle et al., 2013;Huss et al., 2013;Kienholz et al., 2017;Barandun et al., 2018). Snowline altitudes (SLAs), when measured at the end of the melt season, represent the equilibrium line altitude (ELA) of a glacier (Meier, 1962), which is an indicator of seasonal/annual glacier mass balance and its response to climatic variability (Paterson, 1994). ELA is the altitude at which annual glacier mass balance is zero, and is inferred from direct mass balance measurements using the traditional "glaciologic method" with stakes and pits (Meier and Post, 1962). However, in remote areas of High Mountain Asia (HMA), access to glaciers, especially to their accumulation areas, is limited by rugged terrain and difficult logistics. Such measurements are sparse in global glacier databases such as the Randolph Glacier Inventory (RGI) (Pfeffer et al., 2014) or the Global Land Ice Monitoring from Space (GLIMS) database (Raup et al., 2007). In HMA, only a handful of glaciers have been surveyed systematically for mass balance (Bolch et al., 2012;Azam et al., 2018). Field-based ELAs for HMA are reported in Wagnon et al. ( , 2013 for Chhota Shigri and Mera/Pokalde glaciers in the Khumbu region and Acharya and Kayastha (2019) for Yala glacier in Trishuli basin in Nepal. However, long-term records are limited. The relationship between ELA and the end-of-summer SLA can be applied to infer annual mass balance from remotely sensed snowlines in such areas with limited field-based measurements as well as for missing years (Rabatel et al., 2012). However, the SLA/ELA relationship is not straightforward, as it is complicated by the presence of patches of firn (old snow) and superimposed ice in the accumulation area of glaciers due to melting of snow and refreezing of water (Paterson, 1994). Superimposed ice causes the ELA to be situated below the end-of-season SLA (Llibutry, 1998), though this may not always be the case . Conversely, late in the ablation season, under climate warming conditions, the snowline may retreat beyond the firnline, but optical satellite imagery in this case detects the firnline rather than the snowline.
Separating snow and ice surfaces also allows inferring the accumulation-area-ratio (AAR) of a glacier, defined as the ratio of the accumulation area to the entire area of the glacier (Meier, 1962). AAR fluctuations at the end of a hydrological year are an indicator of glacier mass balance at local or regional scales (Dyurgerov et al., 2009). ELA/AAR have been used to estimate glacier mass balances at regional scales using remote sensing methods (Kulkarni, 1992;Dyurgerov, 1996;Rabatel et al., 2008). The "AAR/ELA method" developed by Kulkarni (1992) for the Western Himalaya allows inferring glacier mass balance at regional scales from satellite imagery (Kulkarni et al., 2004). A "template" method based on the relationship between AAR and glacier mass balance was developed by Dyurgerov (1996), and Khalsa et al. (2004). However, the wide application of these methods in HMA is limited by lack of field-based SLA/ELA measurements needed to develop the AAR-mass balance relationship.
In lack of any field-based measurements, ELA can be inferred using indirect methods such as Toe-to-Headwall Altitude Ratios, the Area-Altitude-Ratio, the Area Altitude Balance Ratio and the Area Altitude Balance Index methods and Median Elevation. These methods were reviewed in detail in other studies (Kaser and Osmaston, 2002;Osmaston, 2005); here we only highlight a few of these. For example, the Toe-to-headwall Altitude Ratio method was used to reconstruct ELAs since the Little Ice Age in the Nepalese Himalaya (Kayastha and Harrison, 2008). The median elevation approach is fast, and only requires glacier outlines and a digital elevation model (DEM). Braithwaite and Raper (2009) used this method to estimate ELAs using glacier inventory data from the 21st century and considered the uncertainties to be acceptable (±82 m) (Braithwaite and Raper, 2009). However, they point out that ELAs are underestimated when glaciers are in a state of negative mass balance. King et al. (2017) also used the median elevation method for glacier mass balance in central Himalaya in the last decade and likewise pointed out potential underestimates due to negative mass balance of glaciers in this area (Bolch et al., 2012;Brun et al., 2017). The topographic map-based method estimates ELAs from the inflections of contours on topographic maps (Leonard and Fountain, 2017). Zhao et al. (2016) found a good correlation between ELAs derived using this method and median glacier elevations for Chinese glaciers (R 2 = 0.92). Pellicciotti et al. (2015) assumed an AAR of 0.66 for the Langtang region in Central Himalaya to infer regional ELAs for their melt model. However, this is problematic for debris covered glaciers where AAR values are lower (0.2-0.4) (Kulkarni, 1992). Modeling approaches rely on empirical curves based on temperature and precipitation gradients to estimate modern-day steady-state ELA values (ELA 0 ), as described in Ohmura and Boettcher (2018). Using a similar method, for HMA, Fujita and Nuimura (2011) estimated the "ideal" ELAs for present glaciers based on an energy balance model with downscaled climate variables, validated with several field-based ELAs. All these methods each have advantages and disadvantages, and choosing one or the other may depend on the data available as well as the scale of analysis. Median elevation might be an appropriate proxy for ELA for global applications or when past reconstructions are desired; however, at smaller scales, remote sensing methods might be more accurate. Topographic maps in HMA are not readily available or are outdated, so the topographic method should be used with a DEM. The wide application of empirical methods is hampered by limited climate data at regional or local scales in HMA.
Satellite imagery provides opportunities to improve regional ELA estimates using high temporal and spatial resolution data, provided that systematic approaches are developed. While snow and ice boundaries are well visible on satellite imagery acquired with good contrast or appropriate instrument gains to minimize pixel saturation, distinguishing between snow, firn and ice on the glacier surface is challenging using conventional methods. Huss and Hock (2015) parameterized their melt model by prescribing "firn" above the median glacier elevation easily available in RGI datasets and "bare ice" below it. Thakur et al. (2017) used SAR imagery to separate glacier facies for a small sub-basin in the North West Himalayas, but the large scale application of this method might be limited by the availability of the SAR data and the extensive data processing. Kienholz et al. (2017) defined the limit between ice/firn and snow in Alaska on the basis of Landsat false color composites and a DEM and manually extracted multitemporal SLAs. Huss et al. (2013) and Barandun et al. (2018) extracted SLAs based on ground photographs and a DEM using an innovative method, but this was only applied to a few glaciers. Zhang and Kang (2017) tracked the evolution of snowlines over two decades in the Pamir using Landsat, but no detail is given on the process used for the actual extraction of the SLAs. Guo et al. (2014) extracted SLAs as a single value per year for the Western Himalaya using surface albedo calculated from atmospherically and topographically corrected images overlaid on topographic map for each year studied. In the French Alps and the Andes, SLA/ELA was extracted mainly by manual digitization on aerial photographs or color composites of satellite images (Rabatel et al., 2005(Rabatel et al., , 2008(Rabatel et al., , 2012(Rabatel et al., , 2017. While field-based validation showed good agreement with manually derived remote-sensing ELAs in these studies, manual digitization is time-consuming, and is not applicable over large areas. Klein and Isacks (1999) used spectral unmixing of satellite images to separate the ablation and accumulation areas of glaciers and to extract ELAs on two tropical glaciers and found this to be superior to band ratio techniques.
Despite recent advances in partitioning snow and ice and estimating SLA/ELA, existing optical remote sensing methods are hampered by deep shadows on the glacier accumulation area due to steep topography, icefalls and crevasses. Moreover, the precise date of end-of-ablation season is difficult to define, especially in monsoon-dominated catchments of the eastern Himalaya where snow accumulation and ablation occur concomitantly (Ageta and Higuchi, 1984;Thayyen and Gergan, 2010). SLA extraction methods are not standardized, and guidelines and recommendations such as those established for glacier mapping within the GLIMS project (Racoviteanu et al., 2009) are missing.
In this study, we present an automated method using Landsat imagery to separate various surfaces: exposed glacier ice, snow on ice, snow on land and/or on debris, debris covered ice from "other" (bare, non-glacier/non-snow terrain) and to extract SLAs/ELAs at monthly and annual time scales in two glacierized areas of HMA -the Karakoram and eastern Himalaya. Here we are not distinguishing snow from firn, but rather we consider firn to be included in the snow on ice class. This study was developed to parameterize a temperature index melt model within the framework of the Contribution to High Mountain Asia Runoff from Ice and Snow collaborative project (CHARIS 1 ) (Armstrong et al., 2019). Here we describe the methodology used in the cited study to partition the surface types. The two-step decision-based procedure was implemented in Python, allowing loop-processing series of satellite images and extracting the SLA automatically at each time step. We constructed monthly SLA time series in the Hunza (2013) and Trishuli (2016) subset areas, and estimated ELAs from 2000 to 2016. We automated part of the post-classification correction of problematic areas, and assessed the impact of these corrections on SLA estimates, as an effort toward standardizing snowline mapping procedures using remote sensing and applying them at regional scales in HMA and beyond. 1 nsidc.org/charis

Study Area
This study focuses on two climatically distinct regions of HMA: a subset of the Upper Indus basin in the Karakoram (Hunza sub-basin, 4,861 km 2 and mean elevation 4,863 m a.s.l.) and a subset of the Trishuli basin in the eastern Himalaya (6,086 km 2 , mean elevation 4,699 m a.s.l.) ( Figure 1A). For simplicity, in the paper we refer to these study areas as the "Hunza" and "Trishuli, " and subsets of the two regions. The first study area ("Hunza") is part of the Northern Areas of Pakistan centered on Shimshal Valley, and includes the fast moving surging glacier Khurdopin. The growing terminus lake of Kurdopin glacier poses concerns for hazards and Hispar glacier, which has been undergoing active surges since 2013 (Rashid et al., 2018) (Figure 1B). Climatically, the region is mostly arid, and is primarily influenced by mid-westerly winds originating from the Mediterranean and Caspian Sea regions (Bookhagen and Burbank, 2010). This area receives maximum precipitation as snow in the winter and spring (Fowler and Archer, 2006). Glaciers in the Karakoram are considered of "winter-accumulation-type" (Benn and Owen, 1998;Thayyen and Gergan, 2010) and have been mostly stable or growing, a condition known as the "Karakoram anomaly" (Hewitt, 2005;Minora et al., 2016), most recently attributed to anomalous summer cooling (Forsythe et al., 2017). Fieldbased glacier ELAs are almost non-existent in this rugged area. Remote sensing regional ELA estimates are sparse and variable, ranging from 4,845 m a.s.l. (Scherler et al., 2011) to 4,300-5,500 m a.s.l. based on Landsat imagery (Khan et al., 2015). Estimates vary considerably from one watershed to another due to a strong east-west gradient in precipitation patterns (Mukhopadhyay and Khan, 2016).
The second study area is part of the Narayani basin the Central-Eastern Himalaya of Nepal ( Figure 1C). Climatically, this region is influenced by the south-west Asian summer monsoon circulation system (Yanai et al., 1992;Benn and Owen, 1998). Moist air masses from the Bay of Bengal interact with the topography of the Himalaya and Tibetan Plateau (HTP), causing maximum precipitation on the southern slopes of the Himalaya during the summer months (June-September) (Shrestha, 2000;Bookhagen and Burbank, 2006). Glaciers in this monsoon-influenced part of the Himalaya are of "summer-accumulation" type (Ageta and Higuchi, 1984;Thayyen and Gergan, 2010). The headwaters of glaciers in the Trishuli basin originate from China; our study area spans both the southern and the northern slopes of the Himalaya ( Figure 1C). ELA measurements are scarce in this area as well, and most estimates come from indirect methods (Kayastha and Harrison, 2008) or modeling approaches (Acharya and Kayastha, 2019).

Data Sources
The main satellite data sources for this study are from the Landsat 7 Enhanced Thematic Mapper (ETM +) and the Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS), obtained from the USGS. These two sensors have been acquiring imagery in the visible, near infrared, short wave infrared and thermal bands of the electromagnetic spectrum since 1999 and 2013, respectively [see USGS (2016) for more details]. Monthly Landsat OLI scenes from 2013 (for Hunza) and 2016 (for Trishuli) were used to determine seasonal SLAs fluctuations and to estimate the approximate date of end-ofablation season in each region ( Table 1). The years were selected based on maximum number of cloud-free images with good contrast over snow and ice in each area. The year 2013 had snowy and cloudy conditions in the Trishuli, making it hard to obtain suitable images, so we could not use the same year as for the Hunza for the seasonal analysis. Annual Landsat ETM+ scenes (2000 and OLI (2013 and 2016) were used to map ELAs and to assess ELA fluctuations over the 16-year record. These scenes were selected using a 4-month window in each area (August-October for the Hunza and September-December for the Trishuli) based on the seasonal SLAs and on previous literature (Thayyen and Gergan, 2010). We used digital numbers (DNs) rather than atmospheric reflectance since the latter was not available at the onset of our study.
To extract SLAs, we used elevations from the Shuttle Radar Topography Mission (SRTM), which acquired near-global data in February 2000. We used two SRTM versions: (a) the TABLE 1 | Landsat 7 ETM+ and Landsat 8 OLI remote sensing data used in this study, with the optimized thresholds applied to band ratios and topographic criteria to partition the snow, ice, and debris surfaces. DN ranges represent the (minimum and maximum) values used to map ice and debris cover. Thresholds were used for: A: annual analysis (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) (Jarvis et al., 2008), with void-filling procedures described in Reuter et al. (2007) and (b) the NASA SRTM Version 3.0 (SRTMGL1), at 1-arc second (∼30 m); this was released in 2015 and was voidfilled using elevation data mostly from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model 2 (GDEM2) (NASA-JPL, 2013). The SRTM global DEM targeted a vertical accuracy of ±16 m (Rabus et al., 2003) but in rugged terrain, accuracy can be considerably lower (Berthier et al., 2006;Fujita et al., 2008). In a recent study, Mukul et al. (2017) evaluated the 90 and 30 m versions of the SRTM DEMs and reported an uncertainty (as vertical root mean square error, RMSE z ) of 47.2 and 23.5 m, respectively, based on ground control points in the Himalaya. Bias correction improved the accuracy of both of these DEMs in the cited study. In our study, we accounted for DEM errors in our SLA/ELAs on the basis of the uncertainties reported by Mukul et al. (2017).

Surface Partitioning: Clean Glacier Ice, Snow on Ice, Snow on Land, and Debris Cover
The various snow-ice-debris surfaces were separated using a multi-criteria method implemented in Python as a series of conditional statements, based on band ratios and topographic criteria. Thresholds for various band ratios, elevation, slope, and temperature criteria were selected based on visual inspection of Landsat color composites, previous literature and a priori knowledge. Thresholds for the monthly series were adjusted since they varied seasonally; thresholds for the annual time series were standardized for all images. The criteria used, along with their thresholds, are reported in Table 1; an overview of the surface partition process is presented in Figure 2. The output of surface partition algorithm was a four-class raster when no seasonal snow was present (Figures 2A,B) or a five-class raster for the months when seasonal snow was present. In a preliminary step, we used single band ratios to mask areas under shadow, clouds, and surface water, which pose challenges for the semi-automated band ratio methods (Racoviteanu et al., 2009). Deep shadows are common in mountainous areas on steep slopes especially in the winter due to the low sun angles in the morning around the time of acquisition of the Landsat scenes (∼05:36 GMT or ∼10:36 local time for this particular area). Shadows appear darker than other surfaces in the near-infrared, so we used this wavelength to map shadows from steep terrain, including cast shadow from clouds, which are also an issue (Racoviteanu and Williams, 2012). We applied a threshold to band 5 (ETM+) or band 6 (OLI), which we checked and adjusted for each month in the time series to obtain shadow "masks." Tests for the cloud mapping techniques based on the FMask algorithm (Zhu et al., 2015;Zhu and Helmer, 2018) failed to distinguish between snow and clouds in this area. Therefore, in this study, clouds were mapped using a single band threshold. At near infrared (ETM+ band 5: 1.55-1.75 µm) and shortwave infrared wavelengths (OLI band 6: 1.57-1.65 µm), snow and ice surfaces are dark and clouds are bright, making these two types of surfaces distinguishable. We carefully adjusted the DN thresholds to avoid misclassifying illuminated bare terrain as clouds (Table 1), even though we consider that this would have little/no effect on our surface classifications. Shadow and cloud raster masks were assigned "NoData" values in the final maps ( Figure 2B). Vegetated areas were mapped using the Normalized Difference Vegetation Index (NDVI), defined here as the difference between visible and near infrared bands, i.e., NDVI TM = band 3 − band 4 band 3 + band 4 and The NDVI algorithm results in a raster with values from −1 to 1. This was used to map both the vegetation areas (negative values) and water bodies (positive values) ( Table 1), which were excluded from the potential debris cover map and assigned to the "bare land" class. Mapping the water bodies (lakes and rivers) was only introduced in the Trishuli subset as an improvement of the method, though this did not affect the surface partition on glaciers.
Glacierized areas (clean glacier ice and snow) were mapped using the standard semi-automated "band ratio" technique, which is robust and easy to apply over large areas, and widely used by the glaciological community (Racoviteanu et al., 2009;Bolch et al., 2010;Paul et al., 2013). This technique takes advantage of the spectral difference between snow and ice surfaces and other types of terrain in the visible and near-infrared parts of the electromagnetic spectrum. We calculated the ratio of DNs using two bands ( band 4 band 5 for ETM+ and band 5 band 6 for OLI) and thresholded the resulting raster to obtain a binary map (1 = "snow/ice, " 0 = "other") ( Figure 2A). The segmentation of band ratio images using raw DNs has been found to be superior over band ratios using atmospherically corrected images when cast shadow was present (Paul et al., 2002). We used a band ratio threshold of 1.5 for all our images ( Table 1). We tested various thresholds from the literature (1.2-2) on the basis of false Landsat color composites (ETM+ bands 5, 4, 3 and OLI bands 6, 5, and 3) and we conducted a sensitivity analysis to these thresholds (section "Sensitivity Analysis"). The lower threshold of 1.5 compared to the ones used for glacier mapping (Andreassen et al., 2008;Racoviteanu et al., 2015) allowed mapping the full extent of snow and ice for each time step. The minimum snow and ice area for any given year (2013 for Hunza, 2016 for Trishuli) was considered to represent the glacier surface for that year and was used as a "glacier mask" for the subsequent steps described below.
Snow on ice in the accumulation area of glaciers was separated from exposed glacier ice in the ablation area by thresholding the near-infrared band (ETM+ band 4 at 0.77-0.90 µm and OLI band 5 at 0.851-0.879 µm), along with an elevation and a slope criteria (Figure 2A). At these wavelengths, ice and snow have different brightness temperatures, making them distinguishable from each other. Due to differences in illumination conditions, FIGURE 2 | Surface partition workflow using Landsat bands and a DEM, and its outputs: (A) band ratios for separating exposed glacier ice, snow on ice and debris covered ice from bare land; and (B) band ratios for delineating clouds, shadows, vegetation, and water bodies.
the thresholds varied depending on each image. To pick a threshold, for OLI band 5, for example we defined regions of interest (ROI) on exposed glacier ice surfaces and extracted summary statistics from OLI band 5 over 1 year. On the basis of mean statistics, we tested various thresholds to compute an "optimized" ratio for each month, which was adjusted for each time step based on visual inspection ( Table 1). The sensitivity of the resulting SLAs to these thresholds will be discussed later (section "Sensitivity Analysis"). We concurrently applied a slope criterion (ice <40 • ) and an elevation criterion (ice >2,800 m a.s.l. for the Hunza and ice >4,200 m a.s.l. for the Trishuli) to further constrain the clean ice class, i.e., to exclude steep slopes (rock) and lower elevations (illuminated moraines and bright surface water) from the snow/ice class. The slope thresholds are consistent with previous studies in the Hunza (Hewitt, 2011;Khan et al., 2015), where slopes >35 • represent internal rocks in the accumulation area of glaciers. Based on visual interpretation, here we chose a maximum value of 40 • to constrain the ice class. The outcomes of this step were clean ice and snow on ice classes.
Debris covered ice was mapped by thresholding the band ratio (ETM+ 4/5, OLI 5/6), a slope derived from the DEM, a thermal band (ETM+ band 6, OLI band 11) and an elevation layer. Each threshold was adjusted for each image on the basis of the visual interpretation of color composites; these are reported in Table 1. The thresholds for slope and thermal band were based on previous research (Racoviteanu and Williams, 2012). Minimum and maximum values of thermal bands were adjusted based on visual interpretation; these are generally lower in the fall or winter months, when debris covered ice surfaces are cooler due to less solar heating and lower temperatures, but they vary due to different illumination conditions. Only the exposed (snow-free) part of the debris covered ice was mapped in this step.
Snow on land was estimated for each image using the minimum snow and ice area, i.e., the "glacier mask" described in this section. For any given time step, when the extent of snow and ice was larger than the initial glacier mask area, the additional pixels were classified as seasonal snow and output as snow on land. For the purposes of parameterizing the melt model used in Armstrong et al. (2019), any snow present on the surface of debris covered tongues was also mapped and referred to as snow or land. We performed post-classification corrections to adjust misclassified areas, which included: (a) some glacier edges misclassified as snow as a result of the mixed spectra from glacier ice and rock; (b) shadows in the accumulation areas of glaciers which have a similar spectral signature to ice and were misclassified as exposed glacier ice; and (c) bare illuminated terrain misclassified as debris covered ice. Manual corrections were applied only to the subset study areas for purposes of conducting a sensitivity analysis to manual corrections.
Misclassified areas from (a) and (c) were adjusted manually, with only minimal processing of debris cover, which was not the focus here. Shadows misclassified as ice (b) were corrected automatically by creating a negative buffer (−2,000 m) inside the full glacier ice and snow mask, i.e., "shrinking" the mask. The false snowlines located within the shrunken mask were "erased" automatically with an overlay vector operation for all dates, thus reducing the time needed for manual edits. We refer to the resulting versions over subset areas as "corrected" and the raw, full extent versions as "uncorrected." Automated Snowline Extraction: The "Buffer" Method Snowlines were extracted in an automated way using a buffer applied to the snow and ice areas obtained from the surface partition at each time step (Figures 3A-H). Buffers have been used to extract snowline elevations in the Central Himalaya (Garg et al., 2017) or to estimate uncertainties in mapped glacier area in the Cascades (Granshaw and Fountain, 2017). While Garg et al. (2017) extracted snowlines only over the centerline of glaciers, here we estimated snowlines over the full extent of glaciers. Furthermore, while Garg et al. (2017) fixed the buffer size to 30 m around manually delineated snow areas, in this study we adapted the buffer size to the two DEMs tested (50-m buffer size for the 90-m CGIAR, and a 15-m buffer for the 30 m SRTMGL1). We performed a sensitivity analysis to the DEM and the buffer size used (section "Sensitivity Analysis"). For each time step, we intersected snow buffers to obtain a snowline interval, and then extracted elevations on a pixel-by-pixel basis from the DEM to obtain snowline elevations (Figure 3). Pixel-by-pixel SLA values are sensitive to outliers, so we estimated regional SLAs as the median elevation within the buffer interval over the full or subset image. Using the median rather than the mean accounts for non-normally distributed data in some of the months. The maximum SLA in a given year was considered the ELA for that year (Figure 3). SLAs were output automatically in tables, along with the ID of the scene and the parameters used for each run.

Uncertainty Estimates
Snow and ice mapping derived from remote sensing are subject to uncertainties issued from image classification techniques, mixed pixels and internal rocks, as discussed in previous studies (Racoviteanu et al., 2009;Paul et al., 2013). In this study, we used a standardized band ratio method to distinguish snow and ice from surrounding terrain, with carefully chosen thresholds. On the basis of previous studies (Paul et al., 2013), the accuracy of the glacier outlines derived from remote sensing using the automated technique is estimated as ± 1 pixel size (30 m). Multiple digitizing experiments have also found the accuracy of the automated glacier outlines to be within the variability of those obtained by manual digitization (Paul et al., 2013). Since here we only mapped the entire ice masses, and intersected them to get snowlines, we are not reporting the glacier outline accuracy directly, but we rather focus on the snowline accuracy.
Sources of uncertainty in the SLA estimates come from: (1) the accuracy of the DEM used for the surface partition (90 m CGIAR and 30 m SRTMGL1); (2) the size of the buffer used to extract the snowlines; (3) uncertainty in snow and ice area estimates; and (4) overall uncertainties in the automated method used to partition the surfaces. These were defined and estimated as follows: • ε dem is the vertical error of the SRTM DEM (RMSE z ), defined as 47.2 m for the 90 m CGIAR DEM and ± 23.5 m for the 30 m SRTMGL1 DEM based on Mukul et al. (2017). • ε buffer represents the buffer size used for the snowline extraction, defined as ±15 m for the 30 m SRTMGL1 and ±50 m for CGIAR (see section "Automated Snowline Extraction: The 'Buffer' Method"). • ε outlines is considered to be 1/2 of the pixel size of the Landsat satellite imagery used, i.e., ±15 m. We consider this error to be embedded in the buffer size of ±15 and ±50 m. • ε edit is related to uncertainties in the algorithm itself, which varies with each image and is related to the choice of band ratios thresholds and elevation criteria chosen for each time step. We quantified this error as the difference in SLAs resulting from the uncorrected (raw) and corrected (manually edited) versions of the surface partition for each time step (see section "Surface Partitioning: Clean Glacier Ice, Snow on Ice, Snow on Land and Debris Cover"). Assuming that the individual sources of error are uncorrelated, we estimated SLA accuracy as root mean square error (RMSE z ) for the uncorrected and corrected versions. For the snowlines issued from the uncorrected versions (full and subset extents), the total error (ε uncorr ) is: For the snowlines issued from manually corrected surface partitions, the total error (ε corr ) is:

Snowline Validation
We validated the snowlines in one of our study areas (Trishuli), where PlanetScope high-resolution imagery was available for the same date as the Landsat OLI scene (October 25, 2016) (Planet_Team, 2017). There were no Planet scenes available for the Hunza area, so we could not perform a similar validation there. The PlanetScope sensor acquires data at visible to near-infrared wavelengths (four bands), with a swath width of 24.6 × 16.4 km and 3 m ground resolution. PlanetScope provides image stripes of orthorectified, Top of Atmosphere radiance (at sensor), as analysis-ready data. Image stripes were combined in a single mosaic and bands 4, 3, and 2 were used to create false color composites. We manually digitized a total of 32 visible snowlines as vectors on the PlanetScope image using these color composites for individual glaciers from RGI. We rasterized the lines and extracted elevations along each digitized line based on the SRTMGL1 DEM (30 m), and averaged these to obtain PlanetScope SLAs for each glacier. Given the DEM cell size, a 15 m buffer is automatically embedded in the calculation of PlanetScope SLAs. We compared the OLI-based SLAs with those derived from PlanetScope using parametric statistical tests and basic statistics.

Monthly Surface Partition: Hunza (2013) and Trishuli (2016)
Here we present the monthly fluctuations of the various surfaces for Hunza (May-December 2013) (Figure 4) and Trishuli (June-December 2016) (Figure 5). The areas of "NoData" represent shadows and clouds; while most images are cloud-free over glaciers, deep shadows get progressively larger in the winter (November through December) in both areas due to low sun angles and the effect of rugged terrain. Clouds are more prevalent in the Trishuli, especially during the months of October to December. In Table 2 and Figure 6 we present summaries of total area for surfaces on the glacier and off-glacier in the subset areas (Hunza and Trishuli). In both areas, as snow decreases, exposed glacier ice increases. For example, for the Hunza subset in 2013, the exposed glacier ice increased from a minimum of 4.5% of the glacierized area on May 18 to a maximum of 24.1% of the glacierized area on October 9, and then decreased again to 6.2% in November after snowfall ( Table 2). These seasonal fluctuation patterns are visible in Figure 4, which shows glacier ice progressively being exposed as the season progresses until October 9, and then snow on land increasing to a maximum on November 26. The same patterns are visible in the Trishuli (Figure 5). In the winter months, almost no exposed glacier ice is visible over the entire extent (Figures 4, 5).
In the Trishuli, the exposed glacier ice increased from a minimum of 4.7% of the glacierized area on June 3 to a maximum of 26.8% on December 12, 2016, which was the maximum for this particular year (Table 2 and Figure 6). Minimum exposed glacier ice occurred roughly around the same time of the year (May-June) for both areas, but the maximum occurred later in the season in the Trishuli (December 2016) compared to the Hunza (October 2013) for the years studied ( Table 2). In the Trishuli, the total glacierized area exhibited an apparent decline by about 200 km 2 in exposed glacier ice from spring to winter in 2013; however, this was due to a large number of NoData values due to clouds or shadows obscuring the glacier surface during winter season in this area, and was not a "true" glacier ice loss. Snow on ice in the Hunza decreased from a maximum of 87.7% of the glacierized area on May 18 to a minimum of 63.0% of the glacier area on October 9, when maximum glacier ice was exposed and then increased again to 87.4% on November 26 ( Table 2). In the Trishuli, snow on ice decreased from 81.9% of the glacierized area on June 3 to 55.9% on December 12 (Table 2 and Figure 6) and then increased slightly to 57.6%. The debris covered ice area showed little variability, since debris cover is not expected to change significantly throughout the year ( Table 2). Debris covered surfaces were being progressively exposed as snow at the debris surface melted (Figure 6). Since the snow on land class includes snow on debris as well, the debris covered ice area fluctuation was not a true surface area change. In the Hunza, the debris covered ice area ranged from a minimum of 7.8% of the study area in the spring, when part of the debris covered glacier tongues were covered by snow, to a maximum of 12.9% on October 9 after all snow had melted. Similarly, for the Trishuli, debris cover extent ranged from a minimum of 13.4% on June 3 to a maximum of 17.3% on December 12 ( Table 2).
Outside the glacierized area, snow on land fluctuated from 27.7% of the total subset area in the Hunza on May 18 to a minimum of 0.0% on October 9 in the Hunza and from 7.3% of the total subset area in Trishuli in June to 0.01% on November 26 and December 12. Snow outside glaciers diminished as the melt season progressed, as expected (Table 2 and Figure 6).

Monthly Snowline Fluctuation in the Hunza (2013) and Trishuli (2016)
In the subset area of Hunza (Shimshal valley), the automated surface partition with manual corrections yielded the lowest SLA of 4,727 ± 69 m a.s.l. on May 18 and the highest SLA of 5,171 ± 69 m a.s.l. on October 9 (Table 3 and Figure 7A). After this date, SLA decreased in late fall/winter, which is consistent with the increase in seasonal snow in November and December ( Table 2). The date of highest SLA in the Hunza coincides with the maximum exposed glacier ice on October 9, and therefore we consider this date to be the end of the ablation season in   (Khan et al., 2015) had used July-August as the window of time to estimate maximum SLA for the Hunza, which is earlier than our findings. For the Trishuli, manual corrections of the surface partition yielded the lowest SLA of 5,261 ± 27 m a.s.l. on June 3, and the highest SLAs of 5,495 ± 27 m a.s.l. on December 12, and 5,499 ± 27 m a.s.l. on December 28 (Table 3 and Figure 7B). Since on the December 28 image we detected some seasonal snow (0.8%) (Table 2 and Figure 6), we consider December 12 to be the end of the ablation season for 2016.
The month marking the end of the ablation season was different in the two regions. In the Hunza, the highest SLA in 2013 occurred in October, versus in December in Trishuli in 2016. We cannot elaborate on the potential causes but we speculate that this later-than-usual highest SLA in the Trishuli could be linked with increased summer air temperature trends as noted in the Nepal Himalaya (Fujita and Nuimura, 2011), and potentially drier conditions in 2016. The temporal frequency of images from Landsat does not allow us to determine the precise day of the highest SLA, only the month.   (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) Annual ELAs extracted from images from the end-of-ablation season date established above (October-November) (Figure 8) show little variability from year to year from 2000 to 2016 in both areas (on average within ± 100 roughly with a few exceptions), indicating somewhat stable glacier conditions. In the Hunza, annual ELAs ranged from 5,142 ± 100 m a.s.l. in 2009 to 5,228 ± 97 m a.s.l. in 2013, with an average of 5,176 ± 100 m a.s.l. over the entire period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). This is in general agreement with other ELA estimates in the Karakoram for the same time period, for example Minora et al. (2016), who reported an average ELA of ∼ 5,200 to 5,300 m a.s.l. for the 2001-2010 period in Central Karakoram based on Landsat images; they assumed Julian day 273 or nearest date as the reference for the end of summer. Since they reported the ELA as a single elevation over the entire period, we cannot compare our results with theirs at annual time steps. Our ELA estimates are lower than those of Kääb et al. (2012) (ELA: 5,540 m a.s.l.) for the Karakoram region based on Landsat scenes.

Annual ELAs in the Hunza and Trishuli
In the Trishuli, annual ELAs varied from a minimum of 5,395 ± 48 m a.s.l. in 2000 to a maximum of 5,567 ± 48 m a.s.l. in 2016, with an average of 5,444 ± 48 m a.s.l. over the entire period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). This is in close agreement with regional ELAs of ∼5,400 m a.s.l. used in recent hydrologic studies in the same area (Ragettli et al., 2015). In a previous study, we estimated mean ELA values of 5,468 m a.s.l. in 2003 based on manual digitization on ASTER imagery (Racoviteanu et al., 2013). In the present study, we obtained an ELA of 5,432 m a.s.l. higher ELAs in 2015 and 2016 may be linked with local climatic conditions. Huss and Hock (2015) used a single ELA value of 5,215 m a.s.l. for southwest Asia for their model and projected a rise in the regionally average glacier ELAs by 100 to 500 m from 2010 to 2100. Due to the time scale of our analysis (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and the lack of a clear temporal upward trend in ELA, we cannot make any statements with regards to the suggested trends in the cited study. On the contrary, for the Hunza, our results imply relatively stable conditions, in agreement with cooling trends in summer temperature noted in the last decade (Forsythe et al., 2017).

Comparison With SLAs From High Resolution Imagery
In the Trishuli, the comparison between semi-automated (Landsat 8 OLI) vs. manual (PlanetScope) SLAs (Figure 9) yielded an average difference of 41 m, with a standard deviation of 64 m and a RMSEz of 137 m. SLA values obtained from Landsat and Planet had equal variances based on an F-test at 95% confidence interval. The mean SLA values from PlanetScope (5,301 m a.s.l.) vs. Landsat OLI (5,258 m a.s.l.) were not statistically different based on the t-test for two sample assuming equal variances. On a glacier-by-glacier basis, differences in Planet vs. Landsat SLAs ranged from −53 m (overestimated) to +240 m (under-estimated). In general, the semi-automated algorithm using Landsat 8 OLI underestimated most SLAs compared to PlanetScope, but results show good agreement (R 2 = 0.95) (Figure 9). The underestimation in OLI SLAs is due to the presence of brighter pixels or old snow in the lower ablation area of glaciers. Due to the similarly spectral signal of these ice pixels to snow, the algorithm detects an ice-snow boundary and classifies these pixels as part of the snowlines. Since these patches are situated at lower altitudes, these misclassified pixels introduce negative biases in the SLA estimates, and they need to be filtered out manually.

Spatial Patterns in SLA/ELA and Possible Links With Climate
Spatial patterns in SLAs for the month of October in the Hunza (2013) and in the Trishuli (2016) show a variability in SLA for both areas (Figure 10). In the Hunza (Figure 10A), SLAs decrease by 0.5 m vertical per 1 km eastward and increase by 8 m vertical per 1 km northward. The gradient is oriented slightly in the southeast -northwest direction (176.4 • ), with lower SLA values in the southeast and higher in the northwest. Our results agree with Mukhopadhyay and Khan (2016), who pointed out a strong east-west gradient in precipitation patterns in this area, causing ELA to vary across the Karakoram from about 4,840 to 6,200 m a.s.l. In the Trishuli subset area (Figure 10B), SLAs increase by 11 m vertical per 1 km eastward and 13 m vertical per 1 km change northward, and the trend is northeast (49.6 • ). SLA values in the Trishuli exhibit a larger spatial gradient in SLA than in the Hunza, as they increase toward the drier Tibetan Plateau, on the northeast side of our study region. In the Trishuli, SLA values ranged from 4,414 m a.s.l. on the south side of the divide to a maximum of 6,581 m a.s.l. the north side of the divide (Tibetan Plateau), with an average of 5,425 m a.s.l. This is in agreement with our estimates of mean ELA values of 5,468 m a.s.l. in this area from the previous study based on manual digitization on ASTER imagery (Racoviteanu et al., 2013), as well as other studies mentioned previously Kayastha and Harrison, 2008;Ragettli et al., 2015). The high SLA values of up to 6,581 m a.s.l. in the northern part of the Trishuli area may be due to the drier climate compared to the southern slopes, caused by a decrease in moisture content in the northeast from orographic forcing of monsoon air masses over the Himalaya. The strong gradients in SLA values towards the northeast (Tibetan Plateau) indicate that SLA values calculated over the entire Landsat extent may not be representative of the entire region in this area. To test this, we calculated SLAs separately on the northern slopes of the Himalaya (Tibetan Plateau) and southern slopes (Trishuli basin) of the subset area. Using only the SLA pixels on the southern slopes, we obtained an average SLA of 5,222 m a.s.l., which is about 4% lower than the overall estimates over the entire scene including the Tibetan Plateau (5,425 m a.s.l). On the northern side of the divide, more glacier ice is exposed in late season than on the southern side. The larger expanse of exposed glacier ice on the northern slope, and the higher SLAs appear to bias the regional trends when analyzing the entire Landsat scene in this area. When SLA values are calculated using an entire scene which spans distinct climatic regions, regional SLA values are overestimated by ∼200 m, so a basinby-basin approach would be more appropriate. This observation has important implications for regional ELA estimates intended for mass balance applications, as glacier trends, including ELA trends, do not often follow hydrologic boundaries.

Sensitivity Analysis
While band ratio thresholds are fairly consistent throughout glacierized regions worldwide (e.g., Paul and Andreassen, 2009;  Racoviteanu et al., 2009), it is generally advisable that the thresholds be checked in each region. We found that the SLAs were not very sensitive to the threshold used for band ratios to distinguish snow and ice from surrounding terrain. When using thresholds 1.2 and 2, the differences in resulting SLAs ranged from −73 to +26 m, on average −0.2 to −0.04% compared to SLA using the chosen threshold (1.5) ( Table 4).
With respect to OLI band 5 used to distinguish snow from ice on the glacier surface, we tested two thresholds based on statistics extracted from the ROIs: mean + 1 SD and mean + 3 SD (Figure 11). Since the ROIs digitized on glacier tongues and used to extract band ratio statistic were chosen lower in the (darker) ablation areas of glaciers, we estimated that the adequate DN threshold would be above the mean. The sensitivity analysis for Band 5 (Table 4) showed that the SLAs were fairly sensitive to the choice of threshold, with SLA differences ranging from −244 m (underestimate) to +411 m (overestimated). Band 5 DN values calculated over the ice ROIs show seasonal variability (Figure 11). We found that the optimized DN chosen corresponds closely to mean + 1 SD in the spring and late fall months, and to the mean + 2 SD during the summer season (June to October). The average SLA differences were 1-3% lower than those obtained using the optimized threshold. SLA differences were larger in the spring (May) and late fall (November), when snow covered more of the glacier tongues where ROIs were chosen.
Season snowline altitude estimates were sensitive to the topographic slope threshold during the months of May and November (Table 4). For these months, SLA differences ranged from −271 to +214 m, however these were still within our SLA uncertainty estimates for those months ( Table 3). During the months at the end of ablation season (October and November), SLA differences ranged from −26 to +63 m. Overall, SLA differences due to the various slope thresholds only amount to −1 to 2% compared to those obtained using the optimized value.
To quantify the sensitivity of the snowlines to the DEM used for the surface partition step (90 m CGIAR vs. 30 m SRTMGL1), we compared SLAs obtained from the uncorrected partitioning for the full extent of the Landsat imagery in the Hunza based on the two DEMs, using the same buffer size (50 m). In this area, the uncorrected SLAs over the full extent of the Landsat scene using the 90 m CGIAR ranged from 5,239 ± 575 m a.s.l. (May) to 5,245 ± 86 m a.s.l (October). Using the 30 m SRTM DEM, resulting SLAs ranged from 4,990 ± 575 m a.s.l. (May) to 5,189 ± 106 m a.s.l. for (October) ( Table 4). SLAs obtained using the CGIAR were lower than the values obtained using the SRTMGL1. The differences (CGIAR-SRTMGL1) decreased from the spring and early summer months (May and June) (167 m on average) to the fall months (October-November) (35 m on average) ( Table 4). This difference, calculated as percentage of the CGIAR-based SLA, amounts to a maximum of 4.8% in the spring months and a minimum of 0.6% in October. The two sets of SLA estimates were statistically different based on two sample t-tests at 95% confidence interval (p < 0.05), suggesting that the choice of the DEM significantly impacts the snowline estimates. Regardless of the DEM, we consider the SLA results from the spring months to be less reliable due to poorer image quality, more clouds and more snow. When different DEMs were used only in the last step of the analysis (i.e., extracting elevations along the SLAs) from the same version for partition, differences in resulting SLAs were negligible (∼2 m on average). The size of the buffer (50 m vs. 15 m) used to extract the snowlines after surface partition had little impact on the resulting SLAs based on sensitivity tests for  (brighter), and to mean + 2 SD during July-November when more glacier ice is exposed (darker). The black line shows the "optimized" DN, which was chosen based on the ROI statistics and visual checking using color composites.
the Hunza full extent ( Table 4). Differences in SLA calculated using the 50 and 15 m buffers, respectively, ranged from −1 m in October 9 to 10 m on November 26, with an average of +5 m, or +0.1% difference between the two versions ( Table 4). We conclude that SLAs are not sensitive to the size of the buffer used to extract the snowlines.
To quantify the impact of manual corrections on the SLAs, we compared the SLA resulting from the uncorrected and corrected surface partition in the Hunza and Trishuli. In the Hunza subset area, SLAs from "uncorrected" surface partition ranged from 5,301 ± 578 m a.s.l. in May 18 to 5,261 ± 254 m a.s.l. on November 26, with a month-to-month variability. After manual corrections, SLAs in the subset area ranged from 4,727 m a.s.l. ± 69 (May) to 5,171 ± 69 m a.s.l. (October), with the biggest difference in the spring month (574 m, or 11%) ( Table 3). SLAs resulting from the "corrected" surface partitions were generally lower than the "uncorrected" version (a mean difference of ∼208 m or about 4% of the uncorrected version). The difference in SLA values from uncorrected vs. corrected partition was generally smaller in the summer and fall months (July to October) due to better image quality and less seasonal snow ( Figure 7A). For these months, the average difference in median SLAs between these two versions was on average 68 m, about 75% less than the average difference in the spring and winter months (291 m). In the Trishuli subset area, SLAs from "uncorrected" partition in the Trishuli subset ranged from 5,324 ± 27 m a.s.l. in June to a maximum of 5,534 ± 48 m a.s.l. in December (Table 3). After corrections, SLAs in the subset Trishuli ranged from 5,261 m a.s.l. in June to 5,495 m a.s.l. in December, with the biggest difference between the uncorrected and the corrected version in the summer month (June) (63 m difference, or 1.2%). After the corrections, similarly to the Hunza, SLAs in Trishuli were on average lower than before the corrections (a mean difference of 39 m, or 0.7% compared to the uncorrected values) ( Figure 7B). Most likely, the lower SLAs after corrections are due to removal of shadows in accumulation area, initially misclassified as ice. To assess the effectiveness of the manual corrections, we used a paired t-test on the precorrected and post-corrected versions. The test indicated that SLAs computed after the manual corrections were statistically different than those issued from the uncorrected version at 95% confidence interval (p < 0.05) in both areas, implying that the corrections significantly changed the mean values. Similarly to the Hunza, in the Trishuli subset area, after manual corrections of the surface partition, SLA values were statistically different than the ones obtained from the uncorrected partition (Table 3), based on two-sample t-test at 95% confidence interval, p < 0.05.

Continuity From Landsat 7 ETM+ to Landsat 8 OLI and Impact on Estimated ELA
The temporal trends of ELA may be subject to uncertainties associated with the Landsat sensors used: ETM+ (2000ETM+ ( -2012 and OLI (2013OLI ( -2016, and with respect to the transition from ETM+ to OLI in 2013. In this study, we used Landsat L1T products (pre-collection) and L1TP T1 (collection 1), which are registered, radiometrically calibrated and orthorectified using ground control points GCPs and a DEM (USGS, 2015). These products have stated 50 m global geolocation accuracy, but the actual accuracy was found to be better than the projected. The alignment of the scenes with respect to one another was assessed visually, and further co-registration was not considered necessary. A source of error in the Landsat 7 ETM+ is that scenes acquired after 2003 had the scan line corrector (SLC) turned off due to failure and have data gaps, visible as "stripes" on the images. These were more pronounced toward the edges of the images in our study area, and they were treated as "NoData" in our study. The gaps in the data decrease the total sample size for Landsat 7, when present. The Landsat 8 data used after 2013 has several advantages with respect to ETM+, notably narrower spectral bands, improved calibration and signal-tonoise ratio (SNR), better geometry and 12-bit radiometric resolution compared to 8-bit for Landsat ETM+ (Irons et al., 2012). Landsat OLI data are not saturated, while Landsat ETM+ has pixels that are saturated in bands 1, 2, or 3 (the visible bands), which may appear darker than they really are. In addition, the USGS uses cubic-convolution resampling to smooth the data, not nearest neighbor resampling. As a result, saturated pixels can be mixed with nearby unsaturated ones in this resampling, making it difficult to identify and exclude those from our analysis. Overall, we consider the results from OLI to be more reliable based on these characteristics, with an accuracy of 50 m based on comparison with PlanetScope (section "Uncertainty Estimates"). As no clear ETM+ vs. OLI trend was detected, we did not attempt any inter-calibration, but future studies might consider this issue.

Limitations and Further Work
Delineating the snow/firn-ice boundary using remote sensing in mountainous terrain using optical remote sensing is challenging due to shadows, since snow in shadow appears spectrally similar to clean ice or dirty ice in the ablation area of glaciers. Limitations of the current study include: • Our method does not allow making the distinction between snowline and firnline, but we acknowledge that under warming conditions, the snowline may retreat beyond the firnline. Applying this method on daily high-resolution imagery may facilitate this distinction. • Distinguishing snow on land from snow on ice remains a challenge, since snow has the same spectral signature on and off the glacier; here we relied on a glacier mask at each step. • Band ratio algorithms, while fairly robust, are prone to misclassification errors over frozen lakes and deep/cast shadows over snow and ice (Racoviteanu et al., 2009). Using topographically corrected surface reflectance to reduce variability in illumination instead of top of atmosphere data might mitigate post-classification correction efforts. • Some snowlines are obscured by shadow and they often do not appear as continuous lines, which may slightly bias our regional estimates especially in areas which span different climatic areas. • The thresholds used for surface partition rely on a-priory knowledge to some extent. Using a single threshold is challenging in rough terrain due to differences in topographic illumination or image saturation as well as climatic differences across regions. However, in this study, we used a consistent threshold across all images as much as possible. • Further steps to improve the efficacy of our algorithm include more sophisticated cloud mapping techniques and automated shadow detection using sub-pixel methods (Sirguey et al., 2009). • Post-classification manual corrections of the surface partition were most needed for spring and winter image(s), with minimal corrections and greater accuracy for images acquired later in the season, which were more contrasted and generally cloud-free. We developed tools to automate the post-classification corrections in areas of shadow over snow in the accumulation areas to reduce the time needed for manual edits. • SLA/ELA estimates are limited to monthly temporal resolutions due to the revisit time of the Landsat imagery (16 days), but the method can be easily applied to Sentinel-2 data (5 days). With further testing and refinement and perhaps calibration with field measurements of SLA/ELA, this method can be modified for daily satellite data such as PlanetScope to derive sub-monthly SLAs.

CONCLUSION
Manual digitization of snow and ice on glaciers and subsequent extraction of SLAs is generally a time-consuming process and is difficult to apply over large areas, especially when time series of the snowlines are needed. Here we developed an automated method to separate snow-ice-debris surfaces in two areas of the Himalayas at multi-temporal scales based on Landsat ETM+ and OLI band ratios and topographic criteria. We extracted snowline elevations pixel-by-pixel and estimated ELAs in two areas of HMA: Hunza and Trishuli at monthly time scales (2013 and 2016, respectively) and annual time scales (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). SLA estimates were significantly sensitive to the manual corrections applied to the snow and ice partition results, and fairly sensitive to the topographic slope, the DEM and the band ratio thresholds, particularly during the spring and winter months due to more snow. Snowlines were less sensitive to the size of the buffer used to extract the snowlines. Using this method, we obtained a maximum SLA (∼ELA) of 5,171 ± 27 m a.s.l. in October of 2013 for the Hunza and 5,495 ± 27 m a.s.l. for the month of December of 2016 for the Trishuli., after manual corrections. Over the period studied (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), end-of-the-ablation season annual ELAs fluctuated from 4,917 to 5,336 m a.s.l. for the Hunza, with a 16-year average of 5,177 ± 108 m a.s.l., and 5,395 to 5,565 m a.s.l. for the Trishuli, with an average of 5,444 ± 63 m a.s.l. SLA trends obtained over a smaller subset of the Landsat scenes were representative of the full extent of the image, with an average difference of 100 m. We consider that regional SLA/ELA values obtained using this method are adequate for regional applications such as melt models, when image quality, time of the year, uncertainties due to the DEM used, and band ratio thresholds are considered and assessed. Caution is needed when extracting a "regional" annual ELA value when an image spans various climatic regions, as SLAs may be biased. The time series of snow, ice and debris in two regions provides a valuable training dataset, which may be used in future work to classify images using more sophisticated algorithms such as machine learning classifiers, such as random forests.

DATA AVAILABILITY
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
AR ordered and processed the satellite imagery, developed the python scripts and remote sensing methods to partition the different surfaces and to extract the snowlines, and was also responsible for writing the main body of the manuscript. KR contributed to the development of the algorithm though versioning and feedback, figure generation, and extensive manuscript editing. RA contributed to the theoretical glaciological concepts and manuscript editing. All authors contributed to this work. FUNDING AR, KR, and RA were primarily supported by the Contribution to High Asian Runoff from Ice and Snow (CHARIS) project, funded by the United States Agency for International Development (USAID, Cooperative Agreement No. AID-0AA-A-11-00045). Additional support for AR was from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No 663830. KR was also supported by the NASA Applied Science Grant 80NSSC18K0427 and the NASA HiMAT project.