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

Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I.,
Brun, F., & Kääb, A. (2021). Accelerated global glacier mass loss in the early twenty-first century. Nature,
592(7856), 726–731. https://doi.org/10.1038/s41586-021-03436-z

Link to publication record in Ulster University Research Portal

Published in:
Nature

Publication Status:
Published (in print/issue): 29/04/2021

DOI:
10.1038/s41586-021-03436-z

Document Version
Author Accepted version

General rights
Copyright for the publications made accessible via Ulster University's Research Portal is retained by the author(s) and / or other copyright
owners and it is a condition of accessing these publications that users recognise and abide by the legal requirements associated with these
rights.

Take down policy

The Research Portal is Ulster University's institutional repository that provides access to Ulster's research outputs. Every effort has been
made to ensure that content in the Research Portal does not infringe any person's rights, or applicable UK laws. If you discover content in
the Research Portal that you believe breaches copyright or violates any law, please contact pure-support@ulster.ac.uk.

Download date: 17/05/2024

 Accelerated global glacier mass loss in the early
twenty-first century
Romain Hugonnet1,2,3*, Robert McNabb4,5, Etienne Berthier1, Brian Menounos6,7, Christopher
Nuth5,8,9, Luc Girod5, Daniel Farinotti2,3, Matthias Huss2,3,10, Ines Dussaillant1,11, Fanny Brun12
and Andreas Kääb 5

1
LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, Toulouse, France.
2
Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zürich, Zürich, Switzerland.
3
Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), Birmensdorf, Switzerland.
4
School of Geography and Environmental Sciences, Ulster University, Coleraine, United Kingdom.
5
Department of Geosciences, University of Oslo, Oslo, Norway.
6
Natural Resources and Environmental Studies Institute and Geography, University of Northern British
Columbia, Prince George, British Columbia, Canada.
7
Hakai Institute, Quadra Island, Canada.
8
The Norwegian Defense Research Establishment, Kjeller, Norway.
9
Heimdal Satellite Technologies, Oslo, Norway.
10
Department of Geosciences, University of Fribourg, Fribourg, Switzerland.
11
Department of Geography, University of Zurich, Zurich, Switzerland.
12
IGE, Université Grenoble Alpes, CNRS, IRD, Grenoble INP, Grenoble, France.

*e-mail: romain.hugonnet@gmail.com

Glaciers distinct from the Greenland and Antarctic ice sheets are shrinking rapidly,
altering regional hydrology1, raising global sea-level2 and elevating natural hazards3.
Yet, due to the scarcity of constrained mass loss observations, glacier evolution during
the satellite era is only known as a geographic and temporal patchwork4,5. Here we
reveal the accelerated, albeit contrasted, patterns of glacier mass loss during the early
twenty-first century. By leveraging 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 the first globally complete and consistent estimate of glacier mass change.
We show that, during 2000-2019, glaciers lost 267 ± 16 Gt yr-1, equivalent to 21 ± 3% of
observed sea-level rise6. We identify a mass loss acceleration of 48 ± 16 Gt yr-1 per
decade, explaining 6-19% of the observed acceleration of sea-level rise. Particularly,
thinning rates of glaciers outside ice sheet peripheries doubled over the last two
decades. Glaciers presently lose more mass, and at similar or larger accelerated rates,
than the Greenland or Antarctic ice sheets taken separately7–9. Uncovering the patterns
of mass change in many regions, we find contrasted glacier fluctuations that agree with
decadal variability in precipitation and temperature. Those include a newly-identified
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 foster 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 management of
local water resources and cryospheric risks as well as for regional-to-global sea-level
rise.

