# Health and sustainability of glaciers in High Mountain Asia

Evan Miles, Michael Mccarthy, Amaury Dehecq, Marin Kneib, Stefan Fugger, Francesca Pellicciotti

► To cite this version:

Evan Miles, Michael Mccarthy, Amaury Dehecq, Marin Kneib, Stefan Fugger, et al.. Health and sustainability of glaciers in High Mountain Asia. Nature Communications, 2021, 12 (1), pp.2868. 10.1038/s41467-021-23073-4 . hal-04945318

HAL Id: hal-04945318

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

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

https://doi.org/10.1038/s41467-021-23073-4

OPEN

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

Check for updates

# Health and sustainability of glaciers in High Mountain Asia

Evan Miles ID <fnref n="1" />☒, Michael McCarthy ID <fnref n="1" />,², Amaury Dehecq ID <fnref n="1" />,³, Marin Kneib ID <fnref n="1" />,⁴, Stefan Fugger ID <fnref n="1" />,⁴ & Francesca Pellicciotti ID <fnref n="1" />,⁵

Glaciers in High Mountain Asia generate meltwater that supports the water needs of 250 million people, but current knowledge of annual accumulation and ablation is limited to sparse field measurements biased in location and glacier size. Here, we present altitudinally-resolved specific mass balances (surface, internal, and basal combined) for 5527 glaciers in High Mountain Asia for 2000-2016, derived by correcting observed glacier thinning patterns for mass redistribution due to ice flow. We find that 41% of glaciers accumulated mass over less than 20% of their area, and only 60% ± 10% of regional annual ablation was compensated by accumulation. Even without 21st century warming, 21% ± 1% of ice volume will be lost by 2100 due to current climatic-geometric imbalance, representing a reduction in glacier ablation into rivers of 28% ± 1%. The ablation of glaciers in the Himalayas and Tien Shan was mostly unsustainable and ice volume in these regions will reduce by at least 30% by 2100. The most important and vulnerable glacier-fed river basins (Amu Darya, Indus, Syr Darya, Tarim Interior) were supplied with >50% sustainable glacier ablation but will see long-term reductions in ice mass and glacier meltwater supply regardless of the Karakoram Anomaly.

<footnote n="1">¹Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf, Switzerland. ²British Antarctic Survey, Natural Environment Research Council, Cambridge, UK. ³Laboratory of Hydraulics, Hydrology and Glaciology, ETH Zurich, Zurich, Switzerland. ⁴Institute of Environmental Engineering, ETH Zurich, Zurich, Switzerland. ⁵Department of Geography, Northumbria University, Newcastle, UK. ☒email: evan.miles@wsl.ch</footnote>
Glaciers and snow in High Mountain Asia (HMA) release enough meltwater seasonally to meet the requirements of nearly a quarter of a billion people¹,², and basins fed by these glaciers are the most vulnerable worldwide to ongoing climatic, societal and environmental change³. Assessing the current state and future prevalence of ice and snow reservoirs in HMA is therefore a key priority⁴⁻⁶. However, access to the remote, high-altitude glaciers of HMA can be dangerous and time-consuming, which has restricted field observations of surface mass balance to a sparse coverage of very few sites, mostly at lower-altitude ablation areas⁷,⁸. Mass inputs to glaciers are generally unknown, due to the high uncertainty of reanalysis data at high altitudes and few direct observations⁹⁻¹¹. The current models of glacier change are thus overparameterized and unable to constrain key aspects of glaciers' internal dynamics and interactions with climate, such as the influences of avalanching and debris cover¹². The observational uncertainty in glacier state and varying process representation thus leads to considerable uncertainty in glacier volume change projections within HMA and globally⁶,¹³.

Recent remote sensing studies have advanced regional-scale understanding of the mass change and dynamics of HMA glaciers¹⁴⁻²⁰. However, even high-precision measurements of elevation change cannot resolve the spatial patterns of mass balance across individual glaciers²¹. This specific mass balance (SMB), sometimes called the climatic mass balance, is the combination of accumulation (mass gain) and ablation (mass loss) at a position on the glacier, and combines surface, englacial, frontal, and basal components. Local mass gains and losses due to accumulation and ablation are partially compensated by mass redistribution as a result of ice flow, leading to ice thinning (negative elevation change) or thickening (positive elevation change, Supplementary Figs. 1-2); rather than representing the local signals of SMB, elevation change measurements integrate ice motion²²,²³. In addition, average glacier mass balances estimated from thinning datasets alone typically use an area-averaged glacier surface density which may not be appropriate for glaciers severely out of balance with climate²⁴.

Knowledge of SMB is vital for understanding the regional and local drivers of glacier change⁷,²⁵, and for calibration and validation of numerical models to appropriately represent current and future glacier changes¹²,¹³,²⁶. As the basal and englacial components of SMB are often negligible, while frontal ablation is localized, observations and models often equate SMB with the surface mass balance. New generations of glacier models have incorporated improved spatial representations of glacier surface processes and ice dynamics¹³,²⁷, but without spatially resolved control datasets, these models are forced to calibrate to sparse, biased measurements of surface mass balance or spatially integrated signals of glacier change such as thinning or area change⁵,¹²,²⁸. The scarcity of in situ measurements is particularly problematic in HMA because many of the region's glaciers do not conform to the mass balance patterns assumed by regional-scale models. SMB is typically simplified to linear altitudinal gradients for the ablation and accumulation areas but the prevalence of avalanching²⁹,³⁰ and supraglacial rocky debris⁵,³¹ across the region may lead to distinctive mass balance profiles³², while mass accumulation rates above 6000 m a.s.l. are rarely measured⁹,³³⁻³⁵.

In this study, we provide spatially distributed SMB for glaciers across the entirety of HMA for the 2000–2016 period. We use this dataset to derive glacier mean equilibrium line altitudes (ELAs), the extent of accumulation areas, and the portion of annual ablation that is compensated by annual accumulation as indicators of glacier health in major river basins across the region. Finally, we assess the consequences of the glaciers' contemporary climatic-geometric imbalance in terms of implied changes in ice volume and discharge by 2100.

## Results and discussion

**Altitudinal SMB.** We determine SMB by solving the continuity equation, assuming that englacial and basal mass change is negligible (Eq. 1, “Methods”). We calculate the ice flux divergence by combining estimates of ice thickness³⁶ with observations of ice surface motion¹⁶ and a Monte Carlo-based estimate of depth-averaged correction factor. We use this with elevation change measurements¹⁵ to derive fully distributed and altitudinal SMB, carefully accounting for uncertainty associated with the input data and methods (Methods). Our results correspond to mean annual values for the 2000–2016 period, as constrained by the input elevation change and velocity observations.

The continuity equation has previously been used to determine the SMB on individual glaciers²²,²³,³⁷⁻⁴⁰ but never at a large scale, which requires an automated approach. We first apply our approach to 35 sites to compare to available measurements of surface mass balance⁴¹ (Supplementary Table 1), and another 25 glaciers for which remote sensing-derived estimates of SMB are available³⁸. Our method is consistent with 79% of field measurements to within 0.2 m w.e. a⁻¹ and generally reproduces observed mass balance patterns where glacier velocity is measurable (Supplementary Information). We thus apply this method to the 7341 glaciers in High Mountain Asia with all necessary inputs and an area of 2 km² or greater, for which velocity is generally resolved well. We remove those that are known to be surge-type⁴², and glaciers with inverted or erratic elevation change or mass balance profiles (Supplementary Information), indicative of erroneous input data or undocumented surge behavior, which is common for larger glaciers in this region⁴³. The final set of 5527 glaciers represent 71% of the total volume (56% of total area) of glaciers larger than 2 km² in the region.

