# Ice velocity and thickness of the world's glaciers

Romain Millan, Jérémie Mouginot, Antoine Rabatel, Mathieu Morlighem

► To cite this version:

Romain Millan, Jérémie Mouginot, Antoine Rabatel, Mathieu Morlighem. Ice velocity and thickness of the world's glaciers. Nature Geoscience, 2022, 15, pp.124-129. 10.1038/s41561-021-00885-z. insu-03706390

HAL Id: insu-03706390

https://insu.hal.science/insu-03706390v1

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Ice velocity and thickness of the world's glaciers

Romain Millan ID<fnref n="1" />,²✉, Jérémie Mouginot ID<fnref n="1" />,³, Antoine Rabatel ID⁴ and Mathieu Morlighem ID³,⁴

The effect of climate change on water resources and sea-level rise is largely determined by the size of the ice reservoirs around the world and the ice thickness distribution, which remains uncertain. Here, we present a comprehensive high-resolution mapping of ice motion for 98% of the world's total glacier area during the period 2017-2018. We use this mapping of glacier flow to generate an estimate of global ice volume that reconciles ice thickness distribution with glacier dynamics and surface topography. After reallocating volume connected to the Antarctic ice sheet, the results suggest that the world's glaciers have a potential contribution to sea-level rise of 257 ± 85 mm, which is 20% less than previously estimated. At low latitudes, our findings highlight notable changes in freshwater resources, with 34% more ice in the Himalayas and 27% less ice in the tropical Andes of South America, affecting water availability for local populations. This mapping of glacier flow and thickness redefines our understanding of global ice-volume distribution and has implications for the prediction of glacier evolution around the world, since accurate representations of glacier geometry and dynamics are of prime importance to glacier modelling.

The size and shape of glaciers and ice caps are controlled by changes in the atmosphere and the ocean. Their current shrinkage provides direct and unequivocal evidence of global climate change<fnref n="1" />,². It has become critical to model their future evolution as they provide local populations with freshwater resources<fnref n="1" /> and currently contribute 25–30% to sea-level rise³,⁴, which threatens about 10% of the world’s population living less than 10 m above sea level⁵. While international efforts over the past decades have substantially improved the representation of the ice masses⁶⁻¹⁵, the thickness distribution of the world’s glaciers is still subject to large uncertainties.

In situ measurements of ice thickness exist for only 4,700 of the world's 215,000 glaciers, and most of these glaciers have only partial coverage⁶. A few model-based reconstructions of global glacier thickness have been proposed, with estimates that vary strongly within Randolph Glacier Inventory (RGI) regions, depending on the model used⁹,¹⁰ (Supplementary Table 1), highlighting the lack of confidence that remains in these reconstructions. Indeed, current ice thickness models are mostly using simplified glacier geometry based on width-averaged centreline approaches that are applied glacier by glacier (relying on the glacier's basin boundary), with strong assumptions on the basal shear stress or surface mass balance gradients. While these methods might be reasonable for valley glaciers, they lead to unrealistic bedrock geometries for ice caps and tidewater glaciers, where most of the global ice volume outside Greenland and Antarctica⁶,⁷,⁹,¹⁶,¹⁷ is located. Such inconsistencies in glacier thickness need to be addressed as it has been shown that ice thickness distribution has a greater influence than the total ice volume in projections of glacier evolution for the decades to come⁷,¹⁸,¹⁹. In other words, getting the correct global volume of ice does not help to model a glacier's future if its thickness distribution is not properly captured.

Ice flow defines the transfer of ice within a glacier and largely controls the spatial distribution of the ice volume. While being an essential constraint in glacier models, to obtain realistic simulations of ice dynamics, large-scale mapping of glacier flow remains relatively difficult due to the limited amount of satellite observations with sufficient quality and resolution. This lack of data

hampers our ability to capture a large diversity of glaciers, from small Andean glaciers (<1 km large) flowing at less than 50 m yr⁻<fnref n="1" /> up to large marine terminating glaciers (width of several kilometres), with surface velocity of several kilometres per year. Although some individual regions have been mapped²⁰⁻²⁴, there is to date no global and complete view of glacier flow with sufficient spatial resolution. With the launch during the past decade of the optical sensors Landsat 8 (NASA-US Geological Survey), Sentinel 2A and 2B and the synthetic aperture radar sensors Sentinel-1A and 1B (European Space Agency (ESA) Copernicus), these limitations can be overcome, and the worldwide availability of observations at decametric resolution opens new opportunities for mapping the ice flow of the Earth's glaciers²²,²⁵. The aim of this paper is to (1) provide a new high-resolution coverage of the surface flow velocity of the world's glaciers (contained in the RGI²⁶) and (2) use this new view of glacier dynamics to re-estimate the volume of ice and its sea-level potential.

## Mapping ice flow velocity at global scale

In this article, we present a high-resolution mapping of surface velocity of the world's glaciers derived from about 135,000, 670,000 and 6,684 image pairs acquired between 2017 and 2018 by Landsat 8, Sentinel-2 and Sentinel-1, respectively. For each pair of images acquired up to 400 days apart, glacier motion is tracked using a cross-correlation approach following the displacement of surface features or speckle²⁷,²⁸ (Methods). The resulting displacements are assembled into 62 regional maps at a sampling resolution of 50 m, which is 2–5 times more resolved than available products²⁴, at an accuracy of about 10 m yr⁻<fnref n="1" /> (ref. ²⁵).

Our mapping covers more than 98% of the worldwide glacierized area (Table 1), including regions where no mapping existed before. Examples of such regions include glaciers in Caucasus, New Zealand, the southern cordilleras of South America, sub-Antarctic islands and Europe (compare Fig. 1 with Supplementary Fig. 7). This could be achieved only by merging the data of the three sensors and using a wide range of temporal baselines, allowing comprehensive coverage from the polar regions, where radar provides better correlation, and the mid-latitudes, where optical

