Glacier mass-balance estimates over High Mountain Asia from 2000 to 2021 based on ICESat-2 and NASADEM

Abstract

High Mountain Asia (HMA) glaciers are critical water reserves for montane regions, which are readily influenced by climate change. The glacier mass balance during 2000–2021 over HMA was estimated by comparing the elevations from ICESat-2 and the NASADEM. Radar penetration depth could be one of the intrinsic error sources in estimating glacier mass balance by using NASADEM. Therefore, we doubled elevation differences between the X-band Shuttle Radar Topography Missions (SRTMs) and NASADEM to estimate the potential error. The spatial characteristics of the altitude-dependent penetration depth can be detected in most sub-regions of HMA. Relatively deep penetrations in the Himalaya (2.3–3.7 m) and Hissar Alay (4.3 m) regions and small penetrations in the south-eastern HMA (1.0 m) were observed. The HMA region experienced a significant mass loss at a rate of −0.18 ± 0.12 m w.e. a⁻¹, in which the Hengduan Shan exhibited the highest mass loss of −0.62 ± 0.10 m w.e. a⁻¹, the West Kun Lun experienced a substantial mass gain of 0.23 ± 0.13 m w.e. a⁻¹, and the Karakoram showed a more or less balance. Our results are in agreement with previous studies that assessed the mass balance of HMA glaciers from different methods.

1. Introduction

High Mountain Asia (HMA) contains the densest concentration of alpine glaciers that are sensitive to global warming. Retreat and thinning of alpine glaciers have changed meltwater flow routines and water availability thus affected daily life of the local communities (Immerzeel and others, 2020). The rapid ice mass loss has also contributed to the rise of global sea level. Glacier volume change (mass balance) due to ice thickness change and redistribution can indicate climate change directly. Therefore, it is particularly significant to study glacier mass balance across the HMA.

Glaciological in situ measurements over a limited number of glaciers and subsequent extrapolation to the entire region (WGMS, 2020) may lead to large biases toward low-altitude glaciers (Fujita and Nuimura, 2011). For the alleviation of this problem, geodetic remote-sensing datasets have been introduced to glacier mass-balance analysis to partially compensate for deficiencies of in situ measurements at high elevations. Those datasets include laser altimeter data (Kääb and others, 2012, 2015; Wang and others, 2021), radar altimeter data (Jakob and others, 2021), Digital Elevation Models (DEMs) data (Brun and others, 2017; Zhou and others, 2018; Shean and others, 2020; Bhattacharya and others, 2021; Hugonnet and others, 2021), the Gravity Recovery and Climate Experiment (GRACE) and GRACE-Follow On (GRACE/FO) satellites (Wouters and others, 2019; Ciracì and others, 2020). Accuracy of altimeter-derived mass-balance estimates rely on representativeness of spatial sampling because only a few glaciers are sampled, and the whole glacier area is never fully sampled (Treichler and others, 2019). DEM has advantage of high spatial resolution and temporal coverage. A time series of advanced spaceborne thermal emission and reflection radiometer (ASTER) DEMs can provide detailed glacier elevation information of each DEM pixel. Nevertheless, ASTER-derived DEM may have biased errors or voids in featureless accumulation areas (Wang and Kääb, 2015). Therefore, ASTER DEM stacks for a short period of time were too noisy for mass-balance estimates (Treichler and others, 2019). Gravity is a potential field, and gravity signals can come from any changes in mass, including those from glaciers, permafrost and water storage. However, it is difficult to isolate glacier signals directly from the GRACE measurement (Zhang and others, 2013).

Uncertainties in glacier mass balance can be affected by penetration depth differences of different radar frequencies (Round and others, 2017), different spatial and temporal coverage/resolution of different data sources, as well as data accuracy. Penetration depth of radar signals through ice can be a significant source of uncertainty when using radar-based DEM data such as the C-/X-band Shuttle Radar Topography Mission (SRTM) and the updated NASADEM (NASA JPL, 2021). Radar penetration depths into glaciers can be estimated by different methods: (i) linear extrapolation of ICESat-derived glacier elevation change backdating to the SRTM acquisition date of February 2000 by assuming glacier elevation-change rate that showed the same trend during 2003–2008 and 2000–2003 (Kääb and others, 2012, 2015); (ii) subtracting DEMs from field measurements or optical DEM (Dehecq and others, 2016; Round and others, 2017; Lambrecht and others, 2018); (iii) adding the relative penetration differences between C-/X-band SRTM and the X-band penetration estimation from the assumption of dry snow percentage (Zhou and others, 2019); (iv) doubling the relative penetration differences between SRTM-C and X-band DEMs (Jaber and others, 2019); and (v) ignoring the penetration issues by assuming that SRTM has captured the 1999 autumn elevation profiles (Berthier and others, 2006) and choosing data acquired during autumn season to calculate the mass balance.

