Annual mass change of the world's glaciers from 1976 to 2024 by temporal downscaling of satellite data with in situ observations

Glaciers, distinct from the Greenland and Antarctic ice sheets, play a crucial role in Earth's climate system by affecting global sea levels, regional freshwater availability, nutrient and energy budgets, and local geohazards. Past assessments of regional to global glacier mass changes were limited in spatial coverage, temporal resolution, and/or temporal coverage. Here, we present a new observation-based dataset of glacier mass changes with global coverage and annual resolution from 1976 to 2024. We use geostatistical modeling for the temporal downscaling of decadal glacier-wide elevation change estimates derived from satellite and airborne geodetic data, with glaciological annual in situ observations. In more detail, we spatially interpolate the annual mass balance anomalies from sparse in situ observations and calibrate them to glacier-wide long-term trends from elevation change observations available for individual glaciers for varying time periods and with global glacier coverage from 2000 to 2019. We then extrapolate the results to yearly time series starting between 1915 and 1976, depending on the regional data availability, and extending to 2024. The time series are calculated separately for each of the world's glaciers and then aggregated to gridded (0.5° latitude and longitude), regional, and global estimates of annual glacier mass changes. Since 1976, Earth's glaciers have lost 9179 ± 621 Gt (187 ± 20 Gt per year) of water, contributing 25.3 ± 1.7 mm (0.5 ± 0.2 mm per year) to the global mean sea level rise. About 41 % (∼ 10 mm) of this loss occurred in the last decade, with 6 % (∼ 1.5 mm) occurring in 2023 alone, the record-breaking year of glacier mass loss. We review the strengths and limitations of our new dataset, validate and discuss related uncertainty estimates in a leave-one-out/block-out cross-validation exercise, and compare our results to earlier assessments. The annual mass change time series for individual glaciers and the derived global gridded annual mass change product are available from the World Glacier Monitoring Service (WGMS) at https://doi.org/10.5904/wgms-amce-2025-02 (Dussaillant et al., 2025).

1 Introduction

Glacier monitoring has evolved quickly since its beginning in the late 19th century, with in situ and remotely sensed techniques allowing one to observe detailed changes in area, elevation, volume, and mass, at first only for single glaciers and recently entire regions and the world (Berthier et al., 2023; The GlaMBIE Team, 2025; Thomson et al., 2021; Zemp et al., 2015). In situ, glaciological observations provide extremely valuable information on the annual-to-seasonal temporal variability of glacier changes, reflecting the impact of atmospheric conditions which can be correlated over several hundred kilometers (Braithwaite and Hughes, 2020; Cogley et al., 2011; Fernández and Somos-Valenzuela, 2022; Kaser et al., 2003; Oerlemans, 2001; Østrem and Brugman, 1991). At present, and thanks to the coordination of the World Glacier Monitoring Service (WGMS) and its global network of contributors, glaciological in situ observations exist in nearly all glacierized regions of the Randolph Glacier Inventory (RGI, Pfeffer et al., 2014; RGI Consortium, 2017). Most glaciers have been monitored continuously for periods longer than 10 years, with some of the earliest observations reaching back to the early 20th century. While irreplaceable, one major limitation of the glaciological method lies in the logistical hurdles of maintaining continuous field campaigns. At present, in situ observations are limited to approximately 500 glaciers worldwide, representing less than 0.2 % of the world's glaciers (WGMS, 2025). Secondly, it is challenging to represent the complex mass balance pattern with individual in situ point measurements such that potential sampling biases can accumulate in time when interpolating to glacier-wide estimates. For this reason, glaciological observations often require reanalysis and calibration, with glacier elevation change rates obtained from high-resolution geodetic surveys (Thibert et al., 2008; Thibert and Vincent, 2009; Zemp et al., 2013).

The geodetic or digital elevation model (DEM) differencing method is powerful at providing glacier elevation change observations with high accuracy over large glacierized areas and long periods of time (multi-annual to decadal, Cogley et al., 2011). DEM differencing was initially applied to individual glaciers with DEMs derived from maps (Joerg and Zemp, 2014) and aerial photographs (Belart et al., 2019; Finsterwalder, 1954; Papasodoro et al., 2015; Thibert et al., 2008), but it has now evolved to include data from airborne lidar (Abermann et al., 2010; Echelmeyer et al., 1996), spaceborne altimetry (Jakob and Gourmelen, 2023; Menounos et al., 2024), and satellite-derived DEMs from multiple sensors (Berthier et al., 2023; Toutin, 2001). Recent advances in post-processing techniques (Hugonnet et al., 2022; McNabb et al., 2019; Nuth and Kääb, 2011; Rolstad et al., 2009), supercomputing capabilities, and automated processing pipelines (Girod et al., 2017; Rupnik et al., 2017; Shean et al., 2016) have further enhanced this methodology, enabling its application over entire mountain ranges (Brun et al., 2017; Braun et al., 2019; Dussaillant et al., 2019; Shean et al., 2020) and, recently, globally (Hugonnet et al., 2021). The major limitations of the geodetic method lie firstly in the relatively short period during which corresponding spaceborne sensors have operated (in general after 2000), sensor-related issues (e.g., radar signal penetration into snow and ice), and, importantly, the inability to capture the annual variability of glacier mass changes due to a low signal-to-noise ratio of the elevation changes and the high uncertainties of the volume-to-mass conversion for periods shorter than 5 years (Huss, 2013).