We present the area-weighted mean profile of SMB relative to each glacier's elevation range for all of HMA in Fig. 1, as well as each glacier's mean SMB resolved from our method. The difference between SMB and thinning patterns (Fig. 1b) strongly underlines the necessity of accounting for ice flow²³, and enables our method to resolve accumulation and ablation areas. By representing density differences in accumulation and ablation areas (Methods), our results reveal a consistent bias of +0.07 m w. e. a⁻¹ in past estimates relying on a single density value (Fig. 1, Supplementary Fig. 11). We calculate more negative mass balances despite our subset of glaciers exhibiting a slight positive bias in terms of volume change, and this bias exceeds the reported uncertainty in many subregions (Table S2).

The subregional SMB profiles (Supplementary Figs. 11, 12) emphasize differences in subregional glacier health. Consistent, distinct mass balance gradients are evident for accumulation and ablation areas. Also apparent for many subregions is a mass balance gradient reversal in the lowest elevations, attributable to inverted SMB gradients in debris-covered areas³²,³⁸ and erroneous input data (Fig. S31). The Nyainqentangla subregion shows the most negative SMB profile, with glaciers accumulating mass over only the upper 20% of their elevation range. The Everest, Spiti Lahaul, and Tien Shan subregions exhibit similar normalized elevation SMB profiles despite occupying very different ranges. Glaciers in these subregions accumulate mass over the upper 40% of their elevation range. In contrast, glaciers in the approximately neutral-balance Kunlun Shan and Karakoram accumulate mass over the upper 60% of their elevation range and exhibit less negative SMB in ablation areas.

**ELAs and accumulation area ratios (AARs).** We leverage the distributed SMBs to resolve ELAs and AARs for each glacier (Fig. 2, “Methods”). The area-weighted ELA for the entire region
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Summary of the altitudinally resolved specific mass balance (SMB) results of this study. **a** Mass balance of glaciers in High Mountain Asia based on our specific mass balance results and compared to that derived from elevation change data (dH/dt)$^{15}$. **b** Regional area-weighted mean SMB and dH/dt for the 2000–2016 period, also indicating mean rates of accumulation and ablation, all shown with elevation normalized to each glacier's elevation range. **c** The spatial pattern of glaciers analyzed and our results for glacier and subregional mass balances (SMB), also compared to the dH/dt results for each subregion and the glacier subset that we process. Uncertainty and subregional profiles are shown in Supplementary Figs. 8, 9. Glaciers we did not process are shown in white, and the background is hillshade of the GTOPO30 dataset sourced from the USGS (https://doi.org/10.5066/F7DF6PQS). **d** Geographic position of the HMA domain.</figcaption></figure>

is 5283 m a.s.l., and the regional AAR is 0.51 (Table S3). Our
glacier-specific ELA results extend the local perspective of pre-
vious studies to the entirety of the region. Studies of individual
glaciers in the Tibetan Plateau and in the Central Himalaya have
shown ablation areas extending to 6000 m a.s.l.³³,⁴⁴, but our
results show that this is true for many glaciers on both sides of the
Himalayan Arc and throughout western Tibet (Fig. 2). This is not
due to a bias in our glacier-specific ELA values, as they agree with
the few reported ELAs based on seasonal snowline elevations and
debris extent (Figs. S83–85). Crucially, our ELA results provide a
unique dataset that can enable a novel calibration and validation
of glacier models¹²,⁴⁵.

ELAs and AARs vary considerably between glaciers across the
region, with standard deviations of 678 m and 0.32, respectively.
ELAs follow broad topographic variations, as glacier extent is
limited by the intersection of climate with topographic avail-
ability. However, we find that the median glacier elevation (a
commonly used proxy for ELA) and ELA clearly differ (median
absolute deviation of 193 m, Fig. S16), emphasizing the
importance of accounting for each glacier's geographic and
climatic setting¹⁰. In the Nyainqentanglha subregion, for
example, many glaciers exhibit low ELAs (Fig. 2a) despite losing
mass rapidly¹⁴,¹⁵,²⁰; these maritime glaciers are sustained by high
annual precipitation⁴⁶, but are highly sensitive to warming⁴⁵.

The AAR indicates the portion of glacier area gaining mass at the surface, and thus reflects the glacier's relative health within its local context. AARs are typically in the range of 0.5 to 0.8 for mountain glaciers roughly in climatic-geometric balance, and are lower for glaciers losing mass³². Numerous glaciers in the Karakoram and Kunlun Shan had large accumulation areas (AAR > 0.5) in the early 21st century, a clear manifestation of the

'Karakoram Anomaly'⁴⁷. Our results show that some glaciers in the neighboring Pamir, Pamir Alai, and Tibetan Plateau subregions also have high AARs (Fig. 2b). Across HMA, 40% of glaciers have AAR > 0.5, but such glaciers are very uncommon beyond the area of the Karakoram Anomaly (Fig. S17). Strikingly, the influence of the Karakoram Anomaly is not discernible in the spatial distribution of ELAs; neither the Karakoram nor the Kunlun Shan shows lower ELAs than adjacent subregions (Table S3). The smooth variations of ELA but the abrupt change of AAR around Karakoram Anomaly glaciers suggests that topographic factors might contribute to the currently stable regional mass balance. That is, the glaciers within this zone may be exceptional in part because there is an extensive high-elevation area available for them to accumulate snow, unlike in other regions (Fig. S19)⁴⁸,⁴⁹. Consequently, recent increases in high-altitude precipitation⁵⁰ would affect a disproportionately large glacier area in this subregion.

Contrasting with the high AARs of the Karakoram Anomaly
glaciers, for 16% of studied glaciers (10% of area) no accumula-
tion area exists and annual net loss occurs at the surface across all
elevations (Fig. 2b). 32% of glaciers have very small accumulation
areas (AAR < 0.1), accounting for 19% of glacier area, while 41%
(23% of area) have AAR < 0.2. Smaller glaciers exhibit lower
AARs in our results (Fig. S18), and some of these glaciers may be
cases where the observed velocity is erroneously low. However,
our comparison to reference measurements (Supplementary
Information) highlights numerous small, nearly stagnant glaciers
where the observed surface lowering corresponds directly to SMB
measurements, suggesting that mass replenishment due to ice
flow is negligible. Considering only the 1982 glaciers larger than
5 km², which are more likely to have measurable surface
<figure><img src="image_5.png" /><figcaption>**Fig. 2** The state of studied glaciers in High Mountain Asia in 2000–2016.
a The spatial pattern of equilibrium line altitudes (ELAs) determined in this study. White dots depict the full distribution of glaciers in the region. Black lines depict the outlines of subregions following¹⁴. b The region-wide histogram of ELA values, also indicating median (black) and area-weighted mean (blue) values. c The spatial pattern of accumulation area ratio (AAR) values determined in this study. d The region-wide histogram of AAR values, also indicating median (black) and area-weighted mean (blue) values. Background is the August 2004 NASA Blue Marble⁷⁵.</figcaption></figure>

velocity¹⁶, we still find that 7.5% of glaciers have an AAR < 0.1. Low-AAR glaciers are most concentrated in Eastern Nyainqentangla, which is the subregion with the highest rates of mass loss. Here, 50% of glaciers have an AAR less than 0.2 (Supplementary Table 3, Supplementary Fig. 18). Such “headless” glaciers are also surprisingly common across the rest of the Himalaya, Tibetan Plateau, and Tien Shan (40% of glaciers in these regions); are less frequent in the Pamirs, Hindu Kush, and Karakoram (30%); and are rare in the Kunlun Shan (8%).

Our AAR results depict a picture of strong imbalance for most glaciers in HMA. The very low AAR values for much of HMA suggest that accumulation areas have been substantially reduced by increasing ELAs. However, we note that AARs can be depressed to 0.3–0.5 for heavily debris-covered glaciers sustained by avalanches³², which are common in parts of HMA. Thus, although the AAR pattern highlights the contrast between glaciers affected by the Karakoram Anomaly and the rest of the region, this metric may be biased in some areas. To assess glacier health in an unbiased manner, we determine an additional indicator of glacier health: the ablation balance ratio.

