Himalayan glaciers supply meltwater to densely populated catchments in South Asia, and regional observations of glacier change over multiple decades are needed to understand climate drivers and assess resulting impacts on glacier-fed rivers. Here, we quantify changes in ice thickness during the intervals 1975–2000 and 2000–2016 across the Himalayas, using a set of digital elevation models derived from cold war–era spy satellite film and modern stereo satellite imagery. We observe consistent ice loss along the entire 2000-km transect for both intervals and find a doubling of the average loss rate during 2000–2016 [−0.43 ± 0.14 m w.e. year−1 (meters of water equivalent per year)] compared to 1975–2000 (−0.22 ± 0.13 m w.e. year−1). The similar magnitude and acceleration of ice loss across the Himalayas suggests a regionally coherent climate forcing, consistent with atmospheric warming and associated energy fluxes as the dominant drivers of glacier change.
INTRODUCTION
The Intergovernmental Panel on Climate Change 5th Assessment Report estimates that mass loss from glaciers contributed more to sea-level rise than the ice sheets during 1993–2010 (0.86 mm year−1 versus 0.60 mm year−1, respectively), yet uncertainties for the glacier contribution are three times greater (1). Glaciers also contribute locally to water resources in many regions and serve as hydrological buffers vital for ecology, agriculture, and hydropower, particularly in High Mountain Asia (HMA), which includes all mountain ranges surrounding the Tibetan Plateau (2, 3). Shrinking Himalayan glaciers pose challenges to societies and policy-makers regarding issues such as changing glacier melt contributions to seasonal runoff, especially in climatically drier western regions (3), and increasing risk of outburst floods due to expansion of unstable proglacial lakes (4). Yet, substantial gaps in knowledge persist regarding rates of ice loss, hydrological responses, and associated climate drivers in HMA (2).
Mountain glaciers are known to respond dynamically to a variety of drivers on different time scales, with faster response times than the large ice sheets (5, 6). In the Himalayas, in situ studies document significant interannual variability of mass balances (7–9) and relatively slower melt rates on debris-covered glacier tongues over interannual time scales (10, 11). Yet, the overall effects of surface debris cover are uncertain, as many satellite observations suggest similar ice losses relative to clean-ice glaciers over similar or longer periods (12, 13). Because of the complex monsoon climate in the Himalayas, dominant causes of recent glacier changes remain controversial, although atmospheric warming, the albedo effect due to deposition of anthropogenic black carbon (BC) on snow and ice, and precipitation changes have been suggested as important drivers (14–16).
Model projections of future Himalayan ice loss and resulting impacts require robust observations of glacier response to past and ongoing climate change. Recent satellite remote sensing studies have made substantial advances with improved spatial coverage and resolution to quantify ice mass changes during 2000–2016 (12, 17, 18), and former records extending back to the 1970s have been presented for several areas using declassified spy satellite imagery (13, 19–22). These long-term records are especially critical for extracting robust mass balance signals from the noise of interannual variability (6). Many studies also report the highly heterogeneous behavior of glaciers in localized regions, with some glaciers exhibiting faster rates of ice loss during the 21st century (20, 22). Independent analyses document simultaneously increasing atmospheric temperatures at high-elevation stations in HMA (23–26). To robustly quantify the regional sensitivity of these glaciers to climate change, a reliable Himalaya-wide record of ice loss extending back several decades is needed.
Here, we provide an internally consistent dataset of glacier mass change across the Himalayan range over approximately the past 40 years. We use recent advances in digital elevation model (DEM) extraction methods from declassified KH-9 Hexagon film (27) and ASTER stereo imagery to quantify ice loss trends for 650 of the largest glaciers during 1975–2000 and 2000–2016. All aspects of the analysis presented here only use glaciers with data available during both time intervals unless specified otherwise. We investigate glaciers along a 2000-km transect from Spiti Lahaul to Bhutan (75°E to 93°E), which includes glaciers that accumulate snow primarily during winter (western Himalayas) and during the summer monsoon (eastern Himalayas), but excludes complications of surging glaciers in the Karakoram and Kunlun regions where many glaciers appear to be anomalously stable or advancing (2). Our compilation includes glaciers comprising approximately 34% of the total glacierized area in the region, which represents roughly 55% of the total ice volume based on recent ice thickness estimates (15, 28). This diverse dataset adequately captures the statistical distribution of large (>3 km2) glaciers, thus providing the first spatially robust analysis of glacier change spanning four decades in the Himalayas. We extract DEMs from declassified KH-9 Hexagon images for the 650 glaciers, compile a corresponding set of modern ASTER DEMs, fit a robust linear regression through every 30-m pixel of the time series of elevations, sum the resulting elevation changes for each glacier, divide by the corresponding areas, and translate the volume changes to mass using a density conversion factor of 850 ± 60 kg m−3(see Materials and Methods).
RESULTS
Glacier mass changes
Our results indicate that glaciers across the Himalayas experienced significant ice loss over the past 40 years, with the average rate of ice loss twice as rapid in the 21st century compared to the end of the 20th century (Fig. 1). We calculate a regional average geodetic mass balance of −0.43 ± 0.14 m w.e. year−1 (meters of water equivalent per year) during 2000–2016, compared to −0.22 ± 0.13 m w.e. year−1 during 1975–2000 (−0.31 ± 0.13 m w.e. year−1 for the full 1975–2016 interval) (see Materials and Methods). A 30-glacier moving average shows a quasi-consistent trend across the 2000-km longitudinal transect during both time intervals (Fig. 1), and subregions have similar means and distributions of glacier mass balance. Some central catchments deviate somewhat from the Himalaya-wide mean during 2000–2016 (by approximately 0.1 to 0.2 m w.e. year−1) in the Uttarakhand (~79.0° to 80.0°E), the Gandaki catchment (~83.5° to 84.5°E), and the Karnali catchment (~81° to 83°E), which has fewer larger glaciers and relatively incomplete data coverage. Similar to previous in situ and satellite-based studies (18, 29), we observe considerable variation among individual glacier mass balances, with area-weighted SDs of 0.1 and 0.2 m w.e. year−1during each respective interval for the 650 glaciers. This variability most likely reflects different glacier characteristics such as sizes of accumulation zones relative to ablation zones, topographic shading, and amounts of debris cover. Yet, we find that, in our survey (using a rough average of 60 glaciers per 7000-km2 subregion), local variations tend to average out and mean values are similar across most catchments.
Contrasting distributions of glacier mass balances are evident when comparing between time intervals. The 1975–2000 distribution has a negative tail extending to −0.6 m w.e. year−1, while the 2000–2016 distribution is more negative, extending to −1.1 m w.e. year−1 (Fig. 2A). During the more recent interval, glaciers are losing ice twice as fast on average (Fig. 2B), though this varies somewhat between subregions. For example, we find that the average rate of ice loss has increased by a factor of 3 in the Spiti Lahaul region, and by a factor of 1.4 in West Nepal. We also compile altitudinal distributions of ice thickness change for the glaciers and create a Himalaya-wide average thickness change profile versus elevation (Fig. 2, C and D). These distributed thinning profiles are a function of both thinning by mass loss and of dynamic thinning due to ice flow. We find that the 2000–2016 thinning rate (m year−1) profile is considerably steeper, which is likely caused by a combination of faster mass loss and widespread slowing of ice velocities during the 21st century (2, 30).
We multiply geodetic mass balances by the full glacierized area in the Himalayas between 75° and 93° longitude to estimate region-wide ice mass changes of −7.5 ± 2.3 Gt year−1 during 2000–2016, compared to −3.9 ± 2.2 Gt year−1 during 1975–2000 (−5.2 ± 2.2 Gt year−1 during the full 1975–2016 interval). Recent models using Shuttle Radar Topography Mission (SRTM) elevation data for ice thickness inversion estimate the total glacial ice mass in our region of study to be approximately 700 Gt in the year 2000 (see Materials and Methods) (15, 28). If this estimate is accurate, our observed annual mass losses suggest that of the total ice mass present in 1975, about 87% remained in 2000 and 72% remained in 2016.
Comparison of clean-ice, debris-covered, and lake-terminating glaciers
We study mass changes for different glacier types by separating glaciers into clean-ice (<33% area covered by debris), debris-covered (≥33% area covered by debris), and lake-terminating categories based on a Landsat band ratio threshold and manual delineation of proglacial lakes (see Materials and Methods). All three categories have undergone a similar acceleration of ice loss (Table 1), and debris-covered glaciers exhibit similar and often more negative geodetic mass balances compared to clean-ice glaciers over the past 40 years (Fig. 3). Altitudinal distributions indicate slower thinning for lower-elevation regions of debris-covered glaciers (glacier tongues where debris is most concentrated) relative to clean-ice glaciers, but comparatively faster thinning in mid- to upper elevations (Fig. 4). Lake-terminating glaciers concentrated in the eastern Himalayas exhibit the most negative mass balances due to thermal undercutting and calving (31), though they only comprise around 5 to 6% of the estimated total Himalaya-wide mass loss during both intervals.Table 1Himalaya-wide geodetic mass balances (m w.e. year−1).View this table:
Approximation of required temperature change
As a first approximation of the consistency between observed glacier mass balances and available temperature records, we estimate the energy required to melt the observed ice losses and conservatively estimate the atmospheric temperature change that would supply this energy via longwave radiation to the glaciers, using a simple energy balance approach (Materials and Methods). We propagate significant uncertainties associated with input from global climate reanalysis data, scaling of temperatures from coarse reanalysis grids to specific glacier elevations, and averaging of climate data over the glacierized region. Results suggest that the observed acceleration of ice loss can be explained by an average temperature ranging from 0.4° to 1.4°C warmer during 2000–2016, relative to the 1975–2000 average. This approximately agrees with the magnitude of warming observed by meteorological stations located throughout HMA, which have recorded air temperatures around 1°C warmer on average during 2000–2016, relative to 1975–2000 (Fig. 5). More comprehensive climate observations and models will be essential for further investigation, but these simple energy constraints suggest that the acceleration of mass loss in the Himalayas is consistent with warming temperatures recorded by meteorological stations in the region.
DISCUSSION
Implications for dominant drivers of glacier change in the Himalayas
The parsing of Himalayan glacier energy budgets is not a straightforward task owing to the scarcity of meteorological data, in combination with the complex climate and topography of the region (2). Furthermore, the Himalayas border hot spots of high anthropogenic BC emissions, which may affect glaciers by direct heating of the atmosphere and decreasing albedo of ice and snow after deposition (14). While improved analyses combining observations and high-resolution atmospheric and glacier energy balance models will be required to quantify these effects, the pattern of ice loss we observe has important implications regarding dominant climate influences on Himalayan glacier mass balances. Our results suggest that any drivers of glacier change must explain the region-wide consistency, the doubling of the average rate of ice loss in the 21st century compared to 1975–2000, and the observation that clean-ice, debris-covered, and lake-terminating glaciers have all experienced a similar acceleration of mass loss.
Some studies have suggested a weakening of the summer monsoon and reduced precipitation as primary reasons for negative glacier mass balances, particularly in the Everest region (16). While decreasing accumulation rates may account for a significant portion of the mass balance signal for some glaciers, an extreme Himalaya-wide decrease in precipitation would be required to explain the extensive ice losses we observe, especially given that monsoon-dominated glaciers with high accumulation rates are known to be much more sensitive to temperature than accumulation changes (5, 32). Regional studies of precipitation trends in the Himalayas do not suggest a widespread decrease in precipitation over the past four decades (Supplementary Materials). A uniform BC albedo forcing across the Himalayas is another possible explanation, although BC concentrations measured in snow and ice in the Himalayas have been found to be spatially heterogeneous (14, 33), and high-resolution atmospheric models also show large spatial variability of deposited BC originating from localized emissions in regions of complex terrain (14, 34). Future analyses focused on quantifying the spatial patterns of BC deposition will reveal further insights, yet given the rather homogeneous pattern of mass loss we observe across the 2000-km Himalayan transect, a strong, spatially heterogeneous mechanism seems improbable as a dominant driver of glacier ice loss in the region.
Debris-covered glaciers
Similar thinning rates of debris-covered (thermally insulated) glaciers relative to clean-ice glaciers have been observed by previous studies and have been not only ascribed to surface melt ponds and associated ice cliffs acting as localized hot spots to concentrate melting but also attributed to declining ice flux causing dynamic thinning and stagnation of debris-covered glacier tongues (2). While we cannot yet directly deconvolve relative contributions from these processes, we find that average thinning rates for debris-covered glaciers are slower than clean-ice glaciers at low elevations (glacier tongues where debris is most concentrated), which agrees with reduced melt rates from field studies. In turn, debris-covered glaciers tend to have comparatively faster thinning at mid-range elevations, where debris cover is sparser and also where the majority of total glacierized area resides (Fig. 4). Models of debris-covered glacier processes suggest that this pattern of thinning may cause a reduction in down-glacier surface gradient, which, in turn, reduces driving stress and ice flux and explains why debris-covered ablation zones become stagnant (35). We also find that clean-ice glaciers exhibit a much more pronounced steepening of the thinning profile over time, compared to debris-covered glaciers. It may be that both glacier types experience a uniform thinning phase followed by a changing terminus flux and retreat phase, but the clean-ice glaciers are in a later phase of response to recent climate change (36).
Comparison with previous studies in the Himalayas
To compare our results with previous remote sensing studies that derive mass changes from various sensors (including Hexagon, SRTM, SPOT5, ICESat, and ASTER), we restrict our results to the approximate geographical regions covered by each corresponding study (12, 13, 17–22) and then compute area-weighted average geodetic mass balances. In addition, we compare individual glacier mass balances for the Everest and Langtang Himal regions, where mass changes were previously estimated using declassified Corona and Hexagon imagery (13, 19, 20). Despite factors such as variable spatial resolutions, distinct void-filling methods, heterogeneous spatial and temporal coverages, and different definitions of glacier boundaries, we find that our average mass balances generally agree with previous analyses and overlap within uncertainties (table S1). However, because of the significant variability of individual glacier mass changes within subregions, our results also highlight the importance of sampling a large number of glaciers to obtain a robust average trend for any given area.
Comparison with benchmark mid-latitude glaciers and global average
To gain perspective on mass losses from these low-latitude glaciers in the monsoonal Himalayas, we compare our results with benchmark mid-latitude glaciers in the European Alps, as well as with a global average mass balance trend (fig. S1) (37). The Alps contain the most detailed long-term glaciological and high-elevation meteorological records on Earth, and the climatic sensitivity and behavior of these European glaciers are well understood compared to glaciers in HMA. Air temperatures in the Alps show an abrupt warming trend beginning in the mid-1980s, and Alpine mass balance records display a concurrent acceleration of ice loss, with a continual strongly negative mass balance after that time. Himalayan weather station data indicate a more gradual warming trend, with the strongest warming beginning in the mid-1990s (fig. S1, A and B). We find that mass balance in the Himalayas is less negative compared to the Alps and the global average, despite close proximity to a known hot spot of increasing BC emissions with rapid growth and accompanying combustion of fossil fuels and biomass in South Asia (38). The concurrent acceleration of ice loss observed in both the Himalayas and Europe over the past 40 years coincides with a distinct warming trend beginning in the latter part of the 20th century, followed by the consistently warmest temperatures through the 21st century in both regions.
Conclusion
Our analysis robustly quantifies four decades of ice loss for 650 of the largest glaciers across a 2000-km transect in the Himalayas. We find similar mass loss rates across subregions and a doubling of the average rate of loss during 2000–2016 relative to the 1975–2000 interval. This is consistent with the available multidecade weather station records scattered throughout HMA, which indicate quasi-steady mean annual air temperatures through the 1960s to the 1980s with a prominent warming trend beginning in the mid-1990s and continuing into the 21st century (23–26). We suggest that degree-day and energy balance models focused on accurately quantifying glacier responses to air temperature changes (including energy fluxes and associated feedbacks) will provide the most robust estimates of glacier response to future climate scenarios in the Himalayas.
MATERIALS AND METHODS
Hexagon
U.S. intelligence agencies used KH-9 Hexagon military satellites for reconnaissance from 1973 to 1980. A telescopic camera system acquired thousands of photographs worldwide, after which film recovery capsules were ejected from the satellites and parachuted back to Earth over the Pacific Ocean. With a ground resolution ranging from 6 to 9 m, single scenes from the mapping camera cover an area of approximately 30,000 km2 with overlap of 55 to 70%, allowing for stereo photogrammetric processing of large regions. Images were scanned by the U.S. Geological Survey (USGS) at a resolution of 7 μm and downloaded via the Earth Explorer user interface (earthexplorer.usgs.gov). Digital elevation models were extracted using the Hexagon Imagery Automated Pipeline methodology, which is coded in MATLAB and uses the OpenCV library for Oriented FAST and Rotated BRIEF (ORB) feature matching, uncalibrated stereo rectification, and semiglobal block matching algorithms (27). The majority of the KH-9 images here were acquired within a 3-year interval (1973–1976), and we processed a total of 42 images to provide sufficient spatial coverage (fig. S2).
ASTER
The ASTER instrument was launched as part of a cooperative effort between NASA and Japan’s Ministry of Economy, Trade and Industry in 1999. Its nadir and backward-viewing telescopes provide stereoscopic capability at 15-m ground resolution, and a single DEM covers approximately 3600 km2. Approximately 26,000 ASTER DEMs were downloaded via the METI AIST Data Archive System (MADAS) satellite data retrieval system (gbank.gsj.jp/madas), a portal maintained by the Japanese National Institute of Advanced Industrial Science and Technology and the Geological Survey of Japan. To use all cloud-free pixels (including images with a high percentage of cloud cover), no cloud selection criteria were applied when downloading the images. We used the Data1.l3a.demzs geotiff product, which has a spatial resolution of 30 m. After downloading, the DEMs were subjected to a cleanup process: For a given scene, any saturated pixels (i.e., equal to 0 or 255) in the nadir band 3 (0.76 to 0.86 μm) image were masked in the DEM. Next, the SRTM dataset was used to remove any DEM values with an absolute elevation difference larger than 150 m (relative to SRTM), which effectively eliminated the majority of errors caused by clouds. While more sophisticated cloud masking procedures are possible, the ASTER shortwave infrared detectors failed in April 2008, making cloud detection after this time impossible using standard methods. We examined existing cloud masks derived using Moderate Resolution Imaging Spectroradiometer images as another option (tonolab.cis.ibaraki.ac.jp/ASTER/cloud/). However, these are not optimized for snow-covered regions and often misclassify glacier pixels as clouds. Instead, our large collection of multitemporal ASTER scenes, the SRTM difference threshold, and our robust linear trend fitting algorithm [see description of Random Sample Consensus (RANSAC) in the “Trend fitting of multitemporal DEM stacks” section] effectively excluded any remaining erroneous cloud elevations after the initial threshold. As a final step, all ASTER DEMs were coregistered to the SRTM using a standard elevation–aspect optimization procedure (39). We did not apply fifth-order polynomial correction procedures to the ASTER DEMs for satellite “jitter” effects and curvature bias as done in some previous studies (18). We found that while these types of corrections may reduce the overall average elevation error in a scene, the polynomial fitting can be unwieldy and may introduce unwanted localized biases. By stacking many ASTER DEMs (with 20.5 being the average number of observations per pixel stack during the ASTER trend fitting, see fig. S3E) and using a robust fitting procedure, we found that biases do not correlate across overlapping scenes, and thus tend to cancel out one another. Furthermore, the elevation change results from this portion of our study overlap within uncertainties with Brun et al. (18) (Supplementary Materials) who did perform polynomial corrections. This suggests that for a large-scale regional study using a high number of overlapping ASTER scenes, the satellite jitter and curvature bias corrections have a relatively minimal impact on the final results.
Glacier polygons
To delineate glaciers during all portions of the analysis, we used manually refined versions of the Randolph Glacier Inventory (RGI 5.0) (40). Starting with the original RGI dataset, we edited the glacier polygons to reflect glacier areas during 1975, 2000, and 2016. For the 1975 edit, we used a combination of Hexagon imagery, the Global Land Survey (GLS) Landsat Multispectral Scanner mosaic (GLS1975), and glacier thickness change maps derived from differencing the Hexagon and modern ASTER DEMs, which are particularly useful for debris-covered glacier termini that often have spectral characteristics indistinguishable from surrounding terrain. Debris-covered areas for each glacier were delineated using a Landsat DN TM4/TM5 band ratio with a threshold of 2.0, and glaciers with ≥33% debris cover were assigned to the debris-covered category. For the 2000 edit, we used the GLS2000 Landsat Enhanced Thematic Mapper Plus mosaic, along with glacier thickness change maps derived from differencing ASTER DEMs. For the 2016 edit, we used a custom mosaic of Landsat 8 imagery with acquisition dates spanning 2014–2016. The individually edited glacier polygons were used for all ice volume change and geodetic mass balance computations. The polygons were also used during alignment of the DEMs, where the shapefiles were converted to raster masks with a dilation (morphological operation) of 250 m on the binary rasters. This effectively excluded unstable terrain surrounding the glaciers during the DEM alignment process, such as steep avalanching slopes and unstable moraines.
Trend fitting of multitemporal DEM stacks
Glacier polygons were processed individually—all DEMs from a given time interval (1975–2000 or 2000–2016) that overlap a polygon were selected and resampled to the same 30-m resolution using linear interpolation. The overlapping DEMs were sampled with a 25% extension around each glacier to include nearby stable terrain for alignment and uncertainty analysis (fig. S4). After ensuring that there is no vertical bias, the aligned DEMs were sorted in temporal order as a three-dimensional matrix, and linear trends were fit to every pixel “stack” (i.e., along the third dimension of the matrix) using the RANSAC method. During each RANSAC iteration, a random set of two elevation pixels per stack were selected. A linear trend was fit to these two values, and then all remaining elevation pixels were compared to the trend. Any elevation pixels within 15 m of the trend line were marked as inliers. This process was repeated for 100 iterations, and the iteration with the greatest number of inliers was selected. A final linear fit was performed using all inliers from the best iteration, and this trend was used for each pixel stack’s thickness change estimate. The thickness change maps were subjected to outlier removal using thresholds for maximum slope, maximum thickness change, minimum count per pixel stack, minimum timespan per pixel stack, maximum SD of inlier elevations per pixel stack, and maximum gradient of the thickness change map (fig. S3). In addition, the thickness change pixels were separated into 50-m elevation bins, and pixels falling outside the 2 to 98% quantile range were excluded. Any bins with less than 100 pixels were removed and then interpolated using the two adjacent bins. Before computing ice volume change for the glaciers, the final thickness change maps were visually inspected, any remaining erroneous pixels (which occurred almost exclusively in low-contrast, snow-covered accumulation zones) were manually masked, and a 5 × 5 pixel median filter was applied. We did not attempt to perform seasonality corrections, as no seasonal snowfall records are available and because nearly all the Hexagon DEMs were acquired during winter, thus minimizing any seasonality offsets between regions. For the 1975–2000 interval, we used the Hexagon DEMs and sampled ASTER thickness change trends at the start of the year 2000. For regions with multiple overlapping Hexagon DEMs, we used the same RANSAC method. During the 1975–2000 interval, only two DEMs were available for most glaciers. In these cases, the RANSAC iterations were unnecessary, and we simply differenced the two available DEMs. We did not use SRTM for any thickness change estimates; thus, no correction for radar penetration was necessary.
Mass changes
To compute (mean annual) ice volume changes for individual glaciers, all thickness change pixels falling within a glacier polygon were transformed to an appropriate projected WGS84 UTM coordinate system (zones 43 to 46, depending on longitude of the glacier). Pixel values (m year−1) were then multiplied by their corresponding areas (pixel width × pixel height) and summed together. The resulting ice volume change was then divided by the average glacier area to obtain a glacier thickness change. We used the average of the initial and final glacier areas for a given time interval and excluded slopes greater than 45° to remove any cliffs and nunataks. Last, the glacier thickness change was multiplied by an average ice-firn density (41) of 850 kg m−3 and then divided by the density of water (1000 kg m−3) to compute glacier geodetic mass balance in m w.e. year−1. Because of cloud cover, shadows, and low radiometric contrast, some glacier accumulation zones had gaps in the DEMs and resulting thickness change maps. This is particularly evident in the Hexagon DEMs for the Spiti Lahaul region owing to extensive cloud cover. To fill these gaps, we tested two different void-filling methods for comparison. In the first method, we defined a circular search area with a radius of 50 km around the center of a given glacier. All thickness change pixels from glaciers in this surrounding area were binned (into 50-m elevation bins, and following the same outlier-removal procedure given in the preceding section), and any missing data in the glacier were set to this “regional bin” mean value at the corresponding elevation. In the second method, we filled data gaps using an interpolation procedure, where voids in an individual glacier were linearly interpolated using bin values at upper and lower elevations relative to the missing data (those belonging to the same glacier), and assumed zero change at the highest elevation bin (headwall). Both methods yielded similar results (table S1). In addition, no obvious trends were apparent when geodetic mass balance was plotted versus percent data coverage or glacier size (fig. S5). While smaller glaciers exhibited more scatter, the average mass balance was similar for all glacier sizes. These observations indicate that our representative sample of glaciers, while biased toward larger glaciers, adequately captures the statistical distribution of glacier mass balances in the Himalayas.
To calculate regional geodetic mass balances, we separated glaciers into four subregions (Spiti Lahaul, West Nepal, East Nepal, and Bhutan) as defined by Brun et al. (18). We then calculated the average mass balance for each of these four subregions, weighted by individual glacier areas. Last, we calculated a final average mass balance for the Himalayas, weighted by the total glacierized area (from the RGI 5.0 database) in each of the four subregions, between 75° to 93° longitude. Because of the relatively homogeneous mass balance distribution, we found that this approach resulted in similar values (well within the uncertainties) compared to simply calculating the area-weighted average mass balance of the 650 measured glaciers. To obtain the total mass changes in Gt year−1, we multiplied each subregion mass balance by its total glacierized area and then summed the results from all subregions to get Himalaya-wide totals of −3.9 Gt year−1 for 1975–2000 and −7.5 Gt year−1 for 2000–2016. To calculate contributions to sea-level rise, we used a global ocean surface area of 361.9 × 106 km2 (fig. S4G).
To estimate the total ice mass present in our region of study, we used ice thickness estimates from Kraaijenbrink et al. (15), who used the Glacier bed Topography version 2 model to invert for ice thickness (28) with input from the SRTM DEM (acquired in February of 2000). The ice thickness estimates from (15) did not include glaciers smaller than 0.4 km2, and to estimate the additional mass contribution from these smallest glaciers (along with any other glaciers that are missing thickness estimates), we fit a second-order polynomial to the logarithm of glacier volumes versus the logarithm of glacier areas and evaluated this fit equation for any glaciers without volume data (fig. S6). We then converted glacier volume to mass using a density value of 850 kg m−3. Over our region of study, the ice volumes from the thickness data amounted to 649 Gt, with an additional contribution of 35 Gt from the fitting procedure, for a total of 684 Gt.
Uncertainty assessment
We quantified statistical uncertainty for individual glaciers using an iterative random sampling approach. For a given glacier, the SD of elevation changes from the surrounding stable terrain (σz) was first calculated. For any missing thickness change pixels within the glacier polygon, we also included an extrapolation uncertainty σe. This accounts for additional error in regions with incomplete data, i.e., those glacier regions filled by extrapolating thickness changes from surrounding glaciers, or linear interpolation assuming zero change at the headwall, as described in the previous section. We found that in the Himalaya-wide altitudinal distributions, the maximum SD of thickness change in any 50-m elevation bin above 5000 m is 0.56 m year−1. Nearly all regions with incomplete data coverage are above this elevation, resulting from poor radiometric contrast for snow-covered glacier accumulation zones. We thus conservatively set σe equal to 0.6 m year−1. We then combined both sources of error to get σp for every individual thickness change pixelσp=σ2z+σ2e−−−−−−√(1)
To account for spatial autocorrelation, we started with a normally distributed random error field (with a mean of 0 and an SD of 1) the same size as the thickness change map and then filtered it using an n-by-n moving window average to add spatial correlation, where n is defined as the spatial correlation range divided by the spatial resolution of the thickness change map. We used 500 m for the spatial correlation range, which is a conservative value based on semivariogram analysis in mountainous regions from previous studies (18, 21, 42). The resulting artificial error field En (now with spatial correlation) is scaled by the σp values and added to the thickness change map ΔH as follows, where (x, y) are pixel coordinatesΔHE(x,y)=ΔH(x,y)+En(x,y)⋅σp(x,y)σn(2)
If thickness change data exist at a given pixel location (x, y) on the glacier, σn is the SD of the set of all En values where data exist (i.e., where σe is equal to zero). Conversely, if thickness change data do not exist at a given pixel location (x, y) on the glacier, σn is the SD of the set of all En values where data do not exist (i.e., where σe is equal to 0.6 m year−1). In this way, the second term of Eq. 2 assigns larger uncertainties to regions with incomplete data. Last, all glacier thickness change pixels in ΔHE were summed together to compute a volume change with the introduced error. This procedure was repeated for 100 iterations, and the volume change uncertainty σΔV was computed as the SD of the resulting distribution (fig. S4). For region-wide volume change estimates, we conservatively assumed total correlation between glaciers and computed region-wide uncertainty as follows, where g is the total number of glaciers (~17,400)σΔV region=∑1gσΔV(3)
For glaciers where thickness change data are not available, a measure of uncertainty is still required to factor into the final regional uncertainty estimate. For these glaciers, we estimated σΔVas (42)σΔV=σ2z region⋅Acor5⋅A−−−−−−−−−−−√(4)Acor=π⋅L2(5)
In this case, σz region is the region-wide SD of elevation change over stable terrain (0.42 m year−1) (fig. S7), Acor is the correlation area, L is the correlation range (500 m), and A is the glacier area. Last, all σΔV and σΔV regjon estimates were combined with an area uncertainty (43) of 10% and a density uncertainty (41) of 60 kg m−3 using standard uncorrelated error propagation.
Sensitivity of region-wide glacier mass change estimates
We further tested the sensitivity of our region-wide estimates to potential biases, including (i) the exclusion of small glaciers, (ii) incomplete data coverage for many glacier accumulation zones during 1975–2000, and (iii) void-filling technique. First, we note that our geodetic mass balance analysis only includes glaciers larger than 3 km2. This is because mass balance uncertainties increase with decreasing glacier size, and we find that uncertainties often exceed the magnitude of mass changes for glaciers smaller than ~3 km2. To test whether the neglected small glaciers appreciably affect the result, we also computed mass balances using all available glaciers (i.e., all glaciers with ≥33% data coverage, including those smaller than 3 km2). We find that including the full set of smaller glaciers changes the region-wide geodetic mass balance estimates by a maximum of 0.04 m w.e. year−1 (fig. S4G). Next, we note that the Hexagon DEMs in particular have poor data coverage over glacier accumulation zones (figs. S8 and S9). However, the vast majority of thinning occurs in glacier ablation zones, and the amount of thinning decreases with elevation in a quasi-linear fashion, especially in mid- to upper regions of the glaciers where data gaps are most common. Thus, we hypothesize that we can extrapolate and interpolate with reasonable confidence over accumulation areas. To test the robustness of this assumption, we used the 2000–2016 glacier change data. The ASTER data over this interval have superior radiometric contrast and adequately capture elevation change trends for most accumulation zones. We first set all 2000–2016 thickness change pixels to be empty where the 1975–2000 data are missing to simulate the same data gaps over accumulation zones as in the 1975–2000 data. We then performed the same geodetic mass balance calculations and found that the region-wide geodetic mass balance only changes by 0.01 m w.e. year−1 (fig. S4G, comparing test 3 to test 1). Last, we performed two separate void-filling methods for all tests (see the “Mass changes” section for descriptions of void-filling methods) and observed a maximum change in geodetic mass balance of 0.04 m w.e. year−1. Overall, the relatively small impact of each test suggests that our results are robust to the exclusion of small glaciers, incomplete data coverage over glacier accumulation zones, and void-filling technique.
SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/6/eaav7266/DC1
Fig. S1. Comparison of Himalayan temperature trends and regional mass balance with benchmark mid-latitude glaciers and a global average trend.
Fig. S2. Coverage of glacierized area in the Himalayas.
Fig. S3. Trend fit examples for two large glaciers using ASTER DEMs during 2000–2016, histograms of ASTER pixel counts and timespans per stack (glacier averages), and outlier thresholds.
Fig. S4. Illustration of uncertainty estimation procedure for a single iteration/glacier and Himalaya-wide sensitivity tests.
Fig. S5. Geodetic mass balances during 1975–2000 and 2000–2016 plotted against various parameters.
Fig. S6. Log-log plot of glacier volumes versus areas used to estimate the total ice mass present in our region of study.
Fig. S7. Analysis of elevation change for nonglacier pixels (stable terrain) during both intervals.
Fig. S8. Thickness change maps used in the analysis.
Fig. S9. Thickness change maps for the three remaining Himalayan regions.
Table S1. Geodetic mass balance comparisons with prior studies.
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.
