A Systematic, Regional Assessment of High Mountain Asia Glacier Mass Balance

High-mountain Asia (HMA) constitutes the largest glacierized region outside of the Earth's polar regions. Although available observations are limited, long-term records indicate sustained HMA glacier mass loss since ~1850, with accelerated loss in recent decades. Recent satellite data capture the spatial variability of this mass loss, but spatial resolution is coarse and some estimates for regional and HMA-wide mass loss disagree. To address these issues, we generated 5,797 high-resolution digital elevation models (DEMs) from available sub-meter commercial stereo imagery (DigitalGlobe WorldView-1/2/3 and GeoEye-1) acquired over HMA glaciers from 2007 to 2018 (primarily 2013–2017). We also reprocessed 28,278 ASTER DEMs over HMA from 2000 to 2018. We combined these observations to generate robust elevation change trend maps and geodetic mass balance estimates for 99% of HMA glaciers between 2000 and 2018. We estimate total HMA glacier mass change of −19.0 ± 2.5 Gt yr−1 (−0.19 ± 0.03 m w.e. yr−1). We document the spatial pattern of HMA glacier mass change with unprecedented detail, and present aggregated estimates for HMA glacierized sub-regions and hydrologic basins. Our results offer improved estimates for the HMA contribution to global sea level rise in recent decades with total cumulative sea-level rise contribution of ~0.7 mm from exorheic basins between 2000 and 2018. We estimate that the range of excess glacier meltwater runoff due to negative glacier mass balance in each basin constitutes ~12–53% of the total basin-specific glacier meltwater runoff. These results can be used for calibration and validation of glacier mass balance models, satellite gravimetry observations, and hydrologic models needed for present and future water resource management.


