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

# The unquantified mass loss of Northern Hemisphere marine-terminating glaciers from 2000–2020

Received: 5 January 2022

Accepted: 7 September 2022

Published online: 11 October 2022

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

Check for updates

William Kochtitzky<fnref n="1" />✉, Luke Copland<fnref n="1" />, Wesley Van Wychen<fnref n="1" />,², Romain Hugonnet³,<fnref n="4" />,⁵, Regine Hock⁶,⁷, Julian A. Dowdeswell⁸, Toby Benham⁹, Tazio Strozzi⁹, Andrey Glazovsky<fnref n="10" />, Ivan Lavrentiev<fnref n="10" />, David R. Rounce¹¹, Romain Millan¹², Alison Cook<fnref n="1" />, Abigail Dalton<fnref n="1" />, Hester Jiskoot¹³, Jade Cooley¹³, Jacek Jania¹⁴ & Francisco Navarro¹⁵

In the Northern Hemisphere, ~1500 glaciers, accounting for 28% of glacierized area outside the Greenland Ice Sheet, terminate in the ocean. Glacier mass loss at their ice-ocean interface, known as frontal ablation, has not yet been comprehensively quantified. Here, we estimate decadal frontal ablation from measurements of ice discharge and terminus position change from 2000 to 2020. We bias-correct and cross-validate estimates and uncertainties using independent sources. Frontal ablation of marine-terminating glaciers contributed an average of 44.47 ± 6.23 Gt a⁻<fnref n="1" /> of ice to the ocean from 2000 to 2010, and 51.98 ± 4.62 Gt a⁻<fnref n="1" /> from 2010 to 2020. Ice discharge from 2000 to 2020 was equivalent to 2.10 ± 0.22 mm of sea-level rise and comprised approximately 79% of frontal ablation, with the remainder from terminus retreat. Near-coastal areas most impacted include Austfonna, Svalbard, and central Severnaya Zemlya, the Russian Arctic, and a few Alaskan fjords.

When glaciers terminate in the ocean, mass is lost by frontal ablation at the ice-ocean interface, including mass loss due to iceberg calving (i.e., the mechanical release of icebergs or smaller pieces of ice to the ocean), submarine frontal melting, and subaerial frontal melting and sublimation<fnref n="1" />. Quantifying frontal ablation is essential to partitioning glacier mass loss components, which informs sea-level budgets, estimates of glacier freshwater contributions to the ocean, impacts on marine ecosystems and iceberg hazards. Recent model results suggest that Northern Hemisphere frontal ablation was around 39 Gt a⁻<fnref n="1" /> from

1980 to 1999². Frontal ablation accounted for 10% of the 882 Gt a⁻<fnref n="1" /> global glacier ablation rate (which includes surface melt), although about half of this occurred in the Southern Hemisphere². The same modeling study projected a Northern Hemisphere frontal ablation rate of 50.6 ± 23.8 Gt a⁻<fnref n="1" /> for 2020–2040². However, observation-based estimates of frontal ablation have been limited, partially because many studies only quantify ice discharge³,<fnref n="4" /> as they do not include mass changes due to terminus retreat or advance, or these are only computed for individual glaciers⁵. Alaska⁶ and Svalbard⁷ are the only

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

<footnote n="1">¹Department of Geography, Environment and Geomatics, University of Ottawa, Ottawa, Ontario K1N 6N5, Canada. ²Department of Geography and Environmental Management, University of Waterloo, Waterloo, Ontario, Canada. ³LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, Toulouse, France.</footnote>

<footnote n="4">⁴Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zürich, Zürich, Switzerland. ⁵Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), Birmensdorf, Switzerland. ⁶Department of Geosciences, University of Oslo, Oslo, Norway. ⁷Geophysical Institute, University of Alaska Fairbanks, Fairbanks, USA. ⁸Scott Polar Research Institute, University of Cambridge, Cambridge, UK. ⁹Gamma Remote Sensing, Gümligen (BE), Switzerland.</footnote>

<footnote n="10">¹⁰Institute of Geography, Russian Academy of Sciences, Moscow, Russia. ¹¹Department of Civil and Environmental Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA. ¹²Institut des Géosciences de l'Environnement, CNES, Grenoble, France. ¹³Department of Geography & Environment, University of Lethbridge, Lethbridge, Alberta, Canada. ¹⁴University of Silesia, Katowice, Poland. ¹⁵Universidad Politécnica de Madrid, Madrid, Spain.</footnote>

<footnote>e-mail: will.kochtitzky@uottawa.ca</footnote>
regions in the Northern Hemisphere with complete observation-based frontal ablation estimates, but these have not been updated since 2013.

Whereas many regional-scale studies have derived regional total glacier mass changes from glaciological and geodetic methods<sup>8–11</sup>, satellite gravimetry<sup>12,13</sup>, and numerical modeling<sup>14–16</sup>, no study has determined the frontal ablation rate for the entire Northern Hemisphere. This has mainly been due to a previous lack of spatially and temporally consistent glacier ice thickness, surface velocity, and terminus position observations, from either field or remote sensing surveys, which is exacerbated at high latitudes.

Here we identify every marine-terminating glacier in the Northern Hemisphere, distinct from the Greenland Ice Sheet, and use measurements or estimates of ice thickness, surface velocity, and terminus position changes, together with a climatic mass balance model, to estimate mean frontal ablation (see “Methods” for details; Fig. S1). We split the observations into two periods: 2000 to 2010 and 2010 to 2020. The former decade is based on terminus position observations from summer 2000 and 2010, and velocity observations from 2000 to 2009. The latter decade is based on terminus position observations from summer 2010 and 2020, and velocity observations from 2010 to 2019. Because directly measuring frontal ablation components (i.e., calving, sublimation, and melting) is difficult, we estimate frontal ablation indirectly, by calculating ice discharge (i.e., mass flux through a flux gate perpendicular to glacier flow upstream of the terminus) and mass change due to retreat or advance (henceforth, terminus mass change), while accounting for the climatic mass balance (i.e., the mass changes due to snow accumulation, surface melt and refreezing) below the flux gate. Our computation of frontal ablation excludes submarine melting below floating tongues, but this is likely negligible as virtually no floating glacier tongues remain in the Arctic<sup>17,18</sup>. The residual examples are primarily within Franz Josef Land, Russia<sup>17</sup>.

## Results and discussion

### Complete hemispheric estimates of discharge and terminus change

Guided by the Randolph Glacier Inventory version 6.0 (RGI6)<sup>10,19</sup>, we identified 1496 marine-terminating glaciers and manually mapped their terminus positions in summer 2000, 2010, and 2020 from satellite imagery (see “Methods”)²⁰. The Greenland Periphery has the most marine-terminating glaciers (537), followed by the Russian Arctic (412), Arctic Canada North (252), Svalbard and Jan Mayen (160 + 6 = 166; hereafter Svalbard), Arctic Canada South (86), Alaska (42), and Iceland (1). However, the length of terminus flux gates in Russia is more than two times that of Greenland Periphery glaciers

(Table 1). North Asia has one small marine-terminating glacier on one of the De Long Islands, and we therefore include it in the Russian Arctic RGI6 region. Between 2000 and 2010, 49 of the marine-terminating glaciers lost contact with the ocean and became land-terminating with another 120 glaciers following in the next decade.

