Seasonal glacier and snow loading in Svalbard recovered from geodetic observations

SUMMARY We processed time series from seven Global Navigation Satellite System (GNSS) stations and one Very Long Baseline Interferometry (VLBI) station in Svalbard. The goal was to capture the seasonal vertical displacements caused by elastic response of variable mass load due to ice and snow accumulation. We found that estimates of the annual signal in different GNSS solutions disagree by more than 3 mm which makes geophysical interpretation of raw GNSS time series problematic. To overcome this problem, we have used an enhanced Common Mode (CM) ﬁl-tering technique. The time series are differentiated by the time series from remote station BJOS with known mass loading signals removed a priori. Using this technique, we have achieved a substantial reduction of the differences between the GNSS solutions. We have computed mass loading time series from a regional Climatic Mass Balance (CMB) and snow model that provides the amount of water equivalent at a 1 km resolution with a time step of 7 days. We found that the entire vertical loading signal is present in data of two totally independent techniques at a statistically signiﬁcant level of 95%. This allowed us to conclude that the remaining errors in vertical signal derived from the CMB model are less than 0.2 mm at that signiﬁcance level. Reﬁning the land water storage loading model with a CMB model resulted in a reduction of the annual amplitude from 2 . 1 mm to 1 . 1 mm in the CM ﬁltered time series, while it had only a marginal impact on raw time series. This provides a strong evidence that CM ﬁltering is essential for revealing local periodic signals when a millimetre level of accuracy is required


INTRODUCTION
The Arctic archipelago Svalbard is exposed to climate change phenomena, the temperature is rising, the permafrost is melting, the sea level is rising, and the glaciers are retreating (Hanssen-Bauer et al., 2019).Consequences of climate change, like sea-level rise or increased land-uplift, can be observed by geodetic techniques in an accurate geodetic reference frame.On the other hand, these changes challenge the stability of the geodetic reference frame itself, e.g., the increased land uplift will deform the reference frame over time.Knowledge about the interaction between geophysical processes, crustal deformations and reference frame is mandatory to achieve the GGOS2020 goal of a reference frame with a stability of 0.1 mm/yr (Plag & Pearlman, 2009).
The geodetic observatory in Ny-Ålesund is one of the core stations in the global geodetic network.It was established during the halfdan.kierulf@kartverket.no 1990s with Global Navigation Satellite System (GNSS) antennas, Very Long Baseline Interferometry (VLBI) telescope, Super Conducting Gravity (SCG), absolute gravity points, and control networks (Kierulf et al., 2009a).
Due to Svalbard's remote location and challenging environmental conditions Ny-Ålesund was for a long time the only location with permanent geodetic equipment on the archipelago.Sato et al. (2006b,a) studied the gravity signal in Ny-Ålesund and the interaction between gravity changes and uplift.The uplift in Ny-Ålesund is not linear.It has a seasonal component, that will be studied in details in this manuscript, and an inter-annual signal induced by the long term (years to decades) evolution of glacier mass balance (e.g.Kierulf et al., 2009a).Kierulf et al. (2009b) showed that the uplift changed from year to year and that these variations are very well explained by the changes in the mass balance of the nearby glaciers.Omang & Kierulf (2011) found that also the gravity rate is changing with time.Mémin et al. (2012) showed that topography of glaciers has a significant effect on the gravity rate.
The visco-elastic response of the last ice age (Auriac et al., 2016) and the visco-elastic response of the glacier retreat after the Little Ice Age (LIA) (Mémin et al., 2014) also contribute to the uplift in Ny-Ålesund.In 2005, the Polish research station in Hornsund installed a new GNSS antenna.Rajner (2018) compared results from the stations in Hornsund and Ny-Ålesund and demonstrated that both locations have non-linear uplift.All these papers focus mainly on glacier related phenomena with time spans ranging from years to decades or thousands of years.
The most prominent variations in snowpack and glacier mass are the annual cycle with accumulation of snow each winter and melting in the short Arctic summer.The crusts elastic response of this seasonal variations results in a seasonal cycle also in the GNSS station coordinates and other geodetic equipment.The crust is also exposed to Non-tidal loading (NTL) from atmosphere, ocean and land water (Petrov & Boy, 2004;Mémin et al., 2020).
The main questions in this paper are: (1) How well do GNSS and VLBI capture the seasonal signal from glaciers and snow in Svalbard?(2) Will refining the Land Water Storage (LWS) models with a Climatic Mass Balance (CMB) model improve the loading predictions?To answer these two questions we have studied GNSS time series from six locations on Svalbard (see Fig. 1) and the VLBI antenna in Ny-Ålesund.We have used different analysis strategies both for the GNSS and the VLBI data sets.We have also filtered our time series for NTL and Common Mode (CM) signals to improve the regional accuracy.The model described in van Pelt et al. (2019) simulates glacier CMB and seasonal snow conditions, from which variations in loading from glaciers and snowpack are extracted.
In Section 2 we describe the different data sets used in this study.We describe the softwares and analysis strategies for geodetic analysis, the time series analysis, the CM filtering and the different models used for loading predictions.In Section 3, we compare the geodetic results with the loading signal from glaciers and snow, Atmospheric (ATM), Non-Tidal Ocean (NTO) and LWS.Based on this we discuss possibilities and limitations in our solutions for revealing the seasonal elastic signal.We also study the effect of refining the hydrological model with the CMB model (Subsection 3.4).