ICESat-2's ATLAS instrument can provide elevation measurements with a specific timestamp and without the effect of penetration. The ICESat-2 ATL06 (Land Ice Elevation) product has been applied to study the ice mass balance of Greenland and Antarctic ice sheets (Smith and others, 2020), small glaciers such as those in Svalbard (Sochor and others, 2021) and glacier elevation changes in steep and rugged regions such as HMA (Wang and others, 2021; Shen and others, 2022; Zhao and others, 2022). The ICESat-2 ATL06 product was also found to have a higher accuracy than other DEM products (Chen and others, 2022). Consequently, it is relevant to use ICESat-2 data collected over the HMA from 2018 to 2021 to study the spatiotemporal glacier mass variations. To this end, we will: (1) correct penetration depths for NASADEM to obtain penetration-corrected elevation data in 2000 by doubling the difference between the X-band SRTM DEM and NASADEM; (2) compare penetration estimations estimated with different methods; and (3) obtain the HMA glacier mass balance from 2000 to 2021 through comparing the elevation differences between ICESat-2 and penetration-corrected NASADEM.

2. Study area, data sources and methods

2.1 Study area

HMA has a glacier area of ~9.76 × 10⁴ km² (RGI Consortium, 2017) and accounts for ~7.6% of the global glacier volume, excluding the Greenland and Antarctic ice sheets (Zemp and others, 2019). The climate of the HMA region is mainly controlled by three atmospheric circulation systems: the westerlies, the Indian summer monsoon (ISM) and the East Asia monsoon (Fig. 1). The climate in southern HMA is controlled by the ISM, and that in the north-western HMA is controlled by the westerlies (Ménégoz and others, 2013). The climate in the south-eastern margin zones is additionally affected by the East Asia monsoon. Because of different atmospheric circulation patterns, glaciers have distinct accumulation periods and are consequently categorized into three main categories: winter-, spring- and summer-accumulated glaciers (Maussion and others, 2014). Different regions exhibit various percentages of debris due to local variations in geology, topographic relief and glacier dynamics (Herreid and Pellicciotti, 2020). Glaciers in the western and southern HMA have a high percentage of debris coverage, among which the highest occurs in the Hindu Kush mountain range (19.3%). However, glaciers in the inner HMA have a low percentage of debris coverage of about only 2% (Scherler and others, 2018).

Fig. 1. Glacier distribution of the High Mountain Asia region. Purple polygons indicate glaciers, sub-region boundaries are marked and the names are notated in black. The glacier and region boundaries were from the Randolph Glacier Inventory 6.0. Blue lines and notations represent the flows and names of major HMA rivers. The arrows show the three atmospheric circulation systems of HMA, including the westerlies, Indian summer monsoon and East Asia monsoon.

2.2 Data sources

2.2.1 ICESat-2 data

The advanced topographic laser altimeter system (ATLAS) instrument on the ICESat-2 satellite provides significantly increased spatial coverage compared to previous altimeters by dividing the transmitting laser pulse into six beams (Markus and others, 2017). The ICESat-2 ATL06 (version 4) product uses a window size ranging from 40 to 80 m to segment photons from global geo-located photon data (ATL03) to estimate the geolocated land ice surface height (Smith and others, 2019). Pair tracks are ~3 km apart in the across-track direction, and the strong and weak beams within one pair are separated by ~90 m. Land ice surface height is then determined after corrections for instrument bias (e.g. transmit pulse shape bias correction and first-photon bias correction) (Smith and others, 2019).

The ATL06 product contains latitude, longitude and height above the WGS84 ellipsoid for each measurement footprint. Land ice heights represent the mean surface height averaged over 40 m segments of the ground track and spaced 20 m apart (Smith and others, 2019). ATL06 product has better than 5 cm height accuracy in the Antarctic, and the biases of two beams in one pair are <2 cm (Brunt and others, 2019). Compared with in situ continuously operating reference system and unmanned aerial vehicle data acquired in Qilian Shan, ATL06 data have very high vertical and horizontal positioning accuracy (Zhang and others, 2021), which is indispensable for glacier mass-balance estimates.

We selected the ATLAS/ICESat-2 data that were acquired around the same season as the NASADEM was acquired (mid-February), and a window of 2 months around 15 February was chosen to obtain the appropriate laser altimetry data and minimize seasonal-related bias. Therefore, the data from 15 January to 15 March in 2019, 2020 and 2021, which included more than 2.53 × 10⁶ footprints on the glaciers over the HMA region, were used to estimate glacier mass balance.

2.2.2 DEM datasets

