Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Accelerated global glacier mass loss in the early twenty-first century

Abstract

Glaciers distinct from the Greenland and Antarctic ice sheets are shrinking rapidly, altering regional hydrology1, raising global sea level2 and elevating natural hazards3. Yet, owing to the scarcity of constrained mass loss observations, glacier evolution during the satellite era is known only partially, as a geographic and temporal patchwork4,5. Here we reveal the accelerated, albeit contrasting, patterns of glacier mass loss during the early twenty-first century. Using largely untapped satellite archives, we chart surface elevation changes at a high spatiotemporal resolution over all of Earth’s glaciers. We extensively validate our estimates against independent, high-precision measurements and present a globally complete and consistent estimate of glacier mass change. We show that during 2000–2019, glaciers lost a mass of 267 ± 16 gigatonnes per year, equivalent to 21 ± 3 per cent of the observed sea-level rise6. We identify a mass loss acceleration of 48 ± 16 gigatonnes per year per decade, explaining 6 to 19 per cent of the observed acceleration of sea-level rise. Particularly, thinning rates of glaciers outside ice sheet peripheries doubled over the past two decades. Glaciers currently lose more mass, and at similar or larger acceleration rates, than the Greenland or Antarctic ice sheets taken separately7,8,9. By uncovering the patterns of mass change in many regions, we find contrasting glacier fluctuations that agree with the decadal variability in precipitation and temperature. These include a North Atlantic anomaly of decelerated mass loss, a strongly accelerated loss from northwestern American glaciers, and the apparent end of the Karakoram anomaly of mass gain10. We anticipate our highly resolved estimates to advance the understanding of drivers that govern the distribution of glacier change, and to extend our capabilities of predicting these changes at all scales. Predictions robustly benchmarked against observations are critically needed to design adaptive policies for the local- and regional-scale management of water resources and cryospheric risks, as well as for the global-scale mitigation of sea-level rise.

Main

About 200 million people live on land that is predicted to fall below the high-tide lines of rising sea levels by the end of the century11, whereas more than one billion could face water shortage and food insecurity within the next three decades4. Glaciers distinct from the ice sheets hold a prominent role in these outcomes as the largest estimated contributor to twenty-first century sea-level rise after thermal expansion2, and as one of the most climate-sensitive constituents of the world’s natural water towers12,13. Current glacier retreat temporarily mitigates water stress on populations reliant on ice reserves by increasing river runoff1, but this short-lived effect will eventually decline14. Understanding present-day and future glacier mass change is thus crucial to avoid water-scarcity-induced sociopolitical instability15, to predict the alteration of coastal areas due to sea-level rise4, and to assess the impacts on ecosystems16 and cryosphere-related hazards3.

Nevertheless, glacier mass change stands out as one of the least-constrained elements of the global water cycle, identified as a critical research gap in the Special Report on the Ocean and Cryosphere in a Changing Climate (SROCC) of the Intergovernmental Panel on Climate Change (IPCC)4. Observational limits stem from the fragmented expanse of glacierized surfaces around the globe. Largely inaccessible, only a few hundred of the more than 200,000 glaciers are monitored in situ17. Notwithstanding recent progress in glacier monitoring from space18, global-scale remote-sensing-based studies have been so far limited by (i) the coarse spatial resolution of satellite gravimetry, which is unable to reliably disentangle glacier mass change signals from those of the ice sheets, solid Earth and hydrology in many regions5,19,20; (ii) the sparse repeat sampling of satellite altimetry that operated over short timespans5,10; and (iii) the uneven coverage of optical and radar surface elevation change estimations that account at most for 10% of the world’s glaciers21.

Spatiotemporally resolved estimation

In this study, we leverage large-scale and openly available satellite and airborne elevation datasets as a means of estimation, reference or validation of Earth’s surface elevation over all glaciers and their vicinity between 1 January 2000 and 31 December 2019 (Extended Data Fig. 1). For observational coverage, we rely mostly on NASA’s 20-year-long archive of stereo images from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). We use modern photogrammetry techniques and specifically developed statistical methods to generate and bias-correct nearly half a million Digital Elevation Models (DEMs) at 30 m horizontal resolution. In total, our repeat DEMs cover more than 20 times Earth’s land area (Extended Data Fig. 2).

Changes in glacier elevation based on DEMs are traditionally quantified by differencing pairs of acquisitions from two distinct epochs. To harness the full potential of the repeat temporal coverage provided by the archives, we introduce an approach to producing continuous elevation time series interpolated from all available DEMs (see Methods, Extended Data Fig. 3). This technique allows us to mitigate the strong effects of seasonality while preserving longer, nonlinear glacier elevation changes through time. In total, we independently compute surface elevation time series for about half a billion pixels at a horizontal resolution of 100 m, covering 97.4% of inventoried glacier areas22, with an average of 39 independent observations per pixel between 2000 and 2019 (Extended Data Table 2). Using glacier-wide hypsometric gap-filling methods, we then extend our estimated elevation changes to nearly 99.9% of glacier areas.

We perform an extensive validation by intersecting our elevation time series with 25 million high-precision measurements from NASA’s Ice, Cloud, and land Elevation Satellite (ICESat) and Operation IceBridge campaigns over glaciers, spanning 2003 to 2019. We thereby confirm the absence of temporal and spatial biases in our elevation change estimates (see Methods, Extended Data Fig. 4). We further utilize ICESat data to constrain the spatiotemporal correlations that are either structural to our interpolated elevation time series or that emerge owing to latent, uncorrected ASTER instrument noise, and we propagate our elevation uncertainties into volume change uncertainties accordingly. We validate the reliability of our uncertainty estimates down to the scale of individual glaciers by comparison with independent, high-precision DEM differences for 588 glaciers around the globe (Extended Data Fig. 5).

Integration of elevation changes over each of the 217,175 inventoried glaciers yields volume change, which is subsequently converted to water-equivalent mass change23. Our analysis includes 200,000 km2 of glaciers located near the coast of Greenland (Greenland Periphery) and in the Antarctic seas (Antarctic and Subantarctic), referred to as peripheral glaciers, that are distinct from the Greenland Ice Sheet (GIS) and the Antarctic Ice Sheet (AIS). We aggregate our estimates over the 19 first-order regions of the Randolph Glacier Inventory 6.0 (RGI 6.0)22 (Fig. 1), and report estimates for periods exceeding five years owing to larger uncertainties at shorter timescales (Extended Data Table 1). Uncertainties, provided at the 95% confidence level (two standard deviations), depend primarily on observational coverage. When converting from volume to mass change, our estimates are largely hampered by a poor knowledge of density conversion23, which constitutes the dominant uncertainty component of our glacier mass change assessment.

Fig. 1: Regional glacier mass changes and their temporal evolution from 2000 to 2019.
figure 1

Regional and global mass change rates with time series of mean surface elevation change rates for glaciers (indigo) of the 19 first-order RGI 6.022 regions (white-delimited indigo polygons; region numbers indicated in parentheses), shown on top of a world hillshade36. Regions 2, 5, 9, 17 are further divided (N, S, E and W indicating north, south, east and west, respectively) to illustrate contrasting temporal patterns. Mass change rates are represented by the area of the disk delimiting the inside wedge, which separates the mass change contribution of land-terminating (light grey) and marine-terminating (light blue) glaciers. Mass change rates larger than 4 Gt yr−1 are printed in blue inside the disk (in units of Gt yr−1). The outside ring discerns between land (grey) and marine-terminating (blue) glacier area. Annual time series of mean elevation change (in m yr−1) and regional data coverage are displayed on time friezes at the bottom of the disks.

Source data.

Global contribution to sea-level rise

From 2000 to 2019, global glacier mass loss totalled 267 ± 16 Gt yr−1 (Extended Data Table 1), a mass loss 47% larger than that of the GIS, and more than twice that of the AIS7,8,9 (Table 1). Assuming that all meltwater ultimately reached the ocean, the contribution to sea-level rise was 0.74 ± 0.04 mm yr−1 or 21 ± 3% of the observed rise24. Global glacier mass loss rapidly accelerated (see Methods) at a rate of 48 ± 16 Gt yr−1 per decade (62 ± 8 Gt yr−1 per decade excluding peripheral glaciers), corresponding to a thinning rate acceleration of 0.10 ± 0.02 m yr−1 per decade (0.16 ± 0.02 m yr−1 per decade). While thinning rates increased steadily, mass loss acceleration slightly attenuated in time owing to the decreasing extent of glacier surfaces caused by glacier retreat. Excluding peripheral glaciers, thinning rates nearly doubled, from 0.36 ± 0.21 m yr−1 in 2000 to 0.69 ± 0.15 m yr−1 in 2019. Observational studies have been unable to discern significant (95% confidence interval does not overlap zero) accelerated glacier mass loss19,21, with the exception of a recent gravimetric study20 that estimated an acceleration of 50 ± 40 Gt yr−1 per decade excluding peripheral glaciers. Despite its large uncertainties, this estimate is in agreement with our results. The observed acceleration of mass loss for glaciers exceeds that of the GIS7 and is similar to that of the AIS8. For the AIS, gravimetric observations indicate a decelerating mass loss since the mid-2010s25. We thereby infer that acceleration of sea-level rise since 2000, which is often attributed to the accelerated loss from both the GIS and AIS, also substantially originates from glaciers. Observed sea-level trends24 place the glacier contribution at 6–19% of the acceleration in global sea-level rise, with a mean estimate at 9%. The large spread of this contribution primarily arises from uncertainties in the observed acceleration of sea-level rise24.

