Debris Emergence Elevations and Glacier Change

Debris-covered glaciers represent potentially significant stores of freshwater in river basins throughout High Mountain Asia (HMA). Direct glacier mass balance measurements are extremely difficult to maintain on debris-covered glaciers, and optical remote sensing techniques to evaluate annual equilibrium line altitudes (ELAs) do not work in regions with summer-accumulation type glaciers. Surface elevation and glacier velocity change have been calculated previously for debris-covered glaciers across the region, but the response of debris cover itself to climate change remains an open question. In this research we propose a new metric, i.e. the debris emergence elevation (Z DE ), which can be calculated from a combination of optical and thermal imagery and digital elevation data. We quantify Z DE for 975 debris-covered glaciers in HMA over three compositing periods (1985–1999, 2000–2010, and 2013–2017) and compare Z DE against median glacier elevations, modelled ELAs, and observed rates of both mass change and glacier velocity change. Calculated values of Z DE for individual glaciers are broadly similar to both median glacier elevations and modelled ELAs, but slightly lower than both. Across the HMA region, the average value of Z DE increased by 70 +/− 126 m over the study period, or 2.7 +/− 4.1 m/yr. Increases in Z DE correspond with negative mass balance rates and decreases in glacier velocity, while glaciers and regions that show mass gains and increases in glacier velocity experienced decreases in Z DE . Regional patterns of Z DE , glacier mass balance, and glacier velocities are strongly correlated, which indicates continued overall increases in Z DEE and expansion of debris-covered areas as glaciers continue to lose mass. Our results suggest that Z DE is a useful metric to examine regional debris-covered glacier changes over decadal time scales, and could potentially be used to reconstruct relative mass and ELA changes on debris-covered glaciers using historical imagery or reconstructed debris cover extents.