The frontal ablation rate from all 1496 glaciers was 44.47 ± 6.23 Gt a⁻¹ (all uncertainties given at the 95% confidence level, i.e., two standard deviations) between 2000 and 2010, and 51.98 ± 4.62 Gt a⁻¹ between 2010 and 2020 (Table 1). While this suggests a slight increase in frontal ablation between the two decades, the estimates are within the uncertainties of each other. Ice discharge accounted for 80 and 78% of the frontal ablation for 2000–2010 and 2010–2020, respectively, with the remaining losses due to terminus retreat. Rates of frontal ablation, ice discharge, and terminus mass loss did not vary significantly (i.e., non-overlapping uncertainties) between 2000–2010 and 2010–2020, with few exceptions (Table 1). In Alaska, terminus mass loss increased nearly tenfold while ice discharge decreased by about 15%. In Svalbard, ice discharge nearly tripled, primarily due to the surge of Basin-3 of Austfonna Ice Cap. Owing to improved data quality during 2010–2020 and smaller uncertainties, we focus primarily on estimates of this decade, unless otherwise specified.

Among the seven glacierized regions studied, the Russian Arctic experienced the highest frontal ablation rate, followed by Svalbard and Alaska (Table 1; Fig. 1). Greenland Periphery and Arctic Canada North had similar frontal ablation rates, while Iceland and Arctic Canada South were markedly lower (Table 1; Fig. 1). From 2000–2010 Hubbard Glacier in Alaska showed the highest frontal ablation rate (4.13 ± 0.05 Gt a⁻¹) for a single glacier, but was surpassed by Basin-3 of Austfonna in Svalbard (6.10 ± 0.24 Gt a⁻¹) during the next decade (Fig. 2).

Even though Arctic Canada North has 40% of the total Northern Hemisphere area of marine-terminating glaciers, Alaska, Svalbard, and the Russian Arctic each produce more frontal ablation. The Russian Arctic contains 22% of the marine-terminating glacier area, but it accounts for 33% of total frontal ablation.

### Uneven distribution of frontal ablation across all regions

Within each region, a few large glaciers have a disproportionate impact on the total frontal ablation over the study period (Fig. 2). In Svalbard and the Russian Arctic, 27 and 32% of glaciers (i.e., 61 and 80% of regional marine-terminating glacier area), respectively, make up 90% of the regional frontal ablation rate. A similar, small number of glaciers make up about 90% of frontal ablation in other regions as well, including 24% of glaciers in Arctic Canada North, 14% in Alaska, 28% in Greenland Periphery, and 28% in Arctic Canada South. In total, merely

<table><caption>**Table 1 | Frontal ablation and components**</caption><thead><tr><th rowspan="3"></th><th colspan="2">Frontal ablation</th><th colspan="2">Ice discharge</th><th colspan="2">Terminus mass loss</th><th rowspan="2">Length of flux gates</th><th rowspan="2">Number of glaciers</th></tr><tr><th colspan="2">(Gt a⁻¹)</th><th colspan="2">(Gt a⁻¹)</th><th colspan="2">(Gt a⁻¹)</th></tr><tr><th>2000–2010</th><th>2010–2020</th><th>2000–2010</th><th>2010–2020</th><th>2000–2010</th><th>2010–2020</th><th>(km)</th><th></th></tr></thead><tbody><tr><td>Alaska</td><td>11.59 ± 0.39</td><td>10.68 ± 0.33</td><td>11.49 ± 0.35</td><td>9.79 ± 0.18</td><td>-0.1 ± 0.17</td><td>-0.89 ± 0.28</td><td>80</td><td>42</td></tr><tr><td>Arctic Canada North</td><td>4.14 ± 1.11</td><td>4.28 ± 1.18</td><td>2.68 ± 0.65</td><td>2.24 ± 0.33</td><td>-1.45 ± 0.91</td><td>-2.03 ± 1.14</td><td>481</td><td>252</td></tr><tr><td>Arctic Canada South</td><td>0.09 ± 0.08</td><td>0.09 ± 0.08</td><td>0.04 ± 0.03</td><td>0.03 ± 0.02</td><td>-0.05 ± 0.07</td><td>-0.06 ± 0.07</td><td>40</td><td>86</td></tr><tr><td>Greenland periphery</td><td>4.31 ± 1.57</td><td>3.18 ± 1.09</td><td>2.25 ± 0.77</td><td>1.88 ± 0.43</td><td>-2.05 ± 1.36</td><td>-1.29 ± 1.0</td><td>697</td><td>537</td></tr><tr><td>Iceland</td><td>0.10 ± 0.10</td><td>0.03 ± 0.03</td><td>0.10 ± 0.10</td><td>0.001 ± 0.07</td><td>-0.003 ± 0.02</td><td>-0.03 ± 0.05</td><td>5</td><td>1</td></tr><tr><td>Svalbard and Jan Mayen</td><td>7.62 ± 2.65</td><td>16.82 ± 2.48</td><td>4.88 ± 1.98</td><td>14.41 ± 1.05</td><td>-2.74 ± 1.77</td><td>-2.4 ± 2.25</td><td>657</td><td>166</td></tr><tr><td>Russia: Franz Josef Land</td><td>10.46 ± 5.05</td><td>7.44 ± 3.29</td><td>8.78 ± 4.73</td><td>4.68 ± 2.21</td><td>-1.68 ± 1.76</td><td>-2.76 ± 2.44</td><td>1377</td><td>328</td></tr><tr><td>Russia: Novaya Zemlya</td><td>2.67 ± 1.02</td><td>4.15 ± 0.93</td><td>2.02 ± 0.84</td><td>3.14 ± 0.37</td><td>-0.65 ± 0.57</td><td>-1.0 ± 0.86</td><td>182</td><td>39</td></tr><tr><td>Russia: Severnaya Zemlya</td><td>3.50 ± 1.14</td><td>5.33 ± 0.88</td><td>3.15 ± 1.06</td><td>4.23 ± 0.59</td><td>-0.36 ± 0.43</td><td>-1.1 ± 0.64</td><td>283</td><td>45</td></tr><tr><td>Total</td><td>44.47 ± 6.23</td><td>51.98 ± 4.62</td><td>35.40 ± 5.41</td><td>40.41 ± 2.61</td><td>-9.07 ± 3.07</td><td>-11.57 ± 3.81</td><td>3802</td><td>1496</td></tr></tbody></table>

Decadal mean mass losses (± uncertainties) by region (from west to east) for all Northern Hemisphere marine-terminating glaciers, due to advance or retreat of terminus position (terminus mass change), ice flow through a terminus flux gate (ice discharge), and the sum of these (frontal ablation).
<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Regional frontal ablation of marine-terminating glaciers.** Bar charts show the decadal-mean frontal ablation rates by Northern Hemisphere region from 2000 to 2020 and are divided into ice discharge and terminus mass change (uncertainties</figcaption></figure>

are black lines). Pie charts show the proportion of marine vs. land-terminating glaciers for each region. The Northern Hemisphere glacier proportions are shown in the legend.

1% (15 out of 1496) of marine-terminating glaciers account for 45% of Northern Hemisphere frontal ablation, and 2% (30 glaciers) account for 55% of frontal ablation.

Most marine-terminating glaciers (87% of the 1496 glaciers) contribute <0.04 Gt a⁻¹ to the ocean, accounting for ~14% of Northern Hemisphere frontal ablation. Hence, further studies should focus on the ~200 marine-terminating glaciers contributing more than 0.04 Gt a⁻¹ (see Supplemental Dataset 1), although future glacier changes (e.g., due to surging²¹; tidewater glacier cycle²²) will modify this list.

