## Summary

The vertical movement of the Earth's surface is the result of a number of internal processes in the solid Earth, tidal forces and mass redistribution in the atmosphere, oceans, terrestrial hydrosphere and cryosphere. Close to ice sheets and glaciers, the changes in the ice loads can induce large vertical motions at intraseasonal to secular timescales. The Global Positioning System (GPS) and Very Long Baseline Interferometry (VLBI) antennas in Ny-Ålesund, Svalbard that started observations in 1991 and 1995, respectively, observe vertical uplift rates on the order of 8 ± 2 mm yr^{−1}, which are considerably larger than those predicted by postglacial rebound (PGR) models (order 2 mm yr^{−1}). The observations also indicate increased uplift rates starting some time in 2000. A local GPS campaign network that has been reoccupied annually since 1998, reveals a tilting away from the neighbouring glaciers. The Svalbard glaciers have been undergoing melting and retreat during the last century, with increased melting since about 2000. We compared the observed vertical motion to the motion predicted by loading models using a detailed ice model with annual time resolution as forcing. The model predictions correlate well with the observations both with respect to the interannual variations and the spatial pattern of long-term trends. The regression coefficients for predicted and observed interannual variations in height is 1.08 ± 0.38, whereas the regression coefficient for the predicted and observed spatial pattern turns out to be 1.26 ± 0.42. Estimates of the predicted secular trend in height due to PGR and present-day melting are on the order of 4.8 ± 0.3 mm yr^{−1} and thus smaller than the observed secular trend in height. This discrepancy between predictions and observations is likely caused by the sum of errors in the secular rates determined from observations (due to technique-dependent large-scale offsets) and incomplete or erroneous models (unaccounted tectonic vertical motion, errors in the ice load history, scale errors in the viscoelastic PGR models and the elastic models for present-day melting).

## Introduction

The Geodetic Observatory in Ny-Ålesund, Svalbard is located at 78.9°N and 11.9°E on the southern coast of Kings Bay (Fig. 1). The Observatory is a geodetic station in the global geodetic network of the Global Geodetic Observing System (GGOS, Plag & Pearlman 2009). The geodetic infrastructure includes a 20-m Very Long Baseline Interferometry (VLBI)-antenna, several Global Positioning System (GPS) and combined GPS/Global Navigation Satellite System (GLONASS) receivers, a tide gauge, a superconducting gravimeter, frequently repeated absolute gravity measurements and a co-located Doppler Orbitography and Radiopositioning Integrated by Satellites (DORIS) station (for a detailed description of the station see, e.g. Plag 1998; Plag *et al.* 2000; Kierulf *et al.* 2009).

The station delivers data to the International GNSS Service (IGS), the International VLBI Service for Geodesy and Astrometry (IVS), the International DORIS Service (IDS), the Global Geodynamics Project (GGP) and Global Sea Level Observing System (GLOSS). The station is one of the fewer than 30 geodetic stations globally that have three or more co-located space geodetic techniques. These stations are crucial for the determination and maintenance of the International Terrestrial Reference Frame (ITRF) made available by the International Earth Rotation and Reference System Service (IERS). Being the northernmost of these stations, the Geodetic Observatory in fact plays a key role in the realization of the ITRF (Plag 2006b). Because of its unique location in the northern Arctic and close to the glaciers and ice caps on Svalbard, the observations are particularly valuable for studies of global change phenomena on the Northern Hemisphere, such as the impact of the warming in the Arctic (e.g. Overland 2006) on ice melting and on changes in sea level (e.g. Sato *et al.* 2006; Kierulf & Kristiansen 2006).

The Geodetic Observatory is in a location only about 150 km away from the Mid-Atlantic ridge, and in the past, the area has undergone significant tectonic activity (e.g. Birkenmajer 1981; Blythe & Kleinspehn 1998). To determine the spatial extent of the footprint of the Geodetic Observatory and the spatial variability of present-day deformation in the Kings Bay area, in 1998 a GPS control network was established extending in east—west and north—south directions approximately 50 km × 30 km (see Fig. 1; Plag *et al.* 2000; Bockmann *et al.* 2002; Kierulf *et al.* 2002). From 1998 to 2007, in total six GPS campaigns have been carried out, each with at least four consecutive complete days of observations.

During the last century the glaciers on Svalbard have gradually decreased (Hagen *et al.* 2003; Nuth 2006). Especially around the beginning of this century, deglaciation seems to have accelerated (Kohler *et al.* 2007). Climate change forcing the deglaciation is obvious, for example, in the variations of Arctic air temperatures (e.g. Polyakov *et al.* 2003; Bengtsson *et al.* 2004; Overland 2006) and changes in the ocean currents (e.g. Polyakov *et al.* 2004).