**Sustainable ablation in major basins.** We use our SMB results to assess the origin of glacier ablation in the major basins draining HMA by partitioning the total annual glacier ablation into “imbalance” and “balance” components (Fig. 3, Supplementary Table 4). The “balance” component is the glacier ablation that is compensated on an annual basis by accumulation (Methods), and which we consider to be sustainable in early 21st-century conditions¹. Crucially, the SMB results allow us to determine these values directly for each glacier, whereas prior available

estimates were obtained only at the basin scale by comparing observed thinning with glacier models¹,²⁰. Our results indicate that 40 ± 11% of ablation from HMA glaciers is unsustainable in early 21st century conditions. Although this indicates that the glacier contribution to streamflows will eventually reduce, there is considerable variability between regions and individual glaciers.

Basins fed by the Karakoram Anomaly glaciers (Indus, Amu Darya, Syr Darya, Tarim Interior) are the most important water towers globally³. The prevalence of surging glaciers in these basins leads to a relatively smaller sample size in our study (Supplementary Table 4), and these basins exhibit higher uncertainty in the absolute mass change of glaciers affected by the Karakoram Anomaly¹⁵,²⁰,⁴⁷. Nonetheless, we show that these basins’ glaciers were much healthier compared to the rest of HMA for our study period, with over 50% of annual glacier ablation balanced by accumulation and numerous individual glaciers exceeding 100% balanced ablation (Fig. 3). The Indus basin is of particular interest due to the high dependence of its downstream populations on snow and ice melt, especially in drought conditions¹,²; its high vulnerability to future societal and environmental change²,³; and the high ice-melt content of its summertime streamflows¹. For this key basin, we find that 65% ± 23% of glacier ablation was balanced by accumulation during 2000–2016 (Supplementary Table 4).

Contrasting with basins supported by the Karakoram Anomaly glaciers, nearly all other basins’ supply of glacier ablation is primarily imbalance ablation. Among these, the Ganges–Brahmaputra basin stands out as a vulnerable basin with important water supplies (Supplementary Fig. 21). Although the Ganges and Brahmaputra rivers are sustained by strong monsoonal precipitation, glaciers provide the majority of streamflow in drought years¹. We show that the majority of ablation from these glaciers is imbalance ablation (the balance ratio is 48% ± 9%, Fig. 3, Supplementary Table 4). The future decline in glacier meltwater supply in this region may seem relatively minor on an annual basis due to precipitation excess in normal years¹² but is crucial seasonally, affecting the growth of cash crops in the water scarce pre-monsoon⁵¹. Given the likelihood of considerable cryospheric and environmental change, the societal pathway of adaptation to change in these basins will directly control downstream communities’ resilience to water resource change⁶,⁵².

A previous study¹ determined a regional balance ratio of 38% for drought years; our results indicate that a small balance ratio (40%) is the recent-period norm, rather than the exception. Our results differ at the subregional level: we find a 32% higher balance ratio for the Ganges and Brahmaputra than ref. ¹, but 6–18% lower balance ratios for the basins fed by the Karakoram Anomaly (Supplementary Table 5). Our basin ablation ratios are less than those of ref. ²⁰ by 15–25% for all basins except the Tarim, where the sign of mass balance is uncertain¹⁵,²⁰. Although our imbalance ablation results are slightly higher due to our representation of density, our total ablation estimates are much lower (Supplementary Table 5). This is likely due to compensation of melt and accumulation errors in the model used by that study, which was overparameterized and did not represent the regionally important influence of debris cover¹². Constrained by thinning rates alone, the model would overestimate total melt and compensate this error by exaggerating accumulation, leading to higher balance ratios.