To investigate which parts of the ocean are most affected by frontal ablation, and thus susceptible to related iceberg hazards and marine ecosystem impacts, we introduce an 'intensity index' that spatially aggregates frontal ablation from all glaciers occurring within 50 km of any point in the ocean on a 10 km grid (Fig. 3). The index depends on both the frontal ablation rates and the number of glaciers in the vicinity of each point, and thus varies in a pattern that is largely consistent with the spatial variation in frontal ablation of glaciers in each ocean region (Fig. 2). Several coastal regions show exceptionally high-intensity indices indicating hotspots of potential iceberg occurrence.

The highest intensity index occurs around northeastern Svalbard, due to the vigorous surge of Basin-3 of Austfonna which started in 2012²³, in addition to substantial contributions from the rest of the ~1200 km² ice cap. High indices are also described in the strait between October Revolution and Komsomolets islands in Severnaya Zemlya (Fig. 3), partly due to recent ice shelf break up¹⁷,²⁴. In Alaska,

high indices close to Hubbard and Columbia glaciers are due to exceptionally high frontal ablation rates of these glaciers; these alone accounted for 63% of regional frontal ablation and are thus the conduits for most ice entering the Gulf of Alaska. Higher intensity indices along the western coast than the eastern coast of Novaya Zemlya is also consistent with a concentration of glaciers with high frontal ablation rates. Despite the relative low frontal ablation rates (mean 0.02 Gt a⁻¹) of glaciers in Franz Josef Land, the large number of glaciers yields a high index.

## Role of ice discharge versus terminus retreat

Quantifying the relative importance of ice discharge versus terminus retreat is key to understanding mass loss mechanisms and highlights where future efforts should be focused to reduce observational uncertainties. Terminus retreat was most important in Arctic Canada South, Arctic Canada North, and Greenland Periphery, where it made up 67%, 48%, and 41%, respectively, of total frontal ablation, while only amounting to 8% in Alaska. Thus, ice discharge is dominant in controlling frontal ablation rates in Alaska.

For a retreating glacier, frontal ablation is larger than ice discharge across the flux gate; conversely, for an advancing glacier, frontal ablation is lower than ice discharge. Two examples in Severnaya Zemlya illustrate the connection between terminus change and frontal ablation. We found that the eight glaciers that fed the now collapsed Matusevich Ice Shelf²⁴ had a frontal ablation rate of 0.20 ± 0.17 Gt a⁻¹ from 2000 to 2010, compared to 1.88 ± 0.23 Gt a⁻¹ in 2010–2020 due to the 18.91 ± 0.28 km² a⁻¹ terminus retreat during that decade,
<figure><img src="image_5.png" /><figcaption>**Fig. 2** | Mean frontal ablation rate by glacier for 2010–2020. Dark gray triangles indicate glaciers that have a frontal ablation <0.02 Gt a⁻¹ (not shown in color bar) and squares indicate glaciers where the uncertainty is more than 50% of total frontal ablation. Glacierized area is marked in white. Latitude and longitude in top left of each panel indicate the center point of that panel and the grid is 1° by 1° for all panels, with 1 degree of latitude equal to 111 km. The only glacier off the color bar</figcaption></figure>

scale is Basin-3 of Austfonna (6.18 Gt a⁻¹; Nordaustlandet). For uncertainties see Fig. S3. Frontal ablation for 2000–2010 is shown in Fig. S2. ACN is Arctic Canada North, ACS is Arctic Canada South, GR is Greenland Periphery, SV is Svalbard, FJL is Franz Josef Land, NZ is Novaya Zemlya, SZ is Severnaya Zemlya; other abbreviations are cardinal directions.

the largest of any Northern Hemisphere glacier. On the other hand, the surging outlet of Vavilov Ice Cap²⁵ discharged 1.24 ± 0.01 Gt a⁻¹, but frontal ablation was only 0.29 ± 0.10 Gt a⁻¹, due to terminus advance. Both these events are unique because the Matusevich Ice Shelf is now gone, and Vavilov Ice Cap is unlikely to surge again in the coming decades as no previous surge events have been identified since at least the 1980s, indicating a very long quiescent phase²⁵,²⁶. For the entire Northern Hemisphere, the total terminus mass loss of 11.57 ± 3.81 Gt a⁻¹ during 2010–2020 consisted of 2.98 ± 0.70 Gt a⁻¹ of terminus advance and 14.54 ± 5.97 Gt a⁻¹ of terminus retreat.

We examined glaciers with advancing termini for morphological indicators²⁷ consistent with dynamic instabilities such as surging

or pulsing (hereafter surge-type glaciers). During the surge phase, these glaciers have elevated velocities, typically by an order of magnitude or more above normal, for months to years, which can lead to increased frontal ablation over a similar period²⁸⁻³⁰. Of the 45 glaciers that advanced during the 2000–2010 period, we found 13 confirmed, two probable, and five possible surge-type glaciers, whereas 25 showed no signs of instability, in total accounting for 6.74 ± 0.40 Gt a⁻¹ of frontal ablation. From 2010–2020, 39 glaciers advanced, of which 14 were confirmed, five probable, and four possible surge-type glaciers, while 16 showed no signs of instability, in total accounting for 11.19 ± 0.71 Gt a⁻¹. Only 13 glaciers, six of which are surge-type, advanced over both 2000–2010 and
<figure><img src="image_6.png" /></figure>

**Fig. 3 |** *Intensity index of frontal ablation in near-coastal ocean environments 2010–2020.* The index is defined as the sum of the frontal ablation rates of all glaciers within 50 km of their terminus, based on a 10 km ocean grid. Text in white indicates values outside the color scale. The pie charts show the proportion of ice discharge (blue) and terminus change (red) to total frontal ablation for each region,

and are scaled to frontal ablation for that region. Glacier-covered area is shown in white, ocean in light blue, and other land surfaces in gray. Individual glaciers or ice caps mentioned in the text are named. The grid shows lines of latitude every 10° and longitude every 5°, or 1106 km between lines of latitude for all panels. The coordinates in the bottom left of each box show the center latitude and longitude.

2010–2020, compared to 745 glaciers that retreated during both time periods.

## Comparisons with regional mass balance

Total glacier mass balance is equal to the climatic-basal balance (the mass changes at the surface, inside, and under a glacier) and frontal ablation¹. We subtract our frontal ablation estimates from total mass balance estimates⁹ derived from elevation change observations of all glaciers in a region to estimate the regional climatic-basal balance. Since these estimates⁹ do not include any mass changes below sea level due to retreat or advance, we correct these geodetic balances assuming two end-member scenarios (see “Methods”) of 5.16 to 9.29 Gt a⁻¹ of ice lost below sea level from 2000 to 2020. This suggests that the Northern Hemisphere glacier net mass loss of 223.7 ± 14.6 Gt a⁻¹ from 2000–2020 by Hugonnet et al.⁹ is currently underestimated by 2–4%, although this portion will not contribute to sea-level rise. We find a total mass loss due to ice discharge equivalent to 2.10 ± 0.22 mm of sea level rise between 2000 and 2020 (accounting for the effect of displacement of ocean water by changes in submarine ice, see “Methods”), with an additional 0.11 ± 0.11 mm coming from terminus mass loss above sea level.

In the Russian Arctic, frontal ablation is greater (by 6.38 ± 5.27 Gt a⁻¹) than the total net mass loss rate (Supplemental Table 1, Fig. S5), suggesting a positive regional climatic-basal balance. We estimate a climatic-basal balance close to zero (+1.71 ± 3.45 Gt a⁻¹) in Svalbard, which has overlapping uncertainties with published estimates³¹,³². All other regions have negative climatic-basal balances. Estimates of climatic-basal balance exist for one other region, Greenland Periphery¹⁶, with which we find agreement.