Table 1 Separating mass losses of glaciers and ice sheets

Marine-terminating glaciers collectively represent 40% of Earth’s total glacierized area, yet only contribute 26% to the global mass loss (Fig. 1). This smaller contribution to sea-level rise is uniform for all maritime regions, except where losses of marine-terminating glaciers are dominated by recent large surge events (for example, Svalbard and Jan Mayen; Extended Data Fig. 6). The delayed and asynchronous response of tidewater glaciers to changes in climate26 may partly explain why most marine-terminating glaciers show reduced mass loss. Despite differing mass loss rates, relative acceleration of land- and marine-terminating glaciers within each maritime region are similar (Extended Data Table 3). Notable exceptions exist for glaciers in the Antarctic and Subantarctic, where few land-terminating glaciers are present, and in regions of strong surge-driven mass losses.

Regionally contrasting mass changes

Seven glacierized regions account for 83% of the global mass loss (Extended Data Table 1): Alaska (25%); the Greenland Periphery (13%); Arctic Canada North and South (10% each); Antarctic and Subantarctic, High Mountain Asia (composed of Central Asia, South Asia West and South Asia East) and the Southern Andes (8% each). From 2000 to 2019, specific-mass change (that is, mass change divided by area) strongly varied in latitudinal belts (Fig. 2). The large, northernmost Arctic regions composed of Arctic Canada North, northern Greenland Periphery, Svalbard and Jan Mayen, and the Russian Arctic, all showed moderate specific-mass change rates, averaging –0.28 ± 0.04 metres water equivalent (w.e.) per year. Further South in the Arctic (at latitudes encompassing Alaska, Arctic Canada South, southern Greenland Periphery, Iceland and Scandinavia) specific-mass change rates were consistently more negative, at a near-triple value of −0.74 ± 0.10 m w.e. yr−1, reaching the world’s most negative regional rate over these two decades of −0.88 ± 0.13 m w.e. yr−1 in Iceland. Non-polar regions also experienced substantial mass loss (−0.69 ± 0.11 m w.e. yr−1 on average) with the exception of High Mountain Asia (−0.22 ± 0.05 m w.e. yr−1). The Antarctic and Subantarctic exhibited the least-negative specific-mass change rate of −0.17 ± 0.04 m w.e. yr−1.

Fig. 2: Spatial distribution of glacier elevation change between 2000 and 2019.
figure 2

Glacier mean elevation change rate is aggregated for tiles of 1° × 1° below 60° latitude, 2° × 1° between 60° and 74° latitude and 2° × 2° above 74° latitude, thus representing similar surface areas of approximately 10,000 km2. Disks scale with the glacierized area of each tile and are coloured according to the mean elevation change rate (coloured in grey if less than 50% of the surface is covered by observations or if the 95% confidence interval is larger than 1 m yr−1; only applies to 0.4% of the glacierized area) shown on top of a world hillshade36. Tiles with glacierized areas less than or equal to 10 km2 are displayed at the same minimum size.

Source data.

Our regional mass change estimates closely match those of a recent gravimetric study19 in remote polar regions (Arctic Canada, Svalbard and Jan Mayen, and the Russian Arctic) in which gravimetric uncertainties are considered small owing to weak competing signals27 (Fig. 3). We note, however, the large discrepancies between the latter gravimetric study19 and a more recent one20 in both Iceland and the Russian Arctic. We find good agreement with the dense in situ measurements of Central Europe, Scandinavia and New Zealand5. In High Mountain Asia and the Southern Andes, where gravimetric and glaciological records are less constrained, our mass change estimates of −0.21 ± 0.05 m w.e. yr−1 and −0.67 ± 0.15 m w.e. yr−1, respectively, are slightly more negative than the −0.19 ± 0.06 m w.e. yr−1 and −0.64 ± 0.04 m w.e. yr−1 reported by recent DEM-based studies28,29. For glaciers located in the tropics (Low Latitudes), our estimate of −0.43 ± 0.12 m w.e. yr−1 is about twice as negative as that of a recent interferometric radar study29, −0.23 ± 0.08 m w.e. yr−1, a difference that plausibly originates from biases associated with the penetration of radar signals into ice and firn30.

Fig. 3: Comparison to previous global estimates.
figure 3

Regional and global specific-mass change rates with 95% confidence intervals. a, Estimates from earlier studies5,19,20,21 (coloured) and from this study (black) over the corresponding time periods. b, Global estimates from earlier studies, and annual and 5-year rates from this study. Annual rates from earlier studies are not shown owing to large uncertainties. Uncertainties from this study are conservative, in particular for regions with most negative specific rates (see Methods).

Source data.

Drivers of temporal variabilities

While global glacier mass loss distinctly accelerated, the loss from glaciers peripheral to the GIS and AIS slightly decelerated, from 65 ± 16 Gt yr−1 in 2000–2004 to 43 ± 13 Gt yr−1 in 2015–2019 (Extended Data Table 1). Variability within the ice sheet peripheries was strong, however (Figs. 2, 4a). The peculiar surface elevation change patterns that we capture for glaciers fringing Greenland, particularly notable around the eastern Greenland sub-regions of mass gain in 2015–2019 (Extended Data Fig. 7), mirror those observed by satellite radar altimetry for the outer parts of the GIS31. Similarly, the elevation change rate patterns of Antarctica’s scattered peripheral glaciers largely agree with mass changes reported for the AIS8. Western Antarctic peripheral glaciers substantially thinned (−0.23 ± 0.06 m yr−1) while those of East Antarctica slowly thickened (0.04 ± 0.05 m yr−1). Ice masses surrounding the Antarctic Peninsula, representing 63% of the glacier area in the Antarctic and Subantarctic, experienced moderate, decelerating thinning (−0.19 ± 0.05 m yr−1), in line with recent gravimetric surveys of the entire peninsula25.

Fig. 4: Decadal patterns of glacier thinning are consistent with decadal variations in precipitation and temperature.
figure 4

ac, Difference between 2010–2019 and 2000–2009 for the mean elevation change rates from this study (a) and the mean annual precipitation (b), and mean annual temperature (c) from ERA5. The region inset numbering corresponds to Fig. 2. Tiles are 1° × 1°, with size inversely scaled to uncertainties in the mean elevation change difference shown on top of a world hillshade36. The minimum tile area is 10% for a 95% confidence interval larger than 1 m yr−1, and tiles are displayed at full size for a 95% confidence interval smaller than 0.4 m yr−1.

Source data.

Only two regions of the world beyond the ice sheet peripheries experienced slowdown of glacier thinning. The record thinning rates of Icelandic glaciers during 2000–2004 (1.21 ± 0.18 m yr−1) were nearly halved during 2015–2019 (0.77 ± 0.13 m yr−1), which coincides with the decelerated thinning of Scandinavian glaciers. Both are well corroborated by in situ observations21. Taken together, the slowdown in mass loss from these two regions, in addition to the one of peripheral glaciers of the southeast Greenland Periphery32, define a regional pattern that we refer to as the North Atlantic anomaly.

Elsewhere on Earth, glacier thinning accelerated. The combined mass loss of these regions with increased loss escalated from 148 ± 19 Gt yr−1 in 2000–2004 to 247 ± 20 Gt yr−1 in 2015–2019. Two-thirds of this increase derives from three regions: Alaska (38%), High Mountain Asia (19%) and Western Canada and USA (9%). Glaciers in the latter region experienced a fourfold increase in thinning rates. Most notably, glaciers in northwestern America (Alaska, Western Canada and USA) are responsible for nearly 50% of the accelerated mass loss. The widespread and strong increase of thinning of glaciers in High Mountain Asia brought a large sub-region of sustained thickening in central–western Asia down to a generalized thinning in the late 2010s (Extended Data Fig. 7), suggesting the end of the so-called Karakoram anomaly10. Smaller glacierized regions also underwent strong, sometimes drastic acceleration of thinning. New Zealand, for example, shows a record thinning rate of 1.52 ± 0.50 m yr−1 in 2015–2019, which is a nearly sevenfold increase compared to 2000–2004.

Analysis of climate data reveals that many of the regional patterns of mass change uncovered by our resolved estimates are consistent with large-scale, decadal changes in annual precipitation and temperature (Fig. 4b, c). Strong dipoles that reflect concordant spatial patterns between precipitation change and mass change are observed notably in northwestern America, the southern Greenland Periphery and the Southern Andes. The southern Andean dipole is consistent with the mega-drought33 of the 2010s that drove increased glacier mass loss in the Central Andes. In the Coast Mountains of western Canada and in southeast Alaska, glaciers were severely deprived of precipitation, which instead benefited neighbouring regions of central Alaska and continental USA, correspondingly showing either stable or reduced mass loss. The North Atlantic anomaly coincided with the cool, wet conditions of the last decade. Weaker dipoles can also be observed within the European Alps or Scandinavia. In both regions, glacier thinning slightly accelerated in the northeast and decelerated in the southwest.