<footnote n="1">¹Université Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, Grenoble, France. ²Department of Geosciences and Natural Resource Management, University of Copenhagen, Copenhagen K, Denmark. ³Department of Earth System Science, University of California Irvine, Irvine, CA, USA. ⁴Department of Earth Sciences, Dartmouth College, Hanover, NH, USA. ✉e-mail: romain.millan@univ-grenoble-alpes.fr</footnote>
<table><caption>Table 1 | Ice volume and SLE of glaciers outside the ice sheets</caption><thead><tr><th>Region name</th><th>RGI region</th><th>Glacierized area (×10<sup>3</sup> km<sup>2</sup>)</th><th>Coverage (%)</th><th>Ice volume<sup>a</sup> (×10<sup>3</sup> km<sup>3</sup>)</th><th>SLE<sup>a</sup> (mm)</th><th>Regional difference (%)</th></tr></thead><tbody><tr><td>Alaska and Western Canada</td><td>1, 2</td><td>101</td><td>96</td><td>19.2 ± 5.6</td><td>46.4 ± 15.3</td><td>-4</td></tr><tr><td>Arctic Canada North</td><td>3</td><td>105</td><td>99</td><td>25.4 ± 7.2</td><td>59.9 ± 19.2</td><td>-10</td></tr><tr><td>Arctic Canada South</td><td>4</td><td>41</td><td>100</td><td>7.0 ± 2.1</td><td>17.7 ± 5.8</td><td>-19</td></tr><tr><td>Greenland Periphery</td><td>5</td><td>90</td><td>100</td><td>11.8 ± 3.7</td><td>26.9 ± 9.5</td><td>-25</td></tr><tr><td>Iceland</td><td>6</td><td>11</td><td>99</td><td>3.7 ± 0.9</td><td>9.4 ± 2.6</td><td>-1</td></tr><tr><td>Svalbard</td><td>7</td><td>34</td><td>97</td><td>7.0 ± 2.3</td><td>15.4 ± 5.7</td><td>-6</td></tr><tr><td>Scandinavia</td><td>8</td><td>3</td><td>84</td><td>0.29 ± 0.1</td><td>0.73 ± 0.3</td><td>-4</td></tr><tr><td>Russian Arctic</td><td>9</td><td>52</td><td>98</td><td>15.5 ± 3.9</td><td>33.7 ± 9.6</td><td>6</td></tr><tr><td>North Asia</td><td>10</td><td>2</td><td>63</td><td>0.11 ± 0.05</td><td>0.28 ± 0.1</td><td>-22</td></tr><tr><td>Central Europe</td><td>11</td><td>2</td><td>99</td><td>0.12 ± 0.05</td><td>0.29 ± 0.1</td><td>-11</td></tr><tr><td>Caucasus and Middle East</td><td>12</td><td>1</td><td>78</td><td>0.05 ± 0.03</td><td>0.14 ± 0.1</td><td>-8</td></tr><tr><td>Asia</td><td>13, 14, 15</td><td>98</td><td>100</td><td>9.4 ± 3.7</td><td>23.8 ± 10.5</td><td>34</td></tr><tr><td>Low Latitudes</td><td>16</td><td>2</td><td>82</td><td>0.07 ± 0.04</td><td>0.18 ± 0.1</td><td>-27</td></tr><tr><td>Southern Andes</td><td>17</td><td>29</td><td>99</td><td>5.9 ± 1.6</td><td>14.6 ± 4.4</td><td>10</td></tr><tr><td>New Zealand</td><td>18</td><td>1</td><td>99</td><td>0.07 ± 0.03</td><td>0.18 ± 0.1</td><td>2</td></tr><tr><td>Antarctic and Islands</td><td>19</td><td>133</td><td>100</td><td>35.1 ± 9.1</td><td>61.3 ± 16.0</td><td>-25</td></tr><tr><td>Total</td><td></td><td>706</td><td>98</td><td>140.6 ± 40.4</td><td>256.8 ± 85<sup>b</sup></td><td></td></tr></tbody></table>

RGI region: number in the Randolph Glacier Inventory region. Area: glacierized area in km². Coverage: percentage of total coverage in the ice velocity mapping. Ice vol.: total ice volume in km³ within the RGI. SLE: sea level equivalent within the RGI in mm (see Methods) calculated by integrating only the ice volume above floatation. Reg. Diff.: Regional difference in percent of the ice volume within the RGI of our study compared to the 2019 consensus⁹. *Error budget accounts for the presence of surges (see Methods). †Total SLE given after removing glacierized area connected to the Antarctic ice sheet (54.1 mm SLE)³⁷.

sensors provide the resolution needed to resolve small glaciers. This multi-temporal and multi-sensor approach is crucial to capture the wide variety of glaciers in terms of size and speed, from narrow alpine-type glaciers to large ice caps and from slow valley glaciers to fast tidewater glaciers. Alpine-type glaciers are found throughout all studied mountain ranges, with surface velocities rarely exceeding 500 m yr⁻¹. There are, however, a few exceptions of unusually fast flow velocities for mountain glaciers such as the Gotley Glacier flowing from Mawson Peak on Heard Island in the southern Indian Ocean (RGI-19) with velocities up to 1,500 m yr⁻¹ or Bualtar Glacier in Karakoram (RGI-14), which reaches 2,300 m yr⁻¹ over the study period (Fig. 1). Our mapping includes large ice caps in the Arctic (for example, Severnaya Zemlya, Franz Josef and Devon ice cap), with comprehensive high-resolution mapping that extends up to the ice divide, where velocities are less than a few metres per year, while marine outlet glaciers can reach hundreds to thousands of metres year along their margins (Fig. 1 and Supplementary Fig. 7). The fastest recorded outlet glacier is a surging one, flowing out of the Austfonna ice cap (basin 3 (ref. ²⁶)), with an ice velocity exceeding 2,900 m yr⁻¹. Finally, marine terminating glaciers flowing out from large ice fields in Alaska, the Canadian Arctic and Patagonia are among the fastest glaciers on Earth. In southern Patagonia, Penguin Glacier flows at a speed of about 12,000 m yr⁻¹ at its terminus (Fig. 1), which is comparable to the fastest glacier in Greenland: Jakobshavn Isbræ.

