High Mountain Asian glacier response to climate revealed by multi-temporal satellite observations since the 1960s

High Mountain Asia (HMA) hosts the largest concentration of glaciers outside the polar regions and acts as the water tower of Asia, supporting major rivers through glacier meltwater. Anticipated temperature increases are expected to drive continued glacier recession, with implications for runoff timing and magnitude. However, projections of glacier change carry large uncertainties due to limited long-term observations and sparse in situ mass balance records (≈30 glaciers with data; only two >30 years, both in the Tien Shan). Geodetic methods using satellite-derived elevation changes have improved understanding since 2000, revealing widespread mass loss but regional heterogeneity (e.g., regions of near balance in the Karakoram, Western Kunlun, Eastern Pamir). Earlier use of Hexagon KH-9 imagery extended records back to the 1970s, but accuracy and data voids limited temporal resolution. High-resolution declassified Corona KH-4 imagery offers an opportunity to extend and densify time series into the 1960s. This study aims to reconstruct multi-temporal glacier mass budgets across seven climatically distinct HMA regions over nearly six decades and to attribute variability to climatic drivers, particularly summer temperature and solid precipitation.

Previous geodetic assessments since 2000 indicated heterogeneous glacier mass changes across HMA consistent with global glacier recession, with exceptions termed the Karakoram anomaly where some regions maintained near-balanced mass budgets. Studies leveraging Hexagon KH-9 imagery provided insights back to the mid-1970s and suggested increased mass loss post-2000, but typically covered only two extended periods and suffered from image contrast limitations and large data voids. Earlier Corona-based applications demonstrated potential but were limited to localized areas due to geometric distortions and complex processing. Recent region-wide analyses (e.g., Brun et al. 2017; Shean et al. 2020) provide contemporary baselines, while localized studies in Tien Shan, Central Himalaya (e.g., Langtang, Everest), and Pamir corroborate acceleration in mass loss and highlight climatic controls and non-climatic influences (surging, proglacial lakes).

Study design: Seven climatically diverse glacierized regions across HMA were analyzed: Northern Tien Shan; Ak-Shirak (Central Tien Shan); Purogangri Ice Cap (central Tibetan Plateau); Western Nyainqentanglha; Poiqu region including the Langtang sub-region (Central Himalaya); Gurla Mandhata (West-Central Himalaya); and Muztagh Ata Massif (Eastern Pamir). Glacier outlines: RGI v6.0 outlines were used as a baseline and manually refined using orthoimages and DEM-derived slope/curvature. Outlines were updated for each epoch. DEM sources and generation: Multi-temporal DEMs were generated from declassified Corona KH-4 (1960s–early 1970s) using the Remote Sensing Software Package Graz with a modified collinearity distortion model and GCPs (~30; Landsat and SRTM references). Hexagon KH-9 (mid-1970s) DEMs were produced in Leica Photogrammetry Suite using a frame camera model, Brown’s physical model for distortions, reseau reconstruction, and GCPs (~40). Contemporary high-resolution optical stereo (Pléiades, SPOT-6/7, GeoEye) were processed in PCI Geomatica using RPCs, GCPs (RMSE ≤1 pixel), tie points, and semi-global matching (1–5 m DEMs). ASTER L1A stereo DEMs were generated with NASA Ames Stereo Pipeline (30 m). TanDEM-X interferometric SAR DEMs were produced with GAMMA software (4×4 multilooking, differential interferograms using global TanDEM-X DEM, MCF phase unwrapping) and vertically adjusted over stable terrain. Co-registration, bias correction, outlier removal: All DEMs were co-registered to void-filled 30 m SRTM following Nuth & Kääb. Elevation-dependent biases were corrected with first-order polynomial surfaces. DEM differences with |Δz| > ±150 m were removed; additional filtering within 100 m elevation bins using µ ± 3σ thresholds was applied regionally or glacier-wise as needed. An elevation-dependent sigmoid constraint was applied to limit unrealistic thickness changes by altitude. Gap filling: Small voids were filled with a 4×4 moving-window mean; larger voids were filled hypsometrically using median values per elevation bin, treating surge-type and non-surge-type glaciers separately. Radar penetration correction: X-band penetration biases in TanDEM-X DEMs were quantified by comparing contemporaneous optical DEMs. Estimated mean penetration and propagated mass balance effects: Gurla Mandhata ~0.68 m (clean ice), ~0.24 m (debris), effect ~0.02 m w.e. a−1 (2013–2016); Purogangri ~1.56 m, effect ~0.04 m w.e. a−1 (2012–2018). For Western Nyainqentanglha, ±0.03 m w.e. a−1 was adopted based on C-band estimates. Seasonality corrections: Acquisition-date offsets were minimized; where needed, small corrections were applied based on regional accumulation data (e.g., Ak-Shirak +0.10 m w.e. for 1964–1973; Muztagh Ata +0.05 and +0.03 m w.e. for specified periods), considering glacier accumulation type and seasonality. Conversion to mass and uncertainties: Volume-to-mass conversion used an ice density of 850 ± 60 kg m−3, with larger density uncertainties for short periods (≤3 years: Δρi 100–150 kg m−3). Elevation change uncertainty was computed from off-glacier σ per altitude band accounting for spatial autocorrelation (Sa ≈ 600 m). Total thickness change uncertainty combined DEM differencing and area uncertainties; overall mass balance uncertainty UM combined thickness and density terms. Climate data and attribution: ERA5-Land (1981–present) at glacierized elevations provided 2 m air temperature and precipitation; solid precipitation (SolP) was defined as precipitation when T < 0 °C. HAR v2 data were used for validation in some regions. In situ station records (since the 1970s where available) were compiled to validate ERA5-Land. Primary variables analyzed were summer temperature (SumT) and solid precipitation (SolP). Correlation analyses were conducted between geodetic mass balance for each time period and climate anomalies, both regionally and by clustering into humid (Poiqu, Gurla Mandhata, Western Nyainqentanglha) and cold/dry (Ak-Shirak, Purogangri, Muztagh Ata) regimes. Additional analyses considered non-climatic factors (surge-type glacier behavior, proglacial lake expansion) using literature and mapped changes. Area and length changes: Glacier area changes were computed with mapping uncertainty propagated from image resolution and rectification RMSE. Terminus length changes were derived using a centerline-based reference point method with uncertainty from image resolution and co-registration.