While it should be possible to compare our results with recent observations of Greenland Ice Sheet ice discharge³³,³⁴, the ice sheet currently lacks discharge estimates for approximately 77 to 211 of the slower-flowing ice sheet outlet glaciers, dependent on how outlet margins are combined or separated. Furthermore, recent studies of Greenland Ice Sheet discharge³³,³⁴ include some peripheral glaciers and do not account for mass changes due to terminus advance or retreat, making the comparison less useful. Even so, our hemispheric ice discharge estimate is an order of magnitude lower than the most recent Greenland Ice Sheet discharge estimates³³,³⁴.

Even though the 1496 marine-terminating glaciers comprise only 0.007% of Northern Hemisphere glaciers by count, our first Northern Hemisphere wide frontal ablation estimate indicates that they play a critical role in global glacier mass balance through mass losses by frontal ablation. Our individual glacier estimates provide valuable information for calibrating global glacier models, as well as for coastal communities affected by icebergs and related navigational hazards, and for assessing glacier impacts on near-coastal marine ecosystems.

## Methods

### Marine-terminating glacier inventory

We identified every glacier that had an interface with the ocean during at least part of a daily tidal cycle by manual examination of Landsat 5, Landsat 7, and ASTER imagery from 1999 to 2020. The only glacier not touching the ocean between 1999 and 2002 was Good Friday Glacier (Arctic Canada North), which advanced to become marine-terminating between 2000 and 2010³⁵, yielding a total of 1496 marine-terminating glaciers in the Northern Hemisphere.
For Greenland Periphery, beyond the ice-sheet margin, we included glaciers that were defined by Rastner et al.³⁶ as having either a weak (CL1) or no (CLO) hydrological connection to the ice sheet. Our total number of glaciers and terminus outlines differs from classification of the RGI due to missing glaciers¹⁰, poorly mapped ones, and misclassified terminus types. Poorly mapped glaciers in RGI6 include those with multiple distinct termini and glaciers that share one terminus (mostly in Svalbard).

## Flux gates

We drew 5016 km of flux gates manually, although after correcting fluxgate width to be perpendicular to flow the total length is 3802 km (Table 1). We drew fluxgates with a preference for near-terminus locations where ice thickness observations exist (Fig. S1). If no thickness observations were available, flux gates were drawn as close to the terminus as practical without entering an area that was lost (median distance from the terminus in 2010 was 630 m). Flux gates are constrained horizontally by a lightly edited (to have one terminus per glacier) version of RGI6¹⁰. Unconstrained ice caps that terminate in the ocean, such as those found in the Russian Arctic, have one flux gate that encircles the ice cap. Flux gates were drawn approximately perpendicular to flow direction based on the orientation of medial moraines, structural patterns of other surface-flow features, and/or of valley walls. We used the x and y component velocity observations to adjust the flux gate width such that the entire width of all flux gates were perpendicular to flow⁴. Flux gates were subdivided into ~25 m segments (hereafter flux gate segments) using a unique orthographic coordinate system for each glacier to minimize the effect of distortion from regional projections. We then used these points to extract values from the velocity and ice thickness datasets.

## Frontal ablation estimation

The velocity, ice thickness, and area change measurements were combined to derive mean frontal ablation rates for every glacier in the Northern Hemisphere over the periods 2000–2010 and 2010–2020. We computed the frontal ablation rate ($A_f$; Gt a⁻¹, defined as positive here) from the ice discharge rate ($\dot{D}_{ice}$, Gt a⁻¹, defined as positive here) across a flux-gate close to the glacier terminus and the mass change rate associated with terminus position change ($M_{term}$; Gt a⁻¹; defined positive for retreat and negative for advance):

$$ \dot{A}_f = \dot{D}_{ice} + \dot{M}_{term} \quad (1) $$

$\dot{D}_{ice}$ is given by:

$$ \dot{D}_{ice} = \left( \rho \left( \sum_{n=1}^{N} (V_n \cdot H_n \cdot d_n) \right) - (S_f \cdot \dot{B}_{clim}) \right) \quad (2) $$

where $\rho$ is the vertically averaged density of the ice column (900 kg m⁻³)³⁷, N is the total number of flux gate segments along the flux gate, $V_n$ is the vertically averaged ice velocity component normal to the flux gate at segment n (assumed to be 95% of the surface velocity at that location; m a⁻¹)³⁷, $H_n$ is the ice thickness (m), $d_n$ is the flux gate segment width, $S_f$ is the area of the region below the flux gate not involved in terminus retreat or advance (m²; area shown in Fig. S1A), and $\dot{B}_{clim}$ is the mean specific climatic mass balance rate (kg m⁻² a⁻¹; defined as positive here).

We calculated the mass change due to terminus advance or retreat by:

$$ \dot{M}_{term} = \left( \rho \cdot \frac{\Delta S_{term}}{\Delta t} \cdot \bar{H} + \frac{\Delta S_{term}}{2} \cdot \dot{B}_{clim} \right) \quad (3) $$

where $\bar{H}$ is the mean thickness of the terminus in the area gained or lost (in m) and $\Delta S_{term}$ is the area (m²) of the glacier that is now replaced by

the ocean for retreating glaciers, or which occupies previous ocean area for advancing glaciers (division by 2 is necessary to account for the average area over the decade as no area was lost at the start of decade), and is computed from the digitized termini positions in 2000, 2010, and 2020 (Fig. S1). $\Delta t$ is the number of years between observations.

## Ice surface velocity

Ice surface velocities were primarily derived from Landsat-based, annual displacement mosaics of the Inter-mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE; Fig. S1C)³⁸ at 240 m resolution³⁹. Data gaps smaller than one pixel were filled from the nearest neighbor. To improve the Greenland Periphery and Russian Arctic observations, we included MEaSURES InSAR⁴⁰ and Sentinel-1 data⁴¹, respectively, when the signal-to-noise ratio was lower than that of ITS_LIVE. For glaciers with no ITS_LIVE observations, such as those above 82.7°N, we used Landsat or Sentinel-2 imagery from 2016–2020 and AutoRIFT⁴² to produce velocity maps. For a small ice cap in the Barents Sea (300 km²) and a glacier in the De Long Islands (10 km²), we could not generate reliable velocity results, and thus assumed a 5 ± 5 m a⁻¹ velocity based on evidence of little ice motion (i.e., no flow bands). For Svalbard and the Russian Arctic, where the optical velocity record from 2000 to 2010 was sparse in ITS_LIVE, we filled this gap using 1992–2012 observations from the JERS1, ERS1, ALOS PALSAR, and TerraSAR-X satellites. We used synthetic-aperture radar (SAR) offset-tracking to compute winter velocity from these scenes at 100 m resolution⁴³. We combined all available velocity observations over each decade to derive a mean decadal velocity along the flux gate for each glacier. We do not account for seasonal variability in ice motion due to a lack of data to quantify this variability. In the few instances for which no velocity data was available for the 2000–2010 period (3.9% of flux gates), mean decadal values from 2010 to 2020 were used to fill the gaps.

## Ice thickness