It can be expected that the unloading associated with the rapid melting of the ice cover induces significant signals in the geodetic time-series. In fact, in Juneau, Alaska, uplift rates of more than 50 mm yr^{−1} that to a large extent are attributed to the elastic response of the Earth's crust to the unloading of nearby melting glaciers are causing practical problems (Kelly *et al.* 2007). Khan *et al.* (2007) report large vertical uplift rates close to the Greenland ice sheet, where similar large melting takes place, and explain these with the elastic response of the crust to glacial unloading. The analyses of the GPS and VLBI observations from Ny-Ålesund also result in uplift rates more than three times those predicted by postglacial rebound (PGR) models (see, e.g. Sato *et al.* 2006). Moreover, the geodetic time-series exhibit a significant nonlinear component with the uplift rates increasing with time. Sato *et al.* (2006) introduce the effect of Present Day Ice Melt (PDIM) as a possible explanation for the large uplift rates observed in Ny-Ålesund. However, they find that explaining the observed absolute gravity changes there requires about 1.5 times the PDIM consistent with the observed vertical uplift rates.

It is our goal to quantify the relation between the observed vertical displacement of the Earth's surface and the ice load changes. In the next section, we first describe the geodetic database for displacements and present the time-series resulting from the analyses of these observations. In Section 3, we review the available information on changes in the glaciers on Svalbard and present an empirical ice-melt model. Section 4 introduces the loading model and presents the predictions of the vertical displacement induced by the present-day ice load changes. In Section 5, we compare these predictions to the observed displacements before we discuss the results in Section 6.

## Geodetic Observations and Analysis

The VLBI antenna at Ny-Ålesund has been in operation since 1994. The VLBI observations were analysed as part of a terrestrial reference frame solution at the Goddard Space Flight Center (Kierulf *et al.* 2009) using the Calc-Solve software package (Ma *et al.* 1990). In this solution, station positions, station velocities and radio source positions are estimated as global parameters from all the data. To generate the Ny-Ålesund position time-series, the same terrestrial reference frame solution was carried out, except that the position of Ny-Ålesund was estimated as a local parameter for each daily experiment session. The VLBI results (Fig. 2) are affected by some disturbances in the time-series from 1998 to 2000 that are not seen in the GPS time-series, and the formal uncertainties for the VLBI series are higher in this period (Kierulf *et al.* 2009). Hence only VLBI results from summer 2000 onward are included in our analysis later.

The two GPS stations at the Geodetic Observatory used in the present study are NYAL and NYA1, which started operation in 1992 and 1997, respectively. The stations have operated continuously with the largest gaps being 36 and 27 days for NYAL and NYA1, respectively. In total, only about 6 per cent and 4 per cent of observation days are missing for NYAL and NYA1, respectively (see, Kierulf *et al.* 2009, for a detailed discussion of the two stations).

The GPS observations of NYAL and NYA1 were analysed with the program package GPS Inferred Poitioning System/Orbit Analysis and Simulation Software (GIPSY) and the program GPS Analysis Software of MIT (GAMIT). GIPSY was used in the Precise Point Positioning (PPP) mode (Zumberge *et al.* 1997) using the satellite orbits and clocks provided by the Jet Propulsion Laboratory (JPL). An elevation cut-off angle of 10 degrees, Niell mapping function, and ocean loading coefficients from http://www.oso.chalmers.se/loading/ were used. No ambiguity resolution was performed. The daily free solutions were transformed into ITRF2005 using the parameters for a seven-parameter Helmert Transformation provided in the so-called *x*-files of the JPL.

We also carried out daily GIPSY—PPP analysis for a regional Arctic network including the IGS-stations CHUR, YELL, REYK, TIXI and TRO1 in addition to NYAL and NYA1. The ITRF positions and velocities for the stations define an Arctic reference frame. The GIPSY—REG solutions were determined through daily 7-parameter transformation of the free solutions to this Arctic reference frame.

The GAMIT solutions are regional solutions (Herring 2003) using an Arctic network of permanent GPS stations which—besides NYAL and NYA1—includes the stations TIXI, TRO1, YELL, WTRZ, THU2, BILI, IRKT, HERS, METS and ALGO. Elevation cut-off angle, tropospheric mapping function, and ocean loading were the same as for the GIPSY solution. We used the IGS-orbits. The solution was combined with the Scrips Orbits and Permanent Array Center (SOPAC) global solutions (the higs* files) and connected to ITRF2005 using the program package Global Kalman-filter VLBI and GPS Analysis program (GLOBK).

The time-series for the GIPSY—PPP and GAMIT-solutions for NYAL and NYA1 are shown in Figs 2 (a—d). The decrease of the noise level with time is visible for all four time-series. The daily differences between NYAL and NYA1 for the PPP and GAMIT solutions are on the order of ± 10 mm with a WRMS of 4.1 and 2.7 mm for the PPP and GAMIT solutions, respectively (Figs 2, g and i). The time-series of the daily differences show very small differential trends (−0.3 mm yr^{−1} for the PPP solutions and 0.1 mm yr^{−1} for the GAMIT solutions) and almost no seasonal signal. Thus, the two stations appear to measure the same vertical surface displacements if analysed with the same software. Daily difference for the time-series of one station analysed with different programs are much larger with the WRMS for NYAL—PPP minus NYAL—GAMIT and NYA1–PPP minus NYA1–GAMIT being 7.4 mm and 6.4 mm, respectively (Figs 2, h and j). The trends for the differences are 3.6 mm yr^{−1} and 3.7 mm yr^{−1}, respectively. Moreover, the differences between PPP and GAMIT solutions show clear seasonal signals. The differences between the GPS solutions and the VLBI solutions have a WRMS of 8.8 mm and 10.3 mm for the GAMIT and PPP solutions, respectively (see, e.g. Fig 2f). The trends of the differences between VLBI and GPS solutions are slightly smaller than the differences of the trends of the original time-series, indicating that the trend differences are partially caused by the higher temporal sampling of the GPS time-series. However, interannual variations in the time-series of differences between the PPP and GAMIT solutions as well as between the GPS and VLBI solutions are small. Thus, the interannual variations appear least affected by the different reference frames implicit in the PPP, GAMIT, and VLBI solutions. The similarity of interannual variations in different GPS solutions and two different techniques indicates that these variations are geophysical signals. Moreover, these signals are not or only to a small extent affected by the regional filtering implicit in the reference frame realization of the different techniques and analyses, which indicates that these signals are dominated by small spatial scales.

