Health and sustainability of glaciers in High Mountain Asia

High Mountain Asia (HMA) glaciers and seasonal snow provide meltwater sufficient for nearly 250 million people, making these basins among the most vulnerable globally to climatic, societal, and environmental change. Direct observations of surface mass balance (SMB) are scarce due to remoteness and danger of access, and are biased toward lower-altitude ablation areas. Accumulation rates at high elevations are poorly constrained because reanalyses have high uncertainty and in situ measurements are rare. Consequently, regional glacier models are overparameterized, lack representation of key processes (avalanching, debris cover), and are typically calibrated to sparse or spatially integrated datasets (elevation change, area change), leading to large uncertainties in current state and future projections. Remote sensing has improved regional mass change estimates, but elevation change alone cannot resolve local SMB because thinning/thickening integrate ice flow. Knowing spatially resolved SMB is critical to understand drivers of glacier change and to calibrate models appropriately, especially in HMA where avalanches and supraglacial debris can invert or modify typical altitudinal SMB gradients. This study addresses these gaps by providing spatially distributed, altitudinal SMB for HMA glaciers (2000–2016), deriving equilibrium line altitudes (ELAs), accumulation area ratios (AARs), and the fraction of annual ablation balanced by accumulation as indicators of glacier health, and by assessing implied committed changes in glacier volume and meltwater discharge by 2100 from the current climatic-geometric imbalance.

Recent remote sensing studies quantified regional mass change and glacier dynamics in HMA, yet cannot resolve SMB spatial patterns because elevation change reflects both SMB and dynamic flux divergence. Traditional models simplify SMB to linear altitudinal gradients and are calibrated to sparse field data or integrated signals, which is problematic in HMA where avalanching and debris cover are common and can produce reversed SMB gradients. Few direct high-altitude accumulation measurements exist, adding uncertainty to reanalyses and models. Earlier works reported heterogeneous regional mass balances, including the Karakoram Anomaly (near-neutral or positive mass balance in parts of Karakoram and Kunlun Shan), acceleration of mass loss in the Himalaya, and the influence of proglacial lakes on frontal ablation. Consensus ice thickness products and multi-year velocity datasets have become available, enabling continuity-based SMB derivation at scale. However, assumptions such as a single density conversion (e.g., 850 kg m−3) for converting volume to mass can bias geodetic mass balance, especially for glaciers strongly out of balance.

The study reconstructs spatially distributed annual SMB for 2000–2016 by solving the glacier continuity equation at pixel scale: b = dH/dt − ∇·q/ρ_eff, where dH/dt is surface elevation change, q is ice flux, and ρ_eff accounts for effective density of elevation change vs ice flux components. Inputs: (1) ASTER-based annual surface lowering trends (30 m) for 2000–2016; (2) ITS_LIVE multi-year mean surface velocity vectors (~240 m); (3) multi-model consensus ice thickness (25 m) tied to RGI 6.0 outlines; (4) supraglacial debris extents; (5) DEM (ASTER GDEM) for hypsometry. Processing grid: local projection and resolution varying by glacier size (50 m <15 km²; 100 m up to 80 km²; 200 m >80 km²). Datasets were reprojected and resampled (cubic spline for dH/dt and velocity endpoints; thickness degraded to analysis grid). Ice flux and flux divergence: Ice flux q = h·ȳ·u_surf, where h is ice thickness, u_surf surface velocity, and ȳ is a depth-averaging correction factor representing relative contributions of basal sliding and internal deformation. For each glacier, ȳ is estimated via a 10,000-run Monte Carlo drawing from distributions of ice thickness (from consensus), Glen’s n ~ N(3, 0.067), and basal sliding fraction U(0,1), yielding glacier-specific mean and uncertainty for ȳ. Flux divergence ∇·q is computed with a centered-difference operator on the flux components. Density correction: Flux divergence is assumed at ice density 900 kg m−3. Effective density assigned to dH/dt depends on sign combinations of dH/dt and ∇·q: accumulation (both positive) uses 600 kg m−3; ablation (both negative) uses 900 kg m−3; near-neutral/ambiguous cases use 850 kg m−3. Density uncertainty is taken as ±60 kg m−3. Uncertainty propagation: A Monte Carlo with 1000 runs per glacier perturbs inputs by (1) random uncorrelated noise on dH/dt (using reported σ); (2) systematic scaling of u (68th percentile of reported velocity errors per glacier); (3) systematic scaling of h (68th percentile inter-model standard error) plus (4) random h noise N(0, 10 m); (5) systematic scaling of ȳ (from its MC); and (6) random density variations. Pixel-scale SMB and uncertainties are aggregated. Quality control and aggregation: Calculations are per-pixel but aggregated to 25 m elevation hypsometric bins after smoothing the DEM to mitigate mismatches between velocity and thickness patterns and issues in accumulation areas, icefalls, and tributary junctions. The study focuses on glaciers ≥2 km² with measurable velocity; known surge-type glaciers and those with erratic or inverted dH/dt/SMB profiles are excluded. Additional criteria ensure reliable profiles (e.g., ELA classification accuracy ≥0.5, detrended SMB SD <3 m w.e. a−1, mean SMB uncertainty <3 m w.e. a−1), resulting in 5527 glaciers (71% of total volume of >2 km² glaciers). Deriving ELA and AAR: Pixels are labeled as accumulation or ablation based on SMB sign. Candidate ELAs at each integer elevation are evaluated as binary classifiers; the optimal ELA maximizes segmentation accuracy (Dice coefficient). For glaciers with net loss/gain at all elevations, a linear SMB trend is fitted and the elevation of SMB=0 is extrapolated as a theoretical climatic ELA. AAR is computed as fraction of glacier area above the derived ELA. Ablation balance ratio: For each glacier, total ablation is the sum of negative SMB over all pixels. Imbalance ablation is glacier-wide specific mass balance times area. Balance ablation = total ablation − imbalance ablation. The balance ratio (%) = 100 × balance ablation / total ablation; values >100% indicate net accumulation glaciers. Basin-scale totals account for sampling by scaling with independent imbalance ablation from regional dH/dt and using subset balance:imbalance ratios. Implied committed change to 2100: A Δh parameterization (after Huss & Hock and related work) redistributes annual mass change across glaciers using observed thinning patterns to emulate dynamic response. For glaciers with clear loss, Δh is derived from observed thinning; for ambiguous cases, a standard Δh curve is used. Advance rules: glaciers advance when terminus longitudinal gradient >10°, appending lowest-elevation marginal pixels and thickening by prior terminus height; advancing fraction limited to 50% of volume gain; remaining gain redistributed altitudinally per Δh. SMB is extrapolated beyond current termini using a median ablation gradient of 0.07 m w.e. m−1. Simulations run annually for 200 years; outcomes for 2100 are emphasized. To reflect primary uncertainty (dH/dt), 30 runs per glacier vary SMB within reported dH/dt uncertainty; regional means and standard deviations are reported. Regional aggregations follow established HMA subregions and major river basins, with area-weighted and volume-weighted statistics and uncertainties as described.