Although decadal changes in precipitation explain some of the observed regional anomalies, the global acceleration of glacier mass loss mirrors the global warming of the atmosphere. Aggregated globally over glacierized terrain, we observe modest trends in precipitation during the period 2000–2019 (0.002 m yr−1, +6.2% in 20 years), whereas we detect a strong increase in air temperature (0.030 K yr−1). Combined with our estimate of accelerated mass change, this warming trend yields an observational global glacier mass balance sensitivity to temperature of −0.27 m w.e. yr−1 K−1, in agreement with modelling-based estimates34. Previous studies35 have indicated large multi-decadal variation in rates of glacier mass change across the 20th century, implying that some of the acceleration that we observe could fall within the range of natural variability. Nonetheless, the strong concordance to the increase in global surface temperatures suggests, indirectly, a considerable response to anthropogenic forcing. Together, the contrasting patterns and global-scale sensitivities consistent with meteorological conditions support the notion of a long-term, temperature-driven acceleration in glacier mass loss13 that is still subject to regional and sub-decadal precipitation-driven fluctuations of large magnitude.

Two decades of observational wealth

Benefiting from the nearly complete spatial coverage afforded by ASTER stereo imagery, our global estimate of recent glacier mass change (−275 ± 17 Gt yr−1 for the 2006–2015 IPCC SROCC reference period) shows strongly reduced uncertainties compared to the latest IPCC report4 (−278 ± 226 Gt yr−1) and a recent global study21 (−335 ± 144 Gt yr−1). We resolve the time-varying nature of this mass change signal for nearly all of Earth’s glaciers, which reveals a significantly accelerated mass loss globally. Decadal rates of glacier mass change remain, however, strongly modulated by regional climatic conditions. We capture the magnitude of such fluctuations, most contrasting for North Atlantic and northwestern American glaciers that evolved in opposing directions. At the end of the 2010s, the North Atlantic anomaly brought a whole sub-region of the eastern Greenland Periphery close to balance, whereas the strong increase in thinning rates of High Mountain Asian glaciers probably marks the end of the Karakoram anomaly.

From the spatiotemporally resolved nature of our assessment, multiple possibilities arise to harness observations of the satellite era. Such resolved estimates are not only instrumental for glaciers, but also hold the potential to constrain recent ice sheet mass balance, in particular for the outlet glaciers that are prone to the highest long-term sea-level rise. The improved ability to deconvolve glacier signals from gravimetric observations might foster the detection of nearly two decades of changes in terrestrial water storage. In time, we expect our observational baseline to help drive the development of the next generation of global glaciological and hydrological models, and to ultimately result in more reliable projections at all scales14. In light of the rapid, ongoing change of the cryosphere, the increasingly reliable projections made possible by accurate, global-scale observations are critical for the design of adaptation strategies, with impacts ranging from further sea-level rise4,11 to changes in water management for some of the most vulnerable regions on Earth12,15.

Methods

We summarize the workflow used to process elevation datasets into estimates of glacier mass change for the period of 1 January 2000 to 31 December 2019 (Extended Data Fig. 1).

Glacier inventories

We used the Randolph Glacier Inventory 6.0 (RGI 6.0)22 outlines for all regions except for Caucasus Middle East (region 12). Owing to the high number of uncharted (‘nominal’) glaciers in that region, we updated our inventory with the latest Global Land Ice Measurements from Space (GLIMS) outlines available37. This increased the number of glacier outlines in region 12 from 1,888 to 3,516, representing an increase in total area from 1,307 km2 to 1,336 km2. In Svalbard and Jan Mayen (region 7), we manually updated glacier outlines to account for advances resulting from major surges38,39,40, increasing mapped areas by 228 km2 (Extended Data Fig. 6). In the Greenland Periphery (region 5), we did not analyse the 955 glaciers strongly connected to the ice sheet (RGI 6.0 connectivity level 2) with an area of 40,354 km2, because these are generally included within the ice sheet by studies on the GIS7,9. Our updated inventory numbers 217,175 glaciers covering a total area of 705,997 km2. For the purpose of co-registering and bias-correcting DEMs, we masked ice-covered terrain using the RGI 6.0 for glaciers, the Greenland Ice Mapping Project41 for the GIS, and Bedmap242 for the AIS.

Digital elevation models

We retrieved all ASTER43, ArcticDEM44 and Reference Elevation Model of Antarctica (REMA)45 data intersecting glaciers worldwide (Extended Data Fig. 2), totalling more than 100 TB of data. Because of the non-negligible effects of radar penetration into snow and ice30, we excluded radar elevation datasets from our analysis except for the TanDEM-X 90 m global DEM46 (TanDEM-X). We used TanDEM-X as a globally homogeneous reference47 for co-registration48 and bias correction over ice-free terrain, keeping only elevations with an error smaller than 0.5 m in the provided TanDEM-X height error map. For all DEMs bilinearly resampled to 30 m, co-registration was performed only if more than 100 valid elevation differences (slope >3°, absolute elevation difference <200 m) were available at each iterative step.

From 440,548 ASTER L1A stereo images43 (each covering 60 km × 60 km), we generated, co-registered and bias-corrected 154,656 ASTER DEM strips (30 m resolution; 180 km × 60 km strip size) using improved techniques of MicMac for ASTER49,50. Improvements were made by adjusting the back-looking image for cross-track biases before stereo calculations, by accounting for the curved along-track angle of the satellite Terra, and by stitching the arbitrarily split 60 km × 60 km archive granules into longer strips. The latter operation mitigates edge effects and increases the amount of ice-free terrain available for improved basin-hopping optimizations51 of along-track undulations and satellite jitter parameters. Further details on the processing of ASTER DEMs are available in Supplementary Information.

From 97,649 release 7 ArcticDEM44 DEMs at 2 m resolution and 13,790 release 1.1 REMA45 DEMs at 2 m and 8 m resolution, we stitched and co-registered 40,391 ArcticDEM and 3,456 REMA longer strips to TanDEM-X. Our stitching of the original DEM segments, generated by the Polar Geospatial Center using the Surface Extraction with TIN-based search-space minimization algorithm52, was performed by a sequential pairwise co-registration between same-day acquisitions over all available terrain. This procedure was necessary to increase the amount of ice-free terrain in the final DEM strip for co-registration to TanDEM-X. We allowed for a maximum standard deviation of co-registered elevation differences of 10 m before stopping the sequential co-registration iteration and starting a new strip, instead of the 1-m root-mean-square error originally used44,45.

Elevation time series

Following co-registration, we excluded all DEMs for which the root-mean-square error of the elevation difference with TanDEM-X on ice-free terrain was larger than 20 m. Using all remaining DEMs, we created three-dimensional arrays (time t, space x and y) of elevation h(txy), divided into 2,106 tiles of 1° × 1° containing glaciers and downsampled to 100 m to decrease computing time.

To filter and interpolate our DEMs into elevation time series, we empirically characterized the spatial and temporal variance of elevation observations. For this, two global-scale statistical modelling steps relying on a large sampling of the data were performed. One was used to assess the vertical precision of elevation observations and the other to assess their pairwise dependency with varying time lags (Extended Data Fig. 3a, b).

Concomitantly to the variance modelling process described further below, a multi-step outlier filtering was performed to iteratively improve the quality of the DEMs (Extended Data Fig. 1), which itself affects the empirical estimation of the variances. The filtering algorithms consist of a spatial filter, removing elevations outside a topographical maximum and minimum from the TanDEM-X elevations in the pixel surroundings, and a temporal filter propagated from the TanDEM-X elevation at a given pixel through a maximum possible glacier elevation change rate (Extended Data Fig. 3c). These maxima were first conditioned by extreme values (for example, the maximum observed absolute glacier elevation change rate of 50 m yr−1 on HPS12 glacier, Southern Patagonian Icefield53). Later, those were refined by estimating a linear glacier elevation change rate in the surroundings through weighted least squares54.

In our first global-scale statistical modelling step, we identified a heteroscedasticity in elevation measurements (that is, non-uniform variance; Extended Data Fig. 3a). We determined that the elevation measurement error σh varied with the terrain slope55,56 α, the quality of stereo-correlation49,57 q and the individual performance of each DEM’s co-registration48 σc(txy). To empirically quantify this elevation variance, we used ice-free terrain, where no changes in elevation are expected through time, as a proxy for ice-covered terrain. We randomly sampled up to 10,000 ice-free pixels without replacement for each bin of a studied category of terrain (for example, slope) in each 1° × 1° tile and computed the difference to TanDEM-X. We used the median as a robust estimator of the mean and the square of the normalized median absolute deviation (NMAD) as a robust estimator of the variance to mitigate the effects of elevation outliers58. We found that the empirical variances for the slope \({\sigma }_{\alpha }^{2}\) and the quality of stereo-correlation \({\sigma }_{q}^{2}\) were consistent among regions, and used them to condition a model at the global scale to account for the measurement error independently for each elevation observation h(txy):

$${\sigma }_{h}^{2}(t,x,y)={\sigma }_{{\rm{c}}}^{2}(t,x,y)+{\sigma }_{\alpha }^{2}(\alpha ,q)+{\sigma }_{q}^{2}(q).$$
(1)

In our second step of global-scale statistical modelling, we determined the temporal covariance of glacier elevation change (Extended Data Fig. 3b), which serves as our best unbiased estimator to interpolate elevation observations into continuous time series through Gaussian process59 (GP) regressions. To empirically quantify this temporal covariance, we sampled median temporal variograms by the time lag between pairwise elevation observations Δt of ice-covered pixels. We found that the covariance structure could be estimated by the sum of a pairwise linear (PL) kernel, a periodic (exponential sine squared, ESS) kernel, a local (radial basis function, RBF) kernel, and the product of a pairwise linear and local (rational quadratic, RQ) kernel. This sum decomposes the differences of elevation observations with varying time lags into: an underlying linear trend (the PL), a seasonality (the ESS), a proximity at short time lags (the RBF) and a nonlinear trend (the RQ times PL). Empirical covariances showed little variability between regions. We thus conditioned the parameters of the kernels (periodicity ϕp and variance \({\sigma }_{{\rm{p}}}^{2}\) for the ESS; length scale Δtl and variance \({\sigma }_{{\rm{l}}}^{2}\) for the RBF; length scale Δtnl, variance \({\sigma }_{{\rm{nl}}}^{2}\) and scale mixture αnl for the RQ) at the global scale on the basis of our empirical variograms, whereas the PL kernel was determined directly from the observations of each pixel (xy), and thereby described the temporal covariance as:

$${\sigma }_{h}^{2}(x,y,\Delta t)={\rm{PL}}(x,y,\Delta t)+{\rm{ESS}}({\varphi }_{{\rm{p}}},{\sigma }_{{\rm{p}}}^{2},\Delta t)+{\rm{RBF}}(\Delta {t}_{{\rm{l}}},{\sigma }_{{\rm{l}}}^{2},\Delta t)+{\rm{RQ}}(\Delta {t}_{{\rm{nl}}},{\sigma }_{{\rm{nl}}}^{2},{\alpha }_{{\rm{nl}}},\Delta t){\rm{PL}}(x,y,\Delta t)+{\sigma }_{h}^{2}(t,x,y).$$
(2)

By applying GP regression, we iteratively removed observations outside the 20σ, 12σ, 9σ, 6σ and 4σ credible intervals (Extended Data Fig. 3d). Within the same process, elevation time series were then derived at a monthly time step independently for each of the 400 million pixels (xy) falling on or within 10 km of an inventoried glacier22 (Extended Data Fig. 3e). Further details on the variance estimation, filtering and time series methods are available in Supplementary Information and build on refs. 60,61,62,63.

Validation of elevation time series

We retrieved all ICESat (GLAH1464) and IceBridge (IODEM365 and ILAKS1B66) laser and optical elevations intersecting glaciers worldwide from the National Snow and Ice Data Center. IceBridge data are dominated by 1,220,494 Ames Stereo Pipeline67 photogrammetric 0.5–2 m resolution DEMs65 with a typical footprint of 500 m × 500 m that we down-sampled to a resolution of 50 m to limit repeat spatial sampling when comparing to the 100 m resolution of our elevation time series. We linearly interpolated our GP elevation time series in space and time to match the date and centre of each ICESat footprint or IceBridge pixel68 (Extended Data Fig. 4a–c).

We found that regional and seasonal vertical shifts (typically below 2 m) of surface elevation exist, and attribute these differences to snow cover in the TanDEM-X global DEM46 and the presence of seasonally varying snow cover in ASTER, ArcticDEM and REMA DEMs. At the global scale, these shifts do not affect our annual estimates once differenced into elevation changes, verified by the absence of elevation change bias over glaciers (0.001 ± 0.011 m yr−1). We additionally demonstrated that the uncertainties in our elevation time series (credible interval of the GP regression) are conservative (that is, too large by a factor of about two). We reached the same conclusions at the scale of individual RGI 6.0 regions, and also performed these verifications with several additional relevant variables (Extended Data Fig. 4d). In particular, the absence of a bias with glacier elevation denotes our ability to adequately resolve low-texture glacier surfaces in the accumulation area, including flat, high-latitude ice caps. Further details on the validation of the elevation time series are available in Supplementary Information and build on ref. 69.

Integration of elevation into volume changes

We differenced all elevations h into elevation change according to their value of h on 1 January 2000. We integrated the elevation change dh into volume change dV independently for each glacier and time step using a weighted version of the mean local hypsometric method70 with 100-m elevation bins. Weights were derived from the GP elevation change uncertainties, thus ensuring that pixels with a lower vertical precision in a given elevation bin have a smaller impact on the mean elevation change of that bin. Pixels with a 20-year elevation change larger than five times the NMAD from the median elevation change of the elevation bin were removed53. If no valid elevation estimate existed within a given bin, the elevation change was linearly interpolated from adjacent bins, or extrapolated from the closest bins. For retreating lake- and ocean-terminating glaciers, we excluded any loss below water level, because DEMs refer to the water surface and not the poorly known bathymetry in the deglaciated terrain. We note that these losses do not contribute to sea-level rise.

Uncertainty analysis of volume changes

We propagated our uncertainties in elevation change into uncertainties in volume change by assuming that the uncertainty in the mean elevation change \({\sigma }_{\overline{{\rm{d}}h}}\) and the uncertainty in the glacier area σA are independent:

$${\sigma }_{{\rm{d}}V}^{2}={({\sigma }_{\bar{{\rm{d}}h}}A)}^{2}+{({\sigma }_{A}\bar{{\rm{d}}h})}^{2}.$$
(3)

The uncertainty in the mean elevation change \({\sigma }_{\overline{{\rm{d}}h}}\) is highly subject to spatial correlations due to instrument resolution (spatial scale of 0–150 m), uncorrected ASTER instrument noise50 (0–20 km) and the interpolated nature of our elevation time series (0–500 km). The latter spatial correlation term arises from the fact that neighbouring pixels of a given region share similar temporal data gaps, and are hence likely to have similar interpolation biases that correspond to long-range correlations. To empirically quantify these three sources of spatial correlations, we drew spatial variograms of elevation differences between ICESat and our GP elevation time series71 at each ICESat acquisition date. We found that the spatial correlations greatly varied with the time lag Δt to the closest ASTER, ArcticDEM or REMA observation. For each time lag, we estimated the partial sill sk (correlated variance) by fitting a sum of seven spherical variogram models S(dskrk), with d the spatial lag, at ranges rk (correlation lengths) of 0.15 km, 2 km, 5 km, 20 km, 50 km, 200 km and 500 km (Extended Data Fig. 5a, b). To propagate these spatial correlations when integrating glacier volumes, we computed the time lag to the closest ASTER, ArcticDEM or REMA observation for each time step of our elevation time series and for each glacier pixel to estimate s1 to s6. We then used the GP elevation change uncertainties of each glacier pixel to derive s0. Finally, we propagated the pixel-wise uncertainties in elevation change into the uncertainty in the mean elevation change \({\sigma }_{\overline{{\rm{d}}h}}\) by circular integration of the sum of variograms72 over the glacier area A (Extended Data Fig. 5c):

$${\sigma }_{\overline{dh}}^{2}=\frac{1}{A}\mathop{\sum }\limits_{k=0}^{6}{\sigma }_{\bar{d{h}_{k}}}^{2},$$
(4)

where \({\sigma }_{\overline{d{h}_{k}}}^{2}\) is the integrated variance component correlated with range rk:

$${\sigma }_{\bar{d{h}_{k}}}^{2}={\int }_{A}[{s}_{k}-S(d,{s}_{k},{r}_{k})]{\rm{d}}A.$$
(5)

The reliability of the sum of short-range correlations used to account for uncorrected ASTER instrument noise (0–20 km) was further verified by applying empirical methods to ice-free terrain73, and found to yield larger and more realistic uncertainty estimates than the single-range variograms of 0.2–1 km used in previous studies28,53,54,74,75,76. Our maximum correlation length of 500 km accords with known spatial correlations of mass balance estimates77. Further details on the spatial correlation methods are available in Supplementary Information and build on refs. 78,79,80,81,82,83.

For each glacier, we estimated an uncertainty in the area σA based on a buffer84 of 15 m corresponding to the typical resolution of the optical imagery used to derive these outlines37,85,86,87. These uncertainties vary from about 0.1% of the area for large icefields (>1,000 km2) to 50% of the area and above for small isolated glaciers (<0.1 km2).

Validation of volume changes

We retrieved high-resolution DEMs from LiDAR74,88, Pléiades54,89, Satellite Pour l’Observation de la Terre90,91 and aerial photographs92,93 acquired in Alaska, Western North America, Central Europe and High Mountain Asia between 2000 and 2019. We derived precise volume change estimates during specific periods for 588 glaciers covering 3,300 km2 and compared these to our volume time series extracted over the same glaciers and periods. We found no statistically significant bias of mean elevation change (0.03 ± 0.03 m yr−1; Extended Data Fig. 5d). We then validated that our uncertainties, derived from spatially integrated variograms calibrated on ICESat measurements, matched the empirical errors deduced from the comparison (~92% of 95% uncertainty ranges intersect the high-precision volume changes; Extended Data Fig. 5d–f). On average, our 5-year uncertainties at the 95% confidence level are lower than 0.5 m yr−1 for glaciers larger than 1 km2 and conservative for smaller glaciers. We thus validated the reliability of our improved uncertainty approaches for volume change estimation down to the scale of individual glaciers.

Aggregation to regions

We summed volume changes of glaciers per region. To propagate correlated uncertainties among glaciers of the same region, we extended the spatial statistics approach used at the glacier scale. For each time step, glacier-wide correlated uncertainties were propagated again to yield an uncertainty in the mean regional elevation change \({\sigma }_{\bar{{\rm{d}}{h}_{R}}}\). Having been integrated once over a spatial support (from pixel to glacier), the glacier-wide uncertainties can be propagated again (from glacier to regions) directly by a double sum of covariances based on the same describing variograms, following Krige’s relation71:

$${\sigma }_{\overline{{\rm{d}}{h}_{R}}}^{2}=\frac{1}{{A}_{{\rm{R}}}^{2}}\sum _{i}\sum _{j}\mathop{\sum }\limits_{k=0}^{6}[{\sigma }_{\overline{{\rm{d}}{h}_{k,i}}}{\sigma }_{\overline{{\rm{d}}{h}_{k,j}}}-S({G}_{i}-{G}_{j},{\sigma }_{\overline{{\rm{d}}{h}_{k,i}}}{\sigma }_{\overline{{\rm{d}}{h}_{k,j}}},{r}_{k})]\,{A}_{i}{A}_{j},$$
(6)

where i, j are indexes for glaciers in the region, \({\sigma }_{\overline{{\rm{d}}{h}_{k,i}}}\) is the uncertainty in the mean elevation change \({\sigma }_{\overline{{\rm{d}}{h}_{k}}}\) with range rk and sill sk for glacier i, Gi − Gj is the pairwise distance (spatial lag d) between glaciers i and j on the basis of their outline centroids, and AR is the sum of areas Ai of glaciers i in the region.

Conversion to mass changes

We converted volume change into mass change by using a density conversion factor23 of 850 kg m−3 and an uncertainty of 60 kg m−3. This density conversion uncertainty was applied at the scale of RGI 6.0 regions, as if correlated for all glaciers in the entire region—an assumption that yields more conservative estimates than earlier studies29,54. We made this conservative assumption owing to the limited knowledge of spatiotemporal correlations in density conversion. Consequently, our mass change uncertainties might be too large, in particular for regions with the most negative specific-mass change rates (Fig. 3, Extended Data Fig. 5g, h).

Aggregation to global

We summed our regional volume and mass change estimates into global volume and mass change. Assuming independence of the uncertainty in volume and mass changes between RGI 6.0 regions, we summed regional uncertainties quadratically. We report uncertainties in mass change for periods shorter than five years solely for the global or near-global estimates (for example, Fig. 3b) by assuming that the aggregation of largely independent regions leaves limited temporal autocorrelation of density conversion factors. We compare our regional and global mass changes results with global and regional studies listed by the latest IPCC assessment4 as well as additional recent studies28,53,94,95,96 (Supplementary Table 4).

Time-evolving glacier areas

We accounted for temporal changes in glacier areas when deriving regional or global time series of specific (area-scaled) mass balances or mean elevation change (specific-volume change). We assumed a linear change through time, calibrated on time-evolving glacier outlines of each RGI 6.0 region21. Over the 20-year study period, these time-evolving glacier areas correspond to a nearly 10% decrease of glacier areas around the globe—a non-negligible change when assessing mean elevation change rates. To account for this, we added an additional uncertainty in the time-evolving glacier area at each time step of 1% of the regional area at that time step.

Observed sea-level rise

We derived global mean sea-level trends from a recent study24 with time series extended to match our period of study of 2000–2019, yielding an estimate of sea-level rise of 3.56 ± 0.4 mm yr−1 with an acceleration of 0.15 ± 0.08 mm yr−2. For conversion, we assumed that 361.8 Gt of water-equivalent mass loss amounted to 1 mm of sea-level rise.

Acceleration

Glacier mass change acceleration and its uncertainties were derived from weighted least squares on the 5-year elevation and mass change rates (that is, 2000–2004, 2005–2009, 2010–2014 and 2015–2019), propagating their related uncertainties as independent. Although shorter timescales and smaller spatial domains are affected by temporal autocorrelation, we assumed the 5-year estimates at the global or near-global scale (that is, excluding peripheral glaciers) as temporally uncorrelated. This assumption is supported by timescales described for density conversion factors23, by the validation of our elevation time series with ICESat and IceBridge, and relies on the billions of globally distributed surface elevation observations, leading to large independent and repeat sampling over 5-year periods (Extended Data Table 2).

Distinction between glaciers and ice sheets

When comparing our results to ice sheet studies, we avoided double-counting contributions from peripheral glaciers by subtracting part of our own estimate for RGI 6.0 regions 5 and 19 to ice sheet estimates from IMBIE7,8. Because IMBIE estimates are a weighted mean of three ensemble estimates where half includes peripheral glaciers, the other half does not (gravimetric studies include peripheral glaciers, altimetric studies exclude peripheral glaciers, and input–output studies do both), we assumed that subtracting half of our estimates for the peripheral glaciers was most adequate. Notably, applying this subtraction leads to better agreement of GIS and AIS estimates between IMBIE and a recent study9 over the period 2003–2018 (Table 1).

Temperature and precipitation analysis

We analysed ERA5 precipitation and temperature97 at both annual and seasonal scales. For the latter scale, we considered only winter precipitation and summer temperature. We found similar decadal patterns at both annual and seasonal scales, and thus present annual changes (Fig. 4) to avoid the latitudinal ambiguity of glaciological definitions of seasons. Temperature change was extracted at 700 hPa (about 3,100 m above sea level) to minimize variations in air temperature affected by differences in land surface class at the 0.125° nominal resolution of the ERA5 reanalysis. To estimate trends of annual precipitation and temperature over 2000–2019, we derived ordinary least-squares trends for each ERA5 grid cell containing glaciers. We then area-weighted the global trend by the glacierized area of each grid cell. We detected a small increase in precipitation at the global scale (4.0% in 20 years) and over glaciers (6.2% in 20 years), coherent with the amplification of the global water cycle in a warming world near the Clausius–Clapeyron rate98. The sensitivity of mass change to air temperature was computed by dividing the specific-mass change acceleration by the temperature increase over glacierized terrain for the period 2000–2019.

Data availability

Global, regional, tile and per-glacier elevation and mass change time series, elevation change maps for 5-, 10- and 20-year periods at 100 m resolution, and tables in this article are publicly available at https://doi.org/10.6096/13. Source data are provided with this paper.

Code availability

The code developed for the global processing and analysis of all data, and to generate figures and tables in this article, is publicly available at https://github.com/rhugonnet/ww_tvol_study. Code concomitantly developed for processing ASTER data is available as the Python package pymmaster at https://github.com/luc-girod/MMASTER-workflows (with supporting documentation at https://mmaster-workflows.readthedocs.io) and for processing DEM time series as the Python package pyddem at https://github.com/iamdonovan/pyddem (with supporting documentation at https://pyddem.readthedocs.io).

