Measuring glacier mass changes from space—a review

The method consists in coregistering and then differencing two DEMs, i.e. two grids representing the elevation of the glacier surface at two epochs in time. Since 2012, the method has also been modified to tackle time series of DEMs, either using a simple linear regression through the elevation time series (Nuimura et al 2012, Willis et al 2012) or using a more advanced statistical framework to resolve non linearities in the glacier elevation change signal (Hugonnet et al 2021).

A key preprocessing step is the relative adjustment of the different DEMs on the stable terrain surrounding the glaciers. The compared DEMs must share the same map projection and, importantly, be referenced to the same elevation datum (geoid or ellipsoid, see https://doi.org/10.5194/tc-2021-197-CC1). A first and mandatory step is to correct the horizontal and vertical shifts between the DEMs using, for example, the coregistration method from Nuth and Kääb (2011). In the case of spy satellite images (presented below), an affine transformation (rotation, scaling) or more complex non-linear shifts may have to be corrected beforehand (Maurer and Rupper 2015, Dehecq et al 2020). After applying this translation vector, some sensor-specific artefacts may remain and need to be modelled to approach an unbiased elevation change estimate (Nuth and Kääb 2011). For example, undulations in the DEMs in the along-track direction, related to the jitter of the platform, are common in optical imagery (Girod et al 2017). In other cases, tilts or along/cross track biases between the DEMs have to be corrected (Shean et al 2016). If DEMs of different original resolutions are compared, some biases related to the terrain morphology may also persist. Terrain curvature proved to be a good proxy for correcting the latter biases (Gardelle et al 2012). After relative adjustment, DEMs can be subtracted (or processed as a time series) to map elevation changes.

There are two final hurdles before obtaining the glacier-wide total mass balance, two issues which are shared with the altimetry-based estimates discussed in sections 4 and 5 of this review. (a) The first one is the need to fill in data gaps in the map of elevation changes. These gaps are often distributed unevenly and thus, if not filled properly, may lead to biased estimates. For example, native stereo sensors from the late 1990s or early 2000s such as ASTER or Satellite Pour l'Observation de la Terre (SPOT) 5 acquired imagery using eight-bit sensors (for which images are made of 2⁸, i.e. 256 digital numbers) leading to poor radiometric dynamics in the bright, textureless accumulation area of glaciers, and hence a concentration of data gaps. As these upper regions often experience attenuated elevation changes than the lower reaches (Schwitter and Raymond 1993, Belart et al 2019), filling these gaps with a glacier-wide average value would bias the mass balance estimate. Several approaches have been developed to fill data gaps (Seehaus et al 2020), one of the best performing and robust being the local hypsometric approach that exploits the fact that elevation changes are often correlated with altitude (McNabb et al 2019). (b) The second hurdle is the conversion of the volume change to mass change. This requires an estimate of the density conversion factor (Cuffey and Paterson 2010). Here again, various strategies have been used. Assuming an unchanged vertical density profile (the so-called Sorge's law), many studies (e.g. Bader 1954, Arendt et al 2002) used a single ice density (900 kg m⁻³). Others (e.g. Haag et al 2004) have favoured separating the ablation area (where ice is lost at a density of 900 kg m⁻³) and accumulation area, where they used an average density for the firn (often 500–600 kg m³), assuming, implicitly, a change in the vertical density profile. Huss (2013) used a firn densification model calibrated with in situ measurements of firn density to show that a single density of 850 ± 60 kg m⁻³ was a good compromise for glacier-wide estimates over a measurement period longer than five years (and ideally longer). Improvements in this volume to mass conversion are still needed as it remains the largest source of errors in a recent estimate of global glacier mass change (Hugonnet et al 2021).

Due to these uncertainties in the volume to mass conversion and the precision of the elevation change measurements, DEM differencing is typically applied over periods of 5–10 years, although dense time series of DEMs can also allow exploring short term, seasonal elevation changes (Seehaus et al 2015, Belart et al 2017, Beraud et al 2022). The temporal resolution and timeliness of the geodetic mass balance record can be improved with ancillary data such as in situ mass balance measurements (Zemp et al 2019, 2020) or snow line satellite observations (Davaze et al 2020, Barandun et al 2021).

The principle of DEM generation from a stereoscopic pair is simple and has been used for decades by national mapping agencies to chart the Earth's land topography using overlapping aerial photographs. It consists of estimating the distortions (also called the parallax) between pair, triplet or more images acquired from different viewpoints. If the acquisition geometry (position, pointing angles) of the images are known, the parallax can be converted to terrain elevation. For an historical background, the reader is referred to a review (Toutin 2001) describing the early days of DEM generation from optical satellite images.

The main drawback of optical (visible or near infra-red) images, compared to radar data, is that they are weather dependent. No DEM can be derived during the (polar) night or when clouds obscure the Earth's surface. Another notorious limitation is that DEMs derived from these images often concentrate data gaps and exhibit more noise in the snow-covered accumulation areas. We will see below that this is however not true anymore for modern stereo-imagery.

In the following, we distinguished two generations of stereo sensors depending on their resolution. We focus on the sensors that have been most used for estimating glacier elevation changes and do not aim at an exhaustive overview.

3.2.1. Stereo sensors with decametric (5–15 m) resolution.

3.2.2. Stereo sensor with metric to submetric resolution.

In the second part of the 2000s and during the 2010s, a novel generation of high resolution satellites was launched, including Worldview (WV, DigitalGlobe, now Maxar) and Pléiades (CNES and Airbus Defence & Space, DS). These satellites have a sub-metre resolution and are characterised by a high agility so that stereo pairs (or triplet) can be formed from along-track acquisitions a few tens of seconds apart. Two characteristics have dramatically increased the quality of the DEMs that can be derived from these modern sensors over snow and ice. First, the radiometric depth has evolved from 8 bits for ASTER, SPOT1-5 (and others) to 11 or 12 bits (2048 or 4096 digital numbers) for WV1-4 and Pléiades (and SPOT6/7). Second, the sub-metre resolution means that fine scale features of the terrain can be resolved. Both characteristics lead to images with a lot of contrast even in the snow-covered accumulation areas and almost no saturation. This implies that (a) the entire archives are useful for glaciology (without any special gain setting) and (b) the correlation between the stereo-pair works efficiently so that almost complete DEMs can be constructed even over ice sheets (Shean et al 2016, Howat et al 2019) or the glacier accumulation areas (Berthier et al 2014). The higher spatial resolution and the good knowledge of the orbits also results in a higher vertical precision of the DEMs. The improvement in resolution and precision between ASTER and Pléiades is illustrated for the Meighen Ice Cap (figure 4). As a rule of thumb, the vertical DEM precision is on the order of the image resolution.

