# Drivers of global glacial erosion rates

Sophie L Norris, John C Gosse, Romain Millan, Jeremie Mouginot, Antoine Rabatel, Mathieu Morlighem, Matthew S M Bolton, Richard B Alley

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

## To cite this version:

Sophie L Norris, John C Gosse, Romain Millan, Jeremie Mouginot, Antoine Rabatel, et al.. Drivers of global glacial erosion rates. Nature Geoscience, 2025, 18 (8), pp.732-739. 10.1038/s41561-025-01747-8. hal-05376955

HAL Id: hal-05376955

https://hal.science/hal-05376955v1

Submitted on 21 Nov 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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Drivers of global glacial erosion rates

Received: 11 February 2022

Accepted: 20 June 2025

Published online: 7 August 2025

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

Check for updates

Sophie L. Norris Ⓓ1,2✉, John C. Gosse Ⓓ1✉, Romain Millan Ⓓ3, Jeremie Mouginot Ⓓ3,8, Antoine Rabatel Ⓓ3, Mathieu Morlighem Ⓓ4, Matthew S. M. Bolton Ⓓ5 & Richard B. Alley Ⓓ6,7

Glacial erosion must be quantified to better constrain numerous geomorphic and orogenic processes; however, accurate models of glacial erosion have been limited by sparse data. Here we use machine learning tools to develop equations that integrate glacial erosion and glaciological, topoclimatic and geological variables based on a global-scale synthesis of 181 contemporary glacier-derived erosion rates. The results reveal environment-specific erosion rate equations for surge-type, marine- and land-terminating glacial settings. We demonstrate that glacial velocity is not the most statistically important predictor of glacial erosion in any environment. Instead, an improved prediction of glacial erosion is attained when velocity is considered with additional glaciological, topoclimatic and geological variables, with the most dominant influences exhibited by precipitation, glacial elevation, length, latitude and the underlying geology. Using these equations, we estimate erosion rates for 85% of contemporary glaciers, with 99% eroding between 0.02 and 2.68 mm yr⁻<fnref n="1" />. Our results suggest a need to adjust how we predict or hindcast glacial erosion rates and highlight their sensitivity not only to changes in glacial sliding velocity but also to additional glaciological, topoclimatic and geological influences.

Throughout the Quaternary period (the past 2.58 Myr), glaciation has acted as a dominant control in sculpting landscapes on Earth<fnref n="1" />,². Understanding of the variability in rates and distribution of glacial erosion is needed to assess its influence on many current interdisciplinary societal questions regarding ocean and riverine sedimentation and alkalinity³, the evolution of mountain topography<fnref n="1" />, and the long-term management of used nuclear fuel in geological repositories⁴. Despite the importance of glacial erosion, the physical processes of quarrying, abrasion and subglacial meltwater incision that govern its magnitude and distribution are understudied⁵⁻⁸. Furthermore, owing to the difficulty of sampling in subglacial systems, observations of these processes are limited⁷,⁸. In recent decades, numerous indirect measurements have supplemented these in situ observations under active glaciers to understand material removal better. These methods commonly

comprise sediment discharge estimates from meltwater streams, the volume of bulk sediment stored proglacially and the depth of bedrock removal constrained by geochronological approaches⁵,⁷⁻¹³. Collectively, estimates of glacial erosion using these methodologies record rates that span multiple orders of magnitude²,⁷,⁹,¹²,¹⁴.

Physical reasoning indicates that glacial erosion rates are influenced by many ice, substrate and climate characteristics, but ice-flow velocity is likely to exert notable control (for example, refs. 6,8). Pioneering efforts to derive erosion models from measurements spanning the past ~10,000 years have found dependence on glacial velocity expressed as linear or exponential 'velocity-erosion rules'²,⁵,⁷,⁸,¹¹,¹⁵. Since the introduction of the first velocity-erosion rule, which was primarily based upon observations of an Alaskan surge-type glacier¹⁵, similar velocity-erosion rules have been independently produced from

<footnote n="1">¹Department of Earth and Environmental Sciences, Dalhousie University, Halifax, Nova Scotia, Canada. ²Department of Geography, University of Victoria, Victoria, British Columbia, Canada. ³Université Grenoble Alpes, CNRS, IRD, INRAE, Grenoble INP, IGE, Grenoble, France. ⁴Department of Earth Sciences, Dartmouth College, Hanover, NH, USA. ⁵Department of Earth and Atmospheric Sciences, University of Alberta, Edmonton, Alberta, Canada. ⁶Earth and Environmental Systems Institute, Pennsylvania State University, University Park, PA, USA. ⁷Department of Geosciences, and Earth and Environmental Systems Institute, Pennsylvania State University, University Park, PA, USA. ⁸Deceased: Jeremie Mouginot. ✉e-mail: sophienorris@uvic.ca; john.gosse@dal.ca</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 1** | Distribution of contemporary glacial erosion rates. **a**, Spatial distribution of synthesized glacial erosion rates (see Supplementary Table 1 for referenced datapoints). Regions and their numbers correspond to the RGI²⁵ regions. The total number (no.) of glacial erosion rate measurements and RGI region mean (μ) are displayed. In cases where a circle represents more than</figcaption></figure>