We have compared our solutions for NYAL and NYA1 also to time-series provided by other groups. The time-series computed at the Nevada Geodetic Laboratory using a combination of GIPSY—PPP and the new Ambizap algorithm for ambiguity resolution in a global network of more than 2800 stations (Blewitt 2008) show trends and interannual variability very similar to our GIPSY—PPP solution, whereas the seasonal constituents exhibit some differences. These difference arise from different station-motion models and different transformation from free to ITRF solutions. The similarity of the interannual variations emphasize the insensitivity of these variations to the specific analysis and the implicit reference frame determination.

The secular trends for the two GPS sites NYAL and NYA1 are consistent on the 0.2 mm yr^{−1}-level for each of the two solutions (GIPSY—PPP and GAMIT), whereas the trends for each station exhibit differences of ∼3 mm yr^{−1} between the two solutions (Table 1). The result for the VLBI time-series (i.e. a different technique and a different analysis method) is in between those of the two different GPS analysis methods. The results for the annual and semi-annual amplitudes and phases exhibit a larger coherence between the two stations analysed with the same method than between different methods. In a study of GPS sites distributed over Europe, Kierulf *et al.* (2008) found similar inter-solution differences in trends and amplitudes and phases of annual and semi-annual constituents. Each analysis method realizes the geodetic reference frame in a different way resulting in specific characteristics of the inherent spatial filtering. The geophysical seasonal signal is dominated by large spatial scales and thus strongly modified by the regional filtering specific to each analysis method. The inter-solution biases in secular trends are due to a combination of different factors resulting from different technique-dependent effects, differences in reference frame realizations for the individual solutions, and differences in the model approaches. Based on the results discussed here, in the comparison of observations and model predictions (Section 5), we will utilize the high spatial coherence of homogeneous solutions and account for the inter-solution biases.

The observations collected during the GPS campaigns on the control network are analysed with GIPSY with the standard PPP-solution described earlier and a subsequent ambiguity resolution. In addition, a regional GAMIT solution was also computed using a double differencing approach. For both analyses, the elevation cut-off angle was set to 10 degrees. The two solutions demonstrate some systematic differences, but internal consistency of the solutions appears to be at the sub-millimetre level. Solutions with elevation cut-off angles of 7 and 15 degrees show slightly larger scatters, but the differences are still at the sub-millimetre level.

The vertical rates for the campaign stations are determined for displacement time-series relative to NYAL and NYA1, thus reducing the effect of intra-seasonal variations on the campaign station rates. Over the spatial scale of the control network, intra-seasonal signals to a large extent are common to all stations (common modes). By considering the daily differences of the campaign stations to one or more continuous base stations, a simple form of the regional filtering introduced by Wdowinski *et al.* (1997) is carried out that reduces the temporal variations common to the stations. The spatial pattern of the vertical secular trends (Fig. 3) is computed as the mean of the two GAMIT and GIPSY solutions each relative to NYAL and NYA1 (in total, four samples). Uncertainties are based on the differences between the solutions. The results for the 2007 campaign were excluded, because the receiver at KAPM did not record due to a cable pulled out by seagulls. The lack of observations from that station decreased the quality of the network campaign for that year considerably.

The spatial pattern of the vertical secular changes indicates E—W tilting (Fig. 3), with the western part subsiding and the eastern part of the network uplifting. As will be discussed further, this pattern is consistent with the expected unloading from the retreating glaciers located to the east of the network.

In the next section, we will describe the model for the ice load changes that can be used to predict the expected signals in vertical motion for the VLBI and GPS stations at the Observatory and the control network with sufficient temporal resolution to compare not only the vertical secular rates but also the interannual changes.

## Glacier and Melting

For the prediction of surface displacements induced by changes in the ice load, a model of the spatial-temporal changes in the ice mass is required. From continuity, changes in ice mass *dM*/*dt* at any arbitrary point **x** on a glacier can be expressed by

*b*

_{n}is the specific mass balance, that is, the mass flux into or out of the glacier surface at that point, and the second term on the right-hand side is the planar divergence of the ice flux

**q**averaged from the glacier bed to the surface at that point (Paterson 2004). Ice gain or loss at the bed is neglected.

For reasons discussed further, it is often easier to set up a model of the spatial pattern of surface elevation changes Δ*h*(**x**) = *h*_{2}(**x**) −*h*_{1}(**x**), where *h*_{i} are ice thicknesses at two different times *t*_{i}, *i* = 1, 2. The mass changes described by eq. (1) are related to surface elevation changes d*h*/d*t* through