Zoom In Zoom Out Reset image size

The sub-metre resolution is obtained at the cost of relatively narrow swaths (∼10–20 km) compared to, for example, ASTER (60 km) or SPOT-HRS (120 km). This means that several images, often acquired over the course of several days or weeks, need to be mosaicked to cover a large ice cap or an entire mountain range. Another drawback is that these sensors are commercial, which complicates data accessibility for the scientific community. They mostly acquire images on-demand and contrary to ASTER, they do not continuously build a vast archive of stereo images. Exceptions to this rule are the numerous, time-stamped WV DEMs freely available in the Arctic (ArcticDEM, Porter et al 2018), Antarctic (REMA DEMs, Howat et al 2019) and in High Mountain Asia (Shean et al 2020). Dedicated acquisition campaigns are also performed by the French Space Agency (CNES) through the Pléiades Glacier Observatory ¹⁷ .

3.2.3. Future missions with stereo capabilities.

Spatially-distributed surface elevation information can also be gained from single-pass across-track interferometric (InSAR) data. In this case, a radar system with spatial separation (baseline) between transmitter and receiver is required for the data acquisition, with one antenna transmitting and both receiving the complex SAR signal. This so-called bistatic instrument configuration enables the generation of high quality interferograms that are not affected by temporal decorrelation and ice motion and with reduced atmospheric effects (Bamler and Hartl 1998, Rosen et al 2000).

The cross-track technique is based on the evaluation of phase differences measured by the two SAR antennae. While in a single SAR image the location of each surface point is reduced to the range distance, this technique is a means to measure also the radar look angle allowing to recover the point's location in space. The interferometric phase difference of every pixel is used to determine the difference between the two range distances to the SAR antennae which in turn can be considered a measure of the radar look angle (Bamler and Hartl 1998). A sensitivity measure of the interferometrically derived elevation is the height resulting from a phase change of one fringe (2π) which represents one ambiguity cycle. This height of ambiguity is inversely proportional to the perpendicular component of the baseline, the effective baseline. A larger effective baseline leads to a dense fringe pattern, hence higher precision but also more difficult unwrapping. Conversely, a smaller baseline leads to a noisier DEM but a more robust unwrapping. Additionally, the radar wavelength also has a directly proportional influence, a low frequency SAR (e.g. L-band, 1.25 GHz) generates larger heights of ambiguities and thus less sensitive DEMs than a C- or X-band SAR system (5.3 GHz and 9.6 GHz, respectively).

3.3.1. Elevation data from spaceborne bistatic InSAR systems.

To date, two spaceborne InSAR missions have been launched with the primary objective of measuring Earth's topography.

The first single-pass interferometer, SRTM, was flown from 11 to 22 February 2000 and was a unique experiment with two antennae on the same platform, the second antenna located on the end of a 60 m long mast (Farr et al 2007). Thus, the SAR data were acquired with a known, constant baseline during the mission. The coverage of the resulting DEM is, however, principally limited to a latitude range from 56°S to 60°N due to the inclined orbit of the Space Shuttle and its mapping geometry. The good vertical accuracy (better than 9 m) and the short temporal baseline provided a snapshot of the glacierised areas within these latitudes. The SRTM C-band DEM (processed by the National Aeronautics and Space Administration, NASA), obtained by combining overlapping ascending and descending strip data, has been widely used due to its completeness and it is available for free download in various versions, the latest being released in February 2020 by the Land Processes Distributed Active Archive Center (LPDAAC). The SRTM X-band DEM, processed by the German Aerospace Center (DLR), is also available but has wide gaps due to the narrower swath and is therefore less used.

The second mission, TanDEM-X (TerraSAR-X add-on for digital elevation measurements), is currently in orbit and consists of the coordinated operation of two almost identical satellites—TerraSAR-X (TSX) and TanDEM-X (TDX), the latter named similarly as the mission—flying in close helix formation (Krieger et al 2007). The mission has been operational since December 2010. The primary objective of the TanDEM-X mission is the generation of a worldwide, consistent, timely, and high precision DEM which was released in 2018 (Wessel et al 2018). Operational data acquisition for DEM generation is performed using the bistatic InSAR stripmap mode (30 km swath width): either the TSX or TDX satellite is used as a transmitter to illuminate a common radar footprint on the Earth's surface. The scattered signal is then recorded by both satellites simultaneously. Both satellites act as a large single-pass radar interferometer with the opportunity for flexible baseline selection. Prerequisites for bistatic operation of this fully active system are the pulse repetition frequency synchronisation and the relative phase referencing between the two satellites.

The global DEM product of the TanDEM-X mission (DLR-EOC 2018), whose performances are analysed in Rizzoli et al (2017) and Wessel et al (2018), is the result of the combination of four years (from 12 December 2010 to 16 January 2015) of bistatic data acquisitions with different baselines and geometries. In order to compensate residual offsets and tilts after the interferometric processing, this product is calibrated to selected ICESat points except in Greenland and Antarctica. Hence the global TanDEM-X DEM is not suitable for the derivation of surface elevation changes, but rather as a reference DEM for various purposes. Instead, individual bistatic TanDEM-X acquisitions can be processed to timestamped DEMs and used to generate time series of surface elevations at yearly or multiyear intervals over regions with high mass losses (Abdel Jaber et al 2019, Braun et al 2019). Compared to early volume change assessments for which DEMs from different sensors were combined, the TanDEM-X single-pass interferometer offers significant improvement in terms of spatial resolution (typically 5–10 m) and vertical accuracy. Typically, elevation changes are measured within ±1 m if ice free areas are available for proper adjustment of the DEMs.