DATA AND DATA ANALYSIS 2.1 CMB model
Glacier mass change is primarily the result of surface -atmosphere interactions (affecting snow accumulation and melt), snow processes (affecting melt water retention and runoff) and frontal processes (calving and frontal ablation of tidewater glaciers).Glacier mass changes due to atmosphere -surface -snow interactions are described by the CMB, which describes the mass change of a vertical column of snow/firn/ice, in response to surface mass, and energy exchange and runoff of melt water.The CMB dominates seasonal glacier mass change, with mass gain from snow accumulation during the cold season and melt-driven mass loss during the melt season.Noël et al. (2020) have shown that for all glaciers in Svalbard the mass fluxes of precipitation (+23 Gt/yr) and runoff (-25 Gt/yr) dominate the seasonal climatic mass balance cycle in recent decades , with nearly all runoff concentrated in the summer months (June, July and August) and snow accumulating the rest of the year.These mass fluxes are much larger than the estimated mean ice discharge due to calving and frontal ablation from tidewater glaciers (7 Gt/yr Błaszczyk et al., 2009).Svalbardwide constraints on the seasonality of combined calving and frontal ablation are currently lacking and not considered here.Previous estimates on three glaciers in Svalbard however indicate that frontal ablation is more substantial in summer and early autumn than during winter and spring (Luckman et al., 2015).
Here, we use the CMB model data set, described in van Pelt et al. ( 2019), and extract weekly output for the period 1990-2018.van Pelt et al. ( 2019) used a coupled energy balance -subsurface model (van Pelt et al., 2012) to simulate CMB for all glaciers in Svalbard, as well as seasonal snow conditions in non-glacier terrain.Both the glacier and seasonal snow mass changes are accounted for.They describe weekly mass changes resulting from snow accumulation, surface moisture exchange, melt and rain water refreezing, and retention in snow, and runoff.Runoff estimates are local and no horizontal transport of water is accounted for.

Elastic loading signal
Mass redistribution results in Earth's crust deformation called mass loading (Darwin, 1882).Mass loadings are caused by the ocean water mass redistribution due to gravitational tides and pole tide (ocean tidal loading), by variations of the atmospheric mass (ATM loading), by variations of the bottom ocean pressure due to ocean circulation (NTO loading), and by variations of land water mass stored in soil, snow, and ice (LWS loading).Mass loading crustal deformations have a typical magnitude at a centimetre level (see e.g.Petrov & Boy, 2004).Love (1911) showed that the deformation caused by mass loadings can be found in a form of an expansion into spherical harmonics.Each spherical harmonic of the deformation field is proportional to the spherical harmonic of the surface pressure exerted by loading mass.The proportionality dimensionless coefficients called Love numbers that depend on a harmonic degree are found by solving differential equations.Therefore, when the global pressure field mass redistribution is known, the elastic deformation can be found by expansion of that field into spherical harmonics, scaling the harmonics by Love numbers and performing an inverse spherical harmonics expansion.
Love numbers were computed using the REAR software (Melini et al., 2015) for the Earth reference model STW105 (Kustowski et al., 2008).Time series of NTL from ATM, NTO and LWS have been used in our analysis.Input to the ATM loading is the pressure field from NASA's numerical weather model MERRA2 (Gelaro et al., 2017).The NTO loading uses the model MPIOM06 (Jungclaus et al., 2013), and the LWS loading uses the pressure field of MERRA2 model (Reichle et al., 2011).The MERRA2 model accounts for soil moisture at the depth of 0-2 meter and accumulated snow.3D displacements cause by these loadings were computed using spherical harmonics transform of degree and order 2699 and presented at a global grid 2 × 2 with a time step of 3 or 6 hours.
Then mass loading at a given station is found by interpolation.The time series of these loadings are available at the International Mass Loading Service http://massloading.net(Petrov, 2017).However, the MERRA2 numerical weather model do not adequately describe accumulation and runoff of water, snow, and ice at glaciers.It does not consider all complexity of glacial mass change processes and its resolution, 16x55 km, is insufficient to catch fine details in Svalbard.Here, we test the impact of replacing the above global model component for snow and ice with the regional snow and glacier CMB product with 1 × 1 km resolution.The model is described in Section 2.1.We have re-gridded the 1 × 1 km model to a uniform, regular, latitude-longitude grid with a resolution of 30 × 30 .The model value at a given element of the new grid is where M ab is the model value for element a,b, the r ij,ab is the distance between grid points i,j and a,b, and D is the kernel distance set to 1 km.We have computed mass loading time series, from 1990-08-05 through 2018-08-26 with a step of 7 days, at a 30 × 30 grid from the CMB output using spherical harmonic expansion degree and order 10799.This high resolution was used to correctly model the signal at stations that are located close to the edge of glaciers.
The choice of the degree/order of the expansion is determined by availability of computing resources.The higher degree/order of the spherical harmonic transform, the less errors near the coastal line.Atmospheric, land-water storage, and non-tidal ocean loading are computed with the time resolution of 3 hours and the total computation time using degree/order 2699 is about 4 years per single core for CPUs produced in 2015-2020.Since the glacier model has time resolution of 7 days we can afford to run computation with degree/order 10799 which allows to correctly model the signal at stations that are located close to the edge of glaciers.
However, it is not sufficient to replace the LWS loading computed on the basis of MERRA2 model with the mass loading computed on the basis of the CMB model.Crustal deformation at a given point is affected by mass loading not only from the close vicinity, but also from remote areas.Therefore, in order to account for loading displacement caused by mass redistribution from the area beyond Svalbard archipelago, we computed an additional series of LWS loading using MERRA2 model that was set to zero outside Svalbard archipelago.The total LWS loading displacement is: where Dmerra2 is the displacement from MERRA2 model, D merra2,svalbrard is the loading signal from the MERRA2 model that was set to zero except latitude 76 • < φ < 81 • and longitude 10 • < λ < 34 • (the area including the Svalbard archipelago) and DCMB is the displacement form the CMB model.Fig. 2 shows the high-resolution maps of the rate and amplitude of the annual signal in crustal deformation caused by the water mass change in Svalbard archipelago according to the CMB model.The parameters were estimated in a 4-parameter least square regression (mean value, rate and sine and cosine annual term) for the time series in each grid point.

