# Time-varying uplift in Svalbard-an effect of glacial changes

Halfdan Pascal Kierulf, Jack Kohler, Jean-Paul Boy, Emily Geyman, Anthony

Mémin, Ove Omang, Holger Steffen, Rebekka Steffen

► To cite this version:

Halfdan Pascal Kierulf, Jack Kohler, Jean-Paul Boy, Emily Geyman, Anthony Mémin, et al. Time-varying uplift in Svalbard-an effect of glacial changes. Geophysical Journal International, 2022, 231 (3), pp.1518-1534. 10.1093/gji/ggac264. hal-03745732

HAL Id: hal-03745732

https://hal.science/hal-03745732v1

Submitted on 4 Aug 2022

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Time-varying uplift in Svalbard—an effect of glacial changes

Halfdan Pascal Kierulf,<fnref n="1" />,<fnref n="2" /> Jack Kohler,<fnref n="3" /> Jean-Paul Boy,<fnref n="4" /> Emily C. Geyman,<fnref n="3" />,<fnref n="5" />
Anthony Mémin,<fnref n="6" /> Ove C. D. Omang,<fnref n="1" /> Holger Steffen<fnref n="7" /> and Rebekka Steffen<fnref n="7" />

Accepted 2022 July 6. Received 2022 July 2; in original form 2022 April 21

## SUMMARY

We analyse Global Navigation Satellite System (GNSS) data from Svalbard to understand how uplift rates are controlled by the elastic and viscoelastic response of the solid Earth to changes in glacier mass on annual, interannual, decadal, centennial and millennial timescales. To reveal local patterns of deformation, we filter the GNSS time-series with an enhanced common-mode filtering technique where the non-tidal loading signal is incorporated. This technique reduces the estimated uncertainties for 5-yr time-series from 0.8 to 0.3 mm yr⁻<fnref n="1" />. Analysis of the GNSS data with different software—GAMIT, GipsyX, and GINS—produce consistent results that all indicate large temporal variations in uplift. For example, at the Ny-Ålesund GNSS station, uplift varies between 6 and 12 mm yr⁻<fnref n="1" /> for different 5-yr periods, and also shows a significant increase in the last 15 yr. We show that this increase is due to climate change-related ice mass loss in Svalbard. We constrain recent glacier retreat on Svalbard using a series of digital elevation models, and then correct the GNSS-derived uplift records for the elastic signal from these ice mass changes. The residual uplift signal is relatively constant, confirming the hypothesis that current ice mass changes exert a strong influence on GNSS observations. The relatively constant record of residual uplift can be used to constrain other geophysical signals such as the viscoelastic response of the solid Earth to ice loading during the Little Ice Age and the Last Glacial Period. We review uplift results from previous viscoelastic modelling studies and show that the residual signal cannot yet be fully explained. Our new uplift results thus motivate the need for new viscoelastic modelling of the glacial isostatic adjustment process in Svalbard.

**Key words:** Glaciology; Global change from geodesy; Loading of the Earth; Reference systems; Satellite geodesy; Arctic region; Glacial isostatic adjustment; Present-day ice melting.

## 1 INTRODUCTION

Svalbard is an Arctic archipelago located 600 km north of Norway, between 76.5° and 80.5° north and 10° and 34° east (Fig. 1). Svalbard is currently warming at a rate seven times the global average (Nordli et al. 2020), and is experiencing the effect of climate change such as sea level rise, permafrost thaw, sea ice disappearance and glacial retreat (Hanssen-Bauer et al. 2019). Geodetic measurements within an accurate geodetic reference frame offer an important tool to quantify many of the consequences of climate change. On the other hand, these climate-driven changes to the land, ice and ocean challenge the stability of the geodetic reference frame itself. Knowledge about the interaction between geophysical processes, crustal

deformations and reference frame is critical to achieving the Global Geodetic Observing System 2020 (GGOS2020) goal of a reference frame with a stability of 0.1 mm yr⁻<fnref n="1" /> (Plag & Pearlman 2009).

Svalbard's remote location and extreme weather conditions make geodetic measurements challenging. Despite this, the Ny-Ålesund geodetic observatory—equipped with Very Long Baseline Interferometry (VLBI) antennae, Global Navigation Satellite System (GNSS) receivers, gravimeters and Doppler Orbitography and Radiopositioning Integrated by Satellites (DORIS) beacon—has been providing geodetic data since the 1990s. After Altamimi et al. (2007) questioned the stability of geodetic installations in Ny-Ålesund because GNSS-estimated uplift indicated a substantial increase in 2003, an uplift change was introduced in the International

<footnote n="1">¹Geodetic Institute, Norwegian Mapping Authority, 3507, Hønefoss, Norway. E-mail: halfdan.kierulf@kartverket.no</footnote>

<footnote n="2">²Department of Geosciences, University of Oslo, 0371 Oslo, Norway</footnote>

<footnote n="3">³Norwegian Polar Institute, 9296 Tromsø, Norway</footnote>

<footnote n="4">⁴Institut Terre et Environnement Strasbourg, (ITES, UMR7063: Université de Strasbourg, CNRS, ENGEES), 5 Rue René Descartes, 67081 Strasbourg, France</footnote>

<footnote n="5">⁵Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="6">⁶Université Côte d'Azur, CNRS, OCA, IRD, Géoazur, 06560 Valbonne, France</footnote>

<footnote n="7">⁷Lantmäteriet, Geodetisk Infrastruktur, Lantmäterigatan 2c, 802 64 Gävle, Sweden</footnote>
<table><caption>**Table 1.** List of acronyms.</caption><tr><td>ATM</td><td>Atmospheric</td></tr><tr><td>CF</td><td>Centre-of-figure</td></tr><tr><td>CM</td><td>Common mode</td></tr><tr><td>CMB</td><td>Climatic mass balance</td></tr><tr><td>CNES-CLS</td><td>Centre National d'Études Spatiales-Collecte<br>Localization Satellites</td></tr><tr><td>DEM</td><td>Digital elevation model</td></tr><tr><td>DORIS</td><td>Doppler Orbitography and Radiopositioning<br>Integrated by Satellites</td></tr><tr><td>ECMWF</td><td>European Centre for Medium-Range Weather<br>Forecasts</td></tr><tr><td>EOF</td><td>Empirical orthogonal functions</td></tr><tr><td>GCM</td><td>General circulation model</td></tr><tr><td>GGOS</td><td>Global geodetic observing system</td></tr><tr><td>GRGS</td><td>Groupe de Recherche en Géodésie Spatiale</td></tr><tr><td>GIA</td><td>Glacial isostatic adjustment</td></tr><tr><td>GNSS</td><td>Global navigation satellite system</td></tr><tr><td>IGS</td><td>International GNSS service</td></tr><tr><td>IVS</td><td>International VLBI service</td></tr><tr><td>ITRF</td><td>International terrestrial reference frame</td></tr><tr><td>LGP</td><td>Last glacial period</td></tr><tr><td>LIA</td><td>Little Ice Age</td></tr><tr><td>LWS</td><td>Land water storage</td></tr><tr><td>NNR</td><td>No-Net-Rotation</td></tr><tr><td>NPI</td><td>Norwegian Polar Institute</td></tr><tr><td>NTO</td><td>Non-tidal ocean</td></tr><tr><td>PDIM</td><td>Present-day ice melt</td></tr><tr><td>PREM</td><td>Preliminary Reference Earth Model</td></tr><tr><td>RMS</td><td>Root mean square</td></tr><tr><td>RSL</td><td>Relative sea level</td></tr><tr><td>VLBI</td><td>Very long baseline interferometry</td></tr></table>

Terrestrial Reference Frame 2005 (ITRF2005) catalogue. However, Kierulf *et al.* (2009b) showed that this increase in uplift was associated with a few very warm summers that caused negative glacier mass balance. In later ITRFs, the uplift for Ny-Ålesund has been constant.

Overall, the observation data quality in Svalbard is considered excellent, and crustal deformations and gravity variations on Svalbard have been widely studied (e.g. Sato *et al.* 2006a, b; Kierulf *et al.* 2009a, b; Omang & Kierulf 2011; Mémin *et al.* 2011, 2012, 2014; Auriac *et al.* 2016; Rajner 2018; Kierulf *et al.* 2021b). In particular, many papers have focused on the temporal variations in the Ny-Ålesund uplift record. For example, Kierulf *et al.* (2009b) showed that uplift changed from year to year, and that those variations coincide with the elastic response from annual ice mass variations at nearby glaciers (correlation coefficient of 0.80). Omang & Kierulf (2011) observed similar changes in the gravity record, and Rajner (2018) showed that the uplift variations were also recorded by GNSS antennae in Hornsund in the south of Svalbard (HORN, Fig. 1). Moreover, the recent geodetic data from Svalbard contain a measurable viscoelastic response to ice mass changes from the Little Ice Age (LIA, Mémin *et al.* 2014) and the Last Glacial Period (LGP, e.g. Auriac *et al.* 2016).

All Svalbard stations (Fig. 1) experience large seasonal variations in the vertical component (Kierulf *et al.* 2021b), which are explained at a 95 per cent significance level as an elastic response of seasonal ice and snow mass changes at Svalbard. These mass changes were extracted from a climatic mass balance (CMB) model described in van Pelt *et al.* (2019). However, Kierulf *et al.* (2021b) found that this model was not able to explain interannual variations in uplift. The authors suggest that the main reason for this discrepancy is the dynamic behaviour of the Svalbard glaciers (ice flow and calving).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The GNSS network on Svalbard. The location NYAL (Ny-Ålesund) includes the GNSS stations NYAL and NYA1, the VLBI antenna NYALES20, and the superconducting gravimeter. White areas are covered by ice. AF is the Austfonna ice cap discussed in the text.</figcaption></figure>

The observed interannual variations hamper geophysical interpretation of the estimated long-term uplift rates. While Kierulf *et al.* (2009b) showed correlation between the yearly land uplift in Ny-Ålesund and the annual mass balance of a few local glaciers, they were not able to find the total uplift signal from the melting glaciers, only the temporal variations. With the recent models of ice mass changes (Hugonnet *et al.* 2021; Geyman *et al.* 2022), we have new possibilities to model the interannual uplift signal from the ongoing ice mass changes. This signal can then be separated from the (almost) constant geophysical signals from, for example glacial isostatic adjustment (GIA) and plate tectonics. The crucial hypothesis is *H*: 'The interannual variations in uplift can be explained by changes in Present Day Ice Melt (PDIM)'.

The aims of this paper are to: (1) demonstrate the time-dependent variations of uplift in Svalbard, (2) test the hypothesis *H* and (3) separate the time-varying component due to PDIM from the long-term geophysical signal.

We begin by presenting the processing of new vertical and horizontal velocities for the GNSS stations in Svalbard. We also introduce a new glacier surface elevation change model. We use this model to explore the possible uplift variations caused by the elastic response to temporal variations in PDIM. Furthermore, we review previous work that presented velocities (mainly the uplift component) due to the viscoelastic response to the LGP and LIA. Finally, we explore whether the currently observed GNSS velocity field of
Svalbard can be fully explained by the sum of all these elastic and viscoelastic components.

# 2 DATA AND ANALYSIS

## 2.1 GNSS data and analysis

We use data from five permanent GNSS stations in Svalbard (NYAL, NYA1, LYRS, SVES and HORN), and one station on Bear Island (BJOS) 240 km south of Svalbard (Fig. 1). All these stations are located close to existing settlements with power supply and communication infrastructure. We also use data from HAGN station, located at a nunatak in the middle of the glacier Kongsvegen, 30 km southeast of Ny-Ålesund. This station is powered by solar panels and batteries. In the dark season, data are recorded for 24 hr once a week to save power until the sun returns. Data are downloaded once a year in the spring during a glacier mass balance field campaign. Finally, we analyse data from three GNSS stations on the Norwegian mainland (TRO1, VARS and HONS) to help constrain the contributions from regional processes (e.g. GIA).

Data from the GNSS stations at Svalbard are analysed with the program packages GAMIT (Herring et al. 2018), GipsyX (Bertiger et al. 2020) and GINS (Loyer et al. 2012). The GipsyX solution is directly realized in IGb14 (The International GNSS Service (IGS) realization of ITRF2014 (Altamimi et al. 2016)) through the JPL GipsyX No-Net-Rotation (NNR) products. GINS uses the precise point positioning with integer ambiguity resolution mode (IPPP Katsigianni et al. 2019) using MG3 orbit and clock products. MG3 is the IGS REPRO3 product computed by the Centre National d'Études Spatiales-Collecte Localization Satellites (CNES-CLS) team on behalf of the Géodésie Spatiale (GRGS). To ensure a good global realization in IGb14 of the GAMIT solution a global network of approximately 90 global IGS stations is analysed and combined with the Svalbard stations before transforming to IGb14. Daily coordinate time-series are extracted from these solutions. For more details about the GAMIT and GipsyX analysis strategy, see Kierulf et al. (2021b). More details on the MG3 products and the GINS software can be found in Michel et al. (2021).