About 200 million people live on land predicted to fall below the high-tide lines of rising sea
levels by the end of the century11, while more than one billion could face water shortage and
food insecurity within the next three decades4. Glaciers distinct from the ice sheets play a
prominent role in these repercussions 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 socio-political instability15, to
predict the alteration of coastal areas due to sea-level rise4, and to assess the impacts on
ecosystems16 as well as on cryosphere-related hazards3.
Despite this, glacier mass change stands out as one of the least-constrained elements of the
global water cycle, identified as a critical research gap by 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 out 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
to (i) the coarse spatial resolution of satellite gravimetry, 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 either estimation, reference, or validation of Earth’s surface elevation
over all glaciers and their vicinity between January 1, 2000 and December 31, 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 the land area on Earth (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 produce 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,
non-linear 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 due to latent, uncorrected ASTER instrument noise, and we propagate our
elevation errors into volume change errors accordingly. We validate the reliability of our
uncertainty estimates down to the scale of individual glaciers by comparison to 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, subsequently converted to water-equivalent mass change23. Our analysis includes
200,000 km² of glaciers located in the Greenland Periphery and in the 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 due to larger uncertainties at shorter time scales (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.

Global contribution to sea-level rise

From 2000 to 2019, global glacier mass loss totalled 266 ± 16 Gt yr-1, a mass loss 47% larger
than that of the GIS, and more than twice that of the AIS7–9 (Table 1). Assuming 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 due 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 were yet unable to discern significant
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, the latter 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, 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.
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 (e.g. 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-contrasted 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 (i.e. mass change divided by area) strongly varied in latitudinal belts (Fig. 2).
The large, northernmost Arctic regions composed of northern Arctic Canada and Greenland
Periphery, Svalbard and Jan Mayen, and the Russian Arctic, all showed moderate specific
mass change rates, averaging -0.28 ± 0.04 m w.e. yr-1. 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 -0.74 ± 0.10 m w.e. yr-1, reaching the world’s most negative regional rate over the
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.
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) where
gravimetric uncertainties are considered small due to weak competing signals27 (Fig. 3). We
note, however, 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 or 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 and -0.67 ± 0.15 m w.e. yr-1, respectively, are slightly
more negative than the -0.19 ± 0.06 and -0.64 ± 0.04 m w.e. yr-1 reported from recent
DEM-based studies28,29. In the 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 of -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.

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 (Fig. 2, Fig. 4a). The peculiar surface elevation change patterns we capture for
glaciers fringing Greenland, particularly notable around the eastern Greenland subregions 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 lost mass (-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.
Only two regions of the world beyond the ice sheet peripheries experienced slowdown of
glacier thinning. Icelandic glaciers’ record thinning rates 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 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 US (9%). Glaciers in the latter region
experienced a fourfold increase in thinning rates. Most notably, glaciers in Northwestern
America (Alaska, Western Canada and US) are responsible for nearly 50% of the accelerated
mass loss. The widespread and strong increase of thinning of High Mountain Asian glaciers
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 1.52 ±
0.50 m yr-1 thinning rate in 2015-2019 , which is a nearly sevenfold increase compared to
2000-2004.
Analysis of climate data reveals that much 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, 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 that instead benefited neighbouring
regions of central Alaska and continental US, 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.
While 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 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 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,
globally, reveals a significant accelerated mass loss. Decadal rates of glacier mass change
remain, however, strongly modulated by regional climatic conditions. We capture the
magnitude of such fluctuations, most contrasted 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 subregion of eastern Greenland Periphery close to balance,
while the strong increase in thinning rates of High Mountain Asian glaciers likely marks the
end of the Karakoram anomaly.
From the spatiotemporally resolved nature of our assessment, novel possibilities arise to
harness observations of the satellite era. Not only instrumental for glaciers, such resolved
estimates 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.

1. Pritchard, H. D. Asia’s shrinking glaciers protect large populations from drought stress. Nature

569, 649–654 (2019).

2. WCRP Global Sea Level Budget Group. Global sea-level budget 1993–present. Earth Syst. Sci.

Data 10, 1551–1590 (2018).

3. Stoffel, M. & Huggel, C. Effects of climate change on mass movements in mountain

environments. Progress in Physical Geography: Earth and Environment 36, 421–439 (2012).

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

IPCC Intergovernmental Panel on Climate Change: Geneva, Switzerland (2019).

5. Gardner, A. et al. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009.

Science 340, 852–857 (2013).

6. Nerem, R. S. et al. Climate-change-driven accelerated sea-level rise detected in the altimeter era.

Proc. Natl. Acad. Sci. U. S. A. 115, 2022–2025 (2018).

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

(2020).

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

(2018).

9. Smith, B. et al. Pervasive ice sheet mass loss reflects competing ocean and atmosphere

processes. Science 368, 1239–1242 (2020).

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).

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).

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

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).

14. Huss, M. & Hock, R. Global-scale hydrological response to future glacier mass loss. Nat. Clim.

Chang. (2018) doi:10.1038/s41558-017-0049-x.

15. IPCC. Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part B: Regional Aspects.

(Cambridge University Press, 2014).

16. Cauvy-Fraunié, S. & Dangles, O. A global synthesis of biodiversity responses to glacier retreat.

Nat Ecol Evol 3, 1675–1685 (2019).

17. World Glacier Monitoring Service (WGMS). Fluctuations of Glaciers Database. (2019)

doi:10.5904/WGMS-FOG-2019-12.

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, (2018).

19. Wouters, B., Gardner, A. S. & Moholdt, G. Global Glacier Mass Loss During the GRACE

Satellite Mission (2002-2016). Front. Earth Sci. 7, 535 (2019).

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).

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).

22. RGI Consortium. Randolph Glacier Inventory – A Dataset of Global Glacier Outlines: Version