On temperate mid-latitude glaciers, ELAs can be inferred from the elevation of the transient snowline at the end of the ablation season (Rabatel et al., 2005;Rabatel et al., 2013;Shea et al., 2013) or glacier reflectivity (Brun et al., 2015). In tropical/subtropical regions where accumulation and ablation occur simultaneously, the transient snowline is not correlated to the ELA (Ageta and Higuchi, 1984;Brun et al., 2015). Median glacier elevations (Z med ) are highly correlated with glaciologically determined ELAs (Carrivick and Brewer, 2004;Machguth et al., 2012), though the exact relation varies with glacier type and region (Braithwaite, 1984) and the presence of debris cover (Owen and Benn, 2005;Robson et al., 2016).
Debris-covered glaciers are defined as glaciers with a layer of supraglacial debris that covers some portion of their tongue (Cogley et al., 2011). They are common in many mountain regions (Alaska, European Alps, New Zealand, Coast Ranges of British Columbia, and the Andes), and in High Mountain Asia it is estimated that 10% of the total glacierized area is covered by debris (Scherler et al., 2011;Kraaijenbrink et al., 2017;Scherler et al., 2018). Debris is introduced to the glacier through either mass movements onto the glacier surface from adjoining rock walls or lateral moraines (rockfall, debris flow), or subglacial erosion (Kirkbride, 2011). Supraglacial or englacial debris is then transported down-glacier and becomes thickest at lower elevations where ablation rates are typically highest for debris-free glaciers (Benn and Ballantyne, 1994;van Woerkom et al., 2019). Changes in debris-covered area have been used to infer relative glacier status (Herreid et al., 2015) within HMA, but as glaciological mass balance measurements on debriscovered glaciers are notoriously difficult to maintain, there are few direct estimates of ELAs and glacier mass balance on debris-covered glaciers.
Observational and modelling studies demonstrate that the presence of a debris cover leads to longer glacier lengths and lower accumulation area ratios (AAR; (Scherler, Bookhagen and Strecker, 2011;Anderson and Anderson, 2016;Banerjee, 2017). In a warming climate, increased melt and decreased accumulation will lead to higher ELAs and transient decreases in AAR, reduced ice transport and glacier velocities, and increased ice ablation and mass loss (Shea and Immerzeel, 2016;Dehecq et al., 2019). We hypothesize that debris-covered glaciers with increased mass loss would experience accelerated debris melt-out and locally enhanced ablation rates. This will lead to an up-glacier expansion of the debris-covered area, and Z DE should increase. Conversely, debris-covered glaciers experiencing mass gain (Sakai and Fujita, 2017;de Kok et al., 2018;Farinotti et al., 2020) would exhibit increased ice transport from higher elevations and decreases in Z DE .
Expansion or contraction of Z DE will likely occur over decadal time scales (Anderson and Anderson, 2016) and changes in Z DE could thus be used to infer long-term trends in glacier mass balance as an alternative or complementary approach to geodetic methods based on DEM differencing. Using existing glacier inventories (Pfeffer et al., 2014) and a method for debris cover delineation (Kraaijenbrink et al., 2017), we construct composites of debris cover for three time periods and extract values of Z DE for individual debris-covered glaciers in HMA. We then examine relations between Z DE , median glacier elevations, modelled ELAs, and observed mass balances and velocity changes.

DATA AND METHODS
To examine changes in debris cover and Z DE , debris masks were calculated for three different time periods (1985-1999; FIGURE 1 | Study area map with Randolph Glacier Inventory regions (Table 1) and High Mountain Asia sub-regions (  (Pfeffer et al., 2014). Modelled (Kraaijenbrink et al., 2017) and observed mass balance rates (Brun et al., 2017) and trends in glacier velocity (Dehecq et al., 2019) for individual debris-covered glaciers were used to independently evaluate the trends in Z DE . While the total number of debris-covered glaciers in HMA is over 7,000, our analysis includes only 1) debris-covered glaciers larger than 5 km 2 to match the analysis of Dehecq et al., (2018), and 2) glaciers that have more than 10% of their area below the median elevation covered by debris.

Debris Cover Masks
Debris cover masks were mapped for all HMA glaciers using RGI glacier extents and a normalized difference snow index (NDSI; Hall et al., 1995) calculated from multi-year composites of Landsat imagery (www.usgs.gov). We preferentially select pixels with the highest surface brightness temperatures to construct the NDSI mosaics, which will reduce classification errors due to clouds, shadows, or the presence of snow. An NDSI threshold of 0.25 was then used to classify pixels within RGI glacier extents as either snow/ice (NDSI>0.25) or debris (NDSI<0.25). The classification procedure was performed in Google Earth Engine (GEE) and full details on the classification algorithm and accuracy are given in Kraaijenbrink et al. (2017; https://doi.org/10.5281/zenodo. 2548690). Glacier IDs were assigned to each debris-covered glacier in the analysis using the Randolph Glacier Inventory (RGI V5; RGI Consortium, 2017), and a resulting geodatabase was created and populated with data from this analysis and from other sources (see below). The average date and day of year (DOY) of acquisition for each glacier was also extracted from GEE and assigned to the geodatabase (Supplementary Figures  S1, S2), and used to assess temporal trends in Z DE and the possibility of bias related to the DOY of acquisition.

Median Glacier Elevation, Debris Emergence Elevation, and Modelled ELA
For each glacier retained in the analysis, a median glacier elevation (Z med ) was extracted from the gap-filled, 90 m resolution SRTM V3.0 elevation data that was collected in February 2000 (Farr et al., 2007). While SRTM accuracy decreases with increasing slope angle (Mukul et al., 2017), the debris-covered glaciers we focus on typically have lower angle slopes. To estimate Z DE , the boundary between debris-free and debris-covered ice was mapped as a polyline, and the median elevation of this contour was extracted from SRTM data ( Figure 2). Taking this conservative approach ensures that narrow medial moraines which can extend far up-glacier (or conversely, bare ice patches that extend downglacier) do not affect the determination of Z DE . Alternative approaches for estimation of Z DE (such as using percentiles of elevation for debris-covered pixels) were tested but these produced a wider range of possible values (Supplementary Material), and did not correspond to other ELA metrics (see below). A modelled ELA (ELA mod ) for each glacier, based on mass balance gradient constructed from a degree-day climatology, gridded precipitation (Kraaijenbrink et al., 2017), and the SRTM data, is included as a comparison against median glacier elevations and Z DE . Values of Z med , Z DE , and ELA mod were assigned to each debris-covered glacier in the geodatabase.

Rates of Change: Z DE
Rates of Z DE change (dZ DE /dt; m yr −1 ) were calculated for each glacier using a linear regression between the glacier-averaged dates of acquisition for the debris cover composites and the corresponding value of Z DE . Positive values of dZ DE /dt indicate an increase in Z DE and are likely associated with increased ELAs, more negative net mass balances, and decreased glacier velocities (Dehecq et al., 2019). Negative values of dZ DE /dt indicate decreased ELAs and positive mass balances or glacier surges.
To exclude both surge-type glaciers and other outliers from affecting our analysis, dZ DE /dt values greater than or less than Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 709957 two standard deviations from the mean (5% of all glaciers) were removed prior to analysis and plotting.

Rates of Change: Glacier Mass Balance
Rates of mass change between 2000 and 2016 have been calculated previously for 6,350 glaciers in the HMA region using time series of digital elevation models constructed from re-processed ASTER imagery (Brun et al., 2017). Linear regressions of surface elevation against time were performed for each 30-m glacierized pixel, and outliers were excluded to calculate trends in surface elevation (m yr −1 ). For each glacier, surface elevation trends were aggregated to 100 m elevation bands, and then converted to net glacier mass change with an assumed density of 850 kg m −3 (Huss, 2013). Net glacier mass balances calculated by Brun et al. (2017) were then merged with the debris-covered database of this study on the RGI V5 glacier ID.

Rates of Change: Glacier Velocity
For each debris-covered glacier in our database, we assign a change in glacier velocity using the dataset produced by Dehecq et al. (2019). Systematic surface feature tracking was applied to a database of Landsat-7 pairs between 1999 and 2017 to estimate gridded annual surface velocities at a grid resolution of 240 m. Trends in glacier velocity were then calculated from the annual velocity estimates, and extrapolated to the year 2000 to establish an initial velocity field. The median percent change in glacier velocity from 2000 to 2016 for each RGI glacier was then extracted from all grid points below the median glacier elevation, and we merge this final dataset with our geodatabase on the RGI ID. Full details on the calculation of glacier velocity changes and their uncertainties can be found in Dehecq et al. (2019).

Analyses
The total number of debris covered glaciers mapped for this analysis is N 7,457. We exclude debris-covered glaciers less than 5 km 2 to match the glacier velocity inventory of Dehecq et al. (2018), and we also exclude glaciers where the debris covered area is less than 10% of the area below the ELA, which leaves N 976. RGI regions 13 (South Asia East), 14 (South Asia West) and 15 (Central Asia) are used to broadly classify the studied glaciers, and reflect substantial differences in climate and topography ( Figure 1). We first examine distributions of Z med , Z DE , and ELA mod . Rates of Z DE change are then compared with observed rates of mass change within each region. We also compare the spatial pattern of dZ DE /dt with the spatial pattern of geodetic mass balances (Brun et al., 2017) and trends in glacier velocity (Dehecq et al., 2019). We then subdivide the RGI regions into smaller subregions based on mountain ranges identified in Kraaijenbrink et al. (2017) and compare sub-regional means.

Median Glacier Elevation, Debris Emergence Elevation, and Modelled ELA
Distributions of Z med , Z DE , and ELA mod (Figure 3) show strong similarities within all three RGI regions, which supports the use of dZ DE /dt as a metric of glacier change. Differences between regions highlight the climatic and geomorphic factors that determine glacier size and shape (

Changes in Z DE
In all RGI regions, the average value of Z DE increased from early to late periods (Table 1; Figure 4). To further quantify changes in Z DE through time for each debris-covered glacier, we calculate 1) simple differences in Z DE between periods (early to middle, and middle to late), and 2) the linear trend through time (dZ DE /dt; with units of m/yr) over all three compositing periods, using the average acquisition dates and the observed changes in Z DE at each glacier. In Central Asia, an average Z DE increase of 49 m was observed between the composite midpoints 1990 and 2014. In South Asia (West) and South Asia (East), average Z DE increases of 39 and 110 m were observed, respectively, over the same period. In South Asia (West), average values of Z DE increased 65 m between early and middle periods, and decreased 26 m between the 2005 and 2014 midpoints. This decrease in Z DE corresponds with the decadal scale mass changes and corresponding decreases in temperature and increases in precipitation observed in the region by Hugonnet et al. (2021).
Changes in Z DE between compositing periods ( Figure 4A) show a range of debris-covered glacier responses. For Central Asia and South Asia (East), changes in Z DE are more positive between the middle and late periods. A clear decrease in Z DE can be seen in South Asia (West) in the later period, which suggests either increased ice fluxes or decreased ablation on the glacier termini. Glaciers in this region, which includes the Karakoram and Kunlun Shan ranges, have well-documented anomalous mass gains or neutral mass balances (Hewitt, 2005;Brun et al., 2017;de Kok et al., 2018;Shean et al., 2020). As South Asia (West) contains a large number of surging glaciers (Hewitt, 2005;Quincey et al., 2011) it is also possible that our filtering process does not eliminate surging glaciers completely from the analysis. A glacier surge phase would likely affect Z DE and would be interpreted as a mass gain and velocity acceleration at the terminus. Median values of dZ DE /dt are positive (increasing Z DE over time) in all three HMA regions, with the greatest rates observed in South Asia East ( Figure 4B). Values of dZ DE /dt for individual glaciers range between −10 and +20 m/yr. These trends are in agreement with observed and projected rates of changes in both ELA and seasonal snowline elevations (Pandey et al., 2013;Hasson et al., 2014;Huss and Hock, 2015;Shea et al., 2015;Viste and Sorteberg, 2015;Pieczonka et al., 2018) Figure 5) in each RGI region. Glaciers in Central Asia and South Asia (West) exhibit statistically significant relations between net mass balance and Z DE trends (R 2 0.18-0.35), but this relation is weak or nonexistent in South Asia (East; R 2 0.00; Figure 5). Median glacier elevation and glacier area do not appear to affect the relation between dZ DE /dt and net mass balance. Similar relations are observed between Z DE and changes in glacier velocity ( Figure 5), though velocity changes and Z DE are not as closely linked as mass balance (R 2 0.08 for both Central Asia and South Asia West). Debris-covered glaciers that are stable or gaining mass typically show decreases in Z DE , and a mixed velocity response.
Glaciers in South Asia (East) have considerably thinned and retreated in recent decades, and have relatively fewer glaciers that exhibit mass gains. The lack of significant relations between mass balance and Z DE and glacier velocity and Z DE may simply indicate a sampling bias, or the fact that debriscovered glaciers in advanced state of recession will tend to have more variable changes in Z DE . Alternatively, the debris cover mask in South Asia (East) may be less reliable due to monsoon cloud cover which dominates South Asia (East) more than the other regions. Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 709957

