# Multi-technique estimation of ice mass balance in Greenland: impact of the uncertainties on firn densification and GIA models

Ana Sanchez Lofficial, Laurent Métivier, Luce Fleitout, Kristel Chanard, Marianne Greff-Lefftz, Maylis de la Serve, Louis-Marie Gauer, Emma Gourrion

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

## To cite this version:

Ana Sanchez Lofficial, Laurent Métivier, Luce Fleitout, Kristel Chanard, Marianne Greff-Lefftz, et al. Multi-technique estimation of ice mass balance in Greenland: impact of the uncertainties on firn densification and GIA models. Geophysical Journal International, 2025, 240, pp.1935-1952. 10.1093/gji/ggaf015. insu-05080926

HAL Id: insu-05080926

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

Submitted on 23 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Multi-technique estimation of ice mass balance in Greenland: impact of the uncertainties on firn densification and GIA models

Ana Sanchez Lofficial,<fnref n="1" /> Laurent Métivier,<fnref n="1" />,<fnref n="2" /> Luce Fleitout,<fnref n="3" /> Kristel Chanard,<fnref n="1" />,<fnref n="2" /> Marianne Greff-Lefftz,<fnref n="1" /> Maylis de La Serve,<fnref n="1" />,<fnref n="2" />,<fnref n="4" /> Louis-Marie Gauer<fnref n="1" />,<fnref n="5" />,<fnref n="6" /> and Emma Gourrion<fnref n="1" />

Accepted 2024 December 26. Received 2024 December 24; in original form 2024 April 26

## SUMMARY

We conduct a comprehensive comparison of ice mass balance (IMB) estimates for Greenland derived from satellite observations of ice surface elevation changes (SEC), gravity and global navigation satellite system (GNSS) observations. Our analysis integrates data from the ICESat and CryoSat-2 satellite altimetry missions, augmented by optical stereo-imagery for peripheral glaciers, and GRACE satellite gravimetry mission, spanning the 2003–2008 and 2011–2015 periods. We also consider three firn densification models (FDM) and five glacial isostatic adjustment (GIA) models for correcting the data sets for these effects when necessary. Our results reveal significant differences among FDM corrections applied to SEC observations, with particularly large variations in IMB estimates reaching up to 90 Gt yr⁻<fnref n="1" />. To address this, we develop an innovative method for estimating equivalent firn corrections to the ice elevation observations, based on a least-squares fit of filtered ice SEC observations to GRACE mass-change estimates. This approach is both simple and independent from climate models assumptions and shows minimal sensitivity to GIA model differences. Using this method, we estimate IMBs for Greenland at -217.6 ± 15.7 Gt yr⁻<fnref n="1" /> for 2003–2008 and -253.2 ± 18.8 Gt yr⁻<fnref n="1" /> for 2011–2015. Importantly, these values indicate an acceleration of the thinning rate, not consistently captured by the IMB estimates inferred from the ice SEC observations corrected by FDMs. Finally, we compute elastic ground deformation induced by ice mass change during 2011–2015, using the four proposed mass-variation distributions and compare the predicted vertical velocities with GNSS observations in Greenland, accounting for all GIA models. While all models are consistent with most of the GNSS-derived uplift rates, they cannot fully explain the observed vertical velocities, especially in the South-East Greenland, which confirms the need to refine our understanding of GIA contributions in this region.

**Key words:** Gravity; Geodesy; Glacial rebound.

## 1 INTRODUCTION

Greenland ice mass balance (IMB) estimates are commonly inferred from space-based observations including satellite gravimetry, altimetry, often supplemented with surface mass balance information, and GNSS measurements of the associated ground deformation (e.g. Velicogna 2009; Khan et al. 2016; Zhang et al. 2019; Mouginot et al. 2019; Velicogna et al. 2020; Shepherd et al. 2020). Numerous IMB estimates have been published in recent years, but

they present large uncertainties, driven by the inherent advantages and shortcomings of each spatial observational technique used (e.g. Cazenave & Llovel 2010; Métivier et al. 2010; Shepherd et al. 2012, 2020). Our goal is to integrate various techniques to improve the precision of Greenland ice mass balance estimates over different time periods. A key focus is to analyse the sources of uncertainty inherent in each technique and explore strategies to reduce these uncertainties by leveraging the strengths of combined observational methods.

<footnote n="1">¹Université Paris Cité, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France. E-mail: lalmetiv@ipgp.fr</footnote>

<footnote n="2">²Univ Gustave Eiffel, ENSG, IGN, F-75238 Paris, France</footnote>

<footnote n="3">³Ecole Normale Supérieure and CNRS UMR8538, PSL Research University, F-75005 Paris, France</footnote>

<footnote n="4">⁴Department of Earth and Environmental Sciences, Michigan State University, 48824 East Lansing, MI, USA</footnote>

<footnote n="5">⁵LEGOS, Université de Toulouse, CNES, CNRS, UPS, IRD, F-31400 Toulouse, France</footnote>

<footnote n="6">⁶CNES, Centre Spatial de Toulouse, F-31401 Toulouse Cedex 9, France</footnote>
Satellite altimetry, based on radar (ERS-1/2, Envisat, CryoSat-2) or laser (ICESat 1-2) measurements, measures changes in ice sheet surface height with temporal resolutions of typically a few years. With records extending back to the 1970s, satellite altimetry provides one of the longest records of ice sheet surface elevation changes (SEC), giving precious information on its decadal evolution. Its high spatial resolution data are crucial for modelling ground deformation caused by ice unloading across both short- and long-spatial wavelengths, especially when compared with global navigation satellite system (GNSS) ground deformation observations. However, satellite altimetry does not directly give access to changes in the ice mass variations. The SEC must first be converted to ice mass changes, accounting for processes occurring within the firn layer and the underlying ice layers, which can be challenging (e.g. Huss 2013). The firn layer refers to the upper layer of the total ice/snow column, where snow accumulates and compacts over time, until it becomes ice. These firn layers are characterized by densities that can span a wide range of values, from ~300 to ~917 kg m⁻³, thus introducing potentially significant uncertainties in mass change estimates from measurements of variations in ice elevation. In order to estimate ice mass changes from altimetry observations, it is possible to use a firn densification model (FDM) (e.g. Stevens et al. 2020), that is, a model predicting snow compaction processes within the firn layer, given surface mass balance information inferred from climate models or re-analyses (e.g. Fetweis et al. 2020). In contrast, time-varying satellite gravimetry missions (GRACE and GRACE-FO missions) provide a more 'direct' measurement of ice mass variations, without FDM assumption, offering over 20 yr of observations. Compared to satellite altimetry, space-based gravimetry provides observations with a monthly temporal resolution but at coarser spatial resolution (a few hundred kilometres).

GNSS records of vertical ground displacement also help constrain or validate IMB estimates (e.g. Jiang et al. 2010; Khan et al. 2010, 2013, 2016; Bevis et al. 2012; Wahr & Zhong 2013). Elastic or viscoelastic deformation of the solid Earth can reflect ice mass variations, but GNSS observations are also sensitive to changes in past ice sheets through glacial isostatic adjustment (GIA) (e.g. Peltier et al. 2015; Caron & Ivins 2020). As estimating GIA in Greenland remains challenging today (Spada et al. 2012; Sutterley et al. 2014; Wake et al. 2016), disentangling the impacts of current ice thinning from ice mass loss since the last glacial maximum (LGM) or during the little ice age (LIA) in GNSS time-series is difficult (Kjeldsen et al. 2015; Adhikari et al. 2021).

In this study, we leverage a combination of different geodetic techniques, presented in Section 2, to overcome some of the inherent limitations of each technique in estimating Greenland IMB. In Section 3, we focus on critically evaluating three different FDMs to assess uncertainties in IMB estimates derived from changes in ice elevation and subsequently compare these estimates with those obtained from gravity measurements. We then combine SEC observations and gravity variations to derive an optimized IMB estimate that takes advantage of the high spatial resolution of altimetry and optical stereo-imagery with the more direct mass measurement of satellite gravimetry. This combined estimate allows us to discuss the uncertainties introduced by FDMs in the context of IMB derived from changes in ice elevation in Greenland. Furthermore, in Section 4, we calculate the elastic deformation of the solid Earth resulting from ice mass variations derived from the different estimated IMBs, at very high spatial resolution during the two time periods considered. By examining the consistency between different ice mass variations and GNSS observations across Greenland,

with respect to various glacial isostatic adjustment (GIA) models, we aim to assess the validity of IMB methods against independent observations. Finally, in Section 5, we discuss our results and conclude.

# 2 DATA AND THEIR UNCERTAINTIES

## 2.1 Ice surface elevation changes

We estimate ice SEC over the Greenland ice sheet from satellite altimetry observations during two distinct periods: (1) from October 2003 to March 2008 using data provided by the ICESat satellite mission (Sørensen et al. 2011) and (2) from January 2011 to December 2015 using data provided by the CryoSat-2 satellite mission (Simonsen & Sørensen 2017; Sørensen et al. 2018). To complement these observations and ensure comprehensive coverage, we integrate observations of ice elevation changes over Greenland's peripheral glaciers using the stereophotogrammetry-based global glacier data base provided by Hugonnet et al. (2021) over the same periods. Linear trends in ice SEC, derived from satellite altimetry, were provided as grids with spatial resolutions of 5 × 5 km² for ICESat and 15 × 15 km² for CryoSat-2. Data for the peripheral glaciers were provided at a higher spatial resolution images of 0.1 × 0.1 km². In cases of overlap between satellite altimetry data and peripheral glacier observations, likely due to different definitions of what constitutes peripheral glaciers, we prioritized the latter due to its higher spatial resolution. The resulting linear trends in ice elevation change, reflecting the integrated data set, are shown in Fig. 1 over Greenland for the two periods of interest.

Ice SEC can be converted to ice mass changes using this formula:

$$ \frac{dM_{SEC}}{dt} = \left( \frac{dH_{Obs}}{dt} - \frac{dH_{Elas}}{dt} - \frac{dH_{GIA}}{dt} - \frac{dH_{FAC}}{dt} \right) \times \rho_{ice}, \quad (1) $$

where $\frac{dM_{SEC}}{dt}$ is the local rate of ice mass change per unit of surface, that is, the local IMB, as inferred from SEC data. $\rho_{ice}$ is the density of pure ice (917 kg m⁻³). $\frac{dH_{Obs}}{dt}$ is the observed trend in SEC as shown in Figs 1(a) and (b). $\frac{dH_{Elas}}{dt}$ and $\frac{dH_{GIA}}{dt}$ are rates of changes in the bedrock elevation due, respectively, to elastic deformation and GIA. $\frac{dH_{FAC}}{dt}$ is the elevation rate correction that has to be applied because of presence of air within the ice column (FAC for 'Firn Air Content'). The presence of this last term is a consequence of physical processes occurring within the firn layer, which are discussed below in Section 3.1.

As shown in Figs 1(a) and (b), $\frac{dH_{Obs}}{dt}$ can locally reach a few m yr⁻¹. It should be noted that the colour bars are saturated to improve the readability of the figures. The observed maximum rates (in absolute) reach up to −15.4 m yr⁻¹ over 2003–2008 and −10.7 m yr⁻¹ over 2011–2015. According to errors provided with the data set, and to dedicated publications, the uncertainties in $\frac{dH_{Obs}}{dt}$ estimations are typically of the order of a few tenths of cm yr⁻¹. Considering these values and ice SEC observations, if we assume no ground deformation ($H_{Elas} = H_{GIA} = 0$) and no air in the ice column ($H_{FAC}$), we would obtain total IMB rates over Greenland of −336.0 Gt yr⁻¹ during the 2003–2008 period and −296.4 Gt yr⁻¹ during the 2011–2015 period, with an uncertainty of approximately 25–35 Gt yr⁻¹ over both periods (Sørensen et al. 2011).

The terms $\frac{dH_{Elas}}{dt}$ and $\frac{dH_{GIA}}{dt}$ are very small compared with $\frac{dH_{Obs}}{dt}$. As discussed in Sections 3.2 and 5, these terms induce changes of less than a few cm yr⁻¹ changes to the SEC rate. Consequently, they are generally neglected in the SEC estimations for calculating
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Ice SEC estimated from altimetry data, peripheral glaciers information and load variations in equivalent water height (EWH) as seen by GRACE space gravity mission over Greenland: (a) SEC linear trend (in m yr⁻¹) estimated from ICESat data from October 2003 to March 2008 (Sørensen *et al.* 2011), (b) SEC linear trend (in m yr⁻¹) estimated from CryoSat-2 data from January 2011 to December 2015 (Simonsen & Sørensen 2017; Sørensen *et al.* 2018), (c) and (d) load linear trends in EWH (cm yr⁻¹) from GRACE filtered geoid variations during the same 2003–2008 and 2011–2015 periods, respectively (Gauer *et al.* 2023).</figcaption></figure>

the IMB. In Section 3.2, we estimate that the contribution of $\frac{dH_{GIA}}{dt}$ to be less than 0.5 per cent of Greenland's total IMB solely from $\frac{dH_{Obs}}{dt}$, which is one order of magnitude smaller than uncertainties mentioned before. Therefore, in this study we neglect the $\frac{dH_{GIA}}{dt}$ and $\frac{dH_{Elas}}{dt}$ terms in eq. (1).

## 2.2 Satellite gravimetry data

Our study extends to measurements of the space and time varying Earth's gravity field recorded by the GRACE satellite mission (e.g. Tapley *et al.* 2004), over the same periods. In particular, we use the solution provided by Gauer *et al.* (2023), based on previous
work by Prevost et al. (2019). This solution uses multichannel singular spectrum analysis (M-SSA) (Ghil et al. 2002) to exploit spatiotemporal information of multiple Level-2 GRACE/GRACE-FO solutions, expressed in Stokes coefficients of the gravitational potential spherical harmonics decomposition, filtered using the DDK7 decorrelation and complementary Lobe-Edge filter (Kusche 2007; Kusche et al. 2009; Gauer et al. 2023), to fill observational gaps and mitigate north-south striping artefacts caused by instrumental and processing errors. The solution is corrected for non-tidal high-frequency atmospheric and oceanic mass variation models, namely the Atmosphere and Ocean Dealiasing Level-1B (AOD1B) model (Dobslaw et al. 2017). Gravity linear trends were estimated from the monthly GRACE solutions in order to obtain long-term mass changes over Greenland for periods of interest. Note that the large spatial wavelength GRACE-derived estimates over Greenland might capture signal due to other processes, including ice thinning from adjacent regions like North Canada and Iceland, but also worldwide hydrologic mass redistribution. If we want to compare GRACE and the altimetry IMBs, we need to correct for these contributions in the mass variations derived from GRACE over Greenland. To mitigate the impact of these signals, we evaluated the contribution of ice thinning (and other possible mass variations) from outside Greenland on our GRACE-derived IMB using a mascon solution (GRACE CSR solution; Save et al. 2016). We masked the ice loads over Greenland, then expanded the global geographical grid into spherical harmonics and filtered this external signal to achieve the same resolution as our primary GRACE M-SSA data set. We finally corrected our GRACE M-SSA solution over Greenland for this potential signal leakage. Note that this signal is very small and only slightly reduces IMB estimates over the northwestern regions of Greenland.

The corrected linear trends of GRACE M-SSA geoid changes are presented in Fig. 1. The trends are here expressed in equivalent water height (EWH) per year, that is, as if all the gravity signal was only due to an equivalent water load on the Earth's surface, accounting for ground elastic deformation. This is a classical choice of units, however, in reality, the IMB rate as observed by GRACE measurements is not only due to water or ice mass variation. It also results from GIA perturbations. We may express the local mass trends as observed by GRACE as follows:

$$ \frac{dM_{\text{Grav}}}{dt} = \left( \frac{dH_{\text{GRACE}}^{\text{ewh}}}{dt} - \frac{dH_{\text{GIA}}^{\text{ewh}}}{dt} \right) \rho_{\text{water}}, \quad (2) $$

where $\frac{dM_{\text{Grav}}}{dt}$ is the local rate of ice mass change per unit of surface, that is, the local IMB, as inferred by gravity data. $\frac{dH_{\text{GRACE}}^{\text{ewh}}}{dt}$ is the local GRACE trend expressed in EWH, while $\frac{dH_{\text{GIA}}^{\text{ewh}}}{dt}$ is a correction term due to GIA gravity impact, expressed also in EWH. As we will see below this correction term is small but not negligible (Section 3.2). $\frac{dM_{\text{Grav}}}{dt}$ and $\frac{dM_{\text{SEC}}}{dt}$ quantities would be directly comparable if their spatial resolutions were similar. However, GRACE data have a relatively low spatial resolution, imposing to filter SEC observations to enable a meaningful comparison. This filtering process is described in Section 4.

As mentioned before, the GRACE M-SSA gravity changes investigated in this study were obtained from multiple classical GRACE solutions over the same time period. Discrepancies between these solutions over Greenland provide an estimate of the uncertainty associated with the GRACE data. Integrating all solution discrepancies with respect to the M-SSA solution over Greenland, we obtain GRACE IMB standard deviations of 11.3 Gt yr⁻¹ during the 2003–2008 period, and 14.9 Gt yr⁻¹ during the 2011–2015 period.

## 2.3 GNSS data

Our investigation relies on data from a network of GNSS stations, with their positions shown in Fig. 2. Note that only six GNSS stations were available during the period 2003–2008, while 51 stations recorded the solid Earth deformation during the period 2011–2015. We use vertical daily GNSS station position time-series processed by the Nevada Geodetic Laboratory (NGL) at the University of Nevada (Blewitt et al. 2018). NGL uses the GipsyX 1.0 software (Bertiger et al. 2020) from the Jet Propulsion Laboratory (JPL) to compute position time-series based on single-station precise point positioning (PPP) with carrier phase ambiguity resolution (Zumberge et al. 1997). Satellite orbit and clock data used for the PPP processing are JPL's Repro 3.0 products. All products are based on, and solutions are aligned to, the IGS14 realization (Rebischung et al. 2016) of the ITRF14 reference frame (Altamimi et al. 2016). The VMF1 mapping function with elevation-dependent data weight, and with nominal dry and wet zenith delays from the numerical weather model of the European Centre for Medium-Range Weather Forecasts (Boehm et al. 2006) was used to correct tropospheric effects.

For each GNSS time-series, we adjusted a kinematic model which includes a piece-wise linear function accounting for linear trends and offsets (equipment change, etc.), and two sine functions accounting for annual and semi-annual signals. In some cases, we added discontinuities in velocity to adjust for changes in velocity slopes. In order to obtain realistic uncertainties, the adjustment was performed assuming a variable white plus power-law noise model in our time-series, simultaneously with the kinematic model. Such a noise model has been shown to describe well the background noise in GNSS station position time-series (e.g. Mao et al. 1999; Santamaría-Gómez et al. 2011; Gobron et al. 2021). Noise model parameters were estimated by restricted maximum-likelihood, an unbiased alternative to classical maximum-likelihood estimation (Patterson & Thompson 1971; Harville 1977; Gobron et al. 2022). Finally, we extracted from this analysis vertical velocities with associated realistic uncertainties for each station over each period considered (Fig. 2).

GNSS observation are point-based, scattered and primarily located around the edges of the ice sheet, making it impossible to infer Greenland IMB from GNSS data alone. Instead, we use GNSS observations here to validate and assess IMB inferred from other data sets. Observed GNSS vertical velocities in Greenland are mainly driven by elastic deformation induced by recent changes in the ice load and viscoelastic deformation induced by GIA processes.

# 3 CORRECTING MODELS AND THEIR IMPACTS ON IMB ESTIMATES

We used two types of models for correcting data: FDM and GIA models. Both types of models are complex and based on various assumptions and parameters. It is therefore difficult to evaluate the uncertainties coming from these models. For this reason we chose to test multiple models of various origins. It offers the possibility to evaluate discrepancies and biases that such correcting models typically imply in IMB estimations.

## 3.1 Firn densification models

Within ice sheets and glaciers, an ice column is classically divided into two parts (e.g. Ligtenberg et al. 2011): (1) the underlying ice layer, where density is constant and equal to the density of pure ice
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) GNSS network used in this study (Blewitt *et al.* 2018) and (b) associated estimated vertical velocities over the 2011–2015 period. Data for only six of these stations were available over the 2003–2008 period, and velocities are shown by black boxes. Note that the colours refer only to the location of stations and are consistent with Figs 8 and 9.</figcaption></figure>

$\rho_{\text{ice}} = 917 \text{ kg m}^{-3}$, and (2) the firn layer, at the top, where density is generally lower than $\rho_{\text{ice}}$ and varies with depth and time due to the presence of air and to compaction processes. The thickness of the ice layer can vary due to the basal ice dynamics of the ice sheet, while the firn layer evolves according to surface processes: snowfall, erosion, compaction, etc.

Understanding the variations in density within Greenland's ice column, particularly within the surface firn layer, is necessary for accurately inferring the IMB from changes in ice surface elevation as measured by satellite altimetry or stereophotogrammetry. Indeed, in the firn layer, air content and compaction can induce significant density variations (Cuffey 2008).

We investigated three FDMs, each rooted in different climate and compaction models. The first model, which we have named FDM-IMAU (provided by Institute for Marine and Atmospheric research Utrecht), is forced by the climate model Regional Atmospheric Climate Model (RACMO2) (Noël *et al.* 2016, 2018), and uses the compaction model developed by Kuipers Munneke *et al.* (2015), Ligtenberg *et al.* (2011) and Ligtenberg *et al.* (2018). The second model, named FDM-GSFC (provided by NASA Goddard Space Flight Center), is based on the GSFC-FDMv1 firn densification model developed by Medley *et al.* (2022), which is forced by MERRA-2 atmospheric re-analysis model (Gelaro *et al.* 2017). Finally, the third model, named FDM-JPL (provided by NASA Jet
propulsion Laboratory), has been developed using JPL's Ice Sheet and Sea-level System Model (ISSM) (Larour et al. 2012) and is forced by ERA5 atmospheric re-analysis model (Hersbach et al. 2017). Data from FDM-IMAU were provided as grids of 1 × 1 km² resolution over the period 1958–2017, while FDM-GSFC and FDM-JPL data were provided as 12.5 × 12.5 km² grids during the period 1980–2021.

We estimated local IMB in Greenland from trends in ice elevation changes, assuming the different FDM corrections. Considering assumptions mentioned in Section 2.1 and eq. (1), the IMB at a given location and date, as seen by SEC, can be simply calculated as follows:

$$ \frac{dM_{SEC}}{dt} = \left( \frac{dH_{Obs}}{dt} \right) \times \rho_{ice} + \Delta \dot{M}_{Firn}, \quad (3) $$

where $\Delta \dot{M}_{Firn}$ is a firn correction term accounting for all processes occurring within the firn layer, which is directly related to $\frac{dH_{FAC}}{dt}$ quantity:

$$ \Delta \dot{M}_{Firn} = -\frac{dH_{FAC}}{dt} \times \rho_{ice}. \quad (4) $$

As previously mentioned $H_{FAC}$ is the firn air content expressed in metres. This quantity, sometimes called dip-integrated porosity, represents the volume of air integrated over the depth of a firn column of one square metre base. It can be calculated as follows: (Medley et al. 2022; Stevens et al. 2020):

$$ H_{FAC} = \int_{0}^{H_{F}} \frac{\rho_{ice} - \rho(z)}{\rho_{ice}} dz, \quad (5) $$

where $H_F$ is the thickness of the firn layer and $\rho(z)$ is the density at a given depth $z$ within the firn layer. $H_{FAC}$ is provided by the different FDM. We calculated a trend of firn air content for each period of interest and then calculated the firn correction term for each FDM. In a context of decreasing ice column ($\frac{dH_{Obs}}{dt} < 0$), a positive firn correction $\Delta M > 0$ means that a part of the decrease in ice height affects the firn layer, while a negative firn correction $\Delta M < 0$ implies that the firn layer tends to increase even if the total column decreases. In the latter case, the underlying ice layer is decreasing at a faster rate than the observed change in surface elevation would suggest.

Fig. 3 shows $\Delta \dot{M}_{Firn}$ for our three different FDM during the two periods considered. During the period 2003–2008, all IMB firn correction terms show very high values near the southeast coasts of Greenland (Figs 3a–c). This means that the air content in the firn layer decreases. Consequently, although these regions were subject to a particularly large negative elevation trend (see Fig. 1a), the thinning of the ice layer was less than might have been expected from the SEC data during the period, probably due to significant compaction processes. The three models nevertheless show significantly different spatial distributions and amplitudes of IMB firn corrections, particularly for the FDM-IMAU models close to the western coasts.

The total IMB corrections across Greenland, as determined by the different FDMs, exhibit significant discrepancies over the 2003–2008 period: +116.4 Gt yr⁻¹ for FDM-IMAU, +69.5 Gt yr⁻¹ for FDM-GSFC, +109.8 Gt yr⁻¹ for FDM-JPL. During the 2011–2015 period, the IMB firn correction terms differ significantly from those of the 2003–2008 period, also showing notable discrepancies across various firn models. In particular, FDM-IMAU shows positive IMB correction trends along the entire southern and western coasts, whereas FDM-GSFC and FDM-JPL show negative trends near the southern coasts. While the latter two models exhibit greater

consistency with each other, notable differences in amplitudes are evident, especially along the northwestern coast. The total IMB corrections over Greenland for the various FDM over the period 2011–2015 are +80.8 Gt yr⁻¹ for FDM-IMAU, -2.7 Gt yr⁻¹ for FDM-GSFC, +44.5 Gt yr⁻¹ for FDM-JPL. While the total firn corrections appear negligible using the FDM-GSFC model over this period, it would be significant using the FDM-JPL model or, even more so, the FDM-IMAU model. Thus, while ignoring firn processes would erroneously lead to large negative IMB over Greenland, particularly over the 2003–2008 period, firn models show significant differences that may lead to large uncertainties in IMB derived from ice elevation changes.

To obtain IMB, we combine the firn corrections with SEC derived from satellite altimetry data (ice sheet) and stereophotogrammetry (peripheral glaciers) for the two periods (see Figs 1a and b). The total Greenland IMB obtained are shown in Table 1.

Differences between FDM corrections are therefore extremely significant and can lead to an error of 20–40 per cent in the IMB estimates. This underscores the importance of selecting or developing an accurate FDM when using ice elevation changes to infer IMB.

## 3.2 GIA models

We investigate five GIA models: (1) ICE-6G model developed by Peltier et al. (2015), (2) ANU model from the Australian National University (Lambeck et al. 2014), (3) Caron model from the Jet Propulsion Laboratory (JPL) model and developed by Caron et al. (2018), (4) Huy3 model developed by Lecavalier et al. (2014) and (5) LM17.3 model proposed by Steffen et al. (2021).

The ICE-6G model is a global model that assumes a radially symmetric Earth structure [using the VM5a viscosity profile provided by Peltier et al. (2015)] and has been constrained by paleo-sea level data and recent geodetic observations. The ANU model is also a global model that has been constrained by paleo-sea level data, but this model assumes different viscosity profiles depending on the regions of interest (North America, Europe, Antarctica, remote regions, etc.) (Lambeck et al. 2014). The Caron model provides a scaled version of the global ANU ice history. The scaling factors (according to the regions) and a single viscosity profile have been jointly inverted to better explain the GNSS and gravity observations (Caron et al. 2018). The Huy3 model is a regional model based on glacial and paleo-sea level records in Greenland. A viscosity profile was derived to optimally explain the observations (Lecavalier et al. 2014). Finally, the LM17.3 model is a global/regional model in which the ice history is a combination of previously published ice history models (Steffen et al. 2021). In Greenland, the ice history of the LM17.3 model is the same as that of the Huy3 model, but a different viscosity profile (VM5a profile) is used to calculate the solid Earth’s response to mass redistribution.

Each GIA model provides the solid Earth response to the ice retreat that followed the last glacial period, from which can be extracted present-day surface ground velocity and the associated gravity perturbation that we express here in equivalent water height (Fig. 4). Fig. 4 shows ground vertical velocities in mm yr⁻¹ (top row) and gravity variations expressed in mm yr⁻¹ EWH (bottom row) predicted by the different GIA models. Note that the LM17.3 model only provides vertical velocities.

As expected, the ANU and Caron models show similar patterns in terms of ground vertical velocities and gravitational perturbations (Fig. 4). They show significant positive signals in the northern,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Firn mass correction terms [as defined in eqs (4) and (5), in Section 3.1], during the two periods of interest, 2003–2008 (a–c) and 2011–2015 (d–f), associated with three different firn densification models (FDM): (a) and (d) FDM-IMAU, (b) and (e) FDM-GSFC, (c) and (f) FDM-JPL. The mass balance correction induced by each FDM on the entire Greenland is indicated on top of each figure.</figcaption></figure>

<table><caption>**Table 1.** Total ice mass balance over Greenland estimated from SEC (based on ICESat data for the 2003–2008 period, and CryoSat-2 data for the 2011–2015 period, both completed with stereophotogrammetry data on peripheral glaciers) and assuming different firn densification models (FDM).</caption><thead><tr><th>Period</th><th>FDM-IMAU</th><th>FDM-GSFC</th><th>FDM-JPL</th></tr></thead><tbody><tr><td>2003–2008</td><td>-215.0 Gt yr<sup>-1</sup></td><td>-261.9 Gt yr<sup>-1</sup></td><td>-221.5 Gt yr<sup>-1</sup></td></tr><tr><td>2011–2015</td><td>-212.0 Gt yr<sup>-1</sup></td><td>-295.4 Gt yr<sup>-1</sup></td><td>-248.6 Gt yr<sup>-1</sup></td></tr></tbody></table>

eastern and southwestern regions, although stronger in the Caron model. Similarly, Huy3 and LM17.3 models show similarities in terms of vertical velocities but are far from identical. In particular, Huy3 models show significant negative vertical velocities in the

southwestern region that are accompanied by strong gravity variations. The ICE-6G model shows very strong Earth response signals in the northeastern region, and weaker ones in the western and eastern regions. Unlike the others (except possibly LM17.3), it shows significant ground motion at the southernmost point.

While gravity predictions will be used to correct GRACE data, the vertical velocities predictions will be used to investigate GNSS observations. They may also be used to correct altimetry data ($\frac{dH_{GIA}}{dt}$ term in eq. 1). However, if we integrate the vertical velocities over Greenland and multiply by $\rho_{ice}$, it leads to changes in IMB between -0.9 and 1.6 Gt yr$^{-1}$, depending on the model. These estimations are very small and represent less than 0.5 per cent of the IMB values that can be inferred from altimetry SEC alone (see Section 2.1). For these reasons we will not correct SEC observations for GIA impact.
<figure><img src="image_6.png" /><figcaption>Figure 4. Vertical ground velocity (top row) and gravity perturbation (bottom row) for different GIA models over Greenland: ICE-6G, ANU, Caron, Huy3 and LM17.3 (velocity only for LM17.3). Ground vertical velocities and gravity perturbations (expressed in equivalent water height) are in mm yr⁻¹.</figcaption></figure>

## 4 COMBINING GRAVITY AND SEC MEASUREMENTS FOR IMB ESTIMATION: EQUIVALENT FIRN CORRECTIONS

One way of discriminating between IMB models derived from surface elevation observations and FDMs, is to compare results with independent observations from satellite gravimetry which represent a more direct measurement of mass variations despite their coarser spatial resolution. In this context, we used data from the GRACE mission (see Section 2.2) to derive maps of equivalent water height variations across Greenland over the two periods of interest (Figs 1c and d). The GRACE solution has been corrected for GIA over Greenland, testing the four GIA models that propose gravity trend predictions (ICE-6G, ANU, Caron, Huy3 models, see Section 3.2). Given the significantly coarser resolution of GRACE data compared to that of SEC, we expanded the IMB solutions derived from the SEC into spherical harmonics up to the maximum degree of our GRACE solution (89). We then filtered the resulting spherical harmonic coefficients with a DDK7 (Kusche et al. 2009) and lobe-edge filter, consistent with Gauer et al. (2023), and finally reconstructed the filtered SEC signal in the spatial domain, to ensure comparison with the GRACE M-SSA solution at a similar spatial resolution. This method enables us to validate the different IMB solutions we have obtained, depending on FDM, against the GRACE observations.

In addition to FDM analysis, we developed a method for determining an ‘equivalent firn model’ that optimally reconciles IMBs inferred from SEC and gravimetry observations. We chose to calculate scaling factors, individually tailored to distinct regions across Greenland. We identified 12 different regions based on the drainage basins proposed by Zwally et al. (2012) and Rignot & Mouginot (2012). Since the ice dynamics at the centre of the ice sheet generally

differs from its edges, we divided each basin in two distinct regions, one internal and the other external (coastal). The defined regions are shown in Fig. 5, and are labelled: ‘CW’ for centre west, ‘NC’ for north centre, ‘NE’ for northeast, ‘NW’ for northwest, ‘SW’ for southwest, ‘SE’ for southeast, ‘CE’ for central east, ‘CS’ for central south, ‘CWI’ for central west interior, ‘NEI’ for northeast interior, ‘NWI’ for northwest interior, ‘NCI’ for north central interior. We also considered a few glaciers located outside Greenland, along the Canadian coasts, in order to account for the possible leakage of this signal over Greenland in GRACE observations due to its poor spatial resolution (see Section 2.2 for details on the estimation of this signal). We have defined a single scaling factor for this ‘Outer region’, which is not shown in Fig. 5. For each basin, we used SEC derived from satellite altimetry data and the peripheral glaciers data base, which we multiplied by 1000 (as if the column load were composed of water). We applied the filters mentioned above on this regional ice load to ensure consistency with our GRACE solution. Filtering causes signal from a single basin to smear over other areas, which needs to be accounted for. We use the following relationship for a given point in region j:

$$ \frac{dM_{\text{Grav}}^j}{dt} = \sum_{i=1}^{N_{\text{region}}} \left( \frac{d\widetilde{H}_{\text{SEC}}^i}{dt} \right) \times 1000 \times c_i, \quad (6) $$

where $\frac{dM_{\text{Grav}}^j}{dt}$ is the local mass variation as seen by the GRACE solution in region j, $\frac{d\widetilde{H}_{\text{SEC}}^i}{dt}$ is the signal from the observed ice SEC in region i that modifies the mass trend at the considered point in region j due to filtering and $c_i$ the scaling factor for the region i. As seen above, the ~ symbol indicates that SEC observations have been filtered at the GRACE resolution using the appropriate filters. Each scaling factor can be considered as an equivalent mean ice
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Scaling factors $c_i$ to be applied to SEC to best match GRACE observations: (a) for 2003–2008, and (b) for 2011–2015. Each scaling factor can be interpreted as an effective mean density of the ice variation in the region (g cm⁻³), reflecting the impact of firn processes (see the text for region acronyms). Scaling factors in the interior regions are omitted here due to insignificance.</figcaption></figure>

density expressed in g cm⁻³ (or kg m⁻³ if multiplied by 1000), which should be used in the region to optimally explain the local IMB over Greenland determined from GRACE using regional SEC. It is important to understand that, due to firn processes, this equivalent density can be very different from the true ice/snow density, and can be greater than 1 (water density in g cm⁻³) or even negative. The relationship between the scaling factor and the firn correction term defined in eqs (3) and (4) is relatively straightforward. If we apply the GRACE spatial filtering to eq. (3), we can consider that the mass trends in eqs (3) and (6) are equivalent, and we have:

$$ \frac{\mathrm{d}\tilde{H}_{\mathrm{SEC}}}{\mathrm{d}t} \rho_{\mathrm{ice}} + \Delta \tilde{M}_{\mathrm{Firn}} = \left( \frac{\mathrm{d}\tilde{H}_{\mathrm{SEC}}}{\mathrm{d}t} \right) 1000 c. \qquad (7) $$

Therefore, we have the following relationship:

$$ \begin{aligned} c &= \frac{1}{1000} \left( \rho_{\mathrm{ice}} + \frac{\Delta \tilde{M}_{\mathrm{Firn}}}{\frac{\mathrm{d}\tilde{H}_{\mathrm{SEC}}}{\mathrm{d}t}} \right) \\ &= 0.917 + \frac{1}{1000} \frac{\Delta \tilde{M}_{\mathrm{Firn}}}{\frac{\mathrm{d}\tilde{H}_{\mathrm{SEC}}}{\mathrm{d}t}} \end{aligned} \qquad (8) $$

Note that the value of the scaling factor strongly depends on physical processes occurring within the firn layer. If there are no changes in the firn layer (no variations in thickness or density), then

the scaling factor would be 0.917, that is, equal to the density of the ice in g cm⁻³. If the firn correction term ΔM<sub>Firn</sub> and the SEC trend have opposite signs, the scaling factor will be smaller than 0.917. If it remains larger than 0, it indicates that changes occurred in the firn layer and that thicknesses of the firn layer and of the underlying ice layer have evolved in the same direction as SEC. If SEC and ΔM<sub>Firn</sub> trends have the same sign, the scaling factor will exceed 0.917, and may be greater than 1. This indicates that thicknesses of the firn layer and of the ice underlying layer have evolved in opposition but changes in the ice layer dominate those in the firn layer. If the scaling factor is negative it indicates also that thicknesses of the firn layer and of the ice underlying layer have evolved in opposition but this time the firn layer has changed more dramatically than the underlying ice layer. For example, if we gain one metre of firn with a density of 0.5 g cm⁻³ but lose 90 cm of ice with a density of 0.917 g cm⁻³ at greater depth, there will be a 10 cm increase in ice surface topography but an overall loss of ice mass of −0.325 × 10³ kg m⁻² within the ice column. In this case, the scaling factor would be −3.25.

We estimated regional scaling factors $c_i$ that best explained the GRACE observations by a least-squares inversion for all basins and for the signal outside Greenland, considering various GIA corrections. Table 2 shows results for different GIA models, with uncertainties derived from the GRACE solution uncertainty (see Section 2.2).
<table><caption>**Table 2.** Scaling factor for different Greenland regions (see the text for details of the acronyms), corresponding to equivalent mean densities in $g cm^{-3}$, estimated from a combination of SEC and gravity data and depending on the GIA model used to correct gravity data.</caption><thead><tr><th colspan="6">Scaling factor</th></tr><tr><th colspan="6">2003–2008</th></tr><tr><th>Region</th><th>No GIA correction</th><th>ICE-6G</th><th>GIA model<br>ANU</th><th>Caron</th><th>Huy3</th></tr></thead><tbody><tr><td>NE</td><td>0.280 ± 0.011</td><td>0.366 ± 0.011</td><td>0.342 ± 0.011</td><td>0.317 ± 0.011</td><td>0.316 ± 0.011</td></tr><tr><td>NC</td><td>0.073 ± 0.012</td><td>0.249 ± 0.012</td><td>0.237 ± 0.012</td><td>0.153 ± 0.012</td><td>0.183 ± 0.012</td></tr><tr><td>NW</td><td>0.575 ± 0.004</td><td>0.572 ± 0.004</td><td>0.579 ± 0.004</td><td>0.556 ± 0.004</td><td>0.568 ± 0.004</td></tr><tr><td>CW</td><td>0.752 ± 0.013</td><td>0.720 ± 0.013</td><td>0.767 ± 0.013</td><td>0.710 ± 0.013</td><td>0.758 ± 0.013</td></tr><tr><td>SW</td><td>0.586 ± 0.016</td><td>0.601 ± 0.016</td><td>0.618 ± 0.016</td><td>0.474 ± 0.016</td><td>0.620 ± 0.016</td></tr><tr><td>SE</td><td>0.562 ± 0.004</td><td>0.551 ± 0.004</td><td>0.561 ± 0.004</td><td>0.555 ± 0.004</td><td>0.559 ± 0.004</td></tr><tr><td>CE</td><td>0.839 ± 0.005</td><td>0.845 ± 0.005</td><td>0.852 ± 0.005</td><td>0.850 ± 0.005</td><td>0.850 ± 0.005</td></tr><tr><td>NEI</td><td>0.373 ± 0.067</td><td>0.361 ± 0.067</td><td>0.303 ± 0.067</td><td>0.379 ± 0.067</td><td>0.338 ± 0.067</td></tr><tr><td>NCI</td><td>0.443 ± 0.177</td><td>0.252 ± 0.177</td><td>-0.153 ± 0.177</td><td>0.167 ± 0.177</td><td>0.002 ± 0.177</td></tr><tr><td>NWI</td><td>1.588 ± 0.347</td><td>2.320 ± 0.347</td><td>1.667 ± 0.347</td><td>1.871 ± 0.347</td><td>1.497 ± 0.347</td></tr><tr><td>CWI</td><td>0.474 ± 0.056</td><td>0.767 ± 0.056</td><td>0.518 ± 0.056</td><td>0.481 ± 0.056</td><td>0.492 ± 0.056</td></tr><tr><td>CS</td><td>0.434 ± 0.064</td><td>0.384 ± 0.064</td><td>0.157 ± 0.064</td><td>0.397 ± 0.064</td><td>0.240 ± 0.064</td></tr></tbody><tbody><tr><th colspan="6">2011–2015</th></tr><tr><th>Region</th><th>No GIA correction</th><th>ICE-6G</th><th>GIA model<br>ANU</th><th>Caron</th><th>Huy3</th></tr><tr><td>NE</td><td>0.467 ± 0.011</td><td>0.630 ± 0.011</td><td>0.529 ± 0.011</td><td>0.515 ± 0.011</td><td>0.501 ± 0.011</td></tr><tr><td>NC</td><td>0.450 ± 0.009</td><td>0.559 ± 0.009</td><td>0.564 ± 0.009</td><td>0.507 ± 0.009</td><td>0.533 ± 0.009</td></tr><tr><td>NW</td><td>0.748 ± 0.008</td><td>0.756 ± 0.008</td><td>0.751 ± 0.008</td><td>0.739 ± 0.008</td><td>0.749 ± 0.008</td></tr><tr><td>CW</td><td>0.816 ± 0.013</td><td>0.785 ± 0.013</td><td>0.844 ± 0.013</td><td>0.803 ± 0.013</td><td>0.836 ± 0.013</td></tr><tr><td>SW</td><td>0.965 ± 0.025</td><td>1.022 ± 0.025</td><td>1.001 ± 0.025</td><td>0.908 ± 0.025</td><td>1.003 ± 0.025</td></tr><tr><td>SE</td><td>1.433 ± 0.039</td><td>1.423 ± 0.039</td><td>1.385 ± 0.039</td><td>1.397 ± 0.039</td><td>1.398 ± 0.039</td></tr><tr><td>CE</td><td>0.950 ± 0.027</td><td>0.971 ± 0.027</td><td>0.971 ± 0.027</td><td>0.971 ± 0.027</td><td>0.969 ± 0.027</td></tr><tr><td>NEI</td><td>0.375 ± 0.041</td><td>0.468 ± 0.041</td><td>0.291 ± 0.041</td><td>0.413 ± 0.041</td><td>0.310 ± 0.041</td></tr><tr><td>NCI</td><td>0.426 ± 0.375</td><td>1.088 ± 0.375</td><td>0.116 ± 0.375</td><td>0.350 ± 0.375</td><td>0.216 ± 0.375</td></tr><tr><td>NWI</td><td>-0.404 ± 0.147</td><td>-0.128 ± 0.147</td><td>-0.123 ± 0.147</td><td>-0.158 ± 0.147</td><td>-0.153 ± 0.147</td></tr><tr><td>CWI</td><td>-0.465 ± 0.191</td><td>-0.125 ± 0.191</td><td>-0.580 ± 0.191</td><td>-0.384 ± 0.191</td><td>-0.563 ± 0.191</td></tr><tr><td>CS</td><td>-0.650 ± 0.097</td><td>-0.791 ± 0.097</td><td>-0.537 ± 0.097</td><td>-0.607 ± 0.097</td><td>-0.586 ± 0.097</td></tr></tbody></table>

It should first be noted that all inner regions (NCI, NEI, NWI, CWI and CS regions) present larger uncertainties. We calculated scaling factors for these regions primarily to ensure the overall consistency of the calculations. However, these values are not significant overall, as the SEC rates observed in these regions are very close to zero. The resulting scaling factors may therefore be anomalously large or small, but in any case, this will not affect the total IMB for Greenland. We have chosen to include these values in the table, but we will mostly focus our attention below on the coastal regions (NC, NE, CE, SE, SW, CW, NW regions). In these regions, we find that the impact of the different GIA models is relatively small. Although the GIA models differ, this can be explained by the fact that gravity variations are less sensitive to GIA than to the large mass variations associated with ice thinning. Therefore, we calculated an average scaling factor and a standard deviation from scaling factors obtained with the different GIA models, for each region. These standard deviations reflect biases and uncertainties introduced by the GIA models on the scaling factor. But since we only have four GIA models, we chose to consider that the error of the mean scaling factors as three times the standard deviation in order to encompass the entire discrepancy between GIA models. The mean scaling factors are presented in the first column of Table 3. Uncertainties in this table are the sum of uncertainties coming from gravity data and the error from GIA models. As noted earlier, the inner regions present very large uncertainties compared with coastal regions, illustrating the poor significance of their scaling factor. Fig. 5 shows the value of the mean scaling factors over Greenland coastal regions for both periods of interest.

For the period spanning 2003 to 2008, all scaling factors in coastal regions are significantly lower than ice density. This indicates that all coastal regions have been affected by thinning of the firn layer. Table 3 also presents the scaling factors computed using the various FDMs. Our estimates from GRACE show relatively higher consistency with FDM-JPL and FDM-GSFC models than with FDM-IMAU across most regions, particularly in CW and SW regions where the FDM-IMAU factors exceed ice density. In NE and NC regions, our estimates differ significantly from all FDM (though they are less divergent from FDM-JPL in the NC region). These two regions are particularly affected by GIA, yet none of the considered GIA models can fully explain this discrepancy. Differences are mostly due to the fact that FDM show firn corrections slightly negative (see Fig. 3), especially in NE regions (which is not the case in NC region for FDM-JPL). Our findings from GRACE indicate firn corrections with opposite signs to those suggested by FDM in these regions.

For the period spanning 2011 to 2015, scaling factors are quite different from those of the previous period. Indeed, we find particularly large scaling factors in the southern regions (SW, SE and CE regions), where they are larger than ice density, even greater than 1 in SE region. In these regions, where we mainly observe negative elevation rates, this means that the firn correction $\Delta M_{\text{Firn}}$ is necessarily negative. This implies that, although ice surface elevations have decreased over the period, the firn layer thickness has increased. This is consistent with estimates from FDM-JPL and partly from FDM-GSFC, but not from FDM-IMAU (Fig. 3 and Table 3). In other regions, scaling factors are lower than the ice density. In these regions, the firn layer thickness evolves in the same direction as the
<table><caption>**Table 3.** Equivalent mean densities (in g cm⁻³) inferred from the different FDM and the scaling model from GRACE in the different Greenland regions (see the text for details of the acronyms and for computation of the densities associated with GRACE corrected from four GIA models).</caption><thead><tr><th colspan="5">2003–2008</th></tr><tr><th>Region</th><th>Scaling model</th><th>FDM-IMAU</th><th>FDM-GSFC</th><th>FDM-JPL</th></tr></thead><tbody><tr><td>NE</td><td>0.335 ± 0.042</td><td>0.831</td><td>1.018</td><td>0.9386</td></tr><tr><td>NC</td><td>0.203 ± 0.069</td><td>1.095</td><td>0.936</td><td>0.693</td></tr><tr><td>NW</td><td>0.566 ± 0.014</td><td>0.890</td><td>0.810</td><td>0.755</td></tr><tr><td>CW</td><td>0.739 ± 0.050</td><td>0.996</td><td>0.757</td><td>0.795</td></tr><tr><td>SW</td><td>0.578 ± 0.107</td><td>1.364</td><td>0.730</td><td>0.763</td></tr><tr><td>SE</td><td>0.557 ± 0.010</td><td>0.329</td><td>0.735</td><td>0.628</td></tr><tr><td>CE</td><td>0.849 ± 0.009</td><td>0.304</td><td>0.735</td><td>0.600</td></tr><tr><td>NEI</td><td>0.345 ± 0.109</td><td>3.886</td><td>3.865</td><td>3.406</td></tr><tr><td>NCI</td><td>0.066 ± 0.411</td><td>1.702</td><td>1.280</td><td>1.791</td></tr><tr><td>NWI</td><td>1.839 ± 0.809</td><td>1.238</td><td>1.170</td><td>1.304</td></tr><tr><td>CWI</td><td>0.564 ± 0.233</td><td>3.908</td><td>2.098</td><td>7.087</td></tr><tr><td>CS</td><td>0.294 ± 0.215</td><td>0.654</td><td>0.987</td><td>1.379</td></tr></tbody><tbody><tr><th colspan="5">2011–2015</th></tr><tr><th>Region</th><th>Scaling model</th><th>FDM-IMAU</th><th>FDM-GSFC</th><th>FDM-JPL</th></tr><tr><td>NE</td><td>0.544 ± 0.088</td><td>0.901</td><td>1.173</td><td>0.930</td></tr><tr><td>NC</td><td>0.541 ± 0.043</td><td>0.689</td><td>0.728</td><td>0.535</td></tr><tr><td>NW</td><td>0.749 ± 0.017</td><td>0.774</td><td>0.839</td><td>0.700</td></tr><tr><td>CW</td><td>0.817 ± 0.049</td><td>0.756</td><td>0.903</td><td>0.842</td></tr><tr><td>SW</td><td>0.984 ± 0.091</td><td>0.845</td><td>1.027</td><td>1.037</td></tr><tr><td>SE</td><td>1.401 ± 0.060</td><td>0.530</td><td>1.051</td><td>1.233</td></tr><tr><td>CE</td><td>0.971 ± 0.028</td><td>0.751</td><td>1.213</td><td>0.978</td></tr><tr><td>NEI</td><td>0.371 ± 0.150</td><td>-0.083</td><td>0.694</td><td>0.612</td></tr><tr><td>NCI</td><td>0.443 ± 0.948</td><td>-2.495</td><td>0.115</td><td>-2.647</td></tr><tr><td>NWI</td><td>-0.140 ± 0.170</td><td>0.915</td><td>1.065</td><td>1.279</td></tr><tr><td>CWI</td><td>-0.417 ± 0.465</td><td>2.284</td><td>1.819</td><td>2.905</td></tr><tr><td>CS</td><td>-0.630 ± 0.241</td><td>2.035</td><td>2.706</td><td>4.276</td></tr></tbody></table>

ice SEC. In NW and CW regions, where we observed particularly high negative elevation rates, the scaling factors are closer to that of ice density. As mentioned above, this means that the thinning affects mainly the ice underlying layer.

Comparisons with the GRACE solution for all models, including our scaling approach, are shown in Figs 6 and 7. None of the FDMs appear capable of perfectly reconciling the GRACE observations with SEC. In Fig. 6, we can see significant positive residuals in the north for all FDMs. Residuals are smaller in the south for FDM-JPL and FDM-GSFC, but remain significant for FDM-IMAU. In Fig. 7, residuals appear slightly lower in the north for all FDMs, but deteriorate in the south for FDM-JPL and FDM-GSFC. Residuals are lower for the scaling factor approach for both periods. This is to be expected, since the method was developed to optimally fit GRACE observations. However, it is not perfect and the residual signal would of course be lower if we had used more regions for our least-squares fit to GRACE observations. Nevertheless, this method offers the possibility of explaining GRACE observations particularly well with SEC data. Using this method, we find global IMB for Greenland in the same range as global IMB predicted from FDM models (compare Tables 1 and 4).

## 5 VALIDATION AGAINST GNSS OBSERVATIONS OF ICE-INDUCED GROUND VERTICAL UPLIFT

Here, we examine whether the various IMB estimates proposed in previous sections align with GNSS station observations of motion during the period 2011–2015. We computed the ground vertical deformation induced by changes in ice loading and compared these results with the GNSS observations. We used a Green’s function

approach to calculate the elastic deformations associated with the different IMB solutions. We used the LOADDEF code developed by Martens et al. (2019) for this purpose. In the calculation, we assumed a PREM Earth internal structure (Dziewonski & Anderson 1981), but modified the density and elastic parameters of the crust. The crustal parameters have been fixed to mean values of CRUST1.0 crustal model (Laske et al. 2013) restricted to Greenland. We chose a Green’s function approach because it offers the possibility of accounting for ice loading at shorter wavelengths than a classical approach using spherical harmonic expansions (down to a resolution of 1–2 km in the case of LOADDEF). Indeed, GNSS motions can be very sensitive to local loads (e.g. Khan et al. 2010). It should be noted that all GNSS stations are located close to the coasts outside of the main ice sheet. For this reason, GNSS motions can be particularly sensitive to mass changes occurring in peripheral glaciers. In fact, our calculations show that the signal coming from the peripheral glaciers can account for half of the ground motion signal recorded at stations close to areas of large ice thickness decrease. This is the case for many stations on southeastern regions, for example, PLPK, KUAQ, KSNB, MIK2, VFDG, LYNS, HEL2, etc. For stations in northwestern regions, it can account for up to 30 per cent of the signal, particularly at DSKG or ASKY stations.

Fig. 8 shows the vertical velocities observed by GNSS stations over the period 2011–2015, together with the modelled vertical velocities induced by the different IMB solutions we have studied. No GIA contribution has been introduced into the vertical velocity predictions in this figure. The first observation that can be made from this figure is that the elastic deformation predictions are rather close to each other, given the uncertainties of GNSS velocities. However, there are significant discrepancies between observed and modelled velocities. These discrepancies are particularly marked for stations
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Predicted gravity changes in equivalent water height in Greenland for the period 2003–2008, estimated from ICESat satellite altimetry data, peripheral glaciers information and assuming different FDMs (IMAU, JPL, GSFC, scaling model from GRACE). The predicted gravity has been filtered as the GRACE data presented Fig. 1; top: obtained trends depending on FDM; bottom: differences between GRACE estimate and top trends.</figcaption></figure>

in certain regions, especially in the southeastern (stations VFDG, KUAQ, KSNB, HEL2, KULU, KBUG) and the northern (stations DSKG, THU, MARG, SCBY, KMOR, ALRT, JWLF, KMJP, JGBL, NRD1, LEFN) regions. They are most likely due to GIA processes, particularly in the north where all GIA models tend to show large vertical velocities (Fig. 4).

We then investigated all the GIA models mentioned in Section 3.2 for surface ground motion corrections (ICE-6G, ANU, Caron, Huy3, and LM17.3 GIA models) to see if they are able to explain the discrepancies observed in Fig. 8 between observations and IMB predictions. Fig. 9 shows vertical velocity differences between observed GNSS velocities and model predictions, which here combine elastic deformation due to recent ice thinning and GIA models. The elastic deformations due to recent ice thinning are those derived from our scaling factor approach. Each line in the figure corresponds to a different GIA model. As FDM differences seem to have a relatively small impact on the vertical velocities of GNSS stations, we have chosen to present here only the solution based on our scaling factor approach. But similar figures for each of the FDMs can be found in the supplementary materials. What we can see here is that no GIA model can satisfactorily explain all GNSS vertical velocities. Particularly in the southeast, where significant residuals remain, regardless of the model. Most GIA models

can significantly reduce discrepancies in the northwestern region (not considering NRD1 station whose velocity seems anomalous), except ICE-6G model, which tends to overestimate the vertical velocities there. However ICE-6G, as well as LM17.3 model, better explain vertical velocities at the northeastern region, while Caron and ANU models significantly underestimate the velocities there. Huy3 model predicts also anomalously large residuals for a few stations in southwestern region.

## 6 DISCUSSION AND CONCLUSIONS

One of the main results of our study is that firn corrections can lead to very different IMB solutions when interpreting surface elevation measurements in terms of ice mass changes. We investigated three FDM over two periods and show that there is considerable variability between models. This is actually consistent with previous studies on firn models comparison, that show that models are complex and that it is difficult to propagate uncertainties (Vandcrux et al. 2020). We develop a method, in which we estimate regional scaling factors that offers the possibility of optimally fitting IMB solutions inferred from SEC and estimates based on GRACE data. We show that an adequate division of Greenland surface based on
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Predicted gravity changes in equivalent water height in Greenland for the period 2011–2015, estimated from CryoSat-2 satellite altimetry data, peripheral glaciers information and assuming different FDM (IMAU, JPL, GSFC, scaling model from GRACE). The predicted gravity has been filtered as the GRACE data presented in Fig. 1; top: obtained trends depending on FDM; bottom: differences between GRACE estimate and top trends.</figcaption></figure>

<table><caption>**Table 4.** Total ice mass balance over Greenland estimated from SEC and gravity data, using the scaling approach and the equivalent mean density estimations shown in Table 3.</caption><thead><tr><td>Period</td><td>Scaling factor approach</td></tr></thead><tbody><tr><td>2003–2008</td><td>−217.6 ± 15.7 Gt yr<sup>−1</sup></td></tr><tr><td>2011–2015</td><td>−253.2 ± 18.8 Gt yr<sup>−1</sup></td></tr></tbody></table>

drainage basins (Rignot & Mouginot 2012; Zwally et al. 2012) provides a good fit to GRACE data with a scaling factor per region which has the dimension of a density in $g cm^{-3}$. We tested four different GIA models to correct GRACE (or even SEC) data, before empirically estimating a scaling factor for each region. We found that the GIA in this context had little effect on the scaling factor estimation. We therefore proposed an average scaling factor value per region between all GIA models, with an uncertainty of three times the standard deviation of scaling factors from all GIA models. This method makes it possible to recalculate and propose IMBs for Greenland that are consistent with both SECs and gravimetry, taking into account the uncertainties arising from gravimetry and GIA models.

For the period 2003–2008, we obtain total IMB estimates for Greenland of −215.0 Gt yr$^{-1}$ using FDM-IMAU, −261.9 Gt yr$^{-1}$ using FDM-GSFC, −221.5 Gt yr$^{-1}$ using FDM-JPL, and −217.6 ± 15.7 Gt yr$^{-1}$ using our scaling approach. Our total IMB estimate is smaller in absolute value than with FDM-GSFC and FDM-JPL but larger than with FDM-IMAU. It is relatively close to predictions from FDM-JPL and FDM-IMAU despite regional firn corrections that are quite different (see Table 3), in particular with FDM-IMAU. Our results suggest a significant misestimation of firn correction by FDM in different regions, notably in the southeastern regions. Nevertheless all IMB estimates remain consistent with other studies. Many studies have given a total IMB estimate in Greenland during this period ranging between −190 and −240 Gt yr$^{-1}$ (e.g. Van Den Broeke et al. 2009; Velicogna 2009; Sørensen et al. 2011; Bartletta et al. 2013), but others also propose values around −160 to −180 Gt yr$^{-1}$ (e.g. Wouters et al. 2008; Siemes et al. 2013; Otosaka et al. 2023). For the period 2011–2015, we obtain IMB estimates of −212.0 Gt yr$^{-1}$ for FDM-IMAU, −295.4 Gt yr$^{-1}$ for FDM-GSFC, −248.6 Gt yr$^{-1}$ for FDM-JPL and −253.2 ± 18.8 Gt yr$^{-1}$ using our scaling approach. All results are relatively consistent with those obtained in previous studies (e.g. McMillan et al. 2016; Kappelsberger et al. 2021; Khan et al. 2022). For example, Otosaka et al.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** GNSS vertical velocity observations (in colour) and predictions from elastic deformation due to the current ice loading (grey scale), during the period 2011–2015. The ice load has been inferred from local IMB determined from SEC data and firn corrections, using the three FDM (IMAU, GSFC, JPL) and our scaling approach. No GIA contribution has been added to the velocity prediction here. Colours associated with stations refer to those shown on Fig. 2.</figcaption></figure>

(2023) found an IMB of −213 ± 40 Gt yr⁻¹ over the period 2012–2016, while Mouginot et al. (2019) proposed an IMB of −286 ± 20 Gt yr⁻¹ over the period 2010–2018. Our values indicate an acceleration of the thinning rate, unlike some of the IMB estimates directly deduced from the ice SEC observations corrected by the three FDM. Here, we find that firn corrections can lead to differences in the total Greenland IMB up to almost 90 Gt yr⁻¹, depending on the period. It confirms that firn models are very sensitive to compaction model assumptions and to the choice of a climate model (Vandecrux et al. 2020). Consequently, the choice of a firn model can be a large source of uncertainty, so FDM should be used with caution (Stevens et al. 2020).

In addition to FDM, there are other possible sources of uncertainty in SEC that may lead to inaccurate estimates of IMB, such as the choice of the methods of interpolation (Sørensen et al. 2011), or radar signal penetration for CryoSat-2 (Nilsson et al. 2015; Sømønsen & Sørensen 2017).

Utilizing our 'scaling method' to optimize consistency between altimetry and gravimetry observations presents a promising strategy for mitigating dependence on FDM. In order to possibly bring further constraints on the IMB models, we also analysed GNSS vertical velocity data and compared them with predictions based on satellite altimetry observations, peripheral glaciers information and firn corrections. GNSS stations measure the displacement of bedrock due to elastic deformation induced by present-day ice mass loss, and viscoelastic deformation induced by past ice mass loss. We observe clear differences between GNSS observations and predictions that cannot be explained by FDM, particularly in the northern and

southeastern regions. The missing signal can be related to GIA or also to mass changes in the adjacent regions, especially those in the Canadian Arctic. Adhikari et al. (2021) show that this contribution can be negligible in most of Greenland except for the northwestern regions where it may cause a maximum uplift of ~1 mm yr⁻¹. We have taken into account GIA contributions using five different models in order to improve our predictions and try to explain the signal remaining at GNSS stations. Unfortunately, no GIA models investigated are able to explain all the missing signal. ICE-6G GIA model overestimates the signal in the northeastern regions, while Caron and ANU GIA model tends to overestimate it in the northwestern regions. Huy3 GIA model tends to increase discrepancies in the southwestern regions. LM17.3 GIA model seems to correlate slightly better than other GIA models over Greenland, but, in any case, we found significant residuals on the southeast coast of Greenland, for all GIA models. In particular four GNSS stations (KUAQ, MIK2, KSNB, HEL2) show very large uplift rates, with residuals larger than 7 mm yr⁻¹ (up to 15 mm yr⁻¹). Khan et al. (2016) have already reported a strong uplift in the southeastern part of Greenland, but mostly for the two stations KUAQ, MIK2. Different studies have highlighted the probable presence of lateral heterogeneities beneath the Greenland region (e.g. Darbyshire et al. 2018). Lateral variations in rheological properties, in particular mantle viscosity, may explain why none of the GIA models, all of which assume a 1-D viscosity profile, have succeeded in explaining the GNSS observations (Milne et al. 2018). Considering the particular case of the southeastern region, some studies propose to explain the discrepancies with observations at KUAQ and MIK2 stations by the presence of
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Residuals between GNSS vertical velocity observations and model predictions (in mm yr⁻¹) during the period 2011–2015, at each station location. The model prediction includes the effect of elastic deformation due to recent ice loading (based on IMB from the scaling factor approach) and different GIA models (ICE-6G, Caron, ANU, Huy3 and LM17.3 models, see Section 3.2 for more details). Colours associated with stations refer to those shown in Fig. 2.</figcaption></figure>

a low viscosity zone in the mantle due to the Iceland hotspot and its migration over time (Rogozhina *et al.* 2016; Steffen *et al.* 2018; Steinberger *et al.* 2019; Celli *et al.* 2021). The trace of the hotspot is still debated and we can see in our study that other close by stations, namely KSNB, HEL2, KBUG, PLPK and KULU, also show important, and unexplained uplift rates. Another explanation could be the signal coming from the ice mass loss since 1900, as proposed by Adhikari *et al.* (2021). This study predicts displacements due to

ice mass loss after the little ice age that can reach 5 mm yr⁻¹ in this region. This may explain some of the missing signal at the stations, but not all of it.

In conclusion, we show here that firn corrections can lead to particularly large errors in IMB estimates over Greenland, which is an issue that is probably underestimated today. No FDM appears able to correct for IMB based on SEC and fully explain gravity variation observations. We propose a method for overcoming this
issue, which optimally fits IMB solutions deduced from SEC and GRACE estimates. Using this method we find an IMB for the whole Greenland of −217.6 ± 15.7 Gt yr⁻¹ in 2003–2008 and −253.2 ± 18.8 Gt yr⁻¹ in 2011–2015. This method is particularly interesting as it appears to be not particularly sensitive to differences between GIA models in Greenland. GIA modelling remains however a major issue in Greenland, and the main limitations in explaining GNSS vertical velocity observations. While FDM models impact very significantly IMB estimates from SEC observations (possibly supplemented by gravimetry data), their impact on GNSS velocities appears relatively small. This impact is not negligible, but in view of the uncertainties arising from GIA modelling, it appears to be of secondary importance. Independently of GIA models, we find large discrepancies between GNSS observations and predictions in the southeastern part of Greenland, between Kangerlussaq glacier and Helheim glacier. This confirms previous studies suggesting the possible impact of either a low viscosity zone in the mantle, or the little ice age, on ground deformation in this region. Further analysis and modelling should be carried out to explain specifically the observations made in this region.

## ACKNOWLEDGMENTS

We thank Bert Vermeersen and two anonymous reviewers for their helpful reviews of our manuscript. We thank Louise Sandberg Sorensen and Sebastian Bjerregaard Simonsen for the ICESat and CryoSat-2 solutions. We deeply thank Etienne Berthier and Romain Hugonnet for their helps and advices on peripheral glacier data. We would like also to thank Nicole-Jeanne Schlegel and Eric Larour for discussions and helps on the JPL-FDM, Peter Kuipers Munneke et Brice Noël for discussions and helps on the IMAU-FDM. We also thank Benoit Lecavalier for providing us the Huy3 GIA model solution, and Kurt Lambeck for the ANU GIA solution. We investigated GNSS data downloaded from NGL website http://geodesy.unr.edu. Finally, we thank the Centre National d'Etudes Spatiales (CNES) for their financial support through the TOSCA committee.

## DATA AVAILABILITY

Data are available in a repository at https://doi.org/10.18715/IPGP.2024.lvgxi71n.

## REFERENCES

Adhikari, S., Milne, G., Caron, L., Khan, S., Kjeldsen, K., Nilsson, J., Larour, E. & Ivins, E., 2021. Decadal to centennial timescale mantle viscosity inferred from modern crustal uplift rates in Greenland, *Geophys. Res. Lett.*, **48**(19), e2021GL094040. doi:10.1029/2021GL094040

Altamimi, Z., Rebischung, P., Métivier, L. & Collilieux, X., 2016. ITRF2014: a new release of the international terrestrial reference frame modeling nonlinear station motions, *J. geophys. Res.: Solid Earth*, **121**(8), 6109–131.

Barletta, V.R., Sørensen, L.S. & Forsberg, R., 2013. Scatter of mass changes estimates at basin scale for Greenland and Antarctica, *The Cryosphere*, **7**(5), 1411–32.

Bertiger, W. et al., 2020. GipsyX/RTGx, a new tool set for space geodetic operations and research, *Adv. Space Res.*, **66**(3), 469–89.

Bevis, M. et al., 2012. Bedrock displacements in Greenland manifest ice mass variations, climate cycles and climate change, *Proc. Natl. Acad. Sci.*, **109**(30), 11944–11948.

Blewitt, G., Hammond, W.C. & Kreemer, C., 2018. Harnessing the GPS data explosion for interdisciplinary science, *Eos*, **99**(10.1029), 485, doi:10.1029/2018eo104623

Boehm, J., Werl, B. & Schuh, H., 2006. Troposphere mapping functions for GPS and very long baseline interferometry from European centre for medium-range weather forecasts operational analysis data, *J. geophys. Res.: Solid Earth*, **111**(B2), doi:10.1029/2005JB003629

Caron, L. & Ivins, E.R., 2020. A baseline antarctic gia correction for space gravimetry, *Earth planet. Sci. Lett.*, **531**, 115957. doi:10.1016/j.epsl.2019.115957

Caron, L., Ivins, E., Larour, E., Adhikari, S., Nilsson, J. & Blewitt, G., 2018. Gia model statistics for grace hydrology, cryosphere, and ocean science, *Geophys. Res. Lett.*, **45**(5), 2203–2212.

Cazenave, A. & Llovel, W., 2010. Contemporary sea level rise, *Annu. Rev. Marine Sci.*, **2**(1), 145–173.

Celli, N.L., Lebedev, S., Schaeffer, A.J. & Gaina, C., 2021. The tilted Iceland plume and its effect on the North Atlantic evolution and magmatism, *Earth planet. Sci. Lett.*, **569**, 117048. doi:10.1016/j.epsl.2021.117048

Cuffey, K.M., 2008. A matter of firn, *Science*, **320**(5883), 1596–1597.

Darbyshire, F.A., Dahl-Jensen, T., Larsen, T.B., Voss, P.H. & Joyal, G., 2018. Crust and uppermost-mantle structure of Greenland and the Northwest Atlantic from Rayleigh wave group velocity tomography, *Geophys. J. Int.*, **212**(3), 1546–1569.

Dobslaw, H. et al., 2017. A new high-resolution model of non-tidal atmosphere and ocean mass variability for de-aliasing of satellite gravity observations: AOD1B RL06, *Geophys. J. Int.*, **211**(1), 263–269.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference earth model, *Phys. Earth planet. Inter.*, **25**(4), 297–356.

Fettweis, X. et al., 2020. GrSMBMIP: intercomparison of the modelled 1980–2012 surface mass balance over the Greenland ice sheet, *The Cryosphere*, **14**(11), 3935–3958.

Gauer, L.-m., Chanard, K. & Fleitout, L., 2023. Data-driven gap filling and spatio-temporal filtering of the GRACE and GRACE-FO records, *J. geophys. Res.: Solid Earth*, e2022JB025561. doi:10.1029/2022JB025561

Gelaro, R. et al., 2017. The modern-era retrospective analysis for research and applications, version 2 (MERRA-2), *J. Clim.*, **30**(14), 5419–5454.

Ghil, M. et al., 2002. Advanced spectral methods for climatic time series, *Rev. Geophys.*, **40**(1), 3–1.

Gobron, K., Rebischung, P., Van Camp, M., Demoulin, A. & de Viron, O., 2021. Influence of aperiodic non-tidal atmospheric and oceanic loading deformations on the stochastic properties of global GNSS vertical land motion time series, *J. geophys. Res.: Solid Earth*, **126**(9), e2021JB022370. doi:10.1029/2021JB022370

Gobron, K., Rebischung, P., de Viron, O., Demoulin, A. & Van Camp, M., 2022. Impact of offsets on assessing the low-frequency stochastic properties of geodetic time series, *J. Geod.*, **96**(7), 1–16.

Harville, D.A., 1977. Maximum likelihood approaches to variance component estimation and to related problems, *J. Am. Stat. Assoc.*, **72**(358), 320–338.

Hersbach, H. et al., 2017. Complete ERA5 from 1940: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service (C3S) Data Store (CDS), doi:10.24381/cds.143582cf

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

Huss, M., 2013. Density assumptions for converting geodetic glacier volume change to mass change, *The Cryosphere*, **7**(3), 877–887.

Jiang, Y., Dixon, T.H. & Wdowinski, S., 2010. Accelerating uplift in the North Atlantic region as an indicator of ice loss, *Nat. Geosci.*, **3**(6), 404–407.

Kappelsberger, M.T., Strößenreuther, U., Scheinert, M., Horwath, M., Groh, A., Knöfel, C., Lunz, S. & Khan, S.A., 2021. Modeled and observed bedrock displacements in north-east Greenland using refined estimates of present-day ice-mass changes and densified GNSS measurements, *J. geophys. Res.: Earth Surface*, **126**(4), e2020JF005860. doi:10.1029/2020JF005860

Khan, S.A., Wahr, J., Bevis, M., Velicogna, I. & Kendrick, E., 2010. Spread of ice mass loss into northwest Greenland observed by GRACE and GPS: Northwest Greenland ice loss, *Geophys. Res. Lett.*, **37**(6), doi:10.1029/2010GL042460
Khan, S.A. et al., 2013. Recurring dynamically induced thinning during 1985 to 2010 on Upernavik Isstrøm, West Greenland, *J. geophys. Res.: Earth Surface*, **118**(1), 111–121.

Khan, S.A. et al., 2016. Geodetic measurements reveal similarities between post–last glacial maximum and present-day mass loss from the Greenland ice sheet, *Sci. Adv.* 2(9), e1600931. doi:10.1126/sciadv.1600931

Khan, S.A. et al., 2022. Greenland mass trends from airborne and satellite altimetry during 2011–2020, *J. geophys. Res.: Earth Surface*, **127**(4), e2021JF006505. doi:10.1029/2021JF006505

Kjeldsen, K.K. et al., 2015. Spatial and temporal distribution of mass loss from the Greenland ice sheet since AD 1900, *Nature*, **528**(7582), 396–400.

Kuipers Munneke, P. et al., 2015. Elevation change of the Greenland ice sheet due to surface mass balance and firn processes, 1960–2014, *The Cryosphere*, **9**(6), 2009–2025.

Kusche, J., 2007. Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models, *J. Geod.*, **81**(11), 733–749.

Kusche, J., Schmidt, R., Petrovic, S. & Rietbroek, R., 2009. Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hydrological model, *J. Geod.*, **83**, 903–913.

Lambeck, K., Rouby, H., Purcell, A., Sun, Y. & Sambridge, M., 2014. Sea level and global ice volumes from the last glacial maximum to the holocene, *Proc. Natl. Acad. Sci.*, **111**(43), 15 296–15 303.

Larour, E., Seroussi, H., Morlighem, M. & Rignot, E., 2012. Continental scale, high order, high spatial resolution, ice sheet modeling using the ice sheet system model (ISSM), *J. geophys. Res.: Earth Surface*, **117**(F1), F01022. doi:10.1029/2011JF002140

Laske, G., Masters, G., Ma, Z. & Pasyanos, M., 2013. Update on CRUST1.0
- A 1-degree Global Model of Earth's Crust, *Geophys. Res. Abs.*, **15**, Abstract EGU2013-2658.

Lecavalier, B.S. et al., 2014. A model of Greenland ice sheet deglaciation constrained by observations of relative sea level and ice extent, *Quat. Sci. Rev.*, **102**, 54–84.

Ligtenberg, S., Helsen, M. & Van den Broeke, M., 2011. An improved semi-
empirical model for the densification of Antarctic firn, *The Cryosphere*,
*5*(4), 809–819.

Ligtenberg, S.R., Kuipers Munneke, P., Noël, B.P. & Van Den Broeke, M.R., 2018. Brief communication: improved simulation of the present-day Greenland firn layer (1960–2016), *The Cryosphere*, **12**(5), 1643–1649.

Mao, A., Harrison, C.G. & Dixon, T.H., 1999. Noise in GPS coordinate time series, *J. geophys. Res.: Solid Earth*, **104**(B2), 2797–2816.

Martens, H. R., Rivera, L. & Simons, M., 2019. LoadDef: A python-based toolkit to model elastic deformation caused by surface mass loading on spherically symmetric bodies, *Earth Space Sci.*, **6**(2), 311–323.

McMillan, M. et al., 2016. A high-resolution record of Greenland mass balance, *Geophys. Res. Lett.*, **43**(13), 7002–7010.

Medley, B., Neumann, T.A., Zwally, H.J., Smith, B.E. & Setvens, Max, C., 2022. Simulations of firn processes over the Greenland and Antarctic ice sheets: 1980-2021, *The Cryosphere*, **16**, 3971–4011. doi:10.5194/tc-16-3971-2022

Métivier, L., Greff-Leffitz, M. & Altamimi, Z., 2010. On secular geocenter motion: the impact of climate changes, *Earth planet. Sci. Lett.*, **296**(3–4), 360–366.

Milne, G.A., Latychev, K., Schaeffer, A., Crowley, J.W., Lecavalier, B.S. & Audette, A., 2018. The influence of lateral earth structure on glacial isostatic adjustment in Greenland, *Geophys. J. Int.*, **214**(2), 1252–1266.

Mouginot, J., Rignot, E., Bjørk, A.A., Van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B. & Wood, M., 2019. Forty-six years of Greenland ice sheet mass balance from 1972 to 2018, *Proc. Natl Acad. Sci.*, **116**(19), 9239–9244.

Nilsson, J. *et al.*, 2015. Greenland 2012 melt event effects on cryosat-2 radar altimetry, *Geophys. Res. Lett.*, **42**(10), 3919–3926.

Noël, B., van de Berg, W.J., Machguth, H., Lhermitte, S., Howat, I., Fettweis, X. & Van Den Broeke, M.R., 2016. A daily, 1 km resolution data set of downscaled Greenland ice sheet surface mass balance (1958–2015), *The Cryosphere*, **10**(5), 2361–2377.

Noël, B. *et al.*, 2018. Modelling the climate and surface mass balance of polar ice sheets using RACMO2-part 1: Greenland (1958–2016), *The Cryosphere*, **12**(3), 811–831.

Otosaka, I.N. *et al.*, 2023. Mass balance of the Greenland and arctic ice sheets from 1992 to 2020, *Earth Syst. Sci. Data Discuss.*, **15**, 1597–1616.

Patterson, H.D. & Thompson, R., 1971. Recovery of inter-block information when block sizes are unequal, *Biometrika*, **58**(3), 545–554.

Peltier, W.R., Argus, D. & Drummond, R., 2015. Space geodesy constrains ice age terminal deglaciation: the global ICE-6G_C (VM5A) model, *J. geophys. Res.: Solid Earth*, **120**(1), 450–487.

Prevost, P., Chanard, K., Fleitout, L., Calais, E., Walwer, D., van Dam, T. & Ghil, M., 2019. Data-adaptive spatio-temporal filtering of GRACE data, *Geophys. J. Int.*, **219**(3), 2034–2055.

Rebischung, P., Altamimi, Z., Ray, J. & Garayt, B., 2016. The IGS contribu-
tion to ITRF2014, *J. Geod.*, **90**(7), 611–630.

Rignot, E. & Mouginot, J., 2012. Ice flow in Greenland for the international polar year 2008–2009, *Geophys. Res. Lett.*, **39**(11).

Rogozhina, I. et al., 2016. Melting at the base of the Greenland ice sheet explained by Iceland hotspot history, *Nat. Geosci.*, **9**(5), 366–369.

Santamaria-Gómez, A., Bouin, M.-N., Collilieux, X. & Wöppermann, G., 2011. Correlated errors in GPS position time series: implications for velocity estimates, *J. geophys. Res.: Solid Earth*, **116**(B1). doi:10.1029/2010JB007701

Save, H., Bettadpur, S. & Tapley, B.D., 2016. High-resolution CSR GRACE RL05 mascons, *J. geophys. Res.: Solid Earth*, **121**(10), 7547–7569.

Shepherd, A. et al., 2012. A reconciled estimate of ice-sheet mass balance, *Science*, **338**(6111), 1183–1189.

Shepherd, A. et al., 2020. Mass balance of the Greenland ice sheet from 1992 to 2018, *Nature*, **579**(7798), 233–239.

Siemes, C., Ditmar, P., Riva, R., Slobbe, D., Liu, X. & Farahani, H.H., 2013. Estimation of mass change trends in the earth's system on the basis of grace satellite data, with application to Greenland, *J. Geod.*, **87**(1), 69–87.

Simonsen, S.B. & Sørensen, L.S., 2017. Implications of changing scattering properties on Greenland ice sheet volume change from cryosat-2 altimetry, *Remote Sens. Environ.*, **190**, 207–216.

Sørensen, L.S., Simonsen, S.B., Nielsen, K., Lucas-Picher, P., Spada, G., Adalgeirdsottir, G., Forsberg, R. & Hvidberg, C.S., 2011. Mass balance of the Greenland ice sheet (2003–2008) from icesat data—the impact of interpolation, sampling and firn density, *The Cryosphere*, **5**(1), 173–186.

Sørensen, L.S., Simonsen, S.B., Forsberg, R., Khvorostovsky, K., Meister, R. & Engdahl, M.E., 2018. 25 years of elevation changes of the Greenland ice sheet from ERS, envisat, and cryosat-2 radar altimetry, *Earth planet. Sci. Lett.*, **495**, 234–241.

Spada, G., Ruggieri, G., Sørensen, L.S., Nielsen, K., Melini, D. & Colleoni, F., 2012. Greenland uplift and regional sea level changes from ICEsAt observations and GIA modelling: Greenland uplift and sea level change, *Geophys. J. Int.*, **189**(3), 1457–1474.

Steffen, H., Li, T., Wu, P., Gowan, E.J., Ivins, E., Lecavalier, B., Tarasov, L. & Whitehouse, P.L., 2021. Lm17. 3-a global land vertical land motion model of glacial isostatic adjustment [data set], PANGEA. doi:

Steffen, R., Audet, P. & Lund, B., 2018. Weakened lithosphere beneath Greenland inferred from effective elastic thickness: a hot spot effect?, *Geophys. Res. Lett.*, **45**(10), 4733–4742.

Steinberger, B., Bredow, E., Lebedev, S., Schaeffer, A. & Torsvik, T.H., 2019. Widespread volcanism in the Greenland–North Atlantic region explained by the Iceland plume, *Nat. Geosci.*, **12**(1), 61–68.

Stevens, C.M., Verjans, V., Lundin, J., Kahle, E.C., Horlings, A.N., Horlings, B.I. & Waddington, E.D., 2020. The community firn model (CFM) v1.0, *Geosci. Model Develop.*, **13**(9), 4355–4377.

Sutterley, T.C., Velicogna, I., Csatho, B., van den Broeke, M., Rezvan-
Behbahani, S. & Babonis, G., 2014. Evaluating Greenland glacial iso-
static adjustment corrections using grace, altimetry and surface mass
balance data, *Environ. Res. Lett.*, **9**(1), 014004. doi:10.1088/1748-
9326/9/1/014004

Tapley, B.D., Bettadpur, S., Ries, J.P., Thompson, P.F. & Watkins, M.M., 2004. Grace measurements of mass variability in the earth system, *Science*, **305**(5683), 503–505.
Van Den Broeke, M. et al., 2009. Partitioning recent Greenland mass loss, *Science*, **326**(5955), 984–986.

Vandecrux, B. et al., 2020. The firn meltwater retention model intercomparison project (RetMIP): evaluation of nine firn models at four weather station sites on the Greenland ice sheet, *The Cryosphere*, **14**(11), 3785–3810.

Velicogna, I., 2009. Increasing rates of ice mass loss from the Greenland and Antarctic ice sheets revealed by grace, *Geophys. Res. Lett.*, **36**(19). doi:10.1029/2009GL040222

Velicogna, I. et al., 2020. Continuity of ice sheet mass loss in Greenland and Antarctica from the grace and grace follow-on missions, *Geophys. Res. Lett.*, **47**(8), e2020GL087291.

Wahr, J. & Zhong, S., 2013. Computations of the viscoelastic response of a 3-D compressible earth to surface loading: an application to glacial isostatic adjustment in Antarctica and Canada, *Geophys. J. Int.*, **192**(2), 557–572.

Wake, L.M., Lecavalier, B.S. & Bevis, M., 2016. Glacial isostatic adjustment (GIA) in Greenland: a review, *Curr. Clim. Change Rep.*, **2**(3), 101–111.

Wouters, B., Chambers, D. & Schrama, E., 2008. Grace observes small-scale mass loss in Greenland, *Geophys. Res. Lett.*, **35**(20). doi:10.1029/2008GL034816

Zhang, B. et al., 2019. Geodetic and model data reveal different spatio-temporal patterns of transient mass changes over Greenland from 2007 to 2017, *Earth planet. Sci. Lett.*, **515**, 154–163.

Zumberge, J., Heflin, M., Jefferson, D., Watkins, M. & Webb, F., 1997. Precise point positioning for the efficient and robust analysis of GPS data from large networks, *J. geophys. Res.: Solid Earth*, **102**(B3), 5005–5017.

Zwally, H. J., Giovinetto, M. B., Beckley, M. A. & Saba, J. L., 2012. *Antarctic and Greenland Drainage Systems*, GSFC Cryospheric Sciences Laboratory. Available at: https://earth.gsfc.nasa.gov/cryo/data/polar-altimetry/antarctic-and-Greenland-drainage-systems