GNSS data analysis
In this study we have used 30 seconds daily RINEX data resampled to five minutes, from five permanent GNSS stations on Svalbard (NYAL, NYA1, LYRS, SVES and HORN), and one station on Bear Island (BJOS) 240 km south of Svalbard (see Fig. 1).All stations are located close to existing settlements with infrastructure like power supply and communication.We have also used data from station, HAGN, 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 is recorded for 24 hours once a week to save power until the sun is back.Data is downloaded during a field trip once a year.
GNSS data are analyzed with the program packages Gamit/Globk (Herring et al., 2018) and GipsyX (Bertiger et al., 2020).The GipsyX software is using undifferentiated observations.We are using the Precise Point Positioning (PPP) approach (Zumberge et al., 1997) and the solutions are in the International GNSS Service (IGS) realization of International Terrestrial Reference Frame (ITRF2014) (Altamimi et al., 2016) through the Jet Propulsion Laboratory (JPL) orbit and clock products.We distinguish between the GispyX-FID and GipsyX-NNR solution, whereby either JPL fiducial (FID) or No-Net-Rotation (NNR) orbit and clock products are applied.The NNR products are only constrained via three no-net-rotation parameters to the ITRF2014 solution, whereas the FID products are tied in addition with three translation and one scale parameter to ITRF2014 (Bertiger et al., 2020).Gamit software uses double difference observations.To ensure a good global realization in ITRF2014 of the Gamit solution a global network of approximately 90 global IGS stations was analyzed and combined with the Svalbard stations before transforming to ITRF2014.The global stations were all stable stations with long time series.Daily coordinate time series are extracted from these solutions.
The two stations in Ny-Ålesund belong to the IGS network and are analyzed by several institutions, University of Nevada, Reno (UNR) (Blewitt et al., 2018), JPL (Heflin et al., 2020), and Scripps Orbits and Permanent Array Center (SOPAC) (Bock & Webb, 2012).NYAL and NYA1 are also included in the latest ITRF2014 (Altamimi et al., 2016) solution.Key parameters for the different analysis strategies are given in Table 1.
The time series are analyzed with Hector software (Bos et al., 2008).We have used the following model function: where A is the constant term, B is the rate, Cj is the amplitudes of the sinusoidal constituents, and φj is the corresponding phases.
We have assumed that the temporal correlation in the time series are a combination of white noise and flicker noise.We have used data from 2010-01-01 until 2018-10-01 in all the GNSS results and comparisons.This limited time period ensures that we have the same time period for all the stations (except HAGN which was established in 2013), no breaks due to equipment shift, and the time series overlap with the CMB model (see Section 2.1).
The time series for the vertical component of the Gamit-NMA solution is plotted in Fig 3.