Direct radar-derived measurements, primarily from the Glacier Thickness Database (GlaThiDa) 3.0.3⁴⁴, were used to determine ice thickness along the flux gates and over deglaciated areas due to terminus retreat since 2000. For Arctic Canada North, Center for Remote Sensing of Ice Sheets⁴⁵ radar tomography observations acquired in ~3 km swaths from the 2014 NASA Operation IceBridge mission provided the main ice thickness source. Additional radio-echo sounding observations were incorporated for the Russian Arctic, where available⁴⁶⁻⁵². Whereas most of our observations were acquired during the study period, we used 1288 Svalbard thickness measurements in GlaThiDa from 1980 to 1995 and some observations in Russia that were collected in 1994⁴⁶ and 1997⁴⁷. To ensure thickness data are contemporaneous, we used elevation change estimates from Hugonnet et al.⁹ to correct ice thickness at the fluxgate from the time of the measurement or model to 2005 and 2015 for use in ice discharge calculations. Where point thickness observations exist, we averaged all observations within 100 m of each flux gate segment to derive the thickness at that point (Eq. 2). When only a centerline observation was available, or the thickness observations did not cover the entire length of the flux gate, we assumed a U-shaped valley to determine thickness where no observations are available⁵³:

$$ H_n = \frac{10 - H_{center}}{W_1} * W_2^2 + H_{center} \quad (4) $$

where $W_1$ (m) is distance from the glacier margin and $W_2$ (m) is distance from the glacier centerline with thickness $H_{center}$ (m). We use 10 m as a minimum ice thickness on the glacier margins. Similar approaches can be found in Sánchez-Gámez and Navarro⁵⁴. This may not accurately reflect the glacier thickness if the centerline observation is not the
deepest point on the glacier, or if the glacier valley is not U-shaped, although we have no way to quantify this.

In all, there is at least one ice thickness observation for 268 glaciers, which together contribute 69% of the total frontal ablation from 2010 to 2020. Out of these 268 glaciers, 250 have at least a centerline observation (middle 20% of glacier), 161 have at least 75% coverage across the flux gate, and 120 have at least 90% coverage across the flux gate. When no thickness observations were available, data gaps were filled with modeling results from Millan et al.⁵⁵ (Fig. S1B). We found that the model results overestimated glacier thickness at the fluxgate, even after accounting for elevation changes with an average 135 m bias between model estimates and observations. We empirically removed this bias by modeling it as linear function of the glacier thickness to match existing measurements.

For five small glaciers where no modeled thickness data were available, we assigned a similar thickness to that of a morphologically similar nearby glacier. For small glaciers, modeling occasionally produces ice thickness values that are unrealistically small (e.g., <5 m average). When the model results suggested that the average ice thickness along the flux gate is <30 m, we replaced these values with 30 m and assigned an uncertainty of 20 m (for 144 glaciers during 2000–2010 and 230 glaciers during 2010–2020). For the six marine-terminating glaciers on Jan Mayen we used scarce point measurements⁵⁵ and flow-law theory³⁷ to estimate the thickness.

To determine mean ice thickness H (Eq. 3) over the glacier area that was lost or gained we averaged all available observations. Where no observations were available, we used modeled thickness estimates (debiased in the same way as above) for the terminus area that was gained or lost unless these data were not available, in which case we assumed the thickness to be 60 ± 30% of the mean thickness along the flux gate. On 74 glaciers, primarily in East Greenland, model thickness estimates were unrealistically low (<10 m) in an area that was gained or lost, so we replaced those values with 60 ± 30% of the mean thickness along the flux gate.

### Glacier terminus front positions

To determine $\Delta S_{term}$ (Eq. 3) for the periods 2000–2010 and 2010–2020, glacier terminus positions were digitized manually for every glacier in approximately 2000, 2010, and 2020, primarily from summer, cloud-free, true-color Landsat 5, 7, and 8 imagery (30 m spatial resolution). When Landsat images were not available, we first attempted to fill the gap with true-color ASTER data (15 m spatial resolution), particularly in northern Greenland and Arctic Canada North for 2000 and 2010. Where no cloud-free optical observations were available, we used Radarsat-1 fine beam mode (2000; 8 m pixel resolution; used in Northern Greenland) and ALOS PALSAR fine beam single polarization (2010; 10 m pixel resolution; used in Arctic Canada North) imagery downloaded from the Alaska Satellite Facility. Mapping was carried out in QGIS version 3.10 in a WGS84 Arctic Polar Stereographic projection (EPSG 3995). To measure the area of each polygon we reprojected the dataset into a unique orthographic projection centered on each glacier to eliminate the impact of distortion.

### Climatic mass balance below the flux gate

To account for mass changes unrelated to frontal ablation (e.g., surface melting) between the flux gate and glacier front), we used the Python Glacier Evolution Model (PyGEM)⁵⁶ to model the climatic mass balance $\dot{B}_{clim}$ for each glacier. $\dot{B}_{clim}$ is defined as the mass changes due to snow accumulation, and melt at the surface and refreezing of melt and rain water below the surface¹. PyGEM assumes the basal mass balance and internal ablation is negligible. We calibrated PyGEM using geodetic mass balance estimates⁹. We used monthly ERA5 air temperature and precipitation data⁵⁷ to compute glacier melt with a degree-day model, accumulation with a temperature threshold, and refreezing from annual air temperature. We extracted $\dot{B}_{clim}$ for the

lowest two 10-m elevation bins of each glacier. The glacier area was assumed constant in PyGEM. For two well-investigated glaciers with high frontal ablation rates we used estimates of climatic balance from direct observations: Columbia Glacier, Alaska (8 ± 2 m water equivalent)⁵⁸, and for Basins 2 and 3 of Austfonna, Svalbard (–0.6 ± 0.3 m water equivalent)³²,⁵⁹. For glaciers lacking appropriate input data (7 small glaciers), we used regional (defined by RGI primary regions) averages of the climatic mass balance.

For 490 glaciers from 2000 to 2010 and 602 glaciers from 2010 to 2020, the absolute $\dot{B}_{clim}$ was greater than the ice discharge, which is physically impossible. For 369 (2000–2010) and 406 (2010–2020) of these glaciers, we increased the average discharge and reduced the average climatic balance within uncertainties to constrain each variable to real-world possibilities. For the remaining glaciers, we used ice discharge values neglecting $\dot{B}_{clim}$, but with 100% uncertainty. These 488 and 601 glaciers only made up 0.70 ± 0.44 Gt a⁻¹ (2000–2010) and 0.54 ± 0.46 Gt a⁻¹ (2010–2020) of ice discharge, accounting for 2% of total ice discharge during both time periods. See Dataset S1 for a list of glaciers that have this correction.

### Surge-type glaciers and instabilities

We individually examined glaciers that exhibited a terminus advance from 2000–2010 and/or 2010–2020 to determine if the advance was due to a dynamic instability such as pulsing or surging. We classified the glaciers according to the RGI classification scale (0: no evidence, 1: possible, 2: probable, 3: observed) based on glacier morphology (including tear-drop moraines, pothole fields, crevasse patterns)²⁷ viewed from satellite imagery and literature where available²⁶,⁶⁰–⁶⁴.

For the Nathorstbreen system and an outlet of the Leningradskiy Ice Cap, our measured discharge was greater than terminus mass gain. However, we adjusted these totals within their uncertainties to make physical sense. For six glaciers from 2000 to 2010 and one from 2010 to 2020, the discharge was less than the terminus change; however, in each instance the mass held in terminus advance was <0.02 Gt a⁻¹, so we artificially increased discharge to be within the uncertainty of terminus change while assigning 100% uncertainty to the discharge. These adjusted glaciers accounted for a total frontal ablation of 0.02 Gt a⁻¹.

### Sea-level rise calculations