The time-series are analysed with Hector (Bos et al. 2008) software. We include the trend and annual and semi-annual harmonics in the time-series analysis. None of the stations used in this study had antenna changes during the study period, nor are there any known offsets. We assume that the temporal correlation in the time-series is a combination of white noise and power-law noise, and estimate the spectral index in the time-series analysis.

The time-series for the GAMIT solutions of the stations used in this study are plotted in Fig. 2. The first data point is from January 1995 from the NYAL station. Thus, calculation of PDIM requires glacier elevation data from 1995 onwards.

## 2.2 Glacier elevation data

We constrain recent ice melt on Svalbard from a series of digital elevation models (DEMs) representing Svalbard's ice surface between 1936 and 2018. The DEMs do cover all areas of Svalbard, which requires the usage of earlier time steps than 1995 to obtain a complete set of thickness changes to estimate PDIMs for all GNSS stations. All of the DEMs are constructed through stereo-photogrammetry using aerial photographs or satellite images (Melvær et al. 2014; Noh & Howat 2017; Geyman et al. 2022). If nothing else is explicitly stated, we refer to this ice model as our PDIM-model.

In Section 3.3 we also compare our measured uplift with elastic uplift estimated from the glacier surface elevation change model in Hugonnet et al. (2021). For the seasonal variations in the glacier mass, we use the CMB model in van Pelt et al. (2019). The reader is referred to Kierulf et al. (2021b) for more details about the use of this model.

### 2.2.1 2010

Our reference DEM for glacier change is constructed from aerial photographs acquired during a Norwegian Polar Institute (NPI) mapping campaign conducted from 2008 to 2012, with the bulk of the imagery obtained in the period 2009–2011. DEMs from these campaigns, which were created using the SOCET SET photogrammetry software, are complete for nearly all of Svalbard (Melvær et al. 2014), with only a portion of southern Spitsbergen not yet available. In the latter area, we use elevation data from the ASTER GDEM product, a global compilation of stacked and filtered ASTER DEMs (Nuth et al. 2013) made from images obtained in the period 2000–2013. The data are interpolated onto a 50-m grid (Fig. 3) in the UTM-33N coordinate system, and elevations are relative to the local geoid (ArcGP-2006-sk) defined by the NPI, as are the other DEMs in this paper.

### 2.2.2 1936

The 1936 ice elevation data are from Geyman et al. (2022), and cover nearly the entire archipelago, apart from the eastern half of Austfonna. For the missing area of Austfonna, we simulate elevations using the 2010 DEM minus a difference field derived from $dz_{2010-1936}$, the 1936–2010 elevation changes observed in a buffered area along the missing data edge. The difference field is made using a polynomial fit of $dz_{2010-1936}$ to two predictor variables: (1) the 2010 elevations and (2) the northing coordinate. This simple scheme effectively captures the simultaneous thickening of Austfonna's northern extent as well as the widespread thinning of Austfonna's southern edge and upper elevations during the period 1936–2010.

### 2.2.3 1990

The 1990 ice elevation data are taken from the available archived 1990 NPI DEMs (Melvær et al. 2014), and a 1996 DEM of the Austfonna ice cap made by combining differential SAR interferometry and ICESat laser altimetry (Moholdt & Kääb 2013). In southern Spitsbergen, a hole in the archived 1990 NPI DEM is filled using 1961 and 1970 data, scaled using the assumption that the average mass loss rates from 1936 to 1961/1970 continued to 1990. In the remaining areas with no 1990 data, we add 2/3 of the 1936–2010 difference to the 1936 elevations, reflecting the fact that, elsewhere in Svalbard, roughly 2/3 of the 1936–2010 elevation change occurred in the interval 1936–1990, with the remaining 1/3 occurring from 1990–2010 (Geyman et al. 2022).

### 2.2.4 2013 and 2018

ArcticDEM tiles are constructed from in-track and cross-track high-resolution (~0.5 m) imagery acquired by the DigitalGlobe optical imaging satellites (Porter et al. 2018). Individual DEMs are generated from overlapping pairs of high-resolution optical satellite images (Noh & Howat 2017) to yield 2-m resolution DEM strips.
<table><caption>**Table 2.** GNSS analysis strategies. Elevation dependent site by site functions, based on actual observations are marked with a star (\*). VMF1 is the troposphere mapping functions from Boehm *et al.* (2006). FES2014c is the ocean-tidal loading from Lyard *et al.* (2021).</caption><thead><tr><th></th><th>GAMIT</th><th>GipsyX</th><th>GINS</th></tr></thead><tbody><tr><td>Software version</td><td>10.7</td><td>1.3</td><td>21.2</td></tr><tr><td>Orbit and clock product</td><td>Estimated</td><td>JPL-NNR</td><td>MG3</td></tr><tr><td>Elevation angle cut-off</td><td>10°</td><td>7°</td><td>10°</td></tr><tr><td>Elevation dependent weighting</td><td>a<sup>2</sup> + b<sup>2</sup>/sin(E)<sup>2</sup> (\*)</td><td>1/√sin(E)</td><td>1/[0.15 + 0.85*sin(E)]</td></tr><tr><td>Troposphere mapping function</td><td>VMF1</td><td>VMF1</td><td>VMF1</td></tr><tr><td>Second-order ionosphere model</td><td>IONEX from CODE</td><td>IONEX from JPL</td><td>IONEX from JPL</td></tr><tr><td>Solid Earth tide</td><td>IERS2010</td><td>IERS2010</td><td>IERS2010</td></tr><tr><td>Ocean tidal loading</td><td>FES2004</td><td>FES2004</td><td>FES2014c</td></tr><tr><td>Ocean pole tide model</td><td>IERS2010</td><td>IERS2010</td><td>IERS2010</td></tr><tr><td>Ambiguity</td><td>Resolved</td><td>Resolved</td><td>Resolved</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time-series for the vertical component of the permanent stations in Svalbard. The black curves are the best-fitting model including trend, annual- and semi-annual sinusoidal signals based on the colour-coded time period 2010–2018. This time period is used for the main results of this paper. The location of HORN was shifted by approximately 140 m in 2009, thus the best-fitting curve is not included for the old location. The curves are shifted with respect to each other to improve readability.</figcaption></figure>