NASADEM is an update of the SRTM DEM. It was derived by reprocessing the original SRTM interferometric SAR data using updated interferometric unwrapping algorithms by applying vertical and tilt adjustments based upon GCPs derived from ICESat to improve the vertical accuracy (Crippen and others, 2016). Compared with the C-band SRTM 1 arc-second v3 product, this processing scheme can deal with strong offsets and ramps occurring in regions where the acquisition strips cross each other (Crippen and others, 2016). It has the best performance among the open-access DEMs (Chen and others, 2022), and has few voids remaining. We masked out the void-fill regions in the NASADEM data (Braun and others, 2019), as the fill-in data have a completely different time stamp and may not be inherently consistent with the original radar data on glaciers, or the data used to fill in could have already originated from optical data like ASTER Global Digital Elevation Model (GDEM) (NASA JPL, 2021).

The SRTM X-band DEM was acquired simultaneously with the NASADEM, with the same 30 m resolution; however, the SRTM-X DEM shows an ‘X’ stripe-like coverage due to a smaller swath width. The horizontal datum of the two DEM datasets is the same (WGS84 datum), but the vertical datums are different; the EGM96 geoid height datum was used for the NASADEM, and the WGS84 height datum was used for the X-band SRTM DEM. In this study, we used the X-band SRTM DEM as the reference to estimate the radar penetration depth for NASADEM.

2.2.3 Glacier boundary dataset

The Randolph Glacier Inventory (RGI) 6.0 (RGI Consortium, 2017) used images acquired near year 2000 to generate glacier boundaries all over the world. It included glaciers of all sizes, which are beneficial for glacial mass change evaluation during 2000–2021. The HMA region consists in regions 13, 14 and 15 in RGI 6.0, so we used the glacier boundaries of these three regions to classify the ICESat-2 footprints and DEM datasets into glacier and non-glacier sections.

2.3 Methods

2.3.1 Co-registration of different datasets

The EGM96 datum for the NASADEM was converted to the WGS84 height datum by adding the geoid height difference in the ENVI software. We then divided the HMA into grids of 1° × 1°, and NASADEM was registered to ICESat-2 or X-band SRTM DEM within each grid by applying the method proposed by Nuth and Kääb (2011), which is given by Eqn (1):

(1)$$\displaystyle{{{\rm d}h} \over {\tan ( \alpha ) }} = a \times \cos ( b-\beta ) + \displaystyle{{\overline {{\rm d}h} } \over {\tan ( \overline \alpha ) }}, \;$$

where dh is the elevation difference between the NASADEM and ICESat-2 (or X-band SRTM DEM), α is the slope, β is the aspect, a and b are the magnitude and direction of the shift vector, $\overline {{\rm d}h\;}$ is the mean elevation difference between the two DEM datasets and $\overline {\alpha \;}$ is the mean slope.

We calculated the shift vector by applying Eqn (1) to each grid to correct the displacement. This process was iterated 30 times, and the shift vector that minimizes the product of the std dev. and the median of the elevation difference was selected as the final offset result. Each 1° × 1° NASADEM grid was shifted to the ICESat-2 data during co-registration. The detailed offset vector was shown in Figure S1. The offset in the plane was, at most, 15 m, which was less than the pixel size of the NASADEM, and the offsets in the Z direction of 99.4% (318/320) of the grids were within ±5 m.

For the NASADEM and X-band SRTM DEMs, we corrected for the maximum curvature- and elevation-dependent bias by fitting the second-order polynomial to the elevation difference on stable ground (Gardelle and others, 2012, 2013). Filters of ±150 m in dh were used to remove outliers.

2.3.2 Representativeness of the ICESat-2 footprints

ICESat-2 samples have to match the glaciers within each spatial unit with respect to the glacier hypsometry to get a statistically robust signal (Kääb and others, 2012; Neckel and others, 2014). Glacier hypsometry results of each mountain range derived from ICESat-2 and NASADEM show that there are few differences between them. But for a smaller area of interest like 1° × 1° cells, hypsometry tends to show some discrepancy. If we found that glacier hypsometry between NASADEM and ICESat-2 is larger than 10% at any elevation band for one cell, we iteratively eliminated some footprints at this elevation band until glacier hypsometry difference at all elevations is smaller than 10% (Treichler and others, 2019; Jakob and others, 2021). We generated estimates of elevation change and mass change from the reduced samples. However, if it turns out we cannot control the hypsometry difference between ICESat-2 and NASADEM under 10%, this cell would not be taken into consideration when estimating glacier mass change. In total, 256 cells were retained after this process and 46 cells were eliminated for the 1° × 1° cell mass-balance estimation. Neckel and others (2014) found that glacier elevation change became stable when the ICESat footprint ratio exceeded 60% for each mountain range. Therefore, a bootstrapping analysis was performed by iteratively and randomly selecting ICESat-2 glacier footprints and computing the mean dh value in each iteration, the sampled footprints can represent the entire region since we found the same trend for the ICESat-2 footprints.

2.3.3 Radar penetration-depth correction