The map of surface elevation change rate over the Western Vatnajökull Ice Cap (Iceland) obtained from two TanDEM-X 30 km-wide strips is shown in figure 5. The 2012 DEM is composed of ascending acquisitions in August (western strip) and June (eastern strip). The second coverage was achieved with ascending acquisitions in winter 2016, with the western strip from December and the eastern from November. The DEMs were vertically coregistered over ice free areas. The decrease in elevation from 2012 to 2016 can be clearly observed at the outlet glacier termini, while higher on the ice cap the surface elevation is almost constant or slightly increases. The obvious large spot of surface lowering in the Northwest of the ice cap is caused by the Bárðarbunga caldera collapse in 2014/2015 (Gudmundsson et al 2016). The slight elevation change discrepancy at the boundary between the different DEMs is visible as a linear feature only on the ice cap and may be due to seasonal elevation changes or variable radar penetration.

Zoom In Zoom Out Reset image size

For large basins and on the flat areas at the interior of the ice sheets, TanDEM-X elevation change maps can be complemented with interpolated CryoSat-2 elevation changes at crossover locations (Krieger et al 2020).

3.3.2. Limitations and challenges of InSAR DEMs for elevation change mapping.

3.3.3. Future missions with single pass InSAR capabilities.

Radar altimetry for Earth observation was initially developed for observing ocean topography and is one of the main instrumental records of ocean topography with applications in marine geoid, ocean currents, and sea level change. Its use is now ubiquitous across the cryosphere and inland water. Radar altimeters are nadir-looking instruments recording time between pulse emission and reception, which is then converted to distance. This distance is usually that of the point-of-closest approach (POCA), obtained via retracking of the recorded waveform (Martin et al 1983, Wingham et al 1986, Davis and Moore 1993). Over sloping terrain, POCA is situated off-nadir and requires special procedures to properly relocate the echo (Brenner et al 1983). CryoSat-2 operates a beam-forming radar altimeter which improves the along-track footprint resolution, and its interferometric mode allows accurate echo-location in the across-track direction (figure 6).

Zoom In Zoom Out Reset image size

This revolutionary design led to a paradigm shift in how radar signal is processed for elevation retrieval. In the presence of sloping terrain, CryoSat-2 is akin to the geometry of a side-looking SAR, with the ability to retrieve elevation information beyond the POCA (Wingham et al 2006). The potential of the so-called swath processing has been demonstrated with data acquired by the airborne ASIRAS radar system (Hawley et al 2009) and then by CryoSat-2 (Gray et al 2013, Gourmelen et al 2018), generating between 1 and 2 orders more measurements than POCA from conventional waveform retracking. For a detailed description of swath processing, we refer to Gray et al (2013) and Gourmelen et al (2018). Apart from an increase in the number of elevation measurements retrieved from CryoSat-2, swath processing allows a better sampling of low-lying terrain as the POCA tends to be located in topographic highs (Foresta et al 2016). This leads to a more complete hypsometric coverage of small ice bodies and thus a more representative sampling of ice elevation changes that tends to be strongly correlated with elevation. In addition, since it does not rely on waveform retracking, swath processing allows retrieval of elevation even when the waveform leading edge is difficult to identify such as for high surface slopes and high terrain roughness. The improved spatial coverage coupled with a monthly repeat cycle gives CryoSat-2 unique spatio-temporal capabilities to monitor changes in glacier volume, and ultimately mass.

To assess linear rates of change from swath elevations, the method follows a generic approach developed over ice sheets. The swath elevation data is binned into a grid, with cell resolution typically set around 500 m as a trade-off between spatial resolution and robustness, and then a plane-fit approach is applied in each grid cell (Foresta et al 2016, 2018). The plane-fit algorithm is a simple regression, modelling elevation (z) with the parameters easting (x), northing (y) and time (t). The coefficient for time (t) resolves a constant rate of surface elevation change and is retrieved from the plane-fit model of each individual grid cell, resulting in gridded maps of dh/dt. Several improvements to the method have been made by applying outlier removal, quality filtering and quality weighting methods (Foresta et al 2016, 2018, Morris et al 2020, Jakob et al 2021, Morris et al 2022). The maps of elevation changes typically have coverage of between 60% and 87% (Foresta et al 2016, 2018, Morris et al 2020, Tepes et al 2021). Different gap filling methods have been applied over the years, including spatial interpolation and hypsometric averaging at different spatial scales (Nilsson et al 2015b). Similar to the DEM differencing and laser altimetry methods, the elevation changes are then converted into mass changes, using the glacierised area from glacier masks (section 2) and a volume-to-mass conversion factor (section 3.1).

In regions with highly complex terrain, such as High Mountain Asia and the Gulf of Alaska, the data rate is usually lower due to limitations with closed-loop onboard-tracking (Dehecq et al 2013) and decreasing data quality with higher surface slopes. Larger grid cells are therefore necessary, and since for larger grid cells the plane-fit approach is a too simplistic representation of topography, an auxiliary baseline DEM is required to remove the terrain-related elevation change within each cell (Kääb et al 2012, Jakob et al 2021).

Besides providing linear change over a time-period, CryoSat-2's monthly repeat cycle has enabled multiple studies to resolve time-dependent changes, giving insights into annual and seasonal changes (Gray et al 2015, Foresta et al 2016, 2018, Morris et al 2020, The IMBIE Team 2020, Jakob et al 2021). For monthly change in elevation, additional spatial averaging is needed and resolution is typically set from 2 km for relatively regular surfaces (The IMBIE Team 2020) up to larger glacier units (Jakob et al 2021). Several methods for time series generation exist, from simple monthly averaging (Foresta et al 2018) to more robust statistical weighted mean (Davis and Segura 2001, Gray et al 2015, Malczyk et al 2020).

Limitations related to other elevation change methods (DEM differencing, laser altimetry) are also applicable to radar altimetry-based studies. These include the uncertain volume-to-mass conversion, the need for gap-filling or data interpolation, and reliance on correct glacier masks for generating volume change and mass balance.