**Revealing a different picture of glacier geometry.** This new comprehensive dataset offers an opportunity to constrain the calculation of ice volume contained in the world's glaciers and ice caps on the basis of surface flow velocities. The surface motion and slopes of a glacier can indeed be related to the ice thickness, assuming that the ice deforms in simple shear²⁹. In this relation, we calibrate parameters related to ice rheology and basal sliding (Methods) using in situ ice thickness measurements from the Glacier Thickness Database⁶ (GlaThiDa), the French GLACIOCLIM Observatory³⁰ and other

independent campaigns³¹⁻³⁶. We calculate ice thickness by modelling the shape of glaciers using a three-dimensional approach at regional scales (Methods). Finally, we generate ice thickness maps globally, with corresponding error maps, at the same resolution as the ice velocity maps, with an overall uncertainty of 30% (Supplementary Figs. 2 and 3) and an overall agreement with observed thicknesses of 9 ± 87 m (Supplementary Table 2).

Using the new constraint provided by the surface velocity, the global volume of ice within the RGI is estimated at (140.6 ± 40.4) × 10³ km³ (Table 1). The largest ice reservoirs are located in the Arctic region (including the Greenland periphery, Canadian and Russian Arctic, Iceland and Svalbard), North America and Antarctic periphery with (70.4 ± 20.1) × 10³ km³, (19.2 ± 5.6) × 10³ km³ and (35.1 ± 9.1) × 10³ km³, respectively. By comparison, the ice volume in the high mountains of Asia is (9.4 ± 3.7) × 10³ km³, or 7% of the global estimate. Overall, the total volume above floatation contained in the RGI corresponds to a potential contribution to sea level of 311 ± 99 mm (Table 1). However, we notice that glaciers connected to the Antarctic ice sheet²⁶,³⁷, such as Alexander Island in the Peninsula (Supplementary Fig. 13A), were accounted for in glacier volume estimates, while other similar islands with the same degree of connectivity were not (for example, Berkner and Roosevelt islands in East Antarctica). Thus, when these areas are excluded from the RGI (Supplementary Fig. 13B), as it was done for Greenland in previous studies⁹, the corrected potential contribution of the world's glaciers is 257 ± 85 mm, or 20% lower than previous estimates (Table 1, ref. ⁹).

On the regional level, our results show substantial differences in comparison with previous estimates⁷⁻¹¹ (Supplementary Table 1 and Supplementary Figs. 8–12) (for example, reaching more than ±35% at the regional level, Table 1), which may strongly impact the temporal response of glaciers to climate forcing. These differences are partly explained by the fact that the 2019 consensus⁹ accounts for volume of ice above sea level and not above floatation but, more
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | **Ice velocity of the world's glaciers.** Centre world map illustrates the coverage by RGI region²⁶. The 44 panels are 25×25 km portions of ice velocity, representative of the glacier diversity across all RGI regions. The ice speed is colour coded on a logarithmic scale overlaid on shaded relief. RGI regions are assigned a colour that is also used for the frames of the ice velocity panels.</figcaption></figure>

importantly, by the use here of a regional inversion methodology based on a new dataset of surface flow velocity, which is independent from earlier mass fluxes estimations that were subject to strong assumptions⁷. Volumes 15–25% smaller are found for the Antarctic periphery, the Greenland periphery and the Canadian Arctic, which in our case translates into sea-level equivalents (SLEs) of 61.3 ± 16 mm, 26.9 ± 9.5 mm and 77.6 ± 25 mm, respectively. While our estimates are generally lower, we find in Asia an ice volume 34% larger than that in ref. ⁹, but only 3% larger than that in ref. ⁷. In addition to the new dataset and method presented here, differences in Asia may also come from the use of a more extensive in situ thickness dataset, which provides more data for thicker glaciers³¹,³⁴,³⁵ than previously available, hence improving the calibration regionally. Nevertheless, this result highlights the high level of uncertainty in Asia, where inversions remain poorly constrained.

Due to the limited constraints (for example, surface slope) available at the time to estimate glacier thicknesses, existing models

using flowlines are subject to important biases at the glacier scale, such as discontinuities at the boundaries between adjacent glacier basins or artifacts related to surface slope leading to unrealistic glacier geometries (for example, Vatnojökull Ice Cap, Fig. 2d). Conversely, the use of an inversion of glacier surface flow velocity mapped at regional scale (independent from the glacier inventory catchments) allows us to unveil a new and more physically realistic view of the sub-glacial landscape. Our worldwide ice thickness estimation, hence, reconciles the spatial distribution of the ice masses with the glacier dynamics, morphology and ice divides, highlighting the importance of using ice velocity observations to constrain ice thickness models (Fig. 2a–f). The large differences between our results and existing estimates suggest that similarities in total ice volume are only fortuitous and that no definitive consensus has been reached yet. As a consequence, this new view of glacier geometry calls for a re-evaluation of the evolution of the world's glaciers in numerical models.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Comparison of glacier thickness distribution. **a–f**, Ice thickness for three selected regions: Academy of Sciences Ice Cap (Severnaya Zemlya) (**a,b**), Vatnojökull Ice Cap (Iceland) (**c,d**) and Penny Ice Cap (Baffin Bay) (**e,f**). For each region, two panels show the ice thickness from this study (**a,c,e**) and from the 2019 consensus (**b,d,f**).⁹ Ice thickness is colour coded from white (no ice) to dark blue for ice thickness larger than 700 m. The solid black lines indicate the limits of the RGI glacier basins.</figcaption></figure>