VLBI
The VLBI station NYALES20 participated in 2183 twenty four hour observing sessions from 1994-10-04 to 2020-10-19.We ran several solutions.Solution s1 was obtained using the geodetic analysis software Where (see Kirkvik et al., 2017, for more details).VLBI observing sessions were individually analyzed with the following approach: a priori station coordinates were taken from ITRF2014 including the post-seismic deformation models or VTRF2019d (IVS update of ITRF2014) for newer stations.To define the origin and the orientation of the output station position estimates, tight no-net-translation and no-net-rotation with respect to ITRF2014 were imposed.A priori radio source coordinates were taken from the ICRF3 S/X catalog (Charlot, P. et al., 2020) and corrected for the galactic aberration.The source coordinates were not estimated.A priori Earth orientation parameters were taken from the C04 combined EOP series consistent with ITRF2014.The Earth orientation parameters, polar motion, polar motion rate, UTC-UT1, length of day, and celestial pole offsets, were then estimated for each session.In addition, troposphere and clock parameters was estimated.Key parameters for the VLBI solutions are included in Table 2.
Solution s2 was obtained using VLBI analysis software suite pSolve (http://astrogeo.org/psolve).Source position, station positions, station velocity, sinusoidal position variations at annual, semi-annual, diurnal, semi-diurnal frequencies of all the stations, were estimated as global parameters in a single least square solution using all dual-band ionosphere-free combinations of VLBI group delays from 1980-04-12 to 2020-12-07, in total 14.8 million observations.There are 28 stations that have discontinuities due to seismic events or station repair.These discontinuities and associated non-linear motion was modeled with B-splines with multiple knots, and the B-spline coefficients were treated as global parameters.In addition to global parameters, the Earth orientation parameters, pole coordinates, UT1, their first time derivatives, as well as daily nutation offsets are estimated for each observation session individually.Atmospheric zenith path delay and clock function are modeled with B-splines of the 1st degree with time span 60 and 20 minutes, respectively.A so-called minimum constraints on station positions and velocities and source coordinates were imposed to invert the matrix of the incomplete rank.These constraints require that the net translation and rotation station positions and velocities of a subset of stations be the same as in ITRF2000 catalog and net rotation of the so-called 212 defining sources be the same as in ICRF.It should be noted that s2 solution is independent on the choice of the a priori reference frame, i.e. change in the a priori position does not affect results.The data reduction model included modeling thermal variation of all antennas, oceanic tidal, NTO ATM and LWS loading with one exception, where for station NYALES20 the following LWS model were used Dmerra − D merra,svalbard .Implying that the a priori model totally ignores mass loading exerted by water mass redistribution in Svalbard.
The VLBI network is small and heterogeneous: different stations participate in different experiments.Therefore, the time series of station position should be treated with a great caution: the estimate of the position change of station X affects the position estimate of station Y because of the use of the net translation and net rotation constraints to solve the system of the incomplete rank.An alternative approach to processing time series is estimation of admittance factor.We assume that the time series of the displacement in question d(t) is present in data as a • d(t) where a is a dimensionless parameter called an admittance factor that is assumed constant for the time period of observations.The admittance factor describes what share of the modeled signal is present in observations.
We noticed that seasonal crustal deformations of NYALES20 positions are periodic but not sinusoidal.The shape of these variations is surprisingly stable with time (See Fig. 4).We decomposed the mass loading signal into four components: seasonal, interan-nual, linear trend, and residuals.The decomposition was performed in three steps.First, the mass loading time series were filtered with the low-pass Gaussian filter, which provided a coarse interannual signal (IAV(t)).Second, the time series were folded of the phase in a form p = (t − t0)/∆t, where t is time, t0 is the reference epoch 2000.0,∆t is the period (one year), and then smoothed.That provided a coarse estimate of the seasonal signal (SEA(t), blue curve in Fig. 4).Then we adjusted parameters A, B, ai, si of the decomposition of the loading displacements D(t) described by the Eq 4, using a single least square solution: where where Bi is the basis spline of the 3rd degree with the pivotal knot i. Fig. 4 illustrates the seasonal component of the loading signal at NYALES20.A thin red line at the plot shows result of the best fit of the sinusoidal signal.However, the sinusoidal model provides a poor fit to the data with errors reaching 40% of the seasonal signal.All constituents of this expansion for NYALES20 are shown in Fig. 5.
In solution s3 we did not estimated annual and semi-annual sinusoidal variations of NYALES20 positions, but estimated admittance factors for the up, east, and north components of the IAV(t) + SEA(p(t)) mass loading time series.In contrast to estimation of sinusoidal variations, the shape and phase of the signal remains fixed when we estimate admittance.The adjusted parameter is the scaling factor of the modeled displacement magnitude.The power of this approach is that it allows us to evaluate quantitatively the amount of the modeled signal in data, and test a statistical hypothesis that all model signal is present in the data.
The results of admittance factor estimation are presented in Table 3 in row ADM TOT.Then we estimated the admittance factor for the seasonal SEA(p(t)) and interannual variations IAV(t) separately in the s4 solution.

Gravimetry/SCG
We use gravity measurements from two SCG instruments covering the period 1999 to 2018 to estimate gravity change.Gravity measurements from 1999 to 2013 and 2014 to 2018 are collected with C039 and iGrav012 SCG instrument, respectively.The original gravity measurements have a spacing of 1 second, giving a total of    approximately 620 million measurements.They are re-sampled every minute using a symmetric numerical Finite Impulse Response (FIR) zero phase low-pass filter with a cut-off at 120 seconds (Wenzel, 1996).Data was then cleaned for outliers and earthquakes.We corrected for the effect of air pressure using the value of -0.422 ± 0.004 µGal/hPa found by Sato et al. (2006a).Both the solid earth and ocean tides are removed from the gravity data by estimating a synthetic tide based on Hartmann & Wenzel (1995) tidal model and a set of tidal parameters.The synthetic tide is estimated using ETERNA 3.4 (Wenzel, 1996).
We also estimated and removed the instrumental drift by comparing to absolute gravity measurements.We estimated a linear drift (using unweighted least squares) by comparing to ten AG measurements (2000, 2001, 2002, 2004, 2007, 2010, twice in 2012, 2014, 2017).The estimated value is -2.74107 ± 0.17 µGal/year.Finally, we re-sampled the data first every 5 minutes and then every 1 hour using a symmetric FIR zero phase filter (cutoffs 1250 sec and 2 hours, respectively) and then to daily values using a flat filter.

Filtering of Common Mode and elastic loading signal
It is well known that stations in a region can have a spatially correlated signal, a so-called CM signal (Wdowinski et al., 1997), and that removal of the CM signal can reduce noise in the time series.The CM signal could come from the GNSS analysis strategy and from the strategy for reference frame realization.It could come from mismodeled orbit, clocks or EOPs, or through unmodeled large scale hydrology or atmospheric effects.To remove such signal either CM filtering, Empirical Ortoghonal Functions (EOF) or regional reference frame realization, can be used.All these methods presuppose that we have stations exposed to the same undesirable CM signal.In Arctic areas, we have limited access to nearby stations.All stations on Svalbard are exposed to similar signals from glaciers, using one or several of these stations for removal of the CM signal will not only remove the CM signal, but also the real elastic signal from snow and ice.
The station BJOS at Bear Island is located 240 km south of Svalbard.The Island is small and surrounded by ocean and the local loading signal from ice and snow is approximately 10% of the signal in Ny-Ålesund (see Tabel A2).It is the closest GNSS station outside Svalbard.Time series of the BJOS station are used to estimate the CM signal.Time series for the BJOS station, and hence CM filtering, are only included in the solutions computed by the authors (Gamit-NMA, GipsyX-FID and GipsyX-NNR) and not available in the external solutions (SOPAC, JPL, UNR and ITRF).The CM filtered time series for the i-th station is then: where t is the epoch and H i GN SS and H BJOS GN SS is the time series for station i and BJOS respectively.
The CM filtering removes the common error signal at the stations as well as real measured signal at BJOS.If the station at Bear Island has an unique unmodelled loading signal not present in other Svalbard stations, this unique signal will be erroneously subtracted also from the other stations.
To CM filter a time series where the signal from a loading model is removed, the loading signal for the station(s) used in the CM filtering has to be removed as well.In our case, the loading signal was subtracted both for the BJOS time series before computing the CM signal and for the other Svalbard time series before the CM filtering.The final Svalbard time series are cleaned for both the regional CM signal over Svalbard and Bear Island and the estimated load signal.The CM filtered time series for station i is then: where t is the epoch, H i GN SS is the observed time series, H i L is the estimated loading signal, and CM is the common mode signal.As described earlier, we use the time series from BJOS to estimate the CM signal, but since we remove the estimated loading signal from the time series, we have to remove the loading signal from BJOS time series before computing the CM signal.Therefore, we get ).( 7)