References

  1. 1.

    Pritchard, H. D. Asia’s shrinking glaciers protect large populations from drought stress. Nature 569, 649–654 (2019).

    ADS  CAS  PubMed  Google Scholar 

  2. 2.

    WCRP Global Sea Level Budget Group. Global sea-level budget 1993–present. Earth Syst. Sci. Data 10, 1551–1590 (2018).

    ADS  Google Scholar 

  3. 3.

    Stoffel, M. & Huggel, C. Effects of climate change on mass movements in mountain environments. Prog. Phys. Geogr. 36, 421–439 (2012).

    Google Scholar 

  4. 4.

    IPCC. IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (eds Pörtner, H. O. et al.) (IPCC, 2019).

  5. 5.

    Gardner, A. et al. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009. Science 340, 852–857 (2013).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  6. 6.

    Nerem, R. S. et al. Climate-change-driven accelerated sea-level rise detected in the altimeter era. Proc. Natl Acad. Sci. USA 115, 2022–2025 (2018).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  7. 7.

    IMBIE Team. Mass balance of the Greenland Ice Sheet from 1992 to 2018. Nature 579, 233–239 (2020).

    ADS  Google Scholar 

  8. 8.

    IMBIE team. Mass balance of the Antarctic Ice Sheet from 1992 to 2017. Nature 558, 219–222 (2018).

    ADS  Google Scholar 

  9. 9.

    Smith, B. et al. Pervasive ice sheet mass loss reflects competing ocean and atmosphere processes. Science 368, 1239–1242 (2020).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Kääb, A., Berthier, E., Nuth, C., Gardelle, J. & Arnaud, Y. Contrasting patterns of early twenty-first-century glacier mass change in the Himalayas. Nature 488, 495–498 (2012).

    ADS  PubMed  PubMed Central  Google Scholar 

  11. 11.

    Kulp, S. A. & Strauss, B. H. New elevation data triple estimates of global vulnerability to sea-level rise and coastal flooding. Nat. Commun. 10, 4844 (2019); author correction 10, 5752 (2019).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Immerzeel, W. W. et al. Importance and vulnerability of the world’s water towers. Nature 577, 364–369 (2020).

    CAS  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Marzeion, B., Cogley, J. G., Richter, K. & Parkes, D. Attribution of global glacier mass loss to anthropogenic and natural causes. Science 345, 919–921 (2014).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Huss, M. & Hock, R. Global-scale hydrological response to future glacier mass loss. Nat. Clim. Chang. 8, 135–140 (2018).

    ADS  Google Scholar 

  15. 15.

    IPCC. Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part B: Regional Aspects (Cambridge University Press, 2014).

  16. 16.

    Cauvy-Fraunié, S. & Dangles, O. A global synthesis of biodiversity responses to glacier retreat. Nat. Ecol. Evol. 3, 1675–1685 (2019).

    PubMed  PubMed Central  Google Scholar 

  17. 17.

    World Glacier Monitoring Service (WGMS). Fluctuations of Glaciers Database https://wgms.ch/data_databaseversions/ (2019).

  18. 18.

    Bamber, J. L., Westaway, R. M., Marzeion, B. & Wouters, B. The land ice contribution to sea level during the satellite era. Environ. Res. Lett. 13, 063008 (2018); corrigendum 13, 099502 (2018).

    ADS  Google Scholar 

  19. 19.

    Wouters, B., Gardner, A. S. & Moholdt, G. Global glacier mass loss during the GRACE satellite mission (2002–2016). Front. Earth Sci. 7, 96 (2019).

    ADS  Google Scholar 

  20. 20.

    Ciracì, E., Velicogna, I. & Swenson, S. Continuity of the mass loss of the world’s glaciers and ice caps from the GRACE and GRACE Follow-On missions. Geophys. Res. Lett. 47, 226 (2020).

    Google Scholar 

  21. 21.

    Zemp, M. et al. Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. Nature 568, 382–386 (2019).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  22. 22.

    RGI Consortium. Randolph Glacier Inventory – A Dataset of Global Glacier Outlines. Technical Report https://www.glims.org/RGI/00_rgi60_TechnicalNote.pdf (Global Land Ice Measurements from Space, 2017).

  23. 23.

    Huss, M. Density assumptions for converting geodetic glacier volume change to mass change. Cryosphere 7, 877–887 (2013).

    ADS  Google Scholar 

  24. 24.

    Ablain, M. et al. Uncertainty in satellite estimates of global mean sea-level changes, trend and acceleration. Earth Syst. Sci. Data 11, 1189–1202 (2019).

    ADS  Google Scholar 

  25. 25.

    Velicogna, I. et al. Continuity of ice sheet mass loss in Greenland and Antarctica from the GRACE and GRACE Follow-On missions. Geophys. Res. Lett. 47, L11501 (2020).

    Google Scholar 

  26. 26.

    Larsen, C. F. et al. Surface melt dominates Alaska glacier mass balance. Geophys. Res. Lett. 42, 5902–5908 (2015).

    ADS  Google Scholar 

  27. 27.

    Blazquez, A. et al. Exploring the uncertainty in GRACE estimates of the mass redistributions at the Earth surface: implications for the global water and sea level budgets. Geophys. J. Int. 215, 415–430 (2018).

    ADS  Google Scholar 

  28. 28.

    Shean, D. E. et al. A systematic, regional assessment of High Mountain Asia glacier mass balance. Front. Earth Sci. 7, 363 (2020).

    ADS  Google Scholar 

  29. 29.

    Braun, M. H. et al. Constraining glacier elevation and mass changes in South America. Nat. Clim. Chang. (2019).

  30. 30.

    Dehecq, A. et al. Elevation changes inferred from TanDEM-X data over the Mont-Blanc area: impact of the X-band interferometric bias. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 9, 3870–3882 (2016).

    ADS  Google Scholar 

  31. 31.

    Sandberg Sørensen, L. et al. 25 years of elevation changes of the Greenland Ice Sheet from ERS, Envisat, and CryoSat-2 radar altimetry. Earth Planet. Sci. Lett. 495, 234–241 (2018).

    ADS  Google Scholar 

  32. 32.

    Bevis, M. et al. Accelerating changes in ice mass within Greenland, and the ice sheet’s sensitivity to atmospheric forcing. Proc. Natl Acad. Sci. USA 116, 1934–1939 (2019).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  33. 33.

    Garreaud, R. D. et al. The Central Chile Mega Drought (2010–2018): a climate dynamics perspective. Int. J. Climatol. 40, 421–439 (2020).

    Google Scholar 

  34. 34.

    Raper, S. C. B. & Braithwaite, R. J. Low sea level rise projections from mountain glaciers and icecaps under global warming. Nature 439, 311–313 (2006).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Parkes, D. & Marzeion, B. Twentieth-century contribution to sea-level rise from uncharted glaciers. Nature 563, 551–554 (2018).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Becker, J. J. et al. Global bathymetry and elevation data at 30 arc seconds resolution: SRTM30_PLUS. Mar. Geod. 32, 355–371 (2009).

    Google Scholar 

  37. 37.

    Tielidze, L. G. & Wheate, R. D. The Greater Caucasus glacier inventory (Russia, Georgia and Azerbaijan). Cryosphere 12, 81–94 (2018).

    ADS  Google Scholar 

  38. 38.

    Dunse, T. et al. Glacier-surge mechanisms promoted by a hydro-thermodynamic feedback to summer melt. Cryosphere 9, 197–215 (2015).

    ADS  Google Scholar 

  39. 39.

    McMillan, M. et al. Rapid dynamic activation of a marine-based Arctic ice cap: ice cap dynamic activation. Geophys. Res. Lett. 41, 8902–8909 (2014).

    ADS  Google Scholar 

  40. 40.

    Nuth, C. et al. Dynamic vulnerability revealed in the collapse of an Arctic tidewater glacier. Sci. Rep. 9, 5541 (2019).

    ADS  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Howat, I. M., Negrete, A. & Smith, B. E. The Greenland Ice Mapping Project (GIMP) land classification and surface elevation data sets. Cryosphere 8, 1509–1518 (2014).

    ADS  Google Scholar 

  42. 42.

    Fretwell, P. et al. Bedmap2: improved ice bed, surface and thickness datasets for Antarctica. Cryosphere 7, 375–393 (2013).

    ADS  Google Scholar 

  43. 43.

    NASA/METI/AIST/Japan Spacesystems & U.S./Japan ASTER Science Team. ASTER Level 1A Data Set – Reconstructed, Unprocessed Instrument Data. 2001, NASA EOSDIS Land Processes DAAC, 2001); https://doi.org/10.5067/ASTER/AST_L1A.003.

  44. 44.

    Porter, C. et al. ArcticDEM (Harvard Dataverse, 2018); https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/OHHUKH.

  45. 45.

    Howat, I. M., Porter, C., Smith, B. E., Noh, M.-J. & Morin, P. The reference elevation model of Antarctica. Cryosphere 13, 665–674 (2019).

    ADS  Google Scholar 

  46. 46.

    Rizzoli, P. et al. Generation and performance assessment of the global TanDEM-X digital elevation model. ISPRS J. Photogramm. Remote Sens. 132, 119–139 (2017).

    ADS  Google Scholar 

  47. 47.

    Vassilaki, D. I. & Stamos, A. A. TanDEM-X DEM: comparative performance review employing LIDAR data and DSMs. ISPRS J. Photogramm. Remote Sens. 160, 33–50 (2020).

    Google Scholar 

  48. 48.

    Nuth, C. & Kääb, A. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. Cryosphere 5, 271–290 (2011).

    ADS  Google Scholar 

  49. 49.

    Rupnik, E., Daakir, M. & Pierrot Deseilligny, M. MicMac – a free, open-source solution for photogrammetry. Open Geospat. Data Softw. Stand. 2, 14 (2017).

    Google Scholar 

  50. 50.

    Girod, L., Nuth, C., Kääb, A., McNabb, R. & Galland, O. MMASTER: improved ASTER DEMs for elevation change monitoring. Remote Sens. 9, 704 (2017).

    Google Scholar 

  51. 51.

    Wales, D. J. & Doye, J. P. K. Global optimization by basin-hopping and the lowest energy structures of Lennard–Jones clusters containing up to 110 atoms. J. Phys. Chem. A 101, 5111–5116 (1997).

    CAS  Google Scholar 

  52. 52.

    Noh, M.-J. & Howat, I. M. The surface extraction from TIN based Search-space Minimization (SETSM) algorithm. ISPRS J. Photogramm. Remote Sens. 129, 55–76 (2017).

    ADS  Google Scholar 

  53. 53.

    Dussaillant, I. et al. Two decades of glacier mass loss along the Andes. Nat. Geosci. 12, 802–808 (2019); author correction 13, 711 (2020).

    ADS  CAS  Google Scholar 

  54. 54.

    Brun, F., Berthier, E., Wagnon, P., Kääb, A. & Treichler, D. A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000–2016. Nat. Geosci. 10, 668–673 (2017); author correction 11, 543 (2018).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Toutin, T. Three-dimensional topographic mapping with ASTER stereo data in rugged topography. IEEE Trans. Geosci. Remote Sens. 40, 2241–2247 (2002).

    ADS  Google Scholar 

  56. 56.

    Lacroix, P. Landslides triggered by the Gorkha earthquake in the Langtang valley, volumes and initiation processes. Earth Planets Space 68, 1–10 (2016).

    Google Scholar 

  57. 57.

    Shean, D. E. et al. An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery. ISPRS J. Photogramm. Remote Sens. 116, 101–117 (2016).

    ADS  Google Scholar 

  58. 58.

    Höhle, J. & Höhle, M. Accuracy assessment of digital elevation models by means of robust statistical methods. ISPRS J. Photogramm. Remote Sens. 64, 398–406 (2009).

    ADS  Google Scholar 

  59. 59.

    Williams, C. K. I. & Rasmussen, C. E. Gaussian Processes for Machine Learning Vol. 2 (MIT Press, 2006).

  60. 60.

    Schiefer, E., Menounos, B. & Wheate, R. Recent volume loss of British Columbian glaciers, Canada. Geophys. Res. Lett. (2007).

  61. 61.

    Nuimura, T., Fujita, K., Yamaguchi, S. & Sharma, R. R. Elevation changes of glaciers revealed by multitemporal digital elevation models calibrated by GPS survey in the Khumbu region, Nepal Himalaya, 1992–2008. J. Glaciol. 58, 648–656 (2012).

    ADS  Google Scholar 

  62. 62.

    Willis, M. J., Melkonian, A. K., Pritchard, M. E. & Rivera, A. Ice loss from the Southern Patagonian Ice Field, South America, between 2000 and 2012. Geophys. Res. Lett. 39, L17501 (2012).

    ADS  Google Scholar 

  63. 63.

    Pedregosa, F. et al. Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830 (2011).

    MathSciNet  MATH  Google Scholar 

  64. 64.

    Zwally, H. J., Schutz, R., Hancock, D. & Dimarzio, J. GLAS/ICESat L2 Global Land Surface Altimetry Data (HDF5), Version 34 (NASA Snow and Ice Data Center, 2014); https://nsidc.org/data/GLAH14.

  65. 65.

    Alexandrov, O., McMichael, S. & Beyer., R. A. IceBridge DMS L3 Ames Stereo Pipeline Photogrammetric DEM, Version 1 (accessed 1 June 2019); https://nsidc.org/data/IODEM3/versions/1.

  66. 66.

    Larsen, C. IceBridge UAF Lidar Scanner L1B Geolocated Surface Elevation Triplets, Version 1 (accessed 20 February 2020); https://nsidc.org/data/ILAKS1B/versions/1.

  67. 67.

    Beyer, R. A., Alexandrov, O. & McMichael, S. The Ames Stereo Pipeline: NASA’s open source software for deriving and processing terrain data. Earth Space Sci. 5, 537–548 (2018).

    ADS  Google Scholar 

  68. 68.

    Harding, D. J. ICESat waveform measurements of within-footprint topographic relief and vegetation vertical structure. Geophys. Res. Lett. 32, L21S10 (2005).

    Google Scholar 

  69. 69.

    Gardelle, J., Berthier, E. & Arnaud, Y. Impact of resolution and radar penetration on glacier elevation changes computed from DEM differencing. J. Glaciol. 58, 419–422 (2012).

    ADS  Google Scholar 

  70. 70.

    McNabb, R., Nuth, C., Kääb, A. & Girod, L. Sensitivity of glacier volume change estimation to DEM void interpolation. Cryosphere 13, 895–910 (2019).

    ADS  Google Scholar 

  71. 71.

    Cressie, N. A. C. Statistics for Spatial Data Vol. 4, 613–617 (Wiley, 1993).

  72. 72.

    Rolstad, C., Haug, T. & Denby, B. Spatially integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis: application to the western Svartisen ice cap, Norway. J. Glaciol. 55, 666–680 (2009).

    ADS  Google Scholar 

  73. 73.

    Dehecq, A. et al. Automated processing of declassified KH-9 Hexagon satellite images for global elevation change analysis since the 1970s. Front. Earth Sci. 8, 566802 (2020).

    Google Scholar 

  74. 74.

    Menounos, B. et al. Heterogeneous changes in western North American glaciers linked to decadal variability in zonal wind strength. Geophys. Res. Lett. 46, 200–209 (2018).

    ADS  Google Scholar 

  75. 75.

    Howat, I. M., Smith, B. E., Joughin, I. & Scambos, T. A. Rates of southeast Greenland ice volume loss from combined ICESat and ASTER observations. Geophys. Res. Lett. 35, L17505 (2008).

    ADS  Google Scholar 

  76. 76.

    Wang, D. & Kääb, A. Modeling glacier elevation change from DEM time series. Remote Sens. 7, 10117–10142 (2015).

    ADS  Google Scholar 

  77. 77.

    Cogley, J. G. & Adams, W. P. Mass balance of glaciers other than the ice sheets. J. Glaciol. 44, 315–325 (1998).

    ADS  CAS  Google Scholar 

  78. 78.

    Journel, A. G. & Huijbregts, C. J. Mining Geostatistics Vol. 600 (Academic Press, 1978).

  79. 79.

    Webster, R. & Oliver, M. A. Geostatistics for Environmental Scientists (John Wiley & Sons, 2007).

  80. 80.

    Gräler, B., Pebesma, E. & Heuvelink, G. Spatio-temporal interpolation using gstat. R J. 8, 204 (2016).

    Google Scholar 

  81. 81.

    Mälicke, M. & Schneider, H. D. Scikit-GStat 0.2.6: A Scipy Flavored Geostatistical Analysis Toolbox Written in Python (2019); https://zenodo.org/record/3531816#.YFsJ737Le00.

  82. 82.

    Dussaillant, I., Berthier, E. & Brun, F. Geodetic mass balance of the Northern Patagonian Icefield from 2000 to 2012 using two independent methods. Front. Earth Sci. 6, 8 (2018).

    ADS  Google Scholar 

  83. 83.

    Berthier, E., Scambos, T. A. & Shuman, C. A. Mass loss of Larsen B tributary glaciers (Antarctic Peninsula) unabated since 2002. Geophys. Res. Lett. 39, L13501 (2012).

    ADS  Google Scholar 

  84. 84.

    Granshaw, F. D. & Fountain, A. G. Glacier change (1958–1998) in the North Cascades National Park Complex, Washington, USA. J. Glaciol. 52, 251–256 (2006).

    ADS  CAS  Google Scholar 

  85. 85.

    Pfeffer, W. et al. The Randolph Glacier Inventory: a globally complete inventory of glaciers. J. Glaciol. 60, 537–552 (2014).

    ADS  Google Scholar 

  86. 86.

    Rastner, P. et al. The first complete inventory of the local glaciers and ice caps on Greenland. Cryosphere 6, 1483–1495 (2012).

    ADS  Google Scholar 

  87. 87.

    Bolch, T., Menounos, B. & Wheate, R. Landsat-based inventory of glaciers in western Canada, 1985–2005. Remote Sens. Environ. 114, 127–137 (2010).

    ADS  Google Scholar 

  88. 88.

    Pelto, B. M., Menounos, B. & Marshall, S. J. Multi-year evaluation of airborne geodetic surveys to estimate seasonal mass balance, Columbia and Rocky Mountains, Canada. Cryosphere 13, 1709–1727 (2019).

    ADS  Google Scholar 

  89. 89.

    Wagnon, P. et al. Seasonal and annual mass balances of Mera and Pokalde glaciers (Nepal Himalaya) since 2007. Cryosphere 7, 1769–1786 (2013).

    ADS  Google Scholar 

  90. 90.

    Berthier, E., Schiefer, E., Clarke, G. K. C., Menounos, B. & Rémy, F. Contribution of Alaskan glaciers to sea-level rise derived from satellite imagery. Nat. Geosci. 3, 92–95 (2010).

    ADS  CAS  Google Scholar 

  91. 91.

    Berthier, E., Cabot, V., Vincent, C. & Six, D. Decadal region-wide and glacier-wide mass balances derived from multi-temporal ASTER Satellite Digital Elevation Models. Validation over the Mont-Blanc area. Front. Earth Sci. 4, 63 (2016).

    ADS  Google Scholar 

  92. 92.

    Glacier Monitoring Switzerland. Swiss Glacier Volume Change, Release 2018 (2018); https://doi.glamos.ch/data/volumechange/volumechange_2018_r2018.html.

  93. 93.

    Bauder, A., Funk, M. & Huss, M. Ice-volume changes of selected glaciers in the Swiss Alps since the end of the 19th century. Ann. Glaciol. 46, 145–149 (2007).

    ADS  Google Scholar 

  94. 94.

    Davaze, L., Rabatel, A., Dufour, A., Hugonnet, R. & Arnaud, Y. Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016. Front. Earth Sci. 8, 149 (2020).

    ADS  Google Scholar 

  95. 95.

    Schuler, T. V. et al. Reconciling Svalbard Glacier mass balance. Front. Earth Sci. 8, 523646 (2020).

    Google Scholar 

  96. 96.

    Aðalgeirsdóttir, G. et al. Glacier Changes in Iceland From ~1890 to 2019. Front. Earth Sci. 8, 520 (2020).

    ADS  Google Scholar 

  97. 97.

    Hersbach, H. & Dee, D. ERA5 reanalysis is in production. ECMWF Newsl. 147, 5–6 (2016).

    Google Scholar 

  98. 98.

    Skliris, N., Zika, J. D., Nurser, G., Josey, S. A. & Marsh, R. Global water cycle amplifying at less than the Clausius–Clapeyron rate. Sci. Rep. 6, 38752 (2016).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  99. 99.

    Sakakibara, D., Sugiyama, S., Sawagaki, T., Marinsek, S. & Skvarca, P. Rapid retreat, acceleration and thinning of Glaciar Upsala, Southern Patagonia Icefield, initiated in 2008. Ann. Glaciol. 54, 131–138 (2013).

    ADS  Google Scholar 

  100. 100.

    Farr, T. G. et al. The Shuttle Radar Topography Mission. Rev. Geophys. 45, RG2004 (2007).

    ADS  Google Scholar 