Regional and Sub-Regional Trends
At the sub-regional scale (Supplementary Figures S4, S5, Table 2), relations between dZ DE , glacier mass balance, and glacier velocity for individual glaciers exhibit a greater degree of variability, yet the same basic relations exist. Sub-regions with positive net mass balances and increasing glacier velocities (Karakoram and Kunlun Shan W) notably have stable or decreasing Z DE through time. All other regions, which exhibit mass losses and velocity decreases, show increases in Z DE . When sub-regional means are calculated, strong links between Z DE, glacier mass balance, and glacier velocity are observed ( Figure 6). Spatial patterns of Z DE , glacier mass change, and glacier velocity (Figure 7) further demonstrate the linkages between these metrics. The observed trends in mass balance, velocity, and Z DE are nearly identical across the region, and demonstrate the potential use of Z DE to map historical patterns of glacier change. The most negative mass balance rates and the greatest increases in Z DE are seen in the eastern Himalayas. Regions of anomalous mass gain (Kunlun Shan and Karakoram) that have been identified previously (Hewitt, 2005;Gardelle et al., 2013;Kääb et al., 2015;Brun et al., 2017) coincide with regions of decreases in Z DE . Regional changes in Z DE thus provide a broad indication of glacier mass change, and there are very few outliers to this regional pattern.
Previous work has shown that there are no significant differences in the mass balance of debris-covered and debrisfree glaciers in the region . While debris insulates ice from melt, this is countered by low slope angles, low velocities, and low emergence velocities observed on debriscovered tongues (Brun et al., 2018). And though the glaciological ELA will respond to annual climate fluctuations, the terminus response times of debris-covered glaciers is substantially longer (Banerjee and Shankar, 2013). Thus, both the debris-covered area and changes in Z DE will lag behind annual variations in climate.