- Widespread and accelerating mass loss: All seven regions experienced net mass loss over the multi-decadal record, with acceleration in most since the early 2000s. Regional mass budgets spanned from approximately -0.40 ± 0.07 m w.e. a−1 (Central/Northern Tien Shan) to -0.06 ± 0.07 m w.e. a−1 (Eastern Pamir/Muztagh Ata) over full study periods. - Region-specific trends and rates: • Northern Tien Shan: Mass loss intensified from ~-0.18 ± 0.11 m w.e. a−1 (1964–1971) to -0.49 ± 0.13 m w.e. a−1 (2016–2020). Contemporary (to 2016) loss consistent with earlier studies (~-0.41 ± 0.20 m w.e. a−1). • Ak-Shirak (Central Tien Shan): Persistent loss across the record, peaking at -0.54 ± 0.10 m w.e. a−1 (1980–2002); full-period ~-0.40 ± 0.07 m w.e. a−1. • Purogangri Ice Cap: Moderate long-term loss (-0.15 ± 0.07 m w.e. a−1, 1969–2019) with near-balance during 2000–2012 (-0.03 ± 0.02 m w.e. a−1), then renewed loss (e.g., 2012–2019 ~-0.12 ± 0.05 m w.e. a−1). • Western Nyainqentanglha: Increasing loss from ~-0.26 ± 0.13 m w.e. a−1 (1968–2001) to ~-0.43 ± 0.12 m w.e. a−1 (2001–2019); recent period ~-0.46 ± 0.14 m w.e. a−1 (2012–2019). • Poiqu region (Central Himalaya) including Langtang: Poiqu loss increased from -0.30 ± 0.10 m w.e. a−1 (1974–2004) to -0.42 ± 0.11 m w.e. a−1 (2004–2018). Langtang intensified from -0.20 ± 0.09 m w.e. a−1 (1964–1974) to ~-0.58 ± 0.11 m w.e. a−1 (2015–2019); mean mass budget roughly doubled from -0.23 ± 0.10 (1964–2004) to -0.50 ± 0.11 m w.e. a−1 (2004–2019). • Gurla Mandhata: Near-balance around 2011–2013 (-0.02 ± 0.08 m w.e. a−1) following slight long-term loss (-0.12 ± 0.10 m w.e. a−1, 1966–2011), then increased loss post-2013 (~-0.20 ± 0.11 m w.e. a−1, 2013–2016; -0.24 ± 0.17 m w.e. a−1, 2018–2019). • Muztagh Ata Massif (Eastern Pamir): Little change to slightly positive between 1973 and 2009 (including +0.03 ± 0.10 m w.e. a−1 during 2001–2009), then transition to loss after 2013 (~-0.12 ± 0.09 to -0.11 m w.e. a−1, 2013–2019); full period near balance to slight loss (~-0.06 ± 0.07 m w.e. a−1, 1967–2019). - Area changes mirrored mass changes: Highest area loss in Northern Tien Shan (-0.60 ± 0.07 % a−1) and Western Nyainqentanglha (-0.30 ± 0.02 % a−1). Minimal area change where mass budgets were close to balance for decades (Gurla Mandhata: -0.08 ± 0.03 km2 a−1; Muztagh Ata: -0.07 ± 0.02 km2 a−1). Other regions showed shrinkage of similar magnitude (Ak-Shirak ~-0.19 ± 0.01 % a−1; Purogangri ~-0.17 ± 0.01 % a−1; Poiqu ~-0.18 ± 0.01 % a−1). - Climate attribution: • Summer temperature (SumT) increases are the dominant driver of long-term acceleration in glacier mass loss across regions, with strong correlations in Northern Tien Shan (ERA5-Land: r ≈ 0.97, p = 0.009; in situ: r ≈ 0.92, p = 0.004). • In cold/dry regions (Ak-Shirak, Purogangri, Muztagh Ata), solid precipitation (SolP) is also important (ERA5-Land: R2 ≈ 0.79, p = 0.0007), with decadal near-balance at Purogangri (2000–2012) linked to increased SolP (+42 mm a−1, ~15%) and slight SumT decrease (~-0.24 °C). • In humid/monsoonal regions (Poiqu, Gurla Mandhata, Western Nyainqentanglha), both SumT and precipitation (solid, winter, and annual) correlate with mass balance (e.g., ERA5-Land SumT R2 ≈ 0.71, p = 0.003; SolP R2 ≈ 0.64, p = 0.0001; winter precipitation R2 ≈ 0.66, p = 0.003). • Recent transitions from balanced to negative mass budgets (Muztagh Ata, Purogangri, Gurla Mandhata) coincide with pervasive increases in SumT in the 2010s. - Non-climatic influences: Surge-type glacier dynamics in Ak-Shirak elevated post-surge ablation and sustained higher loss rates after 1980 (surge-type: -0.62 ± 0.10 m w.e. a−1, 1980–2002). Proglacial lake expansion (e.g., Petrov Lake +240% area; Poiqu basin widespread lakes) enhanced calving and subaqueous melt, exacerbating mass loss relative to regional means. - Consistency with prior studies: Contemporary (2000–2019) geodetic mass balances agree closely with regional assessments (differences typically within ±0.02 to ±0.04 m w.e. a−1) and with localized, temporally detailed studies, while extending records back to the 1960s and adding intermediate periods.