Isolating the elastic signal from glacier and snow
The signal in HCM,L(t) (Eq.7) includes all vertical motions not accounted for in the loading models or CM filtering, e.g., unmodeled loading, Glacial Isostatic Adjustment (GIA), tectonics, and noise.Assuming that the GIA and the tectonic component are linear, the left hand side can be written HCM, where LIN is the linear part and ε contains the noise.The noise includes unmodeled loadings, but also station dependent effects like multipath, atmospheric effects, not use of individual antenna calibration, and thermal expansion of antenna monument.Possible unique unmodelled signal from BJOS will also map into the noise term.Splitting the load signal into a signal from glacier and snow, HGS, and other non tidal loadings, HNT L * , we can rewrite Eq. 7 into: i.e. we have isolated the linear part and the elastic signal from glaciers and snow as a sum of known terms.

RESULTS AND DISCUSSION
The vertical component of the different GNSS solutions in Ny-Ålesund and Bear Island as well as the NTL signal are included in Table 4.In addition, the Where results (Solution s1) from the NYALES20 VLBI antenna and the SCG in Ny-Ålesund are included.Some of the time series are plotted in Fig. 6.The horizontal components of the different GNSS solutions are included in Appendix, Table A1.The loading signals for all the different loading models are included in the Appendix, Table A2.

Determination of the loading annual signal
As shown in Kierulf et al. (2009b) the uplift in Ny-Ålesund varies from year to year.Consequently, trends from different time periods can not be compared directly.We have chosen to use the time interval 2010 until 2018 for time series analysis, except for the ITRF2014 time series that ended in 2014.The estimated uplift for the Ny-Ålesund stations agree below the uncertainty level.
The annual signal in Ny-Ålesund varies between the solutions both in phase and amplitude (Table 4, Table A1 and Fig. 7).This implies that the choice of GNSS analysis strategy has a noticeable impact on the estimated seasonal variations.Martens et al. (2020) found similar differences in the estimated annual signal when they compared GNSS time series, in US and Alaska, based on different analysis strategies.Such variations make direct geophysical interpretation of the periodicity in GNSS time series difficult.
The measured vertical annual signal (Table 4) is smaller than the estimated NTL signal for the GipsyX solutions and larger than the estimated NTL signal for the Gamit solutions.The phase of the GNSS solutions are delayed relative to the NTL signal with between 30 and 70 degrees (corresponding to a delay between one and two and a half months) in Ny-Ålesund.
The CM-filtered solutions are closer to the expected signal from NTL and we have less differences between the GipsyX and Gamit solutions, see Fig. 7.The annual signal found with the Where software for VLBI has a smaller amplitude, but the phase is close to the phase estimated from the loading modeling.
The phase of the gravity signal in Ny-Ålesund is close to the phase of the loading models.The annual amplitude in gravity is 3.45 µGal.For a spherical and compressible Earth model the elastic gravity variations can be converted to vertical position variations using a ratio of -0.24 µGal/mm (Mémin et al., 2012).Using this ratio we got an yearly amplitude of 14.4 mm.This is much larger than the estimated annual loading signal of around 4.0 mm in Ny-Ålesund.However, the gravity variations also depends on the direct gravitational attraction.Mémin et al. (2012) discussed how the location of the load affect the ratio between gravity and uplift.Both the distance to and the relative height of the load have an impact.The large annual signal implies mass changes at locations with negative relative heights close to the station.
The SCG-instrument in Ny-Ålesund gives a combined signal from three glacier related factors.The visco-elastic response from past ice mass changes, the immediate elastic response of the ongoing ice mass changes, and the direct gravitational attraction from the ongoing ice mass changes on the glaciers (see Mémin et al., 2014;Breili et al., 2017).The two latter have a clear influence on the annual signal.In addition, soil moisture and accumulated snow close to and mainly below the gravimeter, have a much stronger effect on gravity than on displacements.
Quantifying the gravity signal from these nearby hydrological factors are demanding and out of the scope of this paper.However, they, as well as glaciers, are forced by temperature and precipitation.We assume that they are in phase with the elastic uplift signal.A gravimeter measures gravity changes directly, while VLBI and GNSS evaluate site positions from analysis of observations at a network, and a position estimate of a given station in general depends on measurements at other stations of the network.The phase of the SCG time series is therefore an independent measure of the variations in Ny-Ålesund and the result coincides with the results from the other techniques.

Determination of the loading admittance factors
We found the admittance factor from VLBI solutions for the seasonal vertical displacement does not deviate from 1.0 at a 2σ level, i.e. the LWS signal is fully recovered from the data.At the same time the departure of the admittance factor from 1.0 for the horizontal loading components implies there is a statically significant discrepancy between the computed loading signal and the data.It should be noted that the magnitude of the seasonal signal in North direction is only 0.15 mm and the signal itself is just too small to be detected.The admittance factor for the interannual signal is significantly different from 1.0, which indicates that the loading signal alone cannot explain it.We made an additional analysis to find the admittance factor for the glacier and snow loading signal at the GNSS stations in Svalbard.We computed mass loading for all the GNSS stations in Svalbard and fitted it to the GNSS time series using reciprocal formal uncertainties as weights.Then we computed the χ 2 per degree of freedom of the fit and scaled variances of admittance factor estimates by this amount.Table 5 shows the estimates of the admittance factor from the differenced GNSS time series using Eq. 8. Similar to the VLBI case, the admittance is very close to 1.0 for the vertical seasonal signal (row "All") and it is far away from 1.0 for the interannual signal and the horizontal signal.
Two factors may cause poor modeling of the interannual signal.First, calving and frontal ablation are not included in the CMB model, and therefore, lacks this contribution.Second, other loadings, for instance non-tidal ocean loading may contribute.The seasonal signal has a very specific time dependence pattern, and the approach of admittance factor estimation exploits the uniqueness of this pattern, while the pattern of the interannual signal is more general.
The analysis of the admittance factors give several important results in addition to the very good agreement between the different estimated vertical seasonal components.Analysis of observations shows that the CMB model provides prediction of the vertical mass loading with 1-σ errors of 5%, which corresponds to 0.1 mm.We have a bias wrt to the model of 0.2 ± 0.1 mm, and this bias is not statistically significant at a 95% level (2 sigma).We conclude that analysis of the data from two totally independent techniques, VLBI and GNSS, proves there is no statistically significant deviation at a 95% significance level between the seasonal vertical mass loading signal based on the CMB model and observations of both techniques.