Glacier change observations using glaciological measurements and the geodetic method therefore complement each other by providing different types of information. Zemp et al. (2019) conducted the first study to combine the annual variability from the glaciological observations with the long-term trends of the geodetic method in order to estimate annual mass changes for all 19 RGI glacier regions from 1976 to 2016. Their global estimate was hampered by the limited geodetic observational sample available at the time of the study (only 9 % of Earth's glaciers by number), resulting in high uncertainties. The Hugonnet et al. (2021) global assessment now fills this observational gap by leveraging the repeated acquisitions and global coverage of the Advanced Spaceborne and Thermal Emission and Reflection (ASTER) satellite optical stereo images (Raup et al., 2000). Their assessment provides individual glacier elevation change rates for nearly all glaciers worldwide (97.4 % of the RGI glacier area) from 2000 to 2019.

In this study, we provide a global observation-based assessment to estimate annual glacier mass changes at a glacier-specific level using geostatistical modeling, feasible thanks to the now almost complete coverage with glacier elevation change observations in the latest version of the Fluctuations of Glaciers (FoG) database (WGMS, 2025). To achieve this, we use glaciological observations from approximately 500 glaciers (0.2 % of the world's glaciers) starting between 1915 and 1976 and glacier-wide geodetic observations from approximately 207 000 glaciers (96 % of the world's glaciers covering 96 % of the world's glaciated surface) starting in the 1940s (Table 1 and Fig. 1). The DEM differencing observations used here include the multiple local and regional satellite and airborne geodetic glacier change assessments for 30 000 (14 %) glaciers plus the 20-year estimates from Hugonnet et al. (2021) available for 205 000 (95 %) glaciers in the FoG database. Building upon the methodological foundations laid out in Zemp et al. (2019, 2020), we further develop the approach to spatially interpolate the glaciological annual in situ observations and calibrate them to the long-term trends derived from satellite and airborne geodetic data elevation change observations. The time series are calculated separately for each of the world's glaciers using geostatistical modeling and are then aggregated to regional and global estimates of annual glacier mass changes. Our results include glacier mass changes with the annual temporal resolution for each individual glacier in the RGI (with the starting date between the hydrological years 1915 and 1976 (see the “Data and methods” section), depending on the RGI region) and a global observation-based assessment of annual glacier mass changes since the hydrological year 1976, made available as a global gridded product. Our dataset holds great potential for contributions to internationally coordinated intercomparison exercises such as The GlaMBIE Team (2025).

Table 1 Key characteristics of the glacier elevation and mass change observations used in this study as available from the Fluctuations of Glaciers database (WGMS, 2025).

^* GCOS (2022).

Figure 1 Spatial and hypsometric coverage of glaciological and geodetic observations for each of the 19 first-order regions. The glacier hypsometry from RGI 6.0 (grey) is overlaid (and almost hidden) with the glacier hypsometry of the geodetic elevation changes (Geo, blue) and the glaciological samples (Glac, red) available from the FoG database (WGMS, 2025). Values for the glacier area and the number (N) of glaciers are given for each region, together with the respective percentage area covered and the number of observed and unobserved glaciers (N_unobs). Grey circles overlaid on the map represent the regional glaciated area. The number of used glacier anomalies is noted per region. ^* Glacier outlines in region 12 are from Tielidze and Wheate (2018), as available from GLIMS. ${}^{**}$ Glaciers strongly connected to the Greenland Ice Sheet in region 5 are excluded.

2 Data and methods

2.1 Input datasets

2.1.1 Glacier inventories

We use the digital glacier outlines from the Randolph Glacier Inventory 6.0 (RGI Consortium, 2017) to spatially locate glaciers, attribute their area, and assign the in situ mass balance observations to individual glaciers. RGI 6.0 represents glacier areas near the beginning of the 21st century. It has been preferred to the more recent version 7.0 for two reasons: first, glacier elevation change observations are available only for version 6.0, and second, it is used for comparison with previous observation-based global assessments, also with this version. RGI outlines are available through the RGI portal (https://doi.org/10.7265/N5-RGI-60, GLIMS Consortium, 2005) and Global Land Ice Measurements from Space (GLIMS), an initiative from the early 2000s to improve glacier inventories using satellite data (https://doi.org/10.7265/N5V98602, RGI Consortium, 2017). Due to the high number of incorrectly mapped glaciers in the RGI 6.0 Caucasus and Middle East region (region 12), the Hugonnet et al. (2021) geodetic observations were calculated using the latest available GLIMS outlines (Tielidze and Wheate, 2018). The latter inventory is also used in this study for consistency. For the Greenland Periphery (region 5), we excluded glaciers strongly connected to the ice sheet (RGI 6.0 connectivity level 2), similarly to Hugonnet et al. (2021). To spatially constrain glaciers within the world's climatic regions, we use the 19 first-order glacier regions as defined by the RGI and illustrated in Fig. 1. Glacier regions are implemented directly in the RGI dataset and are accessible via the same DOI. The locations of the glacier and region outlines are represented in Fig. 1. The full regional hypsometric coverage is illustrated in grey bars in Fig. 1.

2.1.2 Glacier elevation and mass change observations

We use the glacier-wide annual mass change observations from the glaciological method and multi-annual trends of elevation changes derived from the geodetic method as available from the latest update of the FoG database. These glacier change observations are collected by the WGMS in annual calls for data through a worldwide network of national correspondents and principal investigators. After integration of the new, homogenized, and corrected observations, a new FoG database version is released. Individual glaciers with available observations are identified in the FoG database with a WGMS-Id. Updated versions of the FoG database can be accessed via the WGMS online portal (https://wgms.ch/data_databaseversions/, last access: May 2025). The results presented here use version WGMS (2025), which is accessible at https://doi.org/10.5904/wgms-fog-2025-02. The almost complete observational coverage of the latest FoG database version WGMS (2025) is depicted in Fig. 1, with the geodetic sample covering 206 554 glaciers or 96 % of the world's glaciers in the glaciers' entire elevation range. More specifically, in this study, we use local and regional satellite and airborne glacier-wide DEM differencing observations available for 29 529 glaciers (14 % of the world's glaciers) plus the 20-year estimates from Hugonnet et al. (2021) available for 205 120 (95 %) glaciers in the FoG database (only glacier-wide estimates from Hugonnet et al. (2021) calculated from elevation change grids covering more than 50 % of the glacier area were incorporated into the FoG). The key characteristics of glacier mass and elevation change observations are summarized in Table 1. For more details on the specific input data, auxiliary data, retrieval algorithms, and uncertainty estimations of the independent FoG glacier elevation and mass change observations, please refer to WGMS (2025). More details on the glaciological method can be found in Østrem and Brugman (1991), Kaser et al. (2003), and Zemp et al. (2013, 2015). For the geodetic method and its error sources, see WMO (2023), and about measuring glacier mass changes from space, see Berthier et al. (2023).

2.2 Methods

To prepare the data for the main processing, glacier-wide records with available mass and elevation change observations identified by a WGMS-Id are selected from the FoG database and related to their corresponding RGI outline identifier (RGIId) using a link-up table. We exclude geodetic records with survey periods shorter than 5 years in view of their large uncertainty for the volume-to-mass conversion (Huss, 2013). For simplicity, throughout this work hydrological years are represented as the last year of the hydrological cycle (e.g., 1976) from 1 October to 30 September in the Northern Hemisphere and from 1 April of the previous year (e.g., 1975) to 31 March of the next year (e.g., 1976) in the Southern Hemisphere. For the Low Latitudes region, we assume the hydrological year to be equal to the calendar year from 1 January to 31 December.

Our processing algorithm is summarized in three key steps described in the following sections and in Fig. 2. First, focusing on a specific glacier in RGI 6.0, we estimate the detrended temporal variability of annual mass change for the glacier, referred to here as the glacier mean annual mass balance anomaly, using the interannual variability of nearby glaciological time series (Sect. 2.2.1 and Fig. 2a). Secondly, we calibrate the mean annual mass balance anomaly to the long-term trends from the geodetic sample available for the respective glacier (Sect. 2.2.2 and Fig. 2b). Third, we integrate all of these calibrated time series into a single weighted average, producing a data-fused mean calibrated annual mass change time series unique to every individual glacier (Sect. 2.2.3 and Fig. 2c). In a fourth step, the mean calibrated annual mass change time series are integrated into larger regions (e.g., grid point, glacier regions, and global) to obtain a region-specific calibrated annual mass change time series, accounting for spatial correlation of error sources. All given uncertainties in the tables, figures, and main text and reported in the dataset files are at the 1σ level (68 % confidence interval), unless stated otherwise.

Figure 2 Methodological steps illustrated at Hintereisferner (Central Europe, upper row), Gulkana (Alaska, middle row), and Mittivakat (Greenland Periphery, lower row). (a) Mean annual mass balance anomaly (black line) and uncertainty (grey). Grey lines correspond to the spatially selected individual glacier annual mass balance anomalies. The mean annual mass balance anomaly is calculated with the kriging spatial correlation function. (b) Calibration of the mean annual mass balance anomaly over geodetic mass balance observations available for each glacier (red and blue lines for negative and positive change rates, respectively). Grey lines correspond to the individual calibrated time series for each geodetic mass change observation. (c) Mean calibrated annual mass change time series (black line) and uncertainty (grey). We note that the full time series for Hintereisferner starts in 1915, and for Gulkana in 1946, for visualization purposes the plots start in 1950. Notations for the plot equations are defined in the text.

2.2.1 Computing the mean annual mass balance anomaly from the neighboring glaciological observations

Direct annual glaciological observations B_a are reported to the FoG database with their relative uncertainties ${\mathit{\sigma }}_{B_{\mathrm{a}}}$ in meters water equivalent (m w.e.) as

$$\begin{array}{}\text{(1)}& B_{\mathrm{a}}\pm {\mathit{\sigma }}_{B_{\mathrm{a}}}.\end{array}$$

In cases where a glaciological series is missing an uncertainty estimate for a given year, we assume it to be equal to the mean of the annual uncertainty estimates within the series. In cases where a glaciological series has no uncertainty estimate for the entire period, we assume it to be equal to the mean annual uncertainty for all glaciological series from glaciers belonging to the same region.

From the individual glacier annual glaciological time series, we estimate an individual glacier annual mass balance anomaly as the glaciological mass change value at year Y minus the mean mass change during the reference period from 2011 and 2020. Choosing this recent reference period allows for exploitation of a larger glaciological sample, thus obtaining a better representativeness of glacier temporal variabilities across all regions. We allow a threshold of at least 8 years of glaciological observations within the 10-year reference period to calculate a glacier annual mass balance anomaly. This means that a glacier needs to have at least 8 years of glaciological in situ observations within the 10-year reference period to calculate their annual mass balance anomaly. In this step we remove low-confidence glacier anomalies from the processing for not being representative of their regional mass balance variability (Table 2). We use a total of 158 individual glacier anomalies for the assessment.

$$\begin{array}{}\text{(2)}& {\mathit{\beta }}_Y=B_{\mathrm{a},Y}-{\stackrel{\mathrm{‾}}{B}}_{\mathrm{a},\text{2011–2020}}\end{array}$$

Starting from a given glacier g belonging to RGI 6.0 (e.g., Hintereisferner, Gulkana, and Mittivakat in Fig. 2), we use the sample of neighboring glacier annual mass balance anomalies to capture the annual temporal variability of its changes. The mean annual mass balance anomaly of glacier g ${\stackrel{\mathrm{‾}}{(\mathit{\beta }}}_{g,Y}$) is then calculated by kriging all glacier annual mass balance anomalies located near the glacier (Fig. 2a):

$$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{g,Y}\pm {\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{g,Y}}=K({\mathit{\beta }}_{i,Y},{\mathit{\rho }}_{\mathit{\beta },Y}(d\left)\right),\end{array}$$

where ρ_β,Y(d) is the spatial correlation of the annual mass balance anomaly (see further below) and K the function applying ordinary kriging to β_i,Y (i.e., kriging with a constant mean, for which we used the PyKrige Python package, https://doi.org/10.5281/zenodo.3738604, Murphy et al., 2024). When interpolating with kriging, the weight between each pair of glaciers is based solely on the distance between them, using a spatial correlation function which can be constrained from an empirical variogram. Additionally, the predicted kriging uncertainty ${\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{g,Y}}$ grows with distance, from the measurement error of the inputs ${\mathit{\sigma }}_{B_{\mathrm{a}}}$ at close distances from a measured glacier to the signal variability (spread of β_Y) at distances far away from any measured glacier, where the prediction is more poorly constrained. Kriging is a core method of spatial statistics (Cressie, 2015) often termed the “best linear unbiased interpolator” due to its nonparametric nature and empirical variance minimization. It emerged in mining applications (Matheron, 1963) and has since become ubiquitous for spatial interpolation across many fields (Webster and Oliver, 2007). In glaciology, kriging has for instance been used to spatially interpolate sparse ablation measurements (Hock and Jensen, 1999) or ice thickness measurements (Fischer, 2009). Recently, the rise of machine learning methods extended the kriging concept to any kind of dimension through Gaussian processes (Rasmussen and Williams, 2006), which have also found applications in glaciology, from remote sensing time series interpolation (Hugonnet et al., 2021) to model error emulation (Edwards et al., 2021).

Table 2 Regional overview of (i) excluded glacier annual mass balance anomalies, (ii) complementary glacier annual mass balance anomalies added in undersampled regions to calculate the mean annual mass balance anomalies, and (iii) complementary annual mass balance glacier anomalies, normalized and used to fill up past years only to extend the series back in time until at least 1976.

In order to estimate the spatial correlation of the annual mass balance anomaly ρ_β,Y(d) to constrain the kriging, we sampled empirical variograms for both local-scale modeled annual mass balance anomalies (Huss and Hock, 2015) and observational 5-year anomalies (Hugonnet et al., 2021), the latter validating the spatial correlation patterns observed in the modeled estimates (Fig. 3). As climatic patterns driving correlations in regional anomalies should have a size that is largely consistent in time, we expect 5-year anomalies to be spatially correlated at distances similar to annual anomalies but with a lesser amplitude due to the canceling of positive and negative anomalies over time. We indeed identify that both anomalies have a significant spatial correlation up to 5000 km, with a smaller amplitude for 5-year anomalies (Fig. 3). The empirical variograms were computed independently for each contiguous region (High Mountain Asia, North America, and South America being considered larger regions) and for each subperiod of 2000–2020 (annual or 5-year), using up to 10 000 glaciers. In total, all variograms sample more than 10 billion pairwise differences between glacier annual mass balance anomalies to obtain this average spatial correlation function. We then modeled the spatial correlation in the annual mass balance anomaly ρ_β,Y as a sum of two exponential models:

$$\begin{array}{}\text{(4)}& {\mathit{\rho }}_{\mathit{\beta },Y}\left(d\right)=s_{\mathrm{1}}{e}^{-\frac{\mathrm{3}d}{r_{\mathrm{1}}}}+s_{\mathrm{2}}{e}^{-\frac{\mathrm{3}d}{r_{\mathrm{2}}}}\text{if}d>\mathrm{0},\text{else}\mathrm{1},\end{array}$$

where d is the distance between two glaciers, s₁=0.37 and s₂=0.59 are the partial sills, and r₁=200 and r₂=5000 km are the correlation ranges. With this correlation function, we estimate that two directly neighboring glaciers have an annual mass balance anomaly correlated at 96 %, 72 % if separated by 60 km, 51 % if separated by 250 km, 32 % if separated by 1000 km, 10 % if separated by 3000 km, and less than 3 % after 5000 km.

Figure 3 Spatial correlation of anomalies in specific mass change as used for kriging. The correlation is a global average estimated from regional empirical correlation sampled from all glaciers worldwide, for 5-year mass balance anomalies with observational estimates based on surface elevation changes (Hugonnet et al., 2021) and for annual mass balance anomalies for modeled estimates from Huss and Hock (2015). The former is used for validation and the latter for kriging in this study. For example, two directly neighboring glaciers have an annual mass balance anomaly correlated at 96 %, 72 % if separated by 60 km, 51 % if separated by 250 km, at 32 % if separated by 1000 km, 10 % if separated by 3000 km, and less than 3 % after 5000 km (blue dotted line).

We note that, because kriging is a nonparametric interpolation method, its prediction primarily depends on the observations themselves, so uncertainties in the correlation function stemming from the modeled estimates of Huss and Hock (2015) have little influence on our results. Furthermore, because our correlations span multiple orders of magnitude (from 10 to 5000 km), the choice of functional form of the correlation has been shown to have a minimal impact on the prediction (Hugonnet et al., 2022). To exemplify this, we compared kriging with inverse-distance weighting, a different interpolation method altogether, and found almost equal regional estimates as those are primarily driven by the input data. Differences between kriging and inverse-distance weighting only showed at the glacier scale, where kriging allows us to refine anomalies further and derive empirical uncertainties.

In undersampled regions (Arctic Canada South, Russian Arctic, Asia South East, Asia South West, and New Zealand), we added complementary glacier anomalies from neighboring regions to calculate the mean glacier anomalies (Table 2, Fig. 4). As a rule, all glacier mean annual mass balance anomalies should cover at least the period between the hydrological years from 1976 to 2024. For glaciers with mean annual mass balance anomalies not going back to 1976, the best-correlated glaciological series from neighboring regions (i.e., climatically similar) are used to fill in the past years only (Table 2, grey sections in Fig. 4, metadata file, decision supported by Zemp et al., 2019, 2020; Braithwaite and Hughes, 2020; Fernández and Somos-Valenzuela, 2022). To reduce the effect of possible climatic differences within the neighboring regions, the amplitude of the complementary glacier anomalies is normalized to the amplitude of the mean glacier anomaly during the reference period.

2.2.2 Calibrating the mean annual mass balance anomaly on the glacier geodetic sample

Geodetic observations are reported to the FoG database with their relative uncertainties as glacier-wide mean rates of elevation change ($\frac{\mathrm{d}h}{\mathrm{d}t}$) (m) during a specific period of record (PoR). Glaciers may contain multiple individual geodetic observations for different time periods, depending on the dates of the DEMs used (red and blue bars in Fig. 2b). To obtain the geodetic mass change rate, we convert elevation changes to glacier-specific mass change rates (m w.e.) by applying a density conversion factor $f_{\mathit{\rho }}\pm {\mathit{\sigma }}_{f_{\mathit{\rho }}}=\mathrm{850}\pm \mathrm{60}$ kg m⁻³ (Huss, 2013). In this step we exclude from the processing unpublished, and therefore low-confidence, DEM differencing estimates available from the FoG database.

$$\begin{array}{}\text{(5)}& {\stackrel{\mathrm{‾}}{B}}_{\mathrm{geo},\mathrm{PoR}}={\stackrel{\mathrm{‾}}{\mathrm{d}h}}_{\mathrm{PoR}}\cdot f_{\mathit{\rho }}\end{array}$$

To calculate the uncertainty in the geodetic mass balance rate (m w.e.), we propagate the uncertainty in the elevation change rate ${\stackrel{\mathrm{‾}}{\mathrm{d}h}}_{\mathrm{PoR}}$ and the uncertainty in the density conversion factor ${\mathit{\sigma }}_{f_{\mathit{\rho }}}$, considering that they are uncorrelated. We justify this by the fact that elevation change errors stem from instrument noise or spatiotemporal prediction, while density conversion errors stem from modeling errors and lack of knowledge of surface conditions, which are independent. We use reported uncertainties for elevation change rates, and we use ${\mathit{\sigma }}_{f_{\mathit{\rho }}}=\mathrm{60}$ kg m⁻³ (Huss, 2013) for the uncertainty in the density conversion factor.

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{B}}_{\mathrm{geo},\mathrm{PoR}}}=\left|{\stackrel{\mathrm{‾}}{B}}_{\mathrm{geo},\mathrm{PoR}}\right|\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathrm{d}h}}_{\mathrm{PoR}}}}{{\stackrel{\mathrm{‾}}{\mathrm{d}h}}_{\mathrm{PoR}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\rho }}}{f_{\mathit{\rho }}}}\right)}^{\mathrm{2}}}\end{array}$$

