# Secular gravity variations at Svalbard (Norway) from ground observations and GRACE satellite data

Anthony Mémin, Yves Rogister, Jacques Hinderer, Ove Omang, Bernard Luck

► To cite this version:

Anthony Mémin, Yves Rogister, Jacques Hinderer, Ove Omang, Bernard Luck. Secular gravity variations at Svalbard (Norway) from ground observations and GRACE satellite data. Geophysical Journal International, 2011, 184, pp.1119-1130. 10.1111/j.1365-246X.2010.04922.x. insu-00675762

HAL Id: insu-00675762

https://insu.hal.science/insu-00675762

Submitted on 23 Aug 2021

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Secular gravity variation at Svalbard (Norway) from ground observations and GRACE satellite data

A. Mémin,<fnref n="1" /> Y. Rogister,<fnref n="1" /> J. Hinderer,<fnref n="1" /> O. C. Omang<fnref n="2" /> and B. Luck<fnref n="1" />

Accepted 2010 December 10. Received 2010 December 9; in original form 2010 March 5

## SUMMARY

The Svalbard archipelago, Norway, is affected by both the present-day ice melting (PDIM) and Glacial Isostatic Adjustment (GIA) subsequent to the Last Pleistocene deglaciation. The induced deformation of the Earth is observed by using different techniques. At the Geodetic Observatory in Ny-Ålesund, precise positioning measurements have been collected since 1991, a superconducting gravimeter (SG) has been installed in 1999, and six campaigns of absolute gravity (AG) measurements were performed between 1998 and 2007. Moreover, the Gravity Recovery and Climate Experiment (GRACE) satellite mission provides the time variation of the Earth gravity field since 2002. The goal of this paper is to estimate the present rate of ice melting by combining geodetic observations of the gravity variation and uplift rate with geophysical modelling of both the GIA and Earth's response to the PDIM. We estimate the secular gravity variation by superimposing the SG series with the six AG measurements. We collect published estimates of the vertical velocity based on GPS and VLBI data. We analyse the GRACE solutions provided by three groups (CSR, GFZ, GRGS). The crux of the problem lies in the separation of the contributions from the GIA and PDIM to the Earth's deformation. To account for the GIA, we compute the response of viscoelastic Earth models having different radial structures of mantle viscosity to the deglaciation histories included in the models ICE-3G or ICE-5G. To account for the effect of PDIM, we compute the deformation of an elastic Earth model for six models of ice-melting extension and rates. Errors in the gravity variation and vertical velocity are estimated by taking into account the measurement uncertainties and the variability of the GRACE solutions and GIA and PDIM models. The ground observations agree with models that involve a current ice loss of 25 km³ water equivalent yr⁻<fnref n="1" /> over Svalbard, whereas the space observations give a value in the interval [5, 18] km³ water equivalent yr⁻<fnref n="1" />. A better modelling of the PDIM, which would include the precise topography of the glaciers and altitude-dependency of ice melting, is necessary to decrease the discrepancy between the two estimates.

**Key words:** Satellite geodesy; Time variable gravity; Global change from geodesy; Arctic region.

## 1 INTRODUCTION

The Svalbard archipelago, Norway, is affected by the post-glacial rebound, or Glacial Isostatic Adjustment (GIA), subsequent to the last Pleistocene deglaciation that started 21 000 years ago and ended 10 000 years ago. Moreover, most of the ice sheets in Svalbard (36 000 km<fnref n="2" />) are presently thinning (Kohler et al. 2007; Dowdeswell et al. 2008; Kääb 2008; Nuth et al. 2010). Both the last deglaciation and present-day ice melting (PDIM) induce deformation of the Earth. These two effects make Svalbard a very interesting zone to study to better understand the geodetic consequences of ice melting at different timescales.

For many years (from 12 to 18 years), positioning observations with different systems (GPS, VLBI, DORIS) have been collected

at the Geodetic Observatory in Ny-Ålesund, Svalbard (e.g. Sato et al. 2006; Kierulf et al. 2009a). Kierulf et al. (2009b) study the ground velocity and find that a predicted melting rate of 37 cm water equivalent yr⁻<fnref n="1" /> (we yr⁻<fnref n="1" />) explains up to 60 per cent of their observed uplift, which is 8.2 ± 0.9 mm yr⁻<fnref n="1" /> in ITRF2005. They attribute the remaining 40 per cent to imperfect modelling and errors in determining the secular rate from observations. Since 1999, a superconducting gravimeter (SG) has been continuously recording the gravity changes and, from 1998 to 2007, six campaigns of absolute gravity (AG) measurements were performed at the SG station, which also contributes to the Global Geodynamics Project (Crossley et al. 1999). Sato et al. (2006), who were interested in the secular geodetic variations, concluded from the first four AG measurements and observed uplift rate that an ice melting rate of

<footnote n="1">¹IPGS-EOST, CNRS/UdS, UMR 7516, 5 rue René Descartes, 67084 Strasbourg Cedex, France. E-mail: anthony.memin@unistra.fr</footnote>

<footnote n="2">²Norwegian Mapping Authority, 3507 Honefoss, Norway</footnote>
75 cm we yr⁻¹ would explain the ground vertical velocity but only half of the gravity rate.

Since its launch in 2002, the Gravity Recovery and Climate Experiment (GRACE) satellite mission has provided a global mapping of the time-varying gravity field of the Earth. These satellite gravimetry data have been used for the estimation of the ice-mass balance over ice-covered lands of various extension, such as Greenland, Antarctica (e.g. Chen et al. 2006c; Luthcke et al. 2006; Velicogna & Wahr 2006a; Velicogna & Wahr 2006b), Alaska and Patagonia (Tamisiea et al. 2005; Chen et al. 2006a; Chen et al. 2007; Luthcke et al. 2008). Their areas are respectively 12 × 10⁶ km², 1.7 × 10⁶ km², 90 000 km² and 17 000 km². Rodell & Famiglietti (1999) compared the modelled variation of water storage to the expected uncertainty associated to the estimate of that variation from future GRACE observations. They concluded that 200 000 km² is the smallest area for which it will be possible to detect seasonal and annual changes. Rodell & Famiglietti (2001) later showed that the seasonal and annual changes of the smaller Illinois basin (United States, 145 800 km²) could be partly—up to 50 per cent—detected by the GRACE mission. However, in spite of the size of the Alaskan and Patagonian ice sheets, which are smaller than the Illinois basin, the GRACE data have revealed a large signal over these regions that has been associated with ice-mass loss of, respectively, 97 ± 13 and 28 km³ we yr⁻¹. These values agree with the estimates of ice-mass change derived from glaciology studies (Tamisiea et al. 2005; Chen et al. 2007; Luthcke et al. 2008). The size of the ice-covered area of the Svalbard archipelago being intermediate between those of Alaska and Patagonia, a sufficiently large amount of ice mass variation over Svalbard can be observed with the GRACE satellites. According to the latest elevation changes provided by ICESat laser altimetry (Moholdt et al. 2010b), the average water equivalent balance of Svalbard is −12 ± 4 cm we yr⁻¹. The ocean that surrounds Svalbard slightly contributes to the gravimetric signal, from 0.8 ± 0.8 to 1.9 ± 0.1 mm we yr⁻¹ (see Quinn & Ponte 2010, for a review). Because the GRACE signal over Svalbard is approximately 60 times larger, most of it can be attributed to land-mass changes.

We aim at providing new estimates of the PDIM over Svalbard. To do so, we have at our disposal updated gravity and precise positioning measurements at the Ny-Ålesund Geodetic Observatory, as well as time-series of the Earth gravity field provided by the GRACE satellite. To separate the GIA signal from the signal coming from PDIM requires both phenomena to be modelled.

The outline of the paper goes as follows. In Section 2, we analyse ground and satellite observations to estimate the gravity variation and ground vertical velocity over Svalbard. In Section 2.1, we first derive the secular variation of gravity using six AG observations. We also calibrate the SG by superimposing its nine-year series and the AG-derived trend. In Section 2.2, we analyse up to 6 years (2003–2009) of GRACE solutions. By considering solutions provided by different centres [Center for Space Research (CSR), GeoforschungsZentrum (GFZ) and Groupe de Recherche en Géodésie Spatiale (GRGS)] and different lengths for the time-series (5 or 6 yr), we provide uncertainties to the gravity variation observed by the GRACE satellites. In Section 2.3, we collect uplift rates found in the literature, they are obtained from GPS and VLBI observations. Next, in Section 3, we model the Earth response to ice melting. We successively consider the ground and space gravity variations in Sections 3.1 and 3.2, respectively. We use the SELEN software developed by Spada & Stocchi (2007) to compute the response of a viscoelastic earth model to the past deglaciation history included in the model ICE-3G (Tushingham & Peltier 1991). To take the PDIM into account, we compute the ground velocity and gravity rate of an

elastic earth model for two models of ice-thinning extension successively combined with three ice-thinning rates. In Section 4, we compare the modelled gravity variation and vertical velocity resulting from both the PDIM and GIA to the observed gravity variation and uplift rate. We consider ground observations in Section 4.1 and GRACE observations in Section 4.2. In Section 5, we discuss various estimates of ice-volume changes. Finally, we summarize our work in Section 6.

## 2 OBSERVATIONS OF GRAVITY VARIATION AND VERTICAL MOTION

### 2.1 Ground gravity measurements

A total of six observations with FG5 absolute gravimeters have been made between 1998 and 2007 at the Ny-Ålesund SG station. Three European institutes participated in the gravity observations: the Bundesamt fuer Kartographie and Geodaesie (BKG, Frankfurt, Germany), Ecole et Observatoire des Sciences de la Terre (EOST, Strasbourg, France) and European Center for Geodynamics and Seismology (ECGS, Walferdange, G.-D. Luxembourg). Four of the six observations (1998–2002) were reported in the paper by Sato et al. (2006). For each measurement, the gravity value at the height of the AG dropping chamber was transferred to the ground by using a constant vertical gradient of −0.3594 µGal mm⁻¹. Moreover, corrections were applied to the raw data, owing to three geophysical phenomena:

(i) the observed polar motion, provided by the International Earth Rotation Service (IERS),

(ii) the observed tides, including the effect of ocean tide loading, and

(iii) the observed atmospheric pressure which is responsible for a gravity variation of −0.42 µGal hPa⁻¹.

We apply the same corrections to the two additional AG observations of 2004 and 2007. All the AG data are shown in Fig. 1 with their 3σ formal errors.

Sato et al. (2006) showed that using the observed parameters for the tide and air pressure corrections, instead of the default ones programmed in the processing package g provided by the manufacturer Micro-g, clearly improves the standard error of the estimated AG value. However, using different FG5 gravimeters can lead to gravity values with an rms error of up to 2 µGal (Robertsson et al. 2001). This was confirmed by comparing 13 absolute gravimeters (including 11 FG5) in 2003 (Francis et al. 2005): a standard deviation of 1.9 µGal was obtained. Consequently, Sato et al. (2006) suggested to add 2 µGal in quadrature to all the formal errors and fit the data with the new error as a weight. They obtained a linear trend of −2.5 ± 0.9 µGal yr⁻¹ (Fig. 1). By applying the same extra correction to the formal errors of the 6 AG observations, we compute the trend with the 3σ formal errors as a weight in the linear least-square fitting of the data and obtained −1.02 ± 0.48 µGal yr⁻¹, which is about 41 per cent of the value found by Sato et al. (2006) who used the first four observations only.

We also analyse 9 yr (2000 January 1 to 2008 December 31) of SG data at Ny-Ålesund. We apply the three geophysical corrections listed earlier. We use a degree-3 polynomial filtering with a 48 h window to remove the residuals of the tidal signal from the time-series. We correct for five offsets larger than 10 µGal mostly due to the refilling of the SG dewar with liquid helium. Finally, we also correct the SG data for a linear drift of 1.81 µGal yr⁻¹ of
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Absolute Gravity (AG) measurements at Ny-Ålesund (after subtraction of a mean value of 983 017 053.76 µGal). The slopes of the dashed (Sato et al. 2006) and solid (this study) lines are respectively -2.5 ± 0.9 and -1.02 ± 0.48 µGal yr⁻¹. The observations were made by BKG (Bundesamt fuer Kartographie and Geodaesie, Frankfurt, Germany), EOST (Ecole et Observatoire des Sciences de la Terre, Strasbourg, France) and ECGS (European Center for Geodynamics and Seismology, Walferdange, G.-D. Luxembourg).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Superposition of relative gravity observations with a Superconducting Gravimeter (black curve) and AG measurements (red cross) at Ny-Ålesund. The red dashed line is the AG linear trend shown as a solid black line in Fig. 1.</figcaption></figure>

instrumental origin by adjusting the SG and AG data. The corrected SG time-series is shown in Fig. 2.

The AG measurements have been made in early 1998 September, late 2000 July, late 2001 and 2002 September, and early 2007 June. At those different times of the year, both the hydrologic and periodic signals have different amplitudes. To make sure that the secular trend deduced from the AG measurements is linear and not affected by the date of the measurements, we have to compare the AG values with the SG data. By correcting the SG drift to fit the secular trend obtained from the AG measurements, we can compare both data sets. If the linear trend deduced from the AG observations is correct, the AG measurements must coincide with the SG data. Otherwise, the AG data are off the SG series. The superposition of

the AG and the SG data shows a very good agreement, confirming the linear trend that we have estimated for the gravity rate in Ny-Ålesund.

## 2.2 Space gravity measurements

The releases 4 and 2 of the GRACE solutions provided respectively by the CSR/GFZ (Bettadpur 2007; Flechtner 2007) and GRGS (Bruinsma et al. 2010) centres are available as the fully normalized spherical harmonic coefficients or Stokes coefficients $C_n^m(t)$ and $S_n^m(t)$ of the monthly (CSR and GFZ) and 10-day (GRGS) gravity potential $V^{\text{sat}}(\theta, \lambda, a, t)$, which is at the surface of the Earth,

$$V^{\text{sat}}(\theta, \lambda, a, t) = \frac{GM}{a} \sum_{n=2}^{N_{\text{max}}} \sum_{m=0}^{n} [C_n^m(t) Y_n^{m,c}(\theta, \lambda) + S_n^m(t) Y_n^{m,s}(\theta, \lambda)], \quad (1)$$

where $\theta$ is the colatitude, $\lambda$ is the longitude, $G$ is the Newtonian constant of gravitation, $M$ is the mass of the Earth, $GM=3.986 \times 10^{14} \text{ m}^3\text{s}^{-2}$, $a=6378$ km is the mean equatorial radius, and $Y_n^{m,c}(\theta, \lambda)$ and $Y_n^{m,s}(\theta, \lambda)$ are the real cosine and sine fully normalized spherical harmonics of degree $n$ and order $m$, respectively. While other studies consider gravity anomalies (e.g. Matsuo & Heki 2010), we consider the gravity disturbance at $a$ given by

$$g^{\text{sat}}(\theta, \lambda, a, t) = \frac{GM}{a^2} \sum_{n=2}^{N_{\text{max}}} \sum_{m=0}^{n} (n+1) \times [C_n^m(t) Y_n^{m,c}(\theta, \lambda) + S_n^m(t) Y_n^{m,s}(\theta, \lambda)]. \quad (2)$$

We estimate the gravity disturbance rates, $\Delta g^{\text{sat}}$, from the variations of $g^{\text{sat}}$ with respect to the appropriate static fields taken over a 6-year period (2003 January–2009 January). We consider 72 monthly solutions from CSR, 70 monthly solutions from GFZ and 210 10-day solutions from GRGS, all the solutions being computed up to the harmonic degree $N_{\text{max}} = 50$. This is the maximum degree
for which the GRGS provides harmonic coefficients and it corresponds to a spatial resolution of ~400 km. The GRGS solutions are stabilized during their generation process (http://bgi.grgs.fr) by gradually constraining the coefficients of degree 2 through degree 50 to the coefficients of the static field (Bruinsma et al. 2010). Besides, we filter the monthly solutions of the CSR and GFZ centres with a low-pass filter in the spectral domain. The low-pass filter is a cosine taper decreasing from 1 at n = 30 to 0 at n = 50 (e.g. de Linage et al. 2009). Finally, gravity variation Δg<sup>sat</sup> is estimated with a least-square fitting of

$$g(t) = \sum_{i=1}^{2} [a_i \cos(i\omega t) + b_i \sin(i\omega t)] + \Delta g^{\text{sat}} t + c, \quad (3)$$

where ω = 2π/T and T = 1 yr. The a<sub>i</sub>'s and b<sub>i</sub>'s give the magnitude of the annual (i = 1) and semi-annual (i = 2) cycles, and c is the static part of g ⋅ Δg<sup>sat</sup> is shown in Fig. 3.

The three solutions show similar geophysical signals, for they contain large patterns over Greenland and Fennoscandia. Over Greenland, the time-decreasing mass signal, implying negative gravity rates, has already been reported in previous studies (e.g. Howat et al. 2007; Barletta et al. 2008; Slobbe et al. 2009). Over Fennoscandia, the post-glacial rebound due to the last deglaciation induces the time-increasing mass signal, implying positive gravity rates, which can also be found in previous studies (e.g. Steffen et al. 2008; Steffen et al. 2009). A signal over Svalbard with a magnitude smaller than −0.5 µGal yr<sup>−1</sup> can be seen in the CSR and GRGS solutions while it is smaller, in absolute value, for the GFZ solutions. To emphasize the geophysical signals that should be common to the three solutions and decrease the magnitude of other uncorrelated signals, we compute the mean of the three solutions (Fig. 4, top panel). The amplitude of the signal over Svalbard remains about −0.5 µGal yr<sup>−1</sup> (Table 1).

We now investigate the influence of the length and variability of the solutions on the gravity trend. We compute the spatial mean of the GRACE signals over Svalbard for each solution and the averages of the three solutions for the two time intervals 2003 January and 2008 January and 2003 January and 2009 January. North-south stripes affect the GRACE solutions and, to a lesser extent, the mean of the gravity variation. To see the influence of the stripes, we make the computation before and after applying the destriping filter of Swenson & Wahr (2006), recently detailed by Duan et al. (2009). We compute the standard deviation of the gravity trend for the three solutions and for the two time intervals.

Regarding the influence of the length of the time-series, Table 1 shows that the gravity trends for the three solutions, either undestriped or destriped, are smaller for the 2003–2008 time interval than for the 2003–2009 time interval. The differences range from 0.15 to 0.37 µGal yr<sup>−1</sup>. The averages of the undestriped and destriped solutions are respectively 0.28 and 0.20 µGal yr<sup>−1</sup> smaller between 2003 and 2008 than between 2003 and 2009. Interannual geophysical phenomena are probably responsible for the difference between the estimated gravity trends.

Table 1 also shows that the solutions provided by the three centres are very different. However, the destriping process decreases the discrepancy and gives much smaller standard deviations for the destriped solutions.

The observation of the gravity field should be as long as possible to determine accurately its long-term non-periodic variation. Consequently, although the standard deviations for both the undestriped and destriped are smaller for the 2003–2008 solutions, we will base our study on the 2003–2009 destriped solutions.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Gravity trends Δg<sup>sat</sup>, in µGal yr<sup>−1</sup>, between 2003 January and 2009 January deduced from the GRACE data. From top to bottom: CSR-RL04, GFZ-RL04 and GRGS-RL02 solutions. The CSR and GFZ solutions are filtered with a cosine taper in the spectral domain for the harmonic degrees between 30 and 50.</figcaption></figure>

## 2.3 Ground velocity measurements

In this study we use the ground velocities provided in the literature. Noticing that a scale factor existed between the VLBI (1994–2004) and GPS data, Sato et al. (2006) corrected the two GPS time-series (1991–2004 and 1998–2004) and computed the average of all the observations in the terrestrial reference frame ITRF2000. They obtained an uplift rate of 5.2 ± 0.6 mm yr<sup>−1</sup>. Rülke et al. (2008) provide a GPS-only Terrestrial Reference Frame (TRF), the Potsdam-Dresden-Reprocessing TRF. In this TRF, the vertical velocities of the two stations NYAL (1994–2005) and NYA1 (1998–2005) at Ny-Ålesund are respectively 7.3 and 7.1 mm yr<sup>−1</sup> which lead to a mean uplift rate of 7.2 mm yr<sup>−1</sup> and a standard deviation of 0.1 mm yr<sup>−1</sup>. In ITRF 2000 this rate approximately reduces to 5.66 ± 0.3 mm yr<sup>−1</sup>, which is close to the rate obtained by Sato
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top panel: mean of the gravity rates deduced from the CSR-RL04, GFZ-RL04 and GRGS-RL02 solutions of the GRACE data. The CSR and GFZ solutions have been previously filtered with a cosine taper in the spectral domain between degrees 30 and 50. Middle panel: Earth gravity response $\Delta g^{GIA}$ to the deglaciation history ICE-3G with the viscosity profile VP used by Sato *et al.* (2006). Bottom panel: mean of the gravity rates deduced from the CSR, GFZ and GRGS solutions that have been destriped and corrected for the GIA contribution. The unit is $\mu\text{Gal yr}^{-1}$.</figcaption></figure>

*et al.* (2006). Recently, using VLBI, GPS and DORIS data, Kierulf *et al.* (2009a) derived a mean uplift rate of $8.2 \pm 0.9 \text{ mm yr}^{-1}$ in ITRF2005. Kierulf *et al.* (2009a), Kierulf *et al.* (2009b) show that the geodetic observations suffer from the processing strategies. Indeed, they obtain vertical velocities of $10.8$ and $7.6 \text{ mm yr}^{-1}$, respectively, using GIPSY-PPP and GAMIT-DD solutions for the NYA1 GPS station between 1997 and 2007. Similar results are obtained for the second GPS station NYAL between 1991 and 2007. In ITRF2000, their mean uplift reduces to $6.43 \pm 1.19 \text{ mm yr}^{-1}$.

Whereas the dispersion of the published horizontal velocities is less than $1 \text{ mm yr}^{-1}$, the vertical velocity takes on different values depending on the TRF, processing strategy and measurement method.

<table><caption>**Table 1.** Mean gravity trend over Svalbard for three GRACE solutions for 5- and 6-year time intervals (from 2003 to either 2008 or 2009). The asterisk indicates that the mean is computed after applying the destriping filter of Swenson & Wahr (2006). Two bottom lines: average and standard deviation $\sigma$ of the three GRACE solutions.</caption><thead><tr><th rowspan="2">Time interval</th><th>2003 January to 2008</th><th>2003 January to 2009</th></tr><tr><th>January</th><th>January</th></tr></thead><tbody><tr><td>Solution</td><td colspan="2">Gravity trend ($\mu\text{Gal yr}^{-1}$)</td></tr><tr><td>GRGS</td><td>-0.99</td><td>-0.80</td></tr><tr><td>CSR</td><td>-0.85</td><td>-0.90*</td></tr><tr><td>GFZ</td><td>-0.55</td><td>-0.95*</td></tr><tr><td>Average</td><td>-0.80</td><td>-0.90*</td></tr><tr><td>$\sigma$</td><td>0.22</td><td>0.05*</td></tr></tbody></table>

According to Tregoning & Lambeck (2010), the origin of ITRF2000 may be closer to the centre of mass of Earth's system than that of ITRF2005. This would lead to a better consistency between GPS uplift rates and GIA model predictions (Tregoning, personal communication, 2010). That was previously mentioned by Argus (2007). Using the errors as a weight, we compute the mean of the three uplift rates in ITRF2000. We obtain $5.64 \pm 1.57 \text{ mm yr}^{-1}$ with a $3\sigma$ error.

## 3 MODELLING OF GIA AND ELASTIC DEFORMATION DUE TO PDIM

### 3.1 Modelling of ground gravity variation

Sato *et al.* (2006) computed the vertical displacement and gravity variation associated to the GIA at Ny-Ålesund up to the harmonic degree 180. They considered a Maxwell Earth with the viscosity profile VP listed in Table 2 and the deglaciation history contained in the model ICE-3G.

The computation by Sato *et al.* (2006) is more accurate than the one based on the publicly available SELEN software (Spada & Stocchi 2007), which will be used in Section 3.2, because they consider a time-dependent ocean function (Milne *et al.* 1999) as described earlier by Okuno & Nakada (2002). With their VP they obtained $1.88 \pm 0.7 \text{ mm yr}^{-1}$ and $-0.31 \pm 0.02 \mu\text{Gal yr}^{-1}$ for the uplift and gravity variation respectively, the uncertainties being estimated from three different deglaciation histories (Sato *et al.* 2006). We obtain a vertical velocity of $1.84 \text{ mm yr}^{-1}$ by running the SELEN software. Using the value of $-0.15 \mu\text{Gal mm}^{-1}$ obtained by Wahr *et al.* (1995) for the ratio between the viscous gravity variation and vertical velocity, we obtain a gravity rate of $-0.28 \mu\text{Gal yr}^{-1}$.

We now turn to the computation of the geodetic consequences of the PDIM. First, we use two models of ice coverage. The first is the SVAL model (Hagedoorn & Wolf 2003) in which the distribution of the 16 major ice-masses is based on the location of the glaciers described by Hagen *et al.* (1993). The thinning ice masses are approximated by co-axial elliptical cylinders that give a simple representation of the ice-covered area and topography. This model is used by Sato *et al.* (2006). The second model is based on the Digital Chart of the World (DCW). It provides a more realistic geographical location of the glaciers. The topography is provided by the GTOPO30 digital elevation model (http://edc.usgs.gov/products/elevation/gtopo30/gtopo30.html). A better location of the glaciers near Ny-Ålesund is given by the Digital Elevation Model (DEM) from the SPIRIT project (Korona *et al.* 2009).
<table><caption>**Table 2.** Viscosity profiles VP, V1 and VM2 used to compute the Earth response to deglaciation histories ICE-3G and ICE-5G. $\eta_u$ and $\eta_l$ are, respectively, the viscosities of the upper and lower mantles, and L is the thickness of the lithosphere. The values of L, $\eta_u$ and $\eta_l$ for the VP, V1 and VM2 profiles can be found, respectively, in the papers by Sato *et al.* (2006), Tushingham & Peltier (1991) and Braun *et al.* (2008).</caption><thead><tr><th>Name</th><th>Deglaciation model</th><th>η<sub>u</sub> (×10<sup>21</sup> Pa s)</th><th>η<sub>l</sub> (×10<sup>21</sup> Pa s)</th><th>L (km)</th></tr></thead><tbody><tr><td>VP</td><td>ICE-3G</td><td>0.5</td><td>10</td><td>100</td></tr><tr><td>V1</td><td>ICE-3G</td><td>1</td><td>2</td><td>120</td></tr><tr><td>VM2</td><td>ICE-5G</td><td>0.5</td><td>2.6</td><td>90</td></tr></tbody></table>

Secondly, we use three ice thinning-rate models. The first model has non-uniform thinning rates, deduced from Kohler *et al.* (2007) for the Ny-Ålesund area and from Dowdeswell *et al.* (2008) for the Austfonna ice cap. For all the other ice basins, we use non-uniform thinning rates given by an average of the long-term estimates of Nuth *et al.* (2010). The other two models are based on Sato *et al.* (2006) who proposed two uniform thinning rates of 47 and 75 cm we yr⁻¹. This provides six combined models of ice coverage, topography, and ice thinning. Their numbering is given in Table 4. Models 1 and 4 are shown in Fig. 5.

And thirdly, we compute the elastic deformation due to the PDIM by convolving the Green functions of the earth model PREM (Dziewonski & Anderson 1981), where we replace the global ocean by a solid crust (Table 3), with the ice-mass variation. We apply the spherical harmonics formalism of Farrell (1972) to compute the Green functions.

To better predict the gravity rates at Ny-Ålesund, we take the topography of the region into account. Indeed, in the vicinity of Ny-Ålesund, Fig. 6 shows large differences between the SVAL and DCW models for the areas located above and under the horizon of the station. We thus use Green's function for the Newtonian gravity variation as defined by Merriam (1992) and Boy *et al.* (2002) and applied to specific glaciers in the Alps by Mémin *et al.* (2009).

The results for the six combined models of ice coverage, topography and ice thinning are listed in Table 4.

## 3.2 Modelling of satellite gravity variation

Using the SELEN software, we compute the Earth response to two past-deglaciation histories, ICE-3G (Tushingham & Peltier 1991) and ICE-5G (Peltier 2004). We consider a spherical, non-rotating, radially stratified and incompressible earth model with a Maxwell viscoelastic rheological constitutive relation. The sea level equation is solved assuming that the shorelines are fixed. These approximations constitute a zeroth-order model (Spada & Stocchi 2007) that can be used for GRACE-like gravity variation estimates. The structure of the model is the following: the P-wave velocity is infinite, the density and S-wave velocity are taken from PREM, and the viscosity profiles VP, V1 and VM2 (Table 2) are radial profiles used by Sato *et al.* (2006), Tushingham & Peltier (1991) and Braun *et al.* (2008). VM2 is an average of the viscosity structure considered by Peltier (2004).

At the output, SELEN provides a spatial grid for the variation of the geoid height $\Delta N(\theta, \lambda, b)$, which can be expanded in spherical harmonics:

$$ \Delta N(\theta, \lambda, b) = \sum_{n=0}^{N_{\max}} \sum_{m=0}^{n} [\Delta N_{n}^{m,c} Y_{n}^{m,c}(\theta, \lambda) + \Delta N_{n}^{m,s} Y_{n}^{m,s}(\theta, \lambda)], \quad (4) $$

where b = 6371 km is the radius of the earth model. To keep consistency with the GRACE solutions, we take $N_{\max} = 50$. We have

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Models of combined non-uniform ice-thinning rate and ice-coverage extension from the Digital Chart of the World (top panel) and from the SVAL model (bottom panel). The unit is m we yr⁻¹.</figcaption></figure>

$$ \begin{aligned} \Delta N(\theta, \lambda, b) &= \frac{\Delta V^{\text{GIA}}(\theta, \lambda, b)}{g_0} \\ &= b \sum_{n=0}^{N_{\max}} \sum_{m=0}^{n} [\Delta V_{n}^{m,c} Y_{n}^{m,c}(\theta, \lambda) + \Delta V_{n}^{m,s} Y_{n}^{m,s}(\theta, \lambda)], \end{aligned} \quad (5) $$
<table><caption>**Table 3.** Density, seismic wave velocities $V_p$ and $V_s$ and quality factors $Q_μ$ and $Q_κ$ of the crust that replaces the ocean layer in the PREM model of Dziewonski & Anderson (1981).</caption><thead><tr><th>Density (kg m<sup>-3</sup>)</th><th>V<sub>p</sub> (m s<sup>-1</sup>)</th><th>V<sub>s</sub> (m s<sup>-1</sup>)</th><th>Q<sub>μ</sub></th><th>Q<sub>κ</sub></th></tr></thead><tbody><tr><td>2600</td><td>5800</td><td>3200</td><td>600</td><td>57823</td></tr></tbody></table>

where $\Delta V^{\text{GIA}}(\theta, \lambda, b)$ is the variation of the gravity potential induced by the GIA at b and $\Delta V_n^{m,c}$ and $\Delta V_n^{m,s}$ are its spherical harmonic coefficients. $g_0$ is the gravity at the surface of the earth model. The gravity variation at $a$, $\Delta g^{\text{GIA}}(\theta, \lambda, a)$, can then be recovered from

$$ \Delta g^{\text{GIA}}(\theta, \lambda, a) = \frac{GM}{a^2} \sum_{n=2}^{N_{\text{max}}} \sum_{m=0}^{n} (n+1) \left(\frac{b}{a}\right)^n \times \left[ \Delta V_n^{m,c} Y_n^{m,c}(\theta, \lambda) + \Delta V_n^{m,s} Y_n^{m,s}(\theta, \lambda) \right] \quad (6) $$

$$ = \frac{GM}{a^3} \sum_{n=2}^{N_{\text{max}}} \sum_{m=0}^{n} (n+1) \left(\frac{b}{a}\right)^{n-1} \times \left[ \Delta N_n^{m,c} Y_n^{m,c}(\theta, \lambda) + \Delta N_n^{m,s} Y_n^{m,s}(\theta, \lambda) \right], \quad (7) $$

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Topography of the DCW-derived (top panel) and SVAL models (bottom panel) split between ice areas above (left panel) and under (right panel) the horizon of Ny-Ålesund (red star). The unit is m.</figcaption></figure>

to which the cosine taper filter described in Section 2.2 is then applied.

The gravity variation associated to the deglaciation history ICE-3G and viscosity profile VP is shown in Fig. 4 (middle panel). The bottom chart in Fig. 4 shows the destriped residual signal of the mean of the three GRACE solutions corrected for the GIA contribution. Comparison with the top chart of Fig. 4, which shows the uncorrected mean GRACE signal, reveals that the mass-loss pattern over Svalbard is enhanced after correction for the GIA.

In Table 5, we list the residual signals for the GRACE solutions and their mean corrected for three GIA models, that is three couples deglaciation history—viscosity profile in the mantle, after application of the destriping filter of Swenson & Wahr (2006). The residual for the mean GRACE solution is lower than −0.97 µGal yr⁻¹ over Svalbard and the standard deviation is 0.11 µGal yr⁻¹. In Section 4, we will estimate the volume of current ice loss separately from space and ground gravity measurements. To compare the two estimates, we will use the same GIA correction, which is obtained by combining the deglaciation history ICE-3G and viscosity profile VP.
<table><caption>**Table 4.** Modelled vertical velocity *u* (in mm yr⁻¹), gravity rate *g*<sub>e</sub> (in µGal yr⁻¹) due to the elastic deformation, and gravity rates *g*<sub>nt</sub> and *g*<sub>n</sub> due to the Newtonian attraction of the load respectively with and without topography taken into account. The ratios of gravity rates and vertical velocity are in µGal mm⁻¹. The thinning rate in brackets is spatially variable.</caption><thead><tr><th>Ice coverage/topography</th><th>Thinning rate (m yr<sup>-1</sup>)</th><th>u</th><th>g<sub>e</sub></th><th>g<sub>nt</sub></th><th>g<sub>n</sub></th><th>(g<sub>e</sub> + g<sub>nt</sub>)/u</th><th>(g<sub>e</sub> + g<sub>n</sub>)/u</th></tr></thead><tbody><tr><td>1 – DCW/GTOPO30</td><td>[-0.7, -0.3]</td><td>1.51</td><td>-0.35</td><td>0.26</td><td>-0.06</td><td>-0.06</td><td>-0.27</td></tr><tr><td>2 – DCW/GTOPO30</td><td>-0.47</td><td>1.76</td><td>-0.39</td><td>0.24</td><td>-0.07</td><td>-0.09</td><td>-0.26</td></tr><tr><td>3 – DCW/GTOPO30</td><td>-0.75</td><td>2.80</td><td>-0.63</td><td>0.38</td><td>-0.09</td><td>-0.09</td><td>-0.26</td></tr><tr><td>4 – SVAL/SVAL</td><td>[-0.7, -0.3]</td><td>1.27</td><td>-0.28</td><td>-0.01</td><td>-0.05</td><td>-0.23</td><td>-0.26</td></tr><tr><td>5 – SVAL/SVAL</td><td>-0.47</td><td>2.03</td><td>-0.46</td><td>-0.02</td><td>-0.06</td><td>-0.24</td><td>-0.26</td></tr><tr><td>6 – SVAL/SVAL</td><td>-0.75</td><td>3.22</td><td>-0.73</td><td>-0.03</td><td>-0.10</td><td>-0.24</td><td>-0.26</td></tr></tbody></table>

<table><caption>**Table 5.** Mean residual signal (in µGal yr⁻¹) over Svalbard of the destriped GRACE-derived gravity variation between 2003 January and 2009 January corrected for three combined deglaciation history-viscosity profile (Table 2).</caption><thead><tr><th></th><th>ICE-3G + V1</th><th>ICE-3G + VP</th><th>ICE-5G + VM2</th><th>σ</th></tr></thead><tbody><tr><td>GRGS</td><td>-0.97</td><td>-1.16</td><td>-1.15</td><td>0.11</td></tr><tr><td>CSR</td><td>-0.91</td><td>-1.10</td><td>-1.09</td><td>0.11</td></tr><tr><td>GFZ</td><td>-1.04</td><td>-1.24</td><td>-1.23</td><td>0.11</td></tr><tr><td>Mean</td><td>-0.97</td><td>-1.17</td><td>-1.16</td><td>0.11</td></tr></tbody></table>

The gravity variation due to the PDIM is modelled with 0.25° × 0.25° grids of ice-thickness variations such that the integrated volumes of ice loss are approximately the same as the models 4 (13 km³ yr⁻¹), 5 (15 km³ yr⁻¹) and 6 (25 km³ yr⁻¹) listed in Table 4. An example of a grid is given on the left panel of Fig. 7. We expand the grids into series of spherical harmonics for which the coefficients are Δσ<sub>n</sub><sup>m,c</sup> and Δσ<sub>n</sub><sup>m,s</sup> and compute the gravity variation induced by the different ice-thinning rates according to

$$ \Delta g^{\text{PDIM}}(\theta, \lambda, b) = 4\pi G \sum_{n=2}^{N_{\text{max}}} \sum_{m=0}^{n} \frac{n+1}{2n+1} (1 + k'_{n}) \times [\Delta\sigma_{n}^{m,c} Y_{n}^{m,c}(\theta, \lambda) + \Delta\sigma_{n}^{m,s} Y_{n}^{m,s}(\theta, \lambda)], \quad (8) $$

where $k'_{n}$ is the load Love number of degree $n$ for the gravity potential. The cosine taper filter is also applied. As an example, the gravity variation due to the PDIM for model 4 is shown in Fig. 7 (right panel). The gravity rate is smaller than −1 µGal yr⁻¹, which

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Left panel: Model 4 of ice thinning given in Table 4. Right panel: gravity variation Δg<sup>PDIM</sup> due to the PDIM for the non-uniform thinning rate. The units are respectively m we yr⁻¹ and µGal yr⁻¹.</figcaption></figure>

is larger in absolute value than the gravity variation accompanying the GIA. As the gravity variations due to the GIA and PDIM are of opposite signs, the removal of the GIA signal from the GRACE solution enhances the remaining signal. Consequently, it really matters to accurately compute the GIA before correcting the GRACE data.

# 4 COMPARISON BETWEEN THEORY AND OBSERVATION

## 4.1 Ground measurements

Fig. 8 is a synthesis of all the computations made earlier. It shows the gravity variation as a function of the secular vertical velocity.

The slope of the black solid line passing trough the origin is −0.15 µGal mm⁻¹, which is the theoretical ratio between the gravity variation and vertical displacement for a viscoelastic deformation (Wahr *et al.* 1995). From the point on this solid line that corresponds to the viscoelastic deformation computed by using the model ICE-3G + VP (Section 3.1), we draw the black dashed line with a slope of −0.26 µGal mm⁻¹, which is the theoretical ratio between the gravity variation and vertical displacement for the elastic deformation generated by a surface loading over a spherical Earth model (de Linage *et al.* 2007). The gravity variation includes the direct Newtonian attraction of the load and the effect owing to the elastic deformation. In the specific case where the topography is neglected,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Observed and computed gravity rate as a function of the uplift rate. The solid black line gives the GIA computed by Wahr *et al.* (1995). Its slope is approximately −0.15 µGal mm⁻¹. The dashed black line corresponds to the theoretical gravity-to-displacement ratio for the elastic deformations computed by de Linage *et al.* (2007). Its slope is −0.26 µGal mm⁻¹. The computations listed in Table 4, to which GIA contribution is added, are represented by the different items combined with colours listed in the legend of the figure. The red and green lines respectively correspond to the ratio deduced for models 2, 3 and 5, 6 listed in Table 4 taking the topography into account. Error bars are that due to GIA modelling (Section 3.1). The rectangle gives the limits of the errors associated to the observations.</figcaption></figure>

the black dashed line is the sum of the GIA and PDIM effects (Mémin *et al.* 2010).

The SVAL rates are given by the green (with topography) and black (without topography) items. The rates of the DCW-derived models are the red (with topography) and blue (without topography) items. The squares and diamonds respectively correspond to the thinning rates of 47 and 75 cm we yr⁻¹. The circles are for a non-uniform thinning rate ranging from 30 to 70 cm we yr⁻¹. The red and green lines respectively correspond to the ratio deduced for models 2, 3 and 5, 6 listed in Table 4, taking the topography into account. They show how the topography provided by the SVAL and DCW models influences the gravity rate. The SVAL model has a smoother topography than the DCW model, which explains why the green line is closer than the red line to the black dashed line.

The green circle with a black contour and the dotted-dash rectangle represent the observations with their error bars. The vertical velocity (5.64 ± 1.57 mm yr⁻¹) was obtained in Section 2.3. The gravity rate (−1.02 ± 0.48 µGal yr⁻¹) was obtained in Section 2.1.

The best agreement between theory and observation is obtained for the modelling based on the SVAL model with a PDIM rate corresponding to a volume change of ~25 km³ we yr⁻¹ (75 cm we yr⁻¹) and the GIA computed with the ICE-3G + VP model.

## 4.2 GRACE data

Table 6 provides the residuals Δg sat − Δg GIA − Δg PDIM, where Δg sat is the gravity variation deduced from the mean of the GRGS, CSR and GFZ solutions (Section 2.2), Δg GIA is the GIA contribution for the ice-mass change history ICE-3G and viscosity profile VP (Table 2) and Δg PDIM is the PDIM contribution for models 4 (Fig. 7), 5 or 6 (Table 4). The residuals are computed as in Section 2.2. In Fig. 9, we show the residuals after applying the destriping filter.

<table><caption>**Table 6.** Mean residual signal (in µGal yr⁻¹) over Svalbard of the destriped GRACE-derived gravity variation between 2003 January and 2009 January corrected for the deglaciation history ICE-3G, viscosity profile VP and PDIM models 4, 5 and 6 (Table 4).</caption><thead><tr><th></th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>GRGS</td><td>−0.18</td><td>0.02</td><td>0.73</td></tr><tr><td>CSR</td><td>−0.12</td><td>0.08</td><td>0.78</td></tr><tr><td>GFZ</td><td>−0.25</td><td>−0.05</td><td>0.65</td></tr><tr><td>Mean</td><td>−0.18</td><td>0.03</td><td>0.72</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Mean gravity rates deduced from the CSR-RL04, GFZ-RL04 and GRGS-RL02 solutions destriped and corrected for the GIA and PDIM contributions. The CSR and GFZ solutions have been previously filtered with a cosine taper in the spectral domain between degrees 30 and 50. From top to bottom, the PDIM rates are respectively non-uniform, −47 and −75 cm we yr⁻¹. The unit is µGal yr⁻¹.</figcaption></figure>
The PDIM models 4 and 5 combined with the GIA model provide residuals that are respectively −0.19 and 0.02 µGal yr⁻¹, the PDIM model 6 and the GIA model give residuals that are larger than 0.65 µGal yr⁻¹. For the time span 2003 January to 2009 January, the GRACE data provide a total volume of ice melting of approximately 15 km³ we yr⁻¹. Given the volume involved in a PDIM model and the residuals, we can estimate the volume of ice loss needed to fit the gravity variation averaged from the GRGS, CSR and GFZ solutions. We obtain 15.5 ± 2.4 km³ we yr⁻¹. That volume loss is closer to the 8.8 ± 3 km³ we yr⁻¹ estimated by Wouters et al. (2008) for the period 2003 February to 2008 January than the 75 km³ we yr⁻¹ proposed by Chen et al. (2006b) for the period 2002 April to 2005 November.

The destriping filter may affect the estimated volume of ice loss. A similar analysis for the undestriped GRACE solutions leads to an ice loss of 9.1 ± 4.2 km³ we yr⁻¹.

## 5 DISCUSSION

As shown in Section 4, the ground observations of gravity variation and vertical displacement are explained by a PDIM of 25 km³ we yr⁻¹, whereas the GRACE-derived gravity variations provide a PDIM comprised between, roughly, 5 and 18 km³ we yr⁻¹, the upper limit being reached by the destriped solutions. This rather large interval includes estimates of ice loss by glaciology studies (Table 7).

Indeed, Nuth et al. (2010) obtained a total ice loss of 9.7 ± 0.55 km³ we yr⁻¹ over 27 000 km² by comparing altimetry from ICESat (2003–2007) to digital elevation models of 1965 and 1990. Dowdeswell et al. (2008) proposed a volume decrease comprised between −2.5 and −4.5 km³ we yr⁻¹ for the Austfonna ice cap, which represents 8000 km² of ice coverage. Using elevation measurements (GNSS surface profiles, airborne and ICESat laser altimetry) for the period 2002–2008, Moholdt et al. (2010a) proposed for Austfonna an ice melting of 1.3 ± 0.5 km³ we yr⁻¹. Therefore, the total volume loss is comprised between 11 and 14.2 km³ we yr⁻¹. Moholdt et al. (2010b) obtained from the analysis of ICESat laser altimetry data over the time interval 2003–2008 a total loss of 4.3 ± 1.4 km³ we yr⁻¹, excluding calving front retreat or advance.

Ground gravimeters are very sensitive to the location of the changing ice mass and, consequently, to the local effects. The ice-thinning rate of 75 cm we yr⁻¹ for the SVAL model is a simple model; a more realistic model such as the DCW model with an ice-thinning rate dependent on the altitude should be considered to improve the estimate of the gravity rate. Actually, the thinning rate is larger in the lowest parts of the glaciers (e.g. Kohler et al.

2007; Kääb 2008; Nuth et al. 2010). Moreover, the mass distribution of numerous surge-type glaciers in Svalbard changes in a very short time. Recently, Sund et al. (2009) reported that the Comfortlessbreen glacier, located south of Ny-Ålesund, started to surge in 2006. Such a process may influence the gravity rate observed at Ny-Ålesund, especially if it is located in the neighbourhood of the station. Contrary to the ground measurements, space gravimetric methods such as GRACE measurements are sensitive to the total mass loss and are not affected by local phenomena. The difference of sensitivity is likely responsible for the discrepancy between the volumes of ice loss deduced from GRACE and ground observations. Another ice-thinning distribution that, when integrated, would give the same ice mass loss as GRACE data could provide vertical velocity and gravity rate similar to the ground observations.

## 6 CONCLUSIONS

We have added two new AG measurements to the previous four observations made at Ny-Ålesund and obtained a new ground-observed gravity rate of −1.02 ± 0.48 µGal yr⁻¹, in agreement with the continuous SG series. Using the GRACE solutions provided by the CSR, GFZ and GRGS groups for the period 2003 January to 2009 January, we have derived the space gravity rate over Svalbard and Northern Europe. After filtering the CSR and GFZ solutions with a low-pass filter and showing similarities with the GRGS solutions, we have computed the mean of the destriped CSR, GFZ and GRGS solutions and obtained a gravity rate of −0.70 ± 0.07 µGal yr⁻¹ over Svalbard. The rate is somewhat lower (in absolute value) if we do not remove the stripes from the solutions: it is −0.52 ± 0.31 µGal yr⁻¹.

We have used the ICE-3G deglaciation history (Tushingham & Peltier 1991) combined with the viscosity profile used by Sato et al. (2006) to model the viscous response of the Earth to the last deglaciation. At this stage, the mean of the GRACE solutions corrected for the GIA contribution shows, over Svalbard, a negative gravity rate, which emphasizes the present ice mass loss: −1.17 ± 0.18 µGal yr⁻¹ for the destriped solutions, −0.99 ± 0.46 µGal yr⁻¹ if the stripes are not removed.

Next, we have computed the elastic response of the Earth to the PDIM. We have considered six models of PDIM by combining two models of ice-thinning extension (SVAL and modified DCW) and three melting rates (spatially non-uniform, −47 cm we yr⁻¹, −75 cm we yr⁻¹). The modelled GIA has been added to the elastic deformation due to the PDIM to compare the theoretical and observed geodetic parameters—the gravity variation given above and ground vertical velocity reported by Sato et al. (2006), Rülke et al.

<table><caption>**Table 7.** Summary of estimates of ice loss in Svalbard from various studies. The asterisk indicates that the destriping filter of Swenson & Wahr (2006) has been applied to the GRACE solutions.</caption><thead><tr><th>Data</th><th>Reference</th><th>Time interval</th><th>Volume loss (km³ yr⁻¹)</th></tr></thead><tbody><tr><td rowspan="4">GRACE</td><td>Chen et al. (2006b)</td><td>04.2002 to 11.2005</td><td>−75</td></tr><tr><td>Wouters et al. (2008)</td><td>02.2003 to 01.2008</td><td>−8.8 ± 3</td></tr><tr><td>This study</td><td>01.2003 to 01.2009</td><td>−9.1 ± 4.2</td></tr><tr><td>This study*</td><td>01.2003 to 01.2009</td><td>−15.5 ± 2.4</td></tr><tr><td>AG and SG</td><td>This study</td><td>1998–2007</td><td>−25</td></tr><tr><td>ICESat and DEM</td><td>Nuth et al. (2010)</td><td>1965/1990 to 2003/2007</td><td>−13.2 ± 1.55</td></tr><tr><td></td><td>Dowdeswell et al. (2008)</td><td>1966–2006</td><td></td></tr><tr><td>ICESat and DEM</td><td>Nuth et al. (2010)</td><td>1965/1990 to 2003/2007</td><td>−11 ± 1.05</td></tr><tr><td>ICESat</td><td>Moholdt et al. (2010a)</td><td>2002–2008</td><td></td></tr><tr><td>ICESat</td><td>Moholdt et al. (2010b)</td><td>2003–2008</td><td>−4.3 ± 1.4</td></tr></tbody></table>
(2008) and Kierulf et al. (2009a). The best agreement is obtained with the model of PDIM made up of the SVAL model and an ice melting rate of 75 cm we yr⁻¹, which gives an annual loss of ice of 25 km³ we. However, this does not agree with the smaller melting rate (from 5 to 18 km³ we yr⁻¹) derived from the GRACE solutions. Glaciology studies (Dowdeswell et al. 2008; Nuth et al. 2010; Moholdt et al. 2010a; Moholdt et al. 2010b) favour a volume of ice loss between 4 and 14.2 km³ we yr⁻¹, which is similar to the interval obtained from the GRACE data. A part of the discrepancy between the ice losses derived from ground and space observations can probably be reduced by taking into account the topography of the glaciers and the altitude dependency of ice melting in the modelling of PDIM.

## ACKNOWLEDGMENTS

We thank T. Sato for providing the SVAL model, E. Berthier for providing the DCW ice coverage and useful suggestions and J. Travelletti for creating an ice coverage mask of the Ny-Ålesund area. We thank two anonymous reviewers for their constructive comments. We acknowledge logistical and financial support from the French Polar Institute (Institut Paul-Emile Victor, IPEV). A. Mémin acknowledges financial support from the Centre National d'Études Spatiales. Part of this work was supported by COST Action ES0701 'Improved constraints on models of Glacial Isostatic Adjustment'.

## REFERENCES

Argus, D.F., 2007. Defining the translational velocity of the reference frame of Earth, *Geophys. J. Int.*, **169**, 830–838.

Barletta, V.R., Sabadini, R. & Bordoni, A., 2008. Isolating the PGR signal in the GRACE data: impact on mass balance estimates in Antarctica and Greenland, *Geophys. J. Int.*, **172**, 18–30.

Bettadpur, S., 2007. UTCSR level-2 processing standards document for level-2 product release 004, GRACE 327–742 (GR-CSR-03-03).

Boy, J.-P., Gegout, P. & Hinderer, J., 2002. Reduction of surface gravity data from global atmospheric pressure loading, *Geophys. J. Int.*, **149**, 534–545.

Braun, A., Kuo, C.-Y., Shum, C.K., Wu, P., van der Wal, W. & Fotopoulos, G., 2008. Glacial isostatic adjustment at the Laurentide ice sheet margin: models and observations in the Great Lakes region, *J. Geodyn.*, **46**, 165–173, doi: 10.1016/j.jog.2008.03.005.

Bruinsma, S., Lemoine, J.-M., Biancale, R. & Valès, N., 2010. CNES/GRGS 10-day gravity field models (release 2) and their evaluation, *Adv. Space Res.*, **45**, 587–601.

Chen, J.L., Tapley, B.D. & Wilson, C.R., 2006a. Alaskan mountain glacial melting observed by satellite gravimetry, *Earth planet. Sci. Lett.*, **248**, 368–378.

Chen, J.L., Wilson, C.R. & Tapley, B.D., 2006b. Satellite gravity measurements confirm accelerated melting of Greenland ice sheet, *Science*, **313**, 1958, doi:10.1126/science.1129007.

Chen, J.L., Wilson, C.R., Blankenship, D.D. & Tapley, B.D., 2006c. Antarctic mass rates from GRACE, *Geophys. Res. Lett.*, **33**, L11502, doi:10.1029/2006GL026369.

Chen, J.L., Wilson, C.R., Tapley, B.D., Blankenship, D.D. & Ivins, E.R., 2007. Patagonia Icefield melting observed by Gravity Recovery and Climate Experiment (GRACE), *Geophys. Res. Lett.*, **34**, L22501, doi:10.1029/2007GL031871.

Crossley, D. et al., 1999. Network of superconducting gravimeters benefits a number of disciplines, *EOS. Am. Geophys. Un.*, **80** (11), 125–126.

de Linage, C., Hinderer, J. & Rogister, Y., 2007. A search for the ratio between gravity variation and vertical displacement due to a surface load, *Geophys. J. Int.*, **171**, 986–994, doi:10.1111/j.1365-246X.2007.03613.x.

de Linage, C., Rivera, L., Hinderer, J., Boy, J.-P., Rogister, Y., Lambotte, S. & Biancale, R., 2009. Separation of coseismic and postseismic gravity

changes for the 2004 Sumatra-Andaman earthquake from 4.6 years of GRACE observations and modelling of the coseismic change by normal-modes summation, *Geophys. J. Int.*, **176**, 695–714.

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

Duan, X.J., Guo, J.Y., Shum, C.K. & van der Wal, W., 2009. On the post-processing removal of correlated errors in GRACE temporal gravity field solutions, *J. Geod.*, **83**, 1095–1106, doi:10.1007/s00190-009-0327-0.

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

Farrell, W.E., 1972. Deformation of the Earth by surface loads, *Rev. Geophys. Space Phys.*, **10**, 761–797.

Flechtner, F., 2007. GFZ level-2 processing standards document for level-2 product release 004, GRACE 327–743 (GR-GFZ-STD-001), Rev. 1.0.

Francis, O. et al., 2005. Results of the international comparison of absolute gravimeters in Walferdange (Luxembourg) of November 2003, in *IAG Symposia, Gravity, Geoid, and Space Missions, Vol. 129*, pp. 272–275, Springer, Berlin.

Hagedoorn, J.M. & Wolf, D., 2003, Pleistocene and recent deglaciation in Svalbard: implications for tide-gauche, GPS and VLBI measurements, *J. Geodyn.*, **35**, 415–423, doi: 10.1016/S0264-3707(03)00004-8.

Hagen, J.O., Liestøl, D., Roland, E. & Jørgensen, T., 1993, *Glacier Atlas of Svalbard and Jan Mayen*. Norsk Polarinstitutt, Oslo.

Howat, I.M., Joughin, I. & Scambos, T.A., 2007. Rapid Changes in Ice Discharge from Greenland Outlet Glaciers, *Science*, **315**, 1559–1561.

Kääb, A., 2008. Glacier volume changes using ASTER satellite stereo and ICESat GLAS laser altimetry. A test study on Edgeø ya, Eastern Svalbard, *IEEE Trans. Geosci. Remote Sens.*, **46**, 10, 2823–2830.

Kierulf, H.P., Pettersen, B.R., MacMillan, D.S. & Willis, P., 2009a. The kinematics of Ny-Ålesund from space geodetic data, *J. Geodyn.*, **48**, 37–46, doi: 10.1016/j.jog.2009.05.002.

Kierulf, H.P., Plag, H.-P. & Kohler, J., 2009b. Surface deformation induced by present-day ice melting in Svalbard, *Geophys. J. Int.*, **179**, 1–13.

Kohler, J., James, T.D., Murray, T., Nuth, C., Brandt, O., Barrand, N.E., Aas, H.F. & Luckman, A., 2007. Acceleration in thinning rate on western Svalbard glaciers, *Geophys. Res. Lett.*, **34**, L18502, doi:10.1029/2007GL030681.

Korona, J., Berthier, E., Bernard, M., Remy, F. & Thouvenot, E., 2009. SPIRIT. SPOT 5 stereoscopic survey of Polar Ice: reference images and topographies during the fourth International Polar Year (2007–2009), *ISPRS J. Photogrammetry Remote Sens.*, **64**(2), 204–212.

Luthcke, S.B., Zwally, H.J., Abdalati, W., Rowlands, D.D., Ray, R.D., Nerem, F.G., McCarthy, J.J. & Chinn, D.S., 2006. Recent Greenland ice mass loss by drainage system from satellite gravity observations, *Science*, **24**, doi: 10.1126/science.1130776.

Luthcke, S.B., Arendt, A.A., Rowlands, D.D., McCarthy, J.J. & Larsen, C.F., 2008. Recent glaciers mass changes in the Gulf of Alaska region from GRACE mascon solutions, *J. Glaciol.*, **54** (188), 2109–2112.

Matsuo, K. & Heki, K., 2010. Time-variable ice loss in Asian high mountains from satellite gravimetry, *Earth planet. Sci. Lett.*, **290**, 30–36.

Milne, G.A., Mitrovica, J.X. & Davis, J.L., 1999. Near-field hydro-isostacy: the implementation of a revised seal-level equation, *Geophys. J. Int.*, **139**, 464–482.

Mémin, A., Rogister, Y., Hinderer, J., Llubes, M., Berthier, E. & Boy, J.-P., 2009. Ground deformation and gravity variations modelled from present-day ice thinning in the vicinity of glaciers, *J. Geodyn.*, **48**(3–5), 195–203, doi: 10.1016/j.jog.2009.09.006.

Mémin, A., Rogister, Y. & Hinderer, J., 2010. Separation of the geodetic consequences of past and present ice melting, *Pure Appl. Geophys.*, submitted.

Merriam, J.B., 1992. Atmospheric pressure and gravity, *Geophys. J. Int.*, **109**, 488–500.

Moholdt, G., Hagen, J.O., Eiken, T. & Schuler, T.V., 2010a. Geometric changes and mass balance of the Austfonna ice cap, Svalbard, *Cryosphere*, **4**, 21–34.
Moholdt, G., Nuth, C., Hagen, J.O. & Kohler, J., 2010b. Recent elevation changes of Svalbard glaciers derived from ICESat laser altimetry, *Remote Sens. Environ.*, **114**(11), 2756–2767, doi:10.1016/j.rse.2010.06.008.

Nuth, C., Moholdt, G., Kohler, J., Hagen, J.O. & Kääb, A., 2010. Svalbard glacier elevation changes and contribution to sea level rise, *J. geophys. Res.*, **115**, F01008, doi:10.1029/2008JF001223.

Okuno, J. & Nakada, M., 2002. Contributions of ineffective ice load on sea-level and free-air gravity, ice sheet, sea level and the dynamic earth, *J. Am. Geophys. Un.*, Geodynamics Series 29, 177–185.

Peltier, W.R., 2004. Global Glacial Isostasy and the Surface of the Ice-Age Earth: the ICE-5G (VM2) Model and GRACE. *Ann. Rev. Earth. Planet. Sci.*, **32**, 111–149.

Quinn, K.J. & Ponte, R.M., 2010. Uncertainty in ocean mass trends from GRACE, *Geophys. J. Int.*, **181**, 762–768, doi:10.1111/j.1365-245X.2010.04508.x.

Robertsson, L. et al., 2001. Results from the fifth international comparison of absolute gravimeters, *ICAG97*, *Metrologia*, **38**, 71–78.

Rodell, M. & Famiglietti, J.S., 1999. Detectability of variations in continental water storage from satellite observations of the time dependent gravity field, *Water Resour. Res.*, **35**, 2705–2723.

Rodell, M. & Famiglietti, J.S., 2001. An analysis of terrestrial water storage variations in Illinois with implications for the gravity recovery and climate experiment (GRACE), *Water Resour. Res.*, **37**, 1327–1340.

Rülke, A., Dietrich, R., Fritsche, M., Rothacher, M. & Steigenberger, P., 2008. Realization of the Terrestrial Reference System by a re-processed global GPS network, *J. geophys. Res.*, **113**, B08409, doi:10.1029/2007JB005231.

Sato, T., Okuno, J., Hinderer, J., MacMillan, D.S., Plag, H.P., Francis, O., Falk, R. & Fukuda, Y., 2006. A geophysical interpretation of the secular displacement and gravity rates observed at Ny-Ålesund, Svalbard in the Arctic-effects of post-glacial rebound and present-day ice melting, *Geophys. J. Int.*, **165**, 729–743.

Slobbe, D.C., Ditmar, P. & Lindenbergh, R.C., 2009. Estimating the rates of mass change, ice volume change and snow volume change in Greenland from ICESat and GRACE data, *Geophys. J. Int.*, **176**, 95–106.

Spada, G. & Stocchi, P., 2007. SELEN: A Fortran 90 program for solving the 'sea level equation', *Comput. Geosci.*, **33**, 538–562.

Steffen, H., Denker, H. & Müller, J., 2008, Glacial isostatic adjustment in Fennoscandia from GRACE data and comparison with geodynamical models, *J. Geodyn.*, **46**, 155–164, doi: 10.1029/j.jog.2008.03.002.

Steffen, H., Gitlein, O., Denker, H., Müller, J. & Timmen, L., 2009, Present rate of uplift in Fennoscandia from GRACE and absolute gravimetry, *Tectonophysics*, **474**, 69–77, doi: 10.1016/j.tecto.2009.01.012.

Sund, M., Eiken, T., Hagen, J.O. & Kääb, A., 2009. Svalbard surge dynamics derived from geometric changes, *Ann. Glaciol.*, **50** (52), doi:10.3189/172756409789624265.

Swenson, S. & Wahr, J., 2006. Post-processing removal of the cor-
related errors in GRACE data, *Geophys. Res. Lett.*, **33**, L08402,
doi:10.1029/2006GL025285.

Tamisiea, M.E., Leuliette, E.W., Davis, J.L. & Mitrovica, J.X., 2005.
Constraining hydrological and cryospheric mass flux in southeastern
Alaska using space-based gravity measurements, *Geophys. Res. Lett.*,
*32*, L20501, doi:10.1029/2005GL023961.

Tregoning, P. & Lambeck, K., 2010. Origin of the ITRF: finding consistency between GPS and GIA models, *EGU General Assembly*, **12**, EGU2010-7352.

Tushingham, A.M. & Peltier, W.R., 1991. ICE-3G: a new global model of Late Pleistocene deglaciation based upon geophysical predictions of postglacial relative sea level change, *J. geophys. Res.*, **96**, 4497–4523.

Velicogna, I. & Wahr, J., 2006a. Measurements of time-variable gravity show mass loss in Antarctica, *Science*, **311**, 1754, doi:10.1126/science.1123785.

Velicogna, I. & Wahr, J., 2006b Acceleration of Greenland ice mass loss in spring 2004, *Science*, **1754**, doi:10.1038/nature05168.

Wahr, J., DaZhong, H. & Trupin, A., 1995. Predictions of vertical uplift caused by changing polar ice volumes on a viscoelastic earth, *Geophys. Res. Lett.*, **22**(8), 977–980.

Wouters, B., Chambers, D. & Schrama, E.J.O., 2008. GRACE observes small-scales mass loss in Greenland, *Geophys. Res. Lett.*, **35**, L20501, doi:10.1029/2008GL034816.