The penetration depth issue of radar signals can lead to underestimation of glacier surface elevation (Kääb and others, 2012). The radar penetration depth into the snow/ice of a glacier needs to be corrected. For C-band DEM, the C-band penetration depth was corrected by doubling the difference between the X-band SRTM DEM and the NASADEM (Jaber and others, 2019) using Eqn (2):

(2)$$k = 2 \times ( {\rm SRT}{\rm M}_{x-{\rm band}}-{\rm NASADE}{\rm M}_{c-{\rm band}}) , \;$$

where k is the C-band penetration depth in each NASADEM elevation bin.

The coordinates of the NASADEM and X-band SRTM DEM were transformed into their corresponding UTM projections to a resolution of 30 m through cubic interpolation and then were co-registered as described in section 2.3.1. We clipped the grid by the glacier boundary of the region and analyzed the penetration depth of each sub-region according to the boundary. Penetration depths exceeding ±15 m were identified as outliers (Neelmeijer and others, 2017). We divided each RGI mountain range into zones with 100 m intervals in elevation and calculated the median value of penetration depth in each zone, and added the value to the NASADEM value of each pixel in the corresponding bin. For areas of higher or lower elevations in which there were not enough X-band SRTM DEM pixels to be used in the penetration-depth correction, the median value of penetration depth of the nearest elevation zone was used instead.

We investigate alternative methods to estimate the penetration depth following the method proposed by Kääb and others (2012, 2015) for comparison of penetration estimate.

2.3.4 Mass-balance estimation

Elevation differences (dh) between the corrected NASADEM and the ICESat-2 elevation were calculated for each on-glacier footprint. The ICESat-2 footprints were enclosed according to the sub-region boundary. Any footprints with dh values exceeding ±200 m were excluded as outliers, assuming that this value was the maximum possible elevation change. First, we removed footprints using the three std dev. (3σ) outlier rejection criteria within each sub-region (Flament and Rémy, 2012), which removed ~5–20% footprints, with the largest removal concentrated in the southeastern HMA and the Himalaya. We then divided each sub-region into elevation zones with 100 m interval and calculated the mean dh (E_i) of each zone. If no data were available for an elevation bin, a zero elevation-change rate was assigned. Volume change was then calculated by multiplying E_i and the glacier area (A_i) of the corresponding zone. The ice density of 850 ± 60 kg m⁻³ was used to translate the volume change into a mass change since the ice density is appropriate for a wide range of conditions and longer-term trends (>5 years) (Huss, 2013), and the final mass-balance change (dB) was obtained using Eqn (3):

(3)$${\rm d}B = \displaystyle{{\rho _{{\rm ice}}} \over {A_{{\rm total}} \times \rho _{{\rm water}}}}\sum\limits_{i = 1}^n {A_i \times E_i} , \;$$

where n is the number of elevation zones, ρ _ice is the glacier ice density, ρ _water is the mass density of water (1000 kg m⁻³) (Huss, 2013) and A _total is the total glacier area in each sub-region.

2.3.5 Mass-balance accuracy assessment

Uncertainties in mass-balance estimates can originate from four sources: residual errors between ICESat-2 DEM and NASADEM after co-registration, uncertainty in penetration depth estimates, glacier area uncertainty and glacier density uncertainty.

Uncertainties of point measurements (i.e. the std dev. of the elevation based on individual grid points) were used to represent lumped uncertainty of geodetic mass balance. In that sense, the uncertainties are implied to be totally correlated (Rolstad and others, 2009). Therefore, decorrelation distance has to be considered. Previous studies have assumed that the ICESat footprints within each 2 km cluster are fully correlated (Moholdt and others, 2010a). In case of few ICESat footprints available, a correlation distance of 5 km was chosen (Moholdt and others, 2010b). We used semi-variogram cloud to simulate correlation distance of ICESat-2 due to its higher resolution and denser distribution (Nuth and others, 2007), and the value obtained was 2 km. The elevation uncertainty (σ _g) that originates from the residual errors between ICESat-2 and NASADEM after co-registration was then calculated by Eqn (5):

(4)$${\rm SE} = \displaystyle{{M_{{\rm ng}}} \over {\sqrt {( N_{{\rm ng}} \times R) /2D} }}, \;$$

(5)$$\sigma _{\rm g} = \sqrt {{\rm S}{\rm E}^2 + {\rm M}{\rm D}^2} , \;$$

where M _ng is the std dev. of the elevation differences in the non-glacier region, N _ng is the number of ICESat-2 footprints in the non-glacier, R is the spatial resolution of the ATL06 product (20 m) and D is the correlation distance (2 km). SE is the standard mean error, and MD is the mean value of the elevation residuals in the non-glacier region.

The penetration correction uncertainty of the NASADEM was calculated from the std dev. of the median penetration depths for all elevation zones below 6500 m (Neelmeijer and others, 2017). We also added the median value of penetration of the non-glacier pixels in calculating penetration correction uncertainty:

(6)$$\sigma _{{\rm penetration}}{\rm} = \sqrt {{\rm std}{( k_i) }^2{\rm} + {\rm penetration}_{{\rm nonglacier}}^2 } , \;$$

where σ _penetration is the penetration correction uncertainty of the NASADEM, and k_i is the median penetration depth estimated for each elevation zone in each sub-region.

The uncertainty of elevation changes in glaciers can be obtained by:

(7)$$\sigma _{\rm h} = \sqrt {\sigma _{\rm g}^2 + \sigma _{{\rm penetration}}^2 }.$$

Regional mass uncertainties can be affected by uncertainty in glacier area σ _A (±10%) and uncertainty in glacier density σ _ice (±60 kg m⁻³). Uncertainty of glacier volume (σ_V) and mass changes (σ_B) will be estimated as follows:

(8)$$\sigma _{\rm V} = \sqrt {{( {\sigma_{\rm h} \times A} ) }^ 2 + {( {\sigma_{\rm A} \times {\rm d}h} ) }^ 2} , \;$$

(9)$$\sigma _{\rm B} = \sqrt {{\left({\displaystyle{{\rho_{{\rm ice}}} \over {\rho_{\rm w} \times A}} \times \sigma_{\rm V}} \right)}^ 2 + {\left({\displaystyle{V \over {\rho_{\rm w} \times A}} \times \sigma_{{\rm ice}}} \right)}^ 2, \;} $$

where σ _V is the uncertainty in glacier volume, σ _h is the uncertainty in elevation changes, A is the glacier area, σ _A is the uncertainty of glacier area, dh is the mean value of glacier elevation-change rate, σ _B is the uncertainty in mass change, ρ _ice is the mean ice density (850 kg m⁻³), ρ _w is the water density, V and σ _V are the glacier volume change and its uncertainty, respectively.

3. Results

3.1 Penetration depths

The radar penetration depths for NASADEM DEM correction were estimated by doubling elevation difference between the X-band SRTM DEM and NASADEM DEM data. Penetration depths in the 15 sub-regions depend strongly on elevation. With increasing elevation, the radar penetration depth into snow and ice for the NASADEM showed a near linear upward trend. Exceptions occurred at the high altitudes of Hissar Alay, Central Himalaya and East Himalaya Mountains, which could be related to the uneven distribution of the X-band SRTM DEM. There were insufficient number of pixels to obtain radar penetration depth estimations at high altitudes. In general, the average penetration depths showed spatial variability in the HMA region (Fig. 2). The lowest penetration depth of ~1.00 m occurred in the South and East Tibetan Mountains, followed by that in the Hindu Kush Mountains at 1.18 m and the Karakoram Mountain at 1.78 m. The largest penetration depth occurred in the Hissar Alay and Qilian Shan mountains, while moderate penetration depths ranging from 2.26 to 4.31 m occurred in the central HMA. Most of the sub-regions exhibited a penetration depth between ±0.1 m for the off-glacier data points.

Fig. 2. NASADEM penetration of the 15 sub-regions of HMA. The bars on the spatial map show the average penetration depth of each region. The line of each subplot donates the penetration depth of each elevation bin (left y-axis), and the shade represents the penetration correction uncertainty of the corresponding bin. The horizontal lines indicate zero line to the penetration estimates. The bars of each subplot display the glacier area of the corresponding bin (right y-axis). Base map was from Esri, USGS, NOAA.

Additionally, the average penetration depths agreed well with estimations by ICESat extrapolation within the error bars in 10 of the 15 regions (Table 1). Regions with larger discrepancies concentrated in the eastern HMA, especially the South and East Tibet Mountains and the Hengduan Shan region.

Table 1. Estimated average NASADEM penetration depth (m) in this study by using method proposed by Jaber and others (2019) and comparisons with estimates based on ICESat extrapolation (error level given is 1 standard error)

Regions with large discrepancies are marked in bold.

3.2 Altitudinal distribution of glacier elevation changes

The altitudinal distributions of the glacier elevation changes in 15 mountain ranges (sub-regions) are shown in Figure 3. The glacier elevation changes in three autumns exhibited similar patterns. For this reason, we merged data from the three years to increase data coverage and statistics. The common pattern of decreasing rates in glacier elevation with increasing altitude can be observed. All the elevation zones in the West Kun Lun Mountain showed a positive trend while thinning occurred in the Himalaya, South and East Tibet, and the Hengduan Shan Mountains, even at the highest elevations.

Fig. 3. Altitudinal distribution of glacier elevation changes in 15 sub-regions in HMA. The blue and orange lines represent the elevation change rates with uncertainty envelopes calculated with footprints of 2019, 2020 and 2021 for each 100 m elevation band (left y-axis). The vertical lines indicate the median glacier elevation, and the horizontal lines indicate zero line to the elevation change. The bars of each subplot display the glacier area of the corresponding bin (right y-axis).