There are two basic approaches to determining the spatial variation of Δ*h* at a regional scale. The first is to measure it directly, for example by taking the difference of two Digital Elevation Model (DEM) of the ice surface made at different dates (e.g. Etzelmüller 2000). In this case, the spatial distribution of Δ*h* is readily obtained and this is indeed the preferred method. In polar areas, such as Svalbard, however, there are not always sufficient data to generate two DEMs of comparable quality. For the Ny-Ålesund area, the most up-to-date topographic map of the region is based on aerial photographs taken in 1966. There are earlier maps, but these are based on oblique terrestrial and aerial photography from the early 20th century and are not of sufficient geodetic quality. There is also aerial photography coverage from 1990, but so far photogrammetric analyses have only been performed over limited areas. However, to a certain extent it is possible to extrapolate to a larger area differences of surface elevation measurements that are spatially more limited. For Svalbard, examples are data obtained with repeated GPS measurements along single or multiple profiles (Hagen *et al.* 2005; Nuth 2006), or the comparison of contemporary GPS profiles to an older DEM. Because glaciers generally correlate to elevation, it is possible to average elevation change measurements within a particular elevation band to obtain Δ*h*(*z*), where *z* is elevation, which then can be applied to a DEM to specify the spatial distribution of loading. Extrapolating Δ*h* measured along profiles on one glacier to neighbouring glaciers or to a larger region can potentially lead to erroneous estimates because ice surface elevation changes dependent not only on local climate but also on glacier geometry and dynamic regime. However, with sufficient spatial sampling of Δ*h* this approach can give reasonable results. Alternately, a spatially averaged value of Δ*h* can be derived for an entire glacier or region utilizing the glacier's hypsometry, that is, its distribution of area *A* as a function of elevation *z*, using

*A*

_{tot}is the total glacier area and the bar implies spatially averaging.

The second approach to determining the spatial distribution of Δ*h* at regional scales is to integrate eq. (2) twice, first over the glacier surface area to eliminate the flux divergence term (where extra loss through calving can be accommodated in an *a*d hoc way by using a larger value of *b*_{n}(**x**) at the appropriate location, so that eq. (1) still holds), and again over time to yield

*B*

_{n}is an areally averaged mass balance (Paterson 2004) defined by Similarly to eq. (3), where the assumption is that Δ

*h*(

*z*) is a simple function of elevation

*z*, for eq. (4) the specific balance

*b*

_{n}is assumed to be a simple function of elevation, but in the case of mass balance, this assumption is better grounded in theory and practice (Oerlemans 2001). Values of

*b*

_{n}cannot simply be applied to a neighbouring glacier, however, because hypsometries can differ from glacier to glacier. A regional specific mass balance curve

*b*

_{n}(

*z*) always has to be applied to the hypsometry of the glacier or region being modelled.

Both approaches described by eqs (3) and (5), respectively, yield a glacier- or region-wide average. The disadvantage of these approaches is that fine-scale loading cannot be modelled. For example, a retreating glacier with enough mass gain in its upper part to compensate for the frontal loss would have an average ice unloading pattern of zero, yet would experience upward deformation at its front. However, as long as the relative differences in distance from the geodetic stations to the lower versus higher altitude of the glacier are small, the spatial variability of the load changes over the glaciers have only a small impact on the predicted uplift.

A further point is that the available observations of mass balance and surface height changes do not have sufficient temporal resolution for an ice load model and the computation of loading predictions with weekly or monthly resolution that could be compared directly to the observed time-series of displacements. However, eq. (5) can be used to derive seasonal mass balance averages because *b*_{n}(*z*) can be broken down into a winter and summer component, *b*_{w}(*z*) and *b*_{s}(*z*), respectively. Here we will focus on year-to-year changes because these appear to be least affected in the observations of surface displacements by technique and/or analysis-dependent effects (see previous section).

For the surface displacements at Ny-Ålesund induced by ice load changes, the contribution of the glaciers around the observatory is most dominant. Therefore, we separate Svalbard into several main regions (Fig. 4). First, we concentrate on the load model for the near-field, that is, Region 1.

Measurements of elevation changes in Region 1 are restricted to a few airborne and ground-based GPS profiles on two of the larger glacier systems in the Ny-Ålesund area, Kongsvegen (Hagen *et al.* 2005; Bamber *et al.* 2005; Raper *et al.* 2005) and Holtedalsfonna (Baumberger 2007). Elevation changes might also be determined from ICEsat data. However, the track spacing is relatively wide, even in Svalbard, and after filtering out data influenced by clouds, the potential data set in the Ny-Ålesund area has been found to be greatly restricted (Christopher Nuth, 2009, personal communication). Instead we use DEMs of the Norwegian Polar Institute, laser altimeter data from a centerline profile flown on Kronebreen in 2007, and ground-based kinematic GPS data from Holtedahlfonna and Kongsvegen in various years from 2000 onward (Hagen *et al.* 2005, Kohler, unpublished data). Due to the sparse temporal sampling, these profiles provide information on secular changes in ice thickness only. Therefore, we can use eq. (3) only for the determination of secular mass changes. Using this equation, we average the available profiles into simple linear functions of elevation *z*:

*z*is given relative to the geoid (i.e. as metres above sea level),

*h*

_{0}is the ice thickness change at sea level,

*h*

_{acc}is the ice thickness gain in the accumulation area of the glaciers above

*z*=

*h*

_{1}. Using a conversion factor of 0.9 mweq m

^{−1}to convert surface elevation changes in mass changes and expressing the elevation changes in mweq, for Kongsvegen α = 0.005 mweq (m yr)

^{−1}, and for Holtedalsfonna α = 0.0013 mweq (m yr)

^{−1}. For the secular changes over the observation interval, for Kongsvegen

*h*

_{0}= 3.04 mweq yr

^{−1},

*h*

_{1}= 600 m and Δ

*h*

_{acc}= 0 mweq yr

^{−1}. For Holtedalsfonna

*h*

_{0}= 1.0 mweq yr

^{−1},

*h*

_{1}= 890 m and Δ

*h*

_{acc}= 0.11 mweq yr

^{−1}. The values for Kongsvegen are constructed by a least-squares (LSQ) fit of eq. (6) to the 1996 and 2004 GPS-profile data for Kongsvegen (Hagen

*et al.*2005). For Holtedalsfonna, the values are obtained in a LSQ fit of eq. (6) to elevation differences averaged in 100 m bins using two DEMs of Kronebreen—Holtedahlfonna constructed from 1990 NPI aerial photographs and 2007 SPOT imagery, respectively.

The retreat on the local glaciers is spatially variable. For land-terminating glaciers, retreat or advance of the glacier can be accounted for by retreating or advancing glacier masks, allowing for a reasonable representation of the spatial distribution of mass loss. On Kronebreen, for example, the front is roughly stable over the past decade, and typical year-to-year changes are typically about 5–15 per cent of the total calving.

For the second approach discussed above (eq. 5), we use several long mass balance records available in northwestern Svalbard: Austre Brøggerbreen (measured since 1966), Midtre Lovénbreen (since 1967), Kongsvegen (since 1986), and Holtedalsfonna (since 2003, Baumberger 2007). The elevation at which *b*_{n} is equal to zero is the equilibrium line altitude (ELA) *z*_{EL}, which varies regionally, increasing with distance to the coast. Yearly changes in mass balance values between the various glaciers are relatively consistent in the area (Hagen *et al.* 2003, Baumberger 2007), indicating that year-to-year changes (not the absolute values) are driven by larger scale climate conditions. We use mass balance values measured on the glaciers mentioned above during the period 1998–2007 to determine a mean annual mass balance.

For Regions 2 to 5, the distance from the geodetic sites to the other glaciated regions are larger and impacts from these regions on the uplift at the geodetic sites are smaller, so a single average value can be used to describe the far-field unloading. Glacier elevation changes are smaller to the east (Nuth *et al.* 2007), and so we assume a mean value of Δ*h* for the more distant regions of Svalbard independent of elevation *z* of

*h*

_{farfield}= −0.25 mweq yr

^{−1}.

## Loading Prediction

### The geophysical model

The loading prediction were calculated using the approach of Farrell (1972), with the elastic earth model being the Gutenberg—Bullen A spherically symmetric model. Sea level changes due to the melt water have not been accounted for. There are several shortcomings of this approach that may lead to errors in the predictions. These include deviations of the crust and upper mantle underneath Svalbard from the radial structure of the earth model, lateral heterogeneities in the viscoelastic properties, viscous contributions from past recent changes in the ice load, and a contribution from sea level changes due to mass redistribution. Realistic deviations from the radial structure of the global model can changes the near-field Green's function by up to 10 per cent. Lateral heteorogeneities in the elastic crustal properties, including fault systems can also contribute significant errors (e.g. Latychev *et al.* 2005) although the effect on vertical displacements is limited. For times scale of up to decades, the viscous contribution is small and can be neglected. Therefore, for the PDIM, the elastic model is sufficient. However, if there were large ice load changes during the last centuries, then a viscous contribution from these changes might be significant. Assuming a steady-state Maxwell rheology, for example, Farrell & Clark (1976) and Clark *et al.* (1978) showed that for time scales of up to 1000 yr, the viscous contribution is on the order of 10 per cent of the elastic signal. However, Yuen *et al.* (1986) and others discussed the effect of transient rheologies on the viscous response to glacial loading and on time scales of a few centuries, transient responds might be significant. Therefore, we cannot exclude a contribution of vertical displacements induced by the unloading particularly during the last century, which has not be accounted for by any of the ice models used to compute the PGR signal. The main effect of this contribution would be to increase the PGR signal. The redistribution of the melt water in the oceans lead to an additional unloading around Svalbard. However, for relatively small ice loads like the one in Svalbard, vertical displacement due to the direct unloading of ice by far exceeds the secondary effect due to redistribution of water in the ocean with the factor being on the order of 20–50. In summary, the model assumption may lead to an error of the prediction on the order of 10–15 per cent, which is acceptable considering the uncertainties in the observations and the ice load model itself.

For the extent and height of the glaciers, we used DEMs and unpublished data provided by the Norwegian Polar Institute (Fig. 4). The loading predictions were calculated with a spatial resolution of the ice model of 500 m.

### Secular trends at the observatory

