# Reversed Surface-Mass-Balance Gradients on Himalayan Debris-Covered Glaciers Inferred from Remote Sensing

Rosie R Bisset, Amaury Dehecq, Daniel N Goldberg, Matthias Huss, Robert

G Bingham, Noel Gourmelen

► To cite this version:

Rosie R Bisset, Amaury Dehecq, Daniel N Goldberg, Matthias Huss, Robert G Bingham, et al. Reversed Surface-Mass-Balance Gradients on Himalayan Debris-Covered Glaciers Inferred from Remote Sensing. Remote Sensing, 2020, 12 (10), pp.1563. 10.3390/rs12101563. hal-04945185

HAL Id: hal-04945185

https://hal.science/hal-04945185v1

Submitted on 13 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Reversed Surface-Mass-Balance Gradients on Himalayan Debris-Covered Glaciers Inferred from Remote Sensing

Rosie R. Bisset <fnref n="1" />,* ID, Amaury Dehecq <fnref n="2" />,<fnref n="3" /> ID, Daniel N. Goldberg ¹ID, Matthias Huss <fnref n="2" />,<fnref n="3" />,<fnref n="4" /> ID, Robert G. Bingham <fnref n="1" /> and Noel Gourmelen ¹ID

Received: 27 March 2020; Accepted: 7 May 2020; Published: 14 May 2020

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** Meltwater from the glaciers in High Mountain Asia plays a critical role in water availability and food security in central and southern Asia. However, observations of glacier ablation and accumulation rates are limited in spatial and temporal scale due to the challenges that are associated with fieldwork at the remote, high-altitude settings of these glaciers. Here, using a remote-sensing-based mass-continuity approach, we compute regional-scale surface mass balance of glaciers in five key regions across High Mountain Asia. After accounting for the role of ice flow, we find distinctively different altitudinal surface-mass-balance gradients between heavily debris-covered and relatively debris-free areas. In the region surrounding Mount Everest, where debris coverage is the most extensive, our results show a reversed mean surface-mass-balance gradient of -0.21 ± 0.18 m w.e. a⁻<fnref n="1" /> (100 m)⁻<fnref n="1" /> on the low-elevation portions of glaciers, switching to a positive mean gradient of 1.21 ± 0.41 m w.e. a⁻<fnref n="1" /> (100 m)⁻<fnref n="1" /> above an average elevation of 5520 ± 50 m. Meanwhile, in West Nepal, where the debris coverage is minimal, we find a continuously positive mean gradient of 1.18 ± 0.40 m w.e. a⁻<fnref n="1" /> (100 m)⁻<fnref n="1" />. Equilibrium line altitude estimates, which are derived from our surface-mass-balance gradients, display a strong regional gradient, increasing from northwest (4490 ± 140 m) to southeast (5690 ± 130 m). Overall, our findings emphasise the importance of separating signals of surface mass balance and ice dynamics, in order to constrain better their contribution towards the ice thinning that is being observed across High Mountain Asia.

**Keywords:** glaciers; surface mass balance; mass continuity; remote sensing; debris cover

## 1. Introduction

The glaciers in High Mountain Asia collectively form the largest glaciated area outside the polar regions, covering an estimated ~118,264 km<fnref n="2" /> [1]. Meltwater from these glaciers feeds into major river basins, including the Indus, Ganges, Brahmaputra, Yellow, and Yangtze, providing water resources to 221 ± 59 million people living in central and southern Asia [2]. Furthermore, the melt rates of the glaciers in High Mountain Asia modify the frequency and magnitude of glacial-lake outburst floods, posing significant threats to downstream communities [3–5]. Quantifying and improving our understanding of glacier surface-mass-balance distribution across High Mountain Asia is critical in the effective prediction and mitigation of these impacts.

<footnote n="1">¹ School of GeoSciences, University of Edinburgh, Edinburgh EH8 9XP, UK; Dan.Goldberg@ed.ac.uk (D.N.G.); r.bingham@ed.ac.uk (R.G.B.); Noel.Gourmelen@ed.ac.uk (N.G.)</footnote>

<footnote n="2">² Laboratory of Hydraulics, Hydrology and Glaciology, ETH Zurich, 8093 Zurich, Switzerland; adehecq@vaw.baug.ethz.ch (A.D.); matthias.huss@unifr.ch (M.H.)</footnote>

<footnote n="3">³ Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), CH-8903 Birmensdorf, Switzerland</footnote>

<footnote n="4">⁴ Department of Geosciences, University of Fribourg, 1700 Fribourg, Switzerland</footnote>

<footnote>* Correspondence: rosie.bisset@ed.ac.uk</footnote>
Glacier surface mass balance (hereafter SMB) is defined as the difference between accumulation and ablation, being negative where there is net melt. Many of the glaciers in High Mountain Asia are characterised by supraglacial debris cover of varying thickness and extent, which plays an important role in modifying SMB through its impact on glacier ablation rates. Field-based studies have shown that a thin layer of supraglacial debris, less than a critical thickness of ~3–8 cm, enhances glacier melt rates (e.g., [6]) through a reduction in the ice-surface albedo, as first demonstrated experimentally [7], and more recently constrained from surface energy balance modelling (e.g., [8,9]). In contrast, debris cover exceeding the critical thickness has the opposite effect, reducing ablation by insulating the ice surface (e.g., [6,10–14]). Other properties of supraglacial debris, such as moisture content, rock type, and grain size, can alter the thermal conductivity of the debris layer, consequently modifying the relationship between debris thickness and surface melt rates (e.g., [15–17]). The ablation rates on heavily debris-covered glaciers are extremely difficult to measure, due to the challenges that are associated with drilling stakes through the debris layer, as well as the large heterogeneity of local ablation rates (e.g., [13,18]).

Previous geodetic studies of ice-surface-elevation change have shown that, collectively, the glaciers in High Mountain Asia are rapidly losing mass [19–22], with a total annual mass change of $-19.0 \pm 2.5 \text{ Gt a}^{-1}$ between 2000 and 2018 [22]. However, despite the well-known importance of debris cover, geodetic studies have identified no clear relationship between supraglacial debris cover and rates of ice thinning in High Mountain Asia [19,20,23]. The presence of supraglacial features on debris-covered ice, such as ice cliffs and meltwater ponds, is one factor that is likely to be contributing to this anomaly. These features enhance localised melt rates and therefore partially offset the effects of reduced surface melting on ice thinning rates in debris-covered areas [24–26]. The reduced emergence velocity of debris-covered glaciers [27,28], leading to greater rates of ice-surface lowering compared to clean-ice glaciers, is another factor that is hypothesised to play a role in the debris anomaly. This surface lowering counteracts against the reduced ice-thinning rate in debris-covered areas, therefore potentially contributing towards the similar rates of thinning that have been observed on both debris-free and debris-covered glaciers [27–29].

The previously observed ice-surface-elevation changes are a result of a combination of both SMB and ice dynamics, as well as other processes, such as basal melting, internal accumulation, and calving (for lake-terminating glaciers) [30]. Therefore, producing regional distributed SMB observations is critical for separating signals of climate and ice flow dynamics and better constraining the region-wide influence of debris cover in glaciological models. Several studies have shown that it is possible to disentangle these contributions towards ice-surface-elevation changes using the principle of mass conservation [25,28,29,31,32]. Here, we present a methodology employing this principle in order to calculate SMB for a larger sample of glaciers across multiple regions, using recently produced spatially extensive remote sensing datasets. All of the mass-conservation terms are estimated from remote-sensing datasets [21,33] and modelled ice-thickness data [34]. After isolating the contribution of SMB towards ice thinning, we investigate the contrasting altitudinal patterns of SMB on debris-covered and clean-ice glaciers at a regional level.

## 2. Materials and Methods

### 2.1. Overview

Using a mass-continuity approach, we combined existing remote-sensing observations and modelled datasets to derive spatially-distributed estimates of mean SMB over a total glaciated area of ~2000 km² within the Pamir-Karakoram-Himalaya. Our computations were performed for 25 glaciers within these regions. We divided each glacier into sections of ~2 km in length (Figure A1). For each section, we computed the ingoing and outgoing ice fluxes from existing feature-tracking-derived ice velocities [33] and modelled ice thickness [34] (see Section 2.2). We calculated the mean elevation change for each section from previously derived geodetic measurements [18] and, based on mass
continuity, constrained the contribution of SMB. Using surface digital elevation models [35] and existing modelled debris-thickness maps [36], we analysed the dependence of SMB on elevation and debris thickness. We used breakpoint analyses in order to determine the elevations of transition points between different altitudinal SMB gradients for each glacier, and applied weighted regression models to determine the altitudinal SMB gradients below and above these elevations.