In a data-fusion step, we calibrate the mean annual mass balance anomaly (obtained from the glaciological sample) of glacier g to a given geodetic mass change rate k (i.e., geodetic sample) belonging to glacier g. We obtain a “k-calibrated” annual mass change time series for every geodetic observation available for glacier g. The k-calibrated annual mass change time series is then calculated as the sum of the geodetic mass change rate k and the mean annual mass balance anomaly over the period of record of the geodetic mass change rate k (grey lines in Fig. 2b). Due to the large uncertainties related to the volume-to-mass conversion factor over short periods of time (Huss, 2013), only geodetic observations longer than 5 years are considered for calibration.

$$\begin{array}{}\text{(7)}& B_{\mathrm{cal},k,Y}={\stackrel{\mathrm{‾}}{B}}_{\mathrm{geo},k,\mathrm{PoR}}+\left({\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{g,Y}-{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{g,\mathrm{PoR}}\right)\end{array}$$

The uncertainty in the k-calibrated annual mass change time series is then propagated from the uncertainty in the individual geodetic mass change rates and the uncertainty in the mean annual anomaly. We consider these two uncertainties to be uncorrelated as they originate from independent measurements (remotely sensed and in situ) that do not share any similar error sources (including density conversion).

2.2.3 Combining the resulting time series into a mean calibrated annual mass change time series

The mean calibrated annual mass change of glacier g (Fig. 2c) is finally calculated as the weighted mean of the K available k-calibrated annual mass change time series (K is equal to the number of geodetic observations longer than 5 years that are available for the glacier). Two different weights are applied. First, there is a weighting relative to the uncertainty in the geodetic mass change ($W_{{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{B}}_{\mathrm{geo},k,\mathrm{PoR}}}}$), where the k-calibrated annual mass change time series are weighted to the inverse of that squared uncertainty. Second, there is a weighting relative to the time lag of the k-calibrated annual mass change time series (W_t) in number of years to the initial (y₀) and final (y₁) years of the period of record of the k geodetic rate. Here, the point is to give more weight to a time series calibrated to a geodetic observation that is temporally close to the given geodetic period of record:

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}={\displaystyle \frac{{\sum }_{k=\mathrm{1}}^K{B}_{\mathrm{cal},k,Y}\cdot W_{{\mathit{\sigma }}_{B_{\mathrm{geo},k,\mathrm{PoR}}}}\cdot W_t}{K}},\end{array}$$

where $W_{{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{B}}_{\mathrm{geo},k,\mathrm{PoR}}}}=\frac{\mathrm{1}}{{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{B}}_{\mathrm{geo},k,\mathrm{PoR}}}^{\mathrm{2}}}$ and $W_t={\left(\frac{\mathrm{1}}{t_Y}\right)}^{p}t=\mathrm{1},y_{\mathrm{0}}<Y<y_{\mathrm{1}}t=Y-y,Y>y_{\mathrm{1}}t=y_{\mathrm{0}}-Y,Y<y$, and p=1 to give temporally distant years a moderate weight.

