Postglacial gravity change in Fennoscandia—three decades of repeated absolute gravity observations

For the ﬁrst time, we present a complete, processed compilation of all repeated absolute gravity (AG) observations in the Fennoscandian postglacial land uplift area and assess their ability to accurately describe the secular gravity change, induced by glacial isostatic adjustment (GIA). The data set spans over more than three decades and consists of 688 separate observations at 59 stations. Ten different organizations have contributed with measurements using 14 different instruments. The work was coordinated by the Nordic Geodetic Commission (NKG). Repre-sentatives from each country collected and processed data from their country, respectively, and all data were then merged to one data set. Instrumental biases are considered and presented in terms of results from international comparisons of absolute gravimeters. From this data set, gravity rates of change (˙ g ) are estimated for all stations with more than two observations and a timespan larger than 2 yr. The observed rates are compared to predicted rates from a global GIA model as well as the state of the art semi-empirical land uplift model for Fennoscandia, NKG2016LU. Linear relations between observed ˙ g and the land uplift, ˙ h (NKG2016LU) are estimated from the AG observations by means of weighted least squares adjustment as well as weighted orthogonal distance regression. The empirical relations are not signiﬁcantly different from the modelled, geophysical relation ˙ g = 0 . 03 − 0 . 163( ± 0 . 016)˙ h . We also present a ˙ g -model for the whole Fennoscandian land uplift region. At many stations, the observational estimates of ˙ g still suffer from few observations and/or unmodelled environmental effects (e.g. local hydrology). We therefore argue that, at present, the best predictions of GIA-induced gravity rate of change in Fennoscandia are achieved by means of the NKG2016LU land uplift model, together with the geophysical relation between ˙ g and ˙ h .