Our mean SMB estimates were generated for the period 2000–2015, which aligns with the periods of the input datasets (see Section 2.2). We carried out our analyses for five regions within High Mountain Asia (Figure 1), which were chosen based on the mutual availability of ice-surface velocity and elevation-change data, and to cover a range of meteorological and environmental conditions. Figure 1 shows the mean percentage debris cover for each region, as computed from an existing global debris-cover-extent dataset [37]. The largest glaciers within each region were included in our analyses (Figure A1), as smaller glaciers (<~10 km in length) do not provide a sufficient number of measurements to compute SMB gradients. The glacier-wide coverage of our SMB calculations was dependent on the availability of ice velocity data, which was often limited in coverage in high-altitude accumulation areas. We used our final SMB results to investigate the influence of supraglacial debris cover on altitudinal SMB gradients at a regional scale.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Regions of study across High Mountain Asia. Coloured boxes show the five regions for which surface mass balance is computed: Pamir (purple), Karakoram (blue), Spiti Lahaul (green), West Nepal (orange), and Everest (red). Pie charts on right show the mean percentage of ice covered by debris (computed from an existing debris cover distribution dataset for 2013–2017 [37]) for the glaciers analysed within each of the five study regions. Turquoise shaded areas show the glaciers from the Randolph Glacier Inventory v6.0 [1]. Blue outlined areas show the major river basins [38], which drain meltwater from the glaciers in High Mountain Asia.</figcaption></figure>

## 2.2. Data

The ice velocities used in this study cover the period 1999–2015 at 120 m resolution and they were produced using semi-automated feature tracking, applied to Landsat multispectral satellite imagery [see 33 for full derivation]. The ice-surface-elevation change dataset used in this study was derived by [21] from differencing of digital elevation models that were produced from ASTER optical satellite stereo imagery, and provides coverage for 2000–2016 at a 30 m resolution. We used this ice-surface-elevation dataset rather than the more recent compilation that was published in [22], because the former aligns more closely with the temporal coverage of the ice-velocity dataset. As a consequence, the SMB dataset that we produced in this study is representative of the mean SMB for the period 2000–2015, which is covered by both the ice-velocity and ice-surface-elevation-change datasets.
The modelled ice thicknesses used in this study are part of a global dataset, which was produced by [34] using an ensemble of models based on ice-thickness inversion from surface characteristics. The surface-elevation data used in this study are from SRTM 1 Arc-Second Global digital elevation data, which were collected in 2000 [35]. The debris-thickness dataset used in this study has a 30 m resolution and was modelled by [36] from thermal infrared Landsat 8 satellite imagery (2013–present). This dataset was selected for use as it is the only regional dataset providing consistent coverage of debris thickness across all five study regions.

## 2.3. Computing Cross-Sectional Ice Fluxes

Cross-sectional ice fluxes were computed from ice velocities [33] and ice thicknesses [34]. Cross-sectional transects between each glacier section were demarcated with two points at the lateral edges of the glacier, with transects being distributed at intervals of approximately 2 km along the length of each glacier. A spacing of approximately 2 km was chosen as a compromise between providing a sufficient number of data points to establish trends, and averaging over a sufficiently long along-flow distance to prevent large correlations and potentially strong correlation of velocity errors. Each transect was divided into 20 segments of equal length. The normal velocity and ice thickness were interpolated to the midpoint of each segment and then multiplied to compute flux. The resulting values were summed to compute the total ice flux perpendicular to each transect. These calculations were carried out twice for every transect, using median velocities for two periods (1999–2003 and 2013–2015), and the results averaged. These periods align with the temporal coverage of Landsat 7 (before the Scan Line Corrector failure) and Landsat 8 [33]. The formula for ice flux $Q_t$, where $t$ is 1999–2003 or 2013–2015, is thus

$$Q_t = \sum_{i=1}^{n=20} u_i h_i l_i \qquad (1)$$

where $i$ is the index of the segment, $u_i$ is the normal velocity interpolated to segment $i$, $h_i$ is the thickness interpolated to the segment, and $l_i$ is the length of the segment. The calculated mean values of $Q_t$ for each transect were used to quantify the ingoing and outgoing ice fluxes ($Q_{in}$ and $Q_{out}$) for each section between every adjacent pair of transects (or group of 3+ transects, where tributaries are present). The depth-averaged velocity depends on the fraction of basal sliding, which is unknown. Since the surface velocities are high for the observed glaciers, we assume that internal deformation makes a negligible contribution and, hence, that ice-surface velocity approximates to depth-averaged velocity. In order to test the impact of this assumption, we produced an additional set of surface-mass-balance estimates for the Gechongkang Glacier, based on a depth-averaged velocity equal to 90% of the surface velocity.

## 2.4. Producing Sectional Surface-Mass-Balance Estimates

The ice-surface area between each pair of transects was digitised and quantified using glacier outlines from the Randolph Glacier Inventory v6.0 [1]. The mean annual ice-surface-elevation change was calculated for each glacier section using glacier-wide elevation-change maps that were derived from digital elevation-model differencing [21]. The mean annual SMB between each pair/group of transects was computed using the mass-continuity method, as used by previous studies [25,28,29,31,32]:

$$b_j = \left( \frac{dh_j}{dt} - \frac{Q_{inj} - Q_{outj}}{A_j} \right) \times \frac{\rho_{ice}}{\rho_{water}} \qquad (2)$$

where $b$ is the mean SMB in glacier section $j$, $\frac{dh}{dt}$ is the mean ice-surface-elevation change, $A$ is the surface area of the section, $\rho_{ice}$ is the density of ice (920 kg m⁻³), and $\rho_{water}$ is the density of water (1000 kg m⁻³).
## 2.5. Statistically Approximating Altitudinal Gradients

For each glacier, we used breakpoint analysis to detect the transition points between contrasting altitudinal SMB gradients and to determine the elevation at which the transition point occurs. We used regression models, weighted by SMB uncertainties (discussed in Section 2.7), to estimate the altitudinal SMB gradients below and above computed breakpoint elevations. These gradients were compared to previously modelled SMB gradients from a regionally calibrated global glacier model used to estimate global glacier runoff changes [39]. We computed mean regional breakpoint elevations as the arithmetic average of breakpoint elevations for each individual glacier within the region. Similarly, we computed the mean regional altitudinal gradients above and below the breakpoint elevation using the arithmetic averages of the gradients for each glacier within each region. Using our collective dataset of sectional SMB values for each region, we used arithmetic averages to approximate the mean regional SMB, which was partitioned by elevation bands. For each glacier section, mean surface elevation and mean debris thickness were computed from digital elevation models [35] and a modelled debris-thickness dataset [36], respectively.

## 2.6. Estimating Equilibrium Line Altitudes

From our generated SMB results, we calculated the equilibrium line altitude (ELA) for each region, which describes the mean elevation at which accumulation and ablation are in balance [30]. Using our elevation-dependent SMB gradients, for each glacier we calculated the ELA value as the elevation at which SMB is equal to zero. We took the arithmetic average of the ELA values for all the glaciers analysed within each region to estimate the mean regional ELAs. We compared our regional ELA values to previous estimates [19,20,40].

## 2.7. Assessing Uncertainties

We assessed uncertainty in SMB through linear error propagation from errors in input data sets, as described by [41]. The error propagation involved an assessment of errors for all flux estimates via Equation (1). The uncertainties associated with the thickness and velocity components are 25% and 5–10 m a⁻¹, respectively. To estimate the error of flux, ε<sub>Qt</sub>, we assumed that both velocity and thickness uncertainty have multivariate Gaussian distributions, and that each have exponentially decaying autocorrelations with a length scale of L = 2 km. This was based on the autocorrelation of elevation differences from a previous geodetic study [20]; but overall, uncertainties were not found to depend strongly on L. Furthermore, thickness and velocity uncertainties were assumed to be independent. Thus, error in Q<sub>t</sub> is given by

$$Q = \sum_{i=1}^{n} (u_i \Delta h_i) l_i + (h_i + \Delta h_i) \Delta u_i l_i \quad (3)$$