Geodynamical interpretation
To study the time series ability to capture the loading signal from glaciers and snowpack changes, we used the CM filtered time series.Other known loading signals were removed using Eq. 8. To have more robust time series, in the following discussion, we have used averaged GNSS time series.The averaged time series are the weighted mean of the daily values from the Gamit-NMA, the GipsyX-NNR and the GipsyX-FID solutions.The annual periodic signal and the linear rate for the time series in Eq. 8 are included in Table 6 together with the elastic signal from glaciers and snow.Detailed results for the individual GNSS solutions are included in the Appendix, Table A3.
The amplitudes of the estimated loading signal from glaciers and snowpack vary with latitude and longitude and depend on the amount of surrounding glaciers and land masses (see Fig. 9).The station HAGN in the middle of the glacier Kongsvegen has the largest estimated annual loading signal, while the westernmost stations NYAL/NYA1 and HORN have the smallest.The GNSS stations SVES and LYRS are located in central parts of Svalbard and here the measured vertical annual signal agrees with the estimated loading signal at the uncertainty level.For the stations closest to the west coast NYAL, NYA1 and HORN the measured amplitudes are slightly larger than expected from the variations in glaciers and snow (∼ 0.7 mm, ∼ 1.0 mm and ∼ 0.5 mm resp.).Although the admittance factor for all stations combined show very good agreement for the seasonal component, the admittance factor for individual stations in Ny-Ålesund and Hornsund are slightly above one, implying that the observed amplitude is somewhat larger than the prediction from the CMB model.
The larger vertical amplitude at NYAL, NYA1 and HORN might be due to lower precision of the CMB models in areas with more variable coastal climate, changes in groundwater and surface hydrology, and seasonal variability in calving/frontal ablation of glaciers.Especially, calving is assumed to be seasonally dependent with higher incidents during summer (when ice flows faster).This may explain a higher observed amplitude in areas like Ny-Ålesund and Hornsund, which have a lot of nearby large calving glaciers.However, the deviation of admittance factors from 1.0 for these stations is still within 2σ of the statistical uncertainty.Longer time series are needed to establish whether there is a statistically significant deviation of observations from the model for these individual stations.
The phase of the vertical loading signal from glaciers and snow varies with only a few days over Svalbard, and corresponds to a maximal value after the end of the melting season, in mid-October.The phase of the GNSS time series agrees with the glaciers and snow signal from the CMB models within a few weeks.
The predicted horizontal seasonal signal is smaller.It is around 0.2 mm in the north component for all locations except HORN.HORN is located in the south of Svalbard with the majority of glaciers located to the north.Consequently the north amplitude is larger, 0.7 mm.The other stations have glaciers both to the north and to the south and the loading signals are cancelled out.The east annual signal varies from 0.7 mm (NYAL, NYA1) to 0.0 mm (SVES), depending on their locations relative to the glaciers.
We see that the estimated annual signal for the GNSS stations NYA1, HAGN, LYRS and HORN agree at 2-sigma level in the east component.HORN is capturing the larger north annual signal for this location.However, the small north signals for the other stations are too small to be detected.LYRS and SVES have a large north annual signal, 1.1 mm resp.2.3 mm larger than expected load-ing signal.We have not established the origin of this discrepancies.Possible reasons are uneven thermal expansion of the antenna monument (steel mast) and artifacts of the atmosphere model.We will investigate these signals in the future.
Our CMB model is limited by Svalbard archipelago.Glacial loading at other islands, such as Iceland and Greenland, can bring a noticeable contribution.Using Green's function from the ocean tide loading provide of Scherneck (1991) † , assuming 100 Gt seasonal ice cycle on Greenland (see figure 6 in Bevis et al., 2012), and the average distance from Greenland to Svalbard of 800 km, we get coarse estimate of the the amplitude of mass loading signal Svalbard due to glacier in Greenland: 0.38 mm.In order to get a more refined estimate of the magnitude of such a contribution, we used Loomis et al. (2019) mascon solution for Greenland for 2008 from processing GRACE mission.The mascon for Greenland at a regular grid 0.5 • × 0.5 • is provided ‡ with a monthly resolution in the height of the column of water equivalent that covers the entire land and is zero otherwise.The mascon excludes the contribution of the atmosphere and ocean but retains the contribution of land water, snow, and ice storage.
We have converted the height of the column of water equivalent to surface pressure and computed the mass loading from the mascon using the same approach as we used for computing atmospheric and land water storage loading using spherical harmonic transform of degree/order 2699.Since MERRA2 land water storage model covers Greenland, we computed mass loading two times: the first time using the total surface pressure from the mascon and the second time after subtracting the pressure anomaly from MERRA2.In the latter case the resulting mass loading signal provides a correction to the mass loading from MERRA2 model for the contribution of glacier derived from GRACE data analysis since MERRA2 has large errors of modeling glacier dynamics.The results are shown in Figure 8 after removal linear trend.The residual signal has amplitude 0.15 mm and is in phase with the mass loading signal from Svalbard glacier while the total signal has amplitude 0.28 mm.This estimate agrees remarkably well with our coarse estimate.Accounting for the contribution of glacier mass loading using GRACE data reduces the admittance factor from 1.10 to 1.04 for the seasonal vertical component of the VLBI solution.
We exercise a caution in results of processing GRACE data, A thorough analysis of systematic errors of mass loading signal from GRACE mascon solution requires significant efforts and is beyond the scope of the present manuscript.However, our estimates shows that the contribution of glaciers in Greenland is not dominated and its accounting improves the agreement of the CMB model with VLBI and GNSS observations.Table 6.Vertical rate and annual signal for GNSS stations in Svalbard.GS are the elastic loading signal from ice and snow.GNSS-CM is the time series using Eq. 8. Max uplift is the date of the maximum value for the annual signal.In the previous section we examined how the different GNSS time series were able to capture the elastic loading signal from local ice and snow changes.In this paragraph we will discuss the effect of removing the loading signal from the time series, both on the unfiltered time series and the CM filtered time series.We will in particular look at the effect of replacing the global hydrological model with a regional CMB model.In the discussion we used an averaged time series from the GNSS solutions; Gamit-NMA, GipsyX-FID and GipsyX-NNR.Due to limited observations during winter the HAGN time series are not directly comparable with the other time series and therefore not included in this discussion.
We have used time series where no loading models where removed and two time series with slightly different loading time series removed (L1, L2).The first loading model, L1, is the sum of the loadings from ATM, NTO and the total LWS signal from the merra2 model.The second loading model, L2, equlas L1 except that the regional LWS signal in merra2 is replaced with the glacier and snow signal in the CMB model using Eq. 2.
HGNSS − Li contains the unmodeled signal in the GNSS time series after removing the modeled loadings, Li.The signal can be presented as a sum of linear trend (from e.g., GIA and tectonics) and noise including unmodeled loading signals.I.e.HGNSS(t) − Li(t) = LIN (t) + ε(t).Also the CM filtered time series (HCM , Eq. 6) and the CM filtered time series where the loading signal is removed (HCM,L i , Eq. 7) can be presented as a sum LIN (t) + ε(t).To examine the quality of the time series using different filtering and loading models we estimate the Root Mean Square (RMS) and annual signal in the noise time series ε(t).The annual signal in ε(t) is the remaining annual signal after removing the loading model Li.A large annual signal in ε(t) indicate that we have remaining unmodelled periodic signal in the time series after removing load Li.The RMS is a measure of the remaining noise in the time series.The results for the averaged time series are included in Table 7. Results for the individual solutions are included in Table A4 in the Appendix.
We see that removal of the loading signal reduce the RMS values on average by 11%, while replacing the regional hydrological signal with a CMB model reduce the RMS with 13%.The improvements for the CM filtered time series are less, 4% and 6%, respectively.The removal of the CM eliminates part of the elastic loading signal, and this may explains the lower reduction for these series.Both for the unfiltered and the CM filtered time series the RMS are reduced with 2-3 % when we replace the regional signal in the merra2 with the glacier and snow signal from the CMB model.The RMSs are very little affected by the CM filtering (4th vers.1st column in Table 7).
Removing the NTL (ATM, NTO and LWS including the glaciers and snow) from the observed time series have an effect on the daily noise scatter (RMS), but very little effect on the annual signal.This implies that removal of the NTL reduces the daily scatter in the GNSS time series.It also implies that the periodic signal is dominated by other factors.As we saw in Section 3 this annual signal depends on the analysis strategy.We conclude that we have an analysis strategy dependent effect in the periodic signal.
The amplitude of the time series is reduced after the CM filtering.The amplitude of the annual signal in the CM filtered time series using load models L1 is reduced to 2.1 mm.The largest effects are when we use load model L2 and the CM filtered time series.For this solution the averaged annual loading signal is 1.1 mm, one third of most other combinations of filtering and loading models.
For the horizontal components including CM filtering and removing the NTL affects the RMS only marginaly.The amplitude of the horizontal seasonal signal reduces from on average 0.6 mm to 0.4 mm for NYAL, NYA1 and HORN when both CM filtering and NTL loadings are removed.The stations LYRS and SVES have much larger amplitudes and no improvements after filtering.
Note, the glacier model used in this study is not able to capture glacier dynamics like continuous flow of ice towards the glacier front, or more dramatic phenomena such as glacier surging (see e.g., Morris et al., 2020;Dunse et al., 2015).These dynamic effects provide a significant contribution to the total glacier mass balance and uplift, especially, on time scales from years and longer (see Kierulf et al., 2009b, for more on the effect of glacier dynamics on the uplift).The linear elastic uplift signal from the CMB models is not sufficient to fully describe the elastic uplift from ice and snow changes over longer time scales.