Further, Ku-band radar signal scatters within the snow and firn volume, with the potential to impact the accuracy of time-dependent elevation, with spectacular examples over the interior of the Greenland Ice Sheet after extreme melt events (Nilsson et al 2015a). This spectacular case study however stems from the combination of a thick and dry firn, and of the absence of repeated and frequent melt events. Studies conducted over a large range of settings have demonstrated the consistency between radar-altimetry and airborne laser or field-based multi-annual elevation change (Gray et al 2015, Foresta et al 2016, Gourmelen et al 2018, Tepes et al 2021), as well as ways to identify change in penetration and to correct for its impact (Nilsson et al 2015a, Slater et al 2019, Gray 2021). Scattering properties can also induce elevation biases at a seasonal timescale (Gray et al 2019, Morris et al 2022), and efforts are underway to understand and address their impact and exploit the potential of radar altimetry for mapping seasonal change.

Over the course of the last decade, CryoSat-2 has demonstrated that beam-forming, interferometric, radar altimetry can monitor glacier change globally, providing a unique blend of spatial and temporal resolution. It offers the possibility to analyse glacier change in response to complex atmospheric and oceanic forcing. The combination of radar altimetry with other available techniques and constantly improving glacier and climate models should lead to more accurate mass balance estimates and better process understanding.

Importantly, the European Commission Copernicus expansion mission concept CRISTAL (Kern et al 2020), a radar altimeter inspired by CryoSat-2 and AltiKa, and due to launch in 2027, will be the first operational mission with a primary objective to monitor glaciers globally for decades to come (figure 2).

The main goal of ICESat was to measure the elevation change of the ice sheets to determine their contribution to sea level change (Schutz et al 2005). Despite design limitations that shortened the operational lifetime of the lasers (Abshire et al 2005), the mission provided global surface elevation measurements of unprecedented accuracy during its seven year mission life (2003–2009). ICESat recorded data during two to three month-long campaigns per year. The data consists of full-waveform 1064 nm laser returns of an ∼65 m elliptical laser footprint (Schutz et al 2005). The mission operated in a 91 day repeat orbit and a 94° inclination that covered latitudes up to ±86°, with precision pointing employed in polar regions to steer the laser footprints to within ±300 m of reference ground tracks. In lower latitudes, the reference ground track pattern was not repeated with as good accuracy. The along-track distance between measurements is 172 m, but the distance between reference ground tracks ranges from ∼6 km close to the poles up to 89 km at the equator. For land ice studies, there exist preprocessed level-2 data products: GLAH12 (optimised for ice sheets), GLAH14 (commonly used for mountain glaciers), and GLAH06 (level 1B, used for smooth polar glaciers). In addition to an elevation value for each footprint, these products include numerous attributes and corrections. The latest release of the data, 2014, is version 34. ICESat's revolutionary data provided multiple glaciological insights, in particular for volume changes of ice sheets (Pritchard et al 2009, Sørensen et al 2011) and glaciers (Gardner et al 2011, Kääb et al 2012). The success of ICESat led to the follow-up mission ICESat-2, launched eight years after the end of the ICESat mission. To bridge the gap between the two missions, NASA'S Operation IceBridge provided annual airborne altimetry data for selected areas in the Arctic and Antarctic (MacGregor et al 2021).

ICESat-2 was launched in September 2018, and had a minimum mission lifetime of three years (Markus et al 2017) met in December 2021. Unlike its predecessor, ICESat-2 uses a digital photon counting detector and 10 kHz micro-pulse laser, providing a greatly increased number of elevation measurements with 0.7 m along-track spacing between pulses (Neumann et al 2019). ICESat-2 operates in a 91 day repeat orbit, like ICESat, but with a latitude coverage up to ±88°. The laser micro-pulse is split into three-beam pairs (six beams total), with each pair consisting of a weak (¼ strength) and strong beam (full strength). Within each pair, beams are separated on the ground by 90 m and the beam pairs by 3.3 km (Markus et al 2017). The surface footprint of each beam is approximately 11 m in diameter (Luthcke et al 2021, Magruder et al 2021). Data is collected nearly continuously globally but the satellite only operates in exact repeat mode in polar regions, with a pointing control of 4.4 m (Luthcke et al 2021). At lower latitudes the mission off-points the instrument to prioritise vegetation coverage over repeat observations. Level 2A Geolocated Photons (ATL03) is the ICESat-2 product of primary interest for studies of detailed glacier features such as crevasses (Herzfeld et al 2021) or supraglacial lakes (Fair et al 2020). Most suitable for glacier change studies are Level 3A Land Ice Height (ATL06) and L3B Annual Land Ice Height (ATL11, corrected for surface slope), that are binned in the along-track in 40 m and 60 m segments, respectively. Level 3B gridded ice height and height change data (ATL14, ATL15) were made public in late 2021.

The 1064 nm laser on board ICESat led to little multiple scattering within snow and ice and therefore no elevation bias (Smith et al 2018). ICESat-2 uses a 532 nm laser that, under certain conditions, can experience a large degree of multiple scattering within the snow and ice. The magnitude of this elevation bias depends on snow conditions and the method of identifying the surface height from the photon cloud (Smith et al 2018). It can be as large as 24 cm in extreme cases (e.g. pure snow/ice with extremely large grains). However, for most applications, the bias is negligible (2–5 cm) and largely cancels for repeat measurements under similar snow conditions.

ICESat was only capable of pointing to the same location on the ground (i.e. reference ground track) to within 10s to 100s of metres and measurements were made every ca. 170 m along-track. This makes it challenging to directly compare near-repeat elevation measurements, as elevation changes will be convolved with differences in elevation due to changes in measurement location (figure 7). This is often referred to as the 'slope error' and must be corrected for. ICESat-2 has far superior pointing control (3–25 m) and closely spaced beam pairs (90 m) that make correcting for 'slope error' robust and straightforward. In steeper or rougher terrain or in the mid-latitudes without quasi-exact repeat tracks, different ground tracks may not be directly comparable altogether without an additional reference elevation dataset and the use of spatial statistics (figure 7). In all cases, the sampling is rather sparse at the scale of individual glaciers or ice caps, and most of the times uneven between the ablation and accumulation areas, such that a careful extrapolation to the entire ice body is required.