• Dataset: Altitudinally resolved SMB maps for 5527 HMA glaciers (2000–2016), representing 71% of ice volume (56% of area) of glaciers >2 km² after quality filtering.
• Accounting for dynamics vs geodetic thinning: Elevation change alone biases SMB; correcting for ice flux divergence reveals a systematic bias of about +0.07 m w.e. a−1 in past mass balance estimates that used a single density.
• Equilibrium line altitude (ELA) and AAR: Area-weighted mean ELA for HMA is 5283 m a.s.l.; regional mean AAR is 0.51. ELA and AAR vary widely (SDs: 678 m and 0.32, respectively). Median glacier elevation is a poor proxy for ELA (median absolute deviation ~193 m).
• Extent of accumulation areas: 16% of studied glaciers (10% of area) have no accumulation area (AAR = 0). 32% have AAR <0.1 (19% of area); 41% have AAR <0.2 (23% of area). Among glaciers >5 km², 7.5% still have AAR <0.1. Low-AAR glaciers are concentrated in Eastern Nyainqentanglha and are common across the Himalaya, Tibetan Plateau, and Tien Shan.
• Karakoram Anomaly context: Many glaciers in Karakoram and Kunlun Shan have AAR >0.5, indicating healthier states, yet ELAs are not anomalously low there relative to neighbors, suggesting topographic availability of high-elevation accumulation areas as a key factor.
• Sustainable vs unsustainable ablation: Regionally, 40 ± 11% of annual glacier ablation is imbalance (unsustainable) ablation; equivalently, about 60% is balanced by accumulation. Using basin scaling, the regional balanced fraction is 60% (subset mean ~50%). Key basins with Karakoram Anomaly influence (Indus, Amu Darya, Syr Darya, Tarim Interior) have >50% balanced ablation; the Indus basin shows 65% ± 23% balanced ablation (2000–2016). The Ganges–Brahmaputra basin has 48% ± 9% balanced ablation, implying vulnerability in drought seasons.
• Glacier health gradients: Subregions like Nyainqentanglha exhibit the most negative SMB, with accumulation confined to the highest ~20% of elevation range; Everest, Spiti Lahaul, and Tien Shan show similar normalized SMB profiles (accumulation over upper ~40%). Karakoram and Kunlun Shan accumulate over upper ~60% and have less negative SMB in ablation areas.
• Committed changes by 2100 from current disequilibrium: Implied regional glacier volume change by 2100 is −23% ± 1% (main analysis; Figure 4). The abstract reports −21% ± 1%. About 25% of glaciers are projected to lose at least 50% of their volume by 2100 without additional warming (the concluding section also notes ~35% as very unhealthy and expected to lose at least half). Subregional outcomes: all Himalayan subregions lose ≥35% volume; Karakoram, Pamir Alai, Pamir, and Hindu Kush lose 10–20%; Kunlun Shan slightly gains volume (+2.1%).
• Committed change in melt supply: Total annual glacier ablation decreases by ~28% by 2100 (−28% ± 6% in body; abstract gives −28% ± 1%), with stronger reductions where volume losses are large; some subregions (Karakoram, Pamir) show ablation decreases stronger than volume losses.
• Water security implications: Even in basins with relatively healthy glaciers (e.g., Indus), significant fractions of ablation are imbalance-derived; tributaries vary (e.g., Chenab and Satluj fed by unhealthy Western Himalayan glaciers vs Karakoram-fed segments). Ganges–Brahmaputra glaciers’ imbalance-dominated ablation threatens pre-monsoon water availability in drought years.