We use ArcticDEM Release seven strips (https://www.pgc.umn.edu/data/arcticdem/) which are available for the period 2011–2017, supplemented with a preliminary data set for 2018. Only the 2015 and 2018 strips have good enough coverage to make a DEM of a sufficiently large part of Svalbard.

Strip data are transformed from their Polar Stereographic projection to the UTM-33N coordinate system, their WGS84 ellipsoidal heights converted to a preliminary geoid elevation using the EGM96 geoid model; these data are interpolated onto the reference 50-m grid. The strips then need to be co-registered to remove planar offsets that arise in the photogrammetric processing. In principle, offsets in the horizontal direction could also be removed (Nuth & Kääb 2011), but the method is more suitable for large scenes, with significant land area. Neglecting horizontal offsets does not increase vertical errors significantly since Svalbard glaciers are low-slope; 90 per cent of the glaciated area is less than 10° (Nuth *et al.* 2013).

We remove vertical strip offsets in three steps: First, individual strips are adjusted by subtracting the median of differences to any land points in the 2010 DEM. Secondly, a provisional DEM is obtained by taking the median of all strip data. Third, strips are adjusted again by subtracting the median of differences to the provisional DEM. A new DEM is then made from the median of the twice-adjusted strip data. There still remains a number of holes in

areas with no strip data, and there are artefacts remaining from the original strips. We filter out data points that are significantly different to the 2010 DEM, and areas that are significantly rougher than the 2010 DEM. Some care is needed in choosing these filtering parameters due to ongoing Svalbard glacier surges (Sevestre & Benn 2015), which lead to both large elevation changes and increased roughness. We then calculate $dz_{201X-2010}$, interpolate the missing areas in that field using the method of Wang *et al.* (2012), which interpolates missing elements in a 2-d array using adjacent real-valued elements, add this field to the 2010 DEM, and use these data to fill holes in the provisional 2015 and 2018 DEMs.

### 2.2.5 Error analysis and uncertainty propagation

Since DEM errors tend to be spatially autocorrelated (Nuth *et al.* 2007; Rolstad *et al.* 2009), we compute spatial variograms of the elevation difference ($\Delta z$) between the two DEMs for each year of interest (1936, 1990, 2013 and 2018) and the ~2010 NPI reference DEM (Melvær *et al.* 2014) over low-slope (<20°), ice-free land area. Specifically, we estimate the sill, $s_k$, by fitting a spherical variogram model, $S(d, s_k, r_k)$, where $d$ is the spatial lag and $r_k$ is the range. As in Rolstad *et al.* (2009), we propagate the pixelwise uncertainties in $\Delta z$ into the total uncertainty in the Svalbard-wide mean elevation change $\sigma$ by circular integration of the variogram over the glaciated area $A$:

$$ \sigma^2 = \frac{1}{A} \int_A [s_k - S(d, s_k, r_k)] dA. \quad (1) $$

The final propagated uncertainties are shown in Table 4, 5 and A1.

## 2.3 Elastic loadings

Earth's crust responds to different loadings as an elastic body (see e.g. Petrov & Boy 2004). Tidal loadings such as ocean-tide loadings and the ocean pole tide typically are included in GNSS analyses. In contrast, non-tidal loadings from Atmospheric (ATM), Non-Tidal Ocean (NTO) and Land Water Storage (LWS), are not.

### 2.3.1 Loading from ice and snow

The LWS in Arctic areas is dominated by local signals from ice and snow. Such loadings are not very well captured in global models. Kierulf *et al.* (2021b) showed that the elastic response on annual changes in ice and snow from a CMB model (van Pelt *et al.* 2019)
<figure><img src="image_5.png" /><figcaption>Figure 3. Glacier changes in Svalbard, left-hand panel from 1990 to 2010, right-hand panel from 2013 to 2018. White colour refers to areas with no data or no glacier coverage on land. The negative mass balance at Austfonna is due to a large glacier surge starting in 2012 (Dunse et al. 2015).</figcaption></figure>

is recovered by observations at 95 per cent significance level. However, the interannual signal, which contains signals from glacier dynamics, frontal ablation and calving, is not clearly present in the data (Kierulf et al. 2021b).

Measurements of glacier changes that cause these interannual variations with sufficient temporal and spatial resolution are not available. However, the PDIM-model described in Section 2.2 gives the changes in glacier masses over different time periods. The average glacier-induced elastic uplift from each time period is estimated from this model using the Farrell (1972) approach. The uplift and horizontal crustal deformation due to PDIM for different time periods are plotted in Fig. 4.

The seasonal time-varying signals from ice and snow are estimated from the CMB model of van Pelt et al. (2019). The signal is added to the non-tidal loadings described in Section 2.3.2. Since the long-term signal from ice is included in the PDIM-model described in Section 2.2, the seasonal signal from ice is detrended before adding it to the other non-tidal loadings.

### 2.3.2 Other elastic loadings

We compute surface displacements due to atmospheric, induced oceanic and hydrological loading using a convolution of general general circulation model (GCM) outputs with the appropriate Green's functions Farrell (1972), describing the Earth's response to surface loads (Petrov & Boy 2004; Mémin et al. 2020). We choose a classical spherically stratified and non-rotating elastic and isotropic Earth model, based on Preliminary Reference Earth Model (PREM, Dziewonski & Anderson 1981). Displacements are computed in the centre-of-figure (CF) reference frame (Blewitt 2003).

Atmospheric loading models are computed using surface pressure fields from the ECMWF (European Centre for Medium-Range Weather Forecasts) Reanalysis v5 (ERA5; Hersbach et al. 2020). The ocean response to pressure can be modelled using either the classical Inverted Barometer (IB) assumption (Petrov & Boy 2004),

or the TUGO-m barotropic ocean model forced by ERA5 pressure and winds updated from Carrère & Lyard (2003). As shown by Mémin et al. (2020), the combination of ECMWF and TUGO-m models better reflects the atmospheric and oceanic loading effects over the entire temporal spectrum, i.e from daily to interannual periods. Hence, we have used TUGO-m model in this study. Hydrological loading is computed using continental water storage (snow, soil-moisture and eventually canopy water) using the GLDAS2.1/Noah model (Rodell et al. 2004).

The total seasonal loading signal for station i is:

$$H_{L}^{i} = H_{\text{ATM}}^{i} + H_{\text{NTO}}^{i} + H_{\text{LWS}}^{i} + H_{\text{Snow}}^{i} + \text{Res}(H_{\text{Ice}}^{i}) - H_{\text{LWS,Svalbard}}^{i} \quad (2)$$

$H_{\text{ATM}}$, $H_{\text{NTO}}$ and $H_{\text{LWS}}$ are the loading time-series from ATM, NTO and LWS, respectively. $H_{\text{Snow/Ice}}$ is the time-series due to snow/ice changes from the CMB model. Res means the residual time-series after removal of the linear trend. Since the long-term evolution of ice mass balance is included in the PDIM-model the linear trend is removed from $H_{\text{Ice}}$. To avoid including the loading signal from ice and snow on Svalbard twice (both from the LWS and CMB model), the local LWS signal from Svalbard ($H_{\text{LWS, Svalbard}}$) is computed separately for the LWS model and removed. $H_{\text{LWS, Svalbard}}$ is the LWS signal for Svalbard, but where the LWS load is set to zero outside the Svalbard archipelago ($76^{\circ} < \phi < 81^{\circ}$ northern latitude and $10^{\circ} < \lambda < 34^{\circ}$ eastern longitude).

## 2.4 Removal of Common Mode Signal and elastic loading signal

The time-series after removing the non-tidal loadings for station i is:

$$H_{\text{GNSS,L}}^{i} = H_{\text{GNSS}}^{i} - H_{L}^{i}, \quad (3)$$

where $H_{\text{GNSS}}$ is the GNSS uplift time-series and $H_L$ is the loading signal time-series defined in Section 2.3.2.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Crustal motions due to PDIM in Svalbard, left-hand panel from 1990 to 2010, right-hand panel from 2013 to 2018. Arrows are horizontal rates and colour vertical rates.</figcaption></figure>

Stations in a region can have a similar artificial behaviour, a so-called Common Mode (CM) signal (Wdowinski et al. 1997). Removal of this CM signal can reduce noise in the time-series. The reason for CM signal is not fully understood yet. It could be an effect of the GNSS analysis strategy and/or the strategy for reference frame realization. Dedicated CM filtering, Empirical Orthogonal Functions (EOF) or different types of local reference frame realization are possible tools to remove this signal. All these methods presuppose that stations are exposed to the same CM signal. In Arctic areas, we have limited access to nearby stations. As all stations on Spitsbergen are exposed to similar glacier change signals, using one or several of these stations to remove the CM signal will also remove the glacier signal we wish to study.

The BJOS station is located on Bear Island, a small island 240 km south of Svalbard. We used the BJOS time-series for CM filtering, assuming that the glacial signal there is negligible. The loading signal is removed from both the station BJOS and the other stations before the CM filtering. We get the CM—and load-filtered time-series:

$$H_{\text{CM,L}}^{i} = H_{\text{GNSS,L}}^{i} - \text{Res}(H_{\text{GNSS,L}}^{\text{BJOS}}), \quad (4)$$

where $H_{\text{GNSS,L}}^{\text{BJOS}}$ are the time-series for BJOS after removal of the loading signal and Res are the residuals after removing the linear trends of the time-series. If no loading signal is included we get the ordinary CM filtering, that is $H_{\text{CM}}^{i} = H_{\text{GNSS}}^{i} - \text{Res}(H_{\text{GNSS}}^{\text{BJOS}})$.

## 2.5 Viscoelastic loadings

The velocity field after removal of contributions from the different processes discussed above should provide a signal that is related to long-term geodynamic processes. These include GIA, erosion/sedimentation, large-scale and/or regional tectonics, dynamic topography and perhaps some other unknown processes (see e.g. Kierulf et al. 2021a). However, the major contribution is expected to be from GIA in response to deglaciation in the Svalbard and

Barents Sea area during the LGP and LIA. We focus on this part in this study, and briefly review the literature in the next section.

### 2.5.1 Last glacial period

A detailed investigation of GIA in Svalbard and the Barents Sea region in response to the last glaciation began about three decades ago. Based on compilations of relative sea level (RSL) data, Breuer & Wolf (1995) and Kaufmann & Wolf (1996) described the land emergence of Svalbard, but did not provide vertical rates that can be compared to GNSS-derived velocities. Fjeldskaar (1994) calculated the viscoelastic response for the whole of northern Europe including Svalbard. His favoured GIA model shows a weak positive GIA-related uplift of less than 0.5 mm yr⁻¹ in Svalbard. The best-fitting earth models of these three studies show a wide spread in lithospheric thickness and mantle structure and viscosity parameters, the reason for which is still under discussion, mainly with respect to the existence of a low-viscosity asthenosphere below a thin (30 km thick) lithosphere, as suggested by Fjeldskaar (1994). Breuer & Wolf (1995) further point to lateral heterogeneities in lithospheric thickness and upper-mantle viscosity ranges in the Barents Sea/Svalbard region. Kaufmann & Wu (1998) investigated lateral heterogeneities of the asthenosphere beneath the Barents Sea, and provided uplift velocities of slightly below 0 mm yr⁻¹ (western Svalbard) to 2 mm yr⁻¹ (eastern Svalbard). Importantly, the ice models used in previous studies have all been of coarse spatial resolution and/or consisted of highly simplified geometries (e.g. circular shapes). More recent and complex ice models offer the opportunity to refine our representation of solid Earth loading and the ensuing viscoelastic response.

Steffen & Kaufmann (2005) used a version of the ANU-ICE model (Lambeck et al. 1998) to determine the radial subsurface structure of the Barents Sea and aid in shedding light on the asthenosphere question with the help of RSL data. While fits to a widely used three-layer model showed no preference for an asthenosphere
below a thin lithosphere of 60 km thickness, a multilayer model with a low-viscosity asthenosphere (10¹⁹ Pa s) at 120–200 km depth improved the fit to RSL data. The three-layer model yields uplift values of 0.7–1.9 mm yr⁻¹ for the GNSS stations on Svalbard, while the multilayer model provided much lower uplift values of just 0.1–0.7 mm yr⁻¹.

Mémin et al. (2014) applied the ICE-3G ice model (Tushingham & Peltier 1991) in their GIA analysis. Together with models for PDIM and the LIA (see Section 2.5.2), they could explain the uplift and gravity change at Ny-Ålesund with a thin lithosphere and low-viscosity asthenosphere model in line with earlier findings. The GIA-related uplift at Ny-Ålesund is in the range of 0.6–1.3 mm yr⁻¹ (Mémin et al. 2014).

Auriac et al. (2016) used six different ice models in their investigation of GIA in the Barents Sea region. Best-fitting three-layer GIA models to RSL and GNSS-derived uplift data disagree in their values for lithospheric thickness and mantle viscosities, and the uplift pattern shows distinct differences. GIA-related uplift velocities for Svalbard are in a range of 0–5 mm yr⁻¹ depending on the GIA model, but again with smaller values in the west, mainly at a level of 0–1 mm yr⁻¹.

Simon et al. (2018) investigated the present-day GIA signal in Europe using two different ice models in their data-driven approach. Although their focus was Fennoscandia and the British Isles, the authors provide a velocity field that also covers Svalbard. Velocities of their model D1 confirm previous GIA uplift ranges of ca. 0–2 mm yr⁻¹, with smaller values in the west.

### 2.5.2 Little Ice Age

To explain the geodetic changes observed at Ny-Ålesund, Mémin et al. (2014) had to consider the deformation induced by the LIA. As the LIA history was not precisely known, they used simple assumptions based on previous works (see reference in Mémin et al. 2014). Their LIA model has a uniform linear thickening from 1700 to 1930 to reach a maximum ice thickness of 23.8 m during the glaciation phase and a uniform linear thinning from 1930 to 2000 during the deglaciation phase, assumed to be total. The LIA-related uplift that explained the geodetic observations at Ny-Ålesund is in the range of 0.2–4.0 mm yr⁻¹ depending on the lithosphere thickness and the upper-mantle viscosities.

A LIA study for Greenland by Adhikari et al. (2021) showed that a lower-mantle viscosity improved the misfit between models and GNSS observations. In addition, the rheological behaviour of the Earth to loads points to an increase in mantle viscosity and lithospheric thickness for increasing loading times (Lau et al. 2021). Thus, different lower-mantle viscosities and lithospheric thicknesses than those obtained from LGP modelling can be expected for LIA results.

We should further note that the material parameters used (shear modulus) and found (asthenosphere viscosity) by Mémin et al. (2014) point to low Maxwell times of 1–2 yr. This would mean that the viscous behaviour of the asthenosphere would clearly affect surface deformations on decadal scales like those due to PDIM investigated in our study. Such an investigations with viscoelastic PDIM modelling was performed, for example, by Nield et al. (2014) for the Antarctic Peninsula, and confirmed a strong viscoelastic behaviour with GNSS measurements. We do not pursue viscoelastic modelling for PDIM here as the viscosity results in Mémin et al. (2014) were retrieved with a simple LIA history model and should be confirmed first, once available, with a sophisticated ice history

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time-series of uplift trends in Ny-Ålesund after application of a 5-yr moving window. The black horizontal line marks the final uplift after subtracting the effect of PDIM (Table 4).</figcaption></figure>

model. For example, a one magnitude larger asthenosphere viscosity would already make an viscous behaviour due to PDIM less likely as Maxwell time would be at least 14 yr.

## 3 RESULTS

### 3.1 Uplift variations in Ny-Ålesund

Fig. 5 shows the uplift rates in Ny-Ålesund after application of a 5-yr moving window to the observations. The uplift varies between 6 mm yr⁻¹ and a bit more than 11 mm yr⁻¹. Except for some anomalous high uplift around 2004 (see Section 1), there is a gradual increase in the uplift (acceleration) since the start of GNSS observations in 1993.

A consequence of these variations is that we must use equal time intervals to compare uplift for different stations. To make them geophysically applicable the time interval must be the same as in the PDIM-model. Fig. 2 illustrates how the uplift rate based on data from 2010 to 2018 does not agree with the data from the complete time span for the stations NYAL and NYA1.

### 3.2 Effect of removing loading and CM signal

The limitation of time-series length increases the uncertainty of the estimated rates. Kierulf et al. (2021b) showed that a careful treatment of the CM signal is necessary to reveal the seasonal signal. Further, they found that removal of CM and loading signals reduced the noise and seasonal signal in the time-series. Hence, removal of these signals is necessary when optimal accuracy for relatively short time-series is desired. In Fig. 6 the average estimated uncertainty and annual amplitude are plotted for the combinations of CM filtering and load removal described in Section 2.4. Time-series processed with GAMIT, GipsyX and GINS were analysed. We see a reduction both in uncertainties and annual amplitude after removing of elastic loadings and the CM signal.

Tables A3, A4 and A5 in the Appendix list more details about the noise characteristics, such as the spectral index for power-law noise or the amount of white and correlated noise for all stations. The GAMIT time-series with no filtering ($H_{GNSS}$) have an estimated average uncertainty of 0.7 mm yr⁻¹ (0.6–0.8 mm yr⁻¹), while the CM filtered time-series including removal of load signal ($H_{CM,L}$)
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Quality parameters for time-series for the different software and filtering techniques. The coloured boxes are the average uncertainties and the black horizontal lines are the average amplitude of the annual signal.</figcaption></figure>

reduce the estimated uncertainty to on average 0.3 mm yr⁻¹ (0.2–0.4 mm yr⁻¹). This tremendous reduction in uncertainty is due to a reduction in the general noise level, but more important, the remaining noise has less temporal correlation. The amount of power-law noise ($\sigma_{PL}$) is reduced in the time-series and the spectral index ($S_i$) is closer to zero, that is the noise is whiter.

The largest effect on the estimated noise is from the CM filtering. However, the combination of CM filtering and load removal reduces the amplitude of the annual signal from on average 5.2 mm to 1.6mm.

For the GipsyX and GINS software we see similar differences. However, for both solutions and especially the GINS solution the unfiltered time-series ($H_{GNSS}$) have a larger uncertainty. After filtering all the time-series are more or less at the same level for all solutions. Note the large improvements for the GINS solution after filtering of elastic loadings ($H_{GNSS, L}$ versus $H_{GNSS}$). This indicates that the larger uncertainties in the $H_{GNSS}$ reflect the true loading signal.

By comparing different GNSS analysis strategies and loading models Kierulf et al. (2021b) concluded that there is ‘an analysis strategy-dependent effect in the periodic signal’. This effect masks the loading signal in the unfiltered time-series. Hence, we will use the uplift rates from CM filtered time-series, including removal of the loading signal, in the rest of this paper.

## 3.3 Uplift in Svalbard filtered for PDIM

Uplift values for PDIM in the different periods can be found in Table A1. The measured uplift for the different software after removal of PDIM can be found in Table A2. In Fig. 7 the differences between measured uplift and the uplift after correction of the elastic response due to PDIM are plotted. The results are differenced relative to the 2010–2018 rate in Table 4. In Fig. 7(a) the residuals are not corrected for PDIM. If there were no time variations in the uplift, these should not significantly differ from zero. We see a clear increase during the 2010s in this figure. In Fig. 7(b) the signals from our PDIM-model, are removed. We see that both the GAMIT solution and the GINS solution are not significantly different from zero in any time periods. This confirms our hypothesis *H*, that the ongoing changes in glacier mass balance in Svalbard explain the observed varying interannual land uplift. The GipsyX solution for 2015–2018 (see Table A2) shows a larger increase in

uplift than can be explained by the PDIM, but for the other time periods the differences are within the uncertainty level. The reason for the discrepancy for GipsyX for this particular time period is unknown.

In Fig. 7(c) the signal from the PDIM-model in Hugonnet et al. (2021) is removed from the GNSS-uplift. This model provides consecutive ice mass changes between 2000, 2005, 2010, 2015 and 2020. For the southernmost stations HORN, SVES and LYRS the periods 2010–2015 and 2015–2010 agree. For the northernmost stations NYAL, NYA1 and HORN the discrepancy between GNSS and PDIM uplift are larger in 2015–2020. This indicates that this glacier model underestimates the increased mass loss after 2015 in this area. The 2000–2010 values for the two Ny-Ålesund stations agree with the later periods.

The agreement between the uplift solutions from the three different GNSS software and the two PDIM-models is listed in Table 3. For each pair of GNSS solution and PDIM-model and for each time period we computed the root mean square (RMS) of the difference between measured uplift and the uplift due to PDIM. We do not have complete overlap in time of our PDIM-model and the PDIM-model of Hugonnet et al. (2021). To have periods as equal as possible we have used 1990–2010, 2010–2015 and 2013–2018 for the first model and 2000–2010, 2010–2015 and 2015–2020 for the second model.

The RMS of the differences between measured uplift and the uplift estimated from our PDIM-model, are around 0.4 mm yr⁻¹ (0.38–0.44 mm yr⁻¹) for all software. The RMS are around 0.7 mm yr⁻¹ (0.63–0.74 mm yr⁻¹) without PDIM-model applied. The inclusion of uplift from PDIM reduces the discrepancies between different periods by approximately 40 per cent. The results using the PDIM-model from Hugonnet et al. (2021) are larger and vary from 0.54 to 0.91 mm yr⁻¹. This is mainly because the large estimated uplift after 2015 is not reflected in a more negative mass balance. For this PDIM-model also the slightly different time periods used can affect the results.

We note that the PDIM from Hugonnet et al. (2021) has in general a less negative mass balance, thus the uplift for this model is lower. For the period 2010–2015, which is common for both models, the uplift in Ny-Ålesund is 0.6 mm yr⁻¹ lower and in Hornsund 1.4 mm yr⁻¹ lower than that corrected with our PDIM-model.

## 3.4 Final results

The uplift from the GAMIT solution and the PDIM-model for the 2010–2018 period is listed in Table 4 and plotted in Fig. A1. The horizontal velocities are included in Table 5. The period 2010–2018 is chosen because this is the longest period for which we have data from all the stations in the network coincident with an epoch of the PDIM-model. According to Table 3 all the GNSS solutions coincide with the variation present in the PDIM-model at the same level. The GipsyX solution gives too large uplift compared to the PDIM-model for 2015–2018. We therefore recommend the GAMIT or GINS solution for further studies. The differences between the GAMIT and the GINS solution is 0.1 mm yr⁻¹ for most stations and does not exceed 1σ (Table A2). We therefore use the GAMIT solution in our analysis.

Note that while both the GNSS velocities and the elastic response to PDIM are varying with time, the differences are not significantly different from constant.
<figure><img src="image_9.png" /></figure>

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

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

**Figure 7.** Residual trends for different time periods. In panel (a) and (b) the values are relative to the 2010–2018 results (Table 4). In panel (c) they are relative to 2010–2020 results. The trends are estimated from *H*<sub>CM,L</sub> (eq. 4) minus the elastic response from PDIM for the same time interval. The colour-coded rectangles represent the 2σ-error bars. The uncertainties from the GNSS time-series and the PDIM-models are added in quadrature. In (a), (b) and (c) the panel-rows are for GAMIT, GipsyX and GINS respectively. In (a) the GNSS velocities are not corrected for PDIM. In (b) the GNSS velocities are relative to our PDIM-model. In (c) the GNSS velocities are relative to the uplift based on Hugonnet *et al.* (2021). For (a) and (b) the periods 1990–2010, 2010–2015 and 2013–2018 are used. In (c) the periods 2000–2010, 2010–2015 and 2015–2020 are used. The first time-periods for NYAL and NYA1 (until 2010, green rectangles) are not filtered for CM and the time intervals are not fully identical to the ones of the different software and the glacier model.
<table><caption>**Table 3.** Agreement of GNSS-derived uplift and PDIM-caused uplift for different periods. The values are the RMS of the difference between measured uplift and uplift estimated from the PDIM-models. For our PDIM-model, we used the time intervals 1990–2010, 2010–2015 and 2013–2018. The same time intervals were applied for the results without PDIM model. In PDIM-Hugonnet [the model using Hugonnet *et al.* (2021)], we used the time intervals 2000–2010, 2010–2015 and 2015–2020. Unit in mm yr⁻¹.</caption><thead><tr><th></th><th>GAMIT</th><th>GipsyX</th><th>GINS</th></tr></thead><tbody><tr><td>No PDIM-model</td><td>0.74</td><td>0.70</td><td>0.63</td></tr><tr><td>Our PDIM-model</td><td>0.43</td><td>0.44</td><td>0.38</td></tr><tr><td>PDIM-Hugonnet</td><td>0.54</td><td>0.91</td><td>0.65</td></tr></tbody></table>

# 4 DISCUSSION

## 4.1 Vertical movements

As stated, land uplift in Ny-Ålesund shows temporal variations, increasing gradually since the late 1990s, with a short period of accelerated uplift around 2004 (Fig. 5). Also for the other Svalbard stations establish around 2010, we see an increase in the uplift.

The uplift after removing the elastic effect of PDIM (Table 4 and for Ny-Ålesund the black line in Fig. 5) is constant through the whole period for each station (due to our hypothesis *H*, that the ongoing changes in glacier mass balance explains the varying interannual land uplift), but varies with location. This line is close to a lower bound for all rates in the figure. The mass balance for the large glacier Kongsvegen 15 km southwest of Ny-Ålesund has had mostly negative mass balance after 1998, except for 2000 and the period 2006–2010 (see https://www.mosj.no/no/klima/land/massebalanse-isbreer.html). For this period the total uplift is slightly below the uplift after removal of the PDIM. In the rest of the period the measured uplift is higher, that is the GNSS measurements support the glaciological findings that Svalbard glaciers are losing mass after 2000.

Geyman *et al.* (2022) predicted an average ice mass loss of –0.73 m yr⁻¹ in Svalbard in 2100 assuming intermediate CO₂ emission scenario RCP4.5. With the extreme scenario RCP8.5, the ice mass loss is predicted to be –0.92 m yr⁻¹. The 2015–2018 ice mass loss (PDIM) in this study is –0.83 m yr⁻¹, a similar ice mass loss rate to that predicted to occur by 2100. In the same period we observe an increased uplift of the GNSS network. This increased uplift requires at least this mass loss to explain the observations. Hence, the geodetic measurements confirm the increased mass loss found by glaciological methods.

## 4.2 Elastic contribution from Greenland

Crustal deformation in Svalbard is also influenced by ice mass variations outside Svalbard, especially Greenland. Kierulf *et al.* (2021a) calculated the elastic uplift signal in Scandinavia and Svalbard due to glacial melting in Greenland and other glaciated areas. The uplift varied between 0.4 and 0.5 mm yr⁻¹ for 2000–2010 increasing to 0.6–0.7 mm yr⁻¹ for 2010–2014. In Ny-Ålesund the vertical uplift signal from Greenland increased by 0.3 mm yr⁻¹ between these two time periods. Shepherd *et al.* (2020) estimated a total ice mass loss of 244 ± 28 Gt yr⁻¹ for Greenland in the period 2012–2017. This results in an uplift in Svalbard of 0.6 ± 0.1 mm yr⁻¹. The estimated ice mass loss in Shepherd *et al.* (2020) is 44 ± 35 Gt yr⁻¹ for 1997–2002 and 174 ± 30 Gt yr⁻¹ for 2002–2007. This gives an average uplift of Svalbard of 0.3 ± 0.1 mm yr⁻¹ from 1997 to 2007. With the ice model from Shepherd *et al.* (2020) the uplift increased with 0.5 mm yr⁻¹ from 1997–2007 to 2012–2017. The time span is

slightly different in these two data sets, but both give an additional increased uplift due to the elastic uplift from ice mass changes on Greenland that could explain most of the small unmodelled increase in the uplift we see in the GNSS time-series in Ny-Ålesund between the period 1990–2010 and later periods, see Fig. 7 and Table A2.

## 4.3 Horizontal movements and baseline length extension

The horizontal deformations are dominated by plate tectonics while the elastic signal is smaller in magnitude when compared to its vertical magnitudes. In Table 5 the measured horizontal velocities, the plate-tectonic velocities using the model in Altamimi *et al.* (2017) and the horizontal signal from PDIM are included for the 2010–2018 period. The elastic signal is also shown in Fig. 4. The main pattern is an outward spreading from the central parts of Svalbard, but masked by local signal in areas with more extreme mass loss, for example from surging. The westernmost stations in Svalbard, Ny-Ålesund and Hornsund, have a PDIM signal of –0.7 and –0.6 mm yr⁻¹ in the east direction, while the more eastern stations have an east-component velocity of around –0.2 mm yr⁻¹. Hornsund in the south is moving south with about 1.0 mm yr⁻¹ while the Ny-Ålesund in the north is moving 0.3 mm yr⁻¹ to north.

We have an expansion of the network (baseline length from HORN to NYA1) of around 1.4 ± 0.1 mm yr⁻¹ measured with GNSS. This agrees well with what is expected from the loading, 1.3 mm yr⁻¹ from 2010 to 2018. The expected expansion from the 1990–2010 loading model is 0.7 ± 0.1 mm yr⁻¹, which means that the expansion has almost doubled now. The observed expansion for the earliest period was 0.2 ± 0.1 mm yr⁻¹ between HORN station and NYA1. However, for HORN we have only data between 2005 and 2010. Our results confirm Rajner (2018), who also observed an acceleration in the baseline length between Hornsund and Ny-Ålesund.

## 4.4 Impact on the reference frame

As noted, a few very warm summers led to increasingly negative mass balance and an increase in GNSS-estimated uplift in 2003 (Kierulf *et al.* 2009b). As a result, ITRF2005 listed an uplift change which later ITRFs have not included. However, with steadily increasing demand for accuracy of the reference frame and given the increasingly extreme climate conditions, we see a need for improved glacier models to model the elastic uplift. And by that provide the needed corrections to ensure the future stability of the geodetic installations in Svalbard, for example for new ITRF realizations.

The challenge with climate-induced deformations in the reference frame was highlighted when ITRF2020 was realized: the scale of VLBI appeared to drift after 2013. The scale and scale rate of ITRF2020 was determined by Satellite Laser Ranging (SLR) solutions from 1997.7 - 2021.0 and VLBI solutions up to 2013.75 (for details see https://itrf.ign.fr/en/solutions/ITRF2020). Recently, the International VLBI Service (IVS)-analysis group was informed about the scale issue and a possible drift in the up-component of the VLBI in Ny-Ålesund (J.M. Gipson, personal communication, April 2022). We have proven that the increased uplift in Ny-Ålesund is explained by variations in PDIM and that improved glacier models are necessary to maximize the benefit of Ny-Ålesund for geodesy and earth observation.
<table><caption>**Table 4.** Uplift rates for the GNSS stations in Svalbard and northern Norway. The rates are based on the GAMIT solution and the time span is 2010–2018. The uplift rates are estimated using eq. (4) except BJOS that was estimated using eq. (3). Res are the residual velocities after subtracting the effect of PDIM. Unit in mm yr⁻¹.</caption><thead><tr><th>Station</th><th>Latitude</th><th>Longitude</th><th>GNSS</th><th>PDIM</th><th>Res</th></tr></thead><tbody><tr><td>NYAL</td><td>78.9296</td><td>11.8651</td><td>9.5 ± 0.1</td><td>2.7 ± 0.2</td><td>6.8 ± 0.2</td></tr><tr><td>NYA1</td><td>78.9296</td><td>11.8653</td><td>9.5 ± 0.2</td><td>2.7 ± 0.2</td><td>6.9 ± 0.3</td></tr><tr><td>LYRS</td><td>78.2288</td><td>15.3973</td><td>7.6 ± 0.2</td><td>2.5 ± 0.2</td><td>5.1 ± 0.3</td></tr><tr><td>HORN</td><td>77.0025</td><td>15.5386</td><td>9.3 ± 0.1</td><td>3.9 ± 0.3</td><td>5.5 ± 0.3</td></tr><tr><td>SVES</td><td>77.8992</td><td>16.7246</td><td>6.3 ± 0.2</td><td>3.3 ± 0.2</td><td>3.0 ± 0.3</td></tr><tr><td>HAGN</td><td>78.8001</td><td>13.1020</td><td>12.7 ± 0.4</td><td>4.1 ± 0.3</td><td>8.5 ± 0.5</td></tr><tr><td>BJOS</td><td>74.5033</td><td>19.0014</td><td>0.0 ± 0.3</td><td>0.3 ± 0.0</td><td>-0.3 ± 0.3</td></tr><tr><td>TRO1</td><td>69.6627</td><td>18.9397</td><td>3.5 ± 0.3</td><td>0.1 ± 0.0</td><td>3.5 ± 0.3</td></tr><tr><td>VARS</td><td>70.3364</td><td>31.0312</td><td>3.3 ± 0.2</td><td>0.1 ± 0.0</td><td>3.2 ± 0.2</td></tr><tr><td>HONS</td><td>70.9771</td><td>25.9649</td><td>2.4 ± 0.1</td><td>0.1 ± 0.0</td><td>2.3 ± 0.1</td></tr></tbody></table>

<table><caption>**Table 5.** Horizontal rates for the GNSS stations in Svalbard and northern Norway. The rates are based on the GAMIT solution and the time span is 2010–2018. The uplift rates are estimated using eq. (3). Plate are the plate-tectonic velocities in ITRF2014 (Altamimi et al. 2017). Res are the residual velocities after subtracting plate motion and the effect of PDIM. Unit in mm yr⁻¹.</caption><thead><tr><th>Station</th><th>GNSS</th><th>Plate</th><th>PDIM</th><th>Res</th></tr></thead><tbody><tr><td colspan="5"><strong>East</strong></td></tr><tr><td>NYAL</td><td>10.0 ± 0.1</td><td>10.4 ± 0.1</td><td>-0.7 ± -0.0</td><td>0.3 ± 0.2</td></tr><tr><td>NYA1</td><td>10.3 ± 0.1</td><td>10.4 ± 0.1</td><td>-0.7 ± -0.0</td><td>0.6 ± 0.1</td></tr><tr><td>LYRS</td><td>12.3 ± 0.1</td><td>11.6 ± 0.1</td><td>-0.2 ± -0.0</td><td>0.9 ± 0.1</td></tr><tr><td>HORN</td><td>11.7 ± 0.1</td><td>12.1 ± 0.1</td><td>-0.6 ± -0.0</td><td>0.2 ± 0.1</td></tr><tr><td>SVES</td><td>15.8 ± 0.5</td><td>12.1 ± 0.1</td><td>-0.3 ± -0.0</td><td>4.0 ± 0.5</td></tr><tr><td>HAGN</td><td>12.8 ± 0.4</td><td>10.8 ± 0.1</td><td>-0.2 ± -0.0</td><td>2.2 ± 0.4</td></tr><tr><td>BJOS</td><td>13.7 ± 0.2</td><td>13.9 ± 0.1</td><td>0.0 ± 0.0</td><td>-0.2 ± 0.3</td></tr><tr><td>TRO1</td><td>14.8 ± 0.1</td><td>15.6 ± 0.1</td><td>0.0 ± 0.0</td><td>-0.8 ± 0.1</td></tr><tr><td>VARS</td><td>17.7 ± 0.1</td><td>18.1 ± 0.1</td><td>0.0 ± 0.0</td><td>-0.4 ± 0.1</td></tr><tr><td>HONS</td><td>16.3 ± 0.1</td><td>16.8 ± 0.1</td><td>0.0 ± 0.0</td><td>-0.5 ± 0.1</td></tr><tr><td colspan="5"><strong>North</strong></td></tr><tr><td>NYAL</td><td>14.8 ± 0.1</td><td>15.5 ± 0.1</td><td>0.3 ± 0.0</td><td>-1.0 ± 0.2</td></tr><tr><td>NYA1</td><td>14.8 ± 0.1</td><td>15.5 ± 0.1</td><td>0.3 ± 0.0</td><td>-0.9 ± 0.1</td></tr><tr><td>LYRS</td><td>14.3 ± 0.2</td><td>15.1 ± 0.1</td><td>-0.0 ± -0.0</td><td>-0.8 ± 0.2</td></tr><tr><td>HORN</td><td>13.2 ± 0.1</td><td>15.1 ± 0.1</td><td>-1.0 ± -0.1</td><td>-0.8 ± 0.1</td></tr><tr><td>SVES</td><td>14.5 ± 0.5</td><td>14.9 ± 0.1</td><td>-0.0 ± -0.0</td><td>-0.4 ± 0.5</td></tr><tr><td>HAGN</td><td>15.0 ± 0.5</td><td>15.3 ± 0.1</td><td>0.1 ± 0.0</td><td>-0.4 ± 0.5</td></tr><tr><td>BJOS</td><td>13.3 ± 0.2</td><td>14.6 ± 0.1</td><td>-0.2 ± -0.0</td><td>-1.1 ± 0.2</td></tr><tr><td>TRO1</td><td>15.1 ± 0.1</td><td>14.6 ± 0.1</td><td>-0.0 ± 0.0</td><td>0.5 ± 0.1</td></tr><tr><td>VARS</td><td>12.5 ± 0.1</td><td>12.7 ± 0.1</td><td>-0.0 ± 0.0</td><td>-0.2 ± 0.2</td></tr><tr><td>HONS</td><td>13.5 ± 0.1</td><td>13.6 ± 0.1</td><td>-0.0 ± 0.0</td><td>-0.0 ± 0.2</td></tr></tbody></table>

## 4.5 Comparison with models of viscoelastic response on LGP and LIA

Mémin et al. (2014) were able to explain multitechnique geodetic rates at Ny-Ålesund with a GIA model that has a thin lithosphere and a low-viscosity asthenosphere (see Section 2.5.2). We compare their modelling results of the 3-D velocity field generated by the LGP and the LIA to our residual (PDIM-corrected and plate tectonics-corrected) velocities. Table 6 lists the results of their best-fitting models including their uncertainty range for the GNSS stations in Svalbard. Given that the observations are corrected for elastic effects, one would expect a good agreement of the modelling result and the corrected observations. However, there is no such agreement nor can one see an apparent consistency or correlation in all three components. Even when considering the uncertainties of the observed and modelled velocities no agreement is achieved.

The SVES station seems to be unstable as its velocities, especially the east component, disagree with the pattern found for other

stations. Although the uplift component fits the modelled results, we assume that it might be 1.5–2.0 mm yr⁻¹ too low.

Considering the uncertainties, and thus the upper limit of the modelled velocities and the lower limit of the observations, the modelled uplift for most stations would still be 0.5–1.0 mm yr⁻¹ too low. The elastic signal from Greenland ice mass loss is around 0.6 mm yr⁻¹ at Svalbard. Without considering uncertainties, differences in the uplift are mainly at a level of 2–3 mm yr⁻¹.

The LGP contribution for all stations is ca. 1 mm yr⁻¹, which agrees with previously published results (see Section 2.5.1). We thus think that this contribution may not change significantly with a new generation of GIA models for the LGP. Hence, we suggest that the modelling of the contribution of the LIA should be improved.

Other geodynamic processes such as erosion and sedimentation may be considered in future as well. However, we note that van der Wal & IJpelaar (2017) showed that such effect might be less than 0.2 mm yr⁻¹ in the uplift component, and therefore will not sufficiently help to explain the velocity field completely. Similarly, effects of dynamic topography in the Svalbard region can also be assumed to be (much) less than 0.2 mm yr⁻¹ (Conrad & Behn 2010).

Finally, we note that improved modelling should not only focus on fitting GNSS velocities. GIA modelling should also consider a good fit to RSL data of the Svalbard archipelago and the Barents Sea region. Additionally, a good fit to gravity observations, both on ground and from space, should be achieved.

## 5 CONCLUSIONS

The Arctic archipelago of Svalbard is experiencing some of the most dramatic modern consequences of climate change. However, the geodetic installations in the Arctic needed to observe these effects of climate change are themselves affected by the climate change phenomena, such as glacier retreat, and therefore the reference frame is potentially degraded. A careful treatment of crustal movements is therefore necessary to sustain a stable and high-quality reference frame over time.

The GNSS stations in Svalbard show large interannual variations, especially in uplift. To ensure the long-term stability of the reference frame, and to understand the processes behind the variations, we pursued three main goals in this study: (1) demonstrate the time-dependent variations of uplift in Svalbard, (2) confirm the hypothesis *H*, that the interannual variations in uplift can be explained by changes in PDIM and (3) separate the time-varying component due to PDIM from the long-term geophysical signal.

We analysed geodetic time-series with three different methods and software, including filtering of CM and removal of
<table><caption>**Table 6.** Estimated surface velocities based on models of Mémin *et al.* (2014) at Svalbard GNSS stations due to glaciation during the LGP and the LIA and their sum, which is compared to the residual GNSS velocities after subtracting plate motion (in the horizontals only) and the effect of PDIM from observed GNSS rates, see Tables 4 and 5 (last columns, respectively). Unit in mm yr⁻¹.</caption><thead><tr><th>Station</th><th>LGP</th><th>LIA</th><th>LGP + LIA</th><th>GNSS-PDIM</th><th>Residuals</th></tr></thead><tbody><tr><td colspan="6"><strong>Up</strong></td></tr><tr><td>NYAL</td><td>1.0 ± 0.4</td><td>2.1 ± 1.9</td><td>3.1 ± 2.2</td><td>6.8 ± 0.2</td><td>3.7 ± 2.2</td></tr><tr><td>NYA1</td><td>1.0 ± 0.4</td><td>2.1 ± 1.9</td><td>3.1 ± 2.2</td><td>6.9 ± 0.3</td><td>3.8 ± 2.2</td></tr><tr><td>LYRS</td><td>0.9 ± 0.3</td><td>2.9 ± 2.5</td><td>3.9 ± 2.7</td><td>5.1 ± 0.3</td><td>1.2 ± 2.7</td></tr><tr><td>HORN</td><td>0.9 ± 0.3</td><td>1.0 ± 1.0</td><td>2.0 ± 1.2</td><td>5.5 ± 0.3</td><td>3.5 ± 1.2</td></tr><tr><td>SVES</td><td>0.9 ± 0.3</td><td>2.7 ± 2.3</td><td>3.6 ± 2.5</td><td>3.0 ± 0.3</td><td>-0.6 ± 2.5</td></tr><tr><td>HAGN</td><td>0.9 ± 0.3</td><td>2.6 ± 2.3</td><td>3.6 ± 2.5</td><td>8.5 ± 0.5</td><td>3.9 ± 2.5</td></tr><tr><td>BJOS</td><td>1.0 ± 0.1</td><td>-0.4 ± 0.5</td><td>0.6 ± 0.5</td><td>-0.3 ± 0.3</td><td>-0.9 ± 0.6</td></tr><tr><td colspan="6"><strong>East</strong></td></tr><tr><td>NYAL</td><td>0.1 ± 0.2</td><td>-0.5 ± 0.4</td><td>-0.4 ± 0.4</td><td>0.3 ± 0.2</td><td>0.7 ± 0.4</td></tr><tr><td>NYA1</td><td>0.1 ± 0.2</td><td>-0.5 ± 0.4</td><td>-0.4 ± 0.4</td><td>0.6 ± 0.1</td><td>1.0 ± 0.4</td></tr><tr><td>LYRS</td><td>0.1 ± 0.2</td><td>-0.3 ± 0.2</td><td>-0.2 ± 0.2</td><td>0.9 ± 0.2</td><td>1.1 ± 0.3</td></tr><tr><td>HORN</td><td>0.1 ± 0.2</td><td>-0.1 ± 0.1</td><td>0.0 ± 0.2</td><td>0.2 ± 0.1</td><td>0.2 ± 0.2</td></tr><tr><td>SVES</td><td>0.1 ± 0.2</td><td>-0.1 ± 0.1</td><td>0.0 ± 0.2</td><td>4.0 ± 0.5</td><td>4.0 ± 0.5</td></tr><tr><td>HAGN</td><td>0.1 ± 0.2</td><td>-0.5 ± 0.3</td><td>-0.4 ± 0.3</td><td>2.2 ± 0.5</td><td>2.6 ± 0.6</td></tr><tr><td>BJOS</td><td>0.0 ± 0.2</td><td>0.0 ± 0.0</td><td>0.0 ± 0.2</td><td>-0.2 ± 0.2</td><td>-0.2 ± 0.3</td></tr><tr><td colspan="6"><strong>North</strong></td></tr><tr><td>NYAL</td><td>0.0 ± 0.1</td><td>0.0 ± 0.0</td><td>0.1 ± 0.1</td><td>-1.0 ± 0.2</td><td>-1.1 ± 0.2</td></tr><tr><td>NYA1</td><td>0.0 ± 0.1</td><td>0.0 ± 0.0</td><td>0.1 ± 0.1</td><td>-0.9 ± 0.1</td><td>-1.0 ± 0.1</td></tr><tr><td>LYRS</td><td>0.0 ± 0.1</td><td>0.3 ± 0.2</td><td>0.3 ± 0.3</td><td>-0.8 ± 0.2</td><td>-1.1 ± 0.4</td></tr><tr><td>HORN</td><td>0.1 ± 0.1</td><td>0.3 ± 0.3</td><td>0.4 ± 0.4</td><td>-0.8 ± 0.1</td><td>-1.2 ± 0.4</td></tr><tr><td>SVES</td><td>0.0 ± 0.1</td><td>0.4 ± 0.3</td><td>0.4 ± 0.3</td><td>-0.4 ± 0.5</td><td>-0.8 ± 0.6</td></tr><tr><td>HAGN</td><td>0.0 ± 0.1</td><td>0.1 ± 0.4</td><td>0.1 ± 0.1</td><td>-0.4 ± 0.5</td><td>-0.5 ± 0.5</td></tr><tr><td>BJOS</td><td>0.0 ± 0.1</td><td>-0.1 ± 0.0</td><td>-0.1 ± 0.1</td><td>-1.1 ± 0.2</td><td>-1.0 ± 0.2</td></tr></tbody></table>

different loadings. The results were compared with the elastic response of PDIM from different periods, using the same time periods for both the loading prediction and time-series analysis. The comparison confirmed our hypothesis *H*, 'the interannual variations in uplift can be explained by changes in PDIM'.

Based on this, we have separated the time-varying component from the underlying long-term geophysical processes. The residual velocities were compared with earlier studies of the viscoelastic response of LGP and LIA, but agreement could not be achieved, nor could we identify a consistent pattern. In the uplift component, there are still differences in most stations of at least 0.5 mm yr⁻¹ when considering both model and observation uncertainties. Without considering uncertainties, differences in the uplift are mainly at a level of 2–3 mm yr⁻¹. We suggest revising currently used GIA models to account for the LGP and LIA, with an aim improving the fit to GNSS, RSL and gravity data. Further, we note that a sophisticated model of the LIA is strongly warranted so that the viscosity of the asthenosphere can be well determined. If low viscosity values in the range of 10¹⁸ Pa s as found by Mémin *et al.* (2014) can be confirmed, then also our PDIM investigation should repeated with viscoelastic modelling instead of a purely elastic one. This is because such low viscosities lead to an early (after a few years only) viscous response of the Earth to PDIM Nield *et al.* (e.g. 2014).

Our results are based on existing PDIM-models. To ensure our availability to maintain a stable geodetic reference frame also in Arctic areas, continuously updated PDIM-models are a necessity. The main achievement in this study is to demonstrate that, with new glacier mass balance models based on satellite remote sensing techniques, it is possible to quantitatively model the time-varying uplift at different GNSS stations.

## ACKNOWLEDGMENTS

We thank the editor Prof Bert Vermeersen, Shfaqat Abbas Khan and one anonymous reviewer for their constructive comments that helped to improve the manuscript. We would like to thank Giorgio Spada for making freely available his TABOO software.

## DATA AVAILABILITY

The GNSS data from Hornsund (HORN) is provided by the Institute of Geophysics of the Polish Academy of Sciences. GNSS time-series are available from the authors. The NTL time-series are available from http://loading.u-strasbg.fr/. Distributed time-series of climatic mass balance were extracted from the data set described in van Pelt *et al.* (2019). The 1990 and 2008–2012 glacier elevation data are publicly available on the NPI website: https://publicdatasets.data.npolar.no/kartdata/S0_Terrengmodell/Delmodell/. The 1936/1938 glacier reconstructions are available on the NPI website (https://www.doi.org/10.21334/npolar.2021.f6afca5c) and on Zenodo (https://doi.org/10.5281/zenodo.5806388).

## REFERENCES

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

Altamimi, Z., Collilieux, X., Legrand, J., Garayt, B. & Boucher, C., 2007. ITRF2005: a new release of the international terrestrial reference frame based on time series of station positions and earth orientation parameters, *J. geophys. Res.*, **112**(B9), doi:10.1029/2007JB004949.
Altamimi, Z., Rebischung, P., Métivier, L. & Collilieux, X., 2016. ITRF2014: a new release of the international terrestrial reference frame modeling nonlinear station motions, *J. geophys. Res.* **121**(8), 6109–6131.

Altamimi, Z., Métivier, L., Rebischung, P., Rouby, H. & Collilieux, X., 2017. ITRF2014 plate motion model, *J. geophys. Int.* **209**(3), 1906–1912.

Auriac, A., Whitehouse, P., Bentley, M., Patton, H., Lloyd, J. & Hubbard, A., 2016. Glacial isostatic adjustment associated with the barents sea ice sheet: a modelling inter-comparison, *Quater. Sci. Rev.* **147**, 122–135.

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

Blewitt, G., 2003. Self-consistency in reference frames, geocenter definition, and surface loading on the solid Earth, *J. geophys. Res.* **108**(B2), doi:10.1029/2002JB002082.

Boehm, J., Werl, B. & Schuh, H., 2006. Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data, *J. geophys. Res.* **111**(B2), doi:10.1029/2005JB003629.

Bos, M., Fernandes, R., Williams, S. & Bastos, L., 2008. Fast error analysis of continuous GPS observations, *J. Geod.*, **82**, 157–166.

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

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

Conrad, C.P. & Behn, M.D., 2010. Constraints on lithosphere net rotation and asthenospheric viscosity from global mantle flow models and seismic anisotropy, *Geochem. Geophys. Geosyst.*, **11**(5), doi:10.1029/2009GC002970.

Dunse, T., Schellenberger, T., Hagen, J.O., Kääb, A., Schuler, T.V. & Reijmer, C.H., 2015. Glacier-surge mechanisms promoted by a hydro-thermodynamic feedback to summer melt, *Cryosphere*, **9**(1), 197–215.

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

Farrell, W.E., 1972. Deformation of the Earth by surface loads., *RGSP, 10*, 761–797.

Fjeldskaar, W., 1994. Viscosity and thickness of the asthenosphere detected from the Fennoscandian uplift, *Earth planet. Sci. Lett.*, **126**, 399–410.

Geyman, E., van Pelt, W., Maloof, A., Aas, H. & Kohler, J., 2022. Historical glacier change on Svalbard predicts doubling of mass loss by 2100, *Nature*, **601**, 374–379.

Hanssen-Bauer, I., Førland, E., Hisdal, H., Mayer, S., Sandø, A. & Sorteberg, A.E., 2019. *Climate in Svalbard 2100 - a knowledge base for climate adaptation*, Tech. rep., NCCS report 1.

Herring, T., King, R., Floyd, M. & McClusky, S., 2018. Introduction to GAMIT/GLOBK release 10.7, Tech. rep, Mass. Instit. of Technol., Cambridge.

Hersbach, H. et al., 2020. The era5 global reanalysis, *Quart. J. R. Meteorol. Soc.* **146**(730), 1999–2049.

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

Katsigianni, G., Perosanz, F., Loyer, S. & Gupta, M., 2019. Galileo millimeter-level kinematic precise point positioning with ambiguity resolution, *Earth, Planets Space*, **71**(76), doi:10.1186/s40623-019-1055-1.

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, *J. geophys. Int.*, **127**(1), 125–140.

Kaufmann, G. & Wu, P., 1998. Lateral asthenospheric viscosity variations and postglacial rebound: a case study for the Barents Sea, *Geophys. Res. Lett.*, **25**(11), 1963–1966.

Kierulf, H.P., Pettersen, B., McMillan, D. & 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. Measuring Surface deformation induced by present-day ice melting in Svalbard, *Geophys. J. Int.*, **179**(1), 1–13.

Kierulf, H.P., Steffen, H., Barletta, V.R., Lidberg, M., Johansson, J., Kristiansen, O. & Tarasov, L., 2021a. A GNSS velocity field

for geophysical applications in Fennoscandia, *J. Geodyn.*, **146**, doi:10.1016/j.jog.2021.101845.

Kierulf, H.P., van Pelt, W.J.J., Petrov, L., Dähnn, M., Kirkvik, A.-S. & Omang, O., 2021b. Seasonal glacier and snow loading in Svalbard recovered from geodetic observations, *J. geophys. Int.*, **229**(1), 408–425.

Lambeck, K., Smither, C. & Ekman, M., 1998. Tests of glacial rebound models for Fennoscandinavia based on instrumented sea- and lake-leve records, *Geophys. J. Int.*, **135**, 375–387.

Lau, H. C.P., Austermann, J., Holtzman, B.K., Havlin, C., Lloyd, A.J., Book, C. & Hopper, E., 2021. Frequency dependent mantle viscoelasticity via the complex viscosity: cases from Antarctica, *J. geophys. Res.*, **126**(11), e2021JB022622, doi:10.1029/2021JB022622.

Loyer, S., Perosanz, F., Mercier, F., Capdeville, H. & Marty, J.-C., 2012. Zero-difference GPS ambiguity resolution at CNES-CLS IGS Analysis Center, *J. Geod.*, **86**(11), 991–1003.

Lyard, F.H., Allain, D.J., Cancet, M., Carrère, L. & Picot, N., 2021. Fes2014 global ocean tide atlas: design and performance, *Ocean Sci.*, **17**(3), 615–649.

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.

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

Mémin, A., Spada, G., Boy, J.-P., Rogister, Y. & Hinderer, J., 2014. Decadal geodetic variations in Ny-Ålesund (Svalbard): role of past and present ice-mass changes, *J. geophys. Int.*, **198**(1), 285–297.

Mémin, A., Boy, J. & Santamaría-Gómez, A., 2020. Correcting GPS measurements for non-tidal loading, *GPS Solut*, **24**(45), doi:10.1007/s10291-020-0959-3.

Michel, A., Santamaría-Gómez, A., Boy, J.-P., Perosanz, F. & Loyer, S., 2021. Analysis of GNSS displacements in europe and their comparison with hydrological loading models, *Remote Sens.*, **13**(22), doi:10.3390/rs13224523.

Moholdt, G. & Kääb, A., 2013. A new DEM of the Austfonna ice cap by combining differential SAR interferometry with ICESat laser altimetry, *Polar Res.*, **31**(1), doi:10.3402/polar.v31i0.18460.

Nield, G.A. et al., 2014. Rapid bedrock uplift in the Antarctic Peninsula explained by viscoelastic response to recent ice unloading, *Earth planet. Sci. Lett.*, **397**, 32–41.

Noh, M.-J. & Howat, I.M., 2017. The surface extraction from TIN based search-space minimization (SETSM) algorithm, *ISPRS J. Photogram. Remote Sens.*, **129**, 55–76.

Nordli, Ø., Wyszyński, P., Gjelten, H., Isaksen, K., Łupikasza, E., Niedźwiedź, T. & Przybylak, R., 2020. Revisiting the extended Svalbard Airport monthly temperature series, and the compiled corresponding daily series 1898–2018, *Polar Res.*, **39**(3619), doi:10.33265/polar.v39.3614.

Melvær Y., Aas, H.F. & Skoglund, A., 2014. Terrengmodell Svalbard (S0 Terrengmodell), Norwegian Polar Institute.

Nuth, C. & Kääb, A., 2011. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change, *Cryosphere*, **5**(1), 271–290.

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

Nuth, C., Kohler, J., König, M., Von Deschwanden, A., Hagen, J., Kääb, A., Moholdt, G. & Pettersson, R., 2013. Decadal changes from a multi-temporal glacier inventory of Svalbard, *Cryosphere*, **7**(5), 1603–1621.

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.

Petrov, L. & Boy, J.-P., 2004. Study of the atmospheric pressure loading signal in VLBI observations, *J. geophys. Res.*, **109**(B3), doi:10.1029/2003JB002500.
Plag, H. & Pearlman, M., 2009. *Global Geodetic Observing System: Meeting the Requirements of a Global Society on a Changing Planet in 2020*, Springer.

Porter, C. et al., 2018. ArcticDEM, https://doi.org/10.7910/DVN/OHHUKH, Harvard Dataverse, V1.

Rajner, M., 2018. Detection of ice mass variation using GNSS measurements at Svalbard, *J. Geodyn.*, **121**, 20–25.

Rodell, M. et al., 2004. The Global Land Data Assimilation System, *Bull. Am. Meteorol. Soc.*, **85**(3), 381–394.

Rolstad, C., Haug, T. & Denby, B., 2009. Spatially integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis: application to the western Svartisen ice cap, Norway, *J. Glaciol.*, **55**(192), 666–680.

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

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

Sevestre, H. & Benn, D.I., 2015. Climatic and geometric controls on the global distribution of surge-type glaciers: implications for a unifying model of surging, *J. Glaciol.*, **61**, 646–662.

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

Simon, K.M., Riva, R.E.M., Kleinherenbrink, M. & Frederikse, T., 2018. The glacial isostatic adjustment signal at present day in northern Europe and

the British Isles estimated from geodetic observations and geophysical models, *Solid Earth*, **9**(3), 777–795.

Steffen, H. & Kaufmann, G., 2005. Glacial isostatic adjustment of Scandinavia and northwestern Europe and the radial viscosity structure of the Earth's mantle, *J. geophys. Int.*, **163**(2), 801–812.

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

van der Wal, W. & IJpelaar, T., 2017. The effect of sediment loading in Fennoscandia and the Barents Sea during the last glacial cycle on glacial isostatic adjustment observations, *Solid Earth*, **8**(5), 955–968.

van Pelt, W. et al., 2019. A long-term dataset of climatic mass balance, snow conditions, and runoff in Svalbard (1957–2018), *Cryosphere*, **13**(9), 2259–2280.

Wang, G., Garcia, D., Liu, Y., De Jeu, R. & Dolman, A.J., 2012. A three-dimensional gap filling method for large geophysical datasets: application to global satellite soil moisture observations, *Environ. Modell. Softw.*, **30**, 139–142.

Wdowinski, S., Bock, Y., Zhang, J., Fang, P. & Genrich, J., 1997. Southern California Permanent GPS Geodetic Array: spatial filtering of daily positions for estimating coseismic and postseismic displacements induced by the 1992 Landers earthquake, *J. geophys. Res.*, **102**(B8), 18 057–18 070.

APPENDIX

<figure><img src="image_12.png" /><figcaption>**Figure A1.** Uplift and horizontal velocities for the GNSS stations in Svalbard. Horizontal velocities are relative to the ITRF2014 plate-motion model (Altamimi et al. 2017). Left-hand panel shows the GNSS observations. Middle panel is the residual velocity after removing the PDIM signal. Right-hand panel is the residual after removing PDIM, LIA and LGP.</figcaption></figure>
<table><caption>**Table A1.** PDIM uplift signal for Svalbard stations for different time intervals. The PDIM signal is described in Section 2.3.1.</caption><thead><tr><th>Station</th><th>1990–2010</th><th>2010–2015</th><th>2013–2018<br>Trend (mm yr<sup>-1</sup>)</th><th>2015–2018</th><th>2010–2018</th></tr></thead><tbody><tr><td>NYAL</td><td>1.8 ± 0.1</td><td>2.4 ± 0.4</td><td>2.8 ± 0.4</td><td>3.1 ± 0.3</td><td>2.7 ± 0.2</td></tr><tr><td>NYA1</td><td>1.8 ± 0.1</td><td>2.4 ± 0.4</td><td>2.8 ± 0.4</td><td>3.1 ± 0.3</td><td>2.7 ± 0.2</td></tr><tr><td>LYRS</td><td>1.6 ± 0.1</td><td>2.2 ± 0.3</td><td>2.5 ± 0.4</td><td>3.0 ± 0.3</td><td>2.5 ± 0.2</td></tr><tr><td>HORN</td><td>2.5 ± 0.2</td><td>3.7 ± 0.6</td><td>4.1 ± 0.7</td><td>4.1 ± 0.4</td><td>3.9 ± 0.3</td></tr><tr><td>SVES</td><td>1.6 ± 0.1</td><td>3.5 ± 0.5</td><td>3.0 ± 0.5</td><td>3.1 ± 0.3</td><td>3.3 ± 0.2</td></tr><tr><td>HAGN</td><td>3.2 ± 0.2</td><td>3.9 ± 0.6</td><td>4.1 ± 0.7</td><td>4.5 ± 0.4</td><td>4.1 ± 0.3</td></tr></tbody></table>

<table><caption>**Table A2.** Uplift for Svalbard stations after removal of PDIM. The uplift is estimated from eq. (4) and subtracted from the elastic response on the PDIM described in Section 2.3.1. Both the GNSS and PDIM results are from the same time interval.</caption><thead><tr><th>Station</th><th>AC</th><th>1990–2010</th><th>2010–2015</th><th>2013–2018<br>Trend (mm yr<sup>-1</sup>)</th><th>2015–2018</th><th>2010–2018</th></tr></thead><tbody><tr><td rowspan="3">NYAL</td><td>GAMIT</td><td>6.1 ± 0.8</td><td>6.7 ± 0.4</td><td>7.2 ± 0.5</td><td>7.5 ± 0.5</td><td>6.8 ± 0.2</td></tr><tr><td>GipsyX</td><td>6.1 ± 0.5</td><td>6.2 ± 0.4</td><td>6.8 ± 0.5</td><td>7.9 ± 0.4</td><td>6.4 ± 0.2</td></tr><tr><td>GINS</td><td>6.7 ± 0.6</td><td>6.9 ± 0.4</td><td>7.1 ± 0.5</td><td>7.3 ± 0.5</td><td>7.0 ± 0.2</td></tr><tr><td rowspan="3">NYA1</td><td>GAMIT</td><td>6.1 ± 0.7</td><td>6.3 ± 0.4</td><td>7.5 ± 0.5</td><td>8.0 ± 0.5</td><td>6.9 ± 0.3</td></tr><tr><td>GipsyX</td><td>5.9 ± 0.5</td><td>6.0 ± 0.4</td><td>7.1 ± 0.5</td><td>8.4 ± 0.4</td><td>6.4 ± 0.2</td></tr><tr><td>GINS</td><td>6.6 ± 0.5</td><td>6.8 ± 0.4</td><td>7.4 ± 0.5</td><td>7.7 ± 0.4</td><td>7.0 ± 0.3</td></tr><tr><td rowspan="3">LYRS</td><td>GAMIT</td><td>—</td><td>5.0 ± 0.4</td><td>5.5 ± 0.5</td><td>5.6 ± 0.5</td><td>5.1 ± 0.3</td></tr><tr><td>GipsyX</td><td>—</td><td>5.2 ± 0.4</td><td>4.8 ± 0.7</td><td>5.9 ± 1.1</td><td>4.6 ± 0.4</td></tr><tr><td>GINS</td><td>—</td><td>4.9 ± 0.4</td><td>6.0 ± 0.5</td><td>5.5 ± 0.5</td><td>5.2 ± 0.2</td></tr><tr><td rowspan="3">HORN</td><td>GAMIT</td><td>—</td><td>5.5 ± 0.6</td><td>5.6 ± 0.7</td><td>6.0 ± 0.5</td><td>5.5 ± 0.3</td></tr><tr><td>GipsyX</td><td>—</td><td>5.0 ± 0.6</td><td>5.0 ± 0.7</td><td>6.1 ± 0.5</td><td>5.0 ± 0.3</td></tr><tr><td>GINS</td><td>—</td><td>5.7 ± 0.6</td><td>5.5 ± 0.7</td><td>5.5 ± 0.6</td><td>5.6 ± 0.3</td></tr><tr><td rowspan="3">SVES</td><td>GAMIT</td><td>—</td><td>2.9 ± 0.6</td><td>3.6 ± 0.6</td><td>3.4 ± 0.6</td><td>3.0 ± 0.3</td></tr><tr><td>GipsyX</td><td>—</td><td>2.4 ± 0.7</td><td>2.4 ± 0.6</td><td>2.9 ± 0.8</td><td>2.1 ± 0.3</td></tr><tr><td>GINS</td><td>—</td><td>3.1 ± 0.6</td><td>3.5 ± 1.0</td><td>2.2 ± 1.4</td><td>3.0 ± 0.4</td></tr><tr><td rowspan="3">HAGN</td><td>GAMIT</td><td>—</td><td>8.3 ± 1.3</td><td>8.6 ± 0.8</td><td>9.2 ± 0.7</td><td>8.5 ± 0.5</td></tr><tr><td>GipsyX</td><td>—</td><td>6.3 ± 1.0</td><td>7.5 ± 0.8</td><td>9.5 ± 0.6</td><td>7.2 ± 0.5</td></tr><tr><td>GINS</td><td>—</td><td>7.4 ± 1.2</td><td>8.0 ± 1.1</td><td>8.6 ± 1.3</td><td>8.1 ± 0.8</td></tr></tbody></table>
<table><caption>**Table A3.** Parameters for GAMIT uplift time-series in Svalbard. TS is different time-series. $H_{\text{GNSS}}$, $H_{\text{CM}}$, $H_{\text{GNSS, L}}$ and $H_{\text{CM, L}}$ are explained in Section 2.4. Amp. is the amplitude of the annual sinusoidal signal, Si is the spectral index, $\sigma_{pL}$ standard deviation for power-law noise and $\sigma_{wh}$ is the standard deviation for white noise.</caption><thead><tr><th>Station</th><th>TS</th><th>Rate<br>(mm yr<sup>-1</sup>)</th><th>Amp.<br>(mm)</th><th>Si</th><th>σ<sub>pL</sub><br>(mm yr<sup>-Si/4</sup>)</th><th>σ<sub>wh</sub><br>(mm)</th></tr></thead><tbody><tr><td rowspan="4">NYAL</td><td>H<sub>GNSS</sub></td><td>9.7 ± 0.6</td><td>5.2 ± 0.6</td><td>-0.8</td><td>11.9</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>10.1 ± 0.7</td><td>3.8 ± 0.6</td><td>-0.8</td><td>11.9</td><td>0.0</td></tr><tr><td>H<sub>CM</sub></td><td>9.7 ± 0.2</td><td>4.1 ± 0.3</td><td>-0.4</td><td>7.3</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>10.0 ± 0.2</td><td>2.0 ± 0.3</td><td>-0.4</td><td>6.5</td><td>0.0</td></tr><tr><td rowspan="4">NYA1</td><td>H<sub>GNSS</sub></td><td>10.1 ± 0.7</td><td>5.5 ± 0.6</td><td>-0.8</td><td>11.7</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>10.5 ± 0.6</td><td>3.9 ± 0.6</td><td>-0.8</td><td>11.2</td><td>0.0</td></tr><tr><td>H<sub>CM</sub></td><td>10.0 ± 0.3</td><td>4.5 ± 0.4</td><td>-0.5</td><td>7.9</td><td>1.4</td></tr><tr><td>H<sub>CM, L</sub></td><td>10.3 ± 0.2</td><td>2.3 ± 0.3</td><td>-0.4</td><td>7.5</td><td>0.0</td></tr><tr><td rowspan="4">LYRS</td><td>H<sub>GNSS</sub></td><td>7.7 ± 0.5</td><td>4.8 ± 0.5</td><td>-0.6</td><td>10.6</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>8.0 ± 0.6</td><td>4.9 ± 0.6</td><td>-0.7</td><td>11.6</td><td>0.0</td></tr><tr><td>H<sub>CM</sub></td><td>7.7 ± 0.3</td><td>3.1 ± 0.4</td><td>-0.7</td><td>7.7</td><td>2.7</td></tr><tr><td>H<sub>CM, L</sub></td><td>7.9 ± 0.3</td><td>1.6 ± 0.3</td><td>-0.6</td><td>6.9</td><td>2.9</td></tr><tr><td rowspan="4">HORN</td><td>H<sub>GNSS</sub></td><td>9.6 ± 0.7</td><td>5.2 ± 0.6</td><td>-0.8</td><td>12.4</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>10.0 ± 0.8</td><td>4.3 ± 0.7</td><td>-0.8</td><td>12.8</td><td>0.0</td></tr><tr><td>H<sub>CM</sub></td><td>9.3 ± 0.3</td><td>3.4 ± 0.3</td><td>-0.6</td><td>6.2</td><td>2.9</td></tr><tr><td>H<sub>CM, L</sub></td><td>9.7 ± 0.2</td><td>1.5 ± 0.2</td><td>-0.4</td><td>5.2</td><td>2.5</td></tr><tr><td rowspan="4">HAGN</td><td>H<sub>GNSS</sub></td><td>11.6 ± 0.8</td><td>5.8 ± 0.9</td><td>-0.9</td><td>11.7</td><td>1.3</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>12.0 ± 0.8</td><td>3.9 ± 0.9</td><td>-0.8</td><td>11.9</td><td>0.5</td></tr><tr><td>H<sub>CM</sub></td><td>12.1 ± 0.6</td><td>4.6 ± 0.6</td><td>-0.9</td><td>7.2</td><td>2.9</td></tr><tr><td>H<sub>CM, L</sub></td><td>12.7 ± 0.4</td><td>1.1 ± 0.5</td><td>-0.8</td><td>6.1</td><td>2.8</td></tr><tr><td rowspan="4">SVES</td><td>H<sub>GNSS</sub></td><td>6.3 ± 0.6</td><td>4.5 ± 0.6</td><td>-0.6</td><td>13.4</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>6.8 ± 0.7</td><td>2.9 ± 0.7</td><td>-0.7</td><td>14.5</td><td>0.0</td></tr><tr><td>H<sub>CM</sub></td><td>6.2 ± 0.3</td><td>4.3 ± 0.4</td><td>-0.4</td><td>10.1</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>6.6 ± 0.3</td><td>1.0 ± 0.4</td><td>-0.4</td><td>9.7</td><td>1.1</td></tr><tr><td rowspan="4">BJOS</td><td>H<sub>GNSS</sub></td><td>-0.7 ± 0.5</td><td>2.9 ± 0.5</td><td>-0.7</td><td>10.0</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>-0.5 ± 0.6</td><td>3.3 ± 0.5</td><td>-0.8</td><td>10.4</td><td>0.0</td></tr></tbody></table>

<table><caption>**Table A4.** Parameters for GipsyX uplift time-series in Svalbard. TS is different time-series. $H_{\text{GNSS}}$, $H_{\text{CM}}$, $H_{\text{GNSS, L}}$ and $H_{\text{CM, L}}$ are explained in the main document. Amp. is the amplitude of the annual sinusoidal signal, Si is the spectral index, $\sigma_{pL}$ is the standard deviation for power-law noise and $\sigma_{wh}$ is the standard deviation for white noise.</caption><thead><tr><th>Station</th><th>TS</th><th>Rate<br>(mm yr<sup>-1</sup>)</th><th>Amp.<br>(mm)</th><th>Si</th><th>σ<sub>pL</sub><br>(mm yr<sup>-Si/4</sup>)</th><th>σ<sub>wh</sub><br>(mm)</th></tr></thead><tbody><tr><td rowspan="4">NYAL</td><td>H<sub>GNSS</sub></td><td>9.3 ± 0.8</td><td>3.6 ± 0.7</td><td>-0.8</td><td>14.3</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>9.7 ± 0.4</td><td>2.0 ± 0.4</td><td>-0.6</td><td>8.9</td><td>1.8</td></tr><tr><td>H<sub>CM</sub></td><td>9.1 ± 0.2</td><td>3.7 ± 0.3</td><td>-0.3</td><td>8.0</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>9.6 ± 0.1</td><td>1.5 ± 0.3</td><td>-0.2</td><td>6.7</td><td>0.0</td></tr><tr><td rowspan="4">NYA1</td><td>H<sub>GNSS</sub></td><td>9.6 ± 0.9</td><td>3.1 ± 0.8</td><td>-0.8</td><td>14.8</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>10.0 ± 0.4</td><td>2.4 ± 0.5</td><td>-0.7</td><td>9.4</td><td>1.9</td></tr><tr><td>H<sub>CM</sub></td><td>9.5 ± 0.2</td><td>3.2 ± 0.3</td><td>-0.5</td><td>7.2</td><td>3.2</td></tr><tr><td>H<sub>CM, L</sub></td><td>9.9 ± 0.2</td><td>1.1 ± 0.3</td><td>-0.3</td><td>6.7</td><td>1.7</td></tr><tr><td rowspan="4">LYRS</td><td>H<sub>GNSS</sub></td><td>6.9 ± 0.7</td><td>2.5 ± 0.8</td><td>-0.6</td><td>17.7</td><td>2.8</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>7.3 ± 0.6</td><td>2.4 ± 0.6</td><td>-0.7</td><td>12.9</td><td>4.9</td></tr><tr><td>H<sub>CM</sub></td><td>6.9 ± 0.6</td><td>3.9 ± 0.7</td><td>-0.7</td><td>13.7</td><td>5.9</td></tr><tr><td>H<sub>CM, L</sub></td><td>7.2 ± 0.6</td><td>2.0 ± 0.6</td><td>-0.7</td><td>13.0</td><td>6.1</td></tr><tr><td rowspan="4">HORN</td><td>H<sub>GNSS</sub></td><td>8.9 ± 0.7</td><td>3.7 ± 0.7</td><td>-0.7</td><td>13.0</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>9.4 ± 0.4</td><td>2.9 ± 0.4</td><td>-0.6</td><td>8.6</td><td>1.8</td></tr><tr><td>H<sub>CM</sub></td><td>8.6 ± 0.2</td><td>3.4 ± 0.3</td><td>-0.3</td><td>8.2</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>9.1 ± 0.2</td><td>1.4 ± 0.3</td><td>-0.2</td><td>6.8</td><td>0.0</td></tr><tr><td rowspan="4">HAGN</td><td>H<sub>GNSS</sub></td><td>11.3 ± 0.8</td><td>4.4 ± 0.9</td><td>-0.9</td><td>11.6</td><td>2.1</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>11.5 ± 0.8</td><td>3.4 ± 0.8</td><td>-1.0</td><td>9.5</td><td>2.8</td></tr><tr><td>H<sub>CM</sub></td><td>11.2 ± 0.6</td><td>3.4 ± 0.8</td><td>-0.7</td><td>10.2</td><td>3.3</td></tr><tr><td>H<sub>CM, L</sub></td><td>11.6 ± 0.4</td><td>1.5 ± 0.6</td><td>-0.6</td><td>8.1</td><td>3.4</td></tr><tr><td rowspan="4">SVES</td><td>H<sub>GNSS</sub></td><td>4.9 ± 0.6</td><td>2.1 ± 0.7</td><td>-0.5</td><td>15.9</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>5.5 ± 0.4</td><td>3.0 ± 0.5</td><td>-0.4</td><td>12.8</td><td>1.8</td></tr><tr><td>H<sub>CM</sub></td><td>4.8 ± 0.4</td><td>2.8 ± 0.5</td><td>-0.4</td><td>12.9</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>5.4 ± 0.4</td><td>1.4 ± 0.5</td><td>-0.4</td><td>12.1</td><td>2.5</td></tr><tr><td rowspan="4">BJOS</td><td>H<sub>GNSS</sub></td><td>-1.0 ± 0.5</td><td>1.9 ± 0.5</td><td>-0.7</td><td>11.6</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>-0.7 ± 0.2</td><td>2.6 ± 0.3</td><td>-0.4</td><td>7.1</td><td>0.0</td></tr></tbody></table>
<table><caption>**Table A5.** Parameters for GINS uplift time-series in Svalbard. TS is different time-series. $H_{GNSS}$, $H_{CM}$, $H_{GNSS, L}$ and $H_{CM, L}$, are explained in the main document. Amp is the amplitude of the annual sinusoidal signal, Si is the spectral index, $\sigma_{pL}$ is the standard deviation for power-law noise and $\sigma_{wh}$ is the standard deviation for white noise.</caption><thead><tr><th>Station</th><th>TS</th><th>Rate<br>(mm yr<sup>-1</sup>)</th><th>Amp.<br>(mm)</th><th>Si.</th><th>σ<sub>pL</sub><br>(mm yr<sup>-5/4</sup>)</th><th>σ<sub>wh</sub><br>(mm)</th></tr></thead><tbody><tr><td rowspan="4">NYAL</td><td>H<sub>GNSS</sub></td><td>9.7 ± 1.2</td><td>4.3 ± 0.9</td><td>-1.0</td><td>16.5</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>10.1 ± 0.5</td><td>0.8 ± 0.4</td><td>-0.7</td><td>10.3</td><td>0.0</td></tr><tr><td>H<sub>CM</sub></td><td>9.5 ± 0.3</td><td>4.0 ± 0.3</td><td>-0.4</td><td>8.3</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>9.9 ± 0.2</td><td>1.8 ± 0.3</td><td>-0.3</td><td>6.8</td><td>0.0</td></tr><tr><td rowspan="4">NYA1</td><td>H<sub>GNSS</sub></td><td>10.1 ± 1.3</td><td>3.7 ± 1.0</td><td>-1.0</td><td>16.9</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>10.4 ± 0.6</td><td>0.8 ± 0.4</td><td>-0.8</td><td>9.9</td><td>1.2</td></tr><tr><td>H<sub>CM</sub></td><td>9.9 ± 0.3</td><td>3.3 ± 0.3</td><td>-0.5</td><td>7.6</td><td>2.1</td></tr><tr><td>H<sub>CM, L</sub></td><td>10.2 ± 0.2</td><td>1.4 ± 0.3</td><td>-0.4</td><td>6.4</td><td>1.6</td></tr><tr><td rowspan="4">LYRS</td><td>H<sub>GNSS</sub></td><td>8.6 ± 0.5</td><td>3.3 ± 0.6</td><td>-0.5</td><td>13.9</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>8.7 ± 0.4</td><td>1.0 ± 0.5</td><td>-0.6</td><td>10.8</td><td>2.3</td></tr><tr><td>H<sub>CM</sub></td><td>8.5 ± 0.3</td><td>3.3 ± 0.5</td><td>-0.4</td><td>10.4</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>8.5 ± 0.2</td><td>0.9 ± 0.3</td><td>-0.3</td><td>8.0</td><td>2.0</td></tr><tr><td rowspan="4">HORN</td><td>H<sub>GNSS</sub></td><td>9.5 ± 1.3</td><td>3.9 ± 0.9</td><td>-1.0</td><td>16.3</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>9.9 ± 0.6</td><td>1.2 ± 0.5</td><td>-0.8</td><td>10.0</td><td>1.0</td></tr><tr><td>H<sub>CM</sub></td><td>9.2 ± 0.3</td><td>3.5 ± 0.3</td><td>-0.5</td><td>8.0</td><td>0.0</td></tr><tr><td>H<sub>CM, L</sub></td><td>9.6 ± 0.2</td><td>1.6 ± 0.3</td><td>-0.3</td><td>6.6</td><td>0.0</td></tr><tr><td rowspan="4">HAGN</td><td>H<sub>GNSS</sub></td><td>12.0 ± 0.9</td><td>3.9 ± 0.9</td><td>-0.8</td><td>13.7</td><td>1.3</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>12.2 ± 0.8</td><td>1.4 ± 0.7</td><td>-0.9</td><td>11.2</td><td>1.7</td></tr><tr><td>H<sub>CM</sub></td><td>11.9 ± 1.1</td><td>3.7 ± 0.9</td><td>-1.1</td><td>11.7</td><td>3.7</td></tr><tr><td>H<sub>CM, L</sub></td><td>12.1 ± 0.8</td><td>1.0 ± 0.5</td><td>-1.0</td><td>9.5</td><td>3.5</td></tr><tr><td rowspan="4">SVES</td><td>H<sub>GNSS</sub></td><td>6.2 ± 1.0</td><td>3.0 ± 0.9</td><td>-0.8</td><td>18.5</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>6.6 ± 1.0</td><td>2.8 ± 0.8</td><td>-0.9</td><td>16.0</td><td>2.9</td></tr><tr><td>H<sub>CM</sub></td><td>6.1 ± 0.7</td><td>2.5 ± 0.7</td><td>-0.7</td><td>13.9</td><td>3.3</td></tr><tr><td>H<sub>CM, L</sub></td><td>6.5 ± 0.8</td><td>1.6 ± 0.7</td><td>-0.8</td><td>14.0</td><td>4.0</td></tr><tr><td rowspan="2">BJOS</td><td>H<sub>GNSS</sub></td><td>-0.7 ± 0.9</td><td>1.0 ± 0.5</td><td>-0.9</td><td>13.8</td><td>0.0</td></tr><tr><td>H<sub>GNSS, L</sub></td><td>-0.5 ± 0.4</td><td>1.3 ± 0.4</td><td>-0.8</td><td>7.3</td><td>2.1</td></tr></tbody></table>