Our preferred ice load model is Model 1, which is based on eq. (6). We assume that the loading for all of north-west Svalbard can be described using the same constants as for Kongsvegen, except for the Holtedalsfonna and Isachsenfonna ice fields, which are in the latitude and longitude ranges of 78°54′N–79°10′N and 12°24′E–13°47′E, respectively, where we use the constants for Holtedalsfonna. Where we refer to Holtedalsfonna below, we refer to this area including the southern part of Isachsenfonna.

To assess the impact of the alternative approaches for the determination of the ice load changes on the loading predictions, we employ several ice load models for Region 1 (Table 2). Model 2 is the same as Model 1 but averaged with eq. (3) for the areas 1–4 in Fig. 5.

The glaciers south of Kongsvegen (south of 78°45′N; green in Fig. 5) cause uplift signals of 0.62 and 0.64 mm yr^{−1} for Models 1 and 2, respectively. For the area north of Holtedalsfonna (north of 79°10′N; blue in Fig. 5) the uplift rates are 0.47 mm yr^{−1} and 0.50 mm yr^{−1} for Models 1 and 2, respectively. However, for the area between 78 and 79°10′N (yellow and red in Fig. 5), including Kongsvegen and Holtedalsfonna, the predicted uplifts are 1.65 mm yr^{−1} and 1.08 mm yr^{−1} for Models 1 and 2, respectively. Thus, the differences between the two approaches are only significant for the middle one of the three latitudinal areas considered. Therefore, we can use either of the two models for the northern and southern areas.

We also compute the spatial fingerprint of each region in vertical displacement for Models 0 and 1. For Model 0, that is, a uniform melting rate of 1 mweq yr^{−1} over the whole of Svalbard, the predicted uplift at the Observatory is 4.6 mm yr^{−1} (Table 3). For Model 1, the predicted vertical displacement at Ny-Ålesund is 3.2 mm yr^{−1}. The total mass loss of Model 1 is equivalent to a mean melting rate over the whole of Svalbard of 0.37 mweq yr^{−1}. For each model listed in Table 2, we can express the uplift rate γ at the Observatory as a simple function of the mean melting rate *r*, that is γ = β*r*, where β has the units mm mweq^{−1}, γ is in mm yr^{−1}, and *r* in mweq yr^{−1}. For Model 0, we have

^{−1}/0.37 mweq yr

^{−1}= 8.7 mm/mweq, and we have Thus, for the same mean melting rate, Model 1 produces an uplift at Ny-Ålesund almost twice as large as Model 0. This is due to the fact that the glaciers have highest melting rates in their lower parts. Especially the glaciers Holtedalsfonna and Kongsvegen have their melting parts closer to Ny-Ålesund than their accumulation areas. Consequently, for Model 1, 88 per cent of the predicted uplift originates from Region 1, whereas for Model 0 only 66 per cent are induced by load changes in that region.

### Predicting interannual variations

The melting during summer time is more variable than snow accumulation during the winter time. Therefore, the separation of the years at the end of the melting season gives a better picture of the actual interannual changes in mass balance. The end of the melting season is approximately at September 2, which corresponds to 0.67 yr.

Annual mass balance values were determined directly from the observations described in Section 3 as the difference between snow and ice accumulation and melting over a full cycle. For example, the net mass balance for 2000 is the difference between accumulation and melting during the period from the end of the melting season in 1999 until the end of the melting season in 2000 (Table 4). To compute a weighted and normalized mean mass balance variation for the three glaciers considered, the mean of each individual mass balance series was removed and the yearly means were computed using the residuals. For the computation of the predicted annual displacements, we have assumed that we can extrapolate the mean mass balance series found for the nearby glaciers to the whole of Svalbard using the geographical pattern of Model 1. The error in this assumption is reduced due to the higher sensitivity of vertical displacements to nearby mass changes as illustrated in Table 3. Eq. (9) can be used for both secular rates and steps in mass balance. Therefore, the mean annual mass balance values for Svalbard determined by using the values in Table 4 and the spatial pattern of Model 1 can be converted into predicted annual displacements at the Observatory using eq. (9).

### The predicted spatial uplift pattern and glacial dynamics

The predictions of the vertical rates induced for the stations of the control network are compiled in Table 5. Some of the control network stations are very close to the glaciers, and therefore dynamic effects might be important.

Kongsvegen is flowing extremely slowly over its entire length, such that it essentially has no dynamics. In this case, eqs (3) and (4) coincide, and eq. (6) is relatively well known. In contrast, Holtedalsfonna has large and uncertain dynamics. It is flowing very rapidly, and the glacier is capable of transporting ice relatively rapidly from the accumulation to melting areas. Here eq. (6) has large uncertainties. To assess the impact of these uncertainties we have tested two alternative models for this glacier (Models 3 and 4; Table 2). Both models are identical to Model 1 except for Holtedalsfonna (Fig. 5, left panel, area 2), where we assume the same mass loss as Model 1, but a different dynamic regime. For Model 3, we assume no dynamics in the glacier, that is Δ*h*(*z*) in eq. (3) coincides with *b*_{n}(*z*) in eq. (5). With the assumptions of , a LSQ fit of eq. (6) to available mass balance values on Holtedalsfonna yields α = 0.00376 mweq (m yr)^{−1}, *h*_{0} = 2.92 mweq yr^{−1}, *h*_{1} = 873 m and *h*_{acc} = 0.36 mweq yr^{−1}. For Model 4, total mass loss of Holtedalsfonna is computed using eq. (6) and then divided by the area of Holtedalsfonna that results in a mean melting rate of −0.29 mweq yr^{−1}.