The study links six decades of geodetic glacier mass balance to climatic variability across diverse HMA regions. The findings indicate that rising summer temperatures are the principal driver of accelerated glacier mass loss, consistent across both humid monsoonal and cold/dry regimes. In cold/dry settings, solid precipitation modulates decadal variability, occasionally offsetting temperature-driven losses (e.g., Purogangri near 2000–2012). In humid regions, precipitation (especially winter and solid precipitation) and summer temperature together influence mass budgets, but recent warming appears to override prior quasi-equilibrium states (e.g., Gurla Mandhata post-2013). The observed transitions from balanced to negative mass budgets in previously stable regions underscore increasing climatic sensitivity, particularly to summer temperature anomalies. Non-climatic processes—surging and glacier–lake interactions—introduce additional spatial/temporal heterogeneity and can amplify ablation beyond climatic forcing. The agreement with independent studies strengthens confidence in the reconstructed trajectories and emphasizes the importance of sub-decadal temporal resolution to capture nonlinear responses. Overall, the results substantiate that recent, pervasive summer warming is driving widespread mass loss and the end of the previously observed ‘Karakoram anomaly’-like stability in several subregions.

This work provides a temporally detailed, multi-sensor geodetic record of glacier mass change across seven HMA regions from the 1960s to 2019 by exploiting declassified Corona KH-4, Hexagon KH-9, and modern optical and SAR data. It demonstrates persistent and accelerating mass loss across most regions, including recent transitions from near-balance to loss in historically stable areas (Eastern Pamir, Purogangri, Gurla Mandhata). Attribution analyses identify summer temperature increases as the dominant driver, with precipitation variability modulating responses, particularly in cold/dry and humid regimes. Non-climatic processes such as surging and proglacial lake expansion further influence regional heterogeneity. The dataset and methods enable improved calibration/validation of glacier mass balance models and highlight the need for sub-decadal monitoring. Future research should prioritize enhanced high-elevation in situ observations (especially precipitation), refined representation of surge and lake-terminating dynamics, and continued development of high-resolution, bias-corrected DEM time series to better constrain projections of glacier recession and downstream hydrological impacts.

- Climate data limitations: ERA5-Land precipitation shows weak agreement with in situ measurements, and station networks are sparse and often at lower elevations than glacierized terrain, limiting representativeness of accumulation. Temperature agreement is stronger but still variable by region. - Remote sensing constraints: Declassified imagery (KH-9) suffers from contrast/dynamic range issues and data voids, particularly at high elevations; Corona KH-4 requires complex distortion modeling. Radar (TanDEM-X) penetration into dry snow/ice introduces elevation biases that require corrections and add uncertainty. DEM co-registration and differencing are sensitive to terrain-dependent biases and spatial autocorrelation. - Conversion and temporal resolution: Assumptions on ice density (and its variability over short periods and surge cycles) affect mass conversion; short observation periods (≤3 years) carry larger density-related uncertainties. Seasonality corrections are approximate where acquisition dates differ. - Spatial scope: Only seven regions are analyzed; while climatically diverse, they may not capture all HMA heterogeneity. Non-climatic processes (surging, lake calving) are not exhaustively quantified across all glaciers. - Precipitation phase partitioning: Solid precipitation definition via temperature thresholds may misclassify precipitation phase under certain meteorological conditions.