INTRODUCTION
Glaciers are sensitive climate indicators that primarily respond to interannual changes in temperature and precipitation (e.g., Bertrand et al., 2012;Harrison, 2013). They constitute important seasonal and long-term hydrologic reservoirs, providing water for hydropower, agriculture, and municipal use (Guido et al., 2016;Ragettli et al., 2016;Milner et al., 2017;Pritchard, 2019) Glaciers can also be a significant natural hazard, especially for regions subject to catastrophic glacier outburst floods (Harrison et al., 2017;Haritashya et al., 2018;Allen et al., 2019).
Many studies have documented changes in HMA glaciers, and their total contribution to global sea level rise. Many earlier studies of HMA glaciers leveraged "traditional" glaciological measurements, while more recent efforts rely on geodetic remote sensing observations, including satellite laser altimetry (e.g., Ice, Cloud, and Land Elevation Satellite [ICESat]), satellite gravimetry (e.g., Gravity Recovery and Climate Experiment [GRACE]), and digital elevation model (DEM) differencing (e.g., DEMs from Pléiades, Satellite Pour l'Observation de la Terre [SPOT], and Advanced Spaceborne Thermal Emission and Reflection Radiometer [ASTER]). These methods have inherent differences in sampling strategy, resolution, and sensitivity, which can lead to discrepancies in results. For detailed reviews of past work, we refer the reader to Bolch et al. (2012), Farinotti et al. (2015), Kääb et al. (2015), Brun et al. (2017), Azam et al. (2018), and Bolch et al. (2019).
Long-term records of glacier length and area suggest that Himalayan glaciers have been retreating since ∼1850 (Bolch et al., 2012). Since the 1960s, traditional and geodetic mass balance observations document mass loss rates of −0.3 to −0.7 m w.e. yr −1 , with greater loss since the mid-1990s (Bolch et al., 2012). Gardner et al. (2013) reported 2003-2009 glacier mass loss rates of −0.24 ± 0.11 m w.e. yr −1 (−29 ± 13 Gt yr −1 ) from ICESat altimetry, −0.16 ± 0.17 m w.e. yr −1 (-19 ± 20 Gt yr −1 ) from GRACE observations, and −0.72 ± 0.22 m w.e. yr −1 (−86 ± 26 Gt yr −1 ) from extrapolating traditional mass balance observations to the full region. The larger mass loss estimate from the traditional data was attributed to a sampling bias of glaciers in accessible areas at lower elevations (Gardner et al., 2013). Kääb et al. (2015) used ICESat data and Shuttle Radar Topography Mission (SRTM) DEM to estimate glacier surface elevation change for 1 • × 1 • cells, with an estimate for total HMA glacier mass loss rate of −0.37 ± 0.1 m yr −1 from 2003 to 2008. While the ICESat and GRACE data provided the first systematic estimates of mass change for the entire HMA region, the sparse sampling of ICESat tracks, coarse resolution of GRACE mascons, and potential mixing of glacier mass change signals with terrestrial water storage and groundwater depletion contributed to high uncertainty estimates. Brun et al. (2017) improved these estimates using time series of ASTER DEMs between 2000 and 2016 to estimate total HMA glacier mass loss of −0.18 ± 0.04 m w.e. yr −1 (−16.3 ± 3.5 Gt yr −1 ). The observed spatial variability across 1 • × 1 • cells was generally consistent with past work: mass change in the Himalayas was negative (−0.33 ± 0.20 to −0.42 ± 0.20 m w.e. yr −1 ), while the Karakoram was near zero (−0.03 ± 0.07 m w.e. yr −1 ), and Kunlun Shan was positive (+0.14 ± 0.08 m w.e. yr −1 ). The near-zero mass change in the Karakoram (a.k.a. the "Karakoram anomaly") is consistent with many other studies (Gardelle et al., , 2013Rankl and Braun, 2016). This situation has persisted for at least ∼50 years (Bolch et al., 2017), and has been attributed to increased winter precipitation (Archer and Fowler, 2004) and decreased mean summer temperature (Fowler and Archer, 2006) since the 1960s (Bashir et al., 2017;Forsythe et al., 2017). Glaciers in the Western Kunlun Shan have also been near a balanced state since the 1970s despite an increasing trend in mean summer and annual temperatures (Wang et al., 2018). de Kok et al. (2018) hypothesized that the positive mass balance in the Kunlun sub-region may be related to increased irrigation in the Tarim basin, which caused an increase in summer snowfall and cloud cover over the Kunlun Shan, with associated decrease in net incoming shortwave radiation.
To assess glacier evolution under future climate scenarios, geodetic glacier mass balance observations can be used to calibrate and validate glacier evolution model projections. These models estimate ∼36-87% mass loss for HMA glaciers by 2100, depending on the glacier evolution model, climate forcing, and representative concentration pathway (RCP) scenario (Kraaijenbrink et al., 2017;Hock et al., 2019). This loss will lead to a significant decrease in seasonal glacier meltwater runoff, with important implications for downstream stakeholders (Huss and Hock, 2018).
While glacier mass balance intercomparison efforts have improved considerably, they are often incomplete, or integrate datasets spanning different time periods with different spatial extent and resolution. Spatially continuous, systematic highresolution observations are needed to capture both the spatiotemporal evolution of glacier mass balance and the processes that influence glacier mass balance (e.g., debris cover evolution, mass redistribution during surging) on a regional scale. To address these issues, we generated a new dataset of high-resolution DEMs from DigitalGlobe/Maxar stereo satellite imagery with unprecedented coverage and accuracy. We also reprocessed the archive of daytime ASTER stereo FIGURE 1 | Composite products for WorldView/GeoEye (A-C) and ASTER (D-F) DEMs after co-registration and quality control. Rows show per-pixel weighted mean, count, and normalized median absolute deviation (NMAD). Note difference in color ramp for count maps in B and E.
imagery acquired between 2000 and 2018. We integrated these observations to derive robust elevation change trends and mass balance estimates for ∼99% of the glacierized area in HMA, with no minimum glacier area threshold. Our results will provide important calibration and validation for mass balance models needed to improve estimates of future HMA glacier change and understand associated downstream impacts in the region.

WorldView/GeoEye (WV/GE) DEMs
We processed all available Level-1B (L1B) panchromatic stereo imagery over HMA from the DigitalGlobe/Maxar archive through late 2017. This included imagery collected by WorldView-1 (0.50 m ground sample distance [GSD] at nadir, beginning in late 2008 for HMA), WorldView-2 (0.46 m GSD, mid 2010), WorldView-3 (0.31 m GSD, late 2014), GeoEye-1 (0.41 m GSD, late 2009), and Quickbird-2 (0.65 m GSD, 2002). A total of 2,562 in-track and 3,235 cross-track stereo pairs were processed using the NASA Ames Stereo Pipeline (ASP) (Shean et al., 2016;Beyer et al., 2018) v2.6.0 (Beyer et al., 2017) following the methods outlined by Shean et al. (2016) (Figure 1). Cross-track pairs were formed by monoscopic images with acquisition time separation of <7 days, convergence angle of 10-70 • , minimum intersection width and height of 6 × 6 km, and cloud cover of <25%. We used the void-filled SRTM-GL1 (Farr et al., 2007) as a seed DEM for initial orthorectification, enabling efficient stereo correlation (Shean et al., 2016). We also used an absolute elevation difference filter of ±200 m relative to the SRTM-GL1 products to remove spurious pixels from the output DEMs. Output DEMs were posted at 2, 8, and 32 m in the appropriate universal transverse mercator (UTM) zone, with elevation values relative to the WGS84 ellipsoid. These data products, detailed processing information, and metadata are available from the National Snow and Ice Data Center (NSIDC) (Shean, 2017a,b,c).

ASTER DEMs
We generated DEMs from 28278 Level-1A (L1A) stereo images collected by the ASTER Visible and Near Infrared (VNIR) instrument (Fujisada, 1998) between June 2000 and June 2018 (Figure 1). We queried the NASA EarthData archive for daytime stereopairs (panchromatic bands 3N and 3B, 15 m GSD) with cloud cover <50%. The ASTER pairs were processed with ASP v2.6.0 using the void-filled SRTM-GL1 product as a seed DEM for initial orthorectification. We used the ASP semi-global matching (SGM) correlator (Hirschmuller, 2008) with default parameters (7 × 7 pixel window), which can offer improved results over the default ASP correlation routines for scenes with limited image resolution and/or texture (i.e., fresh snow-cover and cloudshadowed areas). We used ASP's default SGM disparity map filters (3 × 3 pixel median filter and a 3 × 3 pixel "texture-aware smoothing filter" with scaling factor 0.13) to remove residual artifacts. Finally, isolated clusters of <50 pixels surrounded by missing data ("islands") were removed from the disparity map. We used the rational polynomial coefficient (RPC) model bundled with each L1A scene for stereo triangulation, and culled outliers with large triangulation error (exceeding 3 × 75thpercentile) from the resulting point cloud. Output DEMs were posted at 30 m with elevations relative to the WGS84 ellipsoid. Finally, we performed a 2-pixel erosion of the outer DEM boundaries to remove any residual edge artifacts and filtered the resulting DEMs using an absolute elevation difference filter, removing any DEM pixels with ±100 m offset relative to the voidfilled SRTM-GL1. The ASTER DEMs contain more noise than the WorldView/GeoEye DEMs, so we chose a more conservative threshold to remove spurious DEM pixels (e.g., artifacts near cloud margins), while also preserving the elevation change signals of interest over glacier surfaces.

TanDEM-X Global DEM
We used the publicly available TanDEM-X 1-arcsec (90 m) global DEM as a reference basemap for the region (Rizzoli et al., 2017;Wessel et al., 2018). While relatively coarse, this product offers excellent horizontal and vertical (∼3.5 m absolute and <2 m relative) accuracy. We used the auxiliary products bundled with each DEM tile to mask artifacts/errors observed over some of the relatively steep terrain in HMA, including the theoretical height error map (HEM ≤ 1.5 m), water mask (33 ≤ WAM ≤ 127), and consistency mask (COM = {0,1,2}). Qualitative inspection of DEM tiles around the region suggest that these filters removed most spurious pixels. Further details and code for TanDEM-X processing is available in the tandemx repository (Shean, 2019). The resulting products were used as a reference DEM for co-registration.

Co-registration
All individual 30-m ASTER and 8-m WV/GE DEMs were co-registered to remove any horizontal and vertical offsets relative to the filtered reference TanDEM-X DEM. We identified "static" control surfaces that we assume did not change between the DEM acquisition timestamp and the reference DEM acquisition timestamp (∼2011-2014 for TanDEM-X products). Both input DEMs were clipped to a common intersection and resampled to 30 m via bicubic interpolation. Glacier surfaces were masked with the RGI v6.0 polygons (RGI Consortium, 2017), and outlier elevation difference values over remaining static surfaces were identified and removed. We then used an iterative implementation of the Nuth and Kääb (2011) method, with robust filtering and outlier removal. If the resulting translation exceeded a 200 m threshold or if more than 30 iterations were required, the DEM was not used for subsequent analysis. For all other DEMs, the resulting horizontal translation was applied to each DEM and the median vertical bias was removed. We did not apply any corrections to residual differences (e.g., 5th order polynomials used by Brun et al., 2017), to avoid introducing additional elevation error over glacier surfaces due to overfitting of limited static control surfaces with poor spatial distribution. The dem_align.py script in the demcoreg package  outlines the full co-registration workflow.
We then analyzed the full set of co-registered ASTER and WV/GE DEMs to identify and remove any remaining problematic and/or low-quality DEMs. We estimated population statistics for both DEM sources using metrics from all individual DEMs before and after co-registration (Figure 2). We used these statistics to identify outliers and removed hundreds of DEMs with anomalously high residual bias and spread (Figures 2B,D). The dem_align_post.ipynb notebook in the demcoreg package  outlines the full post-coregistration outlier removal workflow. Figure 2 shows statistics for the final set of co-registered and validated DEMs.

Regional DEM Composites and Mosaics
We generated tiled DEM composite and mosaic products for the HMA region (Shean, 2017a) using the coregistered, filtered DEMs (Figure 1). All products use a custom Albers Equal Area projection that encompasses the RGI glacier polygons (PROJ4 string: +proj=aea +lat_1=25 +lat_2=47 +lat_0=36+lon _0=85 +x_0=0 +y_0=0 +ellps=WGS84 +datum= WGS84 +units=m +no_defs). We generated separate products for WV/GE and ASTER datasets, as well as a combined, blended, gap-free mosaic. Mosaic grid spacing was 8 m for WV/GE products and 30 m for ASTER products, with 100 × 100 km tiles. We also generated 1/3 arcsecond (∼10 m) and 1-arcsecond (∼30 m) resolution products using 1 • × 1 • tiles to conform with previous standards for global DEM mosaics (e.g., SRTM). The ASP dem_mosaic utility was used to generate seamless weighted-average composites and per-pixel composites of count, median, standard deviation, and normalized median FIGURE 2 | Temporal coverage and error metrics for WorldView/GeoEye (A,B) and ASTER (C,D) DEMs before and after co-registration. Each point shows the median elevation offset between an individual DEM and the TanDEM-X Global DEM reference over surfaces assumed to be static. Error bars show the 16-84% spread of these elevation offset values for each DEM. The final set of DEMs used for subsequent analysis and trend-fitting are shown with orange points/bars ("inliers"). Note the apparent drift in median offsets throughout the ASTER mission (C), and zero bias after co-registration.
absolute deviation (NMAD) using all valid DEM elevations at each pixel. We also generated timestamped mosaics using "forward" (most recent DEMs on top) and "reverse" (earliest DEM on top) per-pixel sort order, with associated products showing decimal year of elevation at each pixel. The code for composite and mosaic generation is available in the make_mos.sh and dem_mosaic_validtiles.py scripts in the gmbtools package (Shean and Bhushan, 2019).

Glacier DEM Stack Generation and Elevation Change Trend
We generated "stacks" of combined ASTER and WV/GE DEMs for each glacier polygon (plus a 1 km buffer) at 30 m posting. We computed linear elevation change trend ( dh dt , with units in m yr −1 ) (Figure 3, Figure S1) for each pixel with >5 DEMs and >5 years between the first and last DEM timestamp. Several methods for the trend fitting were considered ( Figure S1), including ordinary least-squares, random sample consensus (RANSAC), and the Theil-Sen estimator (scikit-learn implementation: https://scikit-learn.org/stable/modules/linear_ model.html#theil-sen-regression). We reviewed maps of trend and residual values for a large sample of glaciers across the region, and confirmed that the robust Theil-Sen estimator offered superior output quality for the noisy ASTER+WV/GE DEM stacks ( Figure S1).
The output elevation change trend map for each glacier was then filtered to remove artifacts. For each glacier (plus the buffered extent), a stack per-pixel median grid was created. From this grid, minimum (0.01 percentile) and maximum (99.99 percentile) stack elevation values were determined. These elevation values typically sampled the bottom of a valley floor, and the top of the highest peak in the scene. Then, these minimum and maximum elevation values were padded by and additional −150 and +150 m, respectively. Per-pixel elevation change trends were extrapolated to obtain surface elevation maps for May 31, 2000 and May 31, 2018. If the extrapolated elevations fell outside of the padded bounds, the pixel was masked. A 1-pixel erosion step removed trend values near data gaps, and a 3 × 3pixel median filter followed by a 3 × 3-pixel Gaussian filter were used to remove residual outliers and smooth the output trend maps (Figure 3). Details of this workflow can be found in the stack_interp.py script in the gmbtools package (Shean and Bhushan, 2019).

Glacier Mass Balance
We computed glacier mass balance and mass balance uncertainty at multiple spatial scales: individual glacier polygons, tessellated hexagonal grid cells, glacierized regions, river basins, and the full HMA domain. We use d dt notation for pixel values, and notation for spatially aggregated values, with both representing the annualized rate of change for the period from 2000 to 2018.

Individual Glacier Mass Balance
We used the elevation trend maps to estimate geodetic mass balance for glacier polygons with a minimum of 85% coverage. In practice, 95% of the glaciers in the study area had >97% coverage, eliminating the need to rely on gap-filling approaches (e.g., McNabb et al., 2019). We computed dh dt statistics for each glacier polygon, computed mean glacier elevation change h, and multiplied by polygon area A to estimate total individual glacier volume change V. We used a standard bulk density estimate ρ = 850 kg m −3 (Huss, 2013) to convert this volume change to individual glacier total mass balance M: In contrast to many previous studies, we did not impose a minimum glacier area threshold, and consider all 95536 HMA glacier polygons in RGI v6.0.

Individual Glacier Mass Balance Uncertainty
We calculated individual glacier mass balance uncertainty using approaches similar to those outlined in recent geodetic analyses with similar DEM sources (e.g., Fischer et al., 2015;Berthier et al., 2016;Brun et al., 2017;Menounos et al., 2019).
To estimate glacier elevation change error σ h , we calculated statistics for per-pixel dh dt trend values within a 500 m buffer around each glacier polygon, assuming dh dt values should be 0 for these "static" or "stable" surfaces. We then combined "random" and "systematic" error components to estimate total error.
To estimate random error σ dh dt random (i.e., spread of dh dt noise), we considered both the NMAD and standard deviation of dh dt values over static surfaces. We used the Rolstad et al. (2009) "rule of thumb" approach, and a decorrelation length scale L of 500 m between pixels in the dh dt maps. We scaled the σ dh dt random values for glaciers with A ≥ πL 2 (∼0.78 km 2 ) by an area-dependent coefficient π L 2 5A (Rolstad et al., 2009), which accounts for the fact that larger areas dominated by uncorrelated random error will approach 0 mean error. To account for any residual systematic error σ dh dt systematic (i.e., local dh dt bias), we computed the mean of dh dt values over the same static surfaces. We then combined this σ dh dt systematic value with the scaled σ dh dt random value (computed from the standard deviation of dh dt values) in quadrature to estimate total root-mean-square error (RMSE): Many previous studies report σ h as the scaled σ dh dt random component from the NMAD of dh dt values. These NMAD-only error estimates are more robust to outliers and tend to be relatively small. The combined RMSE estimates from Equation (2) are larger and more conservative, and we use these values for σ h in this study.
Assuming that the three error components (σ h , σ A , σ ρ ) are independent and uncorrelated, the total mass change uncertainty for each glacier polygon is: (using Equation 10 in Joint Committee for Guides in Metrology, 2008). An alternative form of this equation can be written in terms of normalized, dimensionless fractional uncertainty: which is convenient given that our area and density uncertainty are provided as percentages. While Equation (4) appears to be undefined for | h| = 0, if we distribute the | M| term, then as The mb_parallel.py script in the gmbtools package  contains the detailed workflow used for individual glacier mass balance and uncertainty estimation.
As defined above, we assume σ A A = 0.10 and σ ρ ρ = 0.071, so the minimum glacier mass balance uncertainty will be ∼12.2% when σ h | h| is 0. In practice, the σ h | h| term tends to dominate the mass balance uncertainty for individual glaciers, but aggregation over larger areas reduces its influence on total mass balance uncertainty.

Aggregated Glacier Mass Balance
We aggregated the individual glacier mass balance results for the full HMA region and different sub-regions. We computed centroids for each glacier polygon and performed spatial joins to compile statistics for the centroids that fell within larger sub-region polygons. The sub-region aggregation was performed for the commonly used glacierized regions of Kääb et al. (2015) and the HIMAP project , the major hydrologic basins , updated to 6-min spatial resolution), and the Goddard GRACE mascon boundaries (Loomis et al., 2019). We modified three of the Kääb et al. (2015) regions to remove references to country names ("Bhutan" → "East Himalaya, " "East Nepal" → "Central Himalaya" and "West Nepal" → "West Himalaya").
For each sub-region, we computed the sum of glacier mass balance ( M) for all n glac glaciers in cubic meters water equivalent (m 3 w.e. yr −1 ) and gigatons per year (Gt yr −1 ) and divided by aggregated glacier area to obtain mean specific mass balance (m w.e. yr −1 ). We converted Gt to mm sea level rise equivalent using a global ocean area estimate of 3.625 × 10 8 km 2 (Cogley et al., 2011), and considered total contributions from all glaciers and then only those glaciers within exorheic basins.

Aggregated Glacier Mass Balance Uncertainty
Estimating total mass balance uncertainty for larger aggregation areas is more complex than for individual glaciers. Many previous studies compute h for all valid dh dt pixels over glaciers and σ h for all dh dt pixels over stable surfaces in the area of interest (e.g., river basin, 1 • × 1 • cells), often aggregating statistics for elevation bins and scaling using observed glacier hypsometry. The uncertainty is sometimes scaled to account for relatively short autocorrelation length scales of dh dt pixel values over stable areas (e.g., ∼500 m), but rarely for longer autocorrelation length scales (e.g., 10s of km). It is also possible that glacier polygons could be split across aggregation boundaries (e.g., glacier accumulation area in one 1 • × 1 • cell and ablation area in adjacent cell), which will bias aggregated mass balance estimates.
To avoid these issues, we aggregated mass balance uncertainty values for individual glacier polygons, not dh dt pixels. Initially, we used the individual glacier mass balance uncertainty values (σ M ) from Equation (3). Combining these values in quadrature resulted in unrealistically low aggregated uncertainty estimates for the full-HMA region (0.3 Gt yr −1 or 1.7%) due to the large glacier sample size (n glac ≈ 95, 000 for full-HMA region). This approach assumes that the σ M values for all glacier polygons are independent and completely uncorrelated, which is incorrect for several reasons: the density uncertainty σ ρ is constant for all glaciers (perfect spatial autocorrelation), area uncertainty percentage σ A A is also constant for all glaciers, and there is inevitably some large-scale (∼10s of km) spatial autocorrelation of dh dt pixel values over stable terrain that is not accounted for using the scaling described in section 2.5.2. This large-scale σ h spatial autocorrelation is controlled by the input DEM sample counts (Figure 1), DEM dimensions (∼ 60 × 60 km for ASTER), DEM artifact length scales, and slope-and aspect-dependent residual co-registration errors.
If we instead assume that the σ M values are perfectly correlated and compute the sum of all σ M during aggregation, we obtain unrealistically high aggregated uncertainty estimates for the full-HMA region (21.5 Gt yr −1 or 113.4%). A more appropriate estimate of aggregated mass balance uncertainty falls between these two approaches, as the errors are correlated, but only over a limited range of length scales (e.g., Rolstad et al., 2009;Anderson, 2019). To estimate these length scales, we calculated the range of autocorrelation for σ M values using experimental semivariograms with a lag distance equal to the median of the nearest neighbor distances between glacier polygon centroids (∼786 m). We fit a spherical semivariogram to the experimental semivariogram using ordinary least squares, and obtained an autocorrelation length scale estimate of ∼42 km for σ M using the RMSE error metric (Equation 2) and ∼27 km for σ M using the NMAD error metric. See Figure S2 for details and the Mass_balance_cor_working.ipynb notebook in the raster_geostatistics repository (Bhushan, 2020) for additional details.
This approach using total mass balance uncertainty for each glacier (σ M ), however, fails to consider variable spatial autocorrelation length scales for the individual components of σ M in Equation (4). A better strategy is to isolate and consider the spatial autocorrelation of only the σ h values, as the σ ρ and σ A A values are constant across the region, and their inclusion will increase apparent spatial autocorrelation of the σ M values. We repeated our semivariogram analysis for the σ h values and found an autocorrelation length scale of ∼32 km for the σ h RMSE and ∼24 km for the σ h NMAD error metrics ( Figure S2) (dh_dt_sigma_error.ipynb in Bhushan, 2020).
Based on the range of influence for the σ h RMSE values, we defined a tessellated grid of hexagonal cells across the HMA domain to represent independent spatial samples with "radius" of 32 km. The width between opposite sides of these hexagon cells is 55 km, which is roughly equivalent to a quarter-degree cell. We note that aggregation for cells defined in units of decimal degrees (e.g., 1 • × 1 • cells used by many previous HMA-wide inventories) can be problematic, as the cells cover variable total area at different latitudes. For example, a 1 • × 1 • cell at 25 • N covers ∼11,182 km 2 while a 1 • × 1 • cell at 47 • N covers ∼8,455 km 2 , an areal difference of ∼24%. To avoid potential sampling bias, we defined these equal-area hex cells in our custom equalarea projection (section 2.4.2).
We computed the mean glacier elevation change for each hex cell ( h cell ) using the sum of all individual glacier volume change ( V = h·A) estimates for the n glac glacier polygon centroids that fell within each cell: where A cell is the total area of the n glac glacier polygons in the cell. Mean sample size (n glac ) for the 55-km hex cells was 119, with a range of 1-585 ( Figure S3). We substituted h cell and A cell into Equation 1 to compute cell-level total mass balance M cell .
To estimate hex cell elevation change uncertainty (σ h cell ), we computed the sum of all individual glacier volume uncertainty (σ V = σ h ·A) estimates for the n glac glacier polygon centroids that fell within each hex cell, assuming they were correlated (based on the semivariogram analysis above): With this approach, we effectively weight the σ h estimate for each glacier by the relative glacier area A A cell . While the area terms in Equation (6) also have some uncertainty, the common 10% relative area uncertainty will cancel, so it does not directly contribute to the σ h cell estimate. We used Equation (4) to compute total mass balance uncertainty for the cell (σ M cell ), which accounts for this area uncertainty.
For subsequent aggregation over larger HMA sub-regions, we computed mean glacier elevation change ( h agg ) as the sum of all hex cell glacier volume change ( V cell = h cell ·A cell ) estimates for the n cell hex cell centroids that fell within each aggregation polygon containing total glacier area A agg : which is similar to Equation (5). We assumed that σ h cell errors for adjacent hex cells were independent and uncorrelated, and computed aggregated elevation change uncertainty (σ h agg ) in quadrature for all n cell hex cells that fell within each aggregation polygon: Sample sizes for this aggregation (n cell ) were 7-79 for HIMAP regions, 22-327 for Kääb et al. (2015) regions, and 3-117 for hydrologic basins, with a total of 793 hex cells for the full-HMA region. We used these values and Equations (1) and (4) to compute total mass balance M agg and uncertainty σ M agg for the aggregated regions.
To summarize, we estimated uncertainty for each glacier using observed dh dt values on surrounding terrain. We aggregated glacier mass balance and uncertainty values for spatially correlated hex cell samples, then aggregated these independent samples in quadrature for regional and full-HMA uncertainty estimates.

Glacier Meltwater Runoff
While we cannot directly estimate total glacier meltwater runoff contributions from glacierized HMA river basins using our geodetic mass balance estimates, we can estimate the "excess meltwater runoff " (e.g., Radić and Hock, 2014;Brun et al., 2017) or "imbalance component of runoff " (e.g., Pritchard, 2019). This excess glacier meltwater runoff is equal to the glacier mass loss in a given basin (Brun et al., 2017;Pritchard, 2019;Rounce et al., 2020b). In this study, we compute mass loss at multiple scales (individual glaciers, hex cells, basins), and we report excess glacier meltwater runoff for a given basin as the waterequivalent loss from only the glaciers or hex cells with negative mass balance ( M < 0) in that basin. The total excess glacier meltwater runoff can therefore be larger than the total glacier mass loss for the basin, as any glaciers or hex cells with positive mass balance are not included. To avoid thresholding issues (e.g., Anderson, 2019), we conservatively estimate excess glacier meltwater runoff uncertainty as the aggregated total glacier mass balance uncertainty for each basin.
To assess the relative importance of this excess glacier meltwater runoff, we compared with existing model output for the "total glacier meltwater runoff " and "total runoff " in each basin, which provided estimates for excess meltwater "percent of total glacier meltwater runoff " and "percent of total runoff." The total glacier meltwater runoff was estimated from Python Glacier Evolution Model (PyGEM, https://github.com/drounce/ PyGEM), forced with ERA-Interim reanalysis data (Rounce et al., 2020a,b). The monthly meltwater runoff was modeled for a single "gauging station" that moves with the terminus of each glacier and is the sum of modeled balance runoff expected if all glaciers were in equilibrium and any excess (or imbalance) runoff due to mass loss. The monthly values were annualized, and the 2000 to 2018 mean annual total glacier meltwater runoff was computed for each HMA basin, with final values of ∼0.1 to 25 Gt yr −1 (Rounce et al., 2020a,b).
The total runoff for each basin includes all water that moved from the land surface to the river system, including total glacier meltwater runoff, surface runoff (precipitation, snow melt), and groundwater baseflow. We used output from the Water Balance Model (WBM) (Wisser et al., 2010) forced with 1980-2016 ERA-Interim reanalysis data (Grogan, 2016) to estimate mean annual total runoff (at basin mouth for exorheic basins, and total area runoff for endorheic basins) during the 2000-2016 period (D. Grogan, personal communication). These basinspecific values ranged from ∼12 to 1,094 Gt yr −1 , and they show good agreement with available gauge observations (GRDC, 2018) and model output (Huss and Hock, 2018).

Glacier Mass Balance
The total HMA glacier mass balance for the period from 2000 to 2018 was −19.0 ± 2.5 Gt yr −1 (−0.19 ± 0.03 m w.e. yr −1 ). Maps of aggregated values for hex grid cells show the detailed spatial distribution of glacier mass change across the full HMA region (Figure 4, Figure S4). The Himalayan range displays relatively large mass loss, with some spatial variability. Clusters of grid cells with more positive mass balance are observed in the Kunlun and Karakoram ranges. The mass balance gradient across the Karakoram (slightly positive/nearbalance in the west, slightly negative in the east) is consistent with results from recent analyses (Lin et al., 2017;Berthier and Brun, 2019). Hex cell uncertainty values are highest in the Nyainqentanglha and Eastern Tien Shan, where total DEM counts are lowest ( Figure 1E) and individual glacier uncertainty is high.
Regional aggregations show greatest specific mass loss rates ( Figure 5) and total mass loss rates ( Figure S5) across Western, Central, and Eastern Himalayas, Nyainqentanglha, and the Central and Eastern Tien Shan. The Western Kunlun Shan and Eastern Pamir showed slightly positive mass balance, while the Karakoram showed slightly negative mass balance. Our results suggest that the HiMAP sub-region aggregation offers better spatial resolution of mass change patterns than the Kääb et al. (2015) sub-regions. This is especially true for the Tien Shan and Inner Tibetan Plateau, with more negative specific mass balance is observed in the eastern HiMAP sub-regions, though total mass loss for these sub-regions is relatively small (Figure S5).

Mass Balance and Glacier Area
Our results show that larger glaciers dominate total HMA mass balance (Figure 6). In terms of specific mass balance (m w.e. yr −1 ), we found that glaciers between ∼0.1 and ∼10 km 2 clustered around the mean specific mass balance for the full dataset ( Figure 6A). Glaciers smaller than ∼0.1 km 2 had more positive specific mass balance relative to the mean, while glaciers between ∼10 and ∼30 km 2 were more negative ( Figure 6A). Glaciers larger than ∼100 km 2 are concentrated in western HMA (Figure 6C), where more positive mass balance was observed during the period from 2000 to 2018. This indicates that the limited sample of very large glaciers is not necessarily representative of glacier change across the full HMA domain.
Total HMA-wide excess glacier meltwater was −22.71 ± 3.01 Gt yr −1 when calculated using all individual glaciers with negative mass balance, −20.61 ± 2.73 Gt yr −1 when calculated using all hex cells with negative mass balance, and −19.0 ± 2.5 Gt yr −1 when calculated using all aggregated basins with negative mass balance (identical to total mass balance, as all basins are negative). These differences emphasize the importance of aggregration level and methodology for this calculation. The Brahmaputra, Indus, and Ganges river basins had the greatest total excess glacier meltwater runoff, with glaciers contributing 5.23 (4.87) Gt yr −1 , 4.55 (3.92) Gt yr −1 , and 3.26 (3.19) Gt yr −1 above balance glacier meltwater runoff to their respective river systems during the 2000 to 2018 period ( Figure 7A) for all glaciers (hex cells).
This HMA contribution was ∼4-6% of the 0.92± 0.39 mm SLE yr −1 total contribution from all global glaciers and ice caps for this period (Zemp et al., 2019), and a small but non-negligible fraction (∼1.0-1.5%) of total global sea level rise estimates of 3.5

Glacier Mass Balance Methodology
The nearly complete inventory (99% of RGI polygons) of glacier mass balance provides an important snapshot of the detailed spatial pattern of glacier change across the region. We considered all glaciers in HMA (without a minimum area threshold), aggregated these glaciers into independent samples based on observed autocorrelation length scales, and then aggregated these independent samples over larger regions for improved uncertainty estimates. This approach offered a FIGURE 7 | Excess glacier meltwater runoff between 2000 and 2018 for major river basins in HMA: (A) Spatial distribution (Gt yr −1 ), and (B) as a percentage of mean annual PyGEM model total (combined balance and imbalance) glacier meltwater runoff (Rounce et al., 2020a,b). To highlight excess glacier meltwater runoff contributions for highly glacierized basins, bar width is scaled by the percent glacierized area in each basin (total glacier area divided by basin area, with values of ∼0.0-2.1%). Error bars combine basin-specific excess glacier meltwater uncertainty and interannual variability of modeled total glacier meltwater runoff for the 2000-2018 period. The Irrawady, Yellow, and Mekong basins have negligible glacierized area (∼30-200 km 2 ) and are sensitive to basin definitions and aggregation strategy, resulting in less-reliable percentages.
systematic regional assessment, with better coverage over subregions containing relatively small and/or sparse glaciers (e.g., eastern Tibetan Plateau) that could not be measured or were excluded from previous assessments.
A common issue with geodetic glacier mass balance studies is missing data over portions of glacier polygons, especially accumulation areas. These data gaps can be filled using several methods , which can lead to significant differences in estimates for glacier mass balance and uncertainty. Our study involves the most complete coverage for HMA glaciers to date, with mean polygon coverage of 99.5% and <5% missing data for 98% of all RGI glacier polygons, precluding the need for more complex gap-filling approaches.

Glacier Mass Balance Uncertainty
Our results underscore the importance of aggregation methodology for mass balance uncertainty calculations. Individual glacier mass balance uncertainty estimates can be large, often exceeding the magnitude of the apparent change, especially for smaller glaciers. Aggregation significantly reduces the combined mass balance uncertainty, though final estimates can vary greatly depending on methodology.
Unlike many other studies (Rolstad et al., 2009;Nuth and Kääb, 2011;Brun et al., 2017), we did not fit higher-order, 2-D polynomials to dh dt residuals and remove from static surfaces after co-registration. While doing so would significantly reduce apparent σ dh dt systematic bias and large-scale spatial autocorrelation, it would inevitably overfit sparse dh dt samples in many areas and adversely impact glacier mass balance estimates.
We found that large-scale autocorrelation length scales varied depending on the choice of error metric (section 2.5.2), with ∼24 km for σ h from NMAD, ∼32 km for σ h from RMSE, ∼27 km for σ M from NMAD, and ∼42 km for σ M from RMSE ( Figure S2). The increased range of influence for the RMSE metrics is likely related to spatially coherent residual dh dt bias over static surfaces. In addition, we found that the apparent autocorrelation length scales (range of influence) was sensitive to the step and total lag distance used during semivariogram analysis. Finally, we found that common assumptions about circular, isotropic range of influence are likely an oversimplification, as we observe different autocorrelation length scales for directional semivariogram analysis (dh_dt_sigma_error.ipynb in Bhushan, 2020).
Our hex cell approach attempts to account for spatial autocorrelation of σ h error between individual glacier polygon centroids. This is potentially an oversimplification. For one, adjacent polygons may share the same "stable" ground (e.g., ridge between two valley glaciers), resulting in greater apparent spatial autocorrelation than would be found with a regional analysis of σ dh dt random . In addition, rather than the geometric centroid of 2-D glacier polygons, it would be better to use the 3-D "center of mass" for each glacier from gridded ice thickness products (e.g., Huss and Farinotti, 2012), especially when aggregating over GRACE mascons. While this approach will not affect smaller glaciers, it could better represent longer glaciers with curved centerlines.
The hex cell and sub-region aggregation mass balance uncertainty is sensitive to the hex cell size, though full-HMA mass balance values are not affected. We experimented with several hex cell widths (27-55 km) and multiple offsets for the hex cell grid origin. Due to assumptions about spatial dependence (Equation 6), the individual hex cell uncertainty values increase with greater cell area, as do the corresponding uncertainty values for regional aggregation. A related issue involves width-dependent hex cell centroid locations and partitioning amongst different regional aggregation polygons (e.g., adjacent river basins). In practice, we did not observe significant differences in aggregated mass balance values for the range of hex cell widths tested.
Instead of an intermediate hex cell aggregation based on spatial autocorrelation of σ h values assigned to glacier centroids, one could explore a nested spherical semivariogram model (Rolstad et al., 2009) to simultaneously account for small-scale (∼500 m) and large-scale (∼20-30 km) autocorrelation of dh dt values over arbitrary aggregation areas. This nested approach would require a relatively large sample of well-distributed dh dt values over stable surfaces to properly capture the large-scale variability. There are also potential issues with this approach for complex (i.e., non-circular) aggregation polygon geometry (e.g., elongated mountain ranges) and observed anisotropy in spatial autocorrelation.
Our ∼10% glacier polygon area uncertainty assumption is likely overestimating actual uncertainty for the RGI v6.0 inventory, especially when aggregated over the entire region. While inventory error is difficult to assess, the study by Pfeffer et al. (2014) estimated area uncertainty of ∼8% for HMA and ∼5% for all glacier regions in the RGI v3.2 inventory. While an individual glacier may have area uncertainty of 10% or more (e.g., a relatively small, surging glacier), it is unlikely that all glaciers in aggregated regions have area uncertainty of 10%. For larger aggregation areas, the area uncertainty term dominates the total mass balance uncertainty in Equation (4), which justifies future efforts to improve time-variable inventories of glacier outlines.
We are hopeful that future community-driven, opensource software development will continue to standardize these methods, enabling more systematic glacier mass balance and uncertainty estimates for HMA and other glacierized regions.

Comparison With Previous Work
Our results are generally consistent with other full-HMA glacier mass balance assessments, including the comprehensive review by Bolch et al. (2019) and recent analysis by Brun et al. (2017). Considering the range of methods, datasets, time periods and spatial coverage of previous studies, we limit our direct comparison with the "reference" analysis by Brun et al. (2017) (similar methods, similar time period), and refer the reader to Brun et al. (2017) for detailed comparisons with earlier work, and to Bolch et al. (2019) for more complete context.

Full HMA Glacier Mass Balance
Our total HMA-wide mass balance of −19.0 ± 2.5 Gt yr −1 (−0.19 ± 0.03 m w.e. yr −1 ) is larger than the −16.3 ± 3.5 Gt yr −1 (−0.18 ± 0.04 m w.e. yr −1 ) estimate by Brun et al. (2017), although the two are within uncertainty bounds. While this difference in total mass loss rate may be related to differences in methodology, it could also be influenced by the inclusion of an additional 2+ years (2016)(2017)(2018) of DEMs later in the 2000 to 2018 study period. Based on observed long-term glacier mass balance trends for HMA (e.g., Mukherjee et al., 2018;Zhou et al., 2018;Maurer et al., 2019) and glaciers worldwide, we might expect greater mass loss toward the end of the record as opposed to earlier in the record. Despite a more negative total mass balance, we find a smaller sea level rise contribution from exhoreic basins (13.3 ± 2.3 Gt yr −1 ) compared to the 14.6 ± 3.1 Gt yr −1 estimate by Brun et al. (2017).

Sub-region Glacier Mass Balance
We used the recently published HiMAP sub-region definitions  during aggregation, which provide improved geographic partitioning of mass balance for HMA mountain ranges, especially over the Inner Tibetan Plateau and Tien Shan. However, most existing full-HMA geodetic mass balance assessments are limited to the sub-regions defined by Kääb et al. (2015). We generally observe good agreement between our results and those of Brun et al. (2017) for the Kääb et al. (2015) sub-region definitions ( Figure 5D, Figure S5D). The most notable differences for total mass balance in Gt yr −1 (Figure S5D) are over the Kunlun, where our mass gain estimates are smaller, and the Inner Tibetan Plateau where our mass loss estimates are larger. Larger differences in specific mass balance ( Figure 5D) are observed for the Eastern Himalaya and Nyainqentanglha. These differences potentially arise from relatively limited observations for these sub-regions in the study by Brun et al. (2017).

River Basins and Meltwater Runoff
Our estimates for total mass balance in HMA river basins are generally similar to those by Brun et al. (2017) (Figure S6), with notable differences for the Tarim and Inner Tibetan Plateau ( Figure S6D). In both cases, our mass loss estimates are larger than those of Brun et al. (2017), with an apparent sign change for Tarim (−0.87± 0.71 Gt yr −1 in our study vs. +0.4± 1.3 Gt yr −1 in Brun et al., 2017), which is important for full HMA excess glacier meltwater estimates. Some of this discrepancy may be related to basin boundary definitions, but our estimates are consistent with observed differences in glacierized sub-region mass balance estimates over the Kunlun and Tien Shan (Figure S5).
While the total mass balance estimates are important for assessing HMA contributions to sea level rise, the excess glacier meltwater runoff potentially provides information on the longterm mean of seasonal glacier contributions to downstream hydrology. The largest differences between the two metrics are observed for basins with a relatively high percentage of glaciers with positive mass balance (e.g., Indus, Tarim). Our estimates for total excess glacier meltwater runoff are generally similar to those provided by Brun et al. (2017). We note that estimates for excess glacier meltwater runoff vary based on aggregation methodology, with more negative values obtained when aggregating at the glacier polygon level compared with the hex cell or basin level (used by Brun et al., 2017).
To provide context for our observed excess glacier meltwater runoff estimates, we considered the percentage of modeled total glacier meltwater runoff and total runoff in each basin. While these estimates vary based on glacier/hydrologic model, forcing data, and true interannual variability, the long-term mean runoff (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) should capture the relative magnitudes across HMA basins. Our results show that a significant portion of large, endorheic basins on the Tibetan Plateau contain relatively high percentages of excess glacier meltwater (Figure 7). In other words, the long-term mean of seasonal discharge in these river systems may be disproportionately sourced from thinning glaciers that are not being replenished by "excess accumulation" during the accumulation season.
Pritchard (2019) used the excess (imbalance) glacier meltwater estimates by Brun et al. (2017) for a subset of glacierized HMA basins, and estimated that total HMA glacier meltwater production was 1.6x the balance meltwater rates, with higher ratios in drought years. This value would correspond to approximately 38% if plotted on Figure 7B (excess glacier meltwater runoff percentage of total glacier meltwater runoff for each basin), which is slightly higher than our mean value, but comparable. Pritchard (2019) also estimated that the total glacier meltwater runoff contributions were 0.3-3.2% of "gross water inputs" in an average climatic year, though these percentages increased in drought years. The magnitude is similar to our estimates of 0.0-3.3% for the percentage of excess glacier meltwater runoff compared to the 2000-2016 mean of basinspecific total runoff values. While a detailed evaluation is beyond the scope of this study, these percentages offer a sense of the relatively small, but non-negligible contribution of excess glacier meltwater runoff in each basin for a typical year.

Glacier Mass Balance Methodology
Beyond the inherent value of reproducible science, it is worth noting the differences in data and methodology used in our study compared to those used in the study by Brun et al. (2017). We used the more recent, multi-source RGI v6.0 glacier inventory vs. the single-source GAMDAM inventory (Nuimura et al., 2015), processed ASTER DEMs with improved stereo correlation methodology to reduce data gaps (especially over accumulation areas), incorporated an additional 2+ years of ASTER DEMs, integrated an additional 5+ years of high-resolution DEMs from WorldView/GeoEye imagery with improved resolution, dynamic range and geolocation accuracy, used the accurate, self-consistent TanDEM-X global DEM as a reference for co-registration, implemented robust approaches to estimate elevation change trends for each glacier, and relied on a different method to estimate glacier mass balance and uncertainty (without a need to bin dh dt by elevation bands for regional analyses). Considering the range of methodological differences, it is very promising to see good agreement between our results and those of Brun et al. (2017) across the full HMA region and most sub-regions. As a community, we have transitioned from a large spread in estimates for HMA glacier mass balance (factor of ∼2-3), to less than ∼15-20% uncertainty in recent years. The agreement demonstrates value of publicly available remote sensing archives that span multiple decades, and robust approaches involving elevation change trend estimation using long records of potentially noisy observations.

Factors Controlling Glacier Mass Balance
Our assessment of glacier mass balance, which includes smaller glaciers across the Tibetan Plateau, can potentially provide an important constraint for climate reanalysis data evaluation (e.g., Immerzeel et al., 2015). The observed spatial variability in glacier mass balance depends on climatology (forcing) and glacier sensitivity to this forcing. The former involves variables such as temperature, precipitation, seasonality, and insolation. The latter involves factors such as local environmental setting, geometry, ice thickness, bed properties, debris cover properties, and the presence of proglacial lakes. Recent work suggests that glacier sensitivity to climate is the dominant control on HMA mass balance, at least during the 2003-2009 period (Sakai and Fujita, 2017). In general, maritime HMA glaciers are most sensitive to air temperature, while continental HMA glaciers are sensitive to a combination of temperature, precipitation seasonality, and snow/rain differentiation . Our results show an apparent mass gain in the Kunlun Shan and East Pamir, with relatively limited mass loss in the Karakoram and West Pamir (Figure 5), which is consistent with both the observed increases in winter snow water equivalent from 1987 to 2009 (Smith and Bookhagen, 2018), and the higher sensitivity of these continental glaciers to precipitation. Our results also show local spatial variability in glacier mass balance across the region (Figure 4). This is consistent with the hypothesis that glacier sensitivity, controlled for example by variable glacier geometry and terminus type, may play an important role in local glacier mass balance patterns (e.g., Brun et al., 2019). While a detailed analysis of these factors is beyond the scope of this paper, our observations provide a new baseline to interpret glacier mass loss in the context of regional climatology and local geomorphology.

Mass Balance and Glacier Area
The contribution of small glaciers (<0.1 km 2 ) to total HMA mass balance is relatively limited (see annotations in Figure 6C). For example, 46 large glaciers with area >100 km 2 lost more mass than >55,000 glaciers with area <0.3 km 2 ( Figure 6C). However, our results suggest that small glaciers displayed a specific mass balance that was approximately half the specific mass balance of larger glaciers. This is consistent with geodetic mass balance results for the Swiss Alps (Fischer et al., 2015). We now consider several possible explanations and explore further.
There are potential issues with our ability to assess small glacier change. We used the multi-source RGI glacier inventory, which was digitized from medium-resolution imagery (∼30 m) with variable digitization area thresholds used by different operators. As a result, there is likely a bias in the spatial distribution of small glaciers in this inventory. Figure 6C shows that the spatial distribution of small glaciers with area between 0.03 and 0.1 km 2 in the RGI inventory offers a good regional sample. However, many small glaciers with area between 0.01 and 0.03 km 2 are located in the Karakoram and Kunlun. This observation potentially indicates a real increase in the density of small glaciers and/or a bias in the digitization of small glaciers in these sub-regions, which had more positive mass balance than the regional mean for the 2000-2018 period.
It is also possible that the 30 m resolution of our dh dt products is too coarse to properly resolve small glaciers, or that our DEM and trend filtering methodology preferentially smoothed values over small glaciers. For reference, a 0.1 km 2 glacier polygon spans ∼300 × 300 m or ∼10 × 10 pixels for a 30-m elevation change analysis, while a 0.01 km 2 glacier polygon spans ∼100 × 100 m or ∼3 × 3 pixels. In theory, this 9-pixel sample should be sufficient to estimate mass balance, but in practice (and after filtering during DEM generation and post-processing), this is approaching a minimum threshold for significance. To test this, we repeated our analysis without filtering the trend maps. The uncertainty for individual glaciers increased, and the spread of binned mass balance for small glaciers increased, but we observed a bin mean/median distribution that was nearly identical to Figure 6A. Another potential complication involves the static RGI polygon outlines that we used for the 2000-2018 period. Smaller glaciers should experience a greater percent area change, with increased influence of reduced apparent dh dt over surfaces that became ice free between 2000 and 2018.
Many of these small glaciers may actually be perennial snowfields and/or rock glaciers, which we might expect to display reduced elevation change compared to clean ice. This is challenging to assess without high-resolution observations of glacier velocity and debris cover. If we assume that these features are clean ice, the reduced loss rates from small glaciers could be related to the fact that they do not descend to low elevations and hence experience cooler climate conditions. Small glaciers also tend to occupy protected alcoves with limited insolation and more windblown/avalanche snow accumulation. These factors favor the preservation of small glaciers even under conditions of climate warming (DeBeer and Sharp, 2009). Smaller, thinner glaciers should also have shorter response times and will reach equilibrium more quickly than larger glaciers under the same climate forcing (Bahr et al., 1998).
Regardless of the explanation, this result highlights the importance of minimum glacier area threshold selection for geodetic mass balance studies. Small glaciers have received increased attention in recent years (e.g., Huss and Fischer, 2016) and their ∼0.17-0.53 mm SLE yr −1 contribution in the past century is needed to close the total sea level rise budget (Parkes and Marzeion, 2018). It is possible to introduce systematic bias if calculating specific mass balance only from larger glaciers and then multiplying by total glacier area in a region. While the total mass balance and recent/future contribution to sea level rise may not be significantly affected, there are potential implications for water resource management. For example, small glaciers may have increased significance for dry season streamflow in some watersheds, particularly interior endorheic basins where glacier meltwater runoff represents a higher percentage of total runoff. Similarly, future projections from models that were calibrated using only larger glaciers in the region could incorrectly estimate timescales for loss of smaller glaciers.

CONCLUSIONS AND SUMMARY
We generated ∼34000 high-resolution DEMs and produced the first regional high-resolution DEM composites for HMA using stereo imagery acquired by multiple satellite platforms (WorldView-1/2/3, GeoEye-1, and ASTER) between 2000 and 2018. We produced robust elevation change trend maps and geodetic mass balance estimates for 99% of HMA glaciers over the full 18-year period, and aggregated independent samples based on observed spatial autocorrelation of uncertainty estimates. The HMA-wide glacier mass balance was −19.0 ± 2.5 Gt yr −1 (−0.19 ± 0.03 m w.e. yr −1 ) during this period, with potential sea-level rise contribution of up to 0.052 ± 0.007 mm SLE yr −1 (0.037 ± 0.006 mm SLE yr −1 from exorheic basins). We documented the spatial pattern of HMA glacier mass change with greater detail and coverage than previous assessments, with aggregated estimates of glacier mass change for glacierized sub-regions and hydrologic basins. We observed negative mass balance across nearly all sub-regions, with greatest total mass loss across the Himalayas, Nyainqentanglha, and the Tien Shan and positive mass balance in the Western Kunlun Shan and Eastern Pamir. We found that smaller glaciers (<0.1 km 2 ) had less negative specific mass balance than the regional mean, with implications for regional scaling of mass balance estimates from a limited glacier sample. Finally, observed excess glacier meltwater runoff due to negative mass balance was ∼12-53% of total basin-specific glacier meltwater runoff and ∼0.0-3.3% of total basin-specific runoff between 2000 and 2018.
The mass balance results presented here provide a new baseline for future studies of glacier mass balance and landscape evolution in HMA. These results are currently being used for calibration and validation of glacier mass balance models, GRACE trends in total water storage, and hydrologic models used to assess present and future water resources in the region.
Ongoing tasking campaigns for in-track WV/GE sub-meter stereo imagery will fill gaps in existing DEM coverage and provide better repeat high-resolution DEM coverage for HMA glaciers. Combining these growing archives with additional ASTER stereo DEMs, improved ASTER artifact correction algorithms (e.g., Girod et al., 2017), and new satellite laser altimetry archives (NASA ICESat-2 and Global Ecosystem Dynamics Investigation [GEDI]) should enable new elevation change analyses with robust trend fits over shorter intervals. The refined observations of spatiotemporal evolution will enable detailed analysis of the relationships between observed glacier mass balance and evolving climate forcing, surface processes affecting glacier sensitivity (e.g., debris cover evolution), and glacier dynamics-all of which must be better understood to properly assess the present and future role of HMA glaciers for downstream water resources.

DATA AVAILABILITY STATEMENT
The DEM data used in this study are available through the National Snow and Ice Data Center (NSIDC): https://nsidc. org/data/highmountainasia. All tools used for processing and analysis are available on Github, with scripts used for specific processing steps mentioned throughout the text. These datasets and tools will continue to improve and evolve in the future, and we welcome community contributions. A dedicated Github repository contains the glacier mass balance data and copies of processing and analysis scripts/notebooks (https://zenodo.org/ record/3600624).

AUTHOR CONTRIBUTIONS
DS led the effort, processed WV/GE DEMs, developed tools, and performed analysis. SB performed spatial autocorrelation analysis and assisted with uncertainty analysis. PM processed the ASTER DEMs. DR provided PyGEM glacier runoff estimates and assisted with interpretation of regional results. AA coordinated HiMAT project integration and data sharing. BO assisted with the trend-fitting analysis. All authors contributed to discussions and manuscript preparation.

FUNDING
This effort was supported by NASA's High-mountain Asia research program (NNH15ZDA001N-HMA). DS, SB, and AA were supported by NASA award NNX16AQ88G. PM, DR, and BO were supported by NASA award NNX17AB27G. AA was partially supported by the Washington Research Foundation, and by a Data Science Environments project award to the University of Washington eScience Institute from the Gordon and Betty Moore and the Alfred P. Sloan Foundations.