**Implications for water resources.** The contribution of mountain glaciers to water resources is largely uneven across regions due to differences in seasonal precipitation and melting of the snow pack³⁸⁻⁴⁰. To assess the potential impact of the new estimated ice volume on water resources, we focus on six watersheds in the tropical Andes and the high mountains of Asia, where previous studies showed that glaciers contribute more than 5% of the river discharge³⁹,⁴⁰ (Fig. 3). For all considered watersheds, our results further accentuate differences in available regional glacial water resources.

In the upper Indus and Chenab basin in the high mountains of Asia, where more than 8 M people live, glacier meltwater represents more than 50% of the river flow⁴⁰ during the dry season, while the glacier contribution compared with rainfall in other Himalayan regions is relatively small. In this basin, where glaciers are losing mass at a rate of 0.09 ± 0.07 m yr⁻¹ water equivalent (w.e.) (ref. ³), our new estimate shows water reservoirs that are 17–31% larger than previous estimates⁷ (Fig. 3a). Similarly, for sub-basins of the Brahmaputra rivers, which are less densely populated (0.4 M people), our estimates
<figure><img src="image_5.png" /></figure>

**Fig. 3** | Impact of the velocity-derived thicknesses on water resources. a,b, Water resources stored in glaciers for seven watersheds of the Himalayas (a) and South America (b) where glacier contribution to runoff is above 5%³⁹,⁴⁰. Solid black lines represent the limits of the watersheds. Solid blue lines represent main river segments. Number of people per square kilometre inside the watershed is colour coded from white to red. Circle area increases with the watershed population. Circle colour represents available water resources per person per year (assuming a constant mass loss rate until 2100). Percentages indicate the difference in water volume relative to the consensus estimate⁹. Background image is from ArcGIS Online basemaps⁴².

of available water volumes are 30–87% larger (Fig. 3a). These new findings suggest that the timing of the peak water may be delayed in time, thereby temporarily reducing the pressure on freshwater availability in these watersheds.

Conversely, the three catchments considered in the tropical Andes that also provide water for a large share of the population (4M people) have a reduced estimated ice volume by 20% on average, compared with the 2019 consensus (Fig. 3b). In addition, glaciers of these catchments are experiencing among the largest rates of mass loss, at 0.4 m yr⁻¹ w.e. (refs. ³,⁴¹). Our new estimates of glacier volume, therefore, further accentuate the water-pressure needs in the future by bringing the time of the peak water sooner. For example, the drainage basin of La Paz, with more than 2.2 M inhabitants, is strongly dependent on the water from ice melt that represents up to one-third of the water resources during the dry season³⁹. If runoff did not change notably during the past four decades with increasing ice melt compensating the loss of glacierized area³⁹, it seems inevitable that the reduced total reservoir of ice (–23%) will have an impact sooner than anticipated.

**Challenges in constraining ice thickness models.** Continuous observations of glacier surface flow velocity and elevation changes, along with more detailed information on the bedrock topography, are critical to reconstruct the size of ice reservoirs and understand the temporal variability and evolution of the world's glaciers. Our observations of glacier surface flow velocity at global scale allows the use of a physically based approach to reconstruct ice thickness distribution. Using a regional-scale inversion method, our results provide a consistent vision of glacier flow and geometry around the globe, which suggests reconsidering the uncertainties related to current thickness models, and specifically for large ice fields and ice caps. Our reconstruction puts in perspective the high level of uncertainty of glacier thickness models (Supplementary Table 1 and Supplementary Figs. 8–12). Specifically, the accuracy of our estimates still strongly relies on the availability and distribution of in situ thickness measurements, which are lacking for the majority of glaciers⁶. In particular, additional direct observations of ice thicknesses would be beneficial in the Andes, the Himalayas and the Russian Arctic, which are regions hosting 18% of the total volume
and are critical for local populations in terms of water resources. The mapping of glacier surface flow velocity redefines our view of glacier flow dynamics and ice-volume distribution, with important implications for improving reconstructions of past and ongoing changes and especially for modelling the evolution of glacier dynamics and mass balance in a changing climate.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-021-00885-z.

Received: 1 February 2021; Accepted: 16 December 2021; Published online: 7 February 2022

## References

1. Hock, R. et al. in *Special Report on the Ocean and Cryosphere in a Changing Climate* (eds Pörtner, H.-O. et al.) Ch. 2 (IPCC, 2019).

2. Hugonnet, R. et al. Accelerated global glacier mass loss in the early twenty-first century. *Nature* **592**, 726–731 (2021).

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

4. Hock, R. et al. *GlacierMIP—a model intercomparison of global-scale glacier mass-balance models and projections*. *J. Glaciol.* **65**, 453–467 (2019).

5. McGranahan, G., Balk, D. & Anderson, B. The rising tide: assessing the risks of climate change and human settlements in low elevation coastal zones. *Environ. Urban.* **19**, 17–37 (2007).

6. Welty, E. et al. Worldwide version-controlled database of glacier thickness observations. *Earth Syst. Sci. Data* **12**, 3039–3055 (2020).

7. Huss, M. & Farinotti, D. Distributed ice thickness and volume of all glaciers around the globe. *J. Geophys. Res. Earth Surf.* **117**, f04010 (2012).

8. Linsbauer, A., Paul, F. & Haeberli, W. Modeling glacier thickness distribution and bed topography over entire mountain ranges with GlabTop: application of a fast and robust approach. *J. Geophys. Res. Earth Surf.* **117**, F03007 (2012).