We estimated the sea-level equivalent of frontal ablation assuming a global ocean area¹ of 362.5 × 10⁶ km². To account for mass changes below sea level that do not contribute to sea-level change as the ice displaces ocean volume, we subtract the mass losses due to retreat below sea-level or add the gains due to advance. Since the ice fraction of total ice thickness below sea level is unknown, we calculate a higher bound scenario (assuming a fraction of 90% which corresponds to ice close to flotation) and a best estimate lower-bound scenario (50%). Estimates of mass loss below sea level are typically excluded from current global mass balance estimates⁹. Our frontal ablation estimates alone are not able to provide a total sea level rise estimate because frontal ablation is only one component of mass balance.

### Uncertainty in frontal ablation estimates

We estimate the uncertainty in ice discharge $\sigma_D$ assuming each uncertainty component (Eq. 2) is independent, which yields:

$$ \sigma_D^2 = \rho_{ice} \sum_n \left( \sqrt{\left(\frac{\sigma_{V_n}}{V_n}\right)^2 + \left(\frac{\sigma_{H_n}}{H_n}\right)^2 + \left(\frac{\sigma_{d_n}}{d_n}\right)^2} \cdot V_n H_n d_n \right)^2 + (S_f \sigma_{B_{clim}})^2 + (\sigma_{S_f} \dot{B}_{clim})^2 \quad (5) $$

where $\sigma_{V_n}$ is the uncertainty in vertically integrated velocity (assumed to be 95% of the surface velocity at that location)³⁷, $\sigma_{H_n}$ is the
uncertainty in ice thickness, $\sigma_{d_n}$ is the fluxgate width uncertainty, all associated to the flux gate point $n$. At the glacier-scale, $\sigma_{\dot{B}_{clim}}$ is the uncertainty in climatic specific-mass balance, and $\sigma_{S_f}$ is the uncertainty in the area below the flux gate. Similarly, we estimate the uncertainty in terminus mass change $\sigma_{M_{term}}$ assuming each uncertainty component (Eq. 3) is independent, which yields:

$$ \sigma_{M_{term}}^2 = (\sigma_{\Delta S_{term}} H)^2 + (\sigma_H \Delta S_{term})^2 + \left(\sigma_{B_{clim}} \frac{\Delta S_{term}}{2}\right)^2 + \left(\dot{B}_{clim} \frac{\sigma_{\Delta S_{term}}}{2}\right)^2 \quad (6) $$

where $\sigma_{\Delta S_{term}}$ is the uncertainty in the terminus area change and $\sigma_H$ is the uncertainty in the mean thickness over $\Delta S_{term}$. Finally, the uncertainty in frontal ablation $\sigma_A$ is estimated by considering the ice discharge and terminus mass change uncertainties (Eq. 1) as independent:

$$ \sigma_A^2 = \sigma_D^2 + \sigma_M^2 \quad (7) $$

## Uncertainty sources at flux gate point or glacier-scale

We relied on the uncertainties reported by each dataset, whether at the scale of a flux gate point, flux gate, or glacier to propagate the uncertainties through Eqs. 5 and 6.

At flux gate points, we found good agreement (<50% difference) between the average velocity uncertainties $\sigma_{V_n}$ reported in ITS_LIVE and the average dispersion between the ITS_LIVE and MEaSURES estimates at the same flux gates. We found a similar agreement (<30% difference) between measured and modeled ice thicknesses⁵⁵, after accounting for elevation changes⁹ and debiasing (see Supplemental Material for details). The SAR-derived velocity uncertainties average 10 m a⁻¹ for TerraSAR-X and ALOS PALSAR⁴³, 20 m a⁻¹ for JERS-1⁶⁵, and 40 m a⁻¹ for ERS⁶⁶. We derived the uncertainty in flow direction based on the stated uncertainties in the x and y velocity products, translated into a flux gate segment width uncertainty $\sigma_{d_n}$. When these data were unavailable, we assumed a 10° uncertainty in flow direction.

Unless otherwise reported, we assumed a 10% uncertainty in ice thickness $\sigma_{H_n}$. We also assumed a 10% uncertainty for ice thicknesses derived from U-shaped valley-based modeling from a centerline depth measurement³³. When a model depth or thickness observation in the terminus region was not available to calculate terminus mass change, we assumed that the average thickness of the area gained or lost is 60 ± 30% of the average flux gate thickness.

The uncertainty in area gained or lost $\sigma_{S_f}$ was assumed to be one pixel and derived by multiplying the length of the perimeter of the changed terminus area by the width of a pixel (30 m for Landsat and 15 m for ASTER)⁶⁷.

We estimated the uncertainty in the climatic specific-mass balance from the median absolute deviation of 50 simulations from Bayesian calibration of model parameters⁶⁸ forced with ERA5 air temperature and precipitation⁵⁷ calibrated on 2000–2019 glacier mass balances⁹ for each glacier. We assumed the uncertainties in climatic mass balance are independent between glaciers.

## Spatial correlations in ice thickness and velocity uncertainties

We quantified the spatial correlation of velocity and ice thickness uncertainties (Supplemental Fig. 6) by estimating global-scale, empirical variograms⁶⁹,⁷⁰ using the difference between independent sources of estimates at the same flux gates. We performed the spatial correlation analysis at distances covering several orders of magnitude (from 25 m to 1000 km), thereby accounting for biases at glacier and regional scales through long-range correlations⁹.

We compared ITS_LIVE and MEaSURES velocity observations in Greenland (Supplemental Fig. 7) and found fully correlated (100%) variance at short distances (<50 m), and correlated variances at 40–100% within ~700 m. We attribute these short-range correlations to

the resolution of optical imagery and image-matching. Velocity estimates remain correlated at 15–40% within 25 km and 0–15% within 1000 km, highlighting moderate glacier and regional-scale errors due to seasonal differences when estimating yearly velocity. Beyond 1000 km, the velocity estimates are fully decorrelated.

We compared all available ice thickness measurements to debiased model estimates⁵⁵ (Supplemental Fig. 8). Ice thickness variance was correlated at 80–100% within 2 km, the result of short-scale modeling artefacts. Ice thickness variance remains correlated between 50 and 80% within 150 km, and between 0 and 50% within 1000 km, implying large-scale biases in modeled ice thicknesses that are likely owed to the temporal inconsistency of glacier outlines. After 1000 km, ice thickness estimates are fully decorrelated.

More details on spatial correlations are available in the Supplementary Material.

## Uncertainty propagation from pixel to glacier and regional scales

For the ice discharge uncertainty (Eq. 5), we independently considered the ice discharge at the flux gate $\dot{D}_{gate}$ based on velocity and ice thickness estimates (Eq. 1):

$$ \dot{D}_{gate} = \sum_{n} \dot{D}_{n} = \sum_{n} \rho_{ice} V_{n} H_{n} d_{n} \quad (8) $$

We propagated the uncertainties in ice discharge at the flux gate from flux gates $n$ to glaciers, and then repeated the same approach from glacier to region, accounting for spatial correlations in velocity and ice thickness uncertainties⁹,⁷¹:

$$ \sigma_{\dot{D}_{gate}}^2 = \sum_{n} \sigma_{\dot{D}_{n}}^2 + \sum_{n} \sum_{n \neq m} \dot{C}_{\dot{D}}(l_{n-m}) \sigma_{\dot{D}_{n}} \sigma_{\dot{D}_{m}} \quad (9) $$

where $C_{\dot{D}}(l)$ is the spatial correlation in ice discharge at the spatial lag $l$, and $l_{n-m}$ is the distance between flux gates $n$ and $m$ in the glacier, or between glacier $n$ and $m$ in the region. The correlation in ice discharge $C_{\dot{D}}(l)$, ranging from 0 to 1, is estimated either for each glacier or region based on the velocity and ice thickness empirical variograms (Supplemental Figs. 7 and 8).