Zoom In Zoom Out Reset image size

Whatever the method, accurate glacier masks are required to correctly classify individual surface elevation measurements (section 2). Inclusion of measurements on stable terrain outside glaciers (e.g. when an outdated glacier inventory is used) will lead to almost unbiased volume change estimates (the off glacier elevation changes will sum up to almost 0) but an underestimate of glacier thickening/thinning area-averaged rates because the correct total volume change is divided by an area which is too large.

5.2.1. Methods for calculating elevation change from near-repeat observations.

5.2.2. Methods for non-repeat observations and rough surfaces.

New methods take time to develop for new products like ICESat-2, and as of this writing, there are relatively few results in the literature. It is expected that ICESat-2 data will contribute to new insights and a much better understanding of the entire cryosphere, including dynamical changes of the ice sheets (Smith et al 2020) and improved monitoring of mountain glaciers (Wang et al 2021, Fan et al 2022), but also advances in related disciplines such as, snow depths on- and off-glacier (Treichler and Kääb 2017), sea ice and ice shelf thickness (Kacimi and Kwok 2020) or permafrost changes (Michaelides et al 2021). In addition, the accurate elevation measurements of both ICESat and ICESat-2 provide a globally consistent reference elevation dataset that has the potential to greatly improve the quality and accuracy of existing and upcoming DEMs and other elevation data across the globe, which will directly benefit glacier volume change studies.

There are currently no funded Earth observing laser altimetry missions under development. However, the need for the continuation of cryosphere focused laser altimetry is highlighted in the 2018 Decadal Survey for Earth Science and Applications from Space, a guiding document for NASA's future Earth science investment, and has been recommended as one of seven observables to be completed as an Earth System Explorer class mission.

The Earth's gravitational field is described by the geopotential V, i.e. the potential energy per unit mass at a certain point. At a point above the Earth's surface, with spherical coordinates radius r, co-latitude θ and longitude λ, it can be expressed as a sum of Legendre functions:

$$\begin{align*}\begin{array}{l} V(r,\theta ,\lambda ){\mkern 1mu} = {\mkern 1mu} \frac{{GM}}{{{a_e}}}\left\{ {\sum\limits_{l = 0}^\infty {\sum\limits_{m = 0}^l {{{\left( {\frac{{{a_e}}}{r}} \right)}^{l + 1}}} } {P_{lm}}(\cos \theta )} \right.\\ \left. {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \times ({C_{lm}}\cos {\mkern 1mu} m\lambda {\mkern 1mu} + {\mkern 1mu} {S_{lm}}{\mkern 1mu} \sin m\lambda )} \right\} \end{array}\end{align*}$$

where G is the gravitational constant, M the mass of the Earth and a_e denotes its mean equatorial radius. P_lm are the Legendre polynomials of degree l and order m, and C_lm and S_lm are spherical harmonic coefficients. For further details, we refer the reader to Wahr et al (1998). As can be seen from the formula above, coefficients with a higher degree and order are related to increasingly smaller wavelengths. At the same time, the scaling by (a_e /r)^{l +1} implies that the effect of these coefficients have a smaller contribution to the geopotential than the lower-order coefficients. In practice, this means that the orbit of the GRACE satellites are dominated by the large-scale mass distribution of the Earth and that small-scale features, such as glacier mass changes, are more challenging to measure than large-scale features, such as ice-sheet wide mass changes.

By collecting a sufficient amount of range-rate observations, a model of the Earth's gravity field can be derived. The GRACE/GRACE-FO mission provides updated gravity fields typically every month, although weekly or even daily updates are available as well (at the cost of a lower resolution).

At seasonal to decadal time scales, variations in the gravity field are mainly induced by redistribution of water on the Earth's surface, although processes within its interior such as glacial isostatic adjustment (GIA) may contribute as well, as we will see later on. Using an appropriate scaling (see Wahr et al 1998 for details), variations in the geopotential can be expressed as variations in surface water loading. This is usually expressed in units of metre water equivalent height. It should be kept in mind that GRACE can only measure anomalies in surface water loading, not the absolute weight of a glacier or ice sheet.

GRACE data are distributed by the GRACE science data centres, and other parties, as spherical harmonics (the C_lm and S_lm coefficients in the equation above). Since higher-order coefficients are increasingly affected by observational noise, only coefficients up to degrees and orders smaller than 60 or 90 are typically used in the derivation of the surface water loading anomalies, which corresponds to a maximum spatial resolution of approximately 200–300 km. A user-friendly alternative for the spherical harmonics are mascons, i.e. mass concentration blocks. These provide local mass variations on a regularly spaced grid. Compared to the spherical harmonics, they are less affected by noise and therefore require little or no post-processing. Mascons products are available with a grid resolution down to a few tens of kilometres, but it should be kept in mind that the inherent resolution is still in the order of a few hundred kilometres (in other words, the signal in adjacent grid cells will be correlated). Furthermore, once the grid cells are predefined by the processing centres, they do not always cover the area of interest exactly and may contain signals from neighbouring ice bodies (such as the Greenland Ice Sheet when studying the glaciers of the Canadian Archipelago or in Iceland (Sørensen et al 2017)) or other non-glacier signals (e.g. groundwater and soil moisture signals in Patagonia or High Mountain Asia).

Figure 8 illustrates the linear trend in surface water load as observed by GRACE/GRACE-FO from 2002 to 2021. Together with the polar ice sheets, mountain glacier regions stand out with rates in the order of 20–50 cm yr⁻¹ water equivalent mass loss in Alaska, the Arctic Archipelagos and Iceland.

Zoom In Zoom Out Reset image size