The effect of this modification is illustrated in Table 5 for the stations of the GPS control network. The differences are only significant for the control points closest to Holtedalsfonna. This give an indication of the impact of different dynamic behaviour of the glaciers on the spatial pattern of the uplift and demonstrates the possibility to use geodetic techniques to constrain glacial dynamics.

## Comparison of Observed and Predicted Vertical Displacements

### Interannual variations

To estimate interannual variations in observed and predicted surface displacements, we analysed the time-series with a standard LSQ method. For the determination of interannual variations, we use a model function that consists of a sum of Heaviside functions at the approximate end of the melting season in each year. Thus, the model function is given by

where the*t*

_{j}are 1995.67, 1996.67, …, 2006.67 (see Section 4.3 for an explanation of this choice). The amplitude

*A*

_{j}of each Heaviside function gives the difference between the mean before and after each

*t*

_{j}. For example, the offset for the vertical displacement at 2000.67 describes the difference between the mean height for the period 2000.67 to 2001.67 and the mean height for the period 1999.67 and 2000.67. Using eq. (10), we have determined for each individual geodetic solution annual means (Table 6). The individual solutions as well as the combined time-series of annual observed displacements are shown in Fig. 6 together with the predicted displacements (Section 4.3).

A regression between the individual solutions and the predicted annual displacements in height due to PDIM results in regression and correlation coefficients of 0.83 to 1.15 and 0.6 to 0.7, respectively (Fig. 7). As discussed in Section 2, the individual solutions are biased with respect to the secular trend (see Table 1). The regression of observed and predicted annual displacements also results in solution-specific offsets that are approximately 6.6 mm for the two GAMIT solutions, 8 mm for the GIPSY—REG solutions, 9.4 mm for the GIPSY—PPP solutions, and 7.0 mm for the VLBI solution (see row RO in Table 6). These offsets are due to a combination of the technique and analysis-specific factors mentioned in Section 2 as well as a number of common factors such as PGR (≈2 mm yr^{−1}), secular PDIM (≈3.3 mm yr^{−1}), and tectonics. However, since all solutions, independent of technique and analysis method, should contain the same interannual changes induced by PDIM, we can use these offsets to compute an unbiased time-series of average annual displacements (see column COMBINED in Table 6). We consider the resulting time-series of mean annual displacements (Fig. 6) as a good measure of the interannual variations in height at the Geodetic Observatory. A regression of the observed and predicted mean annual displacements results in a regression coefficient of 1.08 ± 0.38 and a correlation coefficient of 0.80 (Fig. 8). These results indicate that the observations very well capture the loading signal induced by the interannual variations in the ice load.

### Spatial fingerprints in secular trends

As described in Section 2, the spatial pattern of observed secular trends in height is computed as the mean of the GIPSY—PPP (plus subsequent ambiguity resolution) and the GAMIT (regional) solutions (Table 7 and Fig. 3). The predicted PDIM is computed for Model 1 in Table 2. For the PGR signal, we have computed the mean and standard deviation of 10 different model predictions of the present day signal in the vertical, using different ice histories and mantle viscosity profiles (see Plag 2006a). The standard deviation of the predictions from the ensemble mean is taken as an indication of the uncertainties of the predictions. The standard deviation decreases from 0.5 mm yr^{−1} at the easternmost network site to 0.3 mm yr^{−1} at the western margin of the network. Similarly, the PGR signal decreases from 1.9 mm yr^{−1} in the eastern part of the network to 1.2 mm yr^{−1} at the western end. We have also compared the spatial pattern in the predicted PRG signal to the 10 000-C_{14}-year shoreline (Bondevik *et al.* 1995) and find a very good agreement between these two patterns.

In Fig. 9, we compare the observed spatial pattern to the sum of the predicted PDIM and PGR. The spatial pattern of the observed trends follows very well the pattern in the predicted trends, although the East—West tilting of the pattern appears to be slightly larger for the observed trends. The optical agreement is confirmed by a regression analysis of the measured and predicted trends for the network that results in a regression coefficient of 1.26 ± 0.42 and correlation coefficient of 0.90 (Fig. 10). However, the predicted spatial pattern appears to be biased to significantly lower uplift rates (see later).

### Secular trends

At Ny-Ålesund, the combined predicted PGR and PDIM trend is 4.8 ± 0.3 mm yr^{−1} (where we did not add an error for PDIM, see Table 7). This is more than 2 mm yr^{−1} smaller than the secular trends determined from the vertical displacements observed by VLBI and GPS for the GAMIT solutions and roughly 5 mm yr^{−1} smaller than the uplift determined in the GIPSY—PPP solution (see Table 1).

## Discussion