DISCUSSION
Debris covered glaciers in the HMA region represent a potentially significant source of meltwater on decadal/century time scales (Kraaijenbrink et al., 2017;Scherler et al., 2018;Herreid and Pellicciotti, 2020). The few studies that have examined the temporal evolution of debris cover have pointed out that a stable debris-covered area corresponds to a neutral or positive mass balance (Herreid et al., 2015).
This work demonstrates that the elevation of the debris-ice transition (Z DE ) has a similar distribution to both median glacier elevations and modelled equilibrium line altitudes (Table 1; Figure 3). As the ELA responds to annual variation in climate, we focus on how Z DE changes through time and compare this with published mass balance and velocity trends. Our discussion here focuses on the implications of our results, the limitations of the research, and opportunities for future research.
Debris covered glacier extents can be difficult to map with remotely sensed imagery (Frey et al., 2012). However, the transition between debris and clean glacier ice can be reliably estimated from cloud-free optical imagery (Rounce et al., 2021). While the elevation of this debris-ice transition will vary over individual glaciers (e.g. Figure 2), the median value of the transition elevation extracted here as Z DE has very similar distributions to both median glacier elevations and modelled ELAs (Figure 3). Limitations to our approach include the use of a static DEM and the duration of the compositing period used to construct debris cover maps. In the absence of accurate DEMs for each composite period, we have relied on SRTM elevations (ca. 2 | Regional and sub-regional trends in Z DE , mass balance, and velocity for debris covered glaciers with a total area greater than 5 km 2 . Values given are mean +/− standard deviation. 2000) for all estimates of Z DE . If we assume that the majority of glaciers in the region have been experiencing mass loss and surface lowering, our use of SRTM data for all periods may result in overestimates of Z DE in the 2013-2017 period, and underestimates in the 1981-2000 period. Nevertheless, we do not expect the reliance on a static DEM to affect our conclusions substantially, as these errors are likely small near the ELA (which is slightly higher in elevation than Z DE , Figure 3). This study relies on compositing periods that cover 7-15 years periods, in order to obtain sufficient snow-free imagery to map the full extent of debris cover accurately. Changes in Z DE at shorter timescales would likely not be detectable with reasonable accuracy, and given the absence of annual mass balance or velocity trends there is no need for finer temporal resolutions. Sub-regional means of Z DE exhibit clear negative relations with mass balance and velocity trends (Figure 7), with negative mass balances and velocity slowdowns associated with increased elevations of debris emergence. However, we do not find strong evidence for these relations at the scale of individual glaciers (Figures 6, 7). Glaciers in Central Asia (Tien Shan, Kunlun Shan, Pamir, and Tibetan ranges) appear to have the strongest links between these three glacier metrics, while the relation is weaker in South Asia (West) and South Asia (East). These differences are possibly the result of sampling bias, the relative health of glaciers in a region, the date of imagery acquisition, debris supply rates, and the density of laketerminating debris-covered tongues.