3.3 Glacier mass-balance changes

3.3.1 Glacier mass balance of HMA mountain ranges

The estimated total mass change rate of glaciers in HMA was ~−0.18 ± 0.12 m w.e. a⁻¹ (−17.53 ± 11.36 Gt a⁻¹) from 2000 to 2021. Figure 4a shows the glacier mass changes per year in HMA aggregated over 1° × 1° grids in detail. The south-eastern HMA experienced a relatively significant thinning rate of up to −1 m w.e. a⁻¹, followed by the southern and north-eastern HMA of ~−0.6 m w.e. a⁻¹. A moderate loss rate of −0.1 m w.e. a⁻¹ was seen in the Hissar Alay, Pamir and West Tien Shan Mountains. In contrast, a cluster of positive trends of ~0.2–0.4 m w.e. a⁻¹ occurred near the West Kun Lun Mountains. The South and East Tibet and Hengduan Mountains exhibited the largest mass-balance uncertainty, followed by the Himalaya, while the Tien Shan, Pamirs and eastern Kun Lun Mountains were subjected to less uncertainty (Fig. 4b).

Fig. 4. Glacier mass-balance changes (a) and uncertainty (b) over the HMA for the period from 2000 to 2021. Data are shown on a 1° × 1° grid. The circle color represents the mass-balance variation, and the circle size is scaled according to the glacier area. Cells that do not match the glacier hypsometry were eliminated. Base map was from Esri, USGS, NOAA.

Mass changes varied among regions, with a decrease in glacier mass occurring in most sub-regions (Fig. 5). The Hengduan Shan experienced the highest rate of mass loss at −0.62 ± 0.10 m w.e. a⁻¹ (−2.67 ± 0.45 Gt a⁻¹), followed by the South and East Tibet Mountains with a mass change rate of −0.47 ± 0.11 m w.e. a⁻¹ (−1.85 ± 0.43 Gt a⁻¹). Moderate loss rates occurred in the central and northern HMA, and a negative mass rate ranging from −0.43 to −0.03 m w.e. a⁻¹ was observed. Positive anomalies of 0.23 ± 0.13 m w.e. a⁻¹ only occurred in the West Kun Lun Mountain.

Fig. 5. Specific glacier mass balance (m w.e. a⁻¹) for the period from 2000 to 2021, aggregated over the RGI mountain ranges. Sub-regions from HiMAP (Bolch and others, 2019) and Kääb and others (2015) can be found in Figure S2.

3.3.2 Glacier mass budget of river basins

The spatial distribution of all river basins in HMA following the basin boundaries of Shean and others (2020) is shown in Figure 6. The most severe mass loss took place in the exorheic Brahmaputra, Ganges and Indus Basins and in the endorheic Ili and Amu Darya Basins. Only Tarim Basin showed a mass gain of +2.91 ± 1.82 Gt a⁻¹.

Fig. 6. Glacier excess melt runoff for the major river basins in HMA during 2000–2021. Colors of polygons denote the glacier runoff in Gt a⁻¹, the heights of bars represent the sizes of the glacier in each basin, and the colors of bars represent the mass balance in units of m w.e. a⁻¹. The red line indicates the dividing line of the exorheic and endorheic basins. Total mass losses were indicated in italics. Base map was from Esri, USGS, NOAA.

To estimate the contribution to sea level rise (SLR), we assume an ocean area of 3.62 × 10⁸ km² (Jakob and others, 2021). The excess discharge is 0 for glaciers with a balanced or positive mass budget (Brun and others, 2017), so we excluded Tarim Basin for further analysis. We considered two conditions for SLR contribution: from all basins and only from exorheic basins. The total potential contribution to SLR from all basins can be equivalated to a SLR of 0.055 ± 0.031 mm a⁻¹. The total contribution from exorheic basins is equivalent to 0.043 ± 0.015 mm a⁻¹.

4. Discussion

4.1 Penetration uncertainties/comparison

The offsets over debris-covered ice and off-glacier terrain are negative in few elevation bands of some sub-regions meaning that we observe a negative penetration. One main reason is the insufficiency of valid samples in the X-band SRTM data in these elevation zones. Another possible reason may be attributed to a biased correction in bedrock areas, due to inaccurate co-registration in high-relief regions (Li and others, 2021). This kind of trend can also be seen in the penetration depths estimated by ICESat backpropagation (Kääb and others, 2012). In addition, the penetration depth is small at the high altitudes of Hissar Alay, Central Himalaya and East Himalaya Mountains, which may be connected to the spatial difference in the glacial distribution.

Secondly, the estimation of penetration depth varies with different estimation methods. Subtracting DEMs from field measurements or optical DEMs can provide accurate penetration depths, but is generally difficult to apply to the entire HMA region since there were no field or optical data available when the NASADEM data were acquired. Therefore, penetration depth can only be estimated through different assumptions. Methods proposed by Jaber and others (2019) and Zhou and others (2019) rely on relative penetration depth differences between the X-band and C-band SRTM or NASADEM. We compare the relative penetration depth differences with previous studies (Table 2). The penetration depths we estimated agree well with those by Li and others (2021) except for the southern and south-eastern HMA. Li and others (2021) used the void-filled SRTM, which may contain fill-in pixels that are composed of different timestamps or from optical data. The void-filled portions are ~20% in the south-eastern HMA, thus the estimates derived from X-band SRTM and NASADEM are more consistent with those obtained by the non-void-filled SRTM (Zhou and others, 2018, 2019). In addition, when we assumed that the average X-band penetration depth of 4 m and the percentage of dry snow area of 40% (i.e. the average SRTM X-band penetration of 1.6 m) (Zhou and others, 2019), the estimates by doubling the relative penetration differences generally agreed with them. However, it is difficult to quantify the percentage of dry snow area for each region in HMA, where each region has its own climatology.

Table 2. Estimated average NASADEM/C-band SRTM penetration depth (m) relative to the X-band SRTM

This study used non-void-filled NASADEM and X-band SRTM to estimate penetration differences. Li and others (2021) and Gardelle and others (2013) used X-band SRTM and the void-filled C-band SRTM, while Zhou and others (2018) masked out void-filled regions in the C-band SRTM to estimate penetration differences.

Furthermore, we compared the average penetration depth estimates of C-band radar proposed by Kääb and others (2012). Compared with ICESat extrapolation, penetration estimates derived from doubling the relative penetration differences agreed well within the error bars in 10 of 15 regions (Table 1). Area with larger discrepancies located in the east HMA, especially the South and East Tibet Mountains and the Hengduan Shan mountains. ICESat only contained ~1000 on-glacier points in each measurement cycle in each sub-region, whereas the SRTM-X DEM can cover 10–40% glaciated areas. Hence, the difference mainly originated from the different sampled areas and elevations. The penetration depth of the C-band radar estimated by ICESat linear extrapolation was ~8–10 m in the south-eastern HMA (in both this study and that by Kääb and others (2012)), which could have been overestimated, as this depth only occurs when the snow density is ~200 kg m⁻³ (Zhou and others, 2018). On the contrary, the penetration depth estimates by ICESat extrapolation are smaller in East Kun Lun and Qilian Shan mountains. Estimates based on ICESat extrapolation may be suitable for areas with relatively consistent interannual elevation variations. Kääb and others (2012) pointed out that the C-band penetration is also expected to be deeper on winter-accumulated glaciers than on summer-accumulated ones. The estimates derived from doubling the relative penetration differences are consistent with this assumption given the wider spatial coverage of X-band SRTM.

We finally considered the method of assuming that SRTM radar captured the 1999 autumn elevation profiles due to a reflective ice layer (Berthier and others, 2006) so that we could use ICESat-2 acquired during autumn season to calculate glacier mass balance. However, such a reflective-ice-layer assumption may work for lower glacier areas but may not be true for the upper accumulation areas in HMA where a distinct firn layer may not exist. C-band radar penetration is deeper in the accumulation areas in HMA, and the assumption may lead to much underestimated mass balance (Table S1).

For errors associated with radar penetration, we have considered the large uncertainties by using conservative error bars with the std dev. of the median penetration depths. Uncertainty in most regions exceeded 50% of the estimated region-mean penetration depths. As one of the significant sources of mass-balance uncertainty in using SAR-based DEMs, we provided a reference data for estimating regional glacier mass balance in HMA based on NASADEM, which may favor the increase of mass-balance accuracy.

4.2 Impact of NASADEM voids

The larger uncertainty in higher altitude regions is associated with the small number of laser footprints collected in glaciated areas. The larger uncertainty may also be attributed to the larger errors in the NASADEM on steeper slopes (Chen and others, 2022) and the relatively larger percentage of voids compared to lower-slope regions.

The exclusion of void-filled areas in the NASADEM may introduce bias in the sampling (Kääb and others, 2012). We matched the glacier hypsometry of ICESat-2 and the void-filled NASADEM so that the area sampled by ICESat-2 is sufficiently representative of the entire region. Mass balance was calculated by comparing ICESat-2 and non-void-filled NASADEM because the void-filled portions may come from completely different data sources. The hypsometry difference between the void-filled and non-void-filled NASADEM may introduce additional uncertainty during mass-balance calculation. Therefore, we calculated the percentage of NASADEM voids in each sub-region. Different regions have different percentage of voids, and the regions that are most affected by voids are the Central Himalaya, followed by Hengduan Shan and South and East Tibet Mountains. Subsequently, we compared the hypsometry of the void-filled NASADEM and the non-void-filled NASADEM since the void-filled portions can serve as an elevation-range reference, and the two hypsometry results agreed well in 14 of 15 regions, except the South and East Tibet Mountain. Nonetheless, due to lack of measurements in the high altitudinal zones, it is difficult to apply the hypsometry filter; thus, we add an additional 10% of mass-balance uncertainty for the South and East Tibet Mountain in the final results. The voids are mostly distributed in the upper elevation zones, which are often above the mass equilibrium line, and any missing elevation measurements may cause the estimated mass balance toward a more negative value.