CONCLUSIONS
In the introduction two questions were asked: (1) How well do GNSS and VLBI capture the seasonal loading signal from glaciers and snow on Svalbard?(2) Will refining the LWS models with a CMB model improve the loading predictions?To answer these questions a network of seven permanent GNSS stations were analyzed with different analysis strategies and softwares.The different time series were studied and compared with loading predictions from ATM, NTO, LWS including glaciers and snow.
We found large discrepancies between the different analysis strategies, both in phase and amplitude, while the estimates of longterm trend were more consistent.This implies that a direct geophysical interpretation of raw GNSS time series is problematic.To overcome this problem, we performed CM filtering utilizing the data from the nearby station at Bear Island.The elastic loading signal was removed from the time series before the CM filtering.The CM filtered time series gave a much better agreement.This confirmed our initial conjecture that the origin of the discrepancies in the raw time series are due to differences in the analysis strategy in the GNSS data processing.The agreement of CM filtered time series strengthened our confidence that we investigated a real geophysical signal, and not artifacts of data analysis.
We have decomposed the LWS signal into the seasonal and interannual signals, and estimated admittance factors from VLBI data and GNSS CM filtered time series.The admittance factors estimates from vertical seasonal constituents for both techniques do not deviate from 1.0 at a 2σ level.Therefore, we conclude that the entire mass loading signal is present in data from totally independent technique at a statistical significance level of 95%.The 1σ uncertainty of admittance factors corresponds to 0.1 mm.This implies that by using the CMB model we can predict seasonal vertical mass loading displacements on Svalbard with the same level of accuracy, and that predicted errors are less than the observation errors.The interannual loading signal was not recovered from observations.Further work is required to explain these discrepancies.However, since calving and frontal ablation are not included in the CMB model, they may contribute to these discrepancies.