9. Farinotti, D. et al. A consensus estimate for the ice thickness distribution of all glaciers on Earth. *Nat. Geosci.* **12**, 168–173 (2019).

10. Huss, M. & Hock, R. A new model for global glacier change and sea-level rise. *Front. Earth Sci.* **3**, 54 (2015).

11. Grinsted, A. An estimate of global glacier volume. *Cryosphere* **7**, 141–151 (2013).

12. Farinotti, D. et al. How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment. *Cryosphere* **11**, 949–970 (2017).

13. Farinotti, D., Huss, M., Bauder, A., Funk, M. & Truffer, M. A method to estimate the ice volume and ice-thickness distribution of alpine glaciers. *J. Glaciol.* **55**, 422–430 (2009).

14. Marzeion, B., Jarosch, A. H. & Hofer, M. Past and future sea-level change from the surface mass balance of glaciers. *Cryosphere* **6**, 1295–1322 (2012).

15. Radić, V. et al. Regional and global projections of twenty-first century glacier mass changes in response to climate scenarios from global climate models. *Clim. Dyn.* **42**, 37–58 (2014).

16. Zorzut, V. et al. Slope estimation influences on ice thickness inversion models: a case study for Monte Tronador glaciers, North Patagonian Andes. *J. Glaciol.* https://doi.org/10.1017/jog.2020.64 (2020).

17. Werder, M. A., Huss, M., Paul, F., Dehecq, A. & Farinotti, D. A Bayesian ice thickness estimation model for large-scale applications. *J. Glaciol.* **66**, 137–152 (2020).

18. Millan, R. et al. Vulnerability of Southeast Greenland glaciers to warm Atlantic water from Operation IceBridge and Ocean Melting Greenland data. *Geophys. Res. Lett.* **45**, 2688–2696 (2018).

19. Morlighem, M. et al. BedMachine v3: complete bed topography and ocean bathymetry mapping of Greenland from multibeam echo sounding combined with mass conservation. *Geophys. Res. Lett.* **44**, 11051–11061 (2017).

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

21. Millan, R., Mouginot, J. & Rignot, E. Mass budget of the glaciers and ice caps of the Queen Elizabeth Islands, Canada, from 1991 to 2015. *Environ. Res. Lett.* **12**, 024016 (2017).

22. Mouginot, J. & Rignot, E. Ice motion of the Patagonian Icefields of South America: 1984–2014. *Geophys. Res. Lett.* **42**, 1441–1449 (2015).

23. Dehecq, A. et al. Twenty-first century glacier slowdown driven by mass loss in High Mountain Asia. *Nat. Geosci.* **12**, 22–27 (2019).

24. Gardner, A. S., Fahnestock, M. & Scambos, T. *ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities* (NSIDC, 2020).

25. Millan, R. et al. Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach. *Remote Sens.* **11**, 2498 (2019).

26. *Randolph Glacier Inventory* 6.0 (RGI Consortium, 2017); https://doi.org/10.7265/N5-RGI-60

27. Michel, R. & Rignot, E. Flow of glaciar Moreno, Argentina, from the repeat-pass shuttle imaging radar images: comparison of the phase correlation method with radar interferometry. *J. Glaciol.* **45**, 93–100 (1999).

28. Mouginot, J., Rignot, E., Scheuchl, B. & Millan, R. Comprehensive annual ice sheet velocity mapping using Landsat-8, Sentinel-1, and RADARSAT-2 data. *Remote Sens.* **9**, 364 (2017).

29. Cuffey, K. & Paterson, W. S. B. *The Physics of Glaciers* (Butterworth-Heinemann/Elsevier, 2010).

30. Gaillardet, J. et al. OZCAR: the French network of critical zone observatories. *Vadose Zone J.* **17**, 180067 (2018).

31. Lambrecht, A., Mayer, C., Aizen, V., Floricioiu, D. & Surazakov, A. The evolution of Fedchenko glacier in the Pamir, Tajikistan, during the past eight decades. *J. Glaciol.* **60**, 233–244 (2014).

32. Björnsson, H. & Pálsson, F. Radio-echo soundings on Icelandic temperate glaciers: history of techniques and findings. *Ann. Glaciol.* https://doi.org/10.1017/aog.2020.10 (2020).

33. Dowdeswell, J. A. et al. Form and flow of the Academy of Sciences Ice Cap, Severnaya Zemlya, Russian High Arctic. *J. Geophys. Res. Solid Earth* https://doi.org/10.1029/2000JB000129 (2002).

34. Azam, M. F. et al. From balance to imbalance: a shift in the dynamic behaviour of Chhota Shigri glacier, western Himalaya, India. *J. Glaciol.* **58**, 315–324 (2012).

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

36. Millan, R. et al. Ice thickness and bed elevation of the northern and southern Patagonian icefields. *Geophys. Res. Lett.* https://doi.org/10.1029/2019GL082485 (2019).

37. Morlighem, M. et al. Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet. *Nat. Geosci.* https://doi.org/10.1038/s41561-019-0510-8 (2019).

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

39. Soruco, A. et al. Contribution of glacier runoff to water resources of La Paz city, Bolivia (16°S). *Ann. Glaciol.* **56**, 147–154 (2015).

40. Lutz, A. F., Immerzeel, W. W., Shrestha, A. B. & Bierkens, M. F. P. Consistent increase in High Asia's runoff due to increasing glacier melt and precipitation. *Nat. Clim. Change* **4**, 587–592 (2014).

41. Dussaillant, I. et al. Two decades of glacier mass loss along the Andes. *Nat. Geosci.* **12**, 802–808 (2019).