The results of the temporal and spatial regressions presented in the previous section are based on relative values in the sense that we compare predicted and observed year to year changes or predicted and observed spatial patterns. The temporal results only depend on the precision of the observed differences in height between two consecutive years, and are less affected by potential problems in the realization of the reference frame. In particular, a secular motion between two reference frames would be absorbed by the regression offset and not affect the regression coefficient between observed and predicted year-to-year changes. The results of the spatial regression only depend on the precision of local vectors in a relatively small network. Again, this eliminates many of the problems arising from the long-term stability of the reference frame.

The good agreement between observed and predicted year-to-year variations in height as well as the observed and predicted spatial pattern in vertical trends shows that space-geodetic observations near changing glaciers or ice sheets provide valuable constraints on mass changes in these ice masses. However, the spatially uniform bias between predicted and observed secular trends of 2 to 5 mm yr^{−1}, depending on technique and analysis method (see previous section), hampers the direct inversion of geodetic observations for mass changes. Assuming that technique and solution-dependent biases have spatial scales larger than the extension of the control network, this bias between predictions and observations is likely caused by a combination of tectonic vertical motion, scale errors in PGR and PDIM and technique-dependent large-scale offsets in the observed rates.

With respect to the observations, the comparison of the time-series resulting from different analyses approaches reveals the sensitivity of the results to the analysis method. The most influential factor in terms of inter-analysis differences is most likely the reference frame, although other factors may contribute, too. The reference frames implicit in the various analyses differ in the extent of spatial filtering, and this affects in particular signals with large spatial scales, such as the seasonal signals. The free solutions are transformed into ITRF, assuming linear motion of all reference sites. However, for some of the reference sites, the ITRF velocities may not be representative due to nonlinear motion of the sites, for example, due to nearby loading or tectonic motion. Such changes in ITRF velocity can be detected through subsequent iterations.

For the secular trend at Ny-Ålesund, the discrepancy between the predicted combined PDIM and PGR signal and the observed trends is between 2 and 5 mm yr^{−1}. The large range of secular trends derived from two independent techniques and different analysis approaches is indicative of the current uncertainty in relating the space-geodetic time-series unanimously to a reference frame fixed to the centre of mass of the whole Earth system. Although an error of 2 mm yr^{−1} could be due to a secular trend of the reference frame with respect to the centre of mass of the whole Earth system (e.g. Plag 2006; Blewitt *et al.* 2006; Willis *et al.* 2006) it is unlikely that this error would be as large as 5 mm yr^{−1}. However, a combination of unaccounted tectonic uplift, reference frame errors, and technique-dependent effects may be sufficient to cause a discrepancy of up to 5 mm yr^{−1}.

Nevertheless, our model predictions may underestimate the effect of PDIM. Sato *et al.* (2006) predicted an uplift of 2.04 mm yr^{−1} at the Observatory for an assumed uniform ice loss of 0.47 m yr^{−1} that corresponds to 4.3 mm yr^{−1} uplift for a uniform melting of 1 mweq yr^{−1}. This is very close to our value of 4.6 mm yr^{−1} for Model 0. However, several factors common to both models may potentially contribute significantly to the error budget. The model derived here for the ice load changes is associated with considerable uncertainties, and this uncertainty can only be reduced through more comprehensive observations of ice surface elevation changes and mass balance variations. Moreover, the loading calculations assume a global spherical-symmetric model, with no lateral heterogeneity taken into account. However, Kaufmann & Wolf (1996) and others have shown that in Svalbard significant lateral heterogeneity of the upper mantle may impact the PGR signal, and this would also be true for the PDIM.

## Conclusions

We have shown that a significant fraction of the large secular uplift (on the order of 8 ± 2 mm yr^{−1}) observed with different space-geodetic techniques and sensors at the Geodetic Observatory in Ny-Ålesund, Svalbard, can be explained by a combination of PGR and increased PDIM. This does, however, not rule out a contribution from tectonic processes on the order of 2 mm yr^{−1}. Year-to-year variations in height, which are on the order of ±7 mm match very well those predicted based on a model of interannual ice mass changes. Similarly, the spatial pattern in secular vertical crustal motion observed on the control network agrees with the fingerprint predicted on the basis of our model for the secular PDIM.

Our study emphasizes the versatility of geodetic measurements close to glaciers to constrain mass changes in the glaciers. The interannual and secular vertical signals induced by ice load changes are significantly above the noise level of the space-geodetic time-series of surface loading at these timescales. The largest uncertainty is likely due to the ice load model. Therefore, space-geodetic observations for sufficiently dense networks could potentially be used to invert for mass changes. Combination with other observations, such as tide gauges, InSAR, satellite altimetry and GRACE, could help to constrain the mass changes to a much better level than would be possible without the space-geodetic point measurements.

### Acknowledgments

The authors would like to thank J. O. Hagen and T. Eiken for helpful discussions about glaciological issues and D. MacMillan for providing the VLBI results. GPS orbits and clocks were taken from the JPL and the International GNSS Service (IGS). Data for IGS stations were down-loaded from the IGS data archive. The work at the University of Nevada, Reno was supported by NASA funding under the Earth Science Program. Part of the work at the Norwegian Mapping Authority was supported by COST Action ES0701 ‘Improved constraints on models of Glacial Isostatic Adjustment’. The authors also thank the reviewers John Wahr and Jeanne Sauber for constructive comments that helped to improved the original manuscript considerably.