The balance ratios for the Amu Darya, Tarim Interior, and Indus basins are particularly affected by the uncertainty and heterogeneity of subregional volume change, demonstrating the need for systematic measurements in this region⁵³, particularly of accumulation rates³⁴. The sustainability of river discharge in these important and vulnerable basins is dependent on the anomalous health of glaciers in the Kunlun Shan and Karakoram ranges⁴⁷. Despite the higher balance ratios in these basins, a
<figure><img src="image_6.png" /><figcaption>**Fig. 3** The quantity and context of glacier ablation for principal drainage basins in High Mountain Asia. Analyzed glaciers are colored according to the portion of total annual ablation that is compensated by accumulation, which is greater than 100% for glaciers gaining mass. The portion of balance ablation derived from our results is shown for major river basins, indicating uncertainty with the dashed lines, and scaled by area according to the total estimated glacier ablation within each ("Methods", Supplementary Tables 4, 5). Basin vulnerability is colored according to the global range³. Background is a hillshade of the GTOPO30 dataset sourced from the USGS (https://doi.org/10.5066/F7DF6PQS).</figcaption></figure>

considerable portion of glaciers are losing mass. Our results show that the major tributaries of the Indus are supplied by glaciers in contrasting health: the Chenab and Satluj are supplied by imbalance ablation from unhealthy glaciers in the Western Himalayas, while the Indus itself is supplied by imbalance ablation from the Ladakh Range and balance ablation from the Karakoram. Combining these distinct signals, our results indicate that despite the near-neutral mass balance in the Karakoram, 19% ± 12% of subregional ablation was imbalance ablation in the early 21st century (Supplementary Table 3).

**Implied glacier and ablation changes.** The considerable imbalance ablation for 2000-2016 indicates climatic-geometric disequilibrium across the region. We determine how glaciers would respond to this disequilibrium, if maintained, and find that early 21st century mass balance regimes imply a change of -23% ± 1% of glacier volume in HMA by 2100 (Fig. 4, Methods). All subregions along the Himalaya lose at least 35% of their present-day volume by 2100, contrasting with volume reductions of 10-20% for the Karakoram, Pamir Alai, Pamir, and Hindu Kush and a volume gain of 2.1% for the Kunlun Shan (Supplementary Table 4). Our results indicate that 25% of glaciers across the region would lose at least 50% of their current ice volume by 2100 without any warming (Supplementary Figs. 22, 23). We calculate a more negative long-term volume change of -34% ± 2% by 2200 than the 27-29% committed loss estimate based on field observations²⁵. By resolving the implied volume change for a large population of glaciers individually, we mitigate against the biases of sparse glaciological measurements⁵⁴. Our estimate is only slightly lower than recent projections of 29% ± 12%¹² to 36% ± 7%⁵ mass loss by 2100 under the RCP2.6 climate scenario,

suggesting that under that climate pathway most glacier loss is already committed by current climatic-geometric disequilibrium.

Associated with the regional losses in glacier volume, we find a regional change of -28% ± 6% in total annual ablation rates by 2100 (Fig. 4). Subregions experience variable changes in glacier ablation largely following the changes in glacier volume, but the ablation changes are stronger than volume loss in the Karakoram and Pamir subregions. As with the volume change estimates, the projected general reduction in glacier ablation should be taken as a baseline estimate of change. Although continued 21st-century warming will lead to a peak in glacier meltwater supply⁵⁵, this will exacerbate glaciers' climatic-geometric imbalance and lead to more severe eventual reductions in glacier ablation.

We note that the mass losses are partly obscured by the Karakoram Anomaly: other than the globally unique mass gains in the Kunlun Shan and parts of the Karakoram, the region has an implied volume loss of 31% by 2100. Considering recent and further climate warming, our results represent minimum estimates of future volume loss; sustained warming would be likely to overcome recent increases in snow accumulation in the Karakoram and Kunlun Shan⁴⁷,⁵⁰, exacerbating regional glacier loss. Current projections for the Karakoram and Kunlun Shan show ice losses of 10-35% by 2100 in response to continued but reduced emissions under RCP2.6, and substantial ice losses of 30-60% for RCP4.5⁵,⁶,¹². Consequently, disentangling the causes of the Karakoram Anomaly⁴⁷ and understanding its resilience to 21st-century warming remain key priorities for scientists and stakeholders alike.

**Implications for glacier modeling and monitoring.** Our baseline estimate of 23% ± 1% volume loss by 2100 will be exceeded given that most climate trajectories indicate continued
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Changes in glacier volume and total ablation by 2100 implied by 2000-2016 mass balance regime, based on a glacier retreat and advance flow parameterization. Regional icons depict the portion of glaciers with positive (blue) or negative (red) glacier-wide mass balance, implied volume change by 2100 as portion of glacier volume (black), and implied change in total annual glacier ablation by 2100 as a portion of current total annual ablation (orange). Background is a hillshade of the GTOPO30 dataset sourced from the USGS (https://doi.org/10.5066/F7DF6PQS).</figcaption></figure>

warming⁵,¹² and the progressive deglaciation of this region will lead to a cascade of changes to ecosystems and society⁶,⁵². We advocate for the improvement of dynamic glacier models (e.g.²⁷) to better reproduce the long-term mass balance of glaciers in HMA by including key unrepresented processes¹² such as localized mass accumulation due to avalanches, reversed mass balance gradients due to supraglacial debris, and frontal ablation due to ice-marginal lakes¹⁷. Spatially explicit glacier models should not be calibrated to glacier thinning datasets alone, which leads to the compensation of SMB and flux divergence errors, and can lead to errors in both melt and accumulation. Instead, models should be constrained with both glacier thinning and surface velocity observations. Our results extend the sparse glaciological measurements in HMA and thus provide the opportunity for novel strategies¹² to calibrate mass balance models directly to each glacier's altitudinally resolved SMB, or to regionally resolved mass balance gradients (Supplementary Information).

Our method has the potential to generate very novel datasets and understanding of SMB patterns worldwide, but it also demonstrates the need for improvements to existing observations. Glacier basal condition and ice rheology are poorly known at all but a few study sites. Novel field measurements of these properties would enable the uncertainties around SMB determined through a continuity approach to be significantly reduced. Robust assessments of elevation change rates are now possible at the regional scale²⁰, but the established average density value²⁴ should be reconsidered for glaciers with small accumulation areas. Problems in the input datasets of velocity and ice thickness forced us to discard results for 25% of the 7341 glaciers analyzed, and prevented application to glaciers smaller than 2 km². Ice thickness is generally the most uncertain input dataset, and additional ice thickness measurements are needed across HMA to constrain ice thickness models⁵⁶, especially for the region's debris-covered areas. New analyses of modern high repetition, high-resolution satellite data are likely to resolve flow patterns in

problematic areas such as small glaciers, tributary junctions, and icefalls⁵⁷.

Despite these challenges, we have resolved multidecadal SMB profiles across HMA, providing detail of the region's heterogeneous glacier health. We show that imbalance ablation, not replenished by annual snow accumulation, dominates the contribution of glaciers into most river basins, with the exception of basins fed by the Karakoram Anomaly glaciers. 35% of glaciers across the region are very unhealthy and are expected to lose at least half their volume by 2100 without additional climate warming. Our results provide a novel, spatially extensive dataset to calibrate and validate a new generation of glacier models capable of resolving glacier mass balance and ice dynamics at high temporal and spatial resolution. This approach paves the way to resolve the SMB across glaciers globally and for multitemporal periods to characterize the trajectory of glacier change.

## Methods

**Continuity approach to mass balance calculation.** Our mass balance reconstruction approach solves the continuity equation (Eq. 1, Supplementary Fig. 1). For any area of the glacier, in Eq. (1), $dH/dt$ is the annual rate of elevation change at the glacier surface, $\dot{b}$ is the annual SMB (surface, internal, and basal mass balance combined, with frontal ablation if relevant) of that area, and $V*\mathbf{q}$ the annual flux divergence (determined below), accounting for the density $\rho$ for each quantity²²,²³,⁵⁸. We aim to calculate $\dot{b}$, which, assuming that frontal ablation and the englacial and subglacial components are negligible, is the surface mass balance. Crucially, however, this does not equate to the glacier's melt, as it is an integrated signal for each pixel, which may experience seasonal accumulation and ablation.

$$ \frac{\rho_{\mathrm{d}H}}{\rho_{\mathrm{H}_2\mathrm{O}}} \frac{\mathrm{d}H}{\mathrm{d}t} = \dot{b} - \frac{\rho_{\mathrm{v}q}}{\rho_{\mathrm{H}_2\mathrm{O}}} \nabla * \mathbf{q} \quad (1) $$

We apply the continuity equation on a fully distributed basis. For this, we use ASTER-based 2000-2016 annual surface lowering trends¹⁵, ITS_LIVE HMA ice surface velocity products¹⁶,⁵⁹ and multi-model consensus ice thickness estimates³⁶ which correspond to the Randolph Glacier Inventory version 6.0 outlines⁴². These datasets are available in different projections and spatial resolutions. For each glacier, we define a grid for our analysis using the local projection used by³⁶, and
we vary the grid resolution based on the size of the glacier: 50 m for small glaciers (<15 km²), 100 m for larger glaciers (up to 80 km²), and 200 m for very large glaciers (>80 km²). We reproject and resample the surface lowering data (provided at 30 m resolution) and its stated uncertainty to this grid using a cubic spline. We used cubic splines to reproject both ends of the surface velocity vectors to preserve true orientation before resampling these data and their stated uncertainty from their 240 m resolution. Finally, we degrade the corresponding ice thickness data (provided at 25 m resolution); possible concerns of circular analysis with this dataset are mitigated by the method's performance for debris-covered areas and with tests using distinct individual ice thickness models (Supplementary Information). To maintain a continuous dataset over each glacier, we do not filter the surface velocity and surface lowering datasets before reprojection, but instead assess the uncertainty through our calculations.

**Ice flux and flux divergence.** We calculate the ice flux vector **q** at each cell according to Eq. (2), where *h* is the ice thickness (m) and **u** is the annual ice surface velocity vector (m a⁻¹).

$$ \mathbf{q} = h\mathbf{y}\mathbf{u}_s \qquad (2) $$

$\mathbf{y}\mathbf{u}_s$ represents the column-average ice velocity, with the constant $\gamma$ representing the relative importance of basal motion and vertical ice shear deformation (Supplementary Fig. 1). We model $\gamma$ for each glacier individually. For this, we use a Monte Carlo analysis to estimate the depth-integrated velocity at a point assuming simple shear with an assumed ratio of ice motion attributable to basal sliding $|\frac{\mathbf{u}_s}{\mathbf{u}_s}|$ and a thickness estimate. For a given ice thickness and basal sliding ratio, we calculate the velocity at each depth following⁶⁰, then determine $\mathbf{y}\mathbf{u}_s$. We perform this calculation for 10,000 sets of randomly drawn values of ice thickness, flow rate factor $n$, and basal sliding for each glacier. For the ice thickness distributions we use the distribution of ice thickness values produced for that glacier by³⁶. For $n$ we note that $n=3$ is appropriate for many glacier modeling situations³⁸ and use a Gaussian distribution with ($\mu=3$, $\sigma=0.067$). For the portion of flow attributable to sliding, this is dependent on both ice rheology and basal state. Neither basal sliding nor ice internal thermal profiles are well constrained for glaciers in HMA, but authors have variously assumed or determined temperate, cold, and polythermal glaciers across the region⁶¹⁻⁶⁶, demonstrating variable thermal regimes and basal conditions across High Asia. We acknowledge that (1) many small, high-altitude glaciers are likely to be cold-bedded⁶⁴, but (2) there is increasing evidence that the lower-elevation tongues are polythermal with temperate beds⁶⁵,⁶⁶. In addition, although there are many large proglacial lakes in the region which are known to affect terminus ice velocities¹⁷,⁶⁷, it is not likely that an extensive portion of glacier ice approaches flotation. Nonetheless, without widespread knowledge of the importance of basal sliding across HMA, we assume a uniform distribution across [0,1] for our basal sliding factors. In addition to providing an estimate of $\gamma$, this Monte Carlo approach allows us to estimate its uncertainty, $\sigma_{\gamma}$.

The flux divergence $\nabla * \mathbf{q}$ represents the vertical component of ice velocity at the glacier surface, which leads to submergence in areas of divergent flow and emergence in areas of convergent flow. We calculate $\nabla * \mathbf{q}$ on a pixel basis using a centered-difference scheme based on the divergence of $\mathbf{q}$ (Eq. 3, Supplementary Fig. 1).

$$ \nabla * \mathbf{q} = \frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y} \qquad (3) $$

**Density correction.** Our continuity approach assumes the mean ice density within the domain does not change with respect to time. This is generally reasonable in the ablation area or over a period that densification processes can be considered constant, leading to uncertainties on the order of 2%⁵⁸. However, the density of snow, firn, and ice at the glacier surface must still be accounted for in order to derive $\dot{b}$. Geodetic studies often use a single value of 850 kg m⁻³ or zonal values for accumulation and ablation areas²⁴.

We first assume that all ice fluxes are composed purely of glacier ice, such that our flux divergence has a density of 900 kg m⁻³. To determine the effective density of our elevation change signal, we consider the physical situation corresponding to the particular values of elevation change and flux divergence (Table S6). Where elevation change and flux divergence both have positive signs, we interpret mass accumulation as occurring and we assign a density of 600 kg m⁻³. If both are negative, we assume this corresponds to ablation of glacier ice with a density of 900 kg m⁻³. There is ambiguity about the state of glacier ice where flux divergence and elevation change are aligned, but this most likely corresponds to dH/dt and SMB values close to zero, and we choose an intermediate density of 850 kg m⁻³ to represent the variable likelihood of elevation change being composed of glacier ice or wetted snow and firn. We assume that the density uncertainty is approximately 60 kg m⁻³ for all values²⁴.

**Uncertainty.** The uncertainty in the simplified continuity equation (Eq. 1) assuming independent errors is given analytically by

$$ \sigma_b = \sqrt{ \left[ \left( \frac{\sigma_{vq}}{\nabla \mathbf{q}} \right)^2 + \left( \frac{\sigma_{pvq}}{\nabla v q} \right)^2 \right] \left( \rho_{vq} \nabla * \mathbf{q} \right)^2 + \left[ \left( \frac{\sigma_{dH}}{dH} \right)^2 + \left( \frac{\sigma_{\rho H}}{\rho dH} \right)^2 \right] \left( \rho_{dH} dH \right)^2 } \qquad (4) $$

In this equation, the flux divergence uncertainty for an individual pixel integrates the uncertainty for each of four fluxes dependent on multiple inputs (Eq. 2, Supplementary Fig. 1). Assuming these inputs to be subjected to completely random error would lead to an unrealistically high uncertainty estimate; given the ice thickness uncertainties, this would effectively assess the change in flux divergence due to a 40% change in ice thickness between adjacent pixels. Instead, we consider the uncertainty of each input dataset in terms of systematic bias and random error at the scale of an individual pixel and its neighbors. We assume that the uncertainties of the input datasets are not correlated to one another, and consider systematic and random errors separately for each.

We therefore derive the normalized ice thickness uncertainty $\sigma_b^N$ for each glacier as the standard error between individual ice thickness estimates on a pixel-by-pixel basis provided by³⁶, which we normalize relative to the consensus thickness estimate. We take the 68th centile value from the empirical distribution of normalized thickness standard errors (ie, 68% of standard errors are below this value; this is equivalent to the standard deviation for a one-sided distribution) as indicative of the glacier-wide systematic ice thickness uncertainty. We additionally consider that ice thickness is likely to have a random error component that the modeled ice thickness datasets do not reproduce, which we estimate to follow a gaussian distribution with ($\mu=0$ m, $\sigma=10$ m).

We use the pixel-wise ITS_LIVE reported error to derive the systematic normalized surface velocity uncertainty $|\frac{\sigma_u}{u}|$ for each glacier. Specifically, we use the reported error in surface velocity magnitude to determine the 68th centile value of $|\frac{\sigma_u}{u}|$ for each glacier, which we consider the systematic uncertainty. There is also a component of random error in the velocity data, but we assume that the random error is negligible at the scale of adjacent pixels in our analysis. We justify this assumption based on two factors. First, the velocity product is a synthesis of multiple years of observations and our target glaciers are non-surging mountain glaciers, which display consistent spatial patterns of velocity. We therefore expect that the flow direction and relative magnitude are generally very accurate, but that the multi-year mean speed is uncertain due to velocity change over the period of analysis and date biases in the data synthesis¹⁶; this is reflected by our systematic error. Second, we note that the x- and y-displacement uncertainty may be random at the scale of the velocity product (240 m) but is not likely to be random at the scale of adjacent pixels in our analysis. Our assumption is that pixel-scale patterns of ice velocity change accurately reflect larger-scale patterns of ice dynamics that are captured by the velocity data.

For the flux calculation, we assess the random error $\sigma_{\gamma}$ as the standard deviation of calculated $\gamma$ values from the 10,000 run Monte Carlo analyses for each glacier, described above. Considering the agreement in dH/dt products from recent studies¹⁵,²⁰, we consider the uncertainty in dH/dt to be limited to random error. We therefore use the reported uncertainty as $\sigma_{dH}$ in our Monte Carlo analysis. Finally, we assume the random uncertainty in density estimates to be 60 kg m⁻³.

We integrate each source of uncertainty by perturbing our input data in a Monte Carlo analysis with 1000 distinct runs for each individual glacier. Using the uncertainties outlined above, we perturb our inputs with (1) random, spatially uncorrelated noise added to the dH/dt data, (2) randomly chosen systematic scaling of the **u** data, (3) random, systematic scaling of the *h* data, (4) random, uncorrelated noise added to the *h* data, (5) random systematic scaling of the $\gamma$ estimates, (6) random variations in the density values $\rho_{dH}$ and $\rho_{vq}$. This enables us to estimate the integrated uncertainty in $\nabla * \mathbf{q}$ and $\dot{b}$. We also use the Monte Carlo results to determine the uncertainty in our derived metrics of ELA, AAR, committed volume loss, and balance ratio as the standard deviation of each metric for the full population of runs.

**Mass balance profiles.** Although our calculations are performed pixel by pixel across each glacier, slight inconsistencies between the observed velocity pattern and modeled ice thickness pattern can lead to an unrealistic pattern in flux divergence. This is due to several factors: (1) systematic decorrelation in the velocity product due to either a lack of identifiable features (particularly in accumulation areas) or rapid ice flow (particularly in icefalls), (2) the necessary use of a shape factor to distribute ice thicknesses across the glacier width⁶⁸,⁶⁹ which can vary from glacier to glacier and even across glaciers, (3) the inability of current ice thickness models to treat glacier tributaries separately³⁶, and (4) the spatial variations of longitudinal stress gradients.

To mitigate this problem and to provide higher-confidence distributions of specific annual mass balance, we segment each glacier into hypometric bins. To remove local undulations at the glacier surface, the ASTER GDEM v3⁷⁰ is resampled to the resolution of our analysis, then smoothed with an 11x11 Gaussian low-pass filter using a 2σ threshold. We segment the resulting DEM into 25 m elevation bins, then intersect the result with a hole-filled version of the debris cover maps provided by³¹. For each segment, we determine the mean values of $\nabla * \mathbf{q}$, dH/dt, and $\dot{b}$, and the uncertainty is assessed for each variable through quadrature of
the distributed estimates. This aggregation step is crucial to reduce the effects of the factors listed above, and to resolve the overall pattern of SMB rather than amplifying noise due to errors in individual datasets. Finally, our SMB results are compared to available surface mass balance measurements from the World Glacier Monitoring Service⁷¹ and other published literature, and to the results of³⁸ (Supplementary Table 1).

Based on the method's performance, we limit our analyses to larger glaciers (>2 km² in area) which are more likely to show a clear velocity signal¹⁶. We also remove surging glaciers from consideration for further processing, which we identify based on the RGI6.0 attributes. We additionally identify glaciers with erratic surface lowering or mass balance patterns, also indicative of surging or lower quality source data. In particular, we limit our glaciers for further analysis to those that satisfy the following conditions: the detrended altitudinal $dH/dt$ profile has a standard deviation of less than 3 m a⁻¹ and the $dH/dt$ profile has a nonnegative correlation with elevation. We consider these characteristics to be indicative of surging behavior. Finally, we only retain glaciers with the following criteria, which we consider to be indicative of higher quality input data and results: the optimized ELA has an Accuracy of at least 0.5, the detrended SMB profile has a standard deviation of less than 3 m w.e. a⁻¹; and the mean SMB uncertainty is less than 3 m w.e. a⁻¹. This leaves a population of 5527 glaciers representing 71% of the total ice volume of RGI regions 13, 14, and 15. Due to the quality controls, subregions are not uniformly sampled (Supplementary Tables 3, 4), which we account for (Regional Results).

**Determination of ELA and AAR.** The ELA is a single elevation contour ideally intended to distinguish between accumulation areas and ablation areas. Given our distributed mass balance dataset, we determine glacier-specific ELAs through an error minimization approach. We first classify pixels as accumulating or ablating mass based on the sign of SMB in our results. We then use each integer elevation within the glacier's elevation range as a binary classifier to produce a segmentation of accumulation and ablation areas. We assess each segmentation relative to our gridded results by determining the confusion matrix and computing its accuracy. We determine the ELA as the elevation that gives the best Dice coefficient⁷² for the segmentation of accumulation and ablation areas (Supplementary Figs. 2–5). For glaciers whose optimal ELA is at either end of the glacier's elevation range (indicating mass loss or gain at all elevations), we fit a linear trend to the SMB and extrapolate significant trends to determine the elevation with SMB=0, which we take as an indicator of the theoretical climatic ELA for the glacier's location. For all glaciers for which we successfully resolve an ELA, we then determine the AAR as the portion of glacier area that lies above the ELA. We compare our ELA results to available datasets in the Supplementary Information.

**Calculation of ablation balance ratio.** Following¹, we calculate the balance portion of ablation (corresponding to the ratio of balance ablation to total ablation) for each glacier, subregion, and river basin. For each glacier, we calculate the total ablation directly based on our distributed SMB results, summing all pixels with a negative SMB. We then calculate the imbalance ablation for each glacier as the specific annual mass balance multiplied by glacier area¹⁵. From the total and imbalance ablation rates, we determine the rate of balance ablation (Supplementary Fig. 2). We express this for each glacier as a ratio of the balance ablation to total ablation, expressed as a percentage (Fig. 3); glaciers experiencing net annual accumulation thus have a balance ratio greater than 100%.

**Calculation of implied volume change.** To assess the volume change implied by our mass balance profiles, we developed a parameterization of glacier retreat and advance similar to²⁸,⁷³. In this framework, the annual mass balance is calculated based on our SMB results, and a $\Delta h$ parameterization is used to redistribute this mass loss or gain across the glacier, updating the ice thickness of the glacier. The SMB dataset only changes based on glacier extent changes, eventually leading to an equilibrium state after numerous iterations. This parameterization approach has been demonstrated to appropriately represent glacier retreat by implicitly representing ice dynamics⁷³.

For each glacier with a clear signal of mass loss (mean mass balance less than -0.1 m w.e. a⁻¹), we develop a $\Delta h$ parameterization based on the thinning rates from¹⁵. For glaciers with ambiguous thinning patterns, we use the $\Delta h$ parameterization from⁷³ directly. For glaciers with a positive mass balance, we found that the 5 m a⁻¹ thickening threshold for advance used by²⁸ did not allow HMA glaciers to advance to a steady state. We therefore instead allow glaciers to advance when the terminus longitudinal gradient exceeds 10 degrees. We determine this longitudinal gradient based on the mean thickness of the lowest $N_t$ on-glacier pixels, where $N_t$ is the number of pixels equating one glacier width in the terminus area, and the size of each pixel. This longitudinal gradient threshold was chosen such that the effective volume-area scaling relationship noted by⁷⁴ holds for our advancing glaciers. If a glacier is allowed to advance, the lowest-elevation $N_t$ for glacier-marginal pixels become appended to the glacier, and are thickened by the prior terminus height, but the advancing fraction is limited to 50% of the glacier's total volume gain. The non-advancing mass accumulation is distributed altitudinally according to the original $\Delta h$ parameterizations²⁸,⁷³. For advancing glaciers, we extrapolate the SMB from the glacier terminus at a rate of 0.07 m w.e. m⁻¹, which is the median observed ablation gradient in our extended database of field measurements (Supplementary Table 1).

We carry out this $\Delta h$ parameterization for our subset of 5527 glaciers, updating ice thickness, glacier extent, SMB, and elevation datasets with an annual timestep. Although our simulations are performed for 200 years, we highlight the volume change results for the year 2100 as the $\Delta h$ parameterization is most robust for multidecadal periods⁷³. Finally, we note that these results depend strongly on the uncertainty of our SMB results, so we run 30 simulations for each glacier, varying the SMB systematically by the $dH/dt$ uncertainty, which is the primary source of uncertainty for glacier-wide mass balance. We therefore report the mean and standard deviation of regional glacier outcomes for 2100 implied by the current mass balance regimes and their uncertainty.

**Regional results.** We perform the above calculations for each glacier within our regional subset. We then aggregate these values to distinct subregions as defined by¹⁴⁻¹⁶, and to major river basins³ to provide a larger-scale perspective on the heterogeneity of glacier health. For ELA and AAR, we determine the area-weighted mean and its uncertainty, as well as the median value within each zone. For ablation balance and implied volume change, we aggregate total ice volumes but do not assume random uncertainty, instead determining the mean normalized uncertainty (Supplementary Tables 3, 4).

For the river basins we also seek to estimate total ablation including glaciers not represented in our regional subset due to their small size or low-quality input data. We therefore determine the total imbalance ablation in each basin from the results of¹⁵ with our basin outlines and correcting for the subregional mass balance biases due to density estimates (Fig. 1). We then used the ratio of balance to imbalance ablation from our subset of glaciers for the basin to estimate the total ablation in the river basin. This leads to a different regional mean ablation balance ratio (60%) than for our subset (50%) by accounting for subregional sampling bias. As the uncertainty of the scaling is the major source of uncertainty for the regional balance ratios (Fig. 3), we do not scale the results of our future simulations, but report the regional aggregated values.

## Data availability

The SMB datasets generated and analyzed during the current study are available in the Zenodo repository, https://doi.org/10.5281/zenodo.3843292. The elevation change data of ref. ¹⁵ are available at https://doi.org/10.1594/PANGAEA.876545. The mean surface velocity data of ref. ¹⁶ are provided by the NASA MEaSUREs ITS_LIVE project and available at https://its-live.jpl.nasa.gov/. The consensus ice thickness dataset of ref. ³⁶ is available at https://doi.org/10.3929/ethz-b-000315707. The glacier outlines of ref. ⁴² are available at https://www.glims.org/RGI/rgi60_dl.html. The supraglacial debris extents of ref. ³¹ are available at https://doi.org/10.5880/GFZ.3.3.2018.005. The WGMS Fluctuations of Glaciers database is available at https://doi.org/10.5904/wgms-fog-2019-12. River basin boundaries used in this study are available at http://www.fao.org/nr/water/aquamaps/. The Global Lakes and Wetlands Database is available at https://www.worldwildlife.org/pages/global-lakes-and-wetlands-database.

## Code availability

The code used to produce specific mass balance across High Mountain Asia and to derive implied volume loss is available on GitHub at https://github.com/miles916/HMA continuity.

Received: 18 December 2020; Accepted: 7 April 2021;

Published online: 17 May 2021

## References

1. Pritchard, H. D. Asia's shrinking glaciers protect large populations from drought stress. *Nature* **569**, 649–654 (2019).

2. Viviroli, D., Kummu, M., Meybeck, M., Kallio, M. & Wada, Y. Increasing dependence of lowland populations on mountain water resources. *Nat. Sustain.* **3**, 917–928 (2020).

3. Immerzeel, W. W. W. et al. Importance and vulnerability of the world's water towers. *Nature* **577**, 364–369 (2020).

4. Bolch, T. et al. The state and fate of Himalayan glaciers. *Science* **306**, 310–314 (2012).

5. 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* **549**, 257–260 (2017).

6. Hock, R. et al. High mountain areas. in: ipcc special report on the ocean and cryosphere in a changing climate. in *PCC Special Report on the Ocean and Cryosphere in a Changing Climate* (eds Pörtner, H.-O. et al.) https://www.ipcc.ch/srocc/ (2019).

7. Azam, M. F. et al. Review of the status and mass changes of Himalayan-Karakoram glaciers. *J. Glaciol.* **64**, 61–74 (2018).

8. Bolch, T. et al. in *The Hindu Kush Himalaya Assessment* (eds Wester, P., Mishra, A., Mukherji, A. & Shrestha, A. B.) 209–255 (Springer International Publishing, 2019). https://doi.org/10.1007/978-3-319-92288-1
9. Immerzeel, W. W., Wanders, N., Lutz, A. F., Shea, J. M. & Bierkens, M. F. P. Reconciling high-altitude precipitation in the upper Indus basin with glacier mass balances and runoff. *Hydrol. Earth Syst. Sci.* **19**, 4673–4687 (2015).

10. Sakai, A. et al. Climate regime of Asian glaciers revealed by GAMDAM glacier inventory. *Cryosphere* **9**, 865–880 (2015).

11. Wang, Q., Yi, S. & Sun, W. Precipitation-driven glacier changes in the Pamir and Hindu Kush mountains. *Geophys. Res. Lett.* **44**, 2817–2824 (2017).

12. Rounce, D. R., Hock, R., Shean, D. & Khurana, T. Quantifying glacier mass change in High Mountain Asia through 2100 using the open-source Python *Glacier Evolution Model* (PyGEM). *Front. Earth Sci.* **7**, 1–20 (2020).

13. Marzeion, B. et al. Partitioning the uncertainty of ensemble projections of global glacier mass change. *Earth's Futur.* **8**, 1–25 (2020).

14. Kääb, A., Treichler, D., Nuth, C. & Berthier, E. Brief Communication: contending estimates of 2003-2008 glacier mass balance over the Pamir-Karakoram-Himalaya. *Cryosphere* **9**, 557–564 (2015).

15. 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.* **10**, 668–673 (2017).

16. Dehecq, A. et al. Twenty-first century glacier slowdown driven by mass loss in High Mountain Asia. *Nat. Geosci.* **12**, 22–27 (2019).

17. King, O., Bhattacharya, A., Bhambri, R. & Bolch, T. Glacial lakes exacerbate Himalayan glacier mass loss. *Sci. Rep.* **9**, 18145 (2019).

18. Maurer, J. M., Schaefer, J. M., Rupper, S. & Corley, A. Acceleration of ice loss across the Himalayas over the past 40 years. *Int. Adv.* **5**, 1–12 (2019).

19. Wouters, B., Gardner, A. S. & Moholdt, G. Global glacier mass loss during the GRACE satellite mission (2002-2016). *Front. Earth Sci.* **7**, 1–11 (2019).

20. Shean, D. E. et al. A systematic, regional assessment of High-Mountain Asia glacier mass balance. *Front. Earth Sci.* **7**, 1–19 (2020). in Press.

21. Davaze, L., Rabatel, A., Dufour, A., Hugonnet, R. & Arnaud, Y. Region-wide annual glacier surface mass balance for the European Alps From 2000 to 2016. *Front. Earth Sci.* **8**, 1–14 (2020).

22. Hubbard, A. et al. Glacier mass-balance determination by remote sensing and high-resolution modelling. *J. Glaciol.* **46**, 491–498 (2000).

23. Berthier, E. & Vincent, C. Relative contribution of surface mass-balance and ice-flux changes to the accelerated thinning of Mer de Glace, French Alps, over 1979-2008. *J. Glaciol.* **58**, 501–512 (2012).

24. Huss, M. Density assumptions for converting geodetic glacier volume change to mass change. *Cryosphere* **7**, 877–887 (2013).

25. Zemp, M. et al. Historically unprecedented global glacier decline in the early 21st century. *J. Glaciol.* **61**, 745–762 (2015).

26. Radić, V. et al. Regional and global projections of twenty-first century glacier mass changes in response to climate scenarios from global climate models. *Clim. Dyn.* **42**, 37–58 (2014).

27. Zekollari, H., Huss, M. & Farinotti, D. Modelling the future evolution of glaciers in the European Alps under the EURO-CORDEX RCM ensemble. *Cryosphere* **13**, 1125–1146 (2019).

28. Huss, M. & Hock, R. A new model for global glacier change and sea-level rise. *Front. Earth Sci.* **3**, 1–22 (2015).

29. Brun, F. et al. Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia. *J. Geophys. Res. Earth Surf.* **124**, 1331–1345 (2019).

30. Wijngaard, R. R. et al. Modeling the response of the langtang glacier and the hintereisferner to a changing climate since the little ice age. *Front. Earth Sci.* **7**, 1–24 (2019).

31. Scherler, D., Wulf, H. & Gorelick, N. Global assessment of supraglacial debris-cover extents. *Geophys. Res. Lett.* **45**, 11,798–11,805 (2018).

32. Benn, D. I. & Lehmkuhl, F. Mass balance and equilibrium-line altitudes of glaciers in high-mountain environments. *Quat. Int.* **65–66**, 15–29 (2000).

33. Sunako, S., Fujita, K., Sakai, A. & Kayastha, R. B. Mass balance of Trambau Glacier, Rolwaling region, Nepal Himalaya: In-situ observations, long-term reconstruction and mass-balance sensitivity. *J. Glaciol.* **65**, 605–616 (2019).

34. Lambrecht, A., Mayer, C., Bohleber, P. & Aizen, V. High altitude accumulation and preserved climate information in the western Pamir, observations from the Fedchenko Glacier accumulation basin. *J. Glaciol.* **66**, 219–230 (2019).

35. Mayer, C. et al. Accumulation studies at a high elevation glacier site in Central Karakoram. *Adv. Meteorol.* **2014**, 1–12 (2014).

36. Farinotti, D. et al. A consensus estimate for the ice thickness distribution of all glaciers on Earth. *Nat. Geosci.* **12**, 168–173 (2019).

37. 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.* **123**, 1094–1115 (2018).

38. Bisset, R. R. et al. Reversed surface-mass-balance gradients on himalayan debris-covered glaciers inferred from remote sensing. *Remote Sens.* **12**, 1–19 (2020).

39. Gao, H. et al. Post-20th century near-steady state of Batura Glacier: observational evidence of Karakoram Anomaly. *Sci. Rep.* **10**, 1–10 (2020).

40. Kääb, A. & Funk, M. Modelling mass balance using photogrammetric and geophysical data: A pilot study at Griesgletscher, Swiss Alps. *J. Glaciol.* **45**, 575–583 (1999).

41. Zemp, M. et al. Global glacier mass balances and their contributions to sea-level rise from 1961 to 2016. *Nature* **568**, 382–386 (2019).

42. The RGI Consortium. *Randolph Glacier Inventory—A Dataset of Global Glacier Outlines* Version 6.0, GLIMS Technical Report (2017).

43. 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.* **61**, 646–662 (2015).

44. Zhao, H. et al. Dramatic mass loss in extreme high-elevation areas of a western Himalayan glacier: observations and modeling. *Sci. Rep.* **6**, 30706 (2016).

45. Sakai, A. & Fujita, K. Contrasting glacier responses to recent climate change in high-mountain Asia. *Sci. Rep.* **7**, 13717 (2017).

46. Maussion, F. et al. Precipitation seasonality and variability over the Tibetan Plateau as resolved by the High Asia Reanalysis. *J. Clim.* **27**, 1910–1927 (2014).

47. Farinotti, D. et al. Manifestations and mechanisms of the Karakoram glacier Anomaly. *Nat. Geosci.* **13**, 8–16 (2020).

48. Rupper, S. & Roe, G. Glacier changes and regional climate: a mass and energy balance approach. *J. Clim.* **21**, 5384–5401 (2008).

49. Scherler, D., Bookhagen, B. & Strecker, M. R. Hillslope-glacier coupling: the interplay of topography and glacial dynamics in High Asia. *J. Geophys. Res. Earth Surf.* **116**, 1–21 (2011).

50. de Kok, R. J., Kraaijenbrink, P. D. A., Tuinenburg, O. A., Bonekamp, P. N. J. & Immerzeel, W. W. Towards understanding the pattern of glacier mass balances in High Mountain Asia using regional climatic modelling. *Cryosphere* **14**, 3215–3234 (2020).

51. Biemans, H. et al. Importance of snow and glacier meltwater for agriculture on the Indo-Gangetic Plain. *Nat. Sustain.* **2**, 594–601 (2019).

52. Roy, J. et al. *In The Hindu Kush Himalaya Assessment: Mountains, Climate Change, Sustainability and People* (eds. Wester, P., Mishra, A., Mukherji, A. & Shrestha, A. B.) 99–125 (Springer International Publishing, 2019). https://doi.org/10.1007/978-3-319-92288-1_4

53. Hoelzle, M. et al. Re-establishing a monitoring programme for glaciers in Kyrgyzstan and Uzbekistan, Central Asia. Geosci. Instrumentation. *Methods Data Syst.* **6**, 397–418 (2017).

54. Gardner, A. S. et al. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009. *Science* **340**, 852–857 (2013).

55. Huss, M. & Hock, R. Global-scale hydrological response to future glacier mass loss. *Nat. Clim. Chang.* **8**, 135–140 (2018).

56. Pritchard, H. D. et al. Towards Bedmap Himalayas: Development of an airborne ice-sounding radar for glacier thickness surveys in High-Mountain Asia. *Ann. Glaciol.* **61**, 35–45 (2020).

57. Altena, B., Scambos, T., Fahnestock, M. & Kääb, A. Extracting recent short-term glacier velocity evolution over southern Alaska and the Yukon from a large collection of Landsat data. *Cryosphere* **13**, 795–814 (2019).

58. Cuffey, K. M. & Paterson, W. S. B. *The Physics of Glaciers* (Elsevier B.V., 2010).

59. Gardner, A. S. et al. Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years. *Cryosphere* **12**, 521–547 (2018).

60. Huss, M., Bauder, A., Werder, M., Funk, M. & Hock, R. Glacier-dammed lake outburst events of Gornersee, Switzerland. *J. Glaciol.* **53**, 189–200 (2007).

61. Mae, S. Ice Temperature of Khumbu Glacier. *J. Jpn. Soc. Snow Ice* **38**, 37–38 (1976).

62. Liu, Y., Hou, S., Wang, Y. & Song, L. Distribution of borehole temperature at four high-altitude alpine glaciers in Central Asia. *J. Mt. Sci.* **6**, 221–227 (2009).

63. Zhang, T. et al. Observed and modelled ice temperature and velocity along the main flowline of East Rongbuk Glacier, Qomolangma (Mount Everest), Himalaya. *J. Glaciol.* **59**, 438–448 (2013).

64. Vincent, C. et al. Reduced melt on debris-covered glaciers: investigations from Changri Nup Glacier. *Nepal. Cryosphere* **10**, 1845–1858 (2016).

65. Miles, K. E. et al. Polythermal structure of a Himalayan debris-covered glacier revealed by borehole thermometry. *Sci. Rep.* **8**, 1–9 (2018).

66. Gilbert, A. et al. The influence of water percolation through crevasses on the thermal regime of Himalayan mountain glaciers. *Cryosphere* **14**, 1273–1288 (2020).

67. King, O., Dehecq, A., Quincey, D. & Carrivick, J. Contrasting geometric and dynamic evolution of lake and land-terminating glaciers in the central Himalaya. *Glob. Planet. Change* **167**, 46–60 (2018).

68. Linsbauer, A., Paul, F. & Haeberli, W. Modeling glacier thickness distribution and bed topography over entire mountain ranges with glabtop: application of a fast and robust approach. *J. Geophys. Res. Earth Surf.* **117**, 1–17 (2012).

69. Huss, M. & Farinotti, D. Distributed ice thickness and volume of all glaciers around the globe. *J. Geophys. Res.* **117**, 1–10 (2012).

70. NASA/METI/AIST/Japan Spacesystems & U.S./Japan ASTER Science Team. ASTER Global Digital Elevation Model V003. (2018). https://doi.org/10.5067/ASTER/ASTGTM.003

71. WGMS. Fluctuations of Glaciers Database (2019). https://doi.org/10.5904/wgms-fog-2019-12
72. Dice, L. R. Measures of the Amount of Ecologic Association Between Species. *Ecology* **26**, 297–302 (2012).

73. Huss, M., Jouvet, G., Farinotti, D. & Bauder, A. Future high-mountain hydrology: a new parameterization of glacier retreat. *Hydrol. Earth Syst. Sci.* **14**, 815–829 (2010).

74. Banerjee, A., Patil, D. & Jadhav, A. Possible biases in scaling-based estimates of glacier change: a case study in the Himalaya. *Cryosphere* **14*, 3235–3247 (2020).

75. Stöckli, R., Vermote, E., Saleous, N., Simmon, R. & Herring, D. The Blue Marble Next Generation—A true color earth dataset including seasonal dynamics from MODIS (2005).

## Acknowledgements

This study would not have been possible without open data policies by key studies<sup>15,16,36,71</sup>, each representing the work of a large body of scientists. We thank Mohd Farooq Azam and two anonymous reviewers whose comments and suggestions improved the study. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement No 772751, RAVEN, "Rapid mass losses of debris-covered glaciers in High Mountain Asia". We acknowledge geospatial data provided by the USGS, NASA, and NASA/METI/AIST/Japan Spacesystems, as well as glaciological data provided by the NSIDC and WGMS.

## Author contributions

E.M., M.M., A.D., and F.P. designed the study and developed the methods. M.K. and S.F. supported with the collection of reference values of SMB and ELA. E.M. performed all calculations and led the writing of the paper, and all authors contributed to the interpretation of results and writing.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-021-23073-4.

**Correspondence** and requests for materials should be addressed to E.M.

**Peer review information** *Nature Communications* thanks Mohd Azam and other, anonymous reviewers for their contributions to the peer review of this work. Peer review reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021