I N T RO D U C T I O N
Glacial isostatic adjustment (GIA) is the response of the Earth to changing loads on its surface due to build-up and ablation of ice sheets and glaciers. The response includes changes in shape (deformation), gravity potential, stress and rotation of the Earth (Wu & Peltier 1982). The effects of GIA that are presently observed result from several glaciations with ice sheets covering large parts of, for example, North America, Northern Europe and Patagonia. this field and Steffen & Wu (2011) review modern observational and modelling efforts in this region. One important observable of GIA, but less used and investigated compared to deformation, is the secular gravity change. Redistribution of masses within the Earth as well as on the surface cause changes in the gravity field. Also the vertical land motion/uplift itself induces changes in gravity on the surface of the Earth. Knowledge about this GIA-induced rate of change of gravity,ġ, is important in many aspects, for example (1) for reduction of terrestrial gravity observations to a certain epoch, (2) as ground truth for satellite gravity missions (e.g. Steffen et al. 2009;Müller et al. 2012), and (3) for constraining and tuning GIA models (e.g. Steffen et al. 2014;Van Camp et al. 2017).
Several models of the GIA-induced vertical displacement rate have been published for Fennoscandia (e.g. Ekman 1996;Lambeck et al. 1998;Milne et al. 2004;Ågren & Svensson 2007;Lidberg et al. 2010). Although several observationalġ-results exist (see Table 1), noġ-model for Fennoscandia has been published so far. This is primarily because terrestrial gravity measurements are time consuming and need on-site manpower. Absolute gravity (AG) observations are consequently more expensive than most other geodetic observations. In addition, combination of gravity measurements is challenging due to sensor-affecting incidents and local gravity effects that may mask the secular trend due to GIA.
The first systematic observations of the GIA-induced gravity change were repeated relative gravity observations along the socalled Fennoscandian land uplift gravity lines. They consist of four east-west high precision relative gravity profiles, approximately following the latitudes 65 • , 63 • , 61 • and 56 • (see Fig. 1). Measurements along the Finnish part of the 63 • line started in 1966 followed by the rest of the lines from the mid-1970s (Kiviniemi 1974;Ekman & Mäkinen 1996;Mäkinen et al. 2005). The work with the Fennoscandian land uplift gravity lines was initiated and coordinated by the Nordic Geodetic Commission (NKG).
From the late 1980s the relative gravity observations along the uplift lines have been complemented and gradually succeeded by repeated AG observations. In 1988, the Finnish Geodetic Institute (FGI) started this work using a free-fall absolute gravimeter JILAgtype (Torge et al. 1987), JILAg#5. This gravimeter was mainly used in Finland but also at some stations in the other Scandinavian and especially the Baltic countries. During the 1990s, the JILAg measurements were complemented by observations with its successor, the FG5 (Niebauer et al. 1995). These first FG5 campaigns were performed by the National Oceanic and Atmospheric Administration (NOAA), USA, in 1993 and1995. Further FG5 campaigns were conducted by the Bundesamt für Kartographie und Geodäsie (BKG) in 1993(BKG) in , 1995(BKG) in , 1998(BKG) in and 2003 on 15 stations distributed in the uplift area. In 2003-2008 comprehensive campaigning was carried out with an FG5 instrument by the Leibniz Universität Hannover (LUH), Germany. During that time also the FGI, the Norwegian University of Life Sciences (NMBU) and Lantmäteriet (the Swedish mapping, cadastral and land registration authority) invested in FG5 gravimeters and started with repeated AG observations. In 2008 the Technical University of Denmark (DTU) started making repeated measurements with their A10 absolute gravimeter (Micro-g La-Coste 2008). Today there are 688 AG observations on 59 stations in the region (Fig. 1), most of them co-located with Global Navigation Satellite Systems (GNSS) reference stations. Two of the stations, Metsähovi (since 1994;Virtanen 2006) and Onsala (since 2008), also house superconducting gravimeters (SG). As in the case of the land uplift lines, the work with absolute observations was and is coordinated by the NKG.
Only parts of the Fennoscandian repeated AG observations have hitherto been published. Gitlein (2009) published the results from the BKG, NOAA and LUH campaigns in 1993-2008 with focus on the LUH data. Ophaug et al. (2016) published all FG5 data on the Norwegian stations. Selected observations have been included in special studies, for example, to address theġ/ḣ ratio (Pettersen 2011). Table 1 gives an overview of publications addressing different parts of the whole data set. This includes some unpublished reports and poster presentations since they, in the absence of better references, sometimes have been cited in the literature.
Besides Fennoscandia, GIA-induced surface deformation and gravity changes can also be observed in North America. Compared to Fennoscandia, both the signal strength and the geographical extent are larger. AG time-series from North America was analysed by, for example, Larson & van Dam (2000) and Lambert et al. (2001Lambert et al. ( , 2006Lambert et al. ( , 2013a. A map of gravity rate of change in North America, but mainly based on relative gravity measurements, was published by Pagiatakis & Salib (2003). In their study, they readjusted the primary Canadian Gravity Standardization Network using relative gravity measurements spanning over 40 yr. The gravity rate of change was introduced as an unknown in the observation equation and AG measurements were used as weighted constraints in the (least squares) adjustment.
The relation betweenġ and the vertical displacement rate of the crust,ḣ, is also an important observable since it (1) is affected by both the vertical movement itself as well as by mass changes beneath the surface and therefore contains information on the underlying geophysics and geodynamics (e.g. Ekman & Mäkinen 1996
In addition, a trustworthy relation betweenġ andḣ also allows us to make transformations between, and combine, the two observables.
As mentioned, in regions like Antarctica and Greenland, the ratio betweenġ andḣ has been used for separating the present-day icemass change signal from the GIA signal, the latter induced by historical ice mass variations (Wahr et al. 1995;James & Ivins 1998;Fang & Hager 2001;Purcell et al. 2011;Memin et al. 2012). From an analytical study with a GIA model for Greenland and Antarctica, Wahr et al. (1995) found the viscous part of the ratio to be ∼−0.154 μGal mm −1 . Using the ice model ICE-3G, James & Ivins (1998) predictedġ andḣ for Antarctica, and found their ratio to be ∼−0.16 μGal mm −1 . These predictions are based on modelling and are difficult to verify by observations, because gravity change due to present-day ice mass variation is superimposed by the viscous gravity signal.
In North America and Fennoscandia the situation is different. Here, the signal is strongly dominated by the past GIA signal and the ice-free conditions make it possible to conduct repeated measurements of both gravity and height changes. Table 2 summarizes published ratios based on observations in these regions. Olsson et al. (2015) investigated the geophysical relation betweeṅ g andḣ in previously glaciated areas (like Fennoscandia and Laurentia) using a GIA model, similar to the one described in Section 2.5, and found that (1) their ratio varies in the spectral domain and is smaller (less negative) in the lower part of the spectrum, implying that for a region where the GIA signal is smooth and has a large geographical extent (Laurentia) the ratio is expected to be smaller than for a region where higher degrees of the spectrum dominate the signal (Fennoscandia), (2) the borderline between the uplift area and the forebulge area (zero line) forġ andḣ does not exactly coincide, which affects their ratio especially where the signal is small, (3) within Fennoscandia the ratio varies laterally in such a way that for practical applications these variations can be neglected, (4) local effects, such as direct attraction and short wavelength elastic deformation from present-day GIA-induced sea level variations do not significantly affect the ratio other than in extreme cases (when the station in question is located very close to and high above the sea).
These conclusions imply that for Fennoscandia it is a reasonable assumption to estimate a single linear relation betweenġ andḣ for the entire region.
For the first time we present estimated gravity rates of change based on all repeated gravity observations, spanning over three decades, in the Fennoscandian land uplift area. All observations are provided and described in detail. Estimatedġ values are assessed by the geophysical relation betweenġ andḣ, found from GIAmodelling, and the uncertainties in these relations are discussed. We also suggest aġ model covering the whole area, based on the state of the art land uplift model and the geophysical relation betweenġ andḣ.
In Section 2, we describe the AG data set, how data from different sources have been processed and merged, known error sources, and uncertainty estimates. We also introduce land uplift data sets and a geophysical GIA model for comparison to our observational gravity rate of change. In Section 3, we estimate observational values oḟ g and compare it with a semi-empirical land uplift model as well as a pure GIA model. The relation betweenġ andḣ is estimated and studied in Section 4 and it is further used for constructing ȧ g-model, covering the whole area. This is followed by a discussion of the results and a summary of conclusions. Detailed information about the stations and all observations are provided as Supporting Information (Tables S2 and S4).

The AG stations
We have used data from 59 stations in the region where repeated AG observations have been conducted ( Fig. 1 and Table S2). Steffen et al. (2012) studied optimal locations for AG observations and concluded that, except for the northwestern part of Russia, these  Ekman & Mäkinen (1996) Fennoscandia −0.16 ± 0.05 to −0.18 ± 0.06 Ekman & Mäkinen (1996) revisited, this time with more observations ofġ as well asḣ (including GNSS) Mäkinen et al. (2005) Fennoscandia −0.163 ± 0.02 Four years of annual AG-observations on eight stations.ḣ from GNSS (Lidberg et al. 2007). For the different stations, the ratios vary between −0.114 ± 0.031 and −0.232 ± 0.059 μGal mm −1 The viscous part of the ratio in an area affected by present-day ice mass change. Different ratios depending on how the present-day signal is corrected for Sato et al. (2012) stations form a complete and adequate network for providing constraints for the study of GIA parameters.
Most of the stations are co-located with permanent GNSS reference stations in the so-called BIFROST (Baseline Inference for Fennoscandian Rebound Observations, Sea level and Tectonics) network (see e.g. Johansson et al. 2002;Lidberg et al. 2010). Many of these stations have GNSS time-series spanning more than 20 yr. The AG stations typically consist of a concrete pillar mounted directly on solid bedrock, housed in the same building as the GNSS station ( Fig. 2). Some of the stations (e.g. Metsähovi, Mårtsbo, Onsala and Trysil) have two or more pillars and are therefore suitable for comparisons of instruments by means of simultaneous observations. Some stations are not dedicated AG stations but rather housed in public, stable buildings.
Metsähovi (MET) and Onsala (ONS) are geodetic fundamental stations in the sense that they host instrumentation for a large variety of observational techniques like AG, superconducting gravity, very long baseline interferometry, satellite laser ranging (MET), tide gauge (ONS) and monitoring of local hydrology.
In addition to the stations discussed above, some hundred other AG stations have also been observed with absolute gravimeters (typically A10 gravimeters). These are more simple stations like a benchmark mounted in a rock, stairs or similar. The purpose of these observations was not to study GIA or other geophysical processes and phenomena but rather to serve as datum points for national gravity reference systems. These stations and observations are therefore not treated here.

The AG observations
During the time period 1988-2015, 688 repeated AG observations were conducted at the stations described above. One observation is here understood to be the mean of a large number of free fall experiments (drops). The drops are normally executed during a time period of ∼12-48 hr and grouped in sets of ∼50-100 drops. If there was more than one consecutive set-up of the instrument (e.g. with different orientations) at one visit of the station, then the results of the different set-ups are merged to one observation.
Many different organizations have contributed with observations (Table 3). Each organization initially processed their own data. One representative for each country (Table S1) then collected, and in some cases reprocessed, all data from stations in his/her country, respectively. Data from all participating countries have then been merged into one database (Table S4).
The bulk of the observations was collected using FG5 gravimeters (Niebauer et al. 1995). These data were processed using the 'g' software (Micro-g LaCoste 2012) with final International Earth Rotation and Reference System Service (IERS) polar coordinates, calibrated rubidium frequencies, and standard modelling of gravitational effects due to earth tides, ocean loading and varying atmospheric pressure, as implemented in the 'g' software [for details concerning e.g. ocean tide loading (OTL) models, see Supporting Information]. There have been attempts to perform a refined modelling of the gravitational effect due to ocean loading, non-tidal ocean loading and global hydrology (Ophaug et al. 2016), as well as the atmosphere (Gitlein 2009;Ophaug et al. 2016). The general conclusions of these studies are that refined modelling does not give any significant improvement with respect to the gravity trends on average. In addition, the lack of corrections for local hydrology, which could dominate the gravity rate at a specific site, is identified as an important issue for further research (see e.g. Van Camp et al. 2016b). Thus, until the refined modelling improves and the effect of local hydrology can be embedded, we stick with the standard processing scheme in this work.
All data are presented in the zero tide system. Some of the first observations (e.g. IMGC from 1976) were originally in the mean tide system but have been reprocessed to the zero tide system (Haller & Ekman 1988), following the IAG resolution from 1983 (IAG 1984).
Details about the data and data processing are given in the Supporting Information.

Instrumental biases
AG observations are in general sensitive to instrumental biases (or offsets). In order to detect such biases the International Bureau of Weights and Measures (BIPM) organized international comparisons of absolute gravimeters on a regular basis between 1981 and Downloaded from https://academic.oup.com/gji/article-abstract/217/2/1141/5304614 by guest on 15 March 2020   (Table 4). For each comparison a Comparison Reference Value (CRV) is determined by the participating instruments and individual instrumental biases relative to the CRV are determined for each instrument. Table 5 summarizes the results for the instruments relevant for this work. The methods for determining CRVs, biases and especially uncertainties have varied through the years. In later years the officially given uncertainties include a systematic component for each instrument which is, in general, not the case for the results of the early comparisons. In order to make the numbers in Table 5 comparable to each other, we have chosen to provide the 2σ uncertainty from the adjustment/estimation of the instrumental biases. Also, since the sign of the reported offset/DoE (degree of equivalence) has changed over the years all values have been converted to DoE (Instrument#XXX-CRV). Table 5 shows that the participating instruments normally agree with the CRV within the uncertainty limits. In a few cases the estimated bias is larger than two times the standard uncertainty and in only two cases (JILAg#5 2001 and FG5#220 2015) the bias is larger than three times the standard uncertainty. As mentioned before, the uncertainties given in Table 5 are taken as two times the standard uncertainty of the estimated biases (from the adjustment), which is how the uncertainties for the first comparisons were reported. The modern way of reporting expanded total uncertainty was not reproducible for these old results. In order to make all results in Table 5 comparable we had to choose this way of giving the uncertainty. From 2009 the officially published uncertainties are found directly from the expanded total measurement uncertainty reported for each instrument combined with the uncertainty of the estimated CRV. This method results in larger uncertainty estimates than those in Table 5, and based on these, none of the instruments in Table 5 was reported to have significant biases compared to the CRV.
Our study includes data from one JILAg instrument (#5). Table 5 indicates that it might have been biased and that the bias might have changed but these results are not significant. Other institutions have also reported on biases for their JILAg instruments. For JILAg#3 of the Hannover group (LUH), an obtained discrepancy to the FG5#220 (LUH) of +9.0 μGal indicates a significant longterm offset between the measuring levels of the two gravimeters (Timmen et al. 2011). Similar discrepancies have also been reported by Torge et al. (1999) when comparing measurements from FG5#101 (BKG) and JILAg#3 performed in the years 1994-1997. These comparisons showed a discrepancy varying between +8.1 and +9.4 μGal. It is interesting that the same long-term bias of +9 Gal was also determined for the JILAg#6 gravimeter (see Pálinkáš Downloaded from https://academic.oup.com/gji/article-abstract/217/2/1141/5304614 by guest on 15 March 2020 Table 4. Overview of official international (ICAG), European (ECAG) and regional (EURAMET) Comparisons of absolute gravimeters, held in Sèvres (France) until the year 2015, Walferdange and Belval (Luxembourg). The standard deviations (1σ ) of all participating instruments' degrees of equivalences (DoEs) are also given for each campaign.   Liard et al. (2003). Some hints are given in Wilmes et al. (2003) that similar offsets may exist in other JILA gravimeters with respect to FG5 meters. Besides these long-term biases, varying biases valid for shorter periods may exists for gravimeters and depend on the experts who re-adjust the instruments from time to time. One major disadvantage of the JILAg design compared with the FG5 instruments is the high sensitivity of JILAg meters to floor tilts occurring during each drop which is triggered by the dropping mechanism similar in all drops of a measuring set. Because the interferometer design is not following the Abbe rule like it is realized in the FG5 instruments (reference and test prism in one vertical line), tilt coupling errors of some microgal could occur at locations with weak floor conditions. That introduces a systematic error in the station determination by the gravimeter and can only be detected by a new set-up of the meter with another orientation. For FG5 gravimeters, the effect has been minimized, see Niebauer et al. (1995).
By assessing local comparisons between some of the instruments relevant for this study also Pettersen et al. (2010) conclude that data from these instruments reveal no systematic biases, but occasional shifts from 1 yr to another are noted. This was also found by Olsson et al. (2016). They showed that time-series from the FG5#233 gravimeter indicated a jump in 2010. The jump occurred during a service of the instrument by the manufacturer, but no real explanation has been found, yet. The effects of that jump could be reduced by introducing a small correction based on the results from the international comparisons.
Based on the results above, data from FG5#233 have been corrected for the suspected jump in this study (see further Section 3) but no other biases between instruments have been considered.

The NKG2016LU land uplift model
NKG2016LU is a successor of the empirical land uplift model NKG2005LU, which has been the official standard model for geodetic land uplift applications in the Nordic countries for the last decade.   Lambeck et al. (1998) was used. For a thorough description of NKG2005LU, seeÅgren & Svensson (2007) and Vestøl (2007).
In 2016, the NKG Working Group on Geoid and Height Systems released the land uplift model NKG2016LU, which is now called semi-empirical in order to emphasize that it, in addition to observations, also includes a GIA modelling component. Notable differences to NKG2005LU include (1) longer GNSS time-series. Vertical velocities from the BIFROST 2015/16 calculation, processed in GAMIT/GLOBK and finalized in 2016 March 1. This is an updated version of Kierulf et al. (2014).
(2) omission of tide gauge data. Spatial and especially temporal variations in the rate of change of mean sea level (e.g. accelerating sea level rise during the last decades) prompted the decision not to include tide gauge data in NKG2016LU.
NKG2016LU comes in two versions, NKG2016LU lev and NKG2016LU abs. NKG2016LU lev is the land uplift as measured with repeated levelling, that is relative to the geoid. NKG2016LU abs (see Fig. 1) is the absolute land uplift in ITRF2008 as observed by GNSS. In the observation points, the mean difference between the BIFROST GNSS solution and the final NKG2016LU abs model is 0.02 ± 0.42 (1σ ) mm yr −1 , which corresponds to ∼−0.003 ± 0.07 (1σ ) μGal yr −1 (see below). As NKG2016LU is given in the same reference frame as the BIFROST GNSS solution, but also includes levelling data, and gives a trustworthy interpolation between the observation points (and thus a value for all gravity points and any other point), we take it rather than the GNSS solution itself as a reference model.
For conversion of the NKG2016LU abs land uplift to gravity rate of change we use the factor C = −0.163 μGal mm from −1 the modelled linear relatioṅ found by Olsson et al. (2015), valid for 1-D geophysical GIA models (normal mode approach) in Fennoscandia. The uncertainty of the factor has been estimated to u(C) ∼ ±0.016 μGal mm −1 (Ophaug et al. 2016).

The geophysical GIA model ICE-6G(VM5a)
In addition to using the state of the art Fennoscandian land uplift model NKG2016LU (based on land uplift observations),ġ is also predicted by means of a standard geophysical GIA model, namely ICE6-G(VM5a), which is widely used throughout the world as a reference for land uplift and gravity observations. The GIA model is based on the viscoelastic normal-mode method, pseudo-spectral approach (Mitrovica et al. 1994;Mitrovica & Milne 1998), with an iterative procedure in the spectral domain and spherical harmonic expansion truncated at degree 192 (Steffen & Kaufmann 2005) and applied using the software ICEAGE (Kaufmann 2004). The ice history is according to the ice model ICE-6G C and earth rheology according to earth model VM5a (Argus et al. 2014;Peltier et al. 2015). The direct attraction term (from present day, GIA-induced sea level variations) in the Green's function for gravity was omitted, following the recommendations from Olsson et al. (2012).

E S T I M AT I O N O F G R AV I T Y T R E N D S F RO M O B S E RVAT I O N S
From the repeated AG observations we estimateġ at all stations with more than two observations and a time span longer than 2 yr. For comparison, we constructed two different data sets (I and II) based on the observations listed in Table S4. Dataset I includes all observations as they are and Dataset II is refined in such way that observations and stations with large uncertainties and suspected errors are removed (see below). These estimated gravity trends are then    compared with NKG2016LU abs (Section 2.4) and a geophysical GIA model (Section 2.5), shown in Table 6. Dataset I consists of all AG observations as they are listed in Table  S4. The gravity rate of change,ġ, and a reference gravity value, g 0 , in the reference epoch, T 0 (mean epoch of all observations), are estimated for each station, i, by means of weighted least squares adjustment (WLSA) with the observation equations where g i j obs is one gravity observation at station i at epoch T j . The observations are weighted with 1/σ 2 tot , where σ tot is the total standard uncertainty as given in Table S4.
In Dataset II only FG5 observations are used, that is IMGC, GABL, JILAg, and A10 observations are omitted and only stations with 5 or more observations spanning over at least 5 yr are considered.
The omission of other absolute observations than those made with FG5 is motivated by the fact that FG5 instruments have a lower observational uncertainty than the other types of instruments. Especially, the internal consistency with this group of AGs is high, which is crucial here when repeatability is more important than the absolute level. Using only one type of instrument decreases the risk of introducing (unknown) offsets between instruments. Since the observations with the omitted instruments in general are concentrated to the earliest part of the time-series (except A10), any offsets would greatly impact trend estimates. Except for the JILAg instrument the omitted instruments have contributed with relatively few observations. In Finland, JILAg#5 was heavily used during the 1990s and early 2000s, especially at the METS station. Fig. 4 shows all observations at METS. Up to 2003 these observations are almost exclusively JI-LAg type, and after 2003 they are only FG5 type. Three different estimates ofġ at METS are shown in Fig. 4; one using only JI-LAg observations (-0.55 ± 0.18 μGal yr −1 ), one using only FG5 (-0.35 ± 0.06 μGal yr −1 ) and one using all available observations (−0.75 ± 0.05 μGal yr −1 ). Using all observations, the estimatedġ agrees very well with the rate predicted by the NKG2016LU abs model. The FG5 trend differs significantly from the trend based on all observations and one reason could be a possible offset between the JILAg#5 and FG5 instruments. Introducing this offset as an unknown in the observation equation (eq. 2) gives an estimate of the offset between JILAg#5 and FG5 of 7.74±0.78 μGal and, at METS,ġ = −0.41±0.06 μGal yr −1 andġ/ḣ = −0.092±0.013 μGal mm −1 . The results from international comparisons (Table 5) indicate that the bias for JILAg#5 might have changed over the years, but these numbers are not significant and the bias for JI-LAg#5 is therefore not taken into account in this work (applies to Dataset I).
Since the FG5 trend (as well as the JILAg trend and the trend corrected for an offset) differs significantly from the land uplift model and because of the problem with the suspected offset between the JILAg and FG5 observations, the METS station is excluded from Dataset II. HONC, TRDA and TROM (Ophaug et al. 2016) and VAAA and KEVO have been pointed out to have gravity trends induced by multiple overlapping processes thus hiding the GIA signal. They are therefore also omitted from Dataset II.
In Dataset II, the shift identified in the FG5#233 time-series (see Section 2.3) is corrected according to method 3c in Olsson et al. (2016), that is, with the DoE reported from the international comparisons (Table 5).
The adjustment of the data in Dataset II is conducted the same way as for Dataset I (eq. 2). Two observations (TRYB 2008.254, MARA 2013) are identified as outliers (deviate more than 3σ tot from the estimated trendline) and are therefore removed. At stations with observations on more than one pillar (METS, MARA, ONSA, TRYC) the observations on the individual pillars have been merged to one, in the adjustment, by assuming the samė g-value on all pillars and estimating an additional parameter for gravity difference between the pillars.
The difference between Dataset I and II (Table 6) can be explained in different ways for different countries. In Finland and the Baltic countries the difference is primarily because of the exclusion of the JILAg data, in Denmark the exclusion of A10 data and in Sweden because of the correction for the identified shift in the FG5#233 time-series. g GIA in Table 6 represents the global GIA model ICE-6G(VM5a), described in Section 2.5. It is included here to show how such a model performs compared to observational data. Table 7 shows the difference between gravity change predicted using the empirical model, NKG2016LU, and the other predictions/estimates ofġ. Foṙ g GIA the standard deviation is smaller compared to the observed rates but on average the AG-observations fit better with the empirical model, that is other types of geodetic observations in the area. g GIA is not specifically tuned to Fennoscandia and modern GIA observations there and systematically underestimates the gravity change (is less negative) compared to both AG-observations and NKG2016LU. Below NKG2016LU will be used as the reference model.

T H E R E L AT I O N B E T W E E Nġ A N Dḣ
For evaluation of the geophysical relation betweenġ andḣ (eq. 1), we apply both WLSA and WODR methods to estimateġ 0 and C iṅ from observations. The first method allows errors in the observations (ġ) to be taken into account, while the latter considers also errors in the regressor (ḣ). In eq. (3),ġ i isġ from Table 6 for station i andḣ i is the corresponding land uplift value from NKG2016LU abs. The standardized residuals, given in Table 6, arē from the WLSA solution. They indicate if the residual between the estimatedġ-value for the station in question (ġ i ) and the trend line (ġ =ġ 0 + C ·ḣ) is smaller (<1) or larger (>1) than the estimated standard uncertainty for thatġ-value. For exampleε i > 3 indicates that the estimatedġ i -value deviates more than 3σ from the trend line.
Using WODR, the minimization problem is defined as (Boggs et al. 1992) where w i and ε i are the weight and the residual ofġ i , and w δ i and δ i are the weight and the residual ofḣ i . For bothġ andḣ the weights were set equal the inverse of the squared standard error. We used the ODR-package (Boggs et al. 1992) of the Python library SciPy to solve the minimization problem defined in eq. (5). Using WODR we circumvent a systematic bias that is introduced if we use WLSA for line fitting when there is uncertainty in the predictor (Pitkänen et al. 2016). Because WLSA aims to minimize the vertical distance between data points and the fitting line, a larger horizontal spread of the predictor will cause the fitting line to accommodate  by sloping (or attenuating) towards zero. This mechanism is known as the attenuation or regression dilution bias (e.g. Hutcheon et al. 2010;Van Camp et al. 2016a). By contrast, WODR is an example of a bivariate regression technique which takes uncertainties of both outcome and predictor into account, and minimizes the shortest distance between data points and the vertical line. As such, the mechanism causing the regression dilution bias never occurs. In Fig. 5, all the estimated gravity rates from Table 6 are plotted against their corresponding land uplift value (NKG2016LU abs), for both data sets. Also the estimated linear relations (WODR) as well as the geophysical relation (eq. 1) are plotted. The bottom panel of Fig. 5 shows the trend for Dataset II forced through the origin, that isġ 0 = 0. It is clear from eq. (1) that the GIA-model predicts a small deviation ofġ 0 from 0. Still, most of earlier studies of this relation (Table 2) have assumedġ 0 = 0 and therefore we include that here for comparison. Collilieux et al. (2014) and Mazzotti et al. (2011), for example, useġ 0 for evaluation of systematic errors inḣ based on the assumption thatġ 0 = 0 would indicate systematic errors in the scale and centre of mass of the GNSS reference frame. It should be noticed that also systematic errors in the gravity rates would result in offsets of the trend line. The identified shift in the time-series of FG5#233 (see Section 3) caused systematically lower estimates of the gravity rates in Sweden. Dataset II is corrected for that shift but Fig. 6 shows the trend line WLSA for Dataset II without this correction.
Although NKG2016LU is our preferred solution forḣ, we have also fit eq. (3) to Dataset II combined withḣ derived from the BIFROST GNSS observations. This implies that the weights forḣ in the WODR algorithm vary between the stations, in contrast toḣ from NKG2016LU which all have the same weights (0.6164 mm yr −1 ). Note that for four of the stations in Dataset IIḣ from GNSS is not available as they are not a part of the BIFROST network and therefore not included in this solution (see Table 6). Table 8 summarizes the results for different combinations of data sets, sub-sets of stations and estimators. The results indicate that the differences between estimates calculated with WLSA and WODR are small, that is, within one sigma for both Dataset I and II.
The empirical results are well within the 95 per cent confidence interval of the geophysical relation and all the empirical relations are smaller (more negative) than the geophysical. The estimates of C based on Dataset II range from −0.163 to −0.177 μGal mm −1 and agree within the geohysical/modelled value's standard error. The agreement between the solutions indicates that the estimates based on Dataset II are quite robust considering weighting strategy and regression method.

D I S C U S S I O N
We have used the complete data sets to estimate homogenous relations betweenġ andḣ for the region. Of course, ratios between  g andḣ can also be estimated station-wise, but the uncertainties in the observations are (still) too large for this to be meaningful, especially whenḣ andġ → 0. Also the uncertainties of estimatedġ (Table 6) are, in general, large compared to the uncertainties of the land uplift model. This is due to the fact that there are still quite few gravity observations at most stations (≤5 for 40 per cent of the stations) and that there are unmodelled local effects, possibly due to local hydrology or sea level variations (for stations very close to the sea), that may introduce both random and systematic errors in the gravity timeseries (Van Camp et al. 2016b). Van Camp et al. (2005) show that with annual or semiannual AG observations we can expect a standard error of ∼0.1 μGal yr after 15-25 yr. This is in agreement with the uncertainty forġ LU in Table 6 but, due to few observations and shorter time spans, only a few of the observational rates are close to this.
Not only is the uncertainty of the land uplift model still smaller than the uncertainty of the observational AG rates, it is also carefully interpolated (and extrapolated) between the points of observations which allows us to predictġ at any location (important e.g. for reduction of gravity observations in general to a certain epoch). Still, we need to choose a relation betweenḣ andġ in order to convert the land uplift model to gravity. The observational and geophysical relations agree within the uncertainty limits, giving us increased confidence in the latter. This suggests it is safe to use the geophysical (modelled) relation betweenġ andḣ in combination with NKG2016LU abs to convert vertical rates to gravity change. We call this model NKG2016LU gdot and it is consequently defined as NKG2016LU gdot = −0.163 · NKG2016LU abs (μGal yr −1 ) (6) Worth noticing about the NKG2016LU gdot is that it is valid in the whole Fennoscandian land uplift area but not on, or very close to the sea. There the relation betweenġ andḣ is different because of the direct attraction from GIA-induced sea level variations, and depends on the physical location of the station relative to the sea. Olsson et al. (2015) show that for stations located closer to the sea than 10 times the height of the station this effect should be considered and requires a local and rigorous treatment (see e.g. Lysaker et al. 2008;Breili 2009;Olsson et al. 2009;Breili et al. 2010;Olsson et al. 2012).
In Fig. 7, NKG2016LU gdot is plotted together with the difference between this model and the observationalġ-values from Dataset II,ġ I I ( Table 6). The deviations of the observational results from the model are well within the 95 per cent uncertainty level of the estimatedġ (cf. Fig. 5). Close to the land uplift maximum there are three stations (SKEL, LYCK and RATA) where the observed value is larger (more negative) than the model. This is partly explained by Olsson et al. (2016) as a consequence of the introduced correction for the jump in combination with few observations after 2013. The large positive anomaly ofġ at ANDO indicates that the observed gravity change signal is dominated by other processes than GIA, for example, tectonics or varying hydrology (Ophaug et al. 2016). Fig. 8 shows the difference between NKG2016LU gdot and NKG2016LU abs converted toġ using the observational relation,ġ = 0.06 − 0.172ḣ (Dataset II, WODR). The difference in Fennoscandia is smaller than ±0.05 μGal yr −1 which means that for 20 yr of epoch reduction, using one or the other model, the difference will be smaller than 1 μGal.
Finally, we make a comparison of the results in Table 6 with the results from land uplift gravity lines (Fig. 1). Mäkinen et al. (2005) presented theġ difference between VAGA and KRAM along the western part of the 63 • line and between VAAB and JOEN along eastern part. On the western part VAGA has been excluded from Dataset II because of too few observations so here the results from Dataset I are used. Table 9 summarizes this comparison and the conclusion is that within the uncertainties all results agree. Since AG observations give different trends at VAAA and VAAB (Table 6) this also confirms that the AG trend at VAAA probably consists of more than the GIA signal.

C O N C L U S I O N S
For the first time, all repeated AG observations  in the Fennoscandian land uplift area were compiled and presented. This means 688 observations at 59 stations across the region. Ten different organizations have contributed with data spanning for more than three decades. The primary application of the observations is to study the GIA-induced gravity rate of change,ġ. This study also clearly demonstrates the possibility to determine theġ/ḣ ratio with sufficient precision to validate corresponding results from, for example, GIA models and to be used for converting absolute land uplift values,ḣ, to surface gravity change,ġ.
For all stations with more than two observations and a time span longer than 2 yr,ġ was estimated and compared to predicted values. Two data sets were derived; Dataset I corresponds to all original data and Dataset II is modified with the intention to reduce effects from known or possible systematic or gross errors and includes only FG5 observations.ġ was also determined at all AG stations using (i) the semi-empirical land uplift model NKG2016LU and (ii) Downloaded from https://academic.oup.com/gji/article-abstract/217/2/1141/5304614 by guest on 15 March 2020 Figure 8. Difference betweenġ-models: NKG2016LU gdot−ġ obs , whereġ obs is NKG2016LU abs converted toġ using the relation-based geodetic AG observations,ġ = 0.06 − 0.172ḣ (Dataset II, WODR). NKG2016LU was chosen as reference model and the mean differences between this model and the empirical values are not significantly deviating from zero (0.03 ± 0.67 and 0.03 ± 0.28 μGal yr −1 for Dataset I and Dataset II, respectively). The standard deviation for the difference between the reference model and the GIA model is smaller, 0.11 μGal yr −1 , but the GIA model systematically underestimates the gravity change.
A linear relation,ġ =ġ 0 + C ·ḣ, valid for the entire region, betweenġ andḣ (NKG2016LU) was determined by means of WLSA and WODR for each data set. The difference between estimates calculated with WLSA and WODR is small, that is within 1 σ . Dataset II results in smaller standard deviations and estimates of C from Dataset II range from −0.163 to −0.177 μGal mm −1 . Estimates of the constant part are not significantly different from zero. All empirical results are smaller than, and well within the 95 per cent confidence interval of, the geophysical relationġ = 0.03 − 0.163ḣ. This implies that using the geophysical relation is a reasonable choice. Just using the simple ratioġ/ḣ = −0.163 (μGal mm −1 ) will differ from using the full relation only by 0.03 μGal yr −1 , that is <1 μGal over 30 yr, and may be a reasonable choice for practical applications. This also exactly coincides with estimates of C (Dataset II) whenġ 0 is assumed to be zero.
The uncertainty ofġ estimated from observations at the gravity stations is relatively high and inhomogeneous (∼0.1−0.6 μGal yr −1 ) when compared to the lower and more homogenous uncertainty obtained by predictingġ from land uplift observations by means of the land uplift model NKG2016LU abs (0.1-0.2 μGal yr −1 ). In addition, the gravity observations are geographically limited to a few discrete points while the land uplift model comes as an interpolated surface (grid) covering the entire region. At present, we therefore recommend using the latter, which we call NKG2016LU gdot (= −0.163·NKG2016LU abs), as the most reliable and suitable method to predict the GIA-induced gravity change in Fennoscandia. A gridded version of NKG2016LU gdot can be downloaded at https://www.lantmateriet.se/en/maps-and-geographic-inf ormati on/GPS-och-geodetisk-matning/Ref erenssystem/Landhojning/.
Continuation of the AG observations at the stations already established is important in order to decrease the uncertainty and enable more accurate determinations of the relation to the land uplift. This will also improve the possibilities to descriminate the GIA signal from other environmental signals.

A C K N O W L E D G E M E N T S
We would like to express our thanks to all organizations and all people that contributed with absolute gravity observations over the years.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Table S1. Providers of the data in Tables S2 and S4.  Table S2. Absolute gravity stations in Fennoscandia. VGG is the vertical gravity gradient 1.2 m above the floor, n obs the number of absolute gravity observations and T the time span between first and last observation. Table S3. Parameters used to reduce JILAg observations on Estonian stations to reference height 1.20 m. Table S4. Complete list of observations. σ set is the set scatter, i.e. the standard deviation of the set values, and σ tot is the total uncertainty taken into account also the uncertainty from possible instrumental systematic effects or biases (Niebauer et al. 1995). JOEN, SODA, VAAA, VAAB observations are at 1.00 m reference height (approximately mean observation height for stations with both JILAg and FG5 observations), the rest at 1.20 m. Different reference heights at different stations only affect the absolute gravity values and not estimates ofġ. NaN means that data were not available for the authors.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.