42. *World Imagery* (Esri, accessed 19 May 2019); https://www.arcgis.com/home/item.html?id=10df2279f9684e4a9f6a7f08febac2a9

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

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2022, corrected publication 2022
## Methods

**Ice velocity mapping.** *Ice velocity calculation and calibration.* Orthorectified optical satellite images from the American US Geological Survey's LANDSAT 8, the European ESA's Sentinel-2 and the French-Israel Centre National d'Etudes Spatiales-Israel Space Agency's Venus and synthetic aperture radar images from the European ESA's Sentinel-1 were used globally to map glacier surface flow velocity in 2017 and 2018. We formed all possible pairs of images, from the same orbits, with repeat cycles ranging from the sensor nominal cycle (2–16 days) to 100 days and from 330 to 400 days as in ref. ¹⁶. For regions with ice caps and ice fields (for example, Patagonia, Iceland, Canadian Arctic, Svalbard and Russian Arctic) with low correlation for optical data in accumulation regions and dense cloud cover, we also processed Sentinel-1 images (6-day to 24-day repeat cycles) using speckle tracking²¹ to increase the spatial coverage and reduce uncertainties. Sentinel-1's acquisitions used here are made with the interferometric wide-swath mode using progressive scans synthetic aperture radar technique. These observations were processed using the same methodology as ref. ²⁵. We use reference and search sub-images for the correlation with a size of 128×32 pixels, search window of 16×4 pixels and sampling spacing of 28×8 pixels. The geocoding of Sentinel-1 displacement maps uses the TanDEM-X 90 m global digital elevation model.

All glacier surface velocity maps from each image pair are finally calibrated to correct for geometric distortions or geolocation errors. Glacierized areas are masked out, and a linear ramp (Landsat 8), quadratic function (Sentinel-1) or constant value (Sentinel-2) is calculated to best fit the value of the displacement on stable ground for the east-west and north-south components independently²⁵. This fit is finally removed from the entire displacement map to obtain zero-centred surface velocity in ice-free regions.

**Mosaicking.** To reduce computational time, more than 62 smaller sub-regions of the RGI regions were defined (Supplementary Information). For each region, we calculate composite velocity mosaics for each sensor as in ref. ²⁵. For each resulting mosaic, we then perform a spatial adaptive filtering using a median filtering (outlier removal) and a Gaussian smoothing. The size of the kernels, for both the median and Gaussian filters, is adapted depending on the distance to the nearest ice margin to preserve the ice flow structure. Thus, no smoothing is applied when a pixel is closer than 150 m from the ice margin. The smoothing then gradually increases to reach a maximum size of 2 km when farther than 6 km from the ice margin. These smoothing seems to provide the most robust results while preserving the smallest ice flow structures.

On the basis of ref. ²⁵, we expect the precision and accuracy on glacier surface velocity to be globally better than 10 m yr⁻¹, although it might locally be larger due to the limited number of observations (for example, cloud, polar night and so on) or geometrical errors in the original images. Regarding glaciers in the periphery of Greenland, we use the mapping of surface velocity performed by ref. ²⁸, which is a weighted average of radar (Sentinel-1A and 1B, Radarsat, ALOS/PALSAR) and optical data (Landsat 8) at a precision of 18 m yr⁻¹ (ref. ²⁸). These mosaics are available from the National Snow and Ice Data Center²⁸.

**Ice thickness estimation.** *Modelling.* We compute the ice thickness on the basis of surface motion and slopes using the shallow-ice approximation (SIA)²⁹. The SIA (with basal sliding) leads to the following expression for surface velocities:

$$v_s = v_b - 2 (\rho g)^n \|\nabla s\|^n A \frac{H^{n+1}}{n+1} \nabla s \quad (1)$$

where $v_s$ and $v_b$ are the surface and basal velocity, respectively, $n$ is the exponent in the Glen's law usually taken as $n=3$, $A$ is the creep parameter, $s$ is the ice surface elevation ($\nabla s$ being the surface slope), $\rho$ is the ice density taken as 917 kg m⁻³ and $g$ is the acceleration due to gravity, which equals 9.81 m s⁻². On the basis of this expression, we calculate the ice thickness, $H$, directly, using the ice velocity and slope map of an entire region (see section 'Mosaicking') as:

$$H = \left( \frac{(v_s - v_b) (n + 1)}{2A (\rho g)^n \|\nabla s\|^n} \right)^{\frac{1}{n+1}} \quad (2)$$

Because the basal velocities are unknown, we introduce $\beta$ in the interval [0,1] that is the ratio between basal and surface velocity such that $v_b = \beta v_s$, which yields:

$$H = \left( \frac{v_s (1 - \beta) (n + 1)}{2A (\rho g)^n \|\nabla s\|^n} \right)^{\frac{1}{n+1}} \quad (3)$$

To be consistent with the consensus estimate, we estimated ice thicknesses for all glaciers with a degree of connectivity to the ice sheet <2, as defined by the RGI.

**Calibration.** The two parameters that need to be calibrated in equation (3) are $A$ and $\beta$. In regions where ice thicknesses have been measured (GlaThiDa,

GLACIOCLIM³⁰ and other sources³²⁻³⁶), we calibrate $A$ using the expression of equation (4) (see section 'Limitations' for $\beta$ values):

$$A = \left( \frac{v_s (1 - \beta) (n + 1)}{2 (\rho g)^n \|\nabla s\|^n H^{n+1}} \right) \quad (4)$$

In most cases, for a specific region, we calculate an average $A$ value using all available observed in situ ice thicknesses within that region and use it to estimate the ice thickness in equation (3) (Supplementary Table 4). Over the European Alps, a large quantity of radar data have been acquired, covering thicknesses ranging from 10 to >400 m. In this region, we calculate a parameterization of $A$ using a quadratic function where $A = f(H)$ because it allows us to better fit large ice thickness data that were underestimated when using a constant $A$ value (Supplementary Fig. 1). Finally, when no data are available to calibrate the ice thickness, we either use the mean $A$ of the closest region or use the standard value of $A = 24 \times 10^{-25} \text{ s}^{-1} \text{Pa}^{-3}$, which corresponds to the typical value for temperate ice²⁹. Concerning the $\beta$ parameter, we use a parameterization based on ice velocity and surface slope (see section 'Limitations').

In the case of glaciers in the periphery of Antarctica (RGI-19), we used the ice thickness reconstruction from BedMachine v.1 that uses a kriging and mass-conservation approach¹⁷. Since Antarctica has been widely and densely covered by ground-penetrating radar surveys since the 1990s, this solution provides a better estimate than our modelling approach. Region RGI-19 was then completed with our own reconstruction of ice thicknesses for islands in the periphery of Antarctica (Joinville Islands, South Georgia, South Orkney, Kerguelen, Heard and McDonald islands). Finally, as the degree of connectivity in Antarctica has not been assigned in the RGI, we remove all glaciers connected directly to the ice sheet or ice shelves (Supplementary Fig. 13), which should correspond to a strong level of connection as defined by the RGI Consortium (value of 2).

**Model validation.** As the in situ data are used to calibrate the model, we also want to test that the resulting mean difference between the model and the in situ measurements is not a result of overfitting. To do so, we have modelled the ice thickness by removing 60% of the thickness lines used for calibration in the Alps, where an important amount of in situ data exists both over thin (10–50 m) and thicker (100–500 m) ice and compare the modelled result with the withdrawn measurements. This comparison is represented in Supplementary Fig. 2. We found a difference of -16 ± 51 m between the model and the withdrawn in situ measurements, which is comparable to the statistics derived from the model that is calibrated using all observations (Supplementary Table 2). This result indicates that our model is not overfitting field observations.

We further analysed the difference with observed thickness data depending on different GlaThiDa thickness bins (Supplementary Fig. 3). On average, errors for ice thickness greater than 100 m are in the range of 25–35%. For ice thicknesses below 100 m, we assume error >50%. This analysis is further used to provide error maps on the modelled ice thickness, which will be provided through a public data portal.

**Limitations.** The model used to calculate ice thickness is strongly dependent on (1) the surface slope and (2) surface flow velocity. Hence, an active surge phase happening during the study period would lead to an erroneous estimation of the ice thickness. Indeed, a tenfold increase in glacier surface flow velocity would lead to an overestimation of the ice thickness by more than 75%, as defined in equation (3). Furthermore, the use of the SIA in the case of a surge is limited since it accounts only for the amount of internal ice deformation, and indirectly of sliding, which is dominant during a surge. In addition, the temporal mismatch between the datasets (see the following), which might not be consistent between the active and quiescent phases of surging glaciers, probably further increases the uncertainty of our modelled thickness distribution. To reflect this uncertainty in ice thickness distribution, we assembled a glacier surge inventory on the basis of the pattern of elevation changes from ref. ² for the period 2015–2019, which encompasses our study period (Supplementary Table 3). This identification is based on the expected expression on the surface elevation of the surge as described by ref. ⁴³, for example, thinning of the accumulation area (upper part of the basin) and thickening of the ablation (lower part) due to the rapid transfer of mass between the zones. We identified 98 glaciers worldwide, which represents ~2% of the global volume and ice area of the RGI. The glaciers are located in Alaska (4), Canadian Arctic (6), High Mountain Asia (73), Russian Arctic (1) and Svalbard (14). For those surging glaciers, we increased the uncertainty assessment to 75% on the basis of the empirical assessment of the impact of a tenfold increase in surface velocity on the ice thickness estimation (see the preceding).

Furthermore, the temporal mismatch between the acquisition of (1) the ice velocity product, (2) the Digital Elevation Model (DEM) used to calculate the surface slopes, (3) the in situ ice thickness measurements and (4) the contours of the RGI outlines (representative of the glacier state in 2003³⁰) are sources of uncertainties that are difficult to quantify accurately. Indeed, while our ice velocity product is precisely time tagged, tracking the date of the different DEM sources
(see section 'DEM and surface slope') and the acquisition of in situ ice thickness measurements for each glacier is far more complex, if not impossible⁶,³⁸.

Since the SIA describes primarily vertical shear, strong limitations are also found when sliding becomes dominant (for example, fast marine terminating glaciers). Further limitations are also encountered when (1) the surface slope of the ice is nearly flat, which tends to overestimate ice thicknesses, and (2) the surface velocity exceeds the surface slope by several orders of magnitude. For fast marine terminating cases, both limitations are encountered, which further increases the ice thickness overestimation.

To mitigate this effect, and to better account for sliding speed $v_b$ versus internal ice deformation, we calculate empirical maps of the $\beta$ parameter on the basis of the observed ratio between the slope and the surface velocity in the interior of ice caps versus in fast-flowing regions (for example, outlet and marine terminating glaciers). In general, ice thicknesses were found to be overestimated when the value of the ratio dropped (typically <0.001 yr m⁻¹) compared with the value of the ratio on an icefield plateau (Patagonia) or the centre part of an ice cap (Russian Arctic) (Supplementary Fig. 4). As a consequence, we set the $\beta$ parameter to be equal to 0.1 (10% of sliding) in areas where the ratio is greater than 0.001 (icefield plateau, centre of ice caps) and modulate this value down to 0.9 (90% of sliding) where the value of the ratio was less than 0.001 (for example, outlet and marine terminating glaciers). An example of a $\beta$ map is represented in Supplementary Fig. 5.

**DEM and surface slope.** Three different sources of DEM were used to estimate the surface slope: (1) the ASTER Global DEM version 3 (GDEM)⁴⁴, (2) the 90-m TanDEM-X global DEM (TDX)⁴⁵ and the Worldview Arctic DEM (WV)⁴⁶.

In the case of the Himalayas, Patagonia, Arctic and pan-Antarctic regions, we first calculated the glacier surface slope from TDX and GDEM and the slope roughness, and we then used the slope value of the DEM that has the lowest roughness. Similarly, WV in combination with GDEM was used for the Canadian Arctic and Alaska.

The SIA is based on the assumption that at the scale of multiple ice thicknesses (10–20 times²⁹,⁴⁷) the flow is linked to the surface slope. To obtain the glacier surface slope at multiple ice thicknesses, we have filtered the surface slope using a spatially dependent Gaussian kernel. As ice thicknesses are unknown, we used the distance from the ice margin as a proxy. As a result, the kernel size varies gradually from 1.5 km at the ice margin to 8.0 km when the distance to the ice margin is larger than 8 km. Ice-free areas are masked during the filtering using the RGI outlines (Supplementary Fig. 6). Finally, we used a lower limit for the surface slope at 1° where ice velocities were less than 10–20 m yr⁻¹, which was set up empirically using observed values of surface slope and velocities at the ice divides of ice caps (for example, Russian, Canadian Arctic, Spitzberg). This allows us to avoid singularities when slope tends towards 0, leading to ice thickness that tends towards infinite values, as defined in equation (3).

**Calculation of SLE.** To calculate SLE, we first calculate the total ice volume in cubic kilometres as:

$$V = \sum \frac{h_i \times \text{psx} \times \text{psy}}{1 \times 10^9} \qquad (5)$$

where psx and psy are the pixel size in metres in the x and y direction, respectively; $h_i$ is the ice thickness at a pixel location. We then account for only the ice above floatation to obtain the volume of water above flotation:

$$V_{\text{af}} = V_{\text{total}} \frac{\rho_{\text{ice}}}{\rho_{\text{water}}} - V_{\text{bsl}} \qquad (6)$$

where $V_{\text{af}}$ is the volume of water above flotation, $V_{\text{total}}$ is the total ice volume, $\rho_{\text{ice}}$ is the density of the ice (0.9167 Gt km⁻³), $\rho_{\text{water}}$ is the density of seawater (1.027 Gt km⁻³) and $V_{\text{bsl}}$ is the volume located below sea level using the geopotential model EGM96 as reference. The error on the total volume is calculated by integrating the error on the ice thickness (see section 'Model Validation') over each RGI region. Finally, we convert the volume to SLE using the ice volume

needed to raise global sea level by 1 mm (361.8 Gt of ice will raise the global sea level by 1 mm):

$$SLE_{(\text{mm})} = V_{\text{af}} \frac{\rho_{\text{water}}}{361.8} \qquad (7)$$

## Data availability

Thickness and ice velocity datasets are publicly available at https://doi.org/10.6096/1007. GLACIOCLIM data available at https://glacioclim.osug.fr/; GlaThiDa data available at https://www.gtn-g.ch/datacatalogue_glathida/; RGI available at https://www.glims.org/RGI/; river basins and population data available at www.worldpop.org and https://www.bafg.de.

## References

43. Meier, F. & Post, A. What are glacier surges? Can. J. Earth Sci. 6, 807–817 (1969).

44. Abrams, M., Crippen, R. & Fujisada, H. ASTER Global Digital Elevation Model (GDEM) and ASTER Global Water Body Dataset (ASTWBD). Remote Sens. 12, 1156 (2020).

45. TanDEM-X - Digital Elevation Model (DEM) - Global, 90 m (German Aerospace Center, 2018); https://doi.org/10.15489/JU28HC7PUI09

46. Porter, C. et al. ArcticDEM (Polar Geospatial Center, 2018); https://doi.org/10.7910/DVN/OHHUKH

47. Bamber, J. L., Hardy, R. J. & Joughin, I. An analysis of balance velocities over the Greenland ice sheet and comparison with synthetic aperture radar interferometry. J. Glaciol. 46, 67–74 (2000).

## Acknowledgements

R.M. acknowledges support from a post-doctoral fellowship awarded by the French Centre National d'Etudes Spatiales (CNES). J.M. and A.R. acknowledge support from the CNES MaiSON project and Labex OSUG@2020 (Investissements d'avenir - ANR10 LABX56). M.M. acknowledges support from MEASURES-3 project (NASA grant 80NSSC18M0083). We gratefully acknowledge CNES, ESA, Copernicus Program, NASA and Deutsches Zentrum für Luft- und Raumfahrt e.V. for providing the observation used in this paper (Landsat 8, ASTER, Sentinel-1 and 2, TanDEM-X, Venus). We thank the RGI and GlaThiDa, the Service National d'Observation GLACIOCLIM, the worldpop and BFG for the glacier basins, ice thickness, river basins and population data. The computing/storage resources (about 4 MCPU-hours and 100 Tb) used for this work were provided by high-performance clusters from GRICAD (Grenoble Alpes Recherche - Infrastructure de Calcul Intensif et de Données). Maps in Fig. 3 were created using ArcGIS software by Esri. ArcGIS and ArcMap are the intellectual property of Esri and are used herein under license. Copyright © Esri. All rights reserved. For more information about Esri software, please visit www.esri.com.

## Author contributions

R.M., J.M., A.R. and M.M. conceived and designed the research. R.M. and J.M. performed the experiments. R.M., J.M. and A.R. analysed the data. R.M., J.M., A.R. and M.M. wrote the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41561-021-00885-z.

**Correspondence and requests for materials** should be addressed to Romain Millan.

**Peer review information** Nature Geoscience thanks Leigh Stearns, Adrian Luckman and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Tom Richardson.

**Reprints and permissions information** is available at www.nature.com/reprints.