6.0: Technical Report, Global Land Ice Measurements from Space, Colorado, USA. (2017)

doi:10.7265/N5-RGI-60.

23. Huss, M. Density assumptions for converting geodetic glacier volume change to mass change.

The Cryosphere 7, 877–887 (2013).

24. Ablain, M. et al. Uncertainty in satellite estimates of global mean sea-level changes, trend and

acceleration. Earth System Science Data 11, 1189–1202 (2019).

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).

26. Larsen, C. F. et al. Surface melt dominates Alaska glacier mass balance. Geophys. Res. Lett. 42,

5902–5908 (2015).

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).

28. Shean, D. E. et al. A Systematic, Regional Assessment of High Mountain Asia Glacier Mass

Balance. Front. Earth Sci. 7, 435 (2020).

29. Braun, M. H. et al. Constraining glacier elevation and mass changes in South America. Nat.

Clim. Chang. (2019) doi:10.1038/s41558-018-0375-7.

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 Journal of Selected Topics in Applied Earth

Observations and Remote Sensing 9, 3870–3882 (2016).

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).

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. U. S. A. 116, 1934–1939 (2019).

33. Garreaud, R. D. et al. The Central Chile Mega Drought (2010–2018): A climate dynamics

perspective. Int. J. Climatol. 40, 421–439 (2020).

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).

35. Parkes, D. & Marzeion, B. Twentieth-century contribution to sea-level rise from uncharted

glaciers. Nature 563, 551–554 (2018).

36. Becker, J. J. et al. Global Bathymetry and Elevation Data at 30 Arc Seconds Resolution:

SRTM30_PLUS. Mar. Geod. 32, 355–371 (2009).

 Mass change rate (Gt yr-1)
Reference
2000-2004 2005-2009 2010-2014 2015-2018* 2000-2018* 2003-2018*

Glaciers This study -227 ± 25 -257 ± 22 -284 ± 23 -292 ± 24 -264 ± 16 -272 ± 16

7
IMBIE minus This study
Greenland (Greenland Periphery) -94 ± 65 -206 ± 56 -267 ± 57 -152 ± 64 -181 ± 31 -205 ± 32
Ice Sheet
Smith et al.9 -200 ± 24
7
IMBIE minus This study
Antarctic (Antarctic and Subantarctic) -36 ± 118 -93 ± 104 -214 ± 94 -157 ± 87 -121 ± 104 -143 ± 104
Ice Sheet
Smith et al.9 -118 ± 48

Table 1. Separating mass losses of glaciers and ice sheets.

Mass losses from glaciers, the GIS and AIS with 95% confidence intervals. Half of the peripheral
glacier component estimated in this study is removed from the ensemble estimates of the ice sheet
mass balance inter-comparison exercise (IMBIE)7,8, see Methods. *The end date is 06/2017 for the
AIS IMBIE estimate.

Fig. 1. Regional glacier mass changes and their temporal evolution from 2000 to 2019.
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 polygons), shown on top of a
world hillshade36. Regions 2, 5, 9, 17 are further divided to illustrate contrasted 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 glaciers (light
blue). Mass change rates larger than 4 Gt yr-1 are printed in blue inside the disk. The outside ring
discerns between land (grey) and marine-terminating (blue) glacier area. Annual time series of mean
elevation change and regional data coverage are displayed on time friezes at the bottom of the disks.