Download references

Acknowledgements

We thank C. Porter for discussions on ArcticDEM and REMA DEMs, B. Meyssignac for comments on sea-level rise and A. Dehecq for input on the presentation of the manuscript. The GLIMS initiative (in particular J. Kargel and B. Raup) allowed the population of a vast archive of ASTER stereo images over glaciers. Hakai Institute and the University of Northern British Columbia provided computational resources for processing ASTER stereo imagery. SPOT6/7 data were obtained from GEOSUD (ANR-10-EQPX-20, programme ‘Investissements d’Avenir’). ArcticDEM DEMs were provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691 and 1542736 and REMA DEMs were provided by the Byrd Polar and Climate Research Center and the Polar Geospatial Center under NSF-OPP awards 1543501, 1810976, 1542736, 1559691, 1043681, 1541332, 0753663, 1548562, 1238993 and NASA award NNX10AN61G. Computer time was provided through a Blue Waters Innovation Initiative. DEMs were produced using data from DigitalGlobe, Inc. R.H. acknowledges a fellowship from the University of Toulouse. E.B. acknowledges support from the French Space Agency (CNES) through ISIS and TOSCA programmes. R.M., C.N., L.G. and A.K. acknowledge support by ESA through Glaciers_cci and EE10 (4000109873/14/I-NB, 4000127593/19/I-NS, 4000127656/19/NL/FF/gp), and by the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC grant agreement number 320816. B.M. acknowledges funding from the National Sciences and Engineering Research Council of Canada, the Canada Research Chairs Program, the Tula Foundation and Global Water Futures. R.H., D.F. and M.H. acknowledge funding from the Swiss National Science Foundation, grant number 184634.