where Δh<sub>i</sub> and Δu<sub>i</sub> are error in the interpolated velocity and thickness at a given segment i in the calculation of Q<sub>t</sub>. As this is a nonlinear expression of Δh<sub>i</sub> and Δu<sub>i</sub>, linear propagation cannot be applied. Therefore, we approximated this expression as

$$\varepsilon_{Q_t} = \sum_{i=1}^{n} (u_i \Delta h_i) l_i + (h_i + 3\sigma_{h_i}) \Delta u_i l_i \quad (4)$$

where σ<sub>h<sub>i</sub></sub> is the standard deviation of the error. Although this replaces a random error (Δh<sub>i</sub>) by a nonrandom term (3σ<sub>h<sub>i</sub></sub>), there is a 99.8% probability that the error term lies within the ±3σ interval [42].

Error was then propagated from the flux and elevation change estimates to SMB errors. The uncertainties in the computed breakpoint elevations were approximated as the standard errors that were computed from the weighted breakpoint analyses. The uncertainties in the computed altitudinal SMB gradients were approximated as the standard errors associated with the coefficients of the weighted regression models.
## 2.8. Validating Against Previous Findings

In order to validate our results, we firstly compared our computed SMB gradients directly to gradients previously calculated from in-situ measurements for the benchmark Chhota Shigri Glacier, in the Spiti Lahaul region. This glacier was selected for validation purposes due to data availability, meeting the minimum size requirements for inclusion in our analyses, and the existence of a record of altitudinally-varying ablation stake measurements covering a similar period to that of our study [10,43]. We are not aware of additional field-based altitudinal SMB gradients for any of the other glaciers included in our analyses and we were unable to compare directly to smaller glaciers for which previous gradients have been estimated from in-situ measurements since small glaciers provide an insufficient number of data points to compute reliable SMB gradients from our approach. Therefore, we made additional comparisons to previously computed gradients for glaciers within the vicinity of the glaciers in our analyses. These gradients were based on in-situ measurements at Abramov Glacier in the Pamir [44], as well as Pokalde Glacier [45,46] and Changri Nup Glacier [45] in the Everest region. For further validation, we compared our computed regional ELA values to previous values [19,20,40], which were estimated from different approaches, e.g., from snowline altitudes [20].

# 3. Results

## 3.1. Regional Surface-Mass-Balance Results

The mean regional SMB results, as partitioned by elevation band, indicate distinctly different elevation-dependent SMB patterns in West Nepal, in comparison to those that were observed in the other four regions (Figure 2). In West Nepal, the lowest mean SMB values (i.e., greatest average melt rates) occur at the lowest elevations (Figure 2d). In this region, the debris cover distribution is very low (Figure 1) and the mean modelled debris thickness is <3.5 cm in all elevation bands (Figure 2i).

In contrast, our results indicate that, in the Pamir, Karakoram, Spiti Lahaul, and Everest regions, the lowest mean SMB values do not occur at the lowest elevations. Instead, the lowest mean SMB values occur at mid-elevations (Figure 2). Within all four of these regions, there is a considerably higher debris-cover distribution when compared to West Nepal (see Figure 1). In addition, the mean debris thickness is the greatest in the lowest elevation band (modelled by [36] as > 30 cm) and decreases consistently with elevation.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Regional elevation-dependent trends in surface mass balance (SMB) and debris thickness. (a–e) Regional mean surface mass balance, binned by elevation, for Pamir (P), Karakoram (K), Spiti Lahaul (SL), West Nepal (WN), and Everest (E). Horizontal bars represent standard error values. (f–j) Corresponding mean debris-thickness values [36].</figcaption></figure>

## 3.2. Altitudinal Surface-Mass-Balance Gradients

We present the individual SMB estimates for every region in Figure A2. At a regional level, West Nepal shows a strong ($R^2$ = 0.65), linear altitudinal SMB gradient (Figure A2a). Meanwhile, in the Everest, Spiti Lahaul, Karakoram, and Pamir regions, no clear linear altitudinal gradients were detected at a regional scale ($R^2$ = 0.14 or less). However, breakpoint analysis of SMB gradients at a glacier-specific scale revealed that breakpoints in altitudinal SMB profiles commonly occur within the Everest, Spiti Lahaul, Karakoram, and Pamir regions (Table A1). Below the elevations where these breakpoints occur, reversed altitudinal SMB gradients (where SMB decreases with increasing elevation) were commonly found. More specifically, convex breakpoints in altitudinal SMB gradients were detected for 15 out of 18 of the analysed glaciers within these four regions. Of the 15 glaciers where breakpoints were detected, 13 show negative altitudinal SMB gradients below their breakpoint elevation. Conversely, in the West Nepal region, where debris cover is very minimal (Figure 1), breakpoints were not detected
for six of the seven analysed glaciers. Positive altitudinal gradients were found for all seven glaciers within this region (Table A1).

Figure 3 shows examples of the segmented altitudinal SMB gradients for the glaciers we analysed in the Everest region, with the corresponding debris-thickness values shown for each SMB data point. All five glaciers transition from negative altitudinal-SMB gradients of between −0.54 ± 0.19 and −0.06 ± 0.10 m w.e. a⁻¹ (100m)⁻¹ at low elevations to positive altitudinal-SMB gradients of between 1.01 ± 0.39 and 1.56 ± 0.49 m w.e. a⁻¹ (100m)⁻¹ at high elevations. The elevations at which convex breakpoints were detected range from 5180 ± 40 m a.s.l. (Ngozumpa Glacier) to 5920 ± 40 m a.s.l. (Gechongkang Glacier). The debris thickness decreases from low elevations to high elevations for every glacier within this region, as shown in Figure 3.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Elevation-dependent trends in SMB and debris thickness in the Everest region. Locations of the five largest glaciers surrounding Mount Everest (blue triangle). (b–f) Altitudinal surface-mass-balance trends for these glaciers. Vertical bars represent (1σ) uncertainties in the computed SMB values, estimated using a linear propagation of errors. Corresponding debris-thickness values [36] are represented on all panels with the colour scale on right. Dotted black lines represent our computed altitudinal SMB gradients, with corresponding R² values of the combined linear trends shown. The values of these gradients and their uncertainties are shown in Table A1. Blue dotted lines represent previously modelled SMB gradients from a regionally-calibrated global glacier model ignoring the effect of supraglacial debris [39].</figcaption></figure>

Table 1 shows the mean regional breakpoint elevations and SMB gradients, which were computed from the glacier-specific results for each region. In the Pamir, Karakoram, Spiti Lahaul, and Everest regions, we found regional mean SMB gradients of between −1.87 ± 2.57 and −0.17 ± 0.30 m w.e. a⁻¹ (100m)⁻¹ below breakpoint elevations, with a transition to mean gradients of between 0.94 ± 0.40
and 1.21 ± 0.41 m w.e. a⁻¹ (100m)⁻¹ above breakpoint elevations (Table 1). The results show that the regional breakpoint elevation varies considerably between regions. The lowest mean regional breakpoint elevation of 3680 ± 280 m a.s.l. is found in the Pamir, while the highest breakpoint elevation of 5520 ± 50 m a.s.l. is found in the Everest region.

<table><caption>**Table 1.** Regional mean altitudinal SMB gradients and breakpoint elevations.</caption><thead><tr><th rowspan="2">Region</th><th rowspan="2">Mean Breakpoint<br>Elevation (m)</th><th colspan="3">Mean Elevation-Dependent SMB Gradient (m w.e. a<sup>-1</sup><br>(100m)<sup>-1</sup>)</th></tr><tr><th>Below Breakpoint<br>Elevation</th><th>Above Breakpoint<br>Elevation</th><th>R<sup>2</sup> Value</th></tr></thead><tbody><tr><td>Pamir</td><td>3680 ± 280</td><td>-1.87 ± 2.57</td><td>1.01 ± 0.88</td><td>0.47</td></tr><tr><td>Karakoram</td><td>3860 ± 220</td><td>-1.04 ± 2.08</td><td>1.05 ± 0.59</td><td>0.22</td></tr><tr><td>Spiti Lahaul</td><td>4890 ± 100</td><td>-0.17 ± 0.30</td><td>0.94 ± 0.40</td><td>0.63</td></tr><tr><td>West Nepal</td><td>Not found</td><td>1.18 ± 0.40</td><td></td><td>0.74</td></tr><tr><td>Everest</td><td>5520 ± 50</td><td>-0.21 ± 0.18</td><td>1.21 ± 0.41</td><td>0.56</td></tr></tbody></table>

For the Pamir, Karakoram, Spiti Lahaul, and Everest regions, where breakpoints were detected for the majority of glaciers (Table A1), we display the mean altitudinal SMB gradients below and above the breakpoint elevation. These gradients represent the arithemetic average of glacier-specific gradients computed from breakpoint analyses (Table A1) within each region, and exclude glaciers where breakpoints were not detected. We also report the mean of the R² values, arising from the breakpoint analyses, from the same glaciers. For West Nepal, where breakpoints were not detected for the majority of glaciers, we report the mean altitudinal gradient (and its associated R² value), computed from the glacier-specific gradients, excluding the single glacier for which a breakpoint was detected. Error values reported are the mean standard error values associated with the gradients computed from linear regression (for West Nepal) and breakpoint analyses (for all other regions).

Our results also indicate considerable intra-regional variability in both the breakpoint elevation and the magnitude of the reversed gradients within each of the Everest, Spiti Lahaul, Karakoram, and Pamir regions, with the latter ranging over as much as two orders of magnitude (Table A1). This variability precludes any clear region-wide elevation-SMB relationships in these regions (as opposed to in West Nepal, Figure A2). However, there is more conformity in above-breakpoint mass-balance gradients (Table A1).

### 3.3. Equilibrium Line Altitudes

From our altitudinal SMB gradients, we estimate regional ELA values of between 4490 ± 140 m (Pamir) and 5700 ± 60 m (West Nepal), as shown in Figure 4. Our ELA values generally increase from the northwest to the southeast of the mountain belt. The exception to this trend is that we observe a slightly higher ELA in West Nepal than in Everest to its southeast. Figure 4 shows a comparison of our ELA values against previous estimates [19,20,40]. The uncertainty ranges of all of our estimated ELA values fall within the error bounds of the ELA values previous estimated from snow line altitudes [20] (Table A2).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Regional equilibrium line altitudes in Pamir-Karakoram-Himalaya. Coloured symbols show regional mean equilibrium line altitudes (ELAs) calculated by previous studies [19,20,40] and our study. Coloured shaded areas enclose the ELA values for each of the five study regions, from northwest (left) to southeast (right).</figcaption></figure>

# 4. Discussion

## 4.1. Influence of Supraglacial Debris Cover

The reversed altitudinal SMB gradients observed at low elevations in the Pamir, Karakoram, Spiti Lahaul, and Everest regions (Table 1; Figure 3) are likely to be attributable to the effects of supraglacial debris cover. Numerous previous studies have shown that supraglacial debris cover significantly influences glacier melt rates, as discussed in Section 1: a thin layer of debris enhances melt rates while a thicker layer of debris (exceeding a critical thickness of ~3–8 cm) reduces melt rates (e.g., [6–9]). Near the snouts of the glaciers in our analyses, where debris cover is thickest, the ice surface is likely to be insulated the most. As a result, the ablation rate is reduced and the SMB is raised (Figure 3). As elevation increases up-glacier from the snout, debris thickness decreases, therefore reducing the insulation effect. As a result, melt rates are enhanced up-glacier from the terminus, therefore contributing towards the reversed SMB gradients observed (Figure 3), which have also been reproduced in models (e.g., [47]). Furthermore, as the debris thins to below the critical thickness, the albedo effect is likely to dominate, therefore further enhancing melt rates and contributing towards the lowest mean SMB values being observed at mid-elevations (Figure 2).

The strong region-wide linear correlation between elevation and SMB for the relatively debris-free glaciers in West Nepal (Table 1; Figure A2) suggests that, where debris cover is largely absent, SMB is dominated by altitudinal climatic gradients. This pattern is also observed on other debris-free glaciers in the Himalayas (e.g., [43,45,46]). As a result, the lowest mean SMB values (i.e., greatest melt rates) occur at the lowest elevations in this region (Figure 2). The absence of negative altitudinal-SMB gradients in West Nepal also provides further evidence that the negative gradients observed in the Pamir, Karakoram, Spiti Lahaul, and Everest regions could be attributed to supraglacial debris cover.

The intra-regional variations in breakpoint elevations and altitudinal gradients observed within the Pamir, Karakoram, Spiti Lahaul, and Everest regions (Figure 3; Table A1) seem to preclude any clear region-wide relationships between elevation and SMB (Figure A2), despite strong correlation at the scale of individual glaciers. This variation could potentially arise from variations in supraglacial debris coverage. Differences in debris-thickness distribution, as well as local debris properties, such as lithology, grain size, and moisture content, may contribute towards the contrasting breakpoint
elevations and magnitudes of reversed altitudinal SMB gradients that were observed within these regions (e.g., [15–17]). Supraglacial ice cliffs and ponds can also influence the SMB of debris-covered glaciers by creating localised areas of enhanced melting [24–27]. Therefore, it is possible that these features could also partially explain the heterogeneity that was observed between debris-covered glaciers. Further detailed investigations of relationships between debris characteristics and altitudinal SMB gradients are required to gain a better understanding of these variations in reversed gradients.

## 4.2. Contribution of Glacier Dynamics

Our results enabled us to separate the relative contributions of SMB and ice dynamics towards ice-surface thinning. There is a considerable difference between our SMB estimates (black dotted line, Figure 5) and the previously observed ice-surface-elevation change (white dotted line). This difference represents the ice emergence velocity that was calculated from the mass convergence/divergence. At low elevations, where we observe a low ablation rate, the emergence velocity is also low (Figure 5), as also indicated by previous studies (e.g., [25,27,48,49]). This can be explained by the low slope and velocity gradient that are characteristic of the stagnant tongues of debris-covered glaciers. Further up-glacier near the breakpoint elevation, where melt rates are high, the emergence velocity is also high (see Figure 5). This is a consequence of the mass convergence at the transition between the stagnant tongue and the steeper debris-free part of the glacier. The reversed SMB gradient in the lower part of the glacier likely helps to maintain this mass convergence, by causing an inflexion of the glacier surface at its middle elevation. Consequently, the surface-elevation change, which is the sum of these two opposite processes of melt and emergence, shows a smooth trend along the length of the glacier, as observed by geodetic studies. Our results demonstrate the importance of ice dynamics in explaining the apparent contradiction between the reversed altitudinal SMB gradients expected on debris-covered ice and the relatively stable or positive altitudinal ice-surface-elevation-change gradients that were observed by geodetic studies.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of surface mass balance and ice-surface-elevation change. (a–e) show the contrasting altitudinal patterns of SMB (black) and ice-surface-elevation change (white) for the five largest glaciers surrounding Mount Everest. Each point represents the mean SMB / ice-surface-elevation change value computed for each ~2 km sub-section of each glacier. Vertical bars represent (1σ) uncertainties in the computed SMB values, estimated using a linear propagation of errors.</figcaption></figure>
### 4.3. Role of Glacier Surging

Referring to a recent inventory of surging glaciers [50], all of the glaciers that we have analysed in the Pamir and Karakoram are surge-type glaciers or have been exhibiting large velocity change during the study period. Surging behaviour can result in temporal variability in ice thickness and velocity (e.g., [51,52]); therefore, it is likely that this behaviour might be complicating the SMB trends and contributing towards the generally weaker $R^2$ values in the Pamir and Karakoram. An implicit assumption of our study is that the ice fluxes are in balance with the SMB gradient, but this is not the case for surging glaciers. Further investigation is required in the future in order to isolate signals of glacier surging, allowing for us to gain more accurate representations of altitudinal SMB gradients in regions where surge-type behaviour is occurring.

### 4.4. Validation of SMB Gradients With Previous In-Situ Measurements

For the benchmark Chhota Shigri Glacier in Spiti Lahaul, we computed an above-breakpoint altitudinal SMB gradient within 7–12% of the gradients previously measured in the field on the debris-free portion of this glacier (see Table 2). For Kangshung and Ngozumpa Glaciers in the Everest region, we computed above-breakpoint altitudinal gradients with differences of between 3% and 26% as compared to previously derived gradients for clean ice on West Changri Nup and Pokalde Glaciers, which are both located in between Kangshung and Ngozumpa Glaciers. For Fortambek Glacier in the Pamir, we computed an above-breakpoint gradient that was 27% lower than the field-based gradient for Abramov Glacier, located approximately 70 km northwest of Fortambek Glacier. The ELA values estimated from our SMB gradients are comparable to those calculated by previous studies [19,20,40], and follow a similar northwest-southeast increasing trend, as shown in Figure 4.

<table><caption>**Table 2.** Comparison of altitudinal SMB gradients with previous field-based findings.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Measurement Period</th><th rowspan="2">Elevation Range (m)</th><th colspan="3">Altitudinal Surface-Mass-Balance Gradient (m w.e. a<sup>-1</sup> (100m)<sup>-1</sup>)</th></tr><tr><th>Spiti Lahaul</th><th>Everest</th><th>Pamir</th></tr></thead><tbody><tr><td rowspan="4">Our study</td><td rowspan="4">2000–2015</td><td>4570–5180</td><td>0.74 ± 0.14 (CS)</td><td>-</td><td>-</td></tr><tr><td>5420–5750</td><td>-</td><td>1.12 ± 0.24 (K)</td><td>-</td></tr><tr><td>5180–5230</td><td>-</td><td>1.56 ± 0.49 (N)</td><td>-</td></tr><tr><td>3400–4020</td><td>-</td><td>-</td><td>0.74 ± 0.90 (F)</td></tr><tr><td>Azam et al. (2016) [43]</td><td>2002–2014</td><td>4400–5200</td><td>0.66 ± 0.09 (CS)</td><td>-</td><td>-</td></tr><tr><td>Wagnon et al. (2007) [10]</td><td>2002–2006</td><td>4400–4900</td><td>0.69 (CS)</td><td>-</td><td>-</td></tr><tr><td rowspan="2">Sherpa et al. (2017) [45]</td><td>2010–2015</td><td>5330–5690</td><td>-</td><td>1.47 (WCN)</td><td>-</td></tr><tr><td>2009–2015</td><td>5430–5690</td><td>-</td><td>1.37 (P)</td><td>-</td></tr><tr><td>Wagnon et al. (2013) [46]</td><td>2009–2012</td><td>5500–5600</td><td>-</td><td>1.51 (P)</td><td>-</td></tr><tr><td>Barandun et al. (2015) [44]</td><td>2011–2015</td><td>3720–4400</td><td>-</td><td>-</td><td>1.02 ± 0.05 (A)</td></tr></tbody></table>

The altitudinal SMB gradients computed above the breakpoint elevation for selected glaciers within the Spiti Lahaul, Everest and Pamir regions are shown. Gradients previously estimated by field-based studies [10,43–46] for debris-free ice in the ablation area are also shown for glaciers within these regions. The measurement periods and elevation ranges associated with each gradient are indicated. The initials provided after each SMB gradient represent Chhota Shigri (CS), Kangshung (K), Ngozumpa (N), West Changri Nup (WCN), Pokalde (P), Fortambek (F), and Abramov (A). WCN is located ~10 km east of N and ~15 km west of K. P is located ~15 km southeast of N and ~15 km southwest of K. A is located ~70 km northwest of F.

### 4.5. Limitations and Future Directions

While the ensemble approach used by [34] to model ice thickness minimises errors in the input ice-thickness dataset, there are still considerable uncertainties that are associated with ice thickness inversion [34]. These uncertainties are likely to have an impact on the accuracy of our results and we have accounted for these errors where possible. However, the improvement of ice thickness estimations in the future will increase the potential of this approach to produce better-resolved estimates of distributed SMB.
Our approach is also partly limited by the inability to sample SMB in high-altitude accumulation areas, due to larger uncertainties in ice velocities that are associated with feature-tracking in snow-covered areas [33]. Additionally, it is difficult to resolve accurately the SMB gradients of small glaciers < ~10 km in length, due to poorer signal-noise relationships and a lack of sufficient data points to accurately predict SMB gradients. The coverage of high-altitude areas and smaller glaciers can be improved in the future with the improvement of optical-satellite-imagery resolution [53] and the enhancement of feature-tracking algorithms.

There are also some uncertainties that are associated with the depth-dependence of ice velocities used to compute SMB, which can vary in space over a single glacier. Testing our depth-averaged velocity assumption for Gechongkang Glacier indicated that changing the depth-averaged velocity from 100% to 90% of the surface velocity resulted in an average change of 0.05 m w.e. $a^{-1}$ in our absolute SMB values (equivalent to an average change of 5%). While the effects of this assumption are relatively small, the associated uncertainties can be reduced in the future with better knowledge of spatially-distributed depth-dependences of mountain glacier velocities.

We recognise that there are some significant uncertainties that are associated with the debris-thickness model used in this study, associated with the coarse resolution of thermal infrared satellite imagery, as well as high temporal and spatial variabilities in surface temperatures and vertical debris-temperature profiles [32,54]. A further limitation is that since we compute the mean debris thickness for each glacier section, we do not account for small-scale variations in debris thickness, which have been demonstrated to often be highly heterogeneous (e.g., [55]). These local-scale variations lead to the formation of supraglacial features, which result in significant variations in melt rates over small spatial scales [25,26]. As previously discussed, it is likely that these localised melt variations may be contributing towards the heterogeneity in our SMB gradients. It is also possible that the debris-cover distribution and thickness may have evolved over the course of our study period [49,56]. However, significant changes generally occur over multi-decadal timescales [49] and, therefore, we consider it to be unlikely that the position of the transition zone between thin debris (below the critical thickness) and thicker debris significantly shifted during our study period.

A significant challenge that is associated with our approach is the highly branching nature of the glaciers that we analysed. Due to the presence of many glacier tributaries, only a semi-automated approach was possible, involving manual matching of ingoing/outgoing ice fluxes with corresponding glacier tributary sections. The development of a fully-automated approach, for example using flow lines to assign fluxes to corresponding tributary sections, would allow for a greater spatial coverage of distributed SMB and an expansion of our approach to further regions across High Mountain Asia.

Due to the uncertainties associated with the datasets produced in this study, we advocate that further field-based data acquisition is critical in order to validate remote-sensing-based observations thoroughly, and to facilitate accurate upscaling of SMB estimates over wider spatial scales.

## 5. Conclusions

In this paper, we have presented an approach for producing spatially distributed estimates of glacial surface mass balance from remote-sensing observations, based on the principle of mass continuity. We applied our approach to the largest glaciers within five key regions of the Pamir-Karakoram-Himalaya. Each glacier was divided into sections of approximately 2 km in length and we computed the ingoing and outgoing ice fluxes for each of these sections using satellite-derived ice velocities and modelled ice thicknesses. Using geodetic measurements, we calculated mean sectional ice-surface-elevation changes and subsequently isolated the contribution of surface mass balance towards the ice thinning rates using mass continuity. Using breakpoint analyses and regression models, we produced estimates of altitudinal surface-mass-balance gradients and equilibrium line altitudes for each of the five study regions. Our results show reversed altitudinal surface-mass-balance gradients in the lower-elevation portions of debris-covered glaciers, with a transition to positive surface-mass-balance gradients at higher elevations. In contrast, our results show continuously
positive altitudinal surface-mass-balance gradients on debris-free glaciers. This demonstrates that there are important differences in altitudinal ablation trends between debris-covered and debris-free glaciers, which were not previously visible from geodetic mass balance datasets. These differences in surface mass balance are likely offset by differences in ice dynamics, leading to similar thinning rates for debris-covered and clean-ice glaciers being observed by remote sensing. Our results show a regional equilibrium-line-altitude spatial gradient, with the values increasing from the northwest to the southeast. In future, with the generalisation and refinement of ice-velocity measurements and glacier thickness datasets, our operational approach can be developed and applied to glacierised mountain regions worldwide, providing the opportunity to uncover regional-scale surface-mass-balance patterns in areas where scale and location create challenges in field-based data acquisition.

**Supplementary Materials:** Supplementary materials are available online at http://www.mdpi.com/2072-4292/12/10/1563/s1.

**Author Contributions:** R.R.B., A.D., N.G. and D.N.G. designed the study. R.R.B. performed the processing, analyses and interpretation, with inputs from A.D., N.G., D.N.G., R.G.B. and M.H. The paper was written by R.R.B., with inputs from all co-authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research funded by the U.K. Natural Environment Research Council, under a studentship hosted by the Edinburgh U.K. Natural Environment Research Council E³ Doctoral Training Partnership (NE/L002558/1).

**Acknowledgments:** We thank Fanny Brun for her insights on debris-covered-glacier processes. We thank Philip Kraaijenbrink, Dirk Scherler, Daniel Farinotti, Fanny Brun and Heidi Sevestre for sharing their datasets. We thank the editor and three anonymous reviewers for their detailed comments that helped to improve this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Data availability:** The data produced in this study are available in the Supplementary Materials.

## Appendix A

<table><caption>Table A1. Glacier-specific altitudinal SMB gradients and breakpoint elevations for all glaciers analysed.</caption><thead><tr><th rowspan="2">Glacier</th><th rowspan="2">RGI 6.0 ID</th><th rowspan="2">Area (km²)</th><th rowspan="2">Elevation Range (m)</th><th rowspan="2">Breakpoint Elevation (m)</th><th colspan="3">Elevation-Dependent SMB Gradient (m w.e. a<sup>-1</sup> (100m)<sup>-1</sup>)</th></tr><tr><th>Below Breakpoint Elevation</th><th>Above Breakpoint Elevation</th><th>R<sup>2</sup> Value</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Pamir</td></tr><tr><td>Fortambek</td><td>13.19750</td><td>40.2</td><td>2930-6283</td><td>3400 ± 210</td><td>-2.04 ± 2.28</td><td>0.74 ± 0.90</td><td>0.27</td></tr><tr><td>Grumm-Grzhimaylo</td><td>13.13574</td><td>152.1</td><td>3615-6814</td><td>4020 ± 400</td><td>-1.77 ± 4.85</td><td>1.17 ± 1.12</td><td>0.24</td></tr><tr><td>Fedchenko</td><td>13.54431</td><td>663.7</td><td>2908-7392</td><td>3620 ± 240</td><td>-1.80 ± 0.58</td><td>1.13 ± 0.61</td><td>0.90</td></tr><tr><td>Garmo</td><td>13.19758</td><td>129.4</td><td>2976-6712</td><td>Not found</td><td>0.32 ± 0.13</td><td></td><td>0.34</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Karakoram</td></tr><tr><td>Hispar</td><td>14.04477</td><td>495.6</td><td>3110-7794</td><td>4230 ± 120</td><td>-0.62 ± 0.34</td><td>1.34 ± 0.72</td><td>0.23</td></tr><tr><td>Biafo</td><td>14.00005</td><td>559.8</td><td>3045-7148</td><td>3560 ± 310</td><td>-1.65 ± 4.00</td><td>1.50 ± 0.91</td><td>0.29</td></tr><tr><td>Baltoro</td><td>14.06794</td><td>809.1</td><td>3385-8569</td><td>3790 ± 230</td><td>-0.84 ± 1.89</td><td>0.30 ± 0.14</td><td>0.15</td></tr><tr><td>Siachen</td><td>14.07524</td><td>1078.0</td><td>3596-7579</td><td>Not found</td><td>-0.48 ± 0.06</td><td></td><td>0.75</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Spiti Lahaul</td></tr><tr><td>Samudra Tapu</td><td>14.15613</td><td>80.0</td><td>4237-6098</td><td>4750 ± 220</td><td>0.27 ± 0.48</td><td>0.97 ± 0.32</td><td>0.48</td></tr><tr><td>Bara Shigri</td><td>14.15447</td><td>112.4</td><td>3931-6309</td><td>5460 ± 90</td><td>0.02 ± 0.03</td><td>1.26 ± 0.75</td><td>0.58</td></tr><tr><td>Tichu Glacier</td><td>14.16068</td><td>24.8</td><td>4154-5875</td><td>4770 ± 50</td><td>-0.02 ± 0.08</td><td>0.78 ± 0.37</td><td>0.49</td></tr><tr><td>Sara Umga</td><td>14.16065</td><td>33.2</td><td>3835-5994</td><td>Not found</td><td>0.31 ± 0.10</td><td></td><td>0.39</td></tr><tr><td>Chhota Shigri</td><td>14.15990</td><td>16.8</td><td>4280-5764</td><td>4570 ± 50</td><td>-0.93 ± 0.59</td><td>0.74 ± 0.14</td><td>0.96</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">West Nepal</td></tr><tr><td>CN5O257E0002</td><td>15.11019</td><td>20.8</td><td>5032-6333</td><td>Not found</td><td>0.59 ± 0.17</td><td></td><td>0.75</td></tr><tr><td>CN5O257D0011</td><td>15.11003</td><td>28.6</td><td>5185-6635</td><td>Not found</td><td>1.41 ± 0.53</td><td></td><td>0.59</td></tr><tr><td>CN5O257D0018</td><td>15.10994</td><td>29.9</td><td>5107-6501</td><td>Not found</td><td>2.36 ± 0.86</td><td></td><td>0.78</td></tr><tr><td>CN5O257D0003</td><td>13.26906</td><td>6.5</td><td>5355-6466</td><td>Not found</td><td>0.71 ± 0.25</td><td></td><td>0.74</td></tr><tr><td>CN5O257D0006</td><td>15.11011</td><td>7.2</td><td>5285-6399</td><td>Not found</td><td>0.47 ± 0.17</td><td></td><td>0.80</td></tr><tr><td>CN5O257D0001</td><td>13.26909</td><td>3.6</td><td>5387-6222</td><td>Not found</td><td>1.54 ± 0.45</td><td></td><td>0.75</td></tr><tr><td>CN5O257D0005</td><td>13.26904</td><td>7.2</td><td>5308-6252</td><td>5590 ± 20</td><td>0.03 ± 0.41</td><td>4.73 ± 1.02</td><td>0.98</td></tr></tbody></table>
<table><caption>**Table A1. Cont.**</caption><thead><tr><th rowspan="2">Glacier</th><th rowspan="2">RGI 6.0 ID</th><th rowspan="2">Area (km²)</th><th rowspan="2">Elevation Range (m)</th><th rowspan="2">Breakpoint Elevation (m)</th><th colspan="3">Elevation-Dependent SMB Gradient (m w.e. a<sup>-1</sup> (100m)<sup>-1</sup>)</th></tr><tr><th>Below Breakpoint Elevation</th><th>Above Breakpoint Elevation</th><th>R<sup>2</sup> Value</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Everest</td></tr><tr><td>Rongbuk</td><td>15.09991</td><td>73.2</td><td>5155-7947</td><td>5680 ± 40</td><td>-0.18 ± 0.13</td><td>1.05 ± 0.20</td><td>0.72</td></tr><tr><td>Ngozumpa</td><td>15.03473</td><td>61.1</td><td>4702-8181</td><td>5180 ± 40</td><td>-0.54 ± 0.19</td><td>1.56 ± 0.49</td><td>0.51</td></tr><tr><td>Kangshung</td><td>15.09803</td><td>64.4</td><td>4587-8799</td><td>5420 ± 50</td><td>-0.06 ± 0.10</td><td>1.12 ± 0.24</td><td>0.70</td></tr><tr><td>Bhote Koshi</td><td>15.03422</td><td>27.2</td><td>4787-6550</td><td>5390 ± 70</td><td>-0.09 ± 0.09</td><td>1.33 ± 0.71</td><td>0.40</td></tr><tr><td>Gechongkang</td><td>15.09921</td><td>47.1</td><td>5328-7927</td><td>5920 ± 40</td><td>-0.20 ± 0.39</td><td>1.01 ± 0.39</td><td>0.47</td></tr></tbody></table>

<table><caption>**Table A2.** Regional ELA estimates.</caption><thead><tr><th rowspan="2"></th><th colspan="4">ELA Estimates (m)</th></tr><tr><th>Our Study</th><th>Gardelle et al. (2013)</th><th>Kääb et al. (2012)</th><th>Scherler et al. (2011)</th></tr></thead><tbody><tr><td>Pamir</td><td>4490 ± 140</td><td>4580 ± 250</td><td>-</td><td>-</td></tr><tr><td>Karakoram</td><td>4720 ± 240</td><td>5030 ± 280</td><td>5540</td><td>4845</td></tr><tr><td>Spiti Lahaul</td><td>5170 ± 140</td><td>5390 ± 140</td><td>5500</td><td>5103</td></tr><tr><td>West Nepal</td><td>5700 ± 60</td><td>5590 ± 138</td><td>-</td><td>-</td></tr><tr><td>Everest</td><td>5690 ± 130</td><td>5840 ± 320</td><td>5550</td><td>5700</td></tr></tbody></table>

Estimated ELAs (and standard errors) computed in this study are shown for each of the five study regions. Other previous estimates for these regions [19,20,40] are also shown.

<figure><img src="image_9.png" /></figure>

**Figure A1.** Detailed breakdown of study sites and glacier sampling method. The coloured areas in (c-g) show the glacier areas which were monitored in each of the five study regions shown in (a). Labels P, K, SL, WN and E denote Pamir, Karakoram, Spiti Lahaul, West Nepal and Everest. (b) shows an example of the flux gates (black cross-sectional lines) used for Baltoro Glacier, shown in subset K1 in d. For the area between each set of two or more flux gates (depending on the presence of tributaries), SMB is computed using the mass continuity method. Background DEMs are from Shuttle Radar Topography Mission [35].
<figure><img src="image_10.png" /></figure>

**Figure A2.** Detailed breakdown of glacier-specific SMB observations within each region. Elevation-dependent SMB observations for every glacier analysed within West Nepal (WN), Everest (E), Spiti Lahaul (SL), Karakoram (K) and Pamir (P) are shown in (a–e). Each glacier within each region is shown in a different colour. RGI codes for each glacier and its corresponding colour are shown within the legend in (f). Error bars show uncertainties associated with SMB estimates. Dashed grey lines show regional simple linear regression trends, which were calculated from all points combined within each region and weighted by SMB uncertainties. The R² values associated with each linear regression trend are also shown.

## References

1. RGI Consortium. Randolph Glacier Inventory—A Dataset of Global Glacier Outlines: Version 6.0. Available online: https://doi.org/10.7265/N5-RGI-60 (accessed on 3 March 2020).

2. Pritchard, H.D.; Arthern, R.J.; Vaughan, D.G.; Edwards, L.A. Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets. *Nature* **2009**, *461*, 971–975. [CrossRef]

3. Richardson, S.D.; Reynolds, J.M. An overview of glacial hazards in the Himalayas. *Quat. Int.* **2000**, *65–66*, 31–47. [CrossRef]

4. Veh, G.; Korup, O.; von Specht, S.; Roessner, S.; Walz, A. Unchanged frequency of moraine-dammed glacial lake outburst floods in the Himalaya. *Nat. Clim. Change* **2019**, *9*, 379–383. [CrossRef]

5. Veh, G.; Korup, O.; Walz, A. Hazard from Himalayan glacier lake outburst floods. *Proc. Natl. Acad. Sci. USA* **2020**, *117*, 907–912. [CrossRef] [PubMed]
6. Kayastha, R.B.; Takeuchi, Y.; Nakawo, M.; Ageta, Y. Practical prediction of ice melting beneath various thickness of debris cover on Khumbu Glacier, Nepal, using a positive degree-day factor. *IAHS Publ.* **2000**, *264*, 71–81.

7. Östrem, G. Ice melting under a thin layer of moraine, and the existence of ice cores in moraine ridges. *Geogr. Ann.* **1959**, *41*, 228–230. [CrossRef]

8. Lejeune, Y.; Bertrand, J.-M.; Wagnon, P.; Morin, S. A physically based model of the year-round surface energy and mass balance of debris-covered glaciers. *J. Glaciol.* **2013**, *59*, 327–344. [CrossRef]

9. Rounce, D.R.; Quincey, D.J.; McKinney, D.C. Debris-covered glacier energy balance model for Imja–Lhotse Shar Glacier in the Everest region of Nepal. *Cryosphere* **2015**, *9*, 2295–2310. [CrossRef]

10. Wagnon, P.; Linda, A.; Arnaud, Y.; Kumar, R.; Sharma, P.; Vincent, C.; Pottakkal, J.G.; Berthier, E.; Ramanathan, A.; Hasnain, S.I.; et al. Four years of mass balance on Chhota Shigri Glacier, Himachal Pradesh, India, a new benchmark glacier in the Western Himalaya. *J. Glaciol.* **2007**, *53*, 603–611. [CrossRef]

11. Pratap, B.; Dobhal, D.P.; Mehta, M.; Bhambri, R. Influence of debris cover and altitude on glacier surface melting: A case study on Dokriani Glacier, Central Himalaya, India. *Ann. Glaciol.* **2015**, *56*, 9–16. [CrossRef]

12. Nicholson, L.; Benn, D.I. Calculating ice melt beneath a debris layer using meteorological data. *J. Glaciol.* **2006**, *52*, 463–470. [CrossRef]

13. Vincent, C.; Wagnon, P.; Shea, J.M.; Immerzeel, W.W.; Kraaijenbrink, P.; Shrestha, D.; Soruco, A.; Arnaud, Y.; Brun, F.; Berthier, E.; et al. Reduced melt on debris-covered glaciers: Investigations from Changri Nup Glacier, Nepal. *Cryosphere* **2016**, *10*, 1845–1858. [CrossRef]

14. Azam, M.F.; Wagnon, P.; Patrick, C.; Ramanathan, A.; Linda, A.; Singh, V.B. Reconstruction of the annual mass balance of Chhota Shigri Glacier, Western Himalaya, India, since 1969. *Ann. Glaciol.* **2014**, *55*, 69–80. [CrossRef]

15. Nakawo, M.; Young, G. Field experiments to determine the effect of a debris layer on ablation of glacier ice. *Ann. Glaciol.* **1981**, *2*, 85–91. [CrossRef]

16. Lambrecht, A.; Mayer, C.; Hagg, W.; Popovnin, V.; Rezepkin, A.; Lomidze, N.; Svanadze, D. A comparison of glacier melt on debris-covered glaciers in the northern and southern Caucasus. *Cryosphere* **2011**, *5*, 525–538. [CrossRef]

17. Nicholson, L.I.; Benn, D.I. Properties of natural supraglacial debris in relation to modelling sub-debris ice ablation. *Earth Surf. Processes Landforms* **2012**, *38*, 409–501. [CrossRef]

18. Mihalcea, C.; Mayer, C.; Diolaiuti, G.; Lambrecht, A.; Smiraglia, C.; Tartari, G. Ice ablation and meteorological conditions on the debris-covered area of Baltoro Glacier, Karakoram, Pakistan. *Ann. Glaciol.* **2006**, *43*, 292–300. [CrossRef]

19. Kääb, A.; Berthier, E.; Nuth, C.; Gardelle, J.; Arnaud, Y. Contrasting patterns of early twenty-first-century glacier mass change in the Himalayas. *Nature* **2012**, *488*, 495–498. [CrossRef]

20. Gardelle, J.; Berthier, E.; Arnaud, Y. Region-wide glacier mass balances over the Pamir-Karakoram-Himalaya during 1999–2011. *Cryosphere* **2013**, *7*, 1263. [CrossRef]

21. Brun, F.; Berthier, E.; Wagnon, P.; Kääb, A.; Treichler, D. A Spatially resolved estimate of high mountain asia glacier mass balances from 2000 to 2016. *Nat. Geosci.* **2017**, *10*, 668–673. [CrossRef]

22. Shean, D.E.; Bhushan, S.; Montesano, P.; Rounce, D.R.; Arendt, A.; Osmanoglu, B. A systematic, regional assessment of High Mountain Asia Glacier mass balance. *Front. Earth Sci.* **2020**, *7*, 363. [CrossRef]

23. Brun, F.; Wagnon, P.; Berthier, E.; Jomelli, V.; Maharjan, S.B.; Shrestha, F.; Kraaijenbrink, P.D.A. Heterogeneous influence of glacier morphology on the mass balance variability in High Mountain Asia. *J. Geophys. Res. Earth Surf.* **2019**, *124*, 1331–1345. [CrossRef]

24. Sakai, A.; Nakawo, M.; Fujita, K. Melt rate of ice cliffs on the Lirung Glacier, Nepal Himalayas, 1996. *Bull. Glaciol. Res.* **1998**, *16*, 57–66.

25. Brun, F.; Wagnon, P.; Berthier, E.; Shea, J.M.; Immerzeel, W.W.; Kraaijenbrink, P.D.A.; Vincent, C.; Reverchon, C.; Shrestha, D.; Arnaud, Y. Ice cliff contribution to the tongue-wide ablation of Changri Nup Glacier, Nepal, Central Himalaya. *Cryosphere* **2018**, *12*, 3439–3457. [CrossRef]

26. Miles, E.S.; Watson, C.S.; Brun, F.; Berthier, E.; Esteves, M.; Quincey, D.J.; Miles, K.E.; Hubbard, B.; Wagnon, P. Glacial and geomorphic effects of a supraglacial lake drainage and outburst event, Everest Region, Nepal Himalaya. *Cryosphere* **2018**, *12*, 3891–3905. [CrossRef]

27. Banerjee, A. Brief communication: Thinning of debris-covered and debris-free glaciers in a warming climate. *Cryosphere* **2017**, *11*, 133–138. [CrossRef]
28. Nuimura, T.; Fujita, K.; Fukui, K.; Asahi, K.; Aryal, R.; Ageta, Y. Temporal changes in elevation of the debris-covered ablation area of Khumbu Glacier in the Nepal Himalaya since 1978. *Arct. Antarct. Alp. Res.* **2011**, *43*, 246–255. [CrossRef]

29. Nuimura, T.; Fujita, K.; Sakai, A. Downwasting of the debris-covered area of Lirung Glacier in Langtang Valley, Nepal Himalaya, from 1974 to 2010. *Quat. Int.* **2017**, *455*, 93–101. [CrossRef]

30. Cuffey, K.M.; Paterson, W.S.B. *The Physics of Glaciers*, 4th ed.; Academic Press: Amsterdam, The Netherlands, 2010.

31. Berthier, E.; Vincent, C. Relative contribution of surface mass balance and ice flux changes to the accelerated thinning of the Mer de Glace (Alps) over 1979–2008. *J. Glaciol.* **2012**, *58*, 1979–2008. [CrossRef]

32. Rounce, D.R.; King, O.; McCarthy, M.; Shean, D.E.; Salerno, F. Quantifying debris thickness of debris-covered glaciers in the Everest Region of Nepal through inversion of a subdebris melt model. *J. Geophys. Res. Earth Surf.* **2018**, *123*, 1094–1115. [CrossRef]

33. Dehecq, A.; Gourmelen, N.; Trouve, E. Deriving large-scale glacier velocities from a complete satellite archive: Application to the Pamir–Karakoram–Himalaya. *Remote Sens. Environ.* **2015**, *162*, 55–66. [CrossRef]

34. Farinotti, D.; Huss, M.; Fürst, J.J.; Landmann, J.; Machguth, H.; Maussion, F.; Pandit, A. A consensus estimate for the ice thickness distribution of all glaciers on earth. *Nat. Geosci.* **2019**, *12*, 168–173. [CrossRef]

35. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The shuttle radar topography mission. *Rev. Geophys.* **2007**, *45*, RG2004. [CrossRef]

36. Kraaijenbrink, P.D.A.; Bierkens, M.F.P.; Lutz, A.F.; Immerzeel, W.W. Impact of a global temperature rise of 1.5 degrees celsius on Asia's glaciers. *Nature* **2017**, *549*, 257–260. [CrossRef]

37. Scherler, D.; Wulf, H.; Gorelick, N. Global assessment of supraglacial debris-cover extents. *Geophys. Res. Lett.* **2018**, *45*, 11798–11805. [CrossRef]

38. The Global Runoff Data Centre. Major River Basins of the World. Available online: https://www.bafg.de/GRDC/EN/02_srvcs/22_gslrs/221_MRB/riverbasins_node.html (accessed on 9 March 2020).

39. Huss, M.; Hock, R. Global-scale hydrological response to future glacier mass loss. *Nat. Clim. Change* **2018**, *8*, 135–140. [CrossRef]

40. Scherler, D.; Bookhagen, B.; Strecker, M.R. Spatially variable response of Himalayan glaciers to climate change affected by debris cover. *Nat. Geosci.* **2011**, *4*, 156–159. [CrossRef]

41. Foresta, L.; Gourmelen, N.; Pálsson, F.; Nienow, P.; Björnsson, H.; Shepherd, A. Surface elevation change and mass balance of icelandic ice caps derived from swath mode CryoSat-2 altimetry. *Geophys. Res. Lett.* **2016**, *43*, 12138–12145. [CrossRef]

42. Abramowitz, M.; Stegun, I.A. *Handbook of Mathematical Functions*; Dover: New York, NY, USA, 1965.

43. Azam, M.F.; Ramanathan, A.L.; Wagnon, P.; Vincent, C.; Linda, A.; Berthier, E.; Sharma, P.; Mandal, A.; Angchuk, T.; Singh, V.B.; et al. Meteorological conditions, seasonal and annual mass balances of Chhota Shigri Glacier, Western Himalaya, India. *Ann. Glaciol.* **2016**, *57*, 328–338. [CrossRef]

44. Barandun, M.; Huss, M.; Sold, L.; Farinotti, D.; Azisov, E.; Salzmann, N.; Usubaliev, R.; Merkushkin, A.; Hoelzle, M. Re-analysis of seasonal mass balance at Abramov Glacier 1968–2014. *J. Glaciol.* **2015**, *61*, 1103–1117. [CrossRef]

45. Sherpa, S.F.; Wagnon, P.; Brun, F.; Berthier, E.; Vincent, C.; Lejeune, Y.; Arnaud, Y.; Kayastha, R.B.; Sinisalo, A. Contrasted surface mass balances of debris-free glaciers observed between the southern and the inner parts of the Everest Region (2007–2015). *J. Glaciol.* **2017**, *63*, 637–651. [CrossRef]

46. Wagnon, P.; Vincent, C.; Arnaud, Y.; Berthier, E.; Vuillermoz, E.; Gruber, S.; Ménégoz, M.; Gilbert, A.; Dumont, M.; Shea, J.M.; et al. The cryosphere seasonal and annual mass balances of Mera and Pokalde Glaciers (Nepal Himalaya) since 2007. *Cryosphere* **2013**, *7*, 1769–1786. [CrossRef]

47. Shea, J.M.; Immerzeel, W.W.; Wagnon, P.; Vincent, C.; Bajracharya, S. Modelling glacier change in the Everest Region, Nepal Himalaya. *Cryosphere* **2015**, *9*, 1105–1128. [CrossRef]

48. Anderson, L.S.; Anderson, R.S. Modeling debris-covered glaciers: Response to steady debris deposition. *Cryosphere* **2016**, *10*, 1105–1124. [CrossRef]

49. Rowan, A.V.; Egholm, D.L.; Quincey, D.J.; Glasser, N.F. Modelling the feedbacks between mass balance, ice flow and debris transport to predict the response to climate change of debris-covered glaciers in the Himalaya. *Earth Planet. Sci. Lett.* **2015**, *430*, 427–438. [CrossRef]

50. Sevestre, H.; Benn, D.I. Climatic and geometric controls on the global distribution of surge-type glaciers: Implications for a unifying model of surging. *J. Glaciol.* **2015**, *61*, 646–662. [CrossRef]
51. Meier, M.F.; Post, A. What are glacier surges? *Can. J. Earth Sci.* **1969**, *6*, 807–817. [CrossRef]

52. Paul, F. A 60-year chronology of glacier surges in the Central Karakoram from the analysis of satellite image time-series. *Geomorphology* **2020**, *352*, 106993. [CrossRef]

53. Millan, R.; Mouginot, J.; Rabatel, A.; Jeong, S.; Cusicanqui, D.; Derkacheva, A.; Chekki, M. Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach. *Remote Sens.* **2019**, *11*, 2498. [CrossRef]

54. Foster, L.A.; Brock, B.W.; Cutler, M.E.J.; Diotri, F. A physically based method for estimating supraglacial debris thickness from thermal band remote-sensing data. *J. Glaciol.* **2012**, *58*, 677–691. [CrossRef]

55. Nicholson, L.I.; McCarthy, M.; Pritchard, H.D.; Willis, I. Supraglacial debris thickness variability: Impact on ablation and relation to terrain properties. *Cryosphere* **2018**, *12*, 3719–3734. [CrossRef]

56. Kirkbride, M.P.; Deline, P. The formation of supraglacial debris covers by primary dispersal from transverse englacial debris bands. *Earth Surf. Process. Landf.* **2013**, *38*, 1779–1792. [CrossRef]

<figure><img src="image_11.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).