Region
As there are no debris-covered glaciers in South Asia (West) or South Asia (East) that are gaining mass, the first issue is sampling bias: we are only able to sample glaciers that are losing mass. This bias is further related to the general health of glaciers in a given region. Sustained negative mass balances in South Asia (West) and South Asia (East) have produced debris-covered glaciers that are in an advanced state of decline, and in some cases are being cut off from their accumulation areas (Naito et al., 2000;Shea et al., 2015;Pieczonka et al., 2018) and in transition to becoming rock glaciers (Anderson et al., 2018). As a result, Z DE will have an upper limit.
Our analysis of the dates of acquisition shows that imagery is typically collected in early August in most regions and composites (Supplementary Figure S1). However, we find a slight bias towards a later average glacier acquisition DOY for RGI15 in the 2013-2017 composite. While debris emergence elevations are expected to change on annual/decadal timescales in response to mass balance forcings, the timing of image acquisition could potentially affect the delineation of debris cover due to snow. Our compositing procedure preferentially selects pixels with the highest brightness temperatures to help minimize this bias. The average date of acquisition also varies between glaciers (Supplementary Figure S2). In using glacier-averaged composite dates for the trend calculations we find greater trend magnitudes than if the arithmetic mid-point of the composites is used.
Additionally, our calculated Z DE metric extends beyond the time periods of the mass balance and velocity data. Looking only at the absolute change in Z DE between the last two debris-cover masks, we observe a similar, albeit weaker, pattern in all three regions (Supplementary Figure S4). Glaciers that experience mass loss and decreases in velocity tend to show increases in Z DE . As debris-covered glaciers thin, there is the possibility of increased debris supply from unstable and de-buttressed lateral moraines (van Woerkom et al., 2019) which would act to increase Z DE . Finally, South Asia (West) and South Asia (East) contain a greater density of lake-terminating debris-covered glaciers (King et al., 2019), which have accelerated rates of mass loss relative to land-terminating glaciers. These factors all introduce some uncertainty in the relation between mass balance and Z DE .
This research is unable to resolve a mechanism for changes in Z DE , though enhanced ablation at the boundary between ice and debris, and a positive albedo feedback (warmer temperatures leading to increased ablation, melt out of entrained debris, FIGURE 6 | Sub-region means of Z DE versus mass balance rate (top panel) and Z DE vs. change in glacier velocity. Error bars are standard deviations for each sub-region and metric, as in Table 2, and points are scaled to the sample size.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 709957 decreased surface albedo, and enhanced melt) is our working hypothesis. Detailed energy balance modelling work at a regional scale would help test this hypothesis. A reduction in ice flux would point towards reduced debris transport, so the correlation between reduced ice fluxes and increase in Z DE may simply be due to the previously observed relations between mass balance and velocity (Dehecq et al., 2019). The conceptual model of supraglacial debris given by (Anderson and Anderson, no date) and (Bozhinskiy et al., 1986) suggest that englacial debris emergence at the surface is governed by the local ablation rate and the concentration of the debris within the ice. When debris begins to emerge and concentrate at the surface, ablation rates are reduced at debris thicknesses greater than 0.05-0.10 m, which leads to an increase in local relief of the emerging medial moraine. The increase in surface slopes promotes dispersal of the debris, and so increased mass loss promotes increased debris emergence.
As the velocity is due primarily to the driving stress (with glacier thinning and mass loss leading to reduced driving stress and reduced velocities), the link between changes in velocity and Z DE may be coincidental, and not causal. Put another way, the rise in Z DE and the reduction in horizontal velocities are consequences of the same mechanism but are not directly linked. The rise in Z DE can be seen as the consequence of an increased ELA, which leaves a larger area with positive vertical velocity (downglacier of the ELA), and thus leads to debris cover expansion up-glacier. The rise in ELA leads to more thinning and thus a reduced driving stress, as highlighted in Dehecq et al (2019).
Debris-covered glaciers respond to ELA shifts at decadal or century-scale timescales (Banerjee, 2017). Given the importance of larger debris-covered glaciers to local water supplies (Fegel et al., 2016), we split the dataset in each RGI region according to glacier size (greater than or less than 10 km 2 ) and re-calculate the mean rates of Z DE change to determine if glacier size affects the rate of Z DE change. In RGI13, larger glaciers have a slightly higher Z DE rate of change (+3.23 m yr −1 , N 81) than smaller glaciers (+2.68 m yr −1 , N 104), though the standard deviations on these estimates are large (5.90 and 4.86 m yr −1 , respectively). In South FIGURE 7 | Spatial distributions of dZ DE /dt (A), mass balance (B), and change in velocity (C) for debris covered glaciers >5 km 2 in HMA. Mass balance data and velocity data from (Brun et al., 2017;Dehecq et al., 2018), respectively.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 709957 Asia (West) and South Asia (East), smaller debris-covered glaciers have higher rates of Z DE change (3.15 and 6.23 m yr −1 , respectively) than larger debris-covered glaciers (2.69 and 5.24 m yr −1 , respectively). Again, large uncertainties prevent any meaningful conclusions about the role of glacier area in the rate of Z DE response. In light of these discussion points, our work indicates that changes in Z DE derived from historical airborne or spaceborne imagery could be used as indicators of past glacier change. Historical regional inventories of Z DE derived from Corona/ Hexagon imagery (Maurer et al., 2016;Bolch et al., 2017) would provide additional information for the development and testing of models of glacier dynamics and the response times of debris-covered glaciers to climatic change. Extension of this work also includes the potential for refinements to reconstructions of debris-cover extents and ELAs from the Last Glacial Maximum (Owen et al., 2002) or early Holocene (Fernández-Fernández et al., 2017) from geomorphological evidence.

CONCLUSION
Debris-covered glaciers are common in High Mountain Asia, and represent a potentially important long term source of meltwater. In this study, we map the elevation of the debrisice transition zone (Z DE ) for debris-covered glaciers in High Mountain Asia over three different compositing periods, identify trends and changes in Z DE between periods, and compare rates of glacier mass change and velocity change with changes in Z DE . Over decadal timescales and subregional spatial scales, Z DE is closely related to both glacier mass balance and glacier velocity. Glacier mass loss and reductions in velocity are coincident with increases in Z DE , and we hypothesize that this is likely a result of increased ablation at the debris-ice transition leading to increased melt-out of englacial debris. Rates of increases in Z DE are of the same magnitude as modelled and observed increases in ELA. Our results provide an avenue for mapping the relative health of debris-covered glaciers from historical imagery, which can be used to improve debris-covered glacier models and projections of future meltwater availability.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JS, WI, PK, and FB designed the study. JS and PK conducted the analyses and prepared the figures, and all authors contributed to the writing.