For the terminus mass change uncertainty, we applied a similar approach to propagate the uncertainty in average ice thickness over the area gained or lost $\sigma_H$, based solely on the ice thickness variogram (Supplemental Fig. 8).

## Data availability

All relevant data outcomes from this study are available in Dataset S1. The python code developed for this project can be downloaded at github.com/willkochtitzky/FrontalAblation. Glacier terminus positions can be downloaded at the Polar Data Catalogue: polardata.ca/pdcsearch/PDCSearch.jsp?doi_id=13257.

## References

1. Cogley, J. G. et al. Glossary of glacier mass balance and related terms. *IHP-VII Tech. Doc. Hydrol. No 86* (2011).

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

3. Sánchez-Gámez, P. et al. Intra- and inter-annual variability in dynamic discharge from the Academy of Sciences Ice Cap, Severnaya Zemlya, Russian Arctic, and its role in modulating mass balance. *J. Glaciol.* **65**, 780–797 (2019).

4. Van Wychen, W. et al. RADARSAT-2 derived glacier velocities and dynamic discharge estimates for the Canadian High Arctic: 2015–2020. *Can. J. Remote Sensing* **46**, 695–714 (2021).

5. Schellenberger, T., Dunse, T., Kääb, A., Kohler, J. & Reijmer, C. H. Surface speed and frontal ablation of Kronebreen and Kongsbreen,
NW Svalbard, from SAR offset tracking. *Cryosph* 9, 2339–2355 (2015).

6. McNabb, R. W., Hock, R. & Huss, M. Variations in Alaska tidewater glacier frontal ablation, 1985-2013. *J. Geophys. Res. Earth Surf.* **120**, 120–136 (2015).

7. Błaszczyk, M., Jania, J. A. & Hagen, J. O. Tidewater glaciers of Svalbard: recent changes and estimates of calving fluxes. *Pol. Polar Res.* **30**, 85–142 (2009).

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

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

10. RGI Consortium. *Randolph Glacier Inventory—A Dataset of Global Glacier Outlines: Version 6.0* (2017).

11. Tepes, P. et al. Changes in elevation and mass of Arctic glaciers and ice caps, 2010–2017. *Remote Sens. Environ.* **261**, 2010–2017 (2021).

12. Ciraci, E., Velicogna, I. & Swenson, S. Continuity of the mass loss of the world's glaciers and ice caps from the GRACE and GRACE Follow-On missions. *Geophys. Res. Lett.* **47**, 9 (2020).

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

14. Noël, B. et al. Six decades of glacial mass loss in the Canadian Arctic Archipelago. *J. Geophys. Res. Earth Surf.* **123**, 1430–1449 (2018).

15. Østby, T. I. et al. Diagnosing the decline in climatic mass balance of glaciers in Svalbard over 1957–2014. *Cryosphere* **11**, 191–215 (2017).

16. Noël, B. et al. A tipping point in refreezing accelerates mass loss of Greenland's glaciers and ice caps. *Nat. Commun.* **8**, 1473 (2017).

17. Dowdeswell, J. A. in *Arctic Ice Shelves and Ice Islands* (eds. Copland, L. & Mueller, D. R.) 55–74 (Springer, 2017).

18. White, A. & Copland, L. Loss of floating glacier tongues from the Yelverton Bay region, Ellesmere Island, Canada. *J. Glaciol.* **65**, 376–394 (2019).

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

20. Kochtitzky, W. & Copland, L. Retreat of Northern Hemisphere marine-terminating glaciers, 2000–2020. *Geophys. Res. Lett.* **49**, 1–10 (2022).

21. Meier, M. F. & Post, A. What are glacial surges? *Can. J. Earth Sci.* **6**, 807–817 (1969).

22. Meier, M. F. & Post, A. Fast tidewater glaciers. *J. Geophys. Res.* **92**, 9051–9058 (1987).

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

24. Willis, M., Melkonian, A. & Pritchard, M. Outlet glacier response to the 2012 collapse of the Matusevich Ice Shelf, Severnaya Zemlya, Russian Arctic. *J. Geophys. Res. Earth Surf.* **120**, 2040–2055 (2015).

25. Willis, M. J. et al. Massive destabilization of an Arctic ice cap. *Earth Planet. Sci. Lett.* **502**, 146–155 (2018).

26. Dowdeswell, J. A. & Williams, M. Surge-type glaciers in the Russian High Arctic identified from digital satellite imagery. *J. Glaciol.* **43**, 489–494 (1997).

27. Jiskoot, H. in *Encyclopedia of Snow Ice and Glaciers* (eds. Singh, V. P., Singh, P. & Haritashya, U. K.) (Springer, 2011).

28. Jiskoot, H., Pedersen, A. K. & Murray, T. Multi-model photogrammetric analysis of the 1990s surge of Sortebræ, East Greenland. *J. Glaciol.* **47**, 677–687 (2001).

29. Mansell, D., Luckman, A. & Murray, T. Dynamics of tidewater surge-type glaciers in northwest Svalbard. *J. Glaciol.* **58**, 110–118 (2012).

30. Van Wychen, W. et al. Characterizing interannual variability of glacial dynamics and dynamic discharge (1999-2015) for the ice masses of Ellesmere and Axel Heiberg Islands, Nunavut, Canada. *J. Geophys. Res. Earth Surf.* **121**, 39–63 (2016).

31. Noël, B. et al. Low elevation of Svalbard glaciers drives high mass loss variability. *Nat. Commun.* **11**, 1–8 (2020).

32. Schuler, T. V. et al. Reconciling Svalbard glacier mass balance. *Front. Earth Sci.* **8**, 156 (2020).

33. Mankoff, K. D. et al. Greenland Ice Sheet solid ice discharge from 1986 through March 2020. *Earth Syst. Sci. Data* **12**, 1367–1383 (2020).

34. King, M. D. et al. Dynamic ice loss from the Greenland Ice Sheet driven by sustained glacier retreat. *Commun. Earth Environ.* **1**, 1–7 (2020).

35. Medrzycka, D., Copland, L., Van Wychen, W. & Burgess, D. Seven decades of uninterrupted advance of Good Friday Glacier, Axel Heiberg Island, Arctic Canada. *J. Glaciol.* **65**, 440–452 (2019).

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

37. Cuffey, K. M. & Paterson, W. S. B. *The Physics of Glaciers* (Academic Press, 2010).

38. Gardner, A. S., Fahnestock, M. A. & Scambos, T. A. *ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities* (Data archived at National Snow and Ice Data Center, 2019).

39. Lei, Y., Gardner, A. & Agram, P. Autonomous repeat image feature tracking (Autorif) and its application for tracking ice displacement. *Remote Sens* **13**, 1–20 (2021).

40. Joughin, I., B. Smith, I. H. & Scambos, T. MEaSURES Greenland ice sheet velocity map from InSAR data, version 2. NASA *Natl Snow Ice Data Cent. Distrib. Act. Arch. Cent.* https://doi.org/10.5067/OC7B04ZM9G6Q (2018).

41. Friedl, P., Seehaus, T. & Braun, M. Global time series and temporal mosaics of glacier surface velocities, derived from Sentinel-1 data. *Earth Syst. Sci. Data* **13**, 4653–4675 (2021).

42. Gardner, A. S. et al. Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years. *Cryosph* **12**, 2 (2018).

43. Paul, P. et al. The glaciers climate change initiative: Methods for creating glacier area; elevation change and velocity products. *Remote Sens. Environ.* **162**, 408 (2015).