Author information

Affiliations

Authors

Contributions

E.B. and R.H. designed the study with contributions from D.F., M.H. and B.M. L.G., C.N., R.M. and A.K. developed ASTER bias-correction methods. R.H. and R.M. developed glacier elevation GP methods. R.H. implemented spatial statistics methods with inputs from F.B. B.M. assembled and analysed ERA5 data. R.H. performed the processing and analysis of all data with main inputs from E.B., as well as R.M., B.M., D.F., M.H., I.D. and F.B. All authors interpreted the results. R.H. led the writing of the paper and all other co-authors contributed.

Corresponding author

Correspondence to Romain Hugonnet.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Peer review information Nature thanks Beata Csatho, Thomas Frederikse, Michael Willis and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Flow chart of the methodology.

Flow diagram describing the processing steps from satellite imagery to global glacier mass change time series. Processing steps correspond to sections in Methods.

Extended Data Fig. 2 Spatial and temporal coverage of ASTER, ArcticDEM and REMA DEMs.

ac, Spatial distribution of DEMs as a strip count for ArcticDEM strips above 50° N (a), ASTER DEM strips (b) and REMA strips below 50° S (c), shown on top of a world hillshade36. 67,986 ArcticDEM and 9,369 REMA strips are counted before co-registration to TanDEM-X. This later reduces their number to 40,391 and 3,456, respectively, owing to the limited stable terrain in polar regions. d, Temporal distribution of the strip count as a bi-mensual histogram from January 2000 to December 2019. We note that ArcticDEM and REMA strip footprints (15 km × 50 km) are generally much smaller than ASTER DEM strip footprints (180 km × 60 km).

Extended Data Fig. 3 Elevation time series estimation.

ae, Empirical and modelled elevation measurement error (a) and temporal covariance of glacier elevation (b) estimated globally. These are used to condition the filtering (c, d) and elevation time series estimation (e) of elevation observations, illustrated here for a 100 m × 100 m pixel on the ablation area of Upsala, where a strong nonlinear elevation loss occurred99. a, Squared measurement error, estimated by the squared NMAD of elevation differences to TanDEM-X on stable terrain as a function of terrain slope and of quality of stereo-correlation. We express the quality of stereo-correlation as a percentage ranging from 0% for poor correlations to 100% for good correlations. b, Variance between pairwise glacier elevations in time, or temporal variogram. The empirical temporal variogram is derived from the aggregated median of variances binned by time lags of 0.25 yr. Here, pixels were selected on glacierized terrain showing a linear trend of elevation change (estimated from weighted least squares) between −1.5 and −1.0 m yr−1. The median of the linear trend at these locations (−1.2 m yr−1) was directly used to derive the linear model (orange), which has a quadratic variance. The other models are calibrated so that their sum (dashed black line) matches the empirical variogram. c, Spatial and temporal filtering by conditioning a maximum linear elevation change rate from the neighbouring TanDEM-X elevations (see Supplementary Information for further details). d, Filtering by successive GP regression fits for credible intervals of size 20σ, 12σ, 9σ, 6σ and 4σ. e, Elevation time series of final GP regression after the removal of outliers.

Extended Data Fig. 4 Validation of elevation time series and uncertainties to ICESat and IceBridge.

ad, ICESat64 and IceBridge65,66 measurements compared to our surface elevation time series over glacierized terrain in the Saint-Elias Mountains, Alaska (ac) and at the global scale (d). b, Absolute z-scores (white to purple) are shown on top of the 2000–2019 surface elevation change. z-scores correspond to elevation differences to ICESat (dashed outlines) or IceBridge (solid outlines), standardized by our time series uncertainty. c, Time series for a 100 m × 100 m pixel extracted on the tongue of Agassiz Glacier with neighbouring ICESat and IceBridge elevation differences for demonstration purposes. d, Summary of global validation statistics for categories of time, season, region, elevation, observation time lag and total elevation change, with density distributions of measurements for ICESat (light grey) and IceBridge (dark grey). Mean elevation differences, subject to snow-cover biases, are shown only by region (summer mean) and by two-month seasonal component (difference to the annual mean) for each hemisphere.

Extended Data Fig. 5 Uncertainty analysis of volume changes and validation using high-resolution DEMs.

ah, Spatial correlation of elevations between the GP time series and ICESat with the time lag to the closest ASTER, ArcticDEM or REMA observation (a, b), propagation of correlations into specific-volume change uncertainties (c), validation of volume change estimates and uncertainties to high-resolution volume changes extracted over the same 588 glaciers and periods (df) and contribution from all uncertainty sources to the 2000–2019 specific-mass change estimates (g, h). a, An empirical spatial variogram is shown and fitted with a sum of spherical models at correlation lengths of 0.15, 2, 5, 20, 50, 200 and 500 km for elevation differences sampled at 720 days (2 years) from the closest observation. b, Spatially correlated variances as a function of the time lag to the closest observation. The model for the variance used during uncertainty propagation is shown in plain lines (sum of quadratic and squared sinusoidal functions optimized by least squares). c, Propagation of elevation change uncertainties to volume change uncertainties with varying glacier area. As this computation is specific to the time lag of each pixel to the closest observation, for each glacier, at each time step, c refers to an example. The spatial correlations are computed for a time lag to the closest observation, representing the average of our study, of 0–1 yr for 50% of observations, 1–2 yr for 20% of observations, 2–3 yr for 20% of observations and 3–4 yr for 10% of observations. We assume a mean pixel-wise uncertainty of 10 m and simplify by considering only the first step of integration over a continuous glacierized area (equation (5)). This assumption leads to slightly larger contributions from short-range correlations than with further propagation to the second propagation step between discontinuous glaciers (equation (6)). Uncertainties are largely dominated by short- to long-range spatial correlations. d, Comparison of specific-volume changes per glacier with 1σ uncertainties. The mean of differences in estimates over all glaciers does not statistically differ from zero. e, f, Theoretical and empirical 1σ uncertainties, and their evolution with glacier size. The theoretical uncertainty is the mean of per-glacier uncertainties derived from spatially integrated variograms and the empirical uncertainty is the NMAD of the difference between high-resolution and GP estimates. g, h, Propagation of uncertainty sources to specific-mass changes for each RGI 6.0 region, and all glaciers with and without the Greenland Periphery and the Antarctic and Subantarctic, which are magnified in h. Uncertainties are largely dominated by the volume-to-mass conversion uncertainties globally, and by uncertainties in glacier outlines for regions with a relevant share of small glaciers.

Extended Data Fig. 6 Two decades of elevation change over various regions.

ah, Elevation change of glaciers between 2000 and 2019 in Coropuna, Peru (a), Pamir Mountains (b), Iceland (c), Karakoram Mountains (d), European Alps (e), Southern Alps, New Zealand (f), West Greenland (note the rotated orientation of map) (g) and Svalbard (h). Except for Svalbard, glacier outlines displayed are from the RGI 6.0. In the background is shown a hillshade derived from several sources36,46,100. In Svalbard, outlines have been updated to include the massive surges of Austfonna Basin 338,39 in the northeast and Nathorstbreen in the southwest40, indicated by blue arrows.

Extended Data Fig. 7 Global evolution of 5-year thinning rates.

ad, Mean elevation change rates aggregated by tiles of 1° × 1° for the periods 2000–2004 (a), 2005–2009 (b), 2010–2014 (c) and 2015–2019 (d). The tile area is inversely scaled to the squared 95% confidence interval of the mean elevation change in the tile, and tiles are coloured with mean elevation change rates, on top of a world hillshade36. The minimum tile area is 10% for a 95% confidence interval larger than 2 m yr−1 and tiles are displayed at full size for a 95% confidence interval smaller than 0.5 m yr−1. Region labelling refers to that of Fig. 2. The acceleration of thinning brings the Karakoram anomaly to its apparent end.

Extended Data Table 1 Regional rates of glacier elevation and mass change from 2000 to 2019
Extended Data Table 2 Regional data coverage of elevation time series from 2000 to 2019
Extended Data Table 3 Regional rates of land- and marine-terminating glaciers in maritime regions

Supplementary information

Supplementary Information

This file contains the Supplementary Methods, Supplementary Discussion, Supplementary Figures 1–9 and Supplementary Tables 1–3.

Supplementary Table 4

Comparison to IPCC SROCC Table 2A.1 with estimates from this study for periods 2006–2015 and 2000–2019 and recent regional studies (blue). An additional decimal is shown for mass balance rates in Gt yr−1 and mm SLE yr−1. Regions are shown in SROCC ordering. SROCC estimates were combined using the most suitable studies (methods ‘xxx’ and ‘x’).

Peer Review File

Source data

Rights and permissions

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Hugonnet, R., McNabb, R., Berthier, E. et al. Accelerated global glacier mass loss in the early twenty-first century. Nature 592, 726–731 (2021). https://doi.org/10.1038/s41586-021-03436-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41586-021-03436-z

Further reading

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing