# Decadal geodetic variations in Ny-Alesund (Svalbard): role of past and present ice-mass changes

Anthony Mémin, Giorgio Spada, Jean-Paul Boy, Rogister Yves, Jacques

Hinderer

► To cite this version:

Anthony Mémin, Giorgio Spada, Jean-Paul Boy, Rogister Yves, Jacques Hinderer. Decadal geodetic variations in Ny-Alesund (Svalbard): role of past and present ice-mass changes. Geophysical Journal International, 2014, 198 (1), pp.285-297. 10.1093/gji/ggu134. hal-01262302

HAL Id: hal-01262302

https://hal.science/hal-01262302

Submitted on 9 Nov 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
# Decadal geodetic variations in Ny-Ålesund (Svalbard): role of past and present ice-mass changes

A. Mémin,<fnref n="1" /> G. Spada,<fnref n="2" /> J.-P. Boy,<fnref n="3" /> Y. Rogister<fnref n="3" /> and J. Hinderer<fnref n="3" />

Accepted 2014 April 3. Received 2014 April 2; in original form 2013 September 22

## SUMMARY

The geodetic rates for the gravity variation and vertical uplift in polar regions subject to past and present-day ice-mass changes (PDIMCs) provide important insight into the rheological structure of the Earth. We provide an update of the rates observed at Ny-Ålesund, Svalbard. To do so, we extract and remove the significant seasonal content from the observations. The rate of gravity variations, derived from absolute and relative gravity measurements, is $-1.39 \pm 0.11 \ \mu\text{Gal yr}^{-1}$. The rate of vertical displacements is estimated using GPS and tide gauge measurements. We obtain $7.94 \pm 0.21$ and $8.29 \pm 1.60 \ \text{mm yr}^{-1}$, respectively. We compare the extracted signal with that predicted by GLDAS/Noah and ECMWF Re-analysis (ERA)-interim hydrology models. We find that the seasonal gravity variations are well-represented by local hydrology changes contained in the ERA-interim model. The phase of seasonal vertical displacements are due to non-local continental hydrology and non-tidal ocean loading. However, a large part of the amplitude of the seasonal vertical displacements remains unexplained. The geodetic rates are used to investigate the asthenosphere viscosity and lithosphere/asthenosphere thicknesses. We first correct the updated geodetic rates for those induced by PDIMCs in Svalbard, using published results, and the sea level change due to the melting of the major ice reservoirs. We show that the latter are at the level of the geodetic rate uncertainties and are responsible for rates of gravity variations and vertical displacements of $-0.29 \pm 0.03 \ \mu\text{Gal yr}^{-1}$ and $1.11 \pm 0.10 \ \text{mm yr}^{-1}$, respectively. To account for the late Pleistocene deglaciation, we use the global ice evolution model ICE-3G. The Little Ice Age (LIA) deglaciation in Svalbard is modelled using a disc load model with a simple linear temporal evolution. The geodetic rates at Ny-Ålesund induced by the past deglaciations depend on the viscosity structure of the Earth. We find that viscous relaxation time due to the LIA deglaciation in Svalbard is more than 60 times shorter than that due to the Pleistocene deglaciation. We also find that the response to past and PDIMCs of an Earth model with asthenosphere viscosities ranging between 1.0 and $5.5 \times 10^{18} \ \text{Pa s}$ and lithosphere (resp. asthenosphere) thicknesses ranging between 50 and 100 km (resp. 120 and 170 km) can explain the rates derived from geodetic observations.

**Key words:** Satellite geodesy; Sea level change; Time variable gravity; Rheology: crust and lithosphere; Rheology: mantle; Arctic region.

## 1 INTRODUCTION

The Svalbard archipelago is located in the Arctic between 11°E and 26°E of longitude and 76°N and 81°N of latitude, close to a continental margin (Fig. 1). Most of the glaciers of this archipelago are thinning (e.g. Kohler et al. 2007; Dowdeswell et al. 2008; Kääb 2008; Moholdt et al. 2010; Nuth et al. 2010) which induces present-day deformation and gravity changes. Svalbard is also subject to the glacial isostatic adjustment (GIA) in response to past ice-mass changes such as those accompanying the late Pleistocene deglaciation.

This deglaciation terminated in Svalbard about 10 000 yr ago letting coasts and fjords ice free (e.g. Landvik et al. 1998).

Post-glacial land emergence data are vertical displacements of the solid Earth relative to sea level cumulated since the end of the Pleistocene deglaciation. They characterize the response of the Earth to past surface mass changes. Using these types of data makes it possible to characterize the rheology of the Earth. For example, post-glacial land emergence observations across Svalbard resulting from the Pleistocene deglaciation cannot be consistently explained by a single three-layer viscoelastic Earth model with a Maxwell

<footnote n="1">¹School of Physical Sciences, University of Tasmania, Private bag 37, TAS 7001, Hobart, Australia. E-mail: anthony.memin@utas.edu.au</footnote>

<footnote n="2">²Dipartimento di Scienze di Base e Fondamenti, Università degli Studi di Urbino ‘Carlo Bo’, I-61029 Urbino (PU), Italy</footnote>

<footnote n="3">³IPG/EOST, UMR 7516 (CNRS/Université de Strasbourg), 5 rue Descartes, F-67084 Strasbourg, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Bathymetry and topography (in metres) in the Svalbard area from the ETOPO1 global relief model (Amante & Eakins 2009). The yellow and red lines show the −1000 and −500 m isobath of the continental margin. The red star shows the location of Ny-Ålesund.</figcaption></figure>

rheology (Breuer & Wolf 1995; Kaufmann & Wolf 1996). This reflects possible lateral variations of the Earth rheological properties. They were attributed to changes in the lithosphere thickness and asthenosphere viscosity. Using a set of only six observations, Breuer & Wolf (1995) showed that the asthenosphere viscosity and lithosphere thickness vary from $3.8 \times 10^{18}$ Pa s and 80 km, close to the continental margin, to $2.0 \times 10^{19}$ Pa s and 140 km, 300 km inland from it. Adding 19 observations, Kaufmann & Wolf (1996) confirmed the existence of lateral variation of the asthenosphere viscosity with values ranging between $10^{18}$ Pa s, near the continental margin, and $10^{21}$ Pa s, further away. A lithosphere thickness above 110 km was suggested while it was not well-resolved by the observations. The existence of an asthenosphere with a low viscosity was also confirmed later by Hagedoorn & Wolf (2003) who compared vertical displacement rates estimated using geodetic techniques to that predicted using models of past and present-day ice-mass changes (PDIMCs).

Several geodetic observations are continuously operating in Ny-Ålesund, Svalbard (Fig. 1). Since 1999, a superconducting gravimeter (SG) has been installed which continuously recorded the gravity variations. Since 1991 and 1994, positioning observations have been collected from two global positioning system (GPS) receivers and one very long baseline interferometry antenna, respectively. Tide gauge data are also available since 1976 and from 1998 to 2012, eight campaigns of absolute gravity (AG) measurements were performed (1998, 2000, 2001, 2002, 2004, 2007, 2010, 2012). Geodetic measurements handled at Ny-Ålesund contain both the secular contributions from the past and PDIMCs. However, up to now, no studies have satisfactorily explained the decadal rates observed either in gravity variations (Mémin *et al.* 2011) or ground vertical motion measurements (Hagedoorn & Wolf 2003; Kierulf *et al.* 2009b; Omang & Kierulf 2011), or both simultaneously (Sato *et al.* 2006b; Omang & Kierulf 2011; Mémin *et al.* 2012). A possible explanation is that most of the previous studies neglected the effects of changes of ice mass on geodetic rates during the last century. For example, two cold periods occurred between 1180 and 1900 corresponding to the Little Ice Age (LIA) in Svalbard (Grove 2001). These cold periods were followed by a global shrinkage of the glaciers (Hagen & Liestøl 1990; Nuth *et al.* 2007). This deglaciation is likely

causing secular geodetic signals as is currently the case in Alaska (e.g. Sato *et al.* 2012). Hagedoorn & Wolf (2003) were the first to consider the LIA deglaciation effects on the vertical displacements at Ny-Ålesund. However, they did not look at the effects on the rate of the ground gravity variations. Geodetic rates induced by the LIA deglaciation can be used to characterize the lithosphere thickness and asthenosphere viscosity (e.g. Sato *et al.* 2012). Because Ny-Ålesund is close to a continental margin, according to Breuer & Wolf (1995) and Kaufmann & Wolf (1996), we expect that the asthenosphere viscosity is $\le 10^{19}$ Pa s.

Most of the glaciers are experiencing shrinking, and this results in a sea level change (SLC; Jacob *et al.* 2012; Gardner *et al.* 2013). SLC is not uniform over the oceans and regional variabilities exist (e.g. Cazenave & Llovel 2010). For example, ice melting across the Greenland ice sheet produces a decrease of the sea level around Svalbard (Spada *et al.* 2012b). This sea level fall should result in an uplift of the archipelago. The SLC induced by the ice-mass change of the major ice reservoirs deforms the Earth surface and its gravity field. Its geodetic effects are, however, usually not taken into account when predicting the geodetic rates due to PDIMC in Svalbard. An estimate of the geodetic rates induced by the SLC is required to assess its integration into the budget of the sources of deformation.

This paper investigates the geodetic consequences of the past and PDIMCs on geodetic rates observed at Ny-Ålesund. Specifically, it focuses on the geodetic contribution of the LIA deglaciation and its implication for determining the asthenosphere viscosity and lithosphere thickness. We first compute the rates of vertical displacement and gravity variations between 2000 and 2013 by extracting and removing seasonal signals from the geodetic time-series (Section 2). We subsequently compare the obtained seasonal signals with those expected from two hydrology models (Section 3). Then, we study the geodetic consequences of the past and PDIMCs (Section 4). We correct observed geodetic rates for PDIMC and SLC contributions (Section 4.1). Finally, we compute the geodetic effects due to the Pleistocene and LIA deglaciations (Section 4.2) and compare them to the residuals (observed rates minus PDIMC and SLC contributions) to investigate the asthenosphere viscosity and lithosphere thickness beneath Ny-Ålesund (Section 4.3). Our concluding remarks follow in Section 5.
## 2 GEODETIC OBSERVATIONS AT NY-ÅLESUND (2000–2013)

Continental hydrology significantly contributes to the seasonal signal contained in the time-series recorded by the SG (Sato *et al.* 2006a; Mémin *et al.* 2011). We expect that it also contributes to the signal recorded by the nearby NYAL and NYA1 GPS stations. To estimate decadal rates of gravity variations and surface velocities we must remove the hydrology signal from SG and GPS time-series.

### 2.1 Rate of gravity variations

Between 2000 and 2013, seven AG measurements (Table 1) have been made at the SG station in Ny-Ålesund using FG5 absolute gravimetre from 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. Luxembourgh). Measurements made in 2000, 2001 and 2002 are reported by Sato *et al.* (2006b), those of 2004 and 2007 by Mémin *et al.* (2011), and the last two are recent observations made in 2010 and 2012. To ensure a homogeneous processing of the raw data, we correct the additional two observations by adopting the same strategy as Sato *et al.* (2006b) and Mémin *et al.* (2011). The corrections account for three geophysical signals: the observed polar motion provided by the International Earth Rotation and Reference Systems Service (IERS), the observed tides, including the effect of ocean tidal loading, and the observed atmospheric pressure using a barometric admittance factor of −0.42 µGal hPa⁻¹ (Sato *et al.* 2006a). Using different absolute gravimetres can induce an uncertainty on the gravity value up to 2 µGal (Robertsson *et al.* 2001; Francis *et al.* 2005). Considering that three instruments were used at Ny-Ålesund, we add 2 µGal in quadrature to all the formal errors (Sato *et al.* 2006b; Mémin *et al.* 2011). The AG measurements (Fig. 2a) have not been made at the same period of the year (Mémin *et al.* 2011). To minimize the effects due to hydrology loading we must estimate the rate of gravity variations by taking the seasonal signal into account (van Camp *et al.* 2013). To do so, we make use of the SG time-series.

We analyse Ny-Ålesund SG data between 2000 January and 2012 June. We apply corrections for the three geophysical phenomena cited above. Contrary to the AG measurements, the change of gravity due to atmospheric pressure is corrected using the ERA-interim model from the European Centre for Medium-Range Weather Forecasts (ECMWF, Dee *et al.* (2011)) and the Toulouse Unstructured Grid Ocean model (TUGO-m). TUGO-m is an update of the previous 2-D gravity wave model MOG2D (Carrère & Lyard 2003). We do not expect that using a model instead of actual observations impacts on our results. Indeed, atmospheric pressure

<table><caption>**Table 1.** Values of absolute gravity measurements, after subtraction of 983 017 000 µGal, handled at Ny-Ålesund between 2000 and 2013.</caption><thead><tr><th>Date</th><th>Value (µGal)</th><th>1σ uncertainty (µGal)</th><th>Reference</th></tr></thead><tbody><tr><td>2000 July</td><td>55.80</td><td>6.10</td><td>Sato et al. (2006b)</td></tr><tr><td>2001 September</td><td>52.20</td><td>1.10</td><td>Sato et al. (2006b)</td></tr><tr><td>2002 September</td><td>51.00</td><td>1.10</td><td>Sato et al. (2006b)</td></tr><tr><td>2004 June</td><td>51.71</td><td>1.19</td><td>Mémin et al. (2011)</td></tr><tr><td>2007 June</td><td>52.10</td><td>1.02</td><td>Mémin et al. (2011)</td></tr><tr><td>2010 August</td><td>36.88</td><td>0.68</td><td>This study</td></tr><tr><td>2012 May</td><td>45.81</td><td>1.15</td><td>This study</td></tr></tbody></table>

Note: 1 µGal = 10 nm s⁻².

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

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Absolute (AG, red squares) and relative (SG, black curve) gravity measurements at Ny-Ålesund. The SG-filtered time-series is the red curve. (b) Seasonal gravity variations from the SG time-series (black curve), GLDAS/Noah (green) and ERA-interim (blue) hydrology models.</figcaption></figure>

corrections using local pressure measurements and a barometric admittance factor account for 90 per cent of the total atmospheric effect (Merriam 1992). Finally, we correct the SG time-series for 9 offsets larger than 10 µGal mostly due to the refilling of the SG dewar with liquid helium (Sato *et al.* 2006b; Mémin *et al.* 2011; Omang & Kierulf 2011).

To remove the seasonal signal from the AG data using the SG time-series, we adapt the method 3 of van Camp *et al.* (2013). The method 3 of Van Camp *et al.* (2013) provides similar uncertainties to their best method within 0.03 and 0.01 µGal yr⁻¹ on average over a 10 and 15 yr time span (see Table 2 therein). We fit

<table><caption>**Table 2.** Rates of gravity variations, $\dot{g}_{ag/sg}$, estimated from Ny-Ålesund gravity measurements.</caption><thead><tr><th>Reference</th><th>Time interval</th><th>ġ<sub>ag/sg</sub><br>(µGal yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Sato et al. (2006b)</td><td>1998–2002</td><td>-2.5 ± 0.9</td></tr><tr><td>Mémin et al. (2011)</td><td>1998–2007</td><td>-1.02 ± 0.48</td></tr><tr><td>Omang & Kierulf (2011)</td><td>2000–2010</td><td>-1.77 ± 0.01</td></tr><tr><td>This study</td><td>2000–2013</td><td>-1.39 ± 0.11</td></tr></tbody></table>
a degree 2 polynomial, solar annual and semi-annual signals from the SG data using the least squares method. We remove the estimated degree 2 polynomial from the SG time-series. The obtained residuals contain high-frequency and seasonal signals. To attenuate the high-frequency signals and estimate the seasonal signal, we run a monthly average on the residual time-series (Fig. 2b). The rate of gravity variations is estimated after correcting AG measurements for the seasonal signal extracted from the SG data. We obtain $\dot{g}_{ag/sg} = -1.39 \pm 0.11 \ \mu\text{Gal yr}^{-1}$, in agreement, within the uncertainties, with the trend estimated by Mémin et al. (2011) but lower than the trend determined by Sato et al. (2006b) and Omang & Kierulf (2011; Table 2). The superposition of SG and AG data shows a very good agreement (Fig. 2a).

## 2.2 Rate of vertical displacements

To determine the rate of vertical displacements at Ny-Ålesund, we use the GPS time-series solutions for stations NYAL and NYA1 from the Jet Propulsion Laboratory (JPL) from 2000 January to 2012 December (Fig. 3). They are publicly available from http://sideshow.jpl.nasa.gov/post/series.html. The GPS analysis is undertaken using the GPS Inferred Positioning System/Orbit Analysis and Simulation Software (GIPSY/OASIS) package in a precise point positioning mode (Zumberge et al. 1997). Standard solid Earth tide, ocean tide loading and pole tide loading corrections are applied according to the IERS2010 standards (Petit & Luzum 2010). In addition to the standard corrections, we correct for atmospheric pressure loading using the ERA-interim surface pressure fields assuming that the ocean responds to changes in the atmosphere as an inverted barometre (e.g. Van Dam & Wahr 1987; Petrov & Boy 2004). The resulting time-series are shown in Fig. 3. The station coordinates are expressed in the reference frame of the International GNSS Service, IGS08, which is derived from the International Terrestrial Reference Frame (ITRF) 2008 (Altamimi et al. 2011) by correcting the position of core sites for receiver antenna calibration updates (Rebischung et al. 2012).

We estimate the rate of vertical displacements performing a similar processing as for the SG measurements. We fit a model including an offset, a trend and a periodic components using the method of weighted least squares, with the vertical displacement uncertainties as weights. The periodic component includes terms to fit solar annual and semi-annual periods, as well as harmonics 1–6 of the

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Variations in the vertical coordinate of the NYAL (grey curve) and NYA1 (black curve) GPS stations, their filtered time-series are respectively the orange and red curves. Data are from the JPL.</figcaption></figure>

<table><caption>**Table 3.** Rates of vertical displacement, $\dot{u}_{gps}$ and $\dot{u}_{tg}$, estimated at Ny-Ålesund.</caption><thead><tr><th>Site</th><th>Time interval</th><th>u<sub>gps/tg</sub><br>(mm yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>NYAL</td><td>2000–2013</td><td>7.21 ± 0.12</td></tr><tr><td>NYA1</td><td>2000–2013</td><td>8.66 ± 0.09</td></tr><tr><td>Mean</td><td>2000–2013</td><td>7.94 ± 0.21</td></tr><tr><td>Tide gauge</td><td>2000–2012</td><td>8.29 ± 1.60</td></tr></tbody></table>

GPS draconitic year (351.4 d), shown to be a source of anomalous energy in GPS time-series (Ray et al. 2008; Tregoning & Watson 2009). Given the time-series length, it is not possible to resolve the solar and the draconitic annual terms simultaneously. However, to estimate a linear trend, we must attenuate any signal at these low frequencies. The vertical velocities of NYAL and NYA1 GPS stations (Table 3) are 7.21 ± 0.12 and 8.66 ± 0.09 mm yr$^{-1}$, respectively, leading to a mean uplift rate $\dot{u}_{gps} = 7.94 \pm 0.21 \ \text{mm yr}^{-1}$ (3σ), in IGS08 (Fig. 3). Our mean uplift rate is close to the estimate by Omang & Kierulf (2011) for NYA1 between late 2000 and early 2011. Expressing their results in ITRF2005, they obtained an uplift rate of 8.50 ± 0.04 mm yr$^{-1}$.

The vertical velocities estimated from GPS time-series depend on the choice of the terrestrial reference frame, processing strategy, measurement method and time interval considered (Sato et al. 2006b; Kierulf et al. 2009a,b; Mémin et al. 2011; Omang & Kierulf 2011). To assess the mean uplift rate estimated from GPS measurements, we use the monthly mean relative sea level (RSL) data (http://www.psmsl.org/) from Ny-Ålesund tide gauge measurements (Fig. 4, black curve). We estimate the rate of RSL changes, $\dot{s}$, between 2000 and 2012, considering a seasonal signal (Fig. 4, red curve). We obtain $\dot{s} = -5.51 \pm 1.34 \ \text{mm yr}^{-1}$. The tide gauge measures SLCs relative to a benchmark fixed to the ground. Therefore, it measures the changes in the difference between the rate of absolute sea level variation, $\dot{n}$, and the rate of vertical coordinate change of the benchmark, $\dot{u}_{tg}$:

$$ \dot{s} = \dot{n} - \dot{u}_{tg}. \qquad (1) $$

Absolute sea level variations are usually obtained from altimetric measurements (see, for example, Nerem et al. 2010). However, the orbit of most of the radar altimetres do not span high latitude areas because of their inclination, which is, for instance, 66° for

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Monthly mean RSL at Ny-Ålesund tide gauge (black) and its estimated model (red), containing a linear trend and a seasonal component. Data are from http://www.psmsl.org.</figcaption></figure>
Topex/Poseidon and its successors. We have to rely on an ocean general circulation model to predict the rate of absolute SLCs; we choose the GLORYS2v1 (GLObal Ocean ReanalYsis and Simula-tion) model (Ferry *et al.* 2012). This model is an eddy-resolving ocean model with spatial resolution of about 25 km with 75 layers, assimilating sea-surface temperature and elevation in ad-dition to temperature and salinity profiles. A sea-ice module is also included in GLORYS. At the location of the tide gauge, we obtain $\dot{n} = 2.78 \pm 0.26 \text{ mm yr}^{-1}$, leading to $\dot{u}_{tg} = \dot{n} - \dot{s} = 8.29 \pm 1.60 \text{ mm yr}^{-1}$, therefore confirming the uplift rate estimated from the GPS time-series.

# 3 SEASONAL GEODETIC VARIATIONS AT NY-ÅLESUND

## 3.1 Hydrology models

We compare the seasonal gravity variations and vertical displacements, observed by the SG and GPS, to global hydrology models. We use the classical high resolution (0.25° and 3 hr) GLDAS/Noah [Global Land Data Assimilation System, Rodell *et al.* (2004)] and the hydrology component of ERA-interim (Dee *et al.* 2011), with about 0.7° and 6 hr spatial and temporal samplings. Both models include 4-layer soil moisture estimates (0–2 m for GLDAS/Noah and 0–2.55 m for ERA-interim) and a snow layer. In addition, GLDAS/Noah also includes the water content stored in the canopy; this contribution is almost negligible compared to the soil-moisture and snow contributions. The forcings of the two models are different: ERA-interim hydrology is coupled to the atmospheric model, whereas the forcing data set of GLDAS/Noah combines various products.

## 3.2 Seasonal variations of gravity

The gravity seasonal signal is extracted from the SG data (black curve, Fig. 2b) as explained in Section 2.1. The phase of this signal is in very good agreement with that from the ERA-interim model (Fig. 2b). Specifically, the SG and ERA-interim seasonal signals have an annual period with amplitudes, from peak to trough, of up to 10–12 µGal. There are few discrepancies of 2–4 µGal in winter 2001, 2004, 2005, 2011 and 2012. These anomalies may be the result of episodic fluctuations in precipitations that are not adequately represented by the hydrology model. They may also be the result of small offsets that remain in the SG time-series because of our processing strategy. The seasonal signal from GLDAS/Noah has lower amplitudes, less than 1 µGal, than the SG data.

We further investigate the discrepancies between the gravity variations induced by the two hydrology models at Ny-Ålesund. Gravity effects can be decomposed into a local and non-local contribution. Assuming hydrology acts as a thin layer loading on a spherical Earth surface, the local contribution is modelled as a Bouguer Plate with a value of $\pm 0.42677 \text{ µGal cm}^{-1}$ (Boy & Hinderer 2006); the sign depends on the relative position of the gravimetre to the surface (above or below). For both hydrology models, we compute the standard deviation of the predicted seasonal gravity variations separating the local hydrology contribution from the non-local (Table 4). In the ERA-interim model, the local contribution is larger than the non-local contribution by more than one order of magnitude. This leads to a total signal with seasonal variations largely driven by local hydrology contribution. In the GLDAS/Noah model the non-local and local contributions have similar amplitudes but opposite signs so that the total contribution is lower than either of those. Gravity

<table><caption>**Table 4.** Standard deviation (µGal) of the gravity variations at Ny-Ålesund predicted by the GLDAS/Noah and ERA-interim hydrology models between 2000 and 2013. We separate the non-local and local contributions.</caption><thead><tr><td>Contribution</td><td>GLDAS/Noah</td><td>ERA-interim</td></tr></thead><tbody><tr><td>Non-Local</td><td>0.31</td><td>0.13</td></tr><tr><td>Local</td><td>0.33</td><td>3.16</td></tr><tr><td>Total</td><td>0.23</td><td>3.22</td></tr></tbody></table>

variations being sensitive to local mass changes (e.g. Mémin *et al.* 2009), and given the agreement between the ERA-interim model and the SG seasonal variations, we conclude that the ERA-interim model adequately represents the local hydrological effects at Ny-Ålesund. In contrast, hydrology signal in the GLDAS/Noah model seems not well represented as it does not fully capture the local contribution at Ny-Ålesund.

## 3.3 Seasonal vertical motions

To estimate the seasonal component of vertical displacements at NYAL and NYA1, we remove the linear trends estimated in Section 2.2 from the time-series plotted in Fig. 3. Contrary to the gravity variations, the residuals for the vertical displacements are dominated by a large high-frequency content requiring a stronger filtering. To retrieve the seasonal signal we attenuate high-frequency signals of the residual time-series using a running average of 3 months (Figs 3 and 5).

Fig. 5(a) shows the seasonal signal extracted both from the GPS time-series (grey and black curves) and from the two hydrology models ERA-interim and GLDAS/Noah (blue and green curves). No agreement is obtained between the amplitude of the observed and predicted seasonal signals. Besides, there is an offset between the phases of the GPS- and model-derived seasonal signals. For NYA1 and NYAL, the correlation coefficients between the observed and predicted seasonal variations are 0.68 and 0.75, for the GLDAS/Noah model, and 0.09 and 0.11, for the ERA-interim model, respectively (Table 5). This suggests that the phases of the seasonal signal predicted from the GLDAS/Noah model are in better agreement with the observations than those from the ERA-interim model. This is confirmed by comparing the variance reduction of the observed seasonal variations corrected for either models. Correcting the observed seasonal variations for the hydrology signal of the GLDAS/Noah model reduces the variance by up to 27–30 per cent while using the ERA-interim model the reduction is less than 1.5 per cent.

The vertical displacements, as recorded by GPS, mainly reflect the deformation induced by loading with a large spatial fingerprint. Since Ny-Ålesund is located on the coast (Fig. 1), GPS are likely measuring the effect of ocean loading. To further investigate the discrepancies mentioned in the previous paragraph, we use TUGO-m to compute the effect of non-tidal ocean loading (NTOL) on NYA1 and NYAL seasonal variations. Adding the vertical displacement due to NTOL to that induced by hydrology leads to a better agreement between the observed and predicted seasonal signals (Fig. 5b). The coefficients of correlation increase to 0.85 and 0.86, for the GLDAS/Noah model, and 0.44 and 0.37, for the ERA-interim model, at NYA1 and NYAL, respectively (Table 5). Considering TUGO-m further reduces the variances of the observed seasonal signals by up to 10 per cent. The total reduction reaches up to 37 and 8–11 per cent for GLDAS/Noah and ERA-interim models,
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Detrended and filtered vertical displacement time-series of NYAL (grey curve) and NYA1 (black curve) GPS stations and seasonal vertical displacements from the GLDAS/Noah (green) and ERA-interim (blue) hydrology models. (b) Same as (a) but the green and blue curves are the sum of the vertical displacements from the GLDAS/Noah and ERA-interim hydrology models with the non-tidal ocean loading computed using TUGO-m.</figcaption></figure>

respectively. It confirms that the phases of the seasonal variations obtained by combining the GLDAS/Noah model and TUGO-m better agree with those of the seasonal vertical movements estimated from the two GPS stations. However, the amplitudes predicted by the models are still much smaller than those observed. We suspect that the remaining discrepancies are partly due to instrumental- and processing-related issues. It may also be due to mis-represented signals in the hydrology models, in particular the snow component (Sato et al. 2006a).

Hydrology and non-tidal ocean general circulation models have some issues at high latitudes. Precipitation data are not as accurate as rainfall data for mid and low latitudes data, due to poorer satellite coverage (the Tropical Rainfall Measuring Mission, for example, is only covering the ±45° area). There is also no radar altimetry from Topex/Poseidon & Jason for latitudes higher than 66°. On the other hand, we assume a spherically symmetric Earth model for the Green's function computation. A different crustal model leads to changes in the Green's function, but in general its effects are small

<table><caption>**Table 5.** Coefficient of correlation between the observed and predicted seasonal vertical displacements at Ny-Ålesund between 2002 and 2011. The variances (mm²) of the observed seasonal vertical displacement after applying the different predictions are shown in parentheses. The variances of NYA1 and NYAL seasonal vertical displacement are 12.6 and 13.6 mm², respectively.</caption><thead><tr><td>Model</td><td>NYA1</td><td>NYAL</td></tr></thead><tbody><tr><td>GLDAS/Noah</td><td>0.68 (9.2)</td><td>0.75 (9.5)</td></tr><tr><td>GLDAS/Noah + TUGO-m</td><td>0.85 (7.9)</td><td>0.86 (8.6)</td></tr><tr><td>ERA-interim</td><td>0.09 (12.5)</td><td>0.11 (13.4)</td></tr><tr><td>ERA-interim + TUGO-m</td><td>0.44 (11.2)</td><td>0.37 (12.5)</td></tr></tbody></table>

on loading estimates, compared to the errors of the hydrology models. The local gravity contribution computation (here the Bouguer approximation) can certainly be improved by taking into account the topography in the close neighbourhood of the station (a few hundred of metres). However, the main limitation still remains the errors in the hydrology model at high latitudes.

The gravity variations mainly reflect local changes while vertical displacements are sensitive to changes on larger spatial scales. Although the amplitudes do not match properly, vertical displacements obtained from GPS estimates are strongly correlated with those predicted by the GLDAS/Noah model combined with TUGO-m. This shows that the GLDAS/Noah better represents the non-local hydrology contribution when combined with TUGO-m than ERA-interim.

## 4 GEODETIC RATES DUE TO ICE-MASS CHANGES

In Section 2, we have derived rates of gravity variations and changes in the vertical coordinates of NYAL and NYA1 GPS stations. Here, we address the issue of their interpretation as a response to ice-mass changes.

In previous studies, the geodetic effects induced by the Pleistocene deglaciation (Hagedoorn & Wolf 2003; Sato et al. 2006b; Kierulf et al. 2009b) were computed. They were combined with estimates of the effects due to PDIMC using simple (Hagedoorn & Wolf 2003; Sato et al. 2006b; Mémin et al. 2011) and realistic (Kierulf et al. 2009b; Mémin et al. 2012) models. In this section, we remove from the updated geodetic rates the effects of PDIMC in Svalbard and SLC resulting from the melting of the main glaciated regions (Section 4.1). Then we investigate the contributions of the deformation induced by the Pleistocene and LIA deglaciations in the residuals (Section 4.2). Estimates of asthenosphere viscosity and lithosphere thickness are finally provided in Section 4.3.

### 4.1 Geodetic consequences of PDIMC

The Earth elastically deforms due to the rapid change of surface loads exerted by the redistribution of mass on its surface. To account for the deformation induced by the PDIMC at Ny-Ålesund, we use the rates of vertical displacements, $\dot{u}_{pdime}$, and gravity variations, $\dot{g}_{pdime}$, computed by Mémin et al. (2012) for an elastic Earth model. Based on published altitude profiles of ice-mass change (Kierulf et al. 2009b; Kohler et al. 2007; Moholdt et al. 2010; Nuth et al. 2010), Mémin et al. (2012) have defined three models (hereafter referred to as A, B and C) which account for the spatial and altitude dependency of the change of mass of ice in Svalbard. The difference between model A and models B and C is the total volume of ice loss of glaciers located between 2 and 110 km from Ny-Ålesund.
<table><caption>**Table 6.** Rates of vertical displacements ($\dot{u}_{pdimec}$) and gravity variations ($\dot{g}_{pdimec}$) at Ny-Ålesund induced by the change of ice mass in Svalbard. $\dot{V}$ is the rate of ice-volume change of PDIMC models A, B and C described in the text.</caption><thead><tr><th>Model</th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td>˙V (km<sup>3</sup> yr<sup>-1</sup>)</td><td>-11.64</td><td>-12.22</td><td>-15.67</td></tr><tr><td>˙u<sub>pdimec</sub> (mm yr<sup>-1</sup>)</td><td>1.88</td><td>2.09</td><td>3.52</td></tr><tr><td>˙g<sub>pdimec</sub> (µGal yr<sup>-1</sup>)</td><td>-0.10</td><td>-0.43</td><td>-0.39</td></tr></tbody></table>

<table><caption>**Table 7.** Vertical displacement rates predicted at Ny-Ålesund, $\dot{u}_{slc}$, due to the sea level change induced by the melting of the main ice sheets and glaciers. The regional mass-loss rate, $\dot{M} = -100$ Gt yr$^{-1}$, is uniformly distributed across the corresponding region. $\dot{M}_{sca}$ is the rescaled mass-loss rate according to Gardner *et al.* (2013) and Shepherd *et al.* (2012). The figures in brackets are the percentage of the total SLC-induced vertical displacement.</caption><thead><tr><th>Melting area</th><th>˙u<sub>slc</sub><br>using M<br>(mm yr<sup>-1</sup>)</th><th>Ṁ<sub>sca</sub><br>(Gt yr<sup>-1</sup>)</th><th>˙u<sub>slc</sub><br>using Ṁ<sub>sca</sub><br>(mm yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Canada</td><td>0.28</td><td>-60 ± 8</td><td>0.17 ± 0.02 (15 per cent)</td></tr><tr><td>Greenland</td><td>0.38</td><td>-232 ± 23</td><td>0.88 ± 0.09 (79 per cent)</td></tr><tr><td>Russian Arctic</td><td>0.35</td><td>-11 ± 4</td><td>0.04 ± 0.01 (4 per cent)</td></tr><tr><td>Alaska</td><td>0.11</td><td>-64 ± 20</td><td>0.07 ± 0.02 (6 per cent)</td></tr><tr><td>Europe</td><td>0.22</td><td>-15 ± 2</td><td>0.03 ± 0.00 (3 per cent)</td></tr><tr><td>Southern Andes</td><td>-0.06</td><td>-29 ± 10</td><td>-0.02 ± 0.01 (2 per cent)</td></tr><tr><td>Antarctica</td><td>-0.08</td><td>-72 ± 43</td><td>-0.06 ± 0.03 (5 per cent)</td></tr><tr><td>Total</td><td>1.19</td><td>-483 ± 54</td><td>1.11 ± 0.10</td></tr></tbody></table>

Models B and C use the same altitude profiles of ice-mass change except for the glaciers between 10 and 110 km from Ny-Ålesund. These three models are in agreement with the change of ice volume, $\dot{V} = -9.1 \pm 4.2$ km$^3$ yr$^{-1}$, estimated by Mémin *et al.* (2011) from the Gravity Recovery And Climate Experiment (GRACE) solutions. However, they lead to different geodetic rates (Table 6). The uplift rate predicted using model C is 13 per cent larger than the estimates by Kierulf *et al.* (2009b) and Omang & Kierulf (2011). However, the gravity rate is 50 per cent larger than the rate they obtained. Indeed, Omang & Kierulf (2011) neglected the effect of the topography of glaciers which has been shown by Mémin *et al.* (2012) to significantly impact on the predicted gravity rates.

The Svalbard archipelago is not the only place where glaciers are losing mass. Most of glaciers on Earth are shrinking (Gardner *et al.* 2013) which results in a spatially non uniform SLC. We estimate the deformation induced by the change in the sea level due to the melting of the major reservoirs of ice. The results for the rate of vertical displacements are given in Table 7. They are obtained assuming a regional mass loss of 100 Gt yr$^{-1}$ from all the glaciated regions considered (this rate of melting is equivalent to a sea-level rise of 0.28 mm yr$^{-1}$). The corresponding surface load is assumed to be uniformly distributed across each region, which is appropriate in view of the large distance separating the mass loss sources from the Svalbard archipelago. In particular, the rates of uplift (or subsidence) are computed by taking the time derivative of

$$u_{slc}(\theta, \lambda, t) = \rho_i G_u *_i I + \rho_w G_u *_o S, \quad (2)$$

where $(\theta, \lambda)$ are colatitude and longitude of Ny-Ålesund, $\rho_i$ and $\rho_w$ are the densities of ice and water, $I(\theta', \lambda', t')$ and $S(\theta', \lambda', t')$ are the ice thickness variation and the relative SLC at point $(\theta', \lambda')$ and time $t'$ and symbols $*_i$ and $*_o$ denote spatio-temporal convolutions over

the melting ice sources and the global oceans, respectively, and $G_u$, Green's function for vertical displacement (see Spada & Stocchi 2007). In eq. (2) above, the first term accounts for the direct effect of ice loading, while the second represents the contribution of the meltwater load. The relative SLC $S$ has been evaluated solving the so-called 'sea level equation; (SLE, see Farrell & Clark 1976) by means of an improved version of the open source code SELEN that includes a module for the elastic deformations that can use user-supplied compressible Love numbers (Spada & Stocchi 2007; Spada *et al.* 2012a). Hence, $S$ is not assumed to be uniform throughout the oceans. In all our computations, the mantle is assumed to be elastic and compressible (the load-deformation coefficients pertain to the Preliminary Reference Earth Model (PREM) of Dziewonski & Anderson 1981), the rotational feedback on sea-level is taken into account following the theory of Milne & Mitrovica (1998), the fixed-shorelines approximation is used and the geodetic variations are expressed in the reference frame of the center of mass of the whole Earth. The maximum harmonic degree of the analysis is $l_{max} = 128$.

The results obtained for an ice-mass loss of 100 Gt yr$^{-1}$ are rescaled according to estimates of the mass balances for each region (Gardner *et al.* 2013; Shepherd *et al.* 2012). The rate of the vertical displacement at Ny-Ålesund due to the SLC resulting from the melting of the major ice reservoirs is $\dot{u}_{slc} = 1.11 \pm 0.10$ mm yr$^{-1}$ (Table 7). Ice-mass changes from Canada and Greenland account for 94 per cent of this uplift. We estimate the rate of gravity variations induced by the SLC at Ny-Ålesund using the gravity variation-to-vertical displacement conversion factor, $-0.26$ µGal mm$^{-1}$, proposed by de Linage *et al.* (2007) and shown to be valid here (Mémin *et al.* 2012). Indeed, Ny-Ålesund is located outside the ocean where variations are assumed to be within a thin layer ($\le 1$ mm yr$^{-1}$) over the surface of a spherical Earth. This leads to a SLC-induced gravity variation rate of $\dot{g}_{slc} = -0.29 \pm 0.03$ µGal mm$^{-1}$.

Table 8 summarizes the predicted rates of vertical displacements and gravity variations induced by PDIMC and SLC at Ny-Ålesund. It also gives the residual rates obtained from correcting the observed displacement and gravity rates for these contributions :

$$\dot{g}_{res} = \dot{g}_{ag/sg} - (\dot{g}_{pdimec} + \dot{g}_{slc}), \quad (3)$$

$$\dot{u}_{res} = \dot{u}_{gps/tg} - (\dot{u}_{pdimec} + \dot{u}_{slc}). \quad (4)$$

The large residual rates range between 2.06 and 6.90 mm yr$^{-1}$ for the vertical displacement and $-1.11$ and $-0.56$ µGal yr$^{-1}$ for the gravity variations, depending on the model. In the next two subsections, we investigate how the present-day geodetic consequences of the past ice-mass changes can explain the derived residual rates.

## 4.2 Geodetic consequences of past ice-mass changes

To account for the deformation induced by the past ice-mass changes, we consider both the late Pleistocene and LIA deglaciations. The space and time evolution of the former ice-sheets during the Pleistocene deglaciation follow the history contained in the global ice model ICE-3G (Tushingham & Peltier 1991). The LIA history in Svalbard is not precisely known (Hagen & Liestøl 1990; Grove 2001; Nuth *et al.* 2007), so we use qualitative descriptions for starting and ending the glaciation and deglaciation phases during the LIA. Grove (2001) states that two cold periods occurred between 1180 and 1500 and between 1700 and 1900. Hagen & Liestøl (1990) mention that the maximum glaciers extension in Svalbard may have been reached by the end of the 19th or the beginning of the 20th
<table><caption>**Table 8.** Residual rates obtained after correcting the observed rates for predicted rates of vertical displacements ($\dot{u}_{res}$ in mm yr$^{-1}$) and gravity variations ($\dot{g}_{res}$ in µGal yr$^{-1}$) at Ny-Ålesund due to present-day ice-mass change (PDIMC) in Svalbard and sea level change (SLC) induced by the melting of the main ice sheets and glaciers.</caption><thead><tr><th rowspan="2">PDIMC model</th><th colspan="2">A</th><th colspan="2">B</th><th colspan="2">C</th></tr><tr><th>u̇<sub>pdimc</sub> + u̇<sub>slc</sub></th><th>ġ<sub>pdimc</sub> + ġ<sub>slc</sub></th><th>u̇<sub>pdimc</sub> + u̇<sub>slc</sub></th><th>ġ<sub>pdimc</sub> + ġ<sub>slc</sub></th><th>u̇<sub>pdimc</sub> + u̇<sub>slc</sub></th><th>ġ<sub>pdimc</sub> + ġ<sub>slc</sub></th></tr></thead><tbody><tr><td>PDIMC+SLC</td><td>2.99 ± 0.10</td><td>-0.38 ± 0.03</td><td>3.20 ± 0.10</td><td>-0.71 ± 0.03</td><td>4.63 ± 0.10</td><td>-0.67 ± 0.03</td></tr><tr><td><strong>Residuals</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GPS</td><td>4.95 ± 0.23</td><td>—</td><td>4.74 ± 0.23</td><td>—</td><td>3.31 ± 0.23</td><td>—</td></tr><tr><td>Tide Gauge</td><td>5.30 ± 1.60</td><td>—</td><td>5.09 ± 1.60</td><td>—</td><td>3.66 ± 1.60</td><td>—</td></tr><tr><td>Gravity</td><td>—</td><td>-1.00 ± 0.11</td><td>—</td><td>-0.67 ± 0.11</td><td>—</td><td>-0.71 ± 0.11</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Discretization of the ice-covered area, bounded by the blue line, during the LIA by means of surface disc loads. Red dots mark the centers of discs with an area of 32 km² each.</figcaption></figure>

century. Referring to Nordli & Kohler (2003), Nuth et al. (2007) characterize the retreat of Svalbard glaciers as starting sometime after 1920. To take into account the estimate of ice-mass changes of the glaciers between 1936–1938 and 1990 by Nuth et al. (2007), we assume that the maximum extension of the glaciers occurred in 1930. We also assume that from 1700 to 1930 the ice thickness has been linearly increasing over the 1091 discs of 32 km² forming the ice-covered area (Fig. 6). In 1930, the glaciers started to shrink and reached their current configuration in 2000, loosing their ice layer of thickness $h_{max}$ gained during the glaciation phase. Nuth et al. (2007) show that the mean annual ice-mass balance over seven regions covering 5123 km² (14 per cent of the total glaciated area in Svalbard) is about –34 cm yr$^{-1}$. It ranges between –16 and –63 cm yr$^{-1}$ depending on the region. Lacking of more details for the ice-height changes across Svalbard, we choose $h_{max}$ according to the study by Nuth et al. (2007). Considering that the deglaciation phase is 70 yr long, we obtain $h_{max} = 23.8$ m.

We use a modified version of the open source TABOO software (Spada 2003; Spada et al. 2003) to compute the current rates of the vertical displacements, $\dot{u}_v$, and gravity variations, $\dot{g}_v$, induced by the viscous deformation due to past ice-mass changes at Ny-Ålesund. The modification consists in including the computation of the gravity variations and their rate using the gravity Green functions (Peltier 1974). Vertical displacement rates are only slightly affected by using either a compressible or an incompressible Earth model (Tanaka et al. 2011). We use an incompressible,

<table><caption>**Table 9.** Earth models used to compute the vertical displacement and gravity variations induced by past-ice mass changes at Ny-Ålesund. Notation: T and $\eta$ correspond to thickness (in km) and viscosity (in $10^{21}$ Pa s), respectively. The core-mantle, upper-lower mantle and mantle-asthenosphere boundaries are at depths of 2891, 670 and 220 km, respectively.</caption><thead><tr><th rowspan="2">Model<br>Parameter</th><th colspan="2">EM-A</th><th colspan="2">EM-B</th><th colspan="2">EM-C</th></tr><tr><th>T</th><th>η</th><th>T</th><th>η</th><th>T</th><th>η</th></tr></thead><tbody><tr><td>Lithosphere</td><td>50–150</td><td>∞</td><td>50–150</td><td>∞</td><td>50–150</td><td>∞</td></tr><tr><td>Asthenosphere</td><td>170–70</td><td>10<sup>-3</sup>–10</td><td>170–70</td><td>10<sup>-3</sup>–10</td><td>170–70</td><td>10<sup>-3</sup>–10</td></tr><tr><td>Upper mantle</td><td>450</td><td>0.5</td><td>450</td><td>1.0</td><td>450</td><td>0.5</td></tr><tr><td>Lower mantle</td><td>2,221</td><td>2.7</td><td>2,221</td><td>2.7</td><td>2,221</td><td>10.0</td></tr></tbody></table>

viscoelastic and spherically symmetric Earth model with four layers between the core-mantle boundary and the surface. The layers correspond to the lithosphere, asthenosphere, upper and lower mantles (Table 9). The radii of the core-mantle and mantle-asthenosphere boundaries are 3480 and 6151 km, respectively. The radius of the boundary between the upper and lower mantles is 5701 km. The thicknesses of the asthenosphere, $T_a$, and the lithosphere, $T_l$, are free parameters that vary simultaneously so that when the lithosphere thickness increases (resp. decreases) the thickness of the asthenosphere decreases (resp. increases). Lithosphere and asthenosphere thicknesses range between 50 and 150 km and 170 and 70 km, respectively. The elastic properties of the Earth models are volume averages of PREM (Dziewonski & Anderson 1981) and their viscoelastic rheology is Maxwellian. We consider three Earth models with similar elastic parameters but different mantle viscosities (Table 9). Models EM-A and EM-B have a lower mantle viscosity of $2.7 \times 10^{21}$ Pa s while for model EM-C it is $10^{22}$ Pa s. The viscosity of the upper mantle of models EM-A and EM-C is $5.0 \times 10^{20}$ Pa s, it is $10^{21}$ Pa s for model EM-B. Mantle viscosities for model EM-A correspond to a volume-averaged estimate of the viscosity structure (VM2) used by Peltier (2004). Upper and lower mantle viscosities of models EM-B and EM-C respectively are perturbations of those of model EM-A. They are employed here to assess the sensitivity of the geodetic rates computed at Ny-Ålesund to changes in the mantle viscosity structure. For all the models, the asthenosphere viscosity, $\eta_a$, ranges between 0.001 and $10 \times 10^{21}$ Pa s.

For a given load history and a set of $\eta_a$, $T_a$ and $T_l$, modifying the mantle viscosities leads to several $\dot{u}_v$ and $\dot{g}_v$. Further variations of $\eta_a$, $T_a$ and $T_l$ result in additional changes of $\dot{u}_v$ and $\dot{g}_v$. Thus for a specific load history, changes in $\dot{u}_v$ and $\dot{g}_v$ reflect the sensitivity of the predicted geodetic rates to the rheological parameters. Therefore, we use the three Earth models listed in Table 9 to compute the variance, $\chi^2(\eta_a, T_a, T_l)$, of $\dot{u}_v$ and $\dot{g}_v$ induced by the Pleistocene and LIA deglaciations. We use the following expression

$$ \chi^2 = \frac{1}{N-1} \sum_{i=1}^{N} (q_i - \bar{q})^2, \qquad (5) $$
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Variance of the vertical displacement rates (mm² yr⁻²) computed using the three Earth models listed in Table 9 for the Pleistocene (a) and LIA (b) deglaciations.</figcaption></figure>

where $N = 3$ is the number of Earth models, $q$ is either $\dot{u}_v$ or $\dot{g}_v$, $q_i$ is the value of $q$ for the $i$th Earth model and $\bar{q}$ is the mean of the $q_i$ values. To simplify the notation, the dependence of $\chi^2$ from $\eta_a$, $T_a$ and $T_l$ has not been made explicit in eq. (5). For the Pleistocene and LIA deglaciations, we have checked that the pattern of $\chi^2$ is broadly similar for both $\dot{u}$ and $\dot{g}_v$. Thus, we only show results for $\dot{u}_v$ (Fig. 7).

For the Pleistocene deglaciation, Fig. 7(a) shows that $\chi^2 \ge 0.1 \text{ mm}^2 \text{ yr}^{-2}$ for $\dot{u}_v$ for two ranges of $\eta_a$, $T_a$ and $T_l$. For $\dot{g}_v$, these ranges give $\chi^2 \ge 1.7 \times 10^{-3} \mu\text{Gal}^2 \text{yr}^{-2}$. These values correspond to standard deviations of 0.32 mm yr⁻¹ and 0.04 µGal yr⁻¹, respectively. In the first range, for which $\chi^2$ reaches its largest value, $\eta_a$ varies between $6.0 \times 10^{20}$ and $3.2 \times 10^{21}$ Pa s and $T_l$ ($T_a$) from 50 to 75 km (145 to 170 km). These values for $\eta_a$ overlap the upper and lower mantle viscosities and would prevent the asthenosphere to be firmly resolved. The second range, with lower $\chi^2$ values, characterizes low asthenosphere viscosity with $\eta_a \le 10^{19}$ Pa s and thick lithosphere (thin asthenosphere) with $T_l \ge 85$ km ($T_a \le 135$ km). Variations of the predicted geodetic rates at Ny-Ålesund, induced by the Pleistocene deglaciation, due to changes in the mantle viscosities are the largest within these two ranges of $\eta_a$, $T_a$ and $T_l$.

A similar analysis is done computing $\chi^2$ for the geodetic rates predicted for the LIA deglaciation. The variances for $\dot{u}_v$ and $\dot{g}_v$ are about three orders of magnitude lower than those computed for the Pleistocene deglaciation (Fig. 7b). This shows that the predicted geodetic rates due to the Pleistocene deglaciation are more sensitive to the Earth structure than those following the LIA. For both $\dot{u}_v$ and $\dot{g}_v$, $\chi^2$ is the largest for $\eta_a \le 10^{19}$ Pa s, $T_l \le 100$ km and $T_a \ge 120$ km. Comparing Figs 7(a) with (b), we see that the sensitivity of the predicted geodetic rates at Ny-Ålesund to changes in the mantle viscosities depends on the asthenosphere viscosity, the lithosphere/asthenosphere thicknesses, the timescale of the deglaciation and the size of the changing glaciated area.

Viscous responses to surface loads are characterized by their relaxation times which approximately scale with the viscosity of layers. For $T_l$ between 50 and 90 km, the changes affect $\eta_a \ge 6.0 \times 10^{20}$ and $\eta_a \le 10^{19}$ Pa s for the Pleistocene and LIA deglaciations respectively. It results that the viscous response to the LIA deglaciation has a relaxation time more than 60 times shorter than the characteristic relaxation time of the viscous response to the Pleistocene deglaciation.

## 4.3 Observed versus predicted geodetic rates

Here, the computations performed for predicting geodetic rates induced by the Pleistocene and LIA deglaciations at Ny-Ålesund (Section 4.2) are used to determine a set of $\eta_a$, $T_a$ and $T_l$ from the residual rates given in Table 8.

Fig. 8 shows $\dot{u}_v$ and $\dot{g}_v$ induced by past ice-mass changes at Ny-Ålesund and computed using the Earth model EM-A. For thick (resp. thin) lithosphere (resp. asthenosphere) and $\eta_a$ between $10^{20}$ and $10^{21}$ Pa s, $\dot{u}_v \simeq 0.1 \text{ mm yr}^{-1}$ and $\dot{g}_v \simeq 0.03 \mu\text{Gal yr}^{-1}$. However, for thin (resp. thick) lithosphere (resp. asthenosphere) and $\eta_a \le 2.0 \times 10^{18}$ Pa s, $\dot{u}_v$ increases to $6.0 \text{ mm yr}^{-1}$ while $\dot{g}_v$ decreases to $-0.80 \mu\text{Gal yr}^{-1}$. For a fixed asthenosphere viscosity, increasing $T_l$ (resp. decreasing $T_a$) monotonically decreases the magnitude of $\dot{u}_v$ and $\dot{g}_v$. A similar conclusion is obtained when fixing the lithosphere or the asthenosphere thickness and increasing $\eta_a$ from 1.0 to $13.0 \times 10^{18}$ Pa s. The same analysis of the geodetic rates induced by the LIA deglaciation only shows that they are monotonically decreasing over the whole range of the asthenosphere viscosity that we have considered. Sato et al. (2012) obtain similar conclusions for the LIA geodetic effects in Alaska. Thus, the amplitude and behaviour of the geodetic rates as a function of $T_a$, $T_l$ and $\eta_a \ge 3.0 \times 10^{19}$ Pa s are mostly driven by the effects of the Pleistocene deglaciation.

The ranges of $\eta_a$, $T_l$ and $T_a$ for which the geodetic rates induced by the past deglaciations match those determined from the observations (Table 8) are shown by the grey patches in Fig. 8. For the vertical displacement, residual rates are those estimated using the tide gauge observations ($\dot{u}_{tg}$). The matching ranges for residual rates obtained from GPS data ($\dot{u}_{gps}$) are bounded by the blue curves on Figs 8(a), (c) and (e). Depending on the PDIMC model, we obtain several possible ranges. However, there are two common features for all the models. Firstly, ranges determined using vertical displacements from tide gauge are systematically broader than those determined using those from GPS or gravity variations. This is due to the fractional uncertainties which are between 30 and 44 per cent for $\dot{u}_{res}$ estimated from tide gauge data and between 11 and 16 per cent for $\dot{g}_{res}$. GPS-derived residual rates have fractional uncertainties of the order 5–7 per cent. Ranges of $\eta_a$, $T_l$ and $T_a$ estimated from $\dot{u}_{res}$ are similar to those determined from $\dot{g}_{res}$. Secondly, two domains exist for the asthenosphere viscosities where geodetic rates due to past deglaciations can explain the residual rates. For the rates of
<figure><img src="image_11.png" /><figcaption>(a) $\dot{u}_{\text{res}}$ obtained using model A</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) $\dot{u}_{\text{res}}$ obtained using model B</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(e) $\dot{u}_{\text{res}}$ obtained using model C</figcaption></figure>

<caption>(f) $\dot{g}_{\text{res}}$ obtained using model C</caption>

**Figure 8.** Rates of vertical displacements (a, c, e) and gravity variations (b, d, f) at Ny-Ålesund due to the Pleistocene and LIA deglaciations for Earth model EM-A. Grey patches show ranges of $\dot{u}_{\text{res}}$ and $\dot{g}_{\text{res}}$ obtained using $\dot{u}_{\text{tg}}$ and $\dot{g}_{\text{ag/sg}}$ and PDIMC models A (a, b), B (c, d) and C (e, f) (Table 8). Black, green and red curves bound the range of residual rates for predicted rates computed using the models EM-A, EM-B and EM-C, respectively. The blue curves bound the range of residual rates computed using the rate of vertical displacements estimated from GPS data. The dash and dash-dotted curves bound the range of residual rates computed using a LIA deglaciation starting 10 yr earlier and later, respectively.
vertical displacements, which have the lowest fractional uncertainties, the first domain includes viscosities ranging from $2.0 \times 10^{20}$ to $5.0 \times 10^{21}$ Pa s. The second includes viscosities lower than $10^{19}$ Pa s. While $T_l$ can reach 85 km for the first domain, it reaches up to 120 km for the second. The first and second domains correspond to geodetic consequences of the Pleistocene and LIA deglaciations, respectively.

In Fig. 8, ranges of $\eta_a$, $T_l$ and $T_a$ obtained using models EM-B and EM-C are bounded by the green and red curves respectively. Changing the upper and lower mantle viscosities results in $\dot{u}_{res}$ and $\dot{g}_{res}$ that lead to values of $\eta_a$, $T_l$ and $T_a$ which are systematically different from those computed for model EM-A. The differences mostly affect the effects related to the late Pleistocene deglaciation which correspond to asthenosphere viscosities of the same order of magnitude as those of the upper and lower mantle (Fig. 7a).

We restrict ourselves to the ranges of $\eta_a$, $T_l$ and $T_a$ that match $\dot{u}_{res}$ estimated from GPS data. Indeed, these ranges are covered by the rates determined using tide gauge and employing gravity observations broadly leads to the same conclusions. Using $\dot{u}_{res}$ obtained from PDIMC models A and B lead to $\eta_a$ ranging between 2.5 and $3.1 \times 10^{18}$ Pa s for $T_l = 50$ km and $T_a = 170$ km. Increasing $T_l$ (resp. decreasing $T_a$) to 72–79 km (resp. 141–148 km), decreases $\eta_a$ to $10^{18}$ Pa s. The PDIMC model C leads to $\dot{u}_{res}$ that matches $\eta_a$ ranging between 4.4 and $5.3 \times 10^{18}$ Pa s and 2.5 and $12.8 \times 10^{20}$ Pa s for $T_l = 50$ km and $T_a = 170$ km. While the former range decreases to $10^{18}$ Pa s when increasing $T_l$ (resp. decreasing $T_a$) to 91–99 km (resp. 121–129 km), the latter reaches $\eta_a = 10^{21}$ Pa s for $T_l$ between 63 and 70 km and $T_a$ between 150 and 157 km. Starting the LIA deglaciation 10 yr earlier and later does not significantly change the range of $\eta_a$, $T_l$ and $T_a$. The largest changes are obtained by using PDIMC models A and B where the lower limit of $\eta_a$ is reduced by up to 28 per cent ($1.8 \times 10^{18}$ Pa s) and its upper limit increased by up to 23 per cent ($3.7 \times 10^{18}$ Pa s). $T_l$ and $T_a$ are changed by up to 13 per cent. Consequently geodetic rates corrected for effects due to SLC and PDIMC can be explained by geodetic effects due to the Pleistocene and LIA deglaciations provided $\eta_a \le 5.5 \times 10^{18}$ Pa s and $T_l \le 100$ km and $T_a \ge 120$ km. Regarding the asthenosphere viscosity, these results agree with those of Breuer & Wolf (1995) and Kaufmann & Wolf (1996). However, the lithosphere thickness is 10–60 km thinner than theirs. Alternatively, if the PDIMC model C is invoked, $\eta_a$ between 2.5 and $12.8 \times 10^{20}$ Pa s and $T_l$ between 50 and 70 km and $T_a$ between 150 and 170 km would also explain the residual rates. This would rule out the existence of a low-viscosity asthenosphere and would favour the existence of an upper mantle with the traditional viscosity ($10^{21}$ Pa s) inferred from the response to the late Pleistocene deglaciation.

## 5 CONCLUDING REMARKS

We estimated a new rate for the gravity variations in Ny-Ålesund using repeated absolute and continuous relative gravity measurements. This rate is $-1.39 \pm 0.11 \ \mu\text{Gal yr}^{-1}$. It is obtained by removing the significant seasonal contribution as suggested by van Camp et al. (2013). We also updated the rate of vertical displacements from the NYAL and NYA1 GPS time-series provided by the JPL. The average vertical velocity is estimated by removing the seasonal solar and draconitic signals. It is $7.94 \pm 0.21 \ \text{mm yr}^{-1}$, in agreement with the vertical uplift rate estimated using the rate of SLC from the Ny-Ålesund tide gauge data and the rate of absolute sea level variations computed using the ocean general circulation model GLORYS.

We extracted the seasonal signals from the gravity variations and vertical displacements time-series and compared them to the GLDAS/Noah and ERA-interim hydrology models. We found that the gravity seasonal signal strongly agrees with that estimated from the ERA-interim model. The seasonal vertical motions better agree with those computed from the combination of the GLDAS/Noah model and the ocean model TUGO-m. This shows that, in Ny-Ålesund, the local and the non-local hydrology contributions are better represented in the ERA-interim and GLDAS/Noah models, respectively. However, the magnitudes of the vertical displacements computed from the hydrology models are smaller than the observations suggesting potential model-related or instrumental errors.

We also investigated the geodetic contributions of the present-day and past ice-mass changes to the observed geodetic rates. We used published rates of gravity variations and vertical displacements derived from realistic models of ice-mass change in Svalbard. We computed the rates of gravity variations and vertical displacements at Ny-Ålesund induced by the SLC due to the melting of the major world ice-sheets and glaciers. We found that these rates are $-0.29 \pm 0.03 \ \mu\text{Gal yr}^{-1}$ and $1.11 \pm 0.10 \ \text{mm yr}^{-1}$, respectively. They represent 21 and 14 per cent of the observed geodetic rates. They were accounted for when estimating rate residuals.

The geodetic contributions from past ice-mass changes in Svalbard were investigated considering both the Pleistocene and the LIA deglaciations. While the geodetic rates induced by the former at Ny-Ålesund are computed using the ICE-3G model, those induced by the LIA are estimated using a surface disc loads model with simple glaciation and deglaciation histories that are linear in time. We investigated the sensitivity to the Earth structure of the predicted geodetic rates at Ny-Ålesund in response to the past deglaciations. It results that, for a Maxwell rheology, the viscous response to the LIA deglaciation can have relaxation time more than 60 times shorter than that characterizing the viscous response to the Pleistocene deglaciation.

Finally, taking into account the effects due to the PDIMC and the Pleistocene and LIA deglaciations, we inferred from the observed geodetic rates that asthenosphere viscosity ranges between $1.0 \times 10^{18}$ Pa s for lithosphere (resp. asthenosphere) thickness ranging between 50 and 100 km (resp. 120 and 170 km). These conclusions agree with previous studies which used land emergence data (Breuer & Wolf 1995; Kaufmann & Wolf 1996). Our lithosphere thickness is thinner than the one obtained by Kaufmann & Wolf (1996). However, they pointed out that the lithosphere is poorly resolved by their data. The Svalbard archipelago is close to a continental margin. To better infer values of the lithosphere thickness and asthenosphere viscosity, laterally heterogeneous models should be employed (Kaufmann et al. 1997), and possible effects from transient rheology should be taken into account (Spada et al. 2011).

## ACKNOWLEDGEMENTS

AM was supported by an Australian Research Council Super Science Fellowship (FS110200045). GS is partly funded by a DiSBeF (Dipartimento di Scienze di Base e Fondamenti, Urbino University) Research grant. The authors thank O. Omang for providing the SG time-series. Financial and logistical support from the French Polar Institute (IPEV) to perform AG measurements at Ny-Ålesund (program GRAVITE 337) is gratefully acknowledged. O. Omang and another anonymous reviewer are acknowledged for their comments.
REFERENCES

Altamimi, Z., Collilieux, X. & Métivier, L., 2011. ITRF2008: an improved solution of the International Terrestrial Reference Frame, *J. Geod.*, **85**, 457-473.

Amante, C. & Eakins, B.W., 2009. ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis, Technical memorandum nesdis ngdc-24, NOAA, 19 pp.

Boy, J.-P. & Hinderer, J., 2006. Study of the seasonal gravity signal in superconducting gravimeter data, *J. Geodyn.*, **41**(1-3), 227-233.

Breuer, D. & Wolf, D., 1995. Deglacial land emergence and lateral upper-mantle heterogeneity in the Svalbard Archipelago-I. First results for simple load models, *Geophys. J. Int.*, **121**, 775-788.

Carrère, C. & Lyard, F., 2003. Modelling the barotropic response of the global ocean to atmospheric wind and pressure forcing – comparisons with observations, *Geophys. Res. Lett.*, **30**(6), 1275, doi:10.1029/2002GL016473.

Cazenave, A. & Llovel, W., 2010. Contemporary sea level rise, *Annu. Rev. Mar. Sci.*, **2**, 145-173.

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.

Dee, D.P. et al., 2011. The ERA-Interim reanalysis: configuration and performance of the data assimilation system, *Quart. J. Royal Meteo. Soc.*, **137**(656), 553-597.

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.

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

Farrell, W.E. & Clark, J.A., 1976. On postglacial sea level, *Geophys. J. R. astr. Soc.*, **46**, 647-667.

Ferry, N. et al., 2012. GLORYS2V1 global ocean reanalysis of the altimetric era (1993–2009) at meso scale, Mercator Ocean – Quaterly Newsletter 44.

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

Gardner, A.S. et al., 2013. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009, *Science*, **340**, 852-857.

Grove, J.M., 2001. The initiation of the “Little Ice Age” in regions round the North Atlantic, *Climate Change*, **48**(1), 53–82.

Hagedoorn, J.M. & Wolf, D., 2003. Pleistocene and Recent deglaciation in Svalbard: implications for tide-gauge, GPS and VLBI measurements, *J. Geodyn.*, **35**, 415-423.

Hagen, J.O. & Liestøl, O., 1990. Long-term mass-balance investigations in Svalbard, 1950-88, *Ann. Glaciol.*, **14**, 102-106.

Jacob, T., Wahr, J., Pfeffer, W.T. & Swenson, S., 2012. Recent contributions of glaciers and ice caps to sea level rise, *Nature*, **482**, 514-518.

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

Kaufmann, G. & Wolf, D., 1996. Deglacial land emergence and lateral upper-mantle heterogeneity in the Svalbard Archipelago-II. Extended results for high-resolution load models, *Geophys. J. Int.*, **127**, 125-140.

Kaufmann, G., Wu, P. & Wolf, D., 1997. Some effects of lateral heterogeneities in the upper mantle on postglacial land uplift close to continental margins, *Geophys. J. Int.*, **128**, 175-187.

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.

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.

Landvik, J.Y. et al., 1998. The last glacial maximum of Svalbard and the Barents sea area: ice sheet extent and configuration, *Quat. Sci. Rev.*, **17**, 43-75.

Mémin, A., Hinderer, J. & Rogister, Y., 2012. Separation of the geodetic consequences of past and present ice-mass change: influence of topography with application to Svalbard (Norway), *Pure appl. Geophys.*, **169**, 1357-1372.

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.

Mémin, A., Rogister, Y., Hinderer, J., Omang, O.C. & Luck, B., 2011. Secular gravity variations at Svalbard (Norway) from ground observations and GRACE satellite data, *Geophys. J. Int.*, **184**(3), 1119-1130.

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

Milne, G.A. & Mitrovica, J.X., 1998. Postglacial sea-level change on a rotating Earth, *Geophys. J. Int.*, **133**, 1-10.

Moholdt, G., Nuth, C., Hagen, J.O. & Kohler, J., 2010. Recent elevation changes of Svalbard glaciers derived from ICESat laser altimetry, *Remote Sens. Environ.*, **114**(11), 2756-2767.

Nerem, R.S., Chambers, D.P., Choe, C. & Mitchum, G.T., 2010. Estimating mean sea level change from the TOPEX and Jason altimeter missions, *Marine Geodesy*, **33**(1), 435-446.

Nordli, P.Ø. & Kohler, J., 2003. The early 30th century warming. Daily observations at Green Harbour, Grønfjorden, Spitsbergen, DNMI KLIMA 12/03, Det Norske Meteorologiske Institut.

Nuth, C., Kohler, J., Aas, H.F., Brandt, O. & Hagen, J.O., 2007. Glacier geometry and elevation changes on Svalbard (1936-90): a baseline dataset, *Ann. Glaciol.*, **46**, 106-116.

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.

Omang, O.C.D. & Kierulf, H.P., 2011. Past and present-day ice mass variation on Svalbard revealed by superconducting gravimeter and GPS measurements, *Geophys. Res. Lett.*, **38**, L22304, doi:10.1029/2011GL049266.

Peltier, W.R., 1974. The impulse response of a Maxwell Earth, *Rev. Geophys. Space Phys.*, **12**(4), 649-669.

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.

Petit, G. & Luzum, B., 2010. IERS Conventions (2010), IERS Technical Note No. 36, Verlag des Bundesamts für Kartographie und Geodäsie Frankfurt am Main, pp. 179.

Petrov, L. & Boy, J., 2004. Study of the atmospheric pressure loading signal in very long baseline interferometry observations, *J. Geophys. Res.-Solid Earth*, **109**(B3), doi:10.1029/2003JB002500.

Ray, J., Altamimi, Z., Collilieux, X. & van Dam, T., 2008. Anomalous harmonics in the spectra of GPS position estimates, *GPS Solut.*, **12**, 55-64.

Rebischung, P., Griffiths, J., Ray, J., Schmid, R., Collilieux, X. & Garayt, B., 2012. IGS08: the IGS realization of ITRF2008, *GPS Solut.*, **16**, 483-494.

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

Rodell, M. et al., 2004. The global land data assimilation system, *Bull. Am. Meteor. Soc.*, **85**(3), 381-394.

Sato, T., Boy, J.-P., Tamura, Y., Matsumoto, K., Asari, K., Plag, H.-P. & Franci, O., 2006a. Gravity tide and seasonal gravity variation at Ny-Ålesund, Svalbard in Arctic, *J. Geodyn.*, **41**, 234-241.

Sato, T., Okuno, J., Hinderer, J., MacMillan, D.S., Plag, H.P., Francis, O., Falk, R. & Fukuda, Y., 2006b. 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.

Sato, T. et al., 2012. Gravity and uplift rates observed in southeast Alaska and their comparison with GIA model predictions, *J. Geophys. Res.*, **117**, B01401, doi:10.1029/2011JB008485.

Shepherd, A. et al., 2012. A reconciled estimate of ice-sheet mass balance, *Science*, **338**(6111), 1183-1189.
Spada, G., 2003. *The theory behind TABOO*, Samizdat Press edn, Available at: http://samizdat.mines.edu (last accessed 24 April 2014).

Spada, G. & Stocchi, P., 2007. SELEN: a Fortran 90 program for solving the “sea level equation”, *Comput. Geosci.*, **33**(4), 538–562.

Spada, G., Colleoni, F. & Ruggieri, G., 2011. Shallow upper mantle rheology and secular ice sheet fluctuations, *Tectonophysics*, **511**, 89–98.

Spada, G., Melini, D., Galassi, G. & Colleoni, F., 2012a. Modeling sea level changes and geodetic variations by glacial isostasy: the improved SELEN code, preprint, arXiv:1212.5061.

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

Spada, G. et al., 2003. *TABOO User Guide*, Samizdat Press edn, Available at: http://samizdat.mines.edu (last accessed 24 April 2014).

Tanaka, Y., Klemann, V., Martinec, Z. & Riva, R.E.M., 2011. Spectral-finite element approach to viscoelastic relaxation in a spherical compressible Earth: application to GIA modelling, *Geophys. J. Int.*, **184**, 220–234.

Tregoning, P. & Watson, C., 2009. Atmospheric effects and spurious signals in GPS analyses, *J. geophys. Res.*, **114**, B09403, doi:10.1029/2009JB006344.

Tushingham, A. & Peltier, W., 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.

van Camp, M.J., de Viron, O. & Warburton, R.J., 2013. Improving the determination of the gravity rate of change by combining superconducting with absolute gravimeter data, *Comp. Geosci.*, **51**, 49–55.

Van Dam, T.M. & Wahr, J.M., 1987. Displacements of the Earth’s surface due to atmospheric loading: effects on gravity and baseline measurements, *J. geophys. Res.*, **92**(B2), 1281–1286.

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