Fig. 2. Spatial distribution of glacier elevation change between 2000 and 2019.
Glacier mean elevation change rate is aggregated for tiles of 1°x1° below 60° latitude, 2°x1° between
60° and 74° latitude and 2°x2° above 74° latitude, thus representing similar surface areas of
approximately 10,000 km². Disks scale with the glacierized area of each tile and are colored according
to the mean elevation change rate (colored in grey if less than 50% 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 below or equal to 10 km² are
displayed at the same minimum size.

Fig. 3. Comparison to earlier global estimates.

Regional and global specific mass change rates with 95% confidence intervals (s.e.m). a, Estimates
from earlier studies5,19–21 (colored) and the estimates from this study over corresponding time periods
(black). b, Global estimates of earlier studies with annual and 5-year rates of this study. Annual rates
of earlier studies are not shown due to large uncertainties. Uncertainties from this study are
conservative, in particular for regions with most negative specific rates (see Methods).
Fig. 4. Decadal patterns of glacier thinning are consistent with decadal variations in
precipitation and temperature.
Difference between 2010-2019 and 2000-2009 for the mean elevation change rates from this study (a)
and the mean annual precipitation (b), mean annual temperature (c) from ERA5. Region inset
numbering corresponds to Fig. 2. Tiles are 1°x1° with size inversely scaled to uncertainties in the
mean elevation change difference shown on top of a world hillshade36. Minimum tile area is 10% for
95% CI larger than 1 m yr-1 and tiles are displayed at full size for 95% CI smaller than 0.4 m yr-1.

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

Glacier inventories
We used the Randolph Glacier Inventory 6.022 (RGI) outlines for all regions except for
Caucasus Middle East (region 12). Due 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 to 1,336 km2. In
Svalbard and Jan Mayen (region 7), we manually updated glacier outlines to account for
advances resulting from major surges38–40, increasing mapped areas by 228 km² (Extended
Data Fig. 6). In the Greenland Periphery (region 5), we did not analyze the 955 glaciers
strongly connected to the ice sheet (RGI connectivity level 2) with an area of 40,354 km², as
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 km². For the purpose of
co-registering and bias correcting DEMs, we masked ice-covered terrain using the RGI 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 Antarctica45
(REMA) data intersecting glaciers worldwide (Extended Data Fig. 2), totalling more than 100
TB of data storage. 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 degrees, absolute elevation difference <200 m) were available
at each iterative step.
From 440,548 ASTER L1A stereo-images43 (each covering 60 km by 60 km), we generated,
co-registered and bias-corrected 154,656 ASTER DEM strips (30 m resolution; 180 km by 60
km strip size) using improved techniques of MicMac for ASTER49,50. Novel 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 x 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 the
Supplementary Information.
From 97,649 release 7 ArcticDEM44 DEMs at 2 m resolution and 13,790 release 1.1 REMA45
DEMs at 2 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 𝑡, space 𝑥 and 𝑦) of elevation
ℎ(𝑡, 𝑥, 𝑦), divided into 2,106 1 x 1° tiles containing glaciers and downsampled to 100 m to
decrease computing time.
In order 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 is 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 maximums were first conditioned by extreme
values (e.g., maximum observed absolute glacier elevation change of 50 m yr-1 on HPS12
glacier, Southern Patagonian Icefield53). Later, those were refined by estimating a linear
glacier elevation change in the surroundings through weighted least-squares54.
In our first global-scale statistical modelling step, we identified an heteroscedasticity in
elevation measurements (i.e., non-uniform variance, Extended Data Fig. 3a). We determined
that the elevation measurement error σℎ varied with the terrain slope55,56 α, the quality of
stereo-correlation49,57 𝑞 and the individual performance of each DEM’s co-registration48
σ𝑐(𝑡, 𝑥, 𝑦). 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 (e.g., slope) in each 1x1° 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
2
effects of elevation outliers58. We found that the empirical variances for the slope σα and the
quality of stereo-correlation σ𝑞² 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 ℎ(𝑡, 𝑥, 𝑦):

2 2 2 2
σℎ(𝑡, 𝑥, 𝑦) = σ𝑐(𝑡, 𝑥, 𝑦) + σα(α, 𝑞) + σ𝑞(𝑞) (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 time lag between pairwise elevation
observations ∆𝑡 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 non-linear trend (the RQ times PL).
Empirical covariances showed little variability between regions. We thus conditioned the
parameters of the ESS, RBF and RQ kernels (ϕ𝑝,σ𝑝,∆𝑡𝑙,σ𝑙,∆𝑡𝑛𝑙,σ𝑛𝑙 and α𝑛𝑙) at the global scale
based on our empirical variograms, while the PL kernel was determined directly from the
observations of each pixel (𝑥, 𝑦):

2
( ) (
σℎ(𝑥, 𝑦, ∆𝑡) = 𝑃𝐿(𝑥, 𝑦, ∆𝑡) + 𝐸𝑆𝑆 ϕ𝑝, σ𝑝², ∆𝑡 + 𝑅𝐵𝐹 ∆𝑡𝑙, σ𝑙², ∆𝑡 + )
2
( 2
)
𝑅𝑄 ∆𝑡𝑛𝑙, σ𝑛𝑙 , α𝑛𝑙, ∆𝑡 · 𝑃𝐿(𝑥, 𝑦, ∆𝑡) + σ (𝑡, 𝑥, 𝑦)
ℎ
(2)

Applying GP regression, we iteratively removed observations outside the 20-, 12-, 9-, 6- and
4-sigma 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 millions
pixels (𝑥, 𝑦) 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
the Supplementary Information and build on additional references60–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 x 500 m that
we downsampled 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 center 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 impact our annual estimates once differenced into elevation
change, 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 Gaussian Process regression) are conservative (i.e. too large by a factor of
about two). We reached the same conclusions at the scale of individual RGI 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 elevation time series are available
in the Supplementary Information and build on an additional reference69.

Integration of elevation into volume changes