To obtain the total mass balance for a region, the local surface mass anomalies need to be summed up over the region of interest. When working with mascons, this is relatively straightforward, but for the spherical harmonics some intermediate steps are required. If the Stokes coefficients were available up to an infinite degree l, we would be able to create a perfect region mask and could simply integrate the mass anomalies over the region. However, as mentioned earlier, the spherical harmonics are only provided up to a limited degree, typically 60 or 90. This, together with the postprocessing often applied to reduce noise in the higher degree coefficients, leads to attenuation and spreading of the signal (Swenson and Wahr 2002). Here, spreading (or leakage) means that due to the limited resolution of the gravimetry data, mass variations from adjacent areas will spill into the region of interest. Likewise, part of the signal in the region of interest will be spread out beyond the boundaries of that region, e.g. in the ocean although recent mascon approaches minimise this effect compared to former spherical harmonic solutions (Wiese et al 2016). Only considering the glacierised areas would therefore result in a biased estimate. Several methods have been proposed to circumvent these limitations. The most common approach is to place one or several unit mass loads in the glacierised areas of interest, represent these as Stokes coefficients with the same characteristics as the GRACE mass anomaly data (e.g. limited to the same maximum degree and applying a similar post processing) and finally applying a scaling to the mass loads so that the difference with the actual GRACE observations is minimised (Wouters et al 2008, 2019, Jacob et al 2012, Chen et al 2013, Gardner et al 2013, Ciracì et al 2020).

Just like any satellite observation, the GRACE data comes with noise, which reflects itself as random fluctuations in the mass anomaly time series in the order of a few gigatons. This means that in regions where the glacier mass change signal is small, the signal-to-noise ratio will be small, and the GRACE results should be interpreted with care. For the larger glacier regions, however, GRACE can capture mass variations with a high degree of confidence.

Even after extracting the mass changes in the glacier regions, there are confounding processes that need to be removed to isolate the glacier signal. Because GRACE has no vertical resolution, it senses the integrated sum of mass redistribution within its footprint, i.e. the sum of the water load changes at the surface and those in the subsurface. Many glaciers are located in regions that were widely covered by thick layers of ice during the Last Glacial Maximum (about 20 kyr ago), which depressed the Earth's surface and led to mass redistribution in the Earth's mantle. This ice has disappeared since, but the Earth's highly-viscous mantle material is still readjusting to the changes in surface load. The associated mass redistribution in the mantle will also be observed by the GRACE satellites. Several models exist to correct for this process, called GIA, which can be used to correct the GRACE data. Whereas models show a wide spread in Antarctica (Martín-Español et al 2016), they are fairly well constrained in mountain glacier regions and the Arctic and GIA mass rates agree in most regions within a few Gt yr⁻¹ (Blazquez et al 2018). Yet, when summed over all glacier regions globally, differences become significant. For example, the GIA corrections used in two of the most recent GRACE-based glacier mass balances studies, Ciracì et al (2020) and Wouters et al (2019) amounted to 58 ± 11 Gt yr⁻¹ and 34 ± 21 Gt yr⁻¹, respectively ca. 18% and 13% of the total glacier mass change signal (figure 9). Furthermore, some regions with a low mantle viscosity are still adjusting to the load changes related to the Little Ice Age (Sørensen et al 2017). A similar model-based correction is applied for this, although these models are underconstrained in some areas and not all studies include this correction. Another solid earth effect, i.e. mass redistribution associated with earthquakes, is harder to correct for, but is mostly of small magnitude and manifests itself as a step change in the mass time series followed by a slow recovery. In the few cases where this affects glacier mass balance estimates, this can often be identified visually.

Zoom In Zoom Out Reset image size

The coarse resolution of the GRACE observations also implies that, even after correcting for GIA, the signal is a mix of mass changes related to glaciers and hydrological processes, such as changes in the surface, canopy and soil moisture water storage, and in the groundwater table. This is practically an issue in regions where glaciers are scattered and small compared to GRACE's footprint, for example in Central Europe and the Caucasus. Furthermore, strong hydrological variations outside the region of interest may affect the glacier mass balance estimates through the leakage effect discussed earlier (An et al 2021). Again, this is accounted for using models, where it should be noted that most models only represent part of the hydrological cycle, e.g. excluding groundwater storage changes and anthropogenic influences, and tend to underestimate long-term changes in hydrological water storage (Scanlon et al 2018). Furthermore, the magnitude of corrections depends strongly on the choice of the model. The global hydrology correction summed up to +6 ± 7 Gt yr⁻¹ in Ciracì et al (2020) vs. −11 ± 6 Gt yr⁻¹ in Wouters et al (2019) using a different model. Because of the disagreement between models, some studies prefer to omit the correction and include it as an uncertainty (Reager et al 2016). Whereas the impact of hydrology is small in the polar regions, it is the major source of uncertainty for most mid-latitude glacier systems. A further complication is that if part of the glacial meltwater runs off in a proglacial lake, a nearby endorheic lake or an aquifer, the net local mass change detected by the GRACE satellites will be smaller than the actual glacier mass loss. Since such processes are currently not included in hydrological models, this will lead to an underestimation of glacier mass balance. Finally, leakage of the glacier mass change signal in the ocean can also be an issue for glacier regions fringing the ocean. This effect is reduced in the mascon-based solutions.

The GRACE-FO mission is currently our only tool to directly monitor glacier mass changes at a global scale, i.e. without a volume to mass conversion. Although it comes with its own limitations, it has proven to be essential to study the response of the larger glacier systems to climate variability. Likewise, the GRACE/GRACE-FO data has become an indispensable source of information in many other fields, and the importance of continued observations has been recognized by the wider scientific community and decision makers (Tapley et al 2019). Several initiatives, including a joint effort of NASA and ESA, have recently emerged to study concepts for future gravity missions (respectively named Mass Change and MAGIC), consisting of one or more pairs of satellites. This would lead to an improved spatial and temporal resolution of the observations, and most importantly, continuity of this unique time series of global glacier mass balance that started in 2002 with the launch of the GRACE satellites.

The different methods presented in this review are applied to the Vatnajökull Ice Cap, a ∼8000 km² ice mass located in Southeast Iceland (figure 10). The goal of this comparison is not to single out the best-performing technique or to provide a reconcile mass change estimate, but rather to illustrate the capabilities and characteristics of all these techniques in a concrete case study.

Zoom In Zoom Out Reset image size

Vatnajökull is large enough so that each technique can (almost) resolve it spatially. Another advantage is that comprehensive glaciological field surveys are available (Björnsson et al 2013, Aðalgeirsdóttir et al 2020). These in situ data measure the SMB and have been performed seasonally since September 1991. In Aðalgeirsdóttir et al (2020), the SMB was corrected for non-SMB components and calving fluxes (Jóhannesson et al 2020) to calculate the total mass balance, the quantity directly comparable to the remote sensing estimates. Applying or not this challenging correction leads to strikingly different results, this is the reason why both time series (i.e. with/without non SMB components) are retained in figure 11.

Zoom In Zoom Out Reset image size

The dh/dt map derived from InSAR TanDEM-X DEMs for West Vatnajökull was shown in an earlier section of this review (figure 5). In principle, a complete coverage of the ice cap could be derived from InSAR DEMs at high spatial resolution. Despite a shorter survey period (3.5 years), the spatial pattern of elevation changes for West Vatnajökull is highly similar to the pattern obtained from stereo-imagery and CryoSat-2 swath altimetry. The spatial resolution of the InSAR dh/dt map (6 m) is an order of magnitude higher than the one from ASTER (100 m) and almost two orders of magnitude higher than for CryoSat-2 (500 m). Fine details of the elevation changes are captured in the InSAR map (e.g. the drainage of the Skaftá cauldrons, Magnússon et al 2021). If a complete map of elevation change was available from TanDEM-X DEMs, the specific challenges to derive an unbiased estimate of the mass balance would be the need (a) to account for temporal changes in X-band radar penetration and (b) to correct the seasonal bias due to the varying acquisition dates of the DEMs.

The dh/dt map based on ASTER DEMs (figure 10(A)) has a 100 m posting to reduce computing time in a recent global scale study (Hugonnet et al 2021) but could potentially be derived at 30 m which is the posting of ASTER DEMs. The map is gap-free thanks to the relatively good abundance of ASTER DEMs over Iceland compared to other regions on Earth. The massive thinning toward glacier fronts and the surface expression of subglacial volcanic activities (such as the Bárðarbunga caldera collapse (Gudmundsson et al 2016)) are well captured. A larger noise level is observed in the accumulation areas (where the thickening is less homogeneous, noisier than for CryoSat-2), a consequence of the limited radiometric dynamic in the eight-bit ASTER stereo-imagery over homogeneous snowfields (figure 4). The time series of cumulative mass change (figure 11) exhibits a mean seasonal cycle of lower amplitude compared to other techniques and only resolves the variability over periods longer than 4–5 years. One should note that this continuous time series, obtained through Gaussian Process regression, is predicted at a monthly time step but is derived from the temporal interpolation of moderate precision and temporally inconsistent ASTER acquisitions (50 in 20 years, on average). Hence, the inter-annual variability of the mass change, e.g. the positive mass balance during glaciological year 2014–2015 or the strongly negative mass balance in 2009–2010 following the Eyjafjallajökull volcanic eruption (Aðalgeirsdóttir et al 2020), is not captured.

CryoSat-2 swath altimetry also provides an almost complete dh/dt map, covering 85% of the ice cap area (figure 10(B)). To calculate an ice cap-wide mass balance, gaps are filled using the regional hypsometric approach (McNabb et al 2019), i.e. using the mean elevation changes of the rest of the ice cap at the same altitude band. CryoSat-2 captures subtle changes in the ice cap interior as well as rapid changes related to the Bárðarbunga caldera collapse or strong thinning at the margins. There is some data loss over the steeper, marginal, parts of the ice cap, including the rapidly thinning margins. CryoSat-2 provides a smooth and continuous time series with monthly estimates thanks to the satellite's repeat period. The monthly repeat allows resolving the seasonal and inter-annual changes, with an amplitude in good agreement with GRACE (before 2016) and with the SMB measurements uncorrected for non surface-mass-balance processes.

Laser altimetry does not allow continuous temporal sampling due to the gap between ICESat (ending in 2009) and ICESat-2 (launched in 2018). Figure 10(C) shows elevation change rates from ICESat (GLAH06 L1B Global Elevation Data v34: Zwally et al 2014) and ICESat-2 (ATL06 L3A Land Ice Height, Version 3: Smith et al 2021) data. The data span the period February 2003 to September 2020. Rates of change are determined by solving least squares fits of offset, rate and seasonality to elevation anomalies within a 480 m radius of solution points that are separated by 480 m which explain the 'stepwise' lines of measurements (figure 10(C)). Elevation anomalies are determined relative to the ArcticDEM 8-m Mosaicked DEM v3 (Porter et al 2018). Despite a rather scarce spatial sampling, laser altimetry is able to capture the main patterns of change, with thinning in the lower reaches and slight thickening in the accumulation area of the two largest outlets of northern Vatnajökull (Dyngjujökull and Brúarjökull). However, ICESat/ICESat-2 has no measurement below 450 m a.s.l. and misses the very strong thinning rate on the tongues of the south and south-east flowing outlets. Given this sparse spatial sampling of a complex pattern of elevation changes, we did not attempt here to estimate the total volume and mass change of Vatnajökull. Hence, laser altimetry data are not plotted in figure 11 and not listed in table 2.

Due to their coarse spatial resolution, gravimetric observations are not able to resolve Vatnajökull mass change alone (figure 10(D)). We corrected the GRACE-based mass change for entire Iceland (updated from Wouters et al 2019), with a scaling factor (0.68) deduced from an analysis of glacier mass change for all major ice masses in Iceland, as 68% of the glacier mass loss in Iceland from 2000/2001 to 2018/2019 originates from Vatnajökull (Aðalgeirsdóttir et al 2020). The processing of GRACE data was detailed in section 6. In Iceland, particular attention needs to be paid to leakage from the Greenland Ice Sheet or to the ocean together with solid Earth processes, related to mass load changes from GIA and little ice age (Sørensen et al 2011). The seasonal cycle is not as smooth as the one from CryoSat-2 with spikes in the time series, due to measurement noise and the leakage from seasonal snow in the surroundings of the ice cap. The GRACE time series also seems to diverge from the SMB-only in situ and CryoSat-2 measurements after 2016, when it gets close to the ASTER-based results and the glaciological measurements corrected for mass losses not occurring at the surface. Reasons for these differences in the time series remain unclear. This may be explained by less reliable measurements towards the end of the GRACE mission and the fact that GRACE/GRACE_FO also captures the mass loss by other Icelandic ice caps. However, a recent study found a good agreement between GRACE and CryoSat-2 mass balance time series for all icelandic ice caps lumped together (Noël et al 2022). Hence, a careful comparison of the CryoSat-2 time series from different research groups is needed.