Figure 1 .
Figure 1.Geodetic network on Svalbard.The location NYAL include the GNSS stations NYAL and NYA1, the VLBI antenna NYALES20 and the SCG instrument.

Figure 2 .
Figure 2. Crustal deformations due to glacier and snow loading according to the CMB model.The panels are: the rate of change (left) and annual signal (right).Important: the CMB model does not account for mass loss due to frontal ablation and calving.

Figure 3 .
Figure 3. GNSS vertical time series for the Gamit-NMA solutions.The color coded time span is the data period used in this study.The black curves are the model function fitted to this period.The HORN station was moved to a new location in 2009.The time series are shifted with respect to each other to improve readability.

Figure 4 .
Figure 4. Folded periodic up LWS mass loading displacements of NYALES20 after removal of the slowly varying constituent.The thick blue line shows the estimate of the seasonal constituent.Green dots show the mass loading signal after removal of the interannual constituent.A red thin line shows a sinusoidal fit in a form a cos 2π p + b sin 2π p, where p is the phase of the seasonal signal in turns. 1990

Figure 5 .
Figure 5. Three constituents of the vertical LWS mass loading at station NYALES20.The thick blue line shows the interannual variation, the green thin line shows the seasonal component, and red dots in the bottom shows the residual signal.The residual signal is artificially shifted by -8mm.The linear trend is removed and not shown.

Figure 6 .
Figure 6.A selected set of detrended time series for Svalbard.The time series are: GipsyX-NNR for NYA1 (upper left), NTL (ATM, NTO, LWS including the glaciers and snow) signal in Ny-Ålesund (middle left), the GipsyX-NNR time series for NYA1 after removal of NTL and CM filtering (lower left), GipsyX-NNR for BJOS (upper right), gravity from the SCG in Ny-Ålesund (middle right) and the NYALES20 Where solution (lower right).The gravity values are converted to millimeter using the convertion ratio -0.24 µGal/mm(Mémin et al., 2012).

Figure 7 .
Figure 7. Seasonal signal in Ny-Ålesund (NYA1).The left panel shows the sum of the annual and semi-annual sinusoidal signal for the time series, the right panel shows the same results relative to the NTL signal (ATM, NTO, LWS including the glaciers and snow).The upper most five curves are from time series analysis of the raw time series.The sixth and seventh curves are CM-filtered time series.The bottom curve of the left panel is the estimated NTL signal.The curves are shifted with respect to each other to improve readability.

3. 4
Figure8.The seasonal contribution of glaciers in Greenland to the vertical displacement of NYALES20 after removal of the slowly varying constituent from processing GRACE data.The green line shows the total contribution, the shadowed blue line shows the residual contribution after subtraction of land water storage pressure from MERRA2 model.The horizontal axis is a phase of the seasonal signal in phase turns.
Annual signal for GNSS stations in Svalbard.The bars are the amplitude and the vectors are the phase.Blue is from the loading prediction from glaciers and snow red is from the GNSS stations.

Table 1 .
GNSS analysis strategies.(*) Elevation dependent site by site functions, where a and b are estimated based on postfit editing of residuals from each station.E is the elevation angle.

Table 3 .
Admittance factors of NYALES20 displacements caused by LWS loading.The first row, ADM TOT shows the admittance factor estimate from s3 solution of the total mass loading signal.Rows ADM SEA and ADM IAV shows estimates of the seasonal and interannual constituents of the loading signal from s4 solution respectively.

Table 4 .
Mémin et al. (2012)emi-signal in Ny-Ålesund and Bear Island.The parameters are estimated trend and annual signal estimated using Eq.3).The results are for different GNSS solutions, VLBI, SCG and NTL in Ny-Ålesund and Bear Island.In the VLBI time series a pure white noise model is assumed.The gravity values ( * ) are converted to millimeter using the convertion ratio -0.24 µGal/mm fromMémin et al. (2012).CM is the CM filtered time series described in Section 2.6.NTL is the sum of non tidal elastic loading signal from ATM, NTO and LWS including the load from snow and glacier from the CMB model.

Table 5 .
Admittance factors for the vertical component of GNSS station in Svalbard caused by the glacier and snow loading.ADM SEA and ADM IAV show estimates of the seasonal and interannual admittance factors.

Table 7 .
Yearly amplitude and RMS in the time series.In each column the three parameters are amplitude of yearly signal in mm, RMS of the time series in mm and changes in RMS relative to the unfiltered time series in percent.The numbers are in mm.H GN SS , L i and H CM,L i are explained in the text.

Table A1 .
Trend and annual signal in Ny-Ålesund and Bear Island.The parameters are estimated trend and annual signal estimated using Eq. 3) for the North and East components.

Table A2 .
NTL vertical variations at GNSS stations in Svalbard.

Table A3 .
Vertical rate and annual signal for GNSS stations in Svalbard.The results are based on the time series using Eq. 8. Max uplift is the date of the maximum value for the annual signal.

Table A4 .
Yearly amplitude and RMS in the time series.In each column the three parameters are amplitude of yearly signal in mm, RMS of the time series in mm and changes in RMS relative to the unfiltered time series in percent.The numbers are in mm.H GN SS , L i and H CM,L i are explained in the text.