4.3 Comparisons with previous studies

Glacier mass balance in this study was generally consistent with other assessments. We found that the mass loss over the entire HMA during 2000–2021 is most similar to those derived from ASTER DEMs (Table 3). The sub-regional mass balances all showed the same spatial patterns with mass balance aggregated over the RGI regions (Table 4), Kääb boundary (Table S2) and HiMAP boundary (Bolch and others, 2019) (Table S3). The glacier mass-balance estimates for the river basins in HMA also conformed to previous studies (Brun and others, 2017; Shean and others, 2020) (Table S4), with the only difference occurring in the Tarim Basin. Results revealed a mass gain of +2.91 ± 1.82 Gt a⁻¹ in Tarim from 2000 to 2021, while a mass loss of −0.87 ± 0.71 Gt a⁻¹ was observed from 2000 to 2018 (Shean and others, 2020).

Table 3. Previously published mass-balance estimates for HMA

^a Represents data obtained by elevation change by assuming the snow density of 850 ± 60 kg m⁻³, which is the same as the density assumption in other studies. NA indicates data not available in the reference.

Table 4. Region-wide mass balance compared with previous studies aggregated over the RGI boundary

Linear-fit of ICESat/ICESat-2 height differences between NASADEM was obtained by following the method of Wang and others (2021). The result of ASTER DEMs during 2000–2019 was derived from the published elevation-change map by Hugonnet and others (2021).

Mass-balance estimation by Hugonnet and others (2021) is the recent research without the influence of radar penetration depth, the glacier mass balance derived from the elevation-change map can serve as a validation of absolute differencing of ICESat-2 and NASADEM. Wang and others (2021) used relative elevation differences between ICESat/ICESat-2 and void-filled SRTM, and performed a linear fit to the relative elevation differences to obtain glacier elevation change. Void-filled portions may bias the mass-balance estimates, so we recalculated the glacier mass balance using the non-void-filled NASADEM to make the comparison more robust. The glacier mass balances derived from absolute differencing of ICESat-2 and NASADEM tend to be less negative than those from the linear-fit (Table 4 and Fig. S3), although mostly within the error envelop, and we attributed the difference to the fewer ICESat footprints. The results show more consistency with the ASTER-based estimates (Table 4), but the elevation changes difference between the NASADEM and ASTER DEMs appear to be systematically biased. Taking the East Himalaya Mountain as an example, the glacier elevation changes in the low altitudes are more negative from the absolute differencing of ICESat-2 and NASADEM (Fig. S4), which might be associated with the penetration depth estimates tending to be slightly overestimated, but we have considered this penetration correction uncertainty in the accuracy assessment.

5. Conclusions

Radar penetration into snow/ice is essential for accurately estimating glacier mass balance based on the radar-derived NASADEM. To resolve the problem, we corrected the NASADEM DEM by doubling elevation differences between the SRTM X-band DEM and NASADEM. The glacier elevation and mass changes from 2000 to 2021 over the HMA region were then estimated by comparing the elevation differences of ICESat-2 and penetration depth corrected NASADEM.

We concluded that: (1) the penetration depths show an elevation-dependent pattern in which the penetration depth into the snow/ice by the C-band radar showed an upward trend with increasing elevation. Relatively larger average penetration depths occurred in the Himalaya and Hissar Alay mountain range, while smaller penetration depths were found in the South and East Tibet and Hengduan Mountains; (2) from 2000 to 2021, the glacier mass balance of HMA was −0.18 ± 0.12 m w.e. a⁻¹ (−17.53 ± 11.36 Gt a⁻¹), where the greatest mass loss occurred in the south-eastern HMA. Compared with previous studies, the West Kun Lun Mountain experienced a more conspicuous mass gain at a rate of 0.19 ± 0.13 m w.e. a⁻¹; (3) the HMA glaciers have experienced mass losses in all basins but the Tarim basin, with the largest loss of −4.99 ± 2.05 Gt a⁻¹ in the exorheic Indus basin, and the excessive glacier melt runoff of HMA was equivalent to a SLR of 0.055 ± 0.031 mm a⁻¹.

The glacier mass balance during 2000–2021 yielded the same spatial patterns as previous studies, and the conformity of DEM-derived mass-balance estimates confirm the necessity of penetration depth correction for the NASADEM. ICESat-2 can be used to estimate glacier mass balance by comparing it with DEMs directly for regions where ICESat-2 data have a reasonable altitudinal distribution.