GRACE excluded, all satellite-based techniques include a conversion of the volume change to mass change, which is done here using an average value of 850 kg m⁻³ (Huss 2013). This value is in principle valid for periods of more than five years (ideally more) and more appropriate for negative mass balance years. This uncertain volume-to-mass conversion remains a major challenge and partly jeopardises the usage of geodetic observations at high temporal resolution. It can easily explain some of the observed differences with GRACE and the in situ measurements, especially during years of balanced or positive mass budgets. Volcanic activity is another source of discrepancies: a striking example is the 65 m collapse (∼2 km³) of the Bárðarbunga caldera in 2014 in the northwest of Vatnajökull (Gudmundsson et al 2016) that can be wrongly attributed to glacier change by elevation change methods.

The average rates of mass change for two different periods are compared in table 2. Over the longest common period (15 September 2002–15 September 2019), the three available estimates show a satisfactory agreement, within ∼10%. However, over the shorter nine year period (15 September 2010–15 September 2019), differences get larger (±25%) with reduced mass loss from CryoSat-2, higher mass loss from ASTER DEMs and intermediate values for GRACE and the 'in situ' measurements corrected for the non-surface components (geothermal melting, volcanic eruptions, the energy dissipation in the flow of water and ice, and calving). Without this correction, in situ measurements are in agreement with CryoSat-2. It is beyond the scope of this review to conclude on the origin of these differences. Understanding them and providing a reconciled estimate of mass change require further work, in particular by examining individual glaciers of the ice cap.

To complement the above case study over a large and flat arctic ice cap (Vatnajökull), we illustrate in figure 12 the ability of the different elevation change methods to sample spatially and temporally glaciers in mountainous areas. The study area covers 100 km by 100 km around Mt Everest (Himalaya), which corresponds roughly to 1/10 of the GRACE resolution. This is why gravimetry data are not shown in this example, they simply cannot resolve glacier mass change at this scale.

Zoom In Zoom Out Reset image size

DEM differencing resolves the changes of individual glaciers and provides an almost 50 year long record albeit with a 5–10 year resolution. There are large data gaps in the elevation change maps due to the lack of contrast, mostly in the upper accumulation areas. This is particularly the case when only two DEMs are compared (figure 12(A)). These gaps are reduced when processing a time series of ASTER DEMs (figure 12(B)). Gaps would almost vanish if DEMs were derived from 11 or 12-bits modern stereo sensors as illustrated by the complete map of elevation differences for Mera Glacier (30 km south of Mt Everest) obtained by comparing Pléiades DEMs from 2012 and 2018 (Wagnon et al 2021).

Altimetry data, from ICESat and CryoSat-2, are not able to resolve individual glaciers. They are only able to provide a region-wide mass balance estimate when aggregated over sufficiently large regions so that glaciers (in particular the different elevations bands) are well sampled. The denser sampling by CryoSat-2 swath altimetry and ICESat-2 allows resolving regions of 100 by 100 km and infer the seasonal elevation changes (Jakob et al 2021, Wang et al 2021). For ICESat, the spatial sampling was sparser and temporally more irregular so that mass loss estimates were mostly restricted to a single value (for the period 2003–2008) over large regions (Kääb et al 2015).

Table 3 summarises the characteristics, strengths and weaknesses of the remote sensing methods to measure glacier mass changes. It includes some general and well-known limitations such as the dependence on cloud-free conditions for optical stereo-images and laser altimetry and, on the other hand, the need to account for time-variable radar penetration in snow and firn for SAR interferometry or radar altimetry.

As illustrated by our case studies in Iceland and Himalaya, the ability of these methods to resolve mass change depends on the glacier-type, their surrounding topography and the size and latitude of the area. Due to their high resolution (here decametric), DEM-based methods excel in mountainous regions hosting, mostly at low to middle latitudes, a myriad of glaciers of varying size, but generally less than 1000 km². Mass balance can be estimated for each individual glacier, albeit with larger uncertainties for smaller ice bodies. In these regions, the generous availability of stable terrain also facilitates the DEMs coregistration and bias correction. Conversely, the rough topography and modest mean glacier size result in a sparse sampling by laser or radar altimetry so that mass balance can be only computed over large regions, typically accounting 10³–10⁴ km² of glaciers. The gravimetric method will also be restricted to region-wide average, and the main challenge is to account for other natural and anthropogenic processes (hydrology, solid earth) influencing its mass change signal.

Over large icefields (Alaska, Patagonia) or Arctic ice caps, DEM differencing is more challenging to apply because of the scarcity of stable terrain, the need to mosaic multiple DEMs of different dates for a full coverage and, specifically for DEMs derived from pre-2010 optical images, the reduced contrasted in the flat accumulation areas. Laser and radar altimetry are well-suited providing a high temporal resolution, down to seasonal time scale. For laser altimetry this is particularly the case for polar areas where ground tracks are repeated exactly. GRACE also works well in principle, although here again the main challenge is to account properly for all non-glacier signals, generally through modelling.

A common limitation of all techniques is that none of them is able to 'see' the loss (or gain) of ice below the waterline (sub-aqueous loss or gain) for water-terminating glaciers experiencing calving front migration. This omission does not matter (or very little) for estimating the sea-level contribution as the solid freshwater of the glacier tongue is displacing salted sea water (Seehaus et al 2015).

All these limitations and sources of uncertainties lead to a relatively large spread between the mass change estimates of the different techniques. This is well illustrated by a compilation for all glaciers and ice caps in the Russian Arctic updated from Tepes et al (2021) (figure 13, table 4). For this large RGI region (covered by 51 500 km² of glaciers), even using the same technique (Gravimetry or DEM differencing), the losses for similar time periods differ by a factor of almost two (table 4). Hence, a first step in an intercomparison exercise will consist in understanding the reasons behind the differences between similar techniques.

Zoom In Zoom Out Reset image size