For later error propagation, we separate the three sources of errors for elevation change, density conversion, and anomaly calculation, given that these are largely independent between themselves, as previously justified, but they have different scales of spatial correlations. To propagate uncertainties to the mean calibrated annual mass change time series ${\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}$, we consider each source to be entirely correlated with itself during the averaging of the different calibrated series. This is justified because there is a single density conversion and anomaly estimation referring to the same period and glacier, which thus have the same errors. For elevation changes, there are sometimes multiple estimations referring to the same glacier that often share errors from similar instruments (e.g., ASTER) or estimation methods, and so we conservatively assume that their errors are fully correlated. Based on previous equations, we separate error propagation for elevation change, density conversion, and anomaly calculation:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{d}h,{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{k=\mathrm{1}}^N{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathrm{d}h}}_{k,\mathrm{PoR}}\cdot }{f}_{\mathit{\rho }},\text{(10)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{f_{\mathit{\rho }},{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{k=\mathrm{1}}^N{\mathit{\sigma }}_{f_{\mathit{\rho }}\cdot }{\stackrel{\mathrm{‾}}{\mathrm{d}h}}_{k,\mathrm{PoR}},\text{(11)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathit{\beta },{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{k=\mathrm{1}}^N{\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{\mathit{\beta }}}_{k,Y}}.\end{array}$$

The total uncertainty in the mean calibrated annual mass change time series for a certain glacier is

$$\begin{array}{}\text{(12)}& {{\mathit{\sigma }}^{\mathrm{2}}}_{{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}={\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{d}h,{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{f_{\mathit{\rho }},{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathit{\beta },{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}^{\mathrm{2}}.\end{array}$$

2.2.4 Integrating glacier mass changes into larger regions

Every glacier with available geodetic observations has a mean calibrated annual mass change (approximately 207 000 glaciers covering 96 % of the world's glaciated area). The remaining unobserved glaciers (g_unobs) are assumed to behave as the regional mean of the observed sample. The individual glacier mean calibrated annual mass change time series can be integrated into any larger region R containing multiple glaciers. The region-specific calibrated annual mass change time series (i.e., grid cell or glacier region) B_cal at year Y and in region R is calculated as the area-weighted mean of the individual calibrated annual mass balance time series of the sample of observed glaciers belonging to region R (or the specific grid point R for the gridded product in Fig. 6ii).

$$\begin{array}{}\text{(13)}& B_{\mathrm{cal},R,Y}={\displaystyle \frac{{\sum }_{g=\mathrm{1}}^N{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}\cdot A_g}{\left({\sum }_{g=\mathrm{1}}^N{A}_{g,Y}\right)}}\end{array}$$

To derive the uncertainty in the region-specific calibrated mass balance, we need to account for spatial correlations between the uncertainties of per-glacier mean calibrated annual mass change. Indeed, our three error sources, i.e., elevation change, density conversion, and mean annual mass balance anomaly calculation, are significantly correlated spatially. For elevation change, we use the spatial correlation in the elevation change error ρ_dh(d) estimated in Hugonnet et al. (2021), as it is the main data source in the FoG database. These spatially correlated elevation errors are largely due to instrument noise and temporal interpolation to match an exact period of estimation:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{{\mathit{\sigma }}^{\mathrm{2}}}_{\mathrm{d}h,B_{\mathrm{cal},R,Y}}& ={\displaystyle \frac{\mathrm{1}}{A_{\mathrm{tot}}}}{\sum }_{g_{\mathrm{1}}}{\sum }_{g_{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{d}h}\left(d_{g_{\mathrm{1}},g_{\mathrm{2}}}\right)\\ & \cdot {\mathit{\sigma }}_{\mathrm{d}h,{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g_{\mathrm{1}},Y}}\cdot {\mathit{\sigma }}_{\mathrm{d}h,{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g_{\mathrm{2}},Y}}\cdot A_{g_{\mathrm{1}}}\cdot A_{g_{\mathrm{2}}},\end{array}\end{array}$$

where $d_{g_{\mathrm{1}},g_{\mathrm{2}}}$ is the distance between glaciers, A_g is the area of glacier g, and $A_{\mathrm{tot}}={\sum }_g{A}_g$ is the regional glacier area.

For density conversion, we estimated a spatial correlation function of the uncertainty in the density conversion ${\mathit{\rho }}_{{\mathit{\sigma }}_{f_{\mathit{\rho }}}}\left(d\right)$ by performing a similar variogram analysis to that detailed for annual anomalies (Sect. 2.2.1) but instead applied to modeled estimates of the annual density of volume change for all glaciers globally. These estimates were obtained by pairing a mass balance model (Huss and Hock, 2015) with a firn densification model (Huss, 2013) calibrated on geodetic mass balances (Hugonnet et al., 2021). We find a spatial correlation function of

$$\begin{array}{}\text{(15)}& {\mathit{\rho }}_{{\mathit{\sigma }}_{f_{\mathit{\rho }}}}\left(d\right)=s_{\mathrm{1}}{e}^{-\frac{\mathrm{3}d}{r_{\mathrm{1}}}}+s_{\mathrm{2}}{e}^{-\frac{\mathrm{3}d}{r_{\mathrm{2}}}}\text{if}d>\mathrm{0},\text{else}\mathrm{1},\end{array}$$

where d is the distance between two glaciers, s₁=0.12 and s₂=0.72 are the partial sills, and r₁=200 and r₂=5000 km are the correlation ranges. These spatially correlated density errors are due to large local and regional variations in precipitation and firn densification, resulting in spatially correlated errors from the average value.

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{{\mathit{\sigma }}^{\mathrm{2}}}_{f_{\mathit{\rho }},B_{\mathrm{cal},R,Y}}& ={\displaystyle \frac{\mathrm{1}}{A_{\mathrm{tot}}}}{\sum }_{g_{\mathrm{1}}}{\sum }_{g_{\mathrm{2}}}{\mathit{\rho }}_{{\mathit{\sigma }}_{f_{\mathit{\rho }}}}\left(d_{g_{\mathrm{1}},g_{\mathrm{2}}}\right)\cdot {\mathit{\sigma }}_{f_{\mathit{\rho }},{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g_{\mathrm{1}},Y}}\\ & \cdot {\mathit{\sigma }}_{f_{\mathit{\rho }},{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g_{\mathrm{2}},Y}}\cdot A_{g_{\mathrm{1}}}\cdot A_{g_{\mathrm{2}}}\end{array}\end{array}$$

For the mean annual mass balance anomalies, we assume that errors in the real values are completely correlated at regional scales and thus propagated as

$$\begin{array}{}\text{(17)}& {{\mathit{\sigma }}^{\mathrm{2}}}_{\mathit{\beta },B_{\mathrm{cal},R,Y}}={\displaystyle \frac{\mathrm{1}}{A_{\mathrm{tot}}}}{\sum }_g{{\mathit{\sigma }}^{\mathrm{2}}}_{\mathit{\beta },{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}\cdot {A^{\mathrm{2}}}_g.\end{array}$$

Finally, following the assumption that correlation between sources has a negligible impact compared to the spatial correlation of errors within the same source, we combine all sources of error propagated at the regional scale as independent:

$$\begin{array}{}\text{(18)}& {{\mathit{\sigma }}^{\mathrm{2}}}_{B_{\mathrm{cal},R,Y}}={{\mathit{\sigma }}^{\mathrm{2}}}_{\mathrm{d}h,B_{\mathrm{cal},R,Y}}+{{\mathit{\sigma }}^{\mathrm{2}}}_{f_{\mathit{\rho }},B_{\mathrm{cal},R,Y}}+{{\mathit{\sigma }}^{\mathrm{2}}}_{\mathit{\beta },B_{\mathrm{cal},R,Y}}.\end{array}$$

The regional mass change (Gt) of water is then obtained by multiplying the specific mass change by the region's (or grid point's) glacierized area S_R, corrected to the year 2000 using the area change rates updated from Zemp et al. (2019). We propagate the uncertainty in the specific regional mass change, the uncertainty in the regional area (Paul et al., 2015), and the uncertainty in the area change, assuming them to be uncorrelated. Errors in the area stem mostly from remote sensing delineation errors, while errors in area change stem from a lack of multi-temporal outlines to constrain area change. They are largely uncorrelated with the error sources described above for elevation change, glaciological measurements, and anomalies. However, elevation change estimates usually already consider errors in area at the scale of each glacier, so we might conservatively be double-counting these:

$$\begin{array}{}\text{(19)}& {\displaystyle }\mathrm{\Delta }{M}_{R,Y}=B_{R,Y}\cdot \left(S_R+\mathrm{\Delta }{S}_{R,Y}\right),\text{(20)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{\Delta }{M}_{R,Y}}=\left|\mathrm{\Delta }{M}_{R,Y}\right|\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{B_{R,Y}}}{B_{R,Y}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{S_R}}{S_R}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{\Delta }{S}_{R,Y}}}{\mathrm{\Delta }{S}_{R,Y}}}\right)}^{\mathrm{2}}},\end{array}$$

where $\frac{{\mathit{\sigma }}_S}{S}=\mathrm{5}\mathit{\%}$ (Paul et al., 2015), and $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }{S}_R}}{\mathrm{\Delta }{S}_R}$ is updated from Zemp et al. (2019).

The global annual (Y) and cumulative mass changes (Gt) and the sea level equivalent for any given period of record are finally calculated as the sum of the regional mass changes, assuming that the regional mass loss uncertainties are independent and uncorrelated between every region. To simplify the combination of annual values into long-term trends or cumulative annual values, we assume the yearly uncertainty to be independent of other years. This is true for glaciological measurement, with an independent uncertainty estimation for each individual year of the time series but not for the elevation change measurements, where uncertainties are correlated over the years of the survey period.

$$\begin{array}{}\text{(21)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{M}_{\mathrm{Glob},Y}& ={\sum }_{R=\mathrm{1}}^{\mathrm{19}}\mathrm{\Delta }{M}_{R,Y}\text{and}\mathrm{\Delta }{M}_{\mathrm{Glob},\mathrm{PoR}}\\ & ={\sum }_{R=\mathrm{1}}^{\mathrm{19}}\mathrm{\Delta }{M}_{R,\mathrm{PoR}}\end{array}\text{(22)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{\mathrm{\Delta }{M}_{\mathrm{Glob},Y}}& =\sqrt{{\sum }_{R=\mathrm{1}}^{\mathrm{19}}{\left({\mathit{\sigma }}_{\mathrm{\Delta }{M}_{R,Y}}\right)}^{\mathrm{2}}}\text{and}{\mathit{\sigma }}_{\mathrm{\Delta }{M}_{\mathrm{Glob},\mathrm{PoR}}}\\ & =\sqrt{{\sum }_{R=\mathrm{1}}^{\mathrm{19}}{\left({\mathit{\sigma }}_{\mathrm{\Delta }{M}_{R,\mathrm{PoR}}}\right)}^{\mathrm{2}}}\end{array}\text{(23)}& {\displaystyle }\begin{array}{rl}{\mathrm{SLE}}_Y& ={\displaystyle \frac{\mathrm{\Delta }{M}_{\mathrm{Glob},Y}}{S_{\mathrm{ocean}}}}\cdot {\mathrm{10}}^{\mathrm{6}}\pm {\mathit{\sigma }}_{{\mathrm{SLE}}_Y}\text{and}{\text{SLE}}_{\mathrm{PoR}}\\ & ={\displaystyle \frac{\mathrm{\Delta }{M}_{\mathrm{Glob},\mathrm{PoR}}}{S_{\mathrm{ocean}}}}\cdot {\mathrm{10}}^{\mathrm{6}}\pm {\mathit{\sigma }}_{{\mathrm{SLE}}_{\mathrm{PoR}}}\end{array}\text{(24)}& {\displaystyle }\begin{array}{rl}& {\mathit{\sigma }}_{{\mathrm{SLE}}_Y}=\sqrt{{\mathit{\sigma }}_{\mathrm{\Delta }{M}_{\mathrm{Glob},Y}}^{\mathrm{2}}+{\mathit{\sigma }}_{S_{\mathrm{ocean}}}^{\mathrm{2}}}\text{and }\\ & {\mathit{\sigma }}_{{\mathrm{SLE}}_{\mathrm{PoR}}}=\sqrt{{\mathit{\sigma }}_{\mathrm{\Delta }{M}_{\mathrm{Glob},\mathrm{PoR}}}^{\mathrm{2}}+{\mathit{\sigma }}_{S_{\mathrm{ocean}}}^{\mathrm{2}}}\end{array}\end{array}$$

$S_{\mathrm{ocean}}=\mathrm{362.5}\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{km}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{S_{\mathrm{ocean}}}=\mathrm{0.1}\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{km}}^{\mathrm{2}}$ (Cogley, 2012).

One strength of producing per-glacier mass change time series is the possibility of integrating them as an area-weighted mean at any given spatial resolution (i.e., regular grids, subregions, regions, and basins). In this study we integrate glacier mass changes at three spatial resolutions: regionally by the 19 RGI first-order regions and globally to allow direct comparison with previous global observation-based assessments by Zemp et al. (2019) and Hugonnet et al. (2021). Further, taking advantage of the per-glacier annual time series, we generate a global gridded product of annual glacier mass changes for the Copernicus Climate Change Service (C3S, https://climate.copernicus.eu/, last access: May 2025) Climate Data Store (CDS). For consistency with other climate observation datasets (e.g., C3S), we provide glacier changes on a global regular grid of 0.5° latitude and longitude. For temporal consistency within all the regions, we extend the global time series only as far as the hydrological year 1976, in contrast to Zemp et al. (2019), who reached back to 1962. This adjustment is due to the absence of annual observations in the Southern Hemisphere regions prior to 1976 (as evidenced in Zemp et al., 2019, Fig. 10). Regional time series start from the date of the first year of mass change records available for the region (see Table 5). Importantly, our fully operational approach allows yearly updates to be produced as soon as new glacier observations are incorporated into the FoG database of the WGMS.

2.2.5 Methodological progress in data fusion of glaciological and geodetic data

The specific methodological improvements in data fusion of glaciological and geodetic data of the present assessment with respect to Zemp et al. (2019) are detailed in Table 3.

Table 3 Specific methodological improvements in data fusion of glaciological and geodetic data with respect to Zemp et al. (2019).

2.2.6 Description of the datasets

The datasets are described in Table 4. The main dataset, Dataset 1, corresponds to individual glacier annual mass change time series provided in .csv files by RGI first-order regions. Glaciers are identified by their RGI-Id, centroid latitude, and centroid longitude corresponding to the RGI60 glacier outline geometry. For Dataset 1, the start of the time series is region-dependent and corresponds to the date of the first year of mass change records available for the region. We also provide for every region and on a glacier-by-glacier basis a .csv file with additional metadata information and a .readme file. The second dataset, Dataset 2, is a byproduct of Dataset 1: it corresponds to an integration of the individual glacier time series on a global grid of 0.5° latitude and longitude. We chose this resolution and a netCDF file format to make the glacier change product consistent and easily usable by other climate observation datasets (e.g., C3S). To ensure global completeness of annual glacier mass changes, this dataset spans the period from 1976 to 2024.

Table 4 Details of the annual glacier mass change output datasets. n/a – not applicable

3 Results

3.1 Global glacier mass changes

Our results provide revised global annual glacier mass changes extending back to the hydrological year 1976 (annual regional glacier mass changes further back in time, depending on the region) at various spatial levels: per-glacier (Fig. 6i), regional, and on a global scale (Figs. 4, 5, and 6; Table 5). Globally, glaciers have lost 9179 ± 621 Gt of water (or 187 ± 20 Gt yr⁻¹), contributing 25.3 ± 1.7 mm of the sea level rise since 1976. About 41 % of the total loss, equivalent to 10 mm of the sea level rise, occurred during the last decade from 2015 to 2024 only (Table 5). The strong acceleration of global glacier mass loss is evident over the last years of observations, where 5 (2019, 2020, 2022, 2023, and 2024) out of the last 6 years present the strongest global glacier mass loss ever recorded, exceeding 430 Gt yr⁻¹ (or 1.2 mm yr⁻¹ of the sea level rise contribution), and the record year 2023 at 540 ± 69 Gt of water. During 2023 alone, the glacier mass loss was about 80 Gt higher than any other year on record, corresponding to 6 % of the total loss since 1975/1976. In only this 1 year, glacier melt raised sea levels by 1.5 ± 0.2 mm.

Figure 4 Specific annual mass change time series for the 19 regions, with Southern Andes separated by second-order regions and with the respective uncertainties. The sizes of the circles in the pie charts represent the mass lost (Gt) by region and the globe since the hydrological year 1976. Colored sections represent the mass lost during the last decade only (2015–2024). The heat maps represent the regional and global glacier mass changes (m w.e.) for every hydrological year over the common period from 1976 to 2024. Global results from 1976 to 2024 are represented as cumulative mass changes (Gt), and the global heat map is in the range ±500 Gt. Regional time series start from the date of the first year of available mass change records. The grey sections show the periods where annual mass balance anomalies from neighboring regions are used to capture the temporal variability.

Figure 5 Cumulative regional glacier mass changes from the hydrological year 1976 to 2024 for the 19 regions. Specific mass changes (m w.e.) indicate the mean height of the water layer lost over a given glacier surface, and large negative values suggest regions where glaciers have suffered the most. By multiplying by the regional glacier area (km²), we obtain the mass change (Gt). Cumulative glacier mass changes (Gt) correspond to the volume of water lost (1 km³ w.e. = 1 Gt) and are related to the regional contributions to global mean sea level rise (mm).

Table 5 Annual rates of regional glacier mass changes (Gt and m w.e.) from 1976 to 2024 and the last decade from 2015 to 2024. Mean area is calculated from the annual changes in area estimated with change rates updated from Zemp et al. (2019).

Regionally, glacier mass losses in 2023 range from −0.3 m w.e. in the less impacted regions to −2.6 ± 0.6 m w.e. in Western North America. Strong glacier mass loss rates (lower than −1.0 m w.e.) were also reported in Central Europe, Alaska, New Zealand, Southern Andes, Iceland, and Scandinavia. In 2024, it was Scandinavian (−2.3 ± 0.4 m w.e.) and Svalbard (−1.6 ± 0.3 m w.e.) glaciers that suffered the most, followed by Russian Arctic, Western North America, New Zealand, Southern Andes, South Asia East, Low Latitudes, Central Europe, and Asia North, with strong mass loss rates lower than −1.0 m w.e. The noteworthy years 2022, 2023, and 2024 are the first – and consecutive – years where all 19 glacier regions experienced average mass loss.

3.2 Regional glacier mass changes

Since the year 2000, all glacier regions have lost ice (Fig. 5). Alaska, Western Canada US, Svalbard, Russian Arctic, North Asia, Caucasus, Central Asia, Asia East, Southern Andes, New Zealand, and Antarctica have experienced increased mass loss during the last decade (2015–2024) compared to the full period since 1976. Extremely negative regional decadal rates reaching lower than −1.2 m w.e. yr⁻¹ are observed in Alaska, Western North America, and Central Europe during the last decade. Record-breaking annual rates (more negative than −2.5 m w.e.) occurred in Western North America (2023), Iceland (2010), Central Andes (1999, 2018), and Central Europe, the latter experiencing in 2022 the most negative regional mass balance ever recorded of −3.1 m w.e. Mass losses larger than 100 Gt of water in a single year are present only in Antarctica (2023) and Alaska (12 years in total, 5 of them during the last decade), with two regional yearly record mass losses of 176 ± 52 and 163 ± 48 Gt in 2004 and 2019, respectively.

4 Discussion

4.1 Multi-spatial dimensions of annual mass change series

The revised dataset claims several strengths primarily related to the enhanced temporal resolution, providing glacier changes at annual resolution and the multi-spatial dimensions for data integration from individual glacier to global. We remind the reader here that, by construction, nearby glaciers share a large fraction of the variance in mass balance variability and are thus not independent. In a visualization example for selected years in Iceland, Fig. 6 depicts the multiple spatial dimensions from the dataset: individual glacier annual time series (Fig. 6i) and annual time series aggregated into any larger-scale region encompassing multiple glaciers, such as regular grid cells of 0.5° latitude and longitude (e.g., Fig. 6ii or any user-specified resolution) and regional to global specific mass changes (Figs. 6iii and 4) and mass changes (Figs. 6iv and 10). This versatility enables identification of individual years marked by significant glacier changes and detection of zones with varying impacts. For instance, it allows us to pinpoint regions and subregions that were affected by specific annual climate variations (e.g., droughts, floods, and heat waves) as well as those with a larger or smaller influence on the yearly contribution to hydrology and annual sea level rise.

Figure 6 Illustration of the multi-spatial dimensions of the global annual mass change series, with an example for Iceland during selected years. (i) Individual glacier annual mass change series, (ii) gridded annual mass change series, (iii) regional specific annual mass change (m w.e.), and (vi) regional mass change (Gt; note that all the previous dimensions i, ii, and iii may also be represented in gigatons).

The represented years in Fig. 6 are chosen arbitrarily: the initial hydrological year considered in the global assessment, the well-known extremely negative mass change hydrological year 2010 attributed to the eruption of the Eyjafjallajökull volcano (Aðalgeirsdóttir et al., 2020; Belart et al., 2020; Möller et al., 2019), the most positive year of the time series in 2015, and the global extreme mass loss year 2023. This example allows us to illustrate the spatially rather homogeneous glacier mass loss of the year 2010, with larger Icelandic glaciers all losing more than 2 m w.e. and some smaller glaciers losing between 1.5 and 2 m w.e., whereas other years show larger variabilities between glaciers and grid points. This example demonstrates the richness of the dataset for interpretation of glacier mass changes at different spatial scales and a deeper analysis of the spatial and temporal impacts of known glaciological trends and anomalies, e.g., the Andes mega-drought (Dussaillant et al., 2019; Garreaud et al., 2017, 2020; Gillett et al., 2006) or the Karakoram anomaly (Farinotti et al., 2020; Gao et al., 2020; Ougahi et al., 2022), at an unprecedented yearly temporal resolution. We note that the annual mass balance anomalies are extracted from a handful of glaciers in each region, and thus, in each region, individual glaciers share a large fraction of these variabilities. This limitation is, however, evidenced well by large uncertainties in undersampled regions and periods.

4.2 Leave-one-out cross-validation

Due to the lack of independent measurements available to compare and validate our glacier change assessment, we applied a leave-one-out cross-validation exercise over selected reference and benchmark glaciers. Reference and benchmark glaciers are selected by considering their fluctuations to be mainly driven by climatic factors. They provide a reliable and well-documented sample of globally distributed long-term observation series, with more than 10 (benchmark) and 30 (reference) years of continuous and ongoing glaciological mass balance measurements. Notably, glaciological time series can be subject to biases inherent to the glaciological method (e.g., Thibert et al., 2008) and are encouraged to be periodically reanalyzed and calibrated with long-term trends derived from high-resolution elevation change surveys (Zemp et al., 2013). To reduce the risk of validating over potentially erroneous “truths”, we do not use all available glaciological time series in this experiment. We select a sample of 74 reference and benchmark glaciers for the leave-one-out cross-validation and then repeat the analysis over a selection of 32 glaciers knowingly reanalyzed.

For each selected glacier, we compare the original “reference” mass balance time series (reference Ba) as available from the FoG database with the estimated leave-one-out mean calibrated annual mass change time series (leave-one-out Ba). The latter is obtained as described in the original methodology by calibrating the mean annual mass balance anomaly of the glacier over its geodetic sample, but this time we exclude the selected glacier anomaly from the processing. Reference and benchmark glaciers are usually highly monitored and contain multiple sources of geodetic observations for different time periods. However, more than 80 % of the world's glaciers present only one source of geodetic observations, i.e., the 20-year elevation change rates from Hugonnet et al. (2021). To make our validation exercise relevant for these undersampled glaciers, we only consider for calibration the elevation change rates from Hugonnet et al. (2021), excluding other sources of geodetic observations. For each glacier, the mean error (ME) and the standard deviation of the residual (S) between the reference and leave-one-out Ba are estimated as metrics to quantify potential systematic errors and the magnitudes of random errors, respectively.

Validation results are shown for all 74 selected reference and benchmark glaciers in Fig. 7. As verified by the low ME of 0.05 m w.e. for annual values and 0.06 m w.e. for long-term trends (Fig. 7a and b), we find no systematic error on the estimated leave-one-out Ba. This means random errors are the largest error source, at 0.51 m w.e. for the annual values and 0.26 m w.e. for long-term trends. Differences in annual values and the long-term trends may come from the different geodetic datasets used for reanalysis. Reference and benchmark glaciers, if reanalyzed, use high-quality local elevation change observations for calibration, whereas here our estimated leave-one-out Ba is calibrated over the Hugonnet et al. (2021) elevation change only.

Figure 7 Leave-one-out cross-validation results and statistics over 74 glaciological time series from reference and benchmark glaciers. Mean error (ME) and standard deviation of the residuals (S) between the estimated leave-one-out Ba cross-validation time series and the 74 reference and benchmark glacier reference Ba time series. (a) Annual values from the estimated leave-one-out Ba against the reference and benchmark glacier reference Ba. Every dot corresponds to a yearly observation. (b) Long-term trends (1976–2024) from the estimated leave-one-out Ba time series against long-term trends for the reference and benchmark glacier reference Ba time series. (c) Amplitude of the annual variability measured as the time series variability SD (not to be confused with the standard deviation of the residuals denoted as S) for the period 1976–2024 from the estimated leave-one-out Ba time series against the reference and benchmark glacier reference Ba time series. (d) ME and (e) S of residuals for each reference and benchmark glacier against the estimated uncertainty of the mean calibrated annual mass change estimate for the same glacier at 1σ. In panels (a), (b), (c), and (d), each value corresponds to 1 of the 74 reference and benchmark glaciers used for cross-validation, and the symbols correspond to the glacier regions to which they belong. The size of the symbol is related to the area of the glacier.

There is a slight underestimation of the glacier annual mass balance variability amplitude, as shown by the slightly negative ME of −0.08 m w.e. between the leave-one-out Ba and the reference Ba time series amplitude. This bias may come from a smoothing of the mean annual mass balance anomaly extreme values when averaging the nearby glaciological time series. This effect is clear for, e.g., the Afoltbreen glacier (Fig. 8g), where the extreme years in the reference mass balance series are less extreme in the leave-one-out series.

Figure 8 Selected examples of the leave-one-out cross-validation results on individual reference glaciers (a–l, more than 30 years of observations) and benchmark glaciers (m–v, more than 10 years of observations). Reference and benchmark glacier reference annual mass balance time series (grey) and the leave-one-out estimated mean calibrated annual mass change time series and uncertainties at the 2σ (95 % confidence interval, light blue) (left plots). Annual values from the estimated leave-one-out Ba against the reference and benchmark glacier reference Ba (right plot). Every dot corresponds to a yearly observation.

In general, regions with a rich glaciological sample, like Central Europe, Scandinavia, Svalbard, Iceland, and Arctic Canada North, perform well. Specifically, Hintereisferner located in the well-sampled Central Europe presents a residual S as low as 0.19 m w.e. and a ME of −0.09 m w.e. The largest random errors are observed in glaciers Brewster in New Zealand (S: 0.94, ME: 0.32 m w.e.), Gulkana (S: 0.82, ME: −0.01 m w.e.), Wolverine in Alaska (S: 1.03, ME: 0.16 m w.e.), and Zongo at low latitudes (S: 0.83, ME: 0.17 m w.e.). The largest systematic errors are observed in glaciers Columbia in North America (S: 0.37, ME: 0.66 m w.e.), Tungnaarjokull in Iceland (S: 0.35, ME: 0.63 m w.e.), Mittivakkat in Greenland (S: 0.54, ME: −0.61 m w.e.), and Sarenne in Central Europe (S: 0.54, ME: −0.95 m w.e.) (Figs. 7 and 8 for selected reference and benchmark glaciers). Still, large differences may occur between reference and predicted values for individual years (outliers in Fig. 7a), also evident in the few annual values where the reference annual mass balance is not included in the uncertainties of the leave-one-out annual estimates (e.g., Mittivakkat, Fig. 8j, and Djankuat, Fig. 8n).

Figure 7d and e show the ME and S of residuals for each reference and benchmark glacier against the uncertainties of the mean calibrated annual mass change time series ${\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{B}}_{\mathrm{cal},g,Y}}$ calculated using our kriging predicting method for the same sample of glaciers. In most cases, the predicted error is larger than both the leave-one-out cross-validation ME and S residuals for all of the glaciers. This means that the kriging method is cautious, predicting larger uncertainties than what is observed in the residuals, thus confirming our uncertainties to be within a conservative range with no significant systematic errors and able to capture random errors introduced by the prediction method well.

The leave-one-out cross-validation shows that our method introduces only small systematic errors and has a per-glacier random error that is conservative. However, most reference and benchmark glaciers are located in regions with a high density of mass balance observations able to compensate for their exclusion. Due to this, our leave-one-out exercise does not necessarily represent a typical glacier in a data-scarce region. To explore the validity of our methodology in undersampled regions, where the mean glacier anomaly may come from distant glaciological time series, below we present a leave-block-out experiment over the full sample of reference and benchmark glaciers.

4.3 Leave-block-out cross-validation

The experiment is similar to the leave-one-out cross-validation, but this time we remove from the processing not one but all of the glacier anomalies surrounding the reference and benchmark glacier for increasing distances: first, we consider all anomalies further than 1 km (i.e., as in leave-one-out cross-validation, but the glacier's mass balance series is excluded), and then we remove the closest glacier's annual mass balance anomalies at different distance thresholds of 60, 120, 250, 500, 800, and 1000 km. At every step a new mean glacier mass balance anomaly is calculated for the reference and benchmark glacier from the evolving sample and calibrated to the 20-year elevation change rates from Hugonnet et al. (2021). The ME and S of the residuals between the estimated leave-block-out mean calibrated annual mass change time series (leave-block-out Ba) and the reference annual mass balance (reference Ba) are estimated at every distance threshold.

Considering all reference Ba values against the leave-block-out Ba for the six different distances, systematic errors (ME) appear to stay stable between 0.06 and 0.08 m w.e. until 500 km and then increase to 0.20 m w.e. above 500 km. Random errors (S) appear to increase gradually as the distance gets larger (Table 6).

Table 6 Leave-block-out cross-validation ME and S residuals for the yearly values (i.e., Fig. 7a) and for the amplitude of the annual variability during the 1976–2024 period (i.e., Fig. 7c) of the estimated leave-block-out Ba calculated from samples at different distance thresholds.

Figure 9 illustrates the leave-block-out glacier-wide residual results as a function of the distance to the closest glacier annual mass balance anomaly considered. There is no apparent influence of the distance on systematic errors in the estimated glacier-wide leave-block-out Ba justified by the absence of trends in Fig. 9a. In these cases, the slight systematic errors will mostly depend on whether the reference series are reanalyzed or not and the quality of the elevation change used for calibration. As expected, random errors (residual S) increase as the glacier mean annual mass balance anomaly is calculated from a more distant sample (Fig. 9b), from 0.5 m w.e. for nearby time series to 1 m w.e. for series located farther than 2000 km. Importantly, in most cases systematic errors are not significant given our sample size, and random errors are captured well by the mean calibrated annual mass change uncertainty at 1σ, independent of the distance of the sample (Fig. 9c and d). This means that our predicted uncertainties reflect the true variability in the residuals and that our model provides realistic confidence intervals for the mean annual mass balance anomaly predictions. On average, S starts to become larger than σ with distances to the closest glacier larger than 500 km, but the large spread suggests that this comes from the randomness of the predictions.

Figure 9 Leave-block-out cross-validation results and statistics over 74 glaciological time series from reference and benchmark glaciers. (a) ME and (b) S as a function of the distance to the closest glacier anomaly used to calculate the mean mass balance anomaly. (c) ME of residuals and (d) S of residuals at different distance thresholds against the estimated uncertainty of the mean calibrated annual mass change estimate for the same glacier at σ. Symbols correspond to the region where the benchmark or reference glacier belongs.

We can conclude that the spatial extrapolation of distant glaciers does not introduce clear systematic errors but increases random errors. However, these increased random errors are predicted well by our uncertainty assessment, showing larger uncertainties over glaciers in undersampled regions. Both the leave-one-out and leave-block-out cross-validations show that our algorithm can capture the annual variability of individual glacier mass changes on glaciers not presenting glaciological time series (99 % of the global glaciers) with realistic and conservative uncertainties.

4.4 Improvements with respect to earlier assessments

A general overview of the improvements with respect to previous assessments by Zemp et al. (2019) and Hugonnet et al. (2021) is provided in Table 7. Regional and global glacier mass change results for the three observation-based estimates are compared in Fig. 10 (Gt).

Table 7 General methodological ameliorations with respect to Zemp et al. (2019) and Hugonnet et al. (2021). n/a – not applicable

Figure 10 Annual glacier mass change (Gt) from this study compared with results from Zemp et al. (2019) and Hugonnet et al. (2021).

Broadly, our new approach effectively corrects the negative bias in the long-term trends observed in Zemp et al. (2019) thanks to the integration of the glacier elevation changes from Hugonnet et al. (2021), which translates into a significant reduction in both regional and global uncertainties that is largely noticeable for the more recent years (Fig. 10, Table 8). Globally, glacier mass change rates between 1976 and 2016 are less negative (−148 ± 20 Gt yr⁻¹) than previously estimated by Zemp et al. (2019) (−204 ± 45 Gt yr⁻¹, Table 7). The resulting 6069 ± 783 Gt of the cumulative mean mass loss for the period is 35 % lower than the 9290 ± 7698 Gt loss predicted by Zemp et al. (2019) for the same period. The differences mainly come from Southern Andes, Alaska, Russian Arctic, Antarctic and Subantarctic, and Greenland Periphery (Table 8) – all regions with limited geodetic coverage in the previous assessment. Both regionally and globally, the years after 2000 are aligned well with the Hugonnet et al. (2021) trends as a consequence of the calibration to their geodetic trends with global coverage.

Table 8 Annual rates of regional glacier mass change for the period 1976–2016 from this study compared with results from Zemp et al. (2019) for the same period and annual rates of regional glacier mass change for the period 2000–2019 from this study compared with results from Hugonnet et al. (2021) for the same period. Mean area is calculated from the annual changes in area estimated with change rates updated from Zemp et al. (2019).

Alaska, Arctic Canada North, Western North America, Russian Arctic, Caucasus, Low Latitudes, New Zealand, and Southern Andes exhibit less negative trends compared to Zemp et al. (2019). In contrast, Asia North and Arctic Canada South trends are more negative in the past. Overall, regional trends agree well with Hugonnet et al. (2021) trends during the overlapping period 2000–2019 (Table 8). Deviations of more than 5 Gt yr⁻¹ are found in Alaska, Greenland Periphery, and Antarctic and Subantarctic.

Differences in Greenland Periphery, Alaska, and Russian Arctic can be explained by the removal of unpublished geodetic observations from our processing chain, which may have biased the trends in Zemp et al. (2019). The new mean calibrated annual mass change time series in these regions rely only on the 2000–2019 Hugonnet et al. (2021) estimates as a calibration reference and the 1985–2000 Huber et al. (2020) estimates for Greenland Periphery's past period. The Antarctic and Subantarctic islands are a region prone to large uncertainties in all studies. Because all estimates agree within (large) uncertainties, no difference can be seen. The Antarctic and Subantarctic region presents no glaciological or geodetic measurements before the year 2000. The signal for the annual variability before 2000 is driven purely by the very distant Echaurren Norte (Central Andes) normalized time series, and past trends are only calibrated over the 2000–2019 Hugonnet et al. (2021) series, which very likely bias the period before 2000 towards more positive values. There is insufficient evidence to support the glacier mass gain observed before 2000 in both our assessment and Zemp et al. (2019). We assume results in the Antarctic and Subantarctic region to very likely be biased by the lack of observations and therefore highly uncertain, as reflected in our large error bars. However, we still include them to provide global glacier mass changes back to 1976.

Most regions display increased interannual variabilities when compared to both previous studies. The large uncertainties associated with the ASTER DEM time series in Hugonnet et al. (2021) prevent the detection of variability at short timescales (typically less than 3 years, Fig. 10). Similarly, in Zemp et al. (2019), the variance decomposition model (Eckert et al., 2011; Krzywinski and Altman, 2014) employed to extract the temporal mass change variability for each region has been shown to contribute to a slight smoothing of the annual amplitude signal (Zemp et al., 2020). Our approach allows us to better represent the interannual variability at the individual glacier level, supported by the leave-one-out cross-validation exercise, and the effect is inherited by the regional and global levels.

Importantly, our assessment reduces global and regional uncertainties compared to Zemp et al. (2019) (Table 8). This reduction is achieved firstly because our approach benefits from the almost complete observational sample of Hugonnet et al. (2021), which greatly reduces uncertainties after the year 2000. Secondly, there are multiple algorithm improvements, such as the spatial search for glacier anomalies and the normalization of neighboring glacier time series with original regional amplitudes. Notably, considering error propagation using empirical functions of spatial correlation for the different sources of error provides more realistic uncertainties than previous work. All these enhancements are further supported by the leave-one-out cross-validation over reference glaciers, confirming that our uncertainties still remain on the conservative side. Compared with Hugonnet et al. (2021), long-term trends agree within uncertainties, but our estimated errors are consistently larger throughout all the regions. This comes from the propagation of the mean mass change anomaly uncertainties, which is the largest source of error in our assessment and a source that was not considered in the Hugonnet et al. (2021) annual results. These differences are likely explained by different assumptions made in the temporal propagation of uncertainties rather than actual errors arising from estimations.

Like Zemp et al. (2019), one of the largest strengths of our method is the ability to provide glacier mass changes back in time to 1976 at a global level and even further back in regions with longer observational records (e.g., 1915 for Central Europe). Thanks to geostatistical modeling for temporal downscaling, we can now capture the annual temporal variability of glacier changes for each of the world's glaciers and then aggregate them into revised regional and global estimates of annual glacier mass changes.

4.5 Known limitations

4.5.1 Scarcity of glaciological in situ observations

The scarcity of glaciological data is the primary limitation in assessing the variability of glacier changes using our methodology. In sparsely observed regions like High Mountain Asia, Southern Andes, Arctic Canada South, Russian Arctic, Greenland, and Antarctica, annual variations in glacier changes depend on limited and distant or neighboring regions' time series, which may not necessarily be representative of the local annual glacier variability. As a consequence, the annual glacier mass change time series exhibit high uncertainties in these regions that are estimated realistically using our method.

We note a significant observational gap in the Southern Hemisphere where the glacier mass change variability before 2000 is driven by the single and very distant Echaurren Norte glacier. This results in a disproportionate importance at the global scale. The best way forward while still maintaining the independence and observation-based nature of the present assessment – crucial for calibration and validation of glacier models – is to bolster in situ glacier monitoring programs in these regions. In the short term, efforts must be directed towards ensuring the continuity of glacier in situ monitoring in the Southern Hemisphere and possible correction of past long-term series.

4.5.2 Availability of past geodetic observations

The lack of geodetic observations for the period before 2000 is consistent for most glacier regions and critical for accurate results of our assessment in less sampled regions, as shown for the Antarctic and Subantarctic islands. The best way to correct possible deviations in past time series is to calibrate them against accurate long-term geodetic glacier elevation changes. Geodetic observations can be temporally enriched in all regions by unlocking historical United States spy satellite archives (e.g., KH-9 Hexagon and Corona declassified satellite imagery) and national historical airborne image archives.

4.5.3 Grid point artifact in polar regions

To assess glacier mass changes, it is crucial to treat a glacier as a single and indivisible entity. Therefore, in order to properly integrate glacier-wide changes into a regular grid system, we consider a glacier to belong to a grid point if its centroid falls within that grid point's boundaries. If the grid cell is sufficiently large, it will encompass multiple glaciers at their full extension within the grid cell and the grid point mean mass change will be determined accordingly. However, in cases where the grid cell is smaller than the glacier's surface area, the grid point containing the glacier's centroid will represent the mass change of the entire glacier, despite not all of its extension being contained within the grid point (Fig. 11a). This discrepancy is particularly evident in polar regions above 60° latitude when integrating mass changes at a 0.5° global grid resolution. Polar grid cells are relatively smaller in area compared to large polar glaciers. Consequently, this leads to a biased estimate of mass change at the grid point containing the glacier centroid, and consequent neighboring glacierized grid points lack a mass change estimate (Fig. 11b).

Figure 11 Mockup example of the grid point artifact in polar regions. (a) Glacier outlines (white) and their centroids (blue points) in the region Arctic Canada North in 1° grid cells. (b) Visual illustration of the 0.5° gridded annual mass change integration.

This issue might be especially critical for deconvolving the glacier signal for gravimetry (GRACE, e.g., Blazquez et al., 2018; Chen et al., 2022) or other applications in polar regions due to coarse resolution of the ancillary datasets (usually not smaller than 0.5°). A potential solution for larger-scale applications with coarser spatial resolution would be an area weight per tile glacier area, but this would bring an additional bias related to the divisibility of the glacier signal.

4.5.4 Calendar years versus hydrological years

Our results present regional glacier mass changes spanning the hydrological years from 1976 to 2024. In glaciological terms, it is widely accepted that the hydrological year starts in winter with the onset of the accumulation season and concludes at the end of the summer or ablation season (Cogley et al., 2011). Consequently, the hydrological year varies across regions (Southern Hemisphere, Northern Hemisphere, and the tropics) and does not align with the calendar year. Gridded annual glacier mass change values for a given hydrological year will not be fully consistent. For grid points located in the Northern Hemisphere, glacier mass changes correspond to the period from 1 October of the previous year to 30 September of the given year, while grid points in the Southern Hemisphere will represent glacier mass changes from 1 April of the previous year to 31 March of the given year. It is important to note that this discrepancy, stemming from the input data, introduces inconsistencies and uncertainties into the gridded global assessments that users should acknowledge. For cumulative values over longer periods, these differences become less significant. Addressing this issue would need an increase in temporal resolution of the input data to monthly observations, which is not feasible at the global level when relying solely on observations.

4.5.5 Glacier-specific area change rates

We took into account the impact of glacier area changes in time for our regional mass change estimates, as done by Zemp et al. (2019). In the former study, the evolution of area change rates is calculated for each first-order glacier region independently represented in two time steps: a past period with no change rates and a linear change rate calculated from a regional sample of glaciers with observations. This assumption is strong since glacier area changes are far from being linear. This is still the best possible guess considering the available observations. In this study, we further assume that the glacier-specific area change rates behave like their regional mean change rates. This may introduce an additional bias since glacier-specific area change rates greatly depend on the size of the considered glaciers, with an observed decreasing mean and increasing variability of relative area changes towards smaller glaciers (Fischer et al., 2014; Paul et al., 2004). This also implies that the average change obtained for a greater region ultimately depends on the glacier size distribution considered in a specific sample, which may or may not be representative of the full regional glacier size range present in the RGI60 inventory used here.

5 Data availability

The annual mass change time series for individual glaciers and the derived global gridded annual mass change product at a spatial resolution of 0.5° latitude and longitude are available from the World Glacier Monitoring Service (Dussaillant et al., 2025; https://doi.org/10.5904/wgms-amce-2025-02). Earlier versions of the gridded product are available from the Copernicus Climate Data Store (CDS) web-based service (Dussaillant et al., 2024; https://doi.org/10.24381/cds.ba597449). The present version will follow in the next C3S phase. As input data we used the FoG database from the World Glacier Monitoring Service (WGMS, 2025; https://doi.org/10.5904/wgms-fog-2025-02) and RGI 6.0 from the National Snow and Ice Data Center (NSIDC, RGI consortium, 2017; https://doi.org/10.7265/4m1f-gd79).

6 Code availability

The code for implementing our approach is available on GitHub (https://github.com/idussa/global_mb_fusion, last access: May 2025) and on Zenodo (https://doi.org/10.5281/zenodo.15344003; Dussaillant and Hugonnet, 2025).

7 Conclusions

Building upon the insights from previous assessments, we present a new dataset of glacier mass changes with global coverage and annual resolution from 1976 to 2024 based on an approach that combines the strengths from glaciological in situ measurements and geodetic satellite observations. Our results offer critical insights into the acceleration of glacier melt at a global scale during recent years. While glaciers featured moderate mass loss rates in the 1970s and 1980s, the rates strongly increased in the 1990s, reaching about 375 ± 16 Gt per year in the last decade (2015–2024) and a record annual loss of 540 ± 69 Gt alone in 2023. Since 1976, glaciers globally have lost 9179 ± 621 Gt (187 ± 20 Gt per year) of water, contributing 25.3 ± 1.7 mm (0.5 ± 0.2 mm per year) to the global mean sea level rise. Nearly half (41 %) of this loss, equivalent to 10 ± 0.2 mm (1.0 ± 0.2 mm per year) of sea level rise, occurred in the last decade (2015–2024), with 6 % (1.5 mm) occurring in the global record year 2023 alone.

Compared to earlier assessments, our approach allows us to combine the global coverage from geodetic satellite observations with the annual variability and temporal coverage from glaciological in situ measurements, globally for the period 1976–2024 and for some regions even further back in time. As a result, our new dataset provides long-term trends and interannual variability of glacier mass changes, both based on glacier observations, which can be aggregated at any regional to global scale. Geostatistical modeling allowed the interpolation and extrapolation of observational gaps and uncertainties in space and time. Our uncertainty assessment combines estimates of observational errors from glaciological and geodetic methods, considering spatial autocorrelation, density conversion, and glacier area changes. We note that the mean annual mass balance anomalies are extracted from a handful of glaciers in each region, and thus, in each region, individual glaciers share a large fraction of these variabilities. Importantly, our glacier-wide uncertainty estimates are discussed in a leave-one-out/leave-block-out cross-validation exercise.

Given the available input data, the primary limitation of the new dataset is the scarcity of glaciological in situ measurements, especially in the Southern Hemisphere, and the strongly reduced spatial coverage of geodetic observations before the year 2000. Our dataset would benefit greatly from the support of in situ monitoring programs in regions with a low density of observations, the exploitation of historical archives such as airborne stereoscopic imagery, and the maintenance of future spaceborne missions for cryosphere observation. Our approach would also benefit from scientific advances in volume-to-mass conversion for geodetic estimates and the mapping or modeling of glacier area changes over time.