44. GlaThiDa Consortium. *Glacier Thickness Database 3.0.1* (World Glacier Monitoring Service, 2019).

45. CReSIS. Multichannel Coherent Depth Sounder (MCoRDS) L3 Grid-ded Data (2016).

46. Dowdeswell, J. A., Gorman, M. R., Glazovsky, A. F. & Macheret, Y. Y. Airborne radio-echo sounding of the ice caps on Franz Josef Land in 1994. *Mater. Glyatsiologicheskikh Issled. Khronika* **80**, 248–254 (1996).

47. Dowdeswell, J. A. et al. Form and flow of the Academy of Sciences Ice Cap, Severnaya Zemlya, Russian High Arctic. *J. Geophys. Res.* **107**, EPM-5 (2002).

48. Macheret, Y. Y., Glazovsky, A. F., Dowdeswell, J. A. & Gorman, M. R. Ice volume change on Franz Josef Land during the last 40 years. *Z. für Gletscherkd. und Glazialgeol.* **35**, 103–116 (2000).

49. Glazovsky, A., Lavrentiev, I. & Vasilenko, E. I. Ice thickness and frontal ablation of nine outlet glaciers on Novaya Zemlya. in *Open Science Conference POLAR 2018 Vol. 2018* (2018).

50. Glazovsky, A., Lavrentiev, I., Vasilenko, E. & Elagina, N. R. N. RES ice thickness and frontal ablation of outlet glaciers in Russian Arctic. in *IASC Workshop on the Dynamics and Mass Budget of Arctic Glaciers & Proglacial Marine Ecosystems 22–24 (2018)*.

51. Sochnev, O. et al. Studies of glaciers in the Russian Arctic for safe marine operations in iceberg waters. *Oil Ind.* **10**, 92–97 (2018).

52. Tarasov, P. et al. Outlet glaciers as iceberg factories: case study for the Kara Sea (2019). in *Proceedings of the Twenty-ninth (2019) International Ocean and Polar Engineering Conference 16–21 (Honolulu, 2019)*.
53. Van Wychen, W. et al. Glacier velocities and dynamic ice discharge from the Queen Elizabeth Islands, Nunavut, Canada. Geophys. Res. Lett. **41**, 484–490 (2014).

54. Sánchez-Gámez, P. & Navarro, F. Ice discharge error estimates using different cross-sectional area approaches: a case study for the Canadian High Arctic, 2016/17. *J. Glaciol.* **64**, 595–608 (2018).

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

56. Rounce, D. R., Hock, R. & Shean, D. E. Glacier mass change in high Mountain Asia through 2100 using the open-source python glacier evolution model (PyGEM). *Front. Earth Sci.* **7**, 1–20 (2020).

57. Hersbach, H. et al. The ERA5 global reanalysis. Q. J. R. Meteorol. Soc. **146**, 1999–2049 (2020).

58. Rasmussen, L. A., Conway, H., Krimmel, R. M. & Hock, R. Surface mass balance, thinning and iceberg production, Columbia Glacier, Alaska, 1948-2007. *J. Glaciol.* **57**, 431–440 (2011).

59. Hagen, J. O., Melvold, K., Pinglot, F. & Dowdeswell, J. A. On the net mass balance of the glaciers and ice caps in Svalbard, Norwegian Arctic. *Arct., Antarct. Alp. Res.* **35**, 264–270 (2003).

60. Clarke, G. K. C., Schmok, J. P., Ommanney, C. S. L. & Collins, S. G. Characteristics of surge-type glaciers. *J. Geophys. Res.* **91**, 7165–7180 (1986).

61. Copland, L., Sharp, M. J. & Dowdeswell, J. A. The distribution and flow characteristics of surge-type glaciers in the Canadian High Arctic. *Ann. Glaciol.* **36**, 73–81 (2003).

62. Grant, K. L., Stokes, C. R. & Evans, I. S. Identification and characteristics of surge-type glaciers on Novaya Zemlya, Russian Arctic. *J. Glaciol.* **55**, 960–972 (2009).

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

64. Jiskoot, H., Murray, T. & Luckman, A. Surge potential and drainage-basin characteristics in East Greenland. *Ann. Glaciol.* **36**, 142–148 (2003).

65. Strozzi, T. et al. Estimation of Arctic glacier motion with satellite L-band SAR data. *Remote Sens. Environ.* **112**, 636–645 (2008).

66. Dowdeswell, J. A., Benham, T. J., Strozzi, T. & Hagen, J. O. Iceberg calving flux and mass balance of the Austfonna ice cap on Nordaustlandet, Svalbard. *J. Geophys. Res. Earth Surf.* **113**, https://doi.org/10.1029/2007JF000905 (2008).

67. Krumwiede, B. S. et al. in *Global Land Ice Measurements From Space* (eds. Kargel, J., Leonard, G., Bishop, M., Kääb, A., R. B.) 481–508 (Springer, 2014).

68. Rounce, D. R. et al. Quantifying parameter uncertainty in a large-scale glacier evolution model using Bayesian inference: application to High Mountain Asia. *J. Glaciol.* **66**, 175–187 (2020).

69. Mälicke, M. & Schneider, H. D. Scikit-GStat 0.2.6: a scipy flavored geostatistical analysis toolbox written in Python. https://doi.org/10.5281/ZENODO.3531816 (2019).

70. Matheron, G. *Les variables régionalisées et leur estimation. Une application de la théorie des fonctions aléatoires aux sciences de la nature*. Masson et Cie (1965).

71. Webster, R. & Oliver, M. A. *Geostatistics for Environmental Scientists: Second Edition*. https://doi.org/10.1002/9780470517277 (2008).

## Acknowledgements

The start of this project and early ideas that informed this work were a result of the “Workshop on the Importance of Calving for the Mass

Balance of Arctic Glaciers” held in Sopot, Poland in 2016, sponsored by the Centre for Polar Studies, University of Silesia and the IASC Network on Arctic Glaciology. We thank Ruitang Yang for her help in improving the error analysis. W.K. acknowledges support from the Vanier Graduate Scholarship. L.C. thanks the Natural Sciences and Engineering Research Council of Canada, University of Ottawa and ArcticNet Network of Centres of Excellence Canada for funding. T.S. acknowledges support from the ESA Glaciers CCI project. A.G. and I.L. thank RG State Contract FMGE-2019-0004, and H.J. the Natural Sciences and Engineering Research Council of Canada. Re.Ho. and D.R.R. acknowledge support from the NASA grant 80NSSC20K1296.

## Author contributions

W.K. did most of the data collection, performed the analyses and calculations, and wrote the initial manuscript draft. L.C. aided in data collection and helped write the initial manuscript. W.W., J.D., A.C., A.D., H.J., J.J., and F.N. helped with initial project planning and contributed to the final manuscript. Re.Ho. contributed substantially to writing the final manuscript draft. J.D., T.B., A.G., and I.L. provided ice thickness data for the Russian Arctic. T.Z. provided velocity observations in Svalbard and the Russian Arctic. H.J. and J.C. provided frontal ablation estimates for Jan Mayen. D.R.R. provided melt model results, R.M. provided glacier thickness data, W.K., Re.Ho., and Ro.Hu. developed the error analysis and RoHu completed it.

## 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/s41467-022-33231-x.

**Correspondence** and requests for materials should be addressed to William Kochtitzky.

**Peer review information** *Nature Communications* thanks Michael Wood and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**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 4.0 International License, which permits use, sharing, adaptation, 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 license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license 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 license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022