one erosion rate and a majority exists, colour coding represents the mode **b**, Latitudinal distribution of synthesized glacial erosion rates. Each point corresponds to an individual measurement of glacial erosion. Basemap in **a** from Natural Earth (https://www.naturalearthdata.com).

marine-terminating glaciers covering a transect from Patagonia to the Antarctic Peninsula⁵, New Zealand's Franz Josef Glacier¹¹, and global datasets that include land, marine and surge-type glaciers⁷,⁸, lending support for a relationship between glacial velocity and erosion. While velocity-erosion rules do not directly reveal subglacial processes⁶,⁸, landscape-evolution models have effectively implemented them to duplicate long-term topographic evolution under glaciers (for example, refs. 1,16–22).

Despite wide implementation of velocity-erosion rules, multiple observations support continuing efforts to expand these rules to include other glaciological, topoclimatic and geological controls⁶. For example, in environments with almost identical glacial sliding velocities, erosion rates have been shown to differ hugely⁵,⁷,¹⁰, with strong positive relationships between glacial erosion and latitude (acting as a proxy for temperature)⁵ and precipitation⁷. More complex thermomechanical ice-sheet and landscape-evolution models support multivariable control of glacial erosion predictions by providing mechanistic links between glacial erosion, lithology, tectonics and climate on geological timescales¹³,²³,²⁴. Additionally, earlier velocity-erosion rules often focused on relatively high-velocity glaciers in surging¹⁵ or marine-terminating⁷ environments likely to be especially geomorphically active⁷; however, recent work has noted limited correlations between velocity and erosion rates, especially in regions that yield minimal rates of erosion⁵,⁷, suggesting other controls, in addition to velocity, may be important.

Here we assemble a global synthesis of 181 published glacial erosion rates (previous velocity-erosion rules used *n* = 3–38) estimated using both in situ and indirect measurements derived from subglacial material removal and sediment export rates. We also compile a broad range of possible glaciological, topoclimatic and geological controlling variables to explore the dependency of glacial erosion rates using machine learning techniques. In line with previous erosion syntheses⁷, we compile erosion measurements from topographically constrained

extant glaciers and ice caps (separated into glacier basins²⁵) characterized by dominantly warm basal thermal regimes (164 erosion rates are averaged over the past ~350 years or less and 17 span the early to late Holocene; Fig. 1 and Supplementary Table 1). Our full dataset shows no strong method-based biases related to measurement technique (Supplementary Figs. 1 and 2) or significant ‘timescale dependence’¹⁴ (coefficient of determination *R*² = 0.0097; see Methods and Supplementary Fig. 3).

Using this synthesis, we first test predictions of previous velocity-erosion rules. We then produce a stratified statistical model that best predicts an erosion rate for each glacier. We use elastic net regression²⁶, a linear-regression approach that reduces variance by combining two regularization techniques (L1 and L2 penalties) to shrink coefficients, select variables and cope with collinearity (Fig. 2). This way, we provide multivariable empirical equations integrating glacier surface and sliding velocity²⁷ absolute latitude²⁵, geothermal heat flux (GHF)²⁸, mean annual precipitation (MAP)²⁹,³⁰, mean annual air temperature (MAAT)²⁹,³⁰, seismicity (as peak ground acceleration (PGA); Methods)³¹, ice thickness²⁷, glacier length²⁵, area²⁵, surface slope²⁵ and elevation²⁵. Where available, measurements for these variables were taken from direct observations or from global datasets that broadly temporally correspond with the collection interval of an erosion rate (Methods); however, we acknowledge that when temporal discrepancies between erosion rate measurements and explanatory variables occur, spuriously low apparent predictor influence could result. Using these new equations, we provide a near-global estimate of contemporary glacial erosion rates.

## Importance of glacial velocity in predicting glacial erosion

We test previous velocity-erosion rules using our global synthesis of 181 glacial erosion rates, between 0.01 and 370 mm yr⁻¹ (Fig. 1 and Supplementary Table 1). These velocity-erosion rules, estimated
a

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

c

<figure><img src="image_7.png" /><figcaption>**Fig. 2 | Correlation matrix (Pearson correlation coefficient, r) of continuous explanatory variables in different glacier environments.** a, Surge-type environment. b, Land-terminating environment. c, Marine-terminating environment. d, All datapoints regardless of the glacier environment. Cells with insignificant *P* values (the probability of obtaining a correlation coefficient as</figcaption></figure>

b

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

d

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

large or larger than we observed by chance if the true correlation was zero (the null hypothesis)) are marked with X at a significance threshold of 0.05. Variables are ordered based on similarity following hierarchical cluster analysis with the complete linkage method. Max., maximum; med., median; precip., precipitation; slid., sliding; surf., surface; temp., temperature; vel., velocity.

from a non-random subset of glaciers included in our dataset, are not especially effective in predicting erosion rates from velocity across the whole dataset, with R² (coefficient of determination) values ≤0.19 and mean absolute error (MAE, the average absolute difference between observed and predicted values) higher for all datasets than for our preferred models (Fig. 3). Note, here and throughout, when we refer to R² or the root mean squared error (RMSE, square root of the average of the squared differences between observed and predicted

values), we refer to the statistic calculated on the scale of model fitting (that is, the natural logarithm of erosion rate). In contrast, MAE is in linear erosion rate units (mm yr⁻¹) to ensure the result is interpretable. Not surprisingly, multiple-variable equations developed here from the whole dataset generally perform much better (Fig. 3). The elastic net regression employed is partly analogous to exponential or power law regression by including both logged and linear versions of our velocity variables, so we are implicitly testing the same formulation as previously
<figure><img src="image_10.png" /><figcaption>**Fig. 3** | **Glaciological, topoclimatic and geological controls on glacial erosion rates.** a–c, Predicted versus measured glacial erosion rate in surge-type (a), land-terminating (b) and marine-terminating (c) environments based on equations presented in this study. These data are compared with previously proposed velocity–erosion rules, based only on previously published datasets<sup>5,7,11,15</sup>. MAE values are the mean absolute error of predictions relative to observations, while R<sup>2</sup> is the coefficient of determination for the natural logarithm of erosion rate measurements and predictions. d, Variable importance in surge-type, marine- and land-terminating environments. Cross-hatched variables are excluded</figcaption></figure>

owing to insufficient data. Scores are the scaled coefficient of each explanatory variable, such that the most important (that is, largest coefficient) is 100 for each bootstrap resample (n = 1,000); bars represent the mean variable importance across all resamples. Thick and thin error bars show the approximate ±1 standard deviation (15.87th–84.13th percentile) and ±2 standard deviation (2.28th–97.72nd percentile) ranges of the bootstrapped distributions, respectively. X markers indicate the median importance from the resampling, while hollow circles indicate variable importances from the final trained model. Each continuous variable is presented as both ‘linear’ and ‘log-transformed’.
published erosion rules. The individual linear trends between the natural log of erosion rate and the predictors, including velocity, are visualized in Supplementary Fig. 4. The leave-one-out cross-validation (LOOCV) performance of our full variable elastic net models against comparable metrics for models trained on only velocity-based explanatory variables reveals that our selected models perform better than the velocity-only counterparts (Supplementary Tables 3 and 4). None of the velocity-only synthesis models yield $LOOCV R^2$ higher than 0.39, while our land-terminating glacier model shows higher coefficients of determination. In land-terminating glaciers we observe $LOOCV R^2 = 0.55$ (RMSE = 1.20; compared with $R^2 = 0.27$ and RMSE = 1.51 for the velocity-only synthesis dataset). LOOCV performance is challenging to evaluate for surge-type and marine-terminating glaciers, given the relatively low predictive skill of all models of these types and the wide range of erosion rates in these subsets. We do not see improvements in $R^2$ in the complete surge and marine models; however, our preferred elastic net regression models exhibit lower MAE (34.19 mm yr$^{-1}$ (surge); 9.21 mm yr$^{-1}$ (marine)) than the velocity-only synthesis model (MAE = 40.23 mm yr$^{-1}$ (surge); 213.08 mm yr$^{-1}$ (marine)). This difference for surge-type and marine-terminating glaciers indicates these glacier types are better predicted on average using our full elastic net regression model (Supplementary Tables 3 and 4).

## Variables controlling glacial erosion rates

We measure collective variable importance for predictors as the sum of final model effects (combining log and linear versions) of variables (dots in Fig. 3d). Our full-predictor elastic net regression analysis finds the strongest predictors of glacial erosion rates in surge-type settings are MAP and glacier length (comprising 46% and 54% of the total variable importance, respectively; Fig. 3d). In contrast, 85% of the erosion signal in land-terminating environments is explained by glacier elevation, median sliding velocity, glacier length, absolute latitude, median surface velocity, MAP, glacier area and seismicity. Interestingly, glacier velocity (maximum and median, sliding and surface velocity) imparts reduced predictive influence in our analysis of surge-type and land-terminating environments (velocity accounts for 27% of the modelled effect in land-terminating and negligible effect in surge-type environments; Fig. 3d). This contrasts with marine-terminating environments, where 83% of the erosion signal results from MAP, GHF and absolute latitude.

### Glaciological controls

Past work has often started with glacial velocity as a predictor of erosion$^{5,7,11,15}$ based on physical reasoning, as noted above. In our analysis, velocity is not the dominant predictor of erosion rate for any environment, although it seems more important for erosion in land-terminating glaciers. Interestingly, our dataset's median sliding velocity is substantially (240%) higher than the global median of median velocities$^{27}$. Similarly, previous studies that presented glacial erosion as explained by glacial velocity were composed of or included multiple marine-terminating glaciers exhibiting high velocities$^{5,7}$ or have focused on high-velocity glaciers$^{11,15}$ that, in some cases, were also undergoing short periodic surge events$^{15}$. Based on these observations, we suggest that while non-velocity predictors appear to strongly influence erosion, given the significant positive correlations between velocity variables and the natural logarithm of erosion rate (Fig. 2), high-velocity environments are also most likely to exhibit high erosion rates, although other variables seem to be more consistent predictors in our analysis.

Including glaciological, geological and topoclimatic variables improves our erosion rate predictions relative to velocity-only erosion rules in all three glacial environments (Fig. 3a–c), echoing results from process-based landscape-evolution modelling$^{13}$. Field data and process understanding document how velocity and erosion

rate can be decoupled$^6$. In particular, a subglacial till layer between ice and bedrock smooths the bed and protects underlying bedrock from further glacial erosion$^{6,22,32,33}$ while allowing glacial movement by subglacial deformation$^{6,22}$. While we are not modelling this explicitly, by including other glaciological, topoclimatic and geological variables, the effects of a subglacial till layer may be partly incorporated by proxy, improving predictions of glacial erosion rate. For example, research has demonstrated that increases in MAP, glacial melt or MAAT may increase the rate of sediment removal by subglacial meltwater$^{34–36}$, which in turn could remove the subglacial till layer, thus enabling incision into previously protected bedrock$^6$. However, to validate this conclusion, direct measurements of glacial erosion rates derived exclusively from bedrock erosion beneath a glacier rather than sediment volumes removed from the glacier, in addition to measurements of the thickness of the subglacial till layer, would be needed.

### Topoclimatic controls

In all three glacial environments, MAP displays a positive correlation with glacial erosion rate (Fig. 2), supporting a connection between glacial erosion and climate$^{7,10,11,37}$. This relationship has been previously highlighted in sedimentary$^{10}$ and contemporary records$^{37}$, and recently empirically linked to glacial erosion$^7$. If MAP effects are separated into linear and log versions of explanatory variables, the share of variable importance for marine-terminating glaciers is strongly weighted towards MAP over log MAP, lending support to previous observations made on a global scale$^7$ (Supplementary Fig. 4). However, log MAP is favoured in land-terminating and surge-type environments, suggesting a power law relationship over the exponential relationship previously proposed$^7$. We thus support the view that higher MAP values may result in an increased surface slope due to snow accumulation and/or a greater subglacial meltwater flux$^7$, both of which would increase the rate of basal sliding and, in turn, abrasion and quarrying$^{5,7,11,23,37}$. Owing to the strong relationship between MAP and erosion, we suggest that in addition to affecting glacial sliding velocity, MAP must also modulate the subglacial erosion rate by other mechanisms. For example, as previously highlighted$^{5,6,11}$, high MAP may increase the sediment removal and transport rate by increasing subglacial water flows, which otherwise would inhibit abrasion and quarrying, and may also facilitate bedrock erosion by carving subglacial channels.

Previous studies have shown that, like MAP, changes in MAAT (and its effects, proxied indirectly through elevation and latitude) can affect a glacier's dynamics and basal thermal conditions and, therefore, its ability to erode$^{5,23}$. On a global scale, we find that MAAT and, to a greater extent, elevation and latitude strongly predict glacial erosion. However, it appears that these topoclimatic variables exert a secondary control compared with MAP. In line with previous studies$^{5,7}$, we speculate that their lesser influence is because we examine erosion in dominantly warm basal thermal regimes. In the majority of cases, these glaciers exceed the thermal boundary conditions needed to facilitate basal ice movement, and thus, an increase in temperature (surface or geothermal) is far less important than a proportional increase in precipitation$^{5,7,23}$.

### Geological controls

The lithology, bedrock structure and seismicity of a catchment together exert influence on erosion rates in all three glacial environments (Fig. 3). In land-terminating environments (and in some boot-strap resamples of marine and surge-type environments), seismicity (active regions defined as those with PGA values >1.7 m s$^{-2}$; Methods) is a useful predictor of erosion, a result for glaciers that mirrors results in fluvial settings$^{38}$. In the fluvial context, it has been proposed that seismically active regions may display weakened bedrock conditions and, therefore, enhanced erosion associated with these fractured
<figure><img src="image_11.png" /></figure>

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

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

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

**Fig. 4 | Predicted contemporary glacial erosion rates displayed on a near-global scale.** *a*, Global variability in predicted glacial erosion rates for RGI²⁵ regions. Owing to incomplete surface velocity and ice thickness data necessary to predict glacial erosion rates, all glaciers in Antarctica (part of RGI region 19) are excluded. **b–d**, Regional variability in predicted glacial erosion rates in the

southern Andes (RGI region 17; **b**), Alaska, western Canada and northwestern USA (RGI zones 1 and 2; **c**), and South and Central Asia (RGI zones 13, 14 and 15; **d**). Grid cells are assigned based on the most-frequent erosion rate value, as classified in the key, within each cell (100-km pixels for **a**; 50-km pixels for **b–d**). Basemaps from Natural Earth (https://www.naturalearthdata.com).

surfaces³⁸⁻⁴⁰. Within a glacial context, studies have linked increased fracture density and favourable orientation to increased erosion by quarrying⁹,⁴¹⁻⁴³, highlighting these as controls on valley formation (for example, refs. 9,44). In line with this work, if PGA is related to fracture density, regions with higher PGA may exhibit higher erosion rates through quarrying. Additionally, we observe that variation in lithology also imparts a control in our largest and most lithologically balanced subset, land-terminating environments. We did not assess the influence of lithology on marine-terminating or surge-type glacier erosion, owing to the unequal distribution of these subsets' lithologies. However, this does not exclude lithology as an erosion control,

and we hypothesize that a similar effect may be observed in surge-type and marine-terminating environments if comprehensive lithology data could be examined. In addition to the roles exerted by lithology and bedrock structure, globally we observe a positive relationship between glacial erosion and GHF using contemporary glacial erosion data. Our results, therefore, support a recently proposed association among these variables documented using artificial erosion surfaces in a landscape-evolution model²⁴ and are consistent with observations in geothermally active regions (for example, ref. 45), which document increased subglacial meltwater production and basal temperatures, which can, in turn, affect ice-flow dynamics²⁴.
## Global distribution of glacial erosion

Using our multiple-variable equations determined from 181 contemporary glacial erosion rates in combination with available glaciological, topoclimatic and geological variables, we estimate erosion rates for 185,081 contemporary glaciers (85% of Randolph Glacier Inventory version 6 (RGI)²⁵; Supplementary Table 2.) On this near-global scale, we estimate that 99% of these glaciers display predicted erosion rates between 0.02 mm yr⁻¹ (95% confidence interval (CI): 0.020–0.021) and 2.68 mm yr⁻¹ (95% CI: 2.619–2.747). Our best estimate suggests that these glaciers collectively erode approximately 23 Gt of bedrock material each year (Supplementary Fig. 5) if all material is evacuated from beneath each glacier. By comparison, fluvial sediment flux from exorheic drainage areas has been estimated to generate -18.5–20 Gt yr⁻¹ (refs. 46–48).

In total volume terms, the majority of glacial erosion occurs above 50°N/S in regions with widespread modern ice cover in Alaska, the Canadian Arctic, Greenland, Scandinavia and the southern Andes (Fig. 4 and Supplementary Figs. 6 and 7). The highest collective rates of glacial erosion are observed in areas of Alaska, Central and South Asia (RGI regions 13–15 (ref. 25)), Caucasus and the Middle East, and New Zealand (regional medians of ~0.5–1.5 mm yr⁻¹; Supplementary Table 5). This global distribution is predominantly driven by the noteworthy effects of precipitation, glacial elevation, length, latitude and the underlying geology, in addition to glacial velocity.

## Challenges in predicting glacial erosion

We acknowledge that the new insights and improved understanding of glacial erosion we develop here are derived from a small population and thus have limitations. First, our erosion rates are derived using a variety of methodologies, many of which reconstruct erosion indirectly from subglacial sediment export rates. While we detect no strong method-based biases in our dataset (Supplementary Figs. 1 and 2), we caution that these methodological differences may influence measured erosion rates. Second, as highlighted by previous studies, erosion rates used to derive our equations have been collected during contemporary accelerated glacial retreat conditions²,⁵,¹⁴,⁴⁹,⁵⁰, which almost certainly increases sediment fluxes⁶,⁴⁹. Therefore, consistent with other studies (for example, refs. ¹⁴,⁴⁹), our equations should not be used to estimate erosion rates on longer timescales or under past or future glaciological or climatic conditions; however, the equations may be pertinent, especially when modelling past and future periods of rapid deglaciation. Lastly, because, to our knowledge, only ref. ⁵¹ has measured glacial erosion where sliding velocity is minimal (< 0.01 m yr⁻¹), we excluded this datapoint and focused our analysis on predominantly warm-based environments. However, as no global database of cold-based glaciers is available, we cannot remove them from our predictions, and in these cases we recommend treating our predicted glacial erosion rates as maximum values only (Methods).

Despite these limitations, our multiple-variable equations indicate that glacial erosion is predicted better using a collection of glaciological, topoclimatic and geological variables, rather than by glacial velocity alone. We provide further evidence to support the idea that MAP influences erosion rates in all glacial environments and that MAAT (and its effect observed through glacial elevation and latitude) is also influential, although to a lesser extent. Moreover, our results highlight the importance of seismicity, lithology and GHF in controlling glacial erosion. In summary, our study emphasizes the need for a holistic assessment of the controls on glacial erosion that includes glaciological, topoclimatic and geological influences.

## Online content

Any methods, additional references, Nature Portfolio 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-025-01747-8.

## References

1. Egholm, D. L., Nielsen, S. B., Pedersen, V. K. & Lesemann, J. E. Glacial effects limiting mountain height. *Nature* **460**, 884–887 (2009).

2. Koppes, M. N. & Montgomery, D. R. The relative efficacy of fluvial and glacial erosion over modern to orogenic timescales. *Nat. Geosci.* **2**, 644–647 (2009).

3. Lehmann, N. et al. Alkalinity responses to climate warming destabilise the Earth's thermostat. *Nat. Commun.* **14**, 1648 (2023).

4. Turner, J. P., Berry, T. W., Bowman, M. J. & Chapman, N. A. Role of the geosphere in deep nuclear waste disposal – an England and Wales perspective. *Earth Sci. Rev.* **242**, 104445 (2023).

5. Koppes, M. et al. Observed latitudinal variations in erosion as a function of glacier dynamics. *Nature* **526**, 100–103 (2015).

6. Alley, R. B., Cuffey, K. M. & Zoet, L. K. Glacial erosion: status and outlook. *Ann. Glaciol.* **60**, 1–13 (2019).

7. Cook, S. J., Swift, D. A., Kirkbride, M. P., Knight, P. G. & Waller, R. I. The empirical basis for modelling glacial erosion rates. *Nat. Commun.* **11**, 759 (2020).

8. Herman, F., De Doncker, F., Delaney, I., Prasicek, G. & Koppes, M. The impact of glaciers on mountain erosion. *Nat. Rev. Earth Environ.* **2**, 422–435 (2021).

9. Hallet, B., Hunter, L. & Bogen, J. Rates of erosion and sediment evacuation by glaciers: a review of field data and their implications. *Glob. Planet. Change* **12**, 213–235 (1996).

10. Fernandez, R. A. et al. Latitudinal variation in glacial erosion rates from Patagonia and the Antarctic Peninsula (46°S–65°S). *Geol. Soc. Am. Bull.* **128**, 1000–1023 (2016).

11. Herman, F. et al. Erosion by an Alpine glacier. *Science* **350**, 193–195 (2015).

12. Delmas, M., Calvet, M. & Gunnell, Y. Variability of Quaternary glacial erosion rates – a global perspective with special reference to the eastern Pyrenees. *Quat. Sci. Rev.* **28**, 484–498 (2009).

13. Patton, H. et al. The extreme yet transient nature of glacial erosion. *Nat. Commun.* **13**, 7377 (2022).

14. Fernandez, R. A., Anderson, J. B., Wellner, J. S. & Hallet, B. Timescale dependence of glacial erosion rates: a case study of Marinelli Glacier, Cordillera Darwin, southern Patagonia. *J. Geophys. Res.* **116**, F01020 (2011).

15. Humphrey, N. F. & Raymond, C. F. Hydrology, erosion and sediment production in a surging glacier: Variegated Glacier, Alaska, 1982–83. *J. Glaciol.* **40**, 539–552 (1994).

16. Harbor, J. M., Hallet, B. & Raymond, C. F. A numerical model of landform development by glacial erosion. *Nature* **333**, 347–349 (1988).

17. Hildes, D. H. D., Clarke, G. K. C., Flowers, G. E. & Marshall, S. J. Subglacial erosion and englacial sediment transport modelled for North American ice sheets. *Quat. Sci. Rev.* **23**, 409–430 (2004).

18. Jamieson, S. S. R., Hulton, N. R. J. & Hagdorn, M. Modelling landscape evolution under ice sheets. *Geomorphology* **97**, 91–108 (2008).

19. Herman, F., Beaud, F., Champagnac, J.-D., Lemieux, J.-M. & Sternai, P. Glacial hydrology and erosion patterns: a mechanism for carving glacial valleys. *Earth Planet. Sci. Lett.* **310**, 498–508 (2011).

20. Anderson, R. S., Molnar, P. & Kessler, M. A. Features of glacial valley profiles simply explained. *J. Geophys. Res.* **111**, F01004 (2006).

21. Herman, F. & Braun, J. Evolution of the glacial landscape of the Southern Alps of New Zealand: insights from a glacial erosion model. *J. Geophys. Res.* **113**, F02009 (2008).

22. Egholm, D. L., Pedersen, V. K., Knudsen, M. F. & Larsen, N. K. Coupling the flow of ice, water, and sediment in a glacial landscape evolution model. *Geomorphology* **141–142**, 47–66 (2012).
23. Lai, J. & Anders, A. M. Climatic controls on mountain glacier basal thermal regimes dictate spatial patterns of glacial erosion. *Earth Surf. Dyn.* **9**, 845–859 (2021).

24. Lai, J. & Anders, A. M. Tectonic controls on rates and spatial patterns of glacial erosion through geothermal heat flux. *Earth Planet. Sci. Lett.* **543**, 116348 (2020).

25. *Randolph Glacier Inventory - A Dataset of Global Glacier Outlines: Version 6.0* GLIMS Technical Report (RGI Consortium, 2017).

26. Zou, H. & Hastie, T. Regularization and variable selection via the elastic net. *J. R. Stat. Soc. B* **67**, 301–320 (2005).

27. Millan, R., Mouginot, J., Rabatel, A. & Morlighem, M. Ice velocity and thickness of the world's glaciers. *Nat. Geosci.* **15**, 124–129 (2022).

28. Davies, J. H. Global map of solid Earth surface heat flow. *Geochem. Geophys. Geosyst.* **14**, 4608–4622 (2013).

29. Karger, D. N. et al. Climatologies at high resolution for the Earth's land surface areas. *Sci. Data* **4**, 170122 (2017).

30. Brown, J. L., Hill, D. J., Dolan, A. M., Carnaval, A. C. & Haywood, A. M. PaleoClim, high spatial resolution paleoclimate surfaces for global land areas. *Sci. Data* **5**, 180254 (2018).

31. Giardini, D., Grünthal, G., Shedlock, K. M. & Zhang, P. The GSHAP Global Seismic Hazard Map. *Ann. Geophys.* **42**, 1225–1228 (1999).

32. Delaney, I. & Adhikari, S. Increased subglacial sediment discharge in a warming climate: consideration of ice dynamics, glacial erosion, and fluvial sediment transport. *Geophys. Res. Lett.* **47**, e2019GL085672 (2020).

33. Delaney, I. & Anderson, L. S. Debris cover limits subglacial erosion and promotes till accumulation. *Geophys. Res. Lett.* **49**, e2022GL099049 (2022).

34. Swift, D. A., Nienow, P. W. & Hoey, T. B. Basal sediment evacuation by subglacial meltwater: suspended sediment transport from Haut Glacier d'Arolla, Switzerland. *Earth Surf. Process.* **30**, 867–883 (2005).

35. Delaney, I., Werder, M. A. & Farinotti, D. A numerical model for fluvial transport of subglacial sediment. *J. Geophys. Res. Earth Surf.* **124**, 2197–2223 (2019).

36. Beaud, F., Flowers, G. E. & Venditti, J. G. Efficacy of bedrock erosion by subglacial water flow. *Earth Surf. Dyn.* **4**, 125–145 (2016).

37. Herman, F. & Brandon, M. Mid-latitude glacial erosion hotspot related to equatorial shifts in southern westerlies. *Geology* **43**, 987–990 (2015).

38. Portenga, E. W. & Bierman, P. R. Understanding Earth's eroding surface with ¹⁰Be. *GSA Today* **21**, 4–10 (2011).

39. Riebe, C. S., Kirchner, J. W., Granger, D. E. & Finkel, R. C. Strong tectonic and weak climatic control of long-term chemical weathering rates. *Geology* **29**, 511–514 (2001).

40. Molnar, P., Anderson, R. S. & Anderson, S. P. Tectonics, fracturing of rock, and erosion. *J. Geophys. Res.* **112**, F03014 (2007).

41. Dühnforth, M., Anderson, R. S., Ward, D. & Stock, G. M. Bedrock fracture control of glacial erosion processes and rates. *Geology* **38**, 423–426 (2010).

42. Hooyer, T. S., Cohen, D. & Iverson, N. R. Control of glacial quarrying by bedrock joints. *Geomorphology* **153–154**, 91–101 (2012).

43. Krabbendam, M. & Glasser, N. F. Glacial erosion and bedrock properties in NW Scotland: abrasion and plucking, hardness and joint spacing. *Geomorphology* **130**, 374–383 (2011).

44. Leith, K., Moore, J. R., Amann, F. & Loew, S. Subglacial extensional fracture development and implications for Alpine valley evolution. *J. Geophys. Res. Earth Surf.* **119**, 62–81 (2014).

45. Fahnestock, M., Abdalati, W., Joughin, I., Brozena, J. & Gogineni, P. High geothermal heat flow, Basal melt, and the origin of rapid ice flow in central Greenland. *Science* **294**, 2338–2342 (2001).

46. Syvitski, J. P. M., Peckham, S. D., Hilberman, R. & Mulder, T. Predicting the terrestrial flux of sediment to the global ocean: a planetary perspective. *Sediment. Geol.* **162**, 5–24 (2003).

47. Peucker-Ehrenbrink, B. Land2Sea database of river drainage basin sizes, annual water discharges, and suspended sediment fluxes. *Geochem. Geophys. Geosyst.* https://doi.org/10.1029/2008GC002356 (2009).

48. Salles, T. et al. Hundred million years of landscape dynamics from catchment to global scale. *Science* **379**, 918–923 (2023).

49. Koppes, M. N. & Hallet, B. Influence of rapid glacial retreat on the rate of erosion by tidewater glaciers. *Geology* **30**, 47–50 (2002).

50. Koppes, M. & Hallet, B. Erosion rates during rapid deglaciation in Icy Bay, Alaska. *J. Geophys. Res.* **111**, F02923 (2006).

51. Cuffey, K. M. et al. Entrainment at cold glacier beds. *Geology* **28**, 351–354 (2000).

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

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025
# Methods

## Synthesis of data

Glacial erosion rates (that is, data that aim to estimate the removal of material from the glacier bed) were synthesized from glacial studies published before January 2021 (Fig. 1 and Supplementary Table 1). Erosion rates from clearly defined glacial catchments for which glaciological, geological and topoclimatic variables are available were collated. We include a considerable number of datapoints, previously synthesized in erosion rate datasets and review studies (for example, refs. 7,9,12,52), which are cited in Supplementary Table 1. Like previous compilations, we synthesized erosion data from studies that utilized four different methods to measure glacial erosion rate, including sediment discharge estimates from meltwater streams, bulk sediment in proglacial environments, the depth of overlying bedrock removal measured using geochronological techniques (terrestrial cosmogenic nuclide dating) and instrumental measurements from extant glaciers (for example, ref. 12). In cases where several studies measured glacier erosion rates for the same glacier, we included all measurements. In line with previous syntheses of glacial erosion rates used to derive erosion rules³, we synthesized short-term erosion rates (164 contemporary erosion rates and 17 spanning the early to late Holocene) derived only from topographically constrained extant glaciers and ice caps (separated into glacier basins as per ref. 25) characterized by dominantly warm basal thermal regimes. Holocene rates included are those used in previous erosion rules or erosion rate compilations.

## Glaciological variables

Where not measured by the original erosion rate study, key glaciological variables including ‘maximum length of an ice mass,’ ‘area of an ice mass,’ ‘median elevation of an ice mass,’ ‘ice surface slope,’ ‘termination type,’ ‘form’ and ‘surge/non surge behaviour,’ were obtained from the RGI²⁵. Where the original author and the RGI had missing variables, data were sourced from the Shuttle Radar Topography Mission digital elevation model (DEM) between 56°S and 60°N (-80% of the globe; https://www2.jpl.nasa.gov/srtm/) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global DEM version 2 (https://asterweb.jpl.nasa.gov/gdem.asp) for all other areas. In cases where the maximum length of an ice mass (longest surface flowline) was not recorded by the RGI²⁵, we used the same algorithm to calculate this variable³². Likewise, ‘median elevation of an ice mass’ was derived for all glaciated DEM cells within their area polygon as in the RGI²⁵.

Following the practices used by previous ‘velocity-erosion rule’ studies⁵,⁷,¹¹,¹⁵, when available, we use measurements of maximum and median sliding and surface velocities and median ice thicknesses (94 of 181 cases for velocity and 20 of 181 cases for ice thickness) that have been measured over multiple years and broadly correspond to the timespan over which erosion rate measurements were made. These were either collected by the authors of the original erosion rate study or another study that took measurements that temporally correspond with measurements taken by the original erosion rate study. For some glaciers, the same value is used for median and maximum surface velocity due to a lack of available measurements. Consistent with ref. 7, we tried to avoid very short-term measurements of glacial velocity, such as those recorded over a few days, because these have been shown to vary greatly⁵⁴; however, in some cases, short-term measurements were included as these were the only velocity measurements available. In some land- and marine-terminating cases, a slight offset exists between glacial velocity and an erosion rate measurement. Additionally, where studies average erosion over longer timescales (for example, >100 years), no velocity measurements span this entire erosion rate interval. Although we have used, to our knowledge, velocity data that are the best available temporal match to erosion rates, we acknowledge these limitations and highlight that mismatches in the interval over which erosion rates and velocities were measured still occur. We also note that surging activity is not well documented for all

erosion rates, and we acknowledge that some surge-type erosion rate periods and velocity measurement periods encompass quiescence and active surge conditions.

In cases where no direct measurement of surface or sliding velocity was recorded, these data were derived from ref. 27. This dataset has a resolution of 50 m, and surface velocity is derived using cross-correlation⁵⁵ between consecutive images taken in 2017 and 2018. Ice thickness was estimated from previously established ice motion and surface slope data²⁷ using the shallow ice approximation and calibrating parameters using field-based measurements³⁶. We extracted median and maximum values of surface velocity and median ice thickness for each glacier in the RGI²⁵. These values are summarized in Supplementary Table 1 and used to estimate the near-global distribution of glacial erosion rates (Supplementary Table 2 and Fig. 4). In cases for which ice thickness and surface velocity variables were not documented corresponding with the period of erosion rate data collection or by ref. 27, we used surface velocities from regional ice velocity datasets⁵⁷⁻⁵⁹ and a global modelled ice thickness dataset⁶⁰.

In all cases where sliding velocity ($U_{sliding}$) is not directly measured, we estimated sliding velocities based on surface speed ($U_{surface}$) and slope ($\beta$)²⁷: $U_{sliding} = \beta U_{surface}$. In this parameterization, we first computed the ratio between surface slopes and surface speed. This ratio was then constrained within the interval (0.001, 0.1). The sliding speeds in regions where this ratio is 0.001 were assumed to be $U_{sliding} = 0.9 \times U_{surface}$, and sliding speeds were set to $U_{sliding} = 0.1 \times U_{surface}$ in regions where this ratio is 0.01. We use a linear interpolation for $\beta$ in between these two regimes.

## Geological variables

We extracted lithologic information from the Global Lithological Map (GLiM)⁶¹ to describe the bed of each glacier from its centroid (glacier area from ref. 25) using the ‘Calculate geometry (centroid)’ tool in ArcGIS 10.7. GLiM data were rasterized at a resolution of 50 m. Within the GLiM polygons labelled as ‘ice and glaciers’, ‘no data’ and ‘water bodies’ (probably proglacial lakes in the context of this study), a lithology value was set based on a nearest neighbour value using the ‘Nibble (spatial analyst)’ tool in ArcGIS 10.7. The lithology specified by the original study was used (if recorded) when erosion rates were recorded from terrestrial cosmogenic nuclide or in situ instrumental studies. Lithologies were then further grouped into three categories: ‘crystal-line’ (inclusive of intrusive igneous and metamorphic), ‘sedimentary’ (inclusive of carbonates, sandstones and shales) and ‘volcanics’ (inclusive of volcanics and basalts).

Following investigations focused on fluvial erosion³⁸ we used PGA (described as the maximum instantaneous ground movement during an earthquake, based on the amplitude of the largest absolute acceleration with a 10% likelihood of exceeding in 50 years³¹,³⁸) as an alternative means to estimate seismicity. We extracted PGA values from the Global Seismic Hazard Assessment Program map³¹. This dataset has a resolution of 10 km. Using this dataset, we assigned locations in Svalbard, Antarctica (omitted from ref. 31) and PGA values <1.7 m s⁻² (a natural break in the global PGA dataset) as seismically inactive. In contrast, all other sites (that is, PGA ≥ 1.7 m s⁻²) were considered active (Fig. 3d).

We assigned an average GHF value to each glacier using a global map of solid earth surface heat flow²⁸. This map has a resolution of 2°, which may limit the differentiation of GHF under neighbouring glaciers where GHF is spatially highly variable, such as in volcanic regions. We extracted median values from glacier centroids (glacier area from ref. 25) using the ‘Calculate geometry (centroid)’ tool in ArcGIS 10.7 and collating the ‘Final_HFafive’ attribute of ref. 28.

## Topoclimatic variables

Glaciers were divided into climate zones according to the Köppen-Geiger classification map⁶². However, these values were excluded from predictor analysis because they were very unequal (>80% were
polar). MAAT and MAP values, which broadly correspond with glacial erosion rates that were on contemporary timescales, were obtained for each glacial erosion rate site from the CHELSA V1.2 climate dataset (http://chelsa-climate.org). These data span the period 1979–2013²⁹ and have a resolution of ~1 km (30 arcsec). MAAT and MAP values for erosion rates that occurred during the late Holocene (4.2–0.3 kyr ago), mid-Holocene (8.326–4.2 kyr ago) and early Holocene (11.7–8.326 kyr ago) were obtained for each erosion rate site from the PaleoClim dataset (http://www.paleoclim.org)³⁰. These datasets have a resolution of ~5 km (2.5 arcsec). For each of these datasets, variables were obtained from glacier centroids (glacier area from ref. 25) using the ‘Calculate geometry (centroid)’ tool in ArcGIS 10.7.

## Dataset and analysis limitations

Glacial erosion rates in our synthesis were collected using four different measurement techniques, and as such, biases may be present due to how each technique constrains an erosion rate (see also refs. 8,12,63). For example, erosion rates calculated using instrumental measurements often only record glacial erosion resulting from abrasion and do not consider the effects of quarrying or subglacial meltwater incision¹². Furthermore, erosion rates derived from bulk sediment in glacial environments can be affected by short-term sediment storage in proglacial lakes or subglacial deposits⁹,¹⁰,¹²,⁶⁴,⁶⁵. Likewise, meltwater stream flux-derived erosion rates are often solely calculated based on suspended sediment loads and do not take into account or measure solute and bedload¹². To minimize method-based biases, we collate all rates where more than one erosion rate measurement exists (collected via different techniques), to reflect the spectrum of potential erosion rates of one glacier.

The distribution of collated erosion rates separated by measurement technique is shown in Supplementary Fig. 1. Although the distribution of erosion rates shown in Supplementary Fig. 1 are from various regions and glacial systems, the median erosion rates separated by measurement technique are within a single order of magnitude. This broad agreement is suggestive that no single measurement technique exerts a strong systematic bias. We explored the measurement technique-specific relationships for each explanatory variable in each of three glacier types investigated (Supplementary Fig. 2). While only a few measurements were collected via instrumental means and terrestrial cosmogenic nuclide dating, making comparison difficult, the trends they exhibit are broadly consistent and variability between methods is comparable to the natural or measurement variability between observations.

The duration over which a glacial erosion rate is recorded has been shown to affect the rate calculated, with erosion rates measured over longer durations often showing lower average values²,⁷,¹²,¹⁴,⁴⁹,⁵⁰,⁶³,⁶⁶,⁶⁷. To minimize this timescale dependence¹⁴ (also observed in a range of geomorphic systems, for example, the Sadler effect⁶³,⁶⁸) we limit our erosion rate synthesis to shorter-term measurements of glacial erosion (164 erosion rates that are averaged over the past ~350 years or less and a further 17 that span the early to late Holocene), rather than including erosion rates that span multiple glaciations. A comparison between our synthesized glacial erosion measurements and the timespan over which they were measured indicates an insignificant correlation ($y = 0.87 \times x^{-0.056}$, $R^2 = 0.0097$; Supplementary Fig. 3), suggesting the timespan over which measurements were recorded imparts a negligible influence on the erosion rates in our dataset.

Temporal discrepancies between when variables are measured and when an erosion rate is measured could affect our results. This is especially important in the case of velocity, which has been shown to vary dramatically on short timescales⁶⁹. This effect may thus affect the relationship observed between a glacier's velocity and the erosion rate measurement⁷,¹⁴,⁶⁷ in this study. To minimize this effect, when available, we utilized measurements of velocity that temporally overlapped with when an erosion rate measurement was made. In addition, to explore

the potential impact of including velocity data from various sources and temporal scopes, we conducted a sensitivity analysis of three data synthesis strategies: (1) favouring remotely sensed data (which often do not temporally overlap with an erosion rate measurement); (2) favouring velocity measurements that temporally overlap with erosion rate measurements; and (3) an approach that uses all available data, but with a preference towards temporally matched author-reported velocities where available. We compared the LOOCV performance of elastic net models given only velocity data as predictors but including all its permutations (that is, median, maximum, surface, sliding, linear velocity ($m \text{ yr}^{-1}$) and their natural logarithms; Supplementary Table 4). We observed that using the full suite of available velocity data sources, with preference to originally reported velocities, provides generally lower RMSE and MAE estimates than other methods, particularly in land-terminating and surge-type environments (Supplementary Tables 3 and 4), supporting our adoption of this data treatment approach.

Less than 15% of erosion rates in our synthesis have been recorded from surge-type glaciers (in comparison with <0.6% of the RGI²⁵,⁷⁰). Erosion rates measured during periods of surging have been suggested to be higher than those under more quiescent conditions¹⁵,⁷¹. This heterogeneity is the motivation behind our stratification of surge-type glaciers in our analyses. We note that in our analysis we attribute RGI categories of ‘possible, probable and observed surge-type glaciers’ as ‘surge-type’, but acknowledge that these glaciers will not always be actively surging.

As stated in the main text, to our knowledge, constraints on erosion from beneath cold-based glaciers that exhibit a minimal sliding velocity have only been recorded by ref. 51 (sliding velocity <0.01 m yr⁻¹). As a result, we cannot statistically assess the controls of glacial erosion in such environments, and we choose to remove this single datapoint and instead focus on glaciers with predominantly warm basal thermal regimes. Our predictive equations are, therefore, not built or tested using data from such extreme cold-based environments and are not suitable for use in these environments. Any predictions made using these equations for cold-based environments are maximum limiting only.

We would have preferred to assess bed roughness and fracture density when developing predictive equations between glaciological, topoclimatic and geological variables and glacial erosion rates. However, we do not include this variable, as, to our knowledge, it has not been estimated globally or empirically determined in enough extant glacial environments to include in our study.

## Statistical analysis

We employed elastic net regression²⁶, a regularized form of linear regression that combines L1 and L2 penalties (as in, Lasso regression⁷² and Ridge regression⁷³, respectively) to predict erosion rate across all cases with complete explanatory covariates. To study environment-specific controls on erosion, we subdivided glaciers into three categories: surge-type, land-terminating (non-surging) and marine-terminating (non-surging). Tuning parameters for the various models were optimized by grid search and exhaustive LOOCV. LOOCV was also used to evaluate predictive model performance once tuned. We use the same metrics to compare both our primary (all-predictor) models and those that explored surface- and sliding-velocity-only predictors ($R^2$ and RMSE for goodness of fit to the natural logarithm of erosion rate, and MAE to assess average residuals in linear erosion rates).

All continuous covariates were input as candidate explanatory variables in two forms: their conventional (that is, linear) representation and the natural logarithm of these values (that is, log-transformed). Erosion rate, the dependent variable, was represented by the natural logarithm of erosion rate. This follows previous studies⁵,⁷¹⁰ and addresses the skewed nature of global erosion rate data over many orders of magnitude, while simultaneously promoting linearity,
a requirement of elastic net regression. The categorical variables, lithology and PGA class, were encoded to preserve each level (that is, non-full-rank parameterization).

Within our synthesis, several datapoints were excluded from statistical analysis based on collection and calculation techniques or glaciological context (for reference). A description outlining the reasoning for the exclusion of an erosion rate is included in Supplementary Table 1. From a statistical perspective, these data were flagged as high-residual points when included in an initial linear model or were otherwise problematic.

We present variable importance as an indication of variable effect in our predictive models (that is, the magnitude of the absolute value of regression coefficients after values for each covariate are converted to z-scores). Variable importance calculation was bootstrapped with 1,000 resamples to assess the stability of effect strength. Variable importance is portrayed as the percent of the maximum resample (as in Fig. 3) or the proportion of total predictive influence (that is, one variable type relative to the sum of all variable importance). The coefficients of non-z-scored variables are relayed in formulas (Supplementary Text 2) and in tabular form (Supplementary Table 6).

Many of the explanatory variables used in this study were highly correlated (Fig. 2). These dependencies include fundamental correlations between linear and log-transformed versions of the same variables. Although our approach is strictly empirical in its perspective, the effects of interrelated and correlated variables must still be evaluated from a mechanistic view. Likewise, those potential drivers with a well-founded physical basis should not be discounted based on their exclusion from, or substantial shrinkage in, our models.

## Erosion rate predictions

Prediction intervals accompanying each glacier's erosion rate were calculated through bootstrap resampling with 2,000 resamples. We employ the 0.632+ approach⁷⁴ to estimate and account for sample noise and bias in our predictions. Within this framework, each bootstrapped prediction is the sum of the model's base prediction and a residual term (which is a weighted combination of training data and out-of-bag residuals, where the weight depends on the relative overfitting rate and the proportion of held-out data for each bootstrap).

Owing to incomplete surface velocity and ice thickness data necessary to predict glacial erosion rates, all glaciers in Antarctica (part of RGI region 19) are excluded from our predictions. We restrict our reported predictions to the range of erosion rates of our training dataset (0–500 mm yr⁻¹). In cases where erosion rates are higher than this range, we do not report erosion rates, but instead, we assign a value of >500 mm yr⁻¹. The total estimated volumetric erosion for the RGI glaciers is the estimated erosion rate for each glacier multiplied by the glacier's area, while censoring exceptionally high erosion rates to 500 mm yr⁻¹. Considering the spatial heterogeneity implicit in glaciized catchment erosion rates (that is, the central point estimate could be meaningfully higher or lower than the across-the-glacier average), including confidence intervals on this value would imply unfounded precision in the global sum. However, we explore the full distribution of bootstrap-derived mass flux (after censoring) in Supplementary Fig. 5.

## Data availability

All data needed to evaluate the conclusions in this paper are present in the paper and/or Supplementary Information. In some available cases, glacier thickness and velocity data products were derived from https://doi.org/10.1038/s41561-021-00885-z (ref. 27); where these data were absent, glacier thickness and velocity data were derived from the ITS_LIVE project (https://its-live.jpl.nasa.gov/), NVE project (https://hdl.handle.net/11250/2828417), https://doi.org/10.5194/tc-12-521-2018 (ref. 58) and https://doi.org/10.3929/ethz-b-000315707. Glacier length, area, elevation, ice surface slope, termination type, form and

surge characteristics were derived from the GLIMS/RGI (https://doi.org/10.7265/N5-RGI-60) and in missing cases was generated using a combination of data from the Shuttle Radar Topography Mission DEM (https://science.jpl.nasa.gov/projects/SRTM/) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global DEM version 2 (https://asterweb.jpl.nasa.gov/gdem.asp). PGA data are available via the Global Seismic Hazard Assessment Program (http://static.seismo.ethz.ch/GSHAP/index.html). GHF was derived from a global map of solid earth surface heat flow (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/ggge.20271)²⁸. Climate zone categories are available as the Köppen–Geiger climate classification map (http://www.gloh2o.org/koppen/). Contemporary and palaeo MAAT and MAP products are available via the CHELSA V1.2 climate dataset (http://chelsa-climate.org) and Paleoclim dataset (http://www.paleoclim.org), respectively.

## Code availability

The R code used to train the machine learning models and the trained R models are available via Zenodo at https://zenodo.org/records/15832805 (ref. 75).

## References

52. Gurnell, A., Hannah, D. & Lawler, D. Suspended sediment yield from glacier basins. *IAHS-AISH Publ.* **236**, 97–104 (1996).

53. Machguth, H. & Huss, M. The length of the world's glaciers – a new approach for the global calculation of center lines. *Cryosphere* **8**, 1741–1755 (2014).

54. Purdie, H. L., Brook, M. S. & Fuller, I. C. Seasonal variation in ablation and surface velocity on a temperate maritime glacier: Fox Glacier, New Zealand. *Arct. Antarct. Alp. Res.* **40**, 140–147 (2008).

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

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

57. Nagy, T. & Andreassen, L. M. *Glacier Surface Velocity Mapping with Sentinel-2 Imagery in Norway Technical Report No. 37* (Norwegian Water Resources and Energy Directorate, 2019).

58. Gardner, A. S. et al. Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years. *Cryosphere* **12**, 521–547 (2018).

59. Gardner, A. S., Fahnestock, M. A. & Scambos, T. A. *MEaSUREs ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities, Version 1* (NASA National Snow and Ice Data Center Distributed Active Archive Center, 2019); https://doi.org/10.5067/6II6VW8LLWJ7

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

61. Hartmann, J. & Moosdorf, N. The new global lithological map database GLiM: a representation of rock properties at the Earth surface. *Geochem. Geophys. Geosyst.* **13**, Q12004 (2012).

62. Beck, H. E. et al. Present and future Köppen–Geiger climate classification maps at 1-km resolution. *Sci. Data* **5**, 180214 (2018).

63. Turowski, J. M. & Cook, K. L. Field techniques for measuring bedrock erosion and denudation. *Earth Surf. Process. Landf.* **42**, 109–127 (2017).

64. Hasholt, B., van As, D., Mikkelsen, A. B., Mernild, S. H. & Yde, J. C. Observed sediment and solute transport from the Kangerlussuaq sector of the Greenland Ice Sheet (2006–2016). *Arct. Antarct. Alp. Res.* **50**, S100009 (2018).

65. Mancini, D. et al. Filtering of the signal of sediment export from a glacier by its proglacial forefield. *Geophys. Res. Lett.* **50**, (2023).
66. Elverhøi, A., Hooke, R. L. & Solheim, A. Late Cenozoic erosion and sediment yield from the Svalbard-Barents Sea region: implications for understanding erosion of glacierized basins. *Quat. Sci. Rev.* **17**, 209–241 (1998).

67. Ganti, V. et al. Time scale bias in erosion rates of glaciated landscapes. *Sci. Adv.* **2**, e1600204 (2016).

68. Sadler, P. M. Sediment accumulation rates and the completeness of stratigraphic sections. *J. Geol.* **89**, 569–584 (1981).

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

70. Sevestre, H. & Benn, D. I. Climatic and geometric controls on the global distribution of surge-type glaciers: implications for a unifying model of surging. *J. Glaciol.* **61**, 646–662 (2015).

71. Knudsen, N. T., Yde, J. C. & Gasser, G. Suspended sediment transport in glacial meltwater during the initial quiescent phase after a major surge event at Kuannersuit Glacier, Greenland. *Geogr. Tidsskr.* **107**, 1–7 (2007).

72. Tibshirani, R. Regression shrinkage and selection via the lasso: a retrospective. *J. R. Stat. Soc. B* **73**, 273–282 (2011).

73. Hoerl, A. E. & Kennard, R. W. Ridge regression: biased estimation for nonorthogonal problems. *Technometrics* **42**, 80–86 (2000).

74. Efron, B. & Tibshirani, R. Improvements on cross-validation: the 632+ bootstrap method. *J. Am. Stat. Assoc.* **92**, 548–560 (1997).

75. Norris, S., Gosse, J. & Bolton, M. Models and code for estimating glacier erosion rate. Zenodo https://doi.org/10.5281/zenodo.12774358 (2025).

## Acknowledgements

This work was undertaken in partnership with and financially supported by the Canadian Nuclear Waste Management Organization (NWMO). J.C.G. and S.L.N. acknowledge funding from the Natural Sciences and Engineering Research Council of Canada (NSERC-DG-06785-19 and NSERC-DG-04031-23, respectively).

R.M., J.M. and A.R. acknowledge support from the French space agency (CNES) and the LabEx OSUG@2020 (Investissements d'avenir—ANR10 LABX56) in establishing worldwide glacier velocities and thicknesses. R.B.A. acknowledges funding from NASA (2021-80NSSC22K0384), NSF (1738934) and the Heising-Simons Foundation (2018-0769). The authors thank J. Fast and D. Evans for assistance with our erosion rate synthesis.

## Author contributions

S.L.N. and J.C.G. designed the project and compiled the erosion rate dataset with input from R.B.A. M.S.M.B. completed statistical analysis. R.M., J.M., A.R. and M.M. conducted surface velocity and ice thickness analysis. S.L.N. analysed the data and wrote the manuscript with input from all co-authors. All authors contributed to address reviewer comments.

## 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-025-01747-8.

**Correspondence and requests for materials** should be addressed to Sophie L. Norris or John C. Gosse.

**Peer review information** *Nature Geoscience* thanks Małgorzata Chmiel, Simon Cook and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Tom Richardson, in collaboration with the *Nature Geoscience* team.

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