The study addresses the key uncertainty in HMA glacier change by resolving spatially distributed SMB, separating local climatic inputs/outputs from dynamic mass redistribution. The findings show that most HMA glaciers have small to negligible accumulation areas, indicating strong climatic-geometric imbalance and implying committed mass loss and reduced meltwater supply even without further warming. The Karakoram Anomaly’s elevated AARs without correspondingly low ELAs suggest that topography (high-elevation area availability) and recent high-altitude precipitation increases disproportionately support glacier health there. Basin-scale sustainable ablation fractions highlight contrasting water security: basins influenced by Karakoram/Kunlun glaciers currently fare better, but imbalance ablation remains substantial, and long-term reductions in both ice mass and meltwater are inevitable. The derived SMB gradients and ELAs provide new, spatially explicit benchmarks for glacier model calibration and evaluation, emphasizing that models should incorporate localized processes (avalanching inputs, debris-modified ablation, frontal/lacustrine ablation) and be constrained by both thinning and velocity observations to avoid compensating errors in melt vs accumulation. The committed losses are already comparable to low-emissions projections, implying that much of the 21st-century mass loss is locked in by current disequilibrium, with continued warming expected to exacerbate declines and shifts in seasonal water availability.

This work delivers the first region-wide, altitudinally resolved SMB dataset for 5527 HMA glaciers (2000–2016), enabling derivation of ELAs, AARs, balanced vs imbalance ablation, and implied committed changes. Key contributions include: (1) demonstrating the necessity of accounting for ice dynamics and effective density to avoid SMB bias; (2) quantifying widespread low AARs and high ELAs indicating unhealthy glacier states outside Karakoram/Kunlun; (3) providing basin-specific balance ratios that reframe water security under current climate; and (4) estimating committed regional volume loss (~21–23% by 2100) and melt reduction (~28%), with stronger losses in Himalayan and Tien Shan subregions. The dataset offers a powerful constraint for next-generation glacier models and monitoring strategies. Future research should prioritize: disentangling mechanisms and resilience of the Karakoram Anomaly; incorporating avalanche accumulation, debris influences, and frontal ablation in models; expanding high-altitude accumulation and ice thickness measurements; improving velocity retrievals over small and complex-flow glaciers; and producing multitemporal SMB reconstructions to track evolving glacier-climate coupling.

• Input data limitations: Ice thickness models are the most uncertain input; velocity products can decorrelate in accumulation areas and across icefalls; elevation change products have uncertainties and potential biases. Small glaciers (<2 km²) were excluded due to weak velocity signals; surge-type and glaciers with erratic profiles were removed, potentially biasing regional sampling.
• Density assumptions: Effective density assignment mitigates bias, but uncertainties remain (±60 kg m−3) and can affect SMB retrievals, especially near-neutral zones.
• Basal conditions and rheology: Poorly constrained across HMA; Monte Carlo assumptions (sliding fraction U[0,1], n≈3) introduce uncertainty in depth-averaged velocity factors.
• Δh parameterization: Simplifies ice dynamics and advance/retreat behavior; results depend on assumed thinning patterns, advance thresholds, and SMB extrapolation beyond termini; robust for multidecadal scales but still a parameterization.
• Regional scaling for basins: Estimating total ablation and balanced fractions requires scaling from the analyzed subset using independent mass-balance and sampling corrections, adding uncertainty, especially in basins with surging glaciers (e.g., Karakoram region).
• Temporal representativeness: SMB represents 2000–2016 mean conditions; committed-change estimates assume persistence of this regime and thus provide baseline/minimum changes, not full climate-driven projections.