We differenced all elevations ℎ into elevation change based on their value of h on January
1st, 2000. We integrated elevation change 𝑑ℎ into volume change 𝑑𝑉 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,
hence assuring 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, elevation
change was linearly interpolated from adjacent bins, or extrapolated from the closest bins.
For retreating lake- and ocean- terminating glaciers, we exclude 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

uncertainties by assuming that the uncertainty in the mean elevation change σ𝑑ℎ and the
uncertainty in the glacier area σ𝐴 are independent:
 2
2
(
σ𝑑𝑉 = σ𝑑ℎ𝐴 )2 + (σ𝐴𝑑ℎ) (3)

The uncertainty in the mean elevation change σ𝑑ℎ 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 neighboring pixels of a given region sharing similar
temporal data gaps, and are hence likely to have similar interpolation biases which
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 ∆𝑡 to the closest ASTER, ArcticDEM or REMA observation.
For each time lag, we estimated the partial sill 𝑠𝑘 (correlated variance) by fitting a sum of
( )
seven spherical variogram models 𝑆 𝑑, 𝑠𝑘, 𝑟𝑘 , with 𝑑 the spatial lag, at ranges 𝑟𝑘 (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 𝑠1 to 𝑠6. We then used the GP
elevation change uncertainties of each glacier pixel to derive 𝑠0. Finally, we propagated the
pixel-wise uncertainties in elevation change into the uncertainty in the mean elevation change
σ𝑑ℎ by circular integration of the sum of variograms72 over the glacier area 𝐴 (Extended Data
Fig. 5c):

6
1
σ𝑑ℎ² = 𝐴
∑ σ𝑑ℎ ² (4)
𝑘=0 𝑘

where σ𝑑ℎ ² is the integrated variance component correlated with range 𝑟𝑘:
𝑘

( (
σ𝑑ℎ ² = ∫ 𝑠𝑘 − 𝑆 𝑑, 𝑠𝑘, 𝑟𝑘 𝑑𝐴
𝑘 𝐴
)) (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–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 the Supplementary
Information and build on additional references78–83.
For each glacier, we estimated an uncertainty in the area σ𝐴 based on a buffer84 of 15 m
corresponding to the typical resolution of the optical imagery used to derive these
outlines37,85–87. These uncertainties vary from about 0.1% of the area for large icefields (>1000
km²) to 50% of the area and above for small isolated glaciers (<0.1 km²).

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 km²
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 95% uncertainties are
lower than 0.3 m yr-1 for glaciers larger than 1 km² 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
between glaciers of the same region, we extend 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 σ𝑑ℎ . Having been integrated
𝑅

once over a spatial support (from pixel to glacier), the glacier-wide errors 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:

6
σ𝑑ℎ ² =
𝑅
1
𝐴𝑅²
𝑖 𝑗 𝑘=0
( 𝑘,𝑖 𝑘,𝑗
(
· ∑ ∑ ∑ σ𝑑ℎ σ𝑑ℎ − 𝑆 𝐺𝑖 − 𝐺𝑗, σ𝑑ℎ σ𝑑ℎ , 𝑟𝑘 𝐴𝑖𝐴𝑗
𝑘,𝑖 𝑘,𝑗
)) (6)

where 𝑖, 𝑗 are indexes for glaciers in the region, σ𝑑ℎ is the uncertainty in the mean elevation
𝑘,𝑖

change σ𝑑ℎ with range 𝑟𝑘 and sill 𝑠𝑘 for glacier i, 𝐺𝑖 − 𝐺𝑗 is the pairwise distance (spatial lag
𝑘

𝑑) between glaciers 𝑖 and 𝑗 based on their outline centroids and 𝐴𝑖 is the area of glacier 𝑖.

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 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 conversative
assumption due 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 regions, we summed regional uncertainties quadratically. We report uncertainties in mass
changes for periods shorter than five years solely for the global or near-global estimates (e.g.
Fig. 3b) by assuming that the aggregation of largely independent RGI 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–96 (Supplementary Table S4).

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 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 (i.e., 2000-2004, 2005-2009,
2010-2014 and 2015-2019), propagating their related uncertainties as independent. While
shorter time scales and smaller spatial domains are affected by temporal autocorrelation, we
assumed the 5-year estimates at the global or near-global scale (i.e., excluding peripheral
glaciers) as temporally uncorrelated. This assumption is supported by time scales 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 regions 5 and 19 to
ice sheet estimates from IMBIE7,8. Because IMBIE estimates are a weighted mean of three
ensemble estimates where half include peripheral glaciers, the other half does not
(gravimetric studies include peripheral glaciers, altimetric studies exclude peripheral glaciers,
and input-output studies can do both), we assumed that subtracting half of our estimates for
the peripheral glaciers was most adequate. Noteworthily, 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 analyzed ERA5 precipitation and temperature97 at both annual and seasonal scales. For
the latter we consider 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 degree 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
detect 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.

37. Tielidze, L. G. & Wheate, R. D. The Greater Caucasus Glacier Inventory (Russia, Georgia and

Azerbaijan). The Cryosphere 12, 81–94 (2018).

38. Dunse, T. et al. Glacier-surge mechanisms promoted by a hydro-thermodynamic feedback to

summer melt. The Cryosphere 9, 197–215 (2015).

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).

40. Nuth, C. et al. Dynamic vulnerability revealed in the collapse of an Arctic tidewater glacier. Sci.

Rep. 9, 5541 (2019).

41. Howat, I. M., Negrete, A. & Smith, B. E. The Greenland Ice Mapping Project (GIMP) land

classification and surface elevation data sets. The Cryosphere vol. 8 1509–1518 (2014).

42. Fretwell, P. et al. Bedmap2: improved ice bed, surface and thickness datasets for Antarctica. The

Cryosphere 7, 375–393 (2013).

43. NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team. ASTER Level

1A Data Set - Reconstructed, unprocessed instrument data. (2001)

doi:10.5067/ASTER/AST_L1A.003.

44. Porter, C. et al. ArcticDEM. (2018) doi:10.7910/DVN/OHHUKH.

45. Howat, I. M., Porter, C., Smith, B. E., Noh, M.-J. & Morin, P. The Reference Elevation Model of

Antarctica. The Cryosphere 13, 665–674 (2019).

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).

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).

48. Nuth, C. & Kääb. Co-registration and bias corrections of satellite elevation data sets for

quantifying glacier thickness change. The Cryosphere 5, 271–290 (2011).

49. Rupnik, E., Daakir, M. & Pierrot Deseilligny, M. MicMac – a free, open-source solution for

photogrammetry. Open Geospatial Data, Software and Standards 2, 14 (2017).

50. Girod, L., Nuth, C., Kääb, A., McNabb, R. & Galland, O. MMASTER: Improved ASTER DEMs

for Elevation Change Monitoring. Remote Sensing 9, 704 (2017).

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).

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).

53. Dussaillant, I. et al. Two decades of glacier mass loss along the Andes. Nat. Geosci. 12, 802–808

(2019).

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).

55. Toutin, T. Three-dimensional topographic mapping with ASTER stereo data in rugged

topography. IEEE Trans. Geosci. Remote Sens. 40, 2241–2247 (2002).

56. Lacroix, P. Landslides triggered by the Gorkha earthquake in the Langtang valley, volumes and

initiation processes the 2015 Gorkha, Nepal, Earthquake and Himalayan Studies: First Results 4.

Seismology. Earth Planets Space 68, 1–10 (2016).

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).

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).

59. Williams, C. K. I. & Rasmussen, C. E. Gaussian processes for machine learning. vol. 2 (MIT

press Cambridge, MA, 2006).

60. Schiefer, E., Menounos, B. & Wheate, R. Recent volume loss of British Columbian glaciers,

Canada. Geophys. Res. Lett. (2007) doi:10.1029/2007GL030780.

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).

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, (2012).

63. Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. J. Mach. Learn. Res. 12,

2825–2830 (2011).

64. Zwally, H. J., Schutz, R., Hancock, D. & Dimarzio, J. GLAS/ICESat L2 Global Land Surface

Altimetry Data (HDF5), Version 34. (2014) doi:10.5067/ICESAT/GLAS/DATA211.

65. Alexandrov, O., McMichael, S. & Beyer., R. A. IceBridge DMS L3 Ames Stereo Pipeline

Photogrammetric DEM, Version 1. (2018).

66. Larsen, C. IceBridge UAF Lidar Scanner L1B Geolocated Surface Elevation Triplets, Version 1.

(2010, updated 2020).

67. Beyer, R. A., Alexandrov, O. & McMichael, S. The Ames Stereo Pipeline: NASA’s Open Source

Software for Deriving and Processing Terrain Data. Earth and Space Science 5, 537–548 (2018).

68. Harding, D. J. ICESat waveform measurements of within-footprint topographic relief and

vegetation vertical structure. Geophys. Res. Lett. 32, 2509 (2005).

69. Gardelle, J., Berthier, E. & Arnaud, Y. Impact of resolution and radar penetration on glacier

elevation changes computed from DEM differencing. J. Glaciol. 77, 121–128 (2012).

70. Mcnabb, R., Nuth, C., Kääb, A. & Girod, L. Sensitivity of glacier volume change estimation to

DEM void interpolation. The Cryosphere c, 895–910 (2019).

71. Cressie, N. A. C. Statistics for spatial data. vol. 4 613–617 (Wiley, 1993).

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).

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, 516 (2020).

74. Menounos, B. et al. Heterogeneous changes in western North American glaciers linked to

decadal variability in zonal wind strength. Geophys. Res. Lett. 2018GL080942 (2018).

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, 1–5

(2008).

76. Wang, D. & Kääb, A. Modeling glacier elevation change from DEM time series. Remote Sensing

7, 10117–10142 (2015).

77. Cogley, J. G. & Adams, W. P. Mass balance of glaciers other than the ice sheets. J. Glaciol. 44,

315–325 (1998).

78. Journel, A. G. & Huijbregts, C. J. Mining geostatistics. vol. 600 (Academic press London, 1978).

79. Webster, R. & Oliver, M. A. Geostatistics for Environmental Scientists. (2007).

80. Gräler, B., Pebesma, E. & Heuvelink, G. Spatio-Temporal Interpolation using gstat. The R

Journal vol. 8 204 (2016).

81. Mälicke, M. & Schneider, H. D. Scikit-GStat 0.2.6: A scipy flavored geostatistical analysis
 toolbox written in Python. (2019). doi:10.5281/zenodo.3531816.

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).

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, 1–6 (2012).

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).

85. Tad Pfeffer, W. et al. The Randolph Glacier Inventory: a globally complete inventory of glaciers.

J. Glaciol. 60, 537–552 (2014).

86. Rastner, P. et al. The first complete inventory of the local glaciers and ice caps on Greenland. The

Cryosphere 6, 1483–1495 (2012).

87. Bolch, T., Menounos, B. & Wheate, R. Landsat-based inventory of glaciers in western Canada,

1985-2005. Remote Sens. Environ. 114, 127–137 (2010).

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. The Cryosphere 13,

1709–1727 (2019).

89. Wagnon, P. et al. Seasonal and annual mass balances of Mera and Pokalde glaciers (Nepal

Himalaya) since 2007. The Cryosphere 7, 1769–1786 (2013).

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).

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).

92. Glacier Monitoring Switzerland. GLAMOS Swiss Glacier Volume Change, release 2018. (2018)

doi:10.18750/volumechange.2018.r2018.

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).

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).

95. Schuler, T. V. et al. Reconciling Svalbard Glacier Mass Balance. Front Earth Sci. 8, 156 (2020).

96. Aðalgeirsdóttir, G. et al. Glacier Changes in Iceland From ∼1890 to 2019. Front Earth Sci. 8, 520

(2020).

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

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).

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).

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

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 of this article are
publicly available at: https://doi.org/10.5281/zenodo.4530314.

Code availability
The code developed for the global processing and analysis of all data, and to generate figures
and tables of 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 (v0.1) at https://pypi.org/project/pymmaster (with supporting documentation at
https://mmaster-workflows.readthedocs.io) and for processing DEM time series as the python
package pyddem (v0.1) at https://pypi.org/project/pyddem (with supporting documentation at
https://pyddem.readthedocs.io).

Acknowledgments
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 University of Northern
British Columbia provided computational resources for processing ASTER stereo imagery.
SPOT6/7 data were obtained thanks to GEOSUD (ANR-10-EQPX-20, program
‘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 provided
through a Blue Waters Innovation Initiative. DEMs produced using data from DigitalGlobe,
Inc. RH acknowledges a fellowship from the University of Toulouse. EB acknowledges
support from the French Space Agency (CNES) through ISIS and TOSCA programs. RM,
CN, LG and AK 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 no. 320816. BM acknowledges funding from the
National Sciences and Engineering Research Council of Canada, the Canada Research Chairs
Program and Global Water Futures. Authors with ETHZ affiliation acknowledge funding
from the Swiss National Science Foundation, Grant Nr. 184634.

Author contribution
EB and RH designed the study with contributions from DF, MH and BM. LG, CN, RM and
AK developed ASTER bias-correction methods. RH and RM developed glacier elevation
Gaussian Process methods. RH implemented spatial statistics methods with inputs from FB.
BM assembled and analyzed ERA5 data. RH performed the processing and analysis of all
data with main inputs from EB, as well as RM, BM, DF, MH, ID and FB. All authors
interpreted the results. RH led the writing of the paper and all other co-authors contributed.

Competing interests
The authors declare no competing interests.

Additional information
Supplementary information is available for this paper.
Correspondence and requests for materials should be addressed to RH.
Reprints and permissions information is available at www.nature.com/reprints.

Extended Data

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 of Methods.

Extended Data Fig. 2. Spatial and temporal coverage of ASTER, ArcticDEM and
REMA DEMs.
Spatial distribution of DEMs as a strip count for ArcticDEM strips above 50 degrees North (a),
ASTER DEM strips (b), and REMA strips below 50 degrees South (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 respectively 40,391 and 3,456 due to limited stable
terrain in polar regions. d, Temporal distribution of the strip count as a bi-mensual histogram from
January 2000 to December 2019. Note that ArcticDEM and REMA strip footprints (15 km by 50 km)
are generally much smaller than ASTER DEM strip footprints (180 km by 60 km).

Extended Data Fig. 3. Elevation time series estimation.

Empirical and modelled elevation measurement error (a) and temporal covariance of glacier elevation
(b) estimated globally. Those are used to condition the filtering (c,d) and elevation time series
estimation (e) of elevation observations, illustrated here for a 100 m by 100 m pixel on the ablation
area of Upsala where a strong nonlinear elevation loss occurred99. a, Squared measurement error,
estimated by the squared Normalized Median Absolute Deviation 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 year. 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 per year. The median of the
linear trend at these locations (-1.2 m per year) 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 Gaussian Process regression fits for credible intervals of
size 20, 12, 9, 6 and 4 standard deviations. e, Elevation time series of final Gaussian Process
regression after the removal of outliers.

Extended Data Fig. 4. Validation of elevation time series and uncertainties to ICESat
and IceBridge.
ICESat64 and IceBridge65,66 measurements compared to our surface elevation time series over
glacierized terrain in the Saint-Elias Mountains, Alaska (a-c), and at the global scale (d). b, Absolute
z-scores (white-to-purple) on top of 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 by 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, seasons, 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 and validation to high-resolution DEMs.

Spatial correlation of elevations between the Gaussian Process (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 (d-f) and
contribution from all uncertainty sources to the 2000-2019 specific mass change estimates (g,h). a,
Empirical spatial variogram is shown and fitted with a sum of spherical models at correlation length
of 0.1, 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 used for the variance used during error 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, panel (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 year for 50% of observations,
1-2-years for 20% of observations, 2-3 years for 20% of observations and 3-4 years for 10% of
observations. We assume a mean pixel-wise error 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-sigma uncertainties. The mean of differences in estimates over all glaciers does not statistically
differ from zero. e,f, Theoretical and empirical uncertainties, and their evolution with glacier size. The
theoretical uncertainty is the mean of glacier errors derived from spatially integrated variograms and
the empirical uncertainty is the Normalized Median Absolute Difference of the difference between
high-resolution and GP estimates. g,h, Propagation of uncertainty sources to specific mass changes for
each RGI region, and all glaciers with and without the Greenland Periphery and the Antarctic and
Subantarctic which are zoomed on panel (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.
Elevation change of glaciers between 2000 and 2019 in Coropuna, Peru (a), the Pamir Mountains (b),
Iceland (c), the Karakoram Mountains (d), the European Alps (e), the Southern Alps, New Zealand (f)
, West Greenland (note rotated orientation of map) (g), and Svalbard (h). Except for Svalbard, glacier
outlines displayed are from the RGI6.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.

Mean elevation change rates aggregated by tiles of 1°x1° for the periods 2000-2004 (a), 2005-2009
(b), 2010-2014 (c) and 2015-2019 (d). Tile area is inversely scaled to the squared 95% CI of the mean
elevation change in the tile and tiles are colored with mean elevation change rates, on top of a world
hillshade36. Minimum tile area is 10% for 95% CI larger than 2 m yr-1 and tiles are displayed at full
size for 95% CI smaller than 0.5 m yr-1. Region labeling 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.
Regional and global mean elevation change and mass change rates over full and 5-year periods of
2000-2019. Mean elevation change is the volume change divided by time-evolving regional glacier
areas (see Methods)21. Areas reported are those of the RGI 6.0 inventory22, except for region 12
(Caucasus Middle East) which was updated with more recent outlines37. Periods are inclusive and
refer to calendar years of 1st January to 31st December. Uncertainties correspond to 95% confidence
intervals. In Greenland, glaciers highly connected to the ice sheet (RGI connectivity level 2) are not
reported.

Extended Data Table 2. Regional data coverage of elevation time series from 2000 to
2019.
Spatial and temporal coverage of our elevation time series after the three steps of elevation outlier
filtering. Nominal glaciers correspond to uncharted glaciers inventoried in the RGI with only an
estimated surface area, present notably in region 10 (North Asia) where they contribute to 3.0% of the
region total glacier area. Those are accounted for in our volume change estimates by applying the
mean elevation change of the region to their reported area. Glaciers without any coverage correspond
to glaciers having no valid, post-filtering elevation change observation within their outline. This
generally occurs when repeat spatial sampling is poor (less than 3 observations in 20 years) for small
glaciers located in high slopes.

Extended Data Table 3. Regional rates of land- and marine-terminating glaciers in

maritime regions.
Uncertainties correspond to 95% confidence intervals. For marine-terminating glaciers, subaqueous
losses are not included (see Methods).