Mapping seasonal glacier melt across the Hindu Kush Himalaya with time series SAR

Current observational data on Hindu Kush Himalayas (HKH) glaciers are sparse and characterizations of seasonal melt dynamics are limited. Time series synthetic aperture radar (SAR) imagery enables detection of reach-scale glacier melt characteristics across continents and these records may be integral to the development and assessment of surface energy balance models of glacier ablation. We analyze C-band Sentinel-1 A/B SAR time series data, comprised of 32,741 Sentinel-1 A/B SAR images, determine the duration of seasonal glacier melting for 76,831 glaciers (65,108 km), defined using optical 15 observations, in the HKH across the calendar years 2017-2019. Signals of melt onset and duration are recorded at 90m spatial resolution and 12-day temporal repeat across 97% (62,907 km) of the HKH cryosphere. Melt signals persist for more than 40% of the year at elevations below 4,000 m a.s.l. and in excess of 15% of each year at elevations exceeding 7,000 m a.s.l. Melt retrievals resemble characteristics of glacio-climatic subregions of the HKH: melt onsets sooner and occurs for a longer portion of each year in the Central and Eastern Himalaya compared to the Western Himalaya and Karakoram regions. 20 Retrievals of seasonal melting span all elevation ranges of significant glacier area in the HKH region, extending greater than 1km above the maximum elevation of an interpolated 0oC summer isotherm. Furthermore, percolation zones are apparent from meltwater retention indicated by signals of delayed refreeze. Time series SAR datasets are suitable to support operational monitoring of glacier surface melt and the development and assessment of surface energy balance models of melt-driven ablation across the global cryosphere. 25


Introduction
Global warming driven by the anthropogenic release of geologic carbon is causing mass wasting of alpine glaciers worldwide Zemp et al., 2006). The Hindu Kush Himalaya (HKH) region, known colloquially as the "Third Pole," has the most ice-covered area on Earth after the high-latitude polar regions (Yao et al., 2012). In contrast to large ice https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
sheets near the poles, these relatively small alpine glaciersperched at some of the highest elevations on Earthare among 30 the most sensitive indicators within the global cryosphere of changes in global climate (Anthwal et al., 2006). Just as the recession of these sensitive mountain ice caps contributes to over 25 percent of global sea level rise (Zemp et al., 2019), disturbances accompanying HKH glacier retreat pose additional and innumerable hazards to humans and natural ecosystems.
Glacier retreat threatens to disturb the dynamics of river systems delivering freshwater resources to nearly 2 billion people across South and Central Asia (Brown et al., 2007;Milner et al., 2017). Outburst floods resulting from glacier mass wasting 35 have killed at least 6,300 people in the Himalayas alone, and have caused extensive damage to property and livelihoods. These outbursts are expected to increase in frequency with continued glacier wasting (Carrivick and Tweed, 2016). Some organisms endemic to alpine aquatic ecosystems may become extinct as they lose biogeochemical regulation from glaciers upstream (Jacobsen et al., 2012). As global temperatures rise and perennial snow and ice cover decreases, societies are faced with difficult decisions around the costs and benefits of adapting to a changing climate within and around the HKH region (Brown 40 et al., 2007). Informed decision-making for successful climate change adaptation will require knowledge of the state of natural systems and how these systems are projected to change alongside future increases in population and global average temperature (Bogardi et al., 2012). Substantial uncertainties exist in the current understanding of projected disturbances associated with a changing climate, environment, and hydrologic regime across the greater Himalayas due to a lack of observations of hydrology and 45 meteorology at high elevations (Litt et al., 2019). Although the general trajectory of changes to the HKH cryosphere is understood (i.e. accelerated glacier mass loss on a decadal scale), a consensus in projecting changes to HKH hydrology is lacking largely because of missing in situ snow and ice monitoring data across these glaciated river basins (Fujita and Nuimura, 2011). However, construction and maintenance of an in situ monitoring station network is costly and labor-intensive because of the complexity of the high-mountain glaciated terrain. Remote sensing observations are commonly used in lieu of weather 50 station measurements. In particular, satellite imaging radar retrieval of alpine glacier melt characteristics has long been proposed as a source of data for hydrologic and glaciologic research (Shi et al., 1994). Advanced climate change projections of glacier wasting require snow property dynamics that describe the retention, refreezing and drainage of liquid water within glacier snow and firn (Pritchard et al., 2020) and therefore observational records of this type are of great value. Surface melting drives accumulation-zone snow-properties, such as percolation and densification, that can increase the seasonal amount of 55 melting (Alexander et al., 2019) and has been linked to increased englacial temperatures resulting in faster ice motion (Miles et al., 2018).
Recent revelations indicate that shortwave radiation drives melting at high elevations, above 5500 m a.s.l., like on Mount Everest , where temperatures never exceed -10ºC (Matthews et al., 2019;Matthews et al., 2020). These finding assert that temperature-indexed melt models are surely underestimating ablation at these elevations through the assumption that air 60 temperatures greater than 0ºC alone drive snow and glacier melt. Further, studies of glacier wasting in High Mountain Asia have shown variability in patterns and magnitude of glacier wasting across sub-regions of the HKH that would be difficult to https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. capture in numerical models using degree-day assumptions (Brun et al., 2017). An observationally based dataset on characteristics of the glacier surface energy balance is necessary to attempt to capture seasonal and regional variability in glacier wasting across the HKH during melt-freeze cycles. 65

Snowmelt Detection and Radar Imaging
This study builds on extensive research into microwave scattering from dry and wet snow and techniques for snowmelt retrieval from imaging radar sensors to present an operational monitor for spatially-resolved glacier surface melt characteristics using synthetic aperture radar (SAR) time series and up-to-date glacier outlines derived from satellite optical imagery across 70 the HKH. Microwave remote sensing is used to reliably monitor melt patterns across glaciers and ice sheets (Abdalati and Steffen, 2001;Bahr et al., 1997;Jezek et al., 1994b;Kimball et al., 2004;Mote et al., 2017;Rawlins et al., 2010). Because imaging radar is independent of solar illumination and largely unaffected by cloud cover and atmospheric convection, the fidelity of radar observations is limited only by the frequency of observational opportunities and the characteristics of the imaging sensor. At C-band frequencies, frozen glacier percolation areas are recognized as one of the brightest radar targets on 75 Earth and glacier surfaces are unambiguous targets for determination of surface melt/freeze characteristics (Jezek et al., 1994a;Rott and Mätzler, 1987). Detection of seasonal melt on ice surfaces at C-Band frequencies (4 -8 GHz) depends on a strong radiometric response at melt onset (MO), when liquid water content introduced to an otherwise frozen snow or firn matrix causes a drastic decrease in the radar backscatter from the medium (Hallikainen et al., 1986). Deep, frozen snow and firn has a high scattering albedo across microwave frequencies (Matzler, 1998), resulting in high radar backscatter intensity over 80 glaciated regions during the frozen months (Winsvold et al., 2018a;Wiscombe and Warren, 1980). The introduction of liquid water in the snow or firn matrix at even hydrologically minimal amounts (<0.1% by volume) causes a pronounced increase in the medium's dielectric constant, increasing radar signal attenuation and diminishing volume scattering. This leads to a pronounced decrease in radar backscatter, usually by half power (-3 dB) or more (Kendra et al., 1998;Shi and Dozier, 1995).
In areas that are seasonally snow-free, like for areas of debris-cover or bare ice, melting conditions are dominated by surface 85 scattering that is significantly darker relative to winter conditions (Lievens et al., 2019). In part because of the relatively strong signal produced by melting, radar-based melt-detection records have been developed across regions of the global cryosphere for several decades using both real and synthetic aperture radar sensors (Bhattacharya et al., 2009;Bindschadler et al., 1987;Koskinen et al., 1997a). Subsequently, snowmelt detection algorithms have been developed using a host of radar sensors to monitor the onset and duration of snowmelt across glaciers and ice sheets (Abdalati and Steffen, 2001;Ashcraft and Long, 90 2007;Bahr et al., 1997;Jezek et al., 1994b;Kayastha et al., 2019;Koskinen et al., 1997b;Winebrenner et al., 1994). Prior applications of SAR mapping of seasonal surface melting over ice sheets and glaciers have been limited by a lack of repeat observations such as those now available from the Sentinel-1 SAR constellation (Lund et al., 2019). Launch of the Sentinel-1 https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
C-band SAR constellation ushered the first openly-available SAR time series dataset for change detection at the Earth surface.
Limited only by the frequency of observations (12-days per orbit direction) 95 Observations from time series SAR data are used to delineate zones of glacier facies and regions of glacier mass balance (Winsvold et al., 2018a). In glacier percolation zones, seasonally wet snow refreezes into ice lenses, pipes, and other percolation-related features that amplify both surface and volume scattering of C-band radar and result in the brightest SAR backscatter being captured during the frozen periods (Jezek et al., 1994a;Rau et al., 2000). Studies have shown that SAR backscatter intensity generally increases with elevation across glacier surfaces during frozen periods, from the glacier terminus, 100 through frozen meltwater percolation zones, and eventually attenuating in zones where dry snow accumulates (Winsvold et al., 2018a). In transitions between these zones there are pronounced backscatter contrasts rather than smooth, gradual transitions. At C-band frequencies, strong radar scattering within glacier percolation areas dominates the backscatter intensity signal during frozen periods (Jezek et al., 1994a). Importantly, the loss of volume scattering during surface melting in areas of meltwater percolation creates a pronounced, and unambiguous, radar signal in time-series observations. The sensitivity of SAR 105 sensors to the introduction of liquid water to an otherwise frozen snowpack or firn structure is a reliable target for the retrieval of glacier surface melt characteristics. In refreezing percolation zones, the upper layers of firn will become frozen first and freeze downward across layers, thus progressivly increasing backscatter and with decreasing total-column liquid water content (Huang et al., 2018). In this way, the timing of refreeze relative to the surface energy balance at the surface provides a direct and spatially resolved indicator of subsurface meltwater storage within the snow or firn and deliniates the percolation zones 110 over mountain glaciers.
This study applies methods used for decades in the canon of research on retrieving melt characteristics from glaciers and ice sheets but using SAR data acquired at a spatiotemporal resolution that can capture variability across a mountain glacier surface and with a frequency of data acquisition that can illustrate seasonal characteristics of melt onset and duration. In this paper we provide a brief history on the use of imaging radar for melt retrieval on glacier surfaces to motivate a simple threshold-115 based method for change detection of surface freeze/thaw state. It is likely that intense incident solar radiation is driving melt at elevations above the 0ºC summer isotherm (Matthews et al., 2019) across the entirety of the HKH, so that the sensitivity of SAR sensors to changes in the melt/freeze condition of glacier surfaces is a viable alternative to temperature elevation lapse rates (Litt et al., 2019) in order to develop and assess surface energy balance models of glacier ablation. Though coarse in temporal resolution relative to a typical meteorological dataset, retrieval of melt status using SAR time series produces 120 mappings with very high spatial resolution and a continuous record of melt timing and duration across glaciated regions. We present an application of this melt retrieval technique at the scale of the HKH with spatiotemporal fidelity adequate to capture seasonal variability in melt timing and duration across individual glacier surfaces.

Setting and Data
The HKH region ( Fig. 1) spans 13 million km 2 , including areas inhabited by 240 million people with nearly 2 billion 125 people relying on the delivery of water resources from catchments that originate within the region (Scott et al., 2019). Glaciers in the HKH region have been wasting rapidly alongside increases in global average temperature with the exception of the Karakoram mountain range at the intersection of South and Central Asia (Gardelle et al., 2012). It has been argued that the Karakoram region has maintained a state of relative equilibrium in its glacier mass balance because of its high latitude and inland setting which shields the region from monsoonal controls on melt that are active in other Himalayan catchments 130 (Kapnick et al., 2014). The wasting of HKH glaciers is thus a spatially and temporally heterogeneous phenomenon where distinct glacio-climatic regimes control ablation across (Bolch et al., 2012).

Figure 1. (Top) Hindu Kush Himalaya (HKH) region (red outline) and 2018 GAMDAM glacier inventory (GGI) (blue). The region
highlighted with the inset map shows GGI glacier areas across the Trishuli Basin in Nepal. GGI and HKH data are overlaid onto a 30m 135 Shuttle Radar Topography Mission (SRTM) DEM hill shade. (Bottom) Sentinel-1 ascending (black) and descending (blue) swath footprints acquired across the study region. Ascending orbit cycle number 56 is highlighted in red to illustrate the SAR image processing approach for time series analysis across distinct orbit cycles.

GAMDAM Glacier inventory (GGI)
The Glacier Area Mapping for Discharge from the Asian Mountains (GAMDAM) glacier inventory (GGI) is a recently 140 updated (July 2019) database on glacier outlines for the region of High Mountain Asia (Fig. 1). These outlines were originally delineated automatically using cloud and snow-free satellite optical imagery in an initial release of the database (Nuimura et https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. al., 2015). As a recent update to the database, each outline was individually inspected for quality control to correct discrepancies where automatic glacier delineation lost accuracy in terrain occluded areas, at debris covered portions of glaciers, and through obstruction under seasonal snowpack. The recently updated dataset covers a much larger area than the previous 145 version (Sakai, 2019). Although these data are the most current, they do not necessarily capture debris-covered portions of glaciers due to confusion with land in optical image classification schemes (Bolch et al., 2019). The 2018 GAMDAM database includes 76,831 distinct glacier outlines spanning a total area of 65,108 km 2 within the HKH (Nuimura et al., 2015).

Sentinel-1 Synthetic Aperture Radar
The Sentinel-1 A/B satellites were launched in April of 2014 and 2016, respectively, and collect C-band (5.405 GHz) 150 SAR data with a combined revisit interval of 6-days over of the majority of the mid-latitude terrestrial Earth. Each Sentinel-1 scene acquired in the interferometric wide-swath (IW) mode has a width of 250 km and a resolution of 5x20 meters in range and azimuth at the equator. This study considered images taken in the IW mode and in both co-polarized and cross-polarized state (VV, VH). Sentinel-1 data were accessed through a cloud-computing platform (discussed below) wherein SAR scenes were radiometrically terrain corrected to backscatter intensity values in decibels (dB) using the European Space Agency's 155 (ESA) Sentinel Application Platform (SNAP) toolbox and the Shuttle Radar Topography Mission (SRTM) 30m digital elevation model (DEM) (Farr, 2007) upon ingestion into the cloud environment. Data from both the ascending and descending orbit nodes were analyzed across the study region for a total consideration of 32,741 individual Sentinel-1 A/B IW scenes across 46 unique orbit cycles captured across the calendar years 2017-2019 (Table 1, Fig. 1b). By combining orbit directions, we utilize observations acquired at day and night. For the purpose of this study we do not attempt to resolve diurnal-scale melt-160 freeze processes and instead focus on retrieving annual characteristics of melt timing and duration.

180
The absolute difference between mean summer and winter VH backscatter from Sentinel-1. Note the larger cross-polarized (VH) response across a greater amount of area compared to co-polarized (VV) data.

Computing Infrastructure
A cloud-computing platform and application programming interface (API) with radiometrically terrain corrected Sentinel-1 A/B data was used to detect melt characteristics across the region (Gorelick et al., 2017). Radiometric terrain 185 correction of Sentinel-1 data was conducted upon ingestion to the cloud server using the ESA's method contained within the Sentinel Applications Platform (SNAP) processing toolbox (ESA). This method for radiometric terrain correction (RTC) https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
utilizes the 30m spatial resolution SRTM DEM. The SAR times series data and API functionality used to derive glacier melting characteristics is available from Google Earth Engine (GEE) (Gorelick et al., 2017).

Melt Classification
We use a threshold-based change detection algorithm applied to time series radar backscatter intensity to classify melt conditions (Ashcraft and Long, 2007). Melt detection is conducted across Sentinel-1 A/B ascending and descending orbit track time series separately and mosaicked into a final image based on a statistical score for seasonal melt magnitude after classification. To classify snowmelt, we conduct a pixel-based temporal classification by comparing each image at interval i 195 to a dry/frozen winter average backscatter value calculated from January to February for each study year. Snowmelt at each image acquisition interval ( ) was classified using Eq. (1): where the ground-range detected backscatter intensity at each image acquisition ( 0 ) within the times series must be less than the difference between the mean winter backscatter ( ̅ 0 ) and a fixed threshold (b). Threshold value (b) have been developed across numerous studies of melt detection with C-band scatterometer and SAR datasets using both ground-based observations and radar scattering model results of changes to backscatter magnitude at the onset of melt. For the purposes of this study we followed previous studies (Baghdadi et al., 1997;Bhattacharya et al., 2009;Engeset et al., 2002;Nagler and 205 Rott, 2000;Oza et al., 2011;Rott and Mätzler, 1987;Steiner and Tedesco, 2014;Trusel et al., 2012) and selected a b value equal to one half of the signal power (3 dB). Due to the lack of high-elevation meteorological data available across the study years, we were limited in our accuracy assessment to only one high elevation (>4,000m a.s.l.) meteorological station, located at the Yala glacier base camp (Fig. 3). Backscatter values averaged across the Yala glacier (4,950m a.s.l) acquired along the Sentinel-1 A/B descending orbit direction are plotted alongside mean daily air temperature recorded at the Yala glacier base 210 https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. camp automatic weather station (Shea, 2016). If we consider air temperature above 0ºC to control glacier surface melt at this location, classification accuracy for melt retrieval using Eq. (1) is 96% in the VH polarization.

Quantifying algorithm performance
Sentinel-1 SAR viewing geometry will vary as the local incidence angle increases with across-track range. At high 220 incidence angles (far range), the sensitivity to volume scatter is diminished and the melting signal is reduced. At C-band frequencies, these effects on volume scatter are strongest only at very high incidence angles (closer to grazing) (Nagler and Rott, 2000). We classified areas as valid for melt detection using a metric of statistical separability for seasonal backscatter intensity across frozen and melt periods, which we interpret as a measure of the strength of the seasonal melt signal Eq. (2): where the score for seasonal separability of backscatter intensity (z) was calculated across each SAR pixel's time series using the difference between the mean winter ̅ 0 (January-February) and summer ̅ 0 (July-August) season backscatter intensities, as compared to the standard deviation of backscatter across the winter months ( 0 ). In computing z, we which orbit direction (ascending or descending) to use for melt classification on a per-pixel basis after applying Eq. (1) across 235 each orbit cycle time series, so as to capture the maximum area of melt signals occurring across the complex terrain.
We consider the seasonal signal valid for time series melt detection if the seasonal backscatter intensities are separated by greater than two standard deviations ( > 2), representing better than 95% confidence in the presence of an annual melt signal. To choose which orbit direction to use for melt classification, we pick the pixel with the greatest z. We find that z generally increases with elevation across subregions of the HKH (Karakoram, western, central, and eastern Himalaya) up to 240 about 5,000m a.s.l. and that, across elevation ranges, the mean z is above the threshold, indicating detection of a seasonal melt signal across all ranges of glacier elevation spanning the four major glacio-climatic subregions of the HKH (Fig. 4). Areas of debris-cover may exhibit radar brightening with snow-free conditions above winter mean ( < 0). These areas occur towards lower elevations where seasonal snow, or firn, does not have significant contribution to the seasonal backscatter response and are not included in our melt classification approach. 245 Sentinel-1 A/B interferometric wide (IW) swath images have a range in viewing angle between 29.1-46.0° (ESA).
Glacier melt retrieval using SAR data commonly begins with a normalization of radar images by viewing angle on a scene by scene basis (Kimball et al., 2004;Man et al., 2014;Winebrenner et al., 1994;Winsvold et al., 2018b). We consider changes for each individual 10x10m pixel time series across distinct, repeating orbit tracks and directions. This approach holds the local incidence angle effectively constant in time. Glacier melt classification and z-score calculation are carried out across 250 images acquired along identical orbit tracks in distinct orbit directions (Fig. 1) and mosaicked into a final dataset for each study year using the greatest z-score observed across each orbit cycle path and in each orbit direction. We thus limit temporal resolution of melt retrievals to 12-days by choosing only observations from the orbit direction with the greater z-score. Time series analysis of SAR acquisitions on distinct orbit tracks eliminates the need to normalize each scene by incidence angle for the purposes of melt retrieval. This method reduces computational cost and eliminates artefacts that may originate from 255 overlapping orbit paths and differences in radar viewing angle. Areas where local incidence angle varies significantly because of complex topography should show little time series variability in backscatter change at the SAR pixel scale when viewed at a distinct and consistent orbit path and direction.
Mean seasonal melt magnitude averaged over 100m elevation bins over all three calendar years of data shows strong (z > 2) melt signals across glacio-climatic sub-regions and across all elevation ranges of significant glaciation (Fig. 4). The 260 occurrence of seasonal melt signals across all ranges of elevation in the HKH is both noteworthy and striking. Until very recently, glacier ablation in the HKH region has been treated in numerical model studies to be controlled by air temperatures greater than 0ºC (Litt et al., 2019). Our finding that SAR backscatter time series observe pronounced melt/freeze cycles across high elevation ranges of HKH glaciers supports other recent findings that intense incident solar radiation is driving glacier melt processes at elevations above the 0ºC monsoon season isotherm (Matthews et al., 2019). 265 https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.

Comparison to temperature elevation lapse rates
Melting on glacier surfaces across the HKH is controlled by the surface energy balance (SEB) between the atmosphere and underlaying snow, firn or ice. We explore the relationship between the Sentinel-1 SAR derived surface melting record and air temperature since measurements of SEB components like radiant energy fluxes are scarcely measured in the HKH and 275 difficult to extrapolate (de Kok et al., 2020). Air temperature across elevations in the Central Himalayas are determined during https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
Temperature-elevation lapse rates were determined using three-day averages of hourly air temperature measurements that were interpolated to fill gaps using methods identical for the calculation of temperature elevation lapse rates in numerical model studies of snowmelt and glacier wasting in the HKH (Baral et al., 2014). We calculated the difference between three-day 280 average air temperatures and divided by the difference in elevation (1,148m) between the two stations in the Langtang River Valley, Nepal. Lapse-rates ranged from 5ºC km -1 in July of 2018 to -13.7ºC km -1 in December of the same year. Temperature elevation lapse rates were used to extrapolate the maximum elevation of three isotherms (-10ºC, -5ºC, and 0ºC) for each day of year in 2018 in order to compare extrapolated temperatures with melt retrievals from Sentinel-1.

Results and Discussion
We report melt onset (MO), freeze onset (FO), and the frequency of melt retrievals across each study year, termed the annual melt fraction (AMF), across the HKH cryosphere for the calendar years 2017-2019. A melting signal is retrieved 290 across 97% (62,907km 2 ) of the mapped glacier area contained in the GAMDAM inventory. Between the Eastern (EH), Central (CH), Western Himalayas (WH) and Karakoram (K) we find similar mean glacier area with elevation, with the EH and CH having slightly higher minimum elevations compared to K and WH (Fig. 4b). In a comparison of the strength of the observed melting signal, we find that the EH and CH show a maximum at close to 6,000m a.s.l. compared to the WH and K at much lower elevations (~5,000m a.s.l.) (Fig. 4c). The maximum z indicates the greatest seasonal contrast in C-band backscatter. 295 Since bright radar signatures and strong seasonal contrasts during melting are indicative areas of percolating snowmelt, these results suggest a percolation zone skewed towards higher elevations in the EH and CH relative to the WH and K where these zones appear distributed over a larger range in elevation.
We detect melt earlier in the year and for a longer portion of each year in the central and eastern Himalaya compared to the western Himalaya and Karakoram regions. Glacio-climatic dynamics characteristic of HKH sub-regions are apparent in 300 summary statistics of melt retrievals: melt onset occurs about one week (6 days) earlier at elevations between 3,500m and 7,000m a.s.l. across the CH and EH relative to the WH and Karakoram (Table 2). Freeze onset in the CH and EH similarly occurs one week (7 days) sooner relative to the WH and Karakoram on average between 3,500m and 7,000 m a.s.l. The AMF https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
is only about 2.5% greater on average in the EH and CH compared to the WH and Karakoram between elevations of 3,500m and 7,000m a.s.l. At elevations below 4,000m a.s.l. across the entire HKH region, melt is retrieved across 79% of mapped 305 glaciers versus 96% of mapped glacier area above 4,000m a.s.l. in elevation. The reduced consistency in glacier melt-area retrievals at elevations where higher temperatures are expected is likely due to a thinner seasonal snowpack and reduced magnitude of the melting signal (reduced z). Here seasonal backscatter changes do not have a large volume scattering component. Considering this, records from these elevations that are missing or indicating anomalous frozen conditions at elevations below 4,000m a.s.l. are likely errors rather than actual non-melting conditions. Similar mischaracterization of 310 seasonal melt onset at low elevations (<4,000m a.s.l.) we interpret as a contribution from the debris-covered surfaces of the glacier exhibiting a radiometric response characteristic of freeze/thaw processes in soils (Entekhabi et al., 2010).
A summary of glacier melt timing with elevation averaged across study years is shown across 100m elevation bins in Figure 5 and tabulated across 1km elevation bins in Table 2. Due to errors in melt classification at elevations below 3,500m a.s.l., we will summarize our observations across elevations >3,500m a.s.l. and assume that the relatively small glaciated area 315 at elevations at or below 3,500m a.s.l. is negligible for the purposes of identifying trends across regions and elevation. We find an earlier MO at lower elevations, as expected from temperature-elevation lapse rates, beginning in mid-March in the eastern Himalaya and progressing to the Western, Central and Karakoram regions into late April (Fig. 5). Above 5,000m a.s.l. there is less variability in MO timing with elevation, but there is distinct regional variability where melt onset occurs on average 5 days sooner in the CH and EH regions compared to the WH and K regions. MO in the EH and CH between elevation ranges 320 of 3,500m -5,000m a.s.l. occurs one week sooner (7 days) than the WH and K regions. Freeze onset shows greatest regional variability at elevation ranges between about 5,000m and 6,500m a.s.l., occurring nine days later on average in the CE and EH compared to the WH and Karakoram. At elevations greater than 7,000m a.s.l., both melt and freeze onset appear to occur contemporaneously across glacio-climate subregions; with less than three days of difference between the WH/Karakoram and the CE/EH. 325 https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.  Table 2. Melt onset (MO), freeze onset (FO), and annual melt fraction (AMF) averaged over all three study years (2017-2019) across glacioclimatic sub-regions identified in the Randolph Glacier Inventory summarized across 1000m elevation bins. Melt retrieval generally increases with elevation due to the homogeneity and strong frozen scattering across regions of glacier percolation. Standard deviation for melt and freeze statistics could not be reported at elevations greater than 7,000m due to limited sampling frequency and are instead reported 335 as missing data (N/A).

Radar Scattering and Glacier Facies
Imaging radar backscatter intensity is linked with glacier facies. In the ablation zone, supra-glacier meltwater features 340 like crevasses, sun-cups, debris-cover, and other heterogeneities are likely to cause highly variable radar scattering mechanisms over short distances (Rott and Mätzler, 1987). Average z is minimum in the HKH across the lowest elevation glacier surfaces (3,000m-4,000m a.s.l.) and we subsequently retrieve valid melt signals at the lowest rate over elevations below 4,000m a.s.l.
There is much greater density of valid retrievals in fields above 4,000m a.s.l. (Table 2). These differences are apparent in the spatial granularity of the S1-SAR product, as shown in Fig. 6. Ablation zone surfaces on valley glaciers are well resolved and 345 show spatial heterogeneity in MO indicative of meltwater features and debris rather than randomly distributed noise. We attribute missing retrievals in these areas to lower sensitivity to melting from a lack of winter-season volume scattering due to differences in snow depth and/or morphology. We find these are more likely to occur in some regions of the HKH like the EH.
Conversely, there are surface-area fractions highly sensitive to surface melting that appear to be consistent with expectations of temperature lapse rates (i.e. earlier melting and later refreeze at lower elevations). We obtain more robust estimates of 350 melting onset and refreeze by spatially aggregating results of the glacier surface melt timing (Eq. 1) using a median window filter of 9x9 pixels after melt classification and z-score validation. In this way, we select signals within those areas that are indicative of seasonal melting. All spatiotemporal characteristics we report herein are after median window filtering of melt retrievals from 10m native resolution to 90m resolution.   We observe that the average MO is found to follow the 0º and -5 o C isotherms for elevations ~4500 to 6500 m a.s.l.
Below and above these elevations, and for FO, we find episodic melting events occurring over a range of elevations. This is 370 especially apparent in the FO around day of year 270 where FO occurs within a roughly two week period across glaciers between 5,000m-7,500m a.s.l. MO and FO signals are retrieved on days and at elevations where temperatures do not exceed -10ºC, which strengthens and expands recent in situ observations on glacier melt at the Khumbu glacier in the Mount Everest region showing longwave incident radiation driving melt at these temperatures and elevations (Matthews et al., 2019). Here we observe that, even at these extreme elevations (>7,000m a.s.l.) melt signals persist for over three months on average across 375 the Central Himalaya, which suggests that liquid water is retained at these high elevations across a seasonal melt cycle and may not be hydrologically negligible. In radar-derived observations we find surface melting to occur at elevations over 1km above the maximum elevation of 0ºC isotherms in the Central Himalaya. We observe melting over the vast majority of glacier surfaces; signals are retrieved across all elevation ranges of the HKH region, indicating that any areas of snow accumulation are experiencing some degree of surface melting during the summer months, purely melt-free, "dry snow" areas may not 380 currently exist in the HKH. https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.

Percolation Zone Meltwater Hydrology
We observe a delayed freeze-up over zones of known and apparent accumulation, indicating percolation and seasonal storage of subsurface meltwater at high elevations. An example of this is given in Figure 8, where refreeze occurs over thirty days later at 6,000m. a.s.l. compared to elevations around 6,600m a.s.l. The time series of mean Sentinel-1 SAR backscatter 385 for descending orbital nodes from two 250m buffered points on the Khumbu glacier show a more rapid brightening for the higher elevation location, whereas backscatter time series extracted from the point of lower elevation show a gradual increase in radar backscatter, indicative a gradually decreasing liquid water content in the snowpack (Fig. 8) (Forster et al., 2014;Miège et al., 2016). Sentinel SAR backscatter time series at the Khumbu Glacier on Mount Everest indicate similar spatial trends at higher (>6,500 m) and lower (<5,300m) elevations, but a refreeze anomaly is apparent within known elevations of regions of 390 meltwater percolation (Matthews et al., 2019). These data show gradual increase in backscatter at a lower elevation point (teal square, ~6,000m a.s.l.) and a point at higher elevation, above the maximum elevation of a three-day mean 0ºC summer isotherm (pink triangle, ~6,600m a.s.l.). This observation indicates that there exists a relationship between frozen percolation zone depth and C-band backscatter intensity across refreeze cycles wherein C-band backscatter gradually increases with frozen percolation zone depth during a refreeze process. 395 https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License. (1) that melt signals are pronounced above the maximum elevation of the 0ºC isotherm 400 extrapolated from meteorological stations within the Central Himalaya and (2) that melt signals persist for over one month at elevations of known meltwater percolation on the Khumbu glacier within ranges of ~5400-6300m a.s.l. relative to both higher and lower elevations. We hypothesize this observed delay in refreeze is due to the retention of liquid water in the glacier percolation zone and a gradual freeze up process.

Summary and Conclusions 405
Synthetic aperture radar time series backscatter images and glacier extent maps derived from optical imagery have long been proposed to inform hydrologic and glaciologic research across the global cryosphere, however a harmonized dataset of glacier surface melt does not exist. We retrieve glacier surface melt timing and duration for the study years 2017-19 across the HKH region using time series C-band SAR from the Sentinel-1 A/B satellites and an inventory of 76,831 glaciers spanning 65,108km 2 of ice-covered area. We quantify the magnitude of the seasonal melt signal by comparing mean summer and winter 410 backscatter using a z-score metric and apply a mask for areas with little confidence in seasonal melt signal (z < 2). We https://doi.org/10.5194/tc-2020-181 Preprint. Discussion started: 3 September 2020 c Author(s) 2020. CC BY 4.0 License.
subsequently retrieve melt characteristics across the vast majority 97% (62,907km 2 ) of HKH glacier area at 90m spatial and 12-day temporal resolution. Across elevations of the HKH region, melt onset generally increases with elevation, beginning in the first week of April at 3,000m a.s.l. and extending with increasing elevation into the second week of June at some of the highest elevations across the HKH (>7,000m a.s.l.). Across glacio-climatic regions of the HKH, melt onset occurs about one 415 week earlier on average in the Central and Western Himalaya compared to melt onset in the Western Himalaya and Karakoram.
Freeze onset similarly occurs one week sooner in the Western Himalaya and Karakoram compared to the Eastern and Central Himalaya at elevations with encompassing the majority of the glaciated area (4,000m-6,000m a.s.l.). Regional variations in melt timing and duration resemble known characteristics of regional glacio-climatic regimes. Importantly, where the presence of liquid water is likely driven by incident solar radiation driving melt at elevations where air temperatures do not exceed 0ºC, 420 the methods presented in this study can provide observational data on the presence of liquid water across data sparse regions for use in development and assessment of surface energy balance models of glacier ablation. Comparison of melt retrievals to temperature-elevation lapse rates calculated using two high-elevation meteorological stations in the Central Himalaya reveals that melt onset persists for over three months at elevations where extrapolated air temperature fields do not exceed -10ºC.
Additionally, we report observations and propose that meltwater retention is detectable in regions of known glacier percolation 425 through retrievals of delayed refreeze and that there may be a scaling relationship between SAR backscatter, frozen percolation zone depth and liquid water content. We produce a geospatial data product of melt onset (DOY), freeze onset (DOY), and fraction of Sentinel-1 observations classified as melt spanning glaciers of the HKH region at 90m spatial resolution for the calendar years 2017-2019 and plan to release annual updates to this database each calendar year. Melt is retrieved across all elevation ranges of HKH glaciers, which suggests that a dry snow accumulation zone in the HKH region is largely absent. The 430 methods presented in this study can provide the basis for an operational monitor of glacier melt dynamics and aid the development and assessment of surface energy balance models of glacier ablation across the global cryosphere.

Acknowledgements
This work was supported by funds provided to The City College of New York by the National Aeronautics and Space Administration Cryosphere Program's High Mountain Asia Team (HiMAT) program, under award number NNX16AQ83G. 435 Portions of this work were conducted at the Jet Propulsion Laboratory, California Institute of Technology, under contract to the National Aeronautics and Space Administration.