1 Introduction

The ongoing Fennoscandian postglacial rebound due to the last glaciation has been investigated in many ways and from different perspectives (e.g. Ekman 1996, 2009; Vestøl 2006; Steffen and Wu 2011). This rebound is most prominent in the Nordic and Baltic countries, but also affects the northern parts of the Netherlands, Germany and Poland as well as north-western Russia. The major visible consequence is the land uplift which has its maximum of slightly more than 1 cm/a near the city of Umeå in Sweden (Steffen and Wu 2011), and leads to about 700 ha of new land in Sweden and Finland emerging from the sea every year (Poutanen and Steffen 2014).

In geodesy, knowledge of the uplift is important since satellite techniques operate in global reference frames, in which this deformation is directly observable. Whenever we transform data referring to such global reference frames into national Earth-fixed reference frames, for instance European Terrestrial Reference System 1989 (ETRS 89) realizations, we must account for land uplift. At the same time, the Continuously Operating Reference Station (CORS) tracking of Global Navigation Satellite System (GNSS) satellites enables accurate measurement of it.

To measure the land uplift in Fennoscandia is one of the major goals of the project BIFROST (Baseline Inferences for Fennoscandian Rebound Observations, Sea Level, and Tectonics). This project, initiated already in 1993, presently uses more than 200 CORS in the Nordic, Baltic and adjacent countries. Results have been published in several publications, e.g. Scherneck et al. (1998, 2003), Milne et al. (2001), Johansson et al. (2002), Mäkinen et al. (2003), Lidberg et al. (2007, 2010) and Kierulf et al. (2014). In this study, we use the latest reprocessing (Kierulf et al. in preparation).

There are different ways to generate land uplift models, but two methods dominate. The first uses geodetic observations (such as tide gauges, precise spirit levelling and/or GNSS) to calculate a model by some estimation method, such as Least Square Collocation (LSC). This method was successfully applied by, for instance, Danielsen (2001) and Vestøl (2006) and gives a strictly empirical model. The second method computes the model in a geophysical meaningful way based on the combination of an Earth model and geologically constrained ice thickness history, selected so that the model best fits selected land uplift observations such as GNSS-derived velocities and/or geological relative sea-level (RSL) changes. This method yields a so-called GIA model. Many GIA models have been developed for Fennoscandia, and an overview can be found in Steffen and Wu (2011). Since then, Zhao et al. (2012), van der Wal et al. (2013), Kierulf et al. (2014), Schmidt et al. (2014), Steffen et al. (2014), Nordman et al. (2015), Root et al. (2015) and Simon et al. (2018) have provided new studies and models.

Here, we combine the two approaches, i.e. combine a strictly empirical model with a GIA model, which gives us what we refer to as a semi-empirical model. For the empirical model, we use geodetic precise levelling and GNSS uplift rates from the BIFROST project together in a LSC approach, while we compute the GIA model by testing many ice-history models and Earth model parameters, and find the GIA model that fits both geological and geodetic GNSS observations the best. The idea is that the GIA model will fill in with uplift information in the gaps between the observation points and extend the empirical model outside the area of observations. The resulting semi-empirical NKG2016LU model replaces the old NKG2005LU model, which was computed in a similar way under the umbrella of the NKG (Vestøl 2006; Ågren and Svensson 2007). This approach has many similarities to the method described by e.g. Hill et al. (2010) or Simon et al. (2018). They also do this combination in a least squares adjustment, but use the observation directly when constraining the a priori GIA model. This is principally the same as our alternative one-step approach, discussed in Sect. 3.4

Our semi-empirical model gives the uplift expressed in a global reference frame, in this case the International Terrestrial Reference Frame 2008 (ITRF2008), since the uplift rates from the CORS are calculated and expressed in this frame. We call this “absolute land uplift”, and the corresponding absolute model is NKG2016LU_abs. It makes sense to say that this uplift is relative to the reference ellipsoid. For some purposes, however, we are more interested in uplift relative to the geoid. This is for instance the case when we correct levelling lines to a common epoch, and subsequently we call such an uplift for levelled uplift and the corresponding model is NKG2016LU_lev.

To convert from “absolute” to “levelled” uplift, we need a model of the geoid rise. We obtain this model in the GIA calculation. We expect this GIA-induced geoid rise model to represent the geophysical geoid change in Fennoscandia reasonably well and to have low uncertainty (see further Sect. 2.3 and the introduction to Sect. 3).

The new model is computed based on similar calculation methods as in 2005. However, there are two important differences: (1) the types of observational data employed and (2) the underlying background model of glacial isostatic adjustment (GIA). In 2005, the observational data consisted of precise levelling, uplift rates from Continuously Operating Reference Stations (CORS) and tide gauge uplift rates. For the new model, we have omitted the latter mainly due to influence of climate-change related effects. Furthermore, the GIA model is now calculated within NKG and replaces the model of Lambeck et al. (1998) used for NKG2005LU.

The NKG2016LU models were officially released in 2016 and have not been changed since then. They are useful for the mapping authorities in the Nordic and Baltic countries to reduce current observations to the reference epochs of the national reference frames, but can also be applied for different types of research. We provide the model for open access download in two versions, as absolute uplift in NKG2016LU_abs with respect to the Earth’s centre of mass in ITRF2008, and as levelled uplift in NKG2016LU_lev with respect to the levelling reference surface (the geoid). Both models, as well as the estimated uncertainty file, are provided in several file formats.

The methods for making GIA and empirical models are well described in the above-mentioned papers and will only be briefly discussed in Sect. 2. More details can be found in the Supplementary Material. The focus in this study is to estimate the uncertainty of the final semi-empirical uplift model NKG2016LU (Sect. 3). After that, in the last section, the conclusions are accompanied by a discussion of the method and resulting models and recommendations for future improvements.

2 The calculation of NKG2016LU

An optimal uplift model would use all available constraint information, including observations as well as the results of GIA modelling. In our case, a natural way is to start with the GIA model and use the observations to calculate correlated corrections to this GIA model as signals in a LSC approach. However, this is not what we have done for the officially released NKG2016LU model. Instead, we have used a two-step procedure where we divide this complex computation into two manageable parts, which historically have been performed by different researchers and groups from different countries. In general, geodesists calculate empirical models and geophysicists make GIA models. Another historical motivation has been to keep the empirical model completely clean from any kind of GIA modelling, relying only on geodetic observations.

In the first step, we calculate an empirical model and a GIA model separately, and in the second step, we combine these two models. To be able to explain the rather elaborate computation of the uncertainty for the final combined uplift model in Sect. 3, we first briefly describe the data and models involved as well as how we combine the models.

2.1 The empirical model

Based on observations with known accuracy, we estimate the land uplift via LSC. Important advantages with an empirical approach are the ability to calculate uncertainty estimates for the result and to detect and remove statistical outliers.

We use geodetic precise levelling observations and estimate GNSS uplift rates from a CORS network. The levelling observations and CORS data cover the Nordic and Baltic countries (Fig. 1a). There are also additional CORS data from some stations in northern Central Europe (Fig. 1b).

Fig. 1
figure 1

a The levelling network used for the empirical model. b All CORS used in the calculation with observed uplift written. The symbols indicate the standard deviation: Red squares < 0.2 mm/a, yellow dots: < 0.5 mm/a and green dots > 0.5 mm/a

Full size image

Re-levelled lines contain direct information of the land uplift as they give the change in the height difference between two epochs and consequently of the relative uplift between the points. In addition, closed loops (or polygons) consisting of lines from different epochs also contain land uplift information (Danielsen 2001). The longer the time span between the epochs, the more accurate and valuable is this type of land uplift information.

For all CORSs, daily solutions are available from a new BIFROST reprocessing (Kierulf et al. in preparation). The GNSS data are analysed with GAMIT/GLOBK, and the results are given in ITRF2008. Uplift rates and uncertainties are estimated using the time-series analysis software Cheetah (Bos et al. 2008) including both white noise and power law noise. The spatial correlation between the stations was neglected in the time-series analysis. The accuracy depends very much on the time span of the data (Kierulf et al. 2014). All stations used in this study have at least three years of data. Following the multiple Student’s t test procedure described in Pelzer (1985) and Vestøl (2006), we have rejected seven out of 177 stations and associated uplift rates.

“Appendix” lists the uplift value and the corresponding standard uncertainty for all 177 stations. Following the variance component estimation method described by, for instance, Caspary (1988), Förstner (1979a, b) and Pelzer (1985), we calculate variance components for all the different datasets (see Table S3 in Supplementary Material). The results of this calculation indicate that the uncertainty estimates listed in the “Appendix” should be scaled by a factor of 1.41. The uncertainty estimates for the CORS in Kierulf et al. (in preparation) are in other words generally too optimistic.

Since we do not use a GIA model in the LSC adjustment for the strictly empirical model, we instead solve deterministically for a polynomial trend surface of the 5th degree to describe the main uplift. In addition, we estimate signals in all observations stochastically. In summary, trend plus signal describes the uplift. For more details, see Vestøl (2006). In between the observation points, we have filled in with points in which we have calculated the signals by using the prediction ability in LSC, and then added the polynomial value. From all these single point values, observed or just predicted, we can interpolate a model as shown in Fig. 2.

Fig. 2
figure 2

The empirical model, uplift in mm/a. The block dots are observed or predicted values from which we generate the contour lines using the method “Triangulation with linear interpolation” in the Surfer version 10.7 software (https://www.goldensoftware.com). The figure shows no extrapolated values, so outside the area with observed or predicted values the figure is grey

Full size image

2.2 The GIA model

The GIA model used for NKG2016LU is computed by Steffen et al. (2016) in the NKG Geodynamics Working Group and is called NKG2016GIA_prel0306. It is part of an NKG activity to generate a new GIA model for Fennoscandia based on a thermodynamically coupled ice-sheet model and a laterally varying Earth model. In the first steps of this long-term activity, a laterally homogeneous (i.e. only varying with depth) Earth model is used. Hence, we apply the commonly used viscoelastic normal-mode method (Peltier 1974; Wu 1978) that is integrated in the software ICEAGE (Kaufmann 2004) and generate a set of spherically symmetric (1D), compressible, Maxwell-viscoelastic earth models that differ in lithospheric thickness, upper and lower mantle viscosity. The surface load is made of a selected ice-sheet history model and the corresponding ocean load that is determined via the sea-level equation (Farrell and Clark 1976). A difference to previous GIA modelling studies is that we do not use one ice model for Fennoscandia, the Barents/Kara seas and the British Isles only. We use a high-variance set of 25 relatively low misfit model runs from an archive called GLAC (Tarasov 2013; Tarasov et al. 2012; Root et al. 2015; Nordman et al. 2015) in our analysis. More than 11,000 different GIA models are then compared to vertical rates of the GNSS data (see “Appendix”) and a large dataset of geological and paleontological relative sea-level (RSL) observations for northern Europe to identify the best-fitting GIA model to the observations. The detailed description of the GIA modelling and the selection of our best-fitting GIA model can be found in the Supplementary Material. The resulting preliminary GIA model is applied for NKG2016LU. Further improvements towards a 3D GIA model will be undertaken in the next years.

Our selected GIA model, called NKG2016LU_prel0306, consists of a 160 km lithospheric thickness, an upper-mantle viscosity of 7 × 1020 Pa s, a lower mantle viscosity of 7 × 1022 Pa s, and uses the GLAC ice-history model #71340. The NKG2016LU_prel0306 uplift field for northern Europe is shown in Fig. 3.

Fig. 3
figure 3

The uplift in northern Europe as calculated with the best-fitting GIA model called NKG2016GIA_prel0306. Unit in mm/a

Full size image

The model is considered as best compromise in fitting data of the whole study area (see Supplementary Material). This is of major importance for our land uplift model, as we are interested in a GIA model that fills the large data gaps naturally existing for instance due to the presence of the Baltic Sea. On land, we have very good control with the observation points that are the base of the empirical part of the land uplift model. The GIA model performs well in most parts of Norway except for larger misfits along the Norwegian Sea coast and an area north of Oslo (Fig. 4). We attribute the first to likely other geodynamic (tectonic), yet not fully defined processes, see also Kierulf (2017). The second misfit north of Oslo also appears in the comparison of the strictly empirical model and the GNSS data and therefore deserves further investigation in future.

Fig. 4
figure 4

Difference between the GNSS uplift rates (observations) and the NKG2016GIA_prel0306 model. Red means that the GNSS uplift rates are larger than the GIA model

Full size image

The fit of the GIA model to observations in Denmark and northern Central Europe is not good, i.e. in the forebulge area south-west of the Tornquist Suture zone (see Supplementary Material). This is because we use laterally homogeneous models. It is known from seismic observations that there is a remarkable jump in lithospheric thickness at the Tornquist Suture zone (Gregersen et al. 2002). Best-fitting GIA models to observation data south-west of this zone (mainly in Denmark, Germany and western Poland) indicate a much lower lithospheric thickness of 100 km only and support the seismic findings. The next GIA model generation should therefore incorporate such lateral variation.

Together with the best-fitting GIA model, we also calculate the corresponding geoid change model illustrated in Fig. 5 (using the ICEAGE software package, Kaufmann 2004).

Fig. 5
figure 5

Geoid change of the best-fitting GIA model NKG2016GIA_prel0306. Unit in mm/a

Full size image

2.3 The combined model

The combined model is created in a remove-interpolate-restore approach, very similar to how the old NKG2005LU model was computed; see Ågren and Svensson (2007).

We summarize the technique in three steps:

  • At the observation points (i.e. in the GNSS CORS and levelling nodal benchmarks), we subtract the GIA model land uplift from the corresponding value of the purely empirical model.

  • The resulting reduced observations (i.e. the difference between the two models) are then gridded by LSC (Moritz 1980) assuming a 1st order Gauss–Markov signal covariance function. The statistical parameters for the signal are chosen based on empirical covariance analysis. We use the estimated standard error of the empirical model at each observation point as observation noise. This step gives the residual surface.

  • Finally, we restore the GIA model to the residual surface grid to get the NKG2016LU_abs grid.

In other words, we take out values from the purely empirical model and handle these as real observations with known uncertainties. In the following, we call these observations model observations. The above procedure implies that the purely empirical model is smoothed. How much depends on the signal-to-noise variance ratio and the density of the reduced observations. The above choice of covariance function further means that NKG2016LU_abs will approach the GIA model with increasing spatial separation from the observations.

Mathematically, the remove-interpolate-restore method can be summarized as follows,

$$ \dot{h}_{{{\text{NKG2016LU\_abs}}}}^{\text{grid}} = \dot{h}_{\text{GIA}}^{\text{grid}} + \overbrace {{{\text{LSC}}\left\{ {\dot{h}_{{{\text{empirical\_abs}}}}^{{{\text{obs}} . {\text{points}}}} - \dot{h}_{\text{GIA}}^{{{\text{obs}} . {\text{points}}}} } \right\}}}^{\text{Residual Surface (grid)}}, $$
(1)

where the notation is obvious from the above discussion. The residual surface for grid node (i, j), \( {\text{RS}}_{ij} \), is then found by a combination of filtering and prediction (interpolation) given by the standard LSC formula (e.g. Moritz 1980),

$$ {\text{RS}}_{ij} = C_{ij}^{T} \left( {C + D} \right)^{ - 1} \left\{ {\dot{h}_{{{\text{empirical\_abs}}}}^{{{\text{obs}} . {\text{ points}}}} - \dot{h}_{\text{GIA}}^{{{\text{obs}} . {\text{ points}}}} } \right\} = A_{ij}^{T} \left\{ {\dot{h}_{{{\text{empirical\_abs}}}}^{{{\text{obs}} . {\text{ points}}}} - \dot{h}_{\text{GIA}}^{{{\text{obs}} . {\text{ points}}}} } \right\}. $$
(2)

Here, C is the signal variance–covariance matrix (in the observation points), D is the variance–covariance matrix for the model observations, and \( C_{ij} \) is a signal cross-covariance column vector expressing the covariance between the signal in grid node (i, j) and in all the observation points. The column vector \( A_{ij}^{{}} = C_{ij}^{T} \left( {C + D} \right)^{ - 1} \) is introduced for use in Sect. 3.3. The propagated variance for grid node (i, j) becomes,

$$ \sigma_{ij}^{2} = C_{0} - C_{ij}^{T} \left( {C + D} \right)^{ - 1} C_{ij} , $$
(3)

where \( C_{0} \) is the (homogeneous) signal variance.

We split the error in the reduced observations into noise and signals and interpret the errors in the model observations as noise described by the variance–covariance matrix D, while the errors of the GIA model are considered as signals to estimate, statistically specified through the variance–covariance matrix C and vector \( C_{ij} \). It follows that the residual surface computed by Eq. (2) then is an estimate of the errors of the GIA model. (Theoretically, other systematic non-GIA effects may also occur. If so, they stem from the empirical model and remain in the final model.) In this section, we assume that the signal variance–covariance function is homogeneous and isotropic.

The variance–covariance matrix D is assumed to be diagonal (uncorrelated observations) containing only variances estimated when calculating the purely empirical model. Assuming uncorrelated observations is a shortcut approximation since these observations, the model observations, are in fact values taken from the empirical model and therefore correlated. The impact of this approximation on the result is not crucial, for the uncertainty estimates however, the effect may be significant and correlation has to be accounted for (see Sects. 3.3 and 3.4).

In the next step, NKG2016LU_lev is derived by removing the geoid rise of the NKG2016GIA_prel0306 model (cf. Fig. 5) from NKG2016LU_abs:

$$ \dot{H}_{{{\text{NKG2016LU\_lev}}}}^{\text{grid}} = \dot{h}_{{{\text{NKG2016LU\_abs}}}}^{\text{grid}} - \dot{N}_{{{\text{NKG}}2016{\text{GIA}}\_{\text{prel}}0306}}^{\text{grid}}$$
(4)

As previously mentioned, NKG2016LU_lev describes the land uplift relative to the geoid. Since the GLAC ice model used to compute NKG2016GIA_prel0306 does not contain any contemporary ice melting and consequent sea-level change, the geoid rise must here be interpreted as an equipotential surface that is still rising due to historical ice melting in the past, through GIA. It can thus be used in present-day sea-level studies to correct postglacial land uplift that is due to old historic deglaciations. Another suitable application of the NKG2016LU_lev model is to reduce precise levelling observations to a chosen reference epoch.

We start by examining the empirical model, the GIA model and the reduced observations at the 1111 observation points (GNSS CORS and levelling nodal benchmarks). The corresponding statistics are given in Table 1, and the reduced observations are illustrated in Fig. 6. The reduced observations show a very clear spatial correlation, which should primarily be due to the errors in the GIA model. We also add statistics for the difference between the 172 cleaned GNSS observations and the GIA model (Fig. 4).

Table 1 Statistics for the empirical model, GIA model and the reduced observations (all in the observation points). The last row contains statistics for the difference between the cleaned GNSS observations and the GIA model
Full size table
Fig. 6
figure 6

Reduced observations (differences between the empirical model and the GIA model in the observation points). The scale is given by the 1 mm/a arrow in the Southern Baltic Sea. Red means that the empirical model gives larger uplift than the GIA model

Full size image

It should be noted that we have chosen not to remove and restore the mean value of the differences between empirical and GIA models in Eq. (1), contrary to the method used for NKG2005LU (Ågren and Svensson 2007). In principle, there might be a bias between the two models. Since the empirical model refers to the reference frame ITRF2008, we may consider this bias as the GIA model’s deviation from ITRF2008, and this must be further checked numerically:

It can be seen in Table 1 that the straightforward mean value is − 0.13 mm/a. However, as the observation points are significantly denser in some areas than in others and we consider the observations to be uncorrelated, this is most likely not a good estimate of how close the GIA model is to ITRF2008. There is for instance a dense net of nodal benchmarks in the southern parts of Norway (blue area), which would highly influence the result since these benchmarks are clearly positively correlated. If we instead study the reduced GNSS observations, which are more evenly distributed, then we get the much smaller mean value − 0.03 mm/a (last row of Table 1). In view of this, and of all the uncertainties involved (cf. Sect. 3), we decided to neglect the mean value remove-restore in Eq. (1). The resulting difference between using and not using mean value remove-restore is mainly that the final model will differ by this amount (i.e. − 0.03 mm/a) from ITRF2008 far away from the observations. Near the observation points, this choice will affect the estimated model negligibly.

The next step in the computation of NKG2016LU_abs is to calculate the residual surface by Eq. (2). This requires that we compute the signal variance C0 and the correlation length (half-length), which we do by fitting the chosen analytical function (1st-order Gauss–Markov) to an empirical covariance function, computed using the cleaned GNSS observations (see Sect. 2.1). We here prefer to study only the 172 GNSS differences from the GIA model, and not all the 1111 reduced observations, as the latter stem from a previous collocation step and are consequently filtered and most likely correlated. Using the reduced observation to find the signal covariance would lead to difficulties in separating the signal correlation from the observation correlation. The empirical covariance function is computed in the way described by Sanso and Sideris (2013, Section 5.8). We couple two and two differences in all combinations and group the couples in classes depending on the distance \( \psi \) between the points in the couple. In all classes, we calculate the empirical covariances as

$$ {\text{COV}}\left( \psi \right) = \sum\limits_{(P,Q)} {\left( {\dot{h}_{{{\text{GNSS}},P}} - \dot{h}_{{{\text{GIA}},P}} - \overline{{\dot{h}_{\text{GNSS}} - \dot{h}_{\text{GIA}} }} } \right)\left( {\dot{h}_{{{\text{GNSS}},Q}} - \dot{h}_{{{\text{GIA}},Q}} - \overline{{\dot{h}_{\text{GNSS}} - \dot{h}_{\text{GIA}} }} } \right)} , $$
(5)

where the summation extends over all pairs of points P and Q in the interval \( \psi - \Delta < \psi_{PQ} < \psi + \Delta \). Notice that, since the mean value \( \overline{{\dot{h}_{\text{GNSS}} - \dot{h}_{\text{GIA}} }} \) is so small (0.03 mm/a; see Table 1), it does not matter in the present case whether it is subtracted or not.

It is assumed that the GNSS observations are uncorrelated, which means that the GNSS noise only affects the empirical variance, not the covariances. The signal variance at zero distance, i.e. C0, is therefore computed by reducing the empirical variance (containing both signal and noise) by the mean variance of the GNSS noise as follows,

$$ C_{0} = \frac{1}{172}\sum\limits_{i = 1}^{172} {\left( {\dot{h}_{{{\text{GNSS}},i}}^{{}} - \dot{h}_{{{\text{GIA,}}i}}^{{}} } \right)^{2} - } \frac{1}{172}\sum\limits_{i = 1}^{172} {\sigma_{{{\text{GNSS-}}\,{\text{noise}},i}}^{2} } = 0.13\,{\text{mm}}^{2} / {\text{year}}^{2} , $$
(6)

where \( \sigma_{{{\text{GNSS-}}\,{\text{noise}},i}}^{{}} \) are the standard uncertainties for the GNSS observations (noise) after rescaling by the variance factor 1.41 estimated in the computation of the empirical model (see Sect. 2.1). The correlation length becomes 150 km. The empirical and analytical covariance functions are plotted in Fig. 7. Using the above parameters, the gridded residual surface in Fig. 8 is obtained.

Fig. 7
figure 7

Empirical signal covariance function estimated using the 172 GNSS observations (red) and the analytical 1st-order Gauss–Markov function computed with variance 0.131 mm2/a2 and correlation length of 150 km (blue)

Full size image
Fig. 8
figure 8

The residual surface, i.e. the gridded difference between the empirical model and the GIA model. The observation points (model observations) are indicated by black dots

Full size image

The final NKG2016LU_abs model (Fig. 9) is obtained by restoring the GIA model to the residual surface [cf. Eq. (1)].

Fig. 9
figure 9

The final NKG2016LU_abs land uplift model. Unit: mm/a

Full size image

The residuals of NKG2016LU_abs are presented in Fig. 10 and their statistics in Table 2. They are very small and almost nothing remains of the systematic effects in Fig. 6. The semi-empirical model thus agrees very well with the empirical model near the observation points, just some minor extra filtering has been made. Farther away from the observations, NKG2016LU_abs approaches the GIA model (as the residual surface approaches zero in Fig. 8). The differences between the original clean GNSS observations and NKG2016LU_abs are presented in Fig. 11 and Table 2, second row. The latter values are considerably larger than the residuals, which is because of the filtering effect in LSC and the fact that the empirical model is computed from both GNSS and levelling observations and that the levelling sometimes has sufficient weight to overrule the GNSS observations. To what extent this effect will happen depends on the relative weighting of the GNSS and the levelling observations and the a priori variance–covariance function of the signals—in other words the key challenges of LSC.

Fig. 10
figure 10

Residuals of the NKG2016LU_abs model, i.e. differences between the model observations (empirical model) and NKG2016LU_abs. Unit: mm/a. Red means that the model observations are higher than the NKG2016LU_abs model

Full size image
Table 2 Statistics for the residuals of the NKG2016LU_abs model (differences between the empirical model in the observation points and NKG2016LU_abs) and the differences between the cleaned GNSS observations and NKG2016LU_abs
Full size table
Fig. 11
figure 11

Difference between the GNSS uplift rates (observations) and the NKG2016LU_abs model. Red indicates that the GNSS uplift rates are higher than the NKG2016LU_abs model

Full size image

The last step is to compute NKG2016LU_lev using Eq. (4).

Compared to NKG2005_abs, there are huge differences to NKG2016LU_abs outside the area of observations as shown in Fig. 12. These are mainly caused by the replacement of the GIA model from Lambeck et al. (1998) with the NKG2016GIA_prel0306 model. In the north-west, it also plays an important role that the values in NKG2005 were deliberately truncated to − 0.72 mm/a here. Inside the area of observations, the differences are much smaller with differences mainly less than 0.5 mm/a. In northern Denmark and the northernmost part of Norway however, the differences reach 1.0 mm/a and have other explanations. For Denmark, a reason may be that no levelling dataset is used in the NKG2005LU-calculation, and for Northern Norway we suspect that the tide gauge dataset used for NKG2005 could be part of the explanation since omitting/including tests show that the tide gauge data have the most effect here.

Fig. 12
figure 12

The difference between NKG2016LU_abs and NKG2005LU_abs. Positive values mean that the former is larger

Full size image

3 Uncertainty estimation

The main purpose of this section is to estimate a realistic standard uncertainty grid for the official semi-empirical model NKG2016LU_abs. The uncertainty of the purely empirical model is determined by standard uncertainty propagation in the method of least squares adjustment and least squares collocation. The challenge is to find good estimates for the combined model where we bring in a geophysical GIA model in the combination step. An important part is to estimate a standard uncertainty grid for the underlying GIA model and then a correlation function that makes it possible to compute variances and covariances for the errors of the GIA model (remember that in our LSC model, the error of the GIA model is viewed as the signal to estimate; cf. Sect. 2.3).

In Sect. 3.1, we calculate preliminary uncertainties for both the empirical and the semi-empirical model in a simple approach assuming a homogenous covariance function describing the GIA model uncertainty. In Sect. 3.2, we elaborate on the technique by developing an uncertainty grid for the GIA model, and in Sect. 3.3, we use this grid together with an empirical correlation function to calculate a final and official uncertainty grid for the already released semi-empirical model NKG2016LU. To check this uncertainty grid, we develop in Sect. 3.4 an alternative and more direct approach for calculating the semi-empirical uplift model as well as the corresponding uncertainties, and further demonstrate that the uncertainty grid obtained this way differs very little from the official one.

Below we often speak of just the uncertainty. We then always mean the one standard error uncertainty (1 σ).

In the following, we discuss the uncertainty of the absolute model NKG2016LU_abs. The uncertainty of NKG2016LU_lev, which gives the uplift relative to the rising geoid, is not specifically addressed. Since the latter is derived from the absolute model by subtracting the geoid rise given as a model (Fig. 5), the uncertainty will increase somewhat. However, as this geoid change rate is quite well known, this contribution to the uncertainty is for practical use negligible in relation to the total value. The magnitude for this rise is very low, less than 0.7 mm/a, and as such even a relative uncertainty as high as 10% contributes very little to the total uncertainty.

3.1 Preliminary uncertainties for the empirical and semi-empirical models

We first present the straightforward uncertainties obtained in standard least squares collocation/adjustment computations of the empirical and semi-empirical models, assuming in both cases the homogenous and isotropic signal covariance function used in the respective computations. The uncertainty grids are computed exactly as is described in Sects. 2.1 (and Section 1 in the Supplementary Material) and 2.3, respectively.

The uncertainty for the purely empirical model is illustrated in Fig. 13. In the area where we have observations, the uncertainty is about the same everywhere (Fig. 13a). As soon as we move away from that area, the uncertainty of the coefficients in the polynomial surface dominates and the estimates blow up (Fig. 13b).

Fig. 13
figure 13

Uncertainty of the empirical model. a For the area of observations only. b For the whole project area. The red dots are either observation points or so-called prediction points (points where the land uplift is predicted in the LSC process)

Full size image

In the combined solution, the accuracy of the model in this extended region will depend very much on the uncertainty of the GIA model. The estimated standard uncertainty grid for the semi-empirical model, using Eq. (3), is given in Fig. 14. We see that the uncertainties gradually approach 0.36 mm/a when one moves away from the observations. This maximum value is equal to the square root of the adopted signal variance, \( C_{0} \), which is a homogeneous measure of the uncertainty of the GIA model. This value is obtained empirically using Eq. (6) and assumed valid also outside the observation area. This is most likely too optimistic as the GNSS observations in question have been used in the selection of the best-fitting GIA model (See Section 2.4 in Supplementary Material).

Fig. 14
figure 14

Preliminary uncertainty of NKG2016LU_abs. The maximum uncertainty far away from the observations is equal to signal standard deviation = 0.36 mm/a. The contour interval is 0.05 mm/a

Full size image

Another problem for the semi-empirical uncertainty grid is that it becomes too low in areas with dense observations. By comparing Figs. 13a and 14, we can see that the uncertainties are almost 2 times lower for the semi-empirical model in the areas with densest observations (like in southern Norway). This is because the correlation between the observations is not considered in the LSC solution in question; see Sect. 2.3 and matrix D in Eq. (3). This does not affect the estimated residual surface much at all, but the uncertainty estimation is more sensitive in this respect.

We conclude that both the collocation uncertainty grids in Figs. 13 and 14 are problematic in that the GIA uncertainty is not well modelled outside the area of observations. The semi-empirical model is best in this respect, but it is still not realistic that the GIA model accuracy is the same everywhere; the uncertainty is more likely to vary and can even be higher than 0.36 mm/a. We should first know the uncertainty of the GIA model, which then should be combined with the uncertainty of the empirical model to obtain the uncertainty of the final semi-empirical model.

3.2 Uncertainty of the GIA model

Here, we derive an uncertainty grid for the NKG2016GIA_prel0306 model (Sect. 2.2) and then estimate a homogeneous and isotropic correlation function, which will be used to compute heterogeneous and non-isotropic variances and covariances in Sects. 3.3 and 3.4. In other words, we will introduce individual variances for all node points, but keep the same isotropic correlation between them.

This best-fitting GIA model NKG2016GIA_prel0306 was selected in a fitting procedure that involved 11,025 different GIA models as described in the Supplementary Material. We compute the standard deviation for a subset of good GIA models (including the best-fitting NKG2016GIA_prel0306) at each grid node. This provides a grid, \( \varGamma_{\text{STD}} \), with higher standard deviation in areas where the GIA model is more sensitive to the choice of Earth and ice models, which is exactly what we want. However, GIA models are also affected by other types of errors due to various assumptions and approximations (e.g. assumed viscoelastic rheology, 1D distribution of rheological parameters, non-modelled effects, presence of non-GIA effects like tectonic motion in the fitted observations, etc.). It is difficult to model the latter type of uncertainty. In this paper, we assume that the total GIA modelling uncertainty is the square root of the squared sum of \( \varGamma_{\text{STD}} \) and a homogeneous uncertainty factor \( E \) that models all “other types of errors”,

$$ \varGamma_{\text{GIA}} = \sqrt {\varGamma_{\text{STD}}^{2} + E^{2} } . $$
(7)

The extra term \( {\rm E} \) is estimated using the differences between the GNSS observations and the GIA model.

To select the subset for \( \varGamma_{\text{STD}} \)-calculation, we compare 11,025 models with the selected GIA model (NKG2016GIA_prel0306) by weighted sum of squared residuals (ψ):

$$ \psi = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^{n} {\left( {\frac{{p_{\text{best}} - p_{i} (a_{j} )}}{{\Delta o_{i} }}} \right)^{2} } } . $$
(8)

Here, Pbest is NKG2016GIA_prel0306, pi(aj) is the predicted value for model aj in point i, Δoi is the uncertainty of a test observation in point i and n the number of points. These points are either the locations of GNSS uplift rates or RSL observations. The points are here restricted to central Fennoscandia where the uplift is higher than 6 mm/a to avoid very unreliable values from some models that rather best fit the periphery of the GIA-affected area. Such unreliable values may ruin the uncertainty estimation. This is because there are naturally models which have very small differences to NKG2016GIA_prel0306 in areas with low uplift or subsidence but large differences in the uplift centre. Nonetheless, the weighted sum can result in a very low ψ-value comparable to that of models that, in turn, provide a very good fit in the uplift centre but not in the outskirts.

By using all models with ψ-value lower than 1.5, we end up with 21 models. This limit, 1.5, is chosen after testing different values from 1 to 2 in 0.125 steps and analysing the number of different ice and Earth models and their statistics (Table 3). The 21 models falling within the 1.5 limit cover 7 different ice models (out of 25) and 14 different Earth models (out of 441) that result in an average difference to NKG2016GIA_prel0306 of 0.66 mm/a and a difference of 0.53 mm/a at the location of the maximum uplift. This is a good balance of model numbers and, in view of the GNSS uncertainty, maximum and average differences to NKG2016GIA_prel0306. A ψ value of 1.75 or higher, for example, yields quite high differences above 0.9 mm/a at the maximum uplift location, which we judge is quite well determined by GNSS within 0.4 mm/a. Hence, the difference of 0.53 mm/a as found for a ψ value of 1.5 seems reliable. Choosing a ψ value of 1.625 instead of 1.5 does not change the difference statistics much, although a larger number of GIA models fall within the 1.625 limit.

Table 3 Tested ψ value limits, resulting number of different GIA models as well as different ice and Earth models that fall within this limit, the maximum difference to NKG2016GIA_prel0306 (mm/a) at the maximum land uplift location and over the whole area and the average difference (mm/a) to NKG2016GIA_prel0306
Full size table

We then compare all these 21 models plus NKG2016GIA_prel0306 and compute the standard deviation (STD) for each grid node, which gives \( \varGamma_{\text{STD}} \) in Eq. (6).

The next step is to compute uncertainty factor E by studying the agreement between the 172 cleaned GNSS observations and the GIA model NKG2016GIA_prel0306.

To compute the uncertainty term E, we consider the variance \( \sigma_{0}^{2} \) of the standardized difference between GNSS and GIA model (i.e. difference between GNSS and GIA, divided by the uncertainty of this difference),

$$ \sigma_{0}^{{}} = \sqrt {\sum\limits_{i = 1}^{172} {\frac{{\left( {\dot{h}_{{{\text{GNSS}},i}}^{{}} - \dot{h}_{\text{GIA}}^{{}} } \right)^{2} }}{{\left( {\sigma_{{{\text{GNSS}},i}}^{2} + \varGamma_{{{\text{GIA}},i}}^{2} } \right)}}} } = \sqrt {\sum\limits_{i = 1}^{172} {\frac{{\left( {\dot{h}_{{{\text{GNSS}},i}}^{{}} - \dot{h}_{\text{GIA}}^{{}} } \right)^{2} }}{{\left( {\sigma_{{{\text{GNSS}},i}}^{2} + \varGamma_{{{\text{STD}},i}}^{2} + E^{2} } \right)}}} } , $$
(9)

where \( \sigma_{{{\text{GNSS}},i}}^{{}} \) is the uncertainty of GNSS observation i multiplied by the variance factor 1.41 estimated in Sect. 2.1. For E equal to zero, we get \( \hat{\sigma }_{0} \) = 1.45. We choose E so that \( \sigma_{0} \) becomes exactly equal to 1, which happens when E = 0.34 mm/a. This yields the final GIA uncertainty grid in Fig. 15.

Fig. 15
figure 15

Uncertainty grid for the GIA model NKG2016GIA_prel0306 derived using Eq. (7) with E = 0.34 mm/a. The contour interval is 0.05 mm/a

Full size image

3.3 Improved uncertainty estimates for the semi-empirical model

In the following, we derive a better uncertainty grid for NKG2016LU_abs than Fig. 14 by making use of the above GIA uncertainty grid in Fig. 15 and derive the corresponding covariance function. For this purpose, we assume the following type of heterogeneous and non-isotropic covariance function COV for the errors in the GIA model,

$$ {\text{COV}}\left( {\psi_{PQ} } \right) = \sigma_{{{\text{GIA}},P}} \sigma_{{{\text{GIA}},Q}} {\text{CORR}}\left( {\psi_{PQ} } \right), $$
(10)

where \( \psi_{PQ} \) is the spherical distance between points P and Q, \( {\text{CORR}} \) is the correlation function, and \( \sigma_{{{\text{GIA}},i}} \) is the standard deviation for point i obtained from the GIA uncertainty grid, either directly or by interpolation.

As it is assumed that the correlation function is homogeneous and isotropic, we can estimate an empirical correlation function using the differences of the 172 clean GNSS observations from the GIA model in the same way as described in Sect. 2.3 and Eq. (5). The only difference is that we now divide by the respective standard uncertainties \( \sigma_{{{\text{GIA}},P}} \) and \( \sigma_{{{\text{GIA}},Q}} \), taken from the GIA uncertainty grid, to obtain correlation instead of variance,

$$ {\text{CORR}}\left( \psi \right) = \sum\limits_{(P,Q)} {\frac{{\left( {\dot{h}_{{{\text{GNSS}},P}} - \dot{h}_{{{\text{GIA}},P}} } \right)\left( {\dot{h}_{{{\text{GNSS}},Q}} - \dot{h}_{{{\text{GIA}},Q}} } \right)}}{{\sigma_{{{\text{GIA}},P}} \sigma_{{{\text{GIA}},Q}} }}} , $$
(11)

where the summation extends over all pairs of points P and Q in the interval \( \psi - \Delta < \psi_{PQ} < \psi + \Delta \) and where the mean value is not subtracted; cf. the comment after Eq. (5). We assume that the GNSS errors are uncorrelated, which means that they do not affect Eq. (11) if \( \psi > 0 \). For the zero distance, we need to divide by the variance \( \left( {\sigma_{{{\text{GNSS}},P}}^{2} + \sigma_{{{\text{GIA}},P}}^{2} } \right) \), where the GNSS variance has been multiplied by 1.412 as before. Note that our choice of E by normalizing Eq. (9) implies that the empirical estimate of \( {\text{CORR}}\left( 0 \right) \) becomes exactly 1. The empirical correlation function for the present case is illustrated in Fig. 16 together with a fitted analytical correlation function (1st order Gauss–Markov) with variance 1 and correlation length 100 km.

Fig. 16
figure 16

Empirical correlation function for the GIA model errors estimated using the GNSS observations and the GIA uncertainty grid in Fig. 15 (red line). The analytical 1st-order Gauss–Markov function computed with unit variance and correlation length

Full size image

Using the estimated analytical correlation function in Fig. 16, Eq. (10) and the GIA model uncertainty grid of Fig. 15, we are now able to compute the covariance function for the GIA model errors and obtain a new and updated variance–covariance matrix \( C_{\text{GIA}} \). The next step is to use this new covariance information to compute an updated uncertainty grid for the semi-empirical NKG2016LU_abs. To do this, we rewrite Eqs. (1) and (2) as

$$ \dot{h}_{{{\text{NKG2016LU\_abs,}}ij}}^{{}} = \dot{h}_{{{\text{GIA}} .ij}}^{{}} + A_{ij}^{T} \left( {\dot{h}_{\text{empirical}}^{{{\text{obs}} . {\text{pts}}}} - \dot{h}_{\text{GIA}}^{{{\text{obs}} . {\text{pts}}}} } \right) = \dot{h}_{{{\text{GIA}} .ij}}^{{}} + A_{ij}^{T} \dot{h}_{\text{empirical}}^{{{\text{obs}} . {\text{pts}}}} - A_{ij}^{T} \dot{h}_{\text{GIA}}^{{{\text{obs}} . {\text{pts}}}} , $$
(12)

or in a rearranged matrix form,

$$ \dot{h}_{{{\text{NKG2016LU\_abs,}}ij}}^{{}} = \left[ {\begin{array}{*{20}c} 1 & { - A_{ij}^{T} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\dot{h}_{{{\text{GIA}} .ij}}^{{}} } & {\dot{h}_{\text{GIA}}^{{{\text{obs}} . {\text{pts}}}} } \\ \end{array} } \right]^{T} + A_{ij}^{T} \dot{h}_{\text{empirical}}^{{{\text{obs}} . {\text{pts}}}} , $$
(13)

where the coefficients in the column vector \( A_{ij}^{{}} \) are exactly those used to compute NKG2016LU_abs in Sect. 2.3. With updated variances and covariances for the GIA model errors, we recalculate the uncertainty. If we assume that the errors of the empirical model observations and GIA model are uncorrelated, and split the variance and covariance of the GIA model into \( \left[ {\begin{array}{*{20}c} {\sigma_{{{\text{GIA}} .ij}}^{2} } & {C_{{{\text{GIA}},ij}}^{T} } \\ {C_{{{\text{GIA}},ij}} } & {C_{\text{GIA}} } \\ \end{array} } \right] \), then straightforward uncertainty propagation yields,

$$ \sigma_{{{\text{NKG2016LU\_abs,}}ij}}^{2} = \sigma_{{{\text{GIA}} .ij}}^{2} + A_{ij}^{T} \left( {C_{\text{GIA}} + D_{\text{empirical}} } \right)A_{ij}^{{}} - 2A_{ij}^{T} C_{{{\text{GIA}},ij}} , $$
(14)

where \( \sigma_{{{\text{GIA}}.ij}}^{2} \) is the variance of the GIA model error in the node point ij, \( C_{\text{GIA}} \) is the variance–covariance matrix of the GIA model errors in the observation points, \( C_{{{\text{GIA}},ij}} \) is the corresponding cross-covariance column vector between grid node (i, j) and the observation points, and \( D_{\text{empirical}} \) is the variance–covariance matrix for the observation noise (see below). Note that before (e.g. in Sect. 2.3) we referred to the GIA model error as signals, now we consider all errors as noise and use the derived variances and covariances and apply the statistical law of uncertainty propagation directly on Eq. (13), keeping \( A_{ij}^{{}} \) fixed to the existing values from Sect. 2.3. Note further that, if the same homogeneous and isotropic variance–covariance matrices as in Sect. 2.3 are used in Eq. (14), then this equation turns into Eq. (3) exactly.

Another problem with the preliminary standard uncertainty grid for the semi-empirical model NKG2016LU_abs in Fig. 14 is that the correlations between the (purely) empirical model observations are not properly considered. We now model these correlations by assuming an analytical 1st-order Gauss–Markov correlation function and set the correlation length to the length used for the computation of the purely empirical model, i.e. 40 km; see Figure S3 in the Supplementary Material. This correlation function is then used to compute the matrix \( D_{\text{empirical}} \) from the individual observation standard uncertainties using Eq. (10). Note that this correlation length is based on residuals with respect to the signals. In reality, the uplift in the model is estimated signals plus trend, and a longer correlation length is very likely. Using the full covariance matrix from the empirical calculation would be the optimal solution. Using the above correlation function is a stepwise improvement, and we refer to Sect. 3.4 for a more optimal handling of the observation noise.

The final NKG2016LU_abs uncertainty grid computed according to the above developments is illustrated in Fig. 17.

Fig. 17
figure 17

Final uncertainty grid of NKG2016LU_abs, based on the GIA model uncertainty grid in Fig. 15, the analytical correlation function in Fig. 16 and the above modelling of correlations for the empirical model observations. The two plots only differ in that they have different colour scales. Plot a has the same colour scale as Fig. 14. Both figures have a contour interval of 0.05 mm/a

Full size image

Now, we have replaced the homogeneous and isotropic covariance function used in Sect. 2.3 (in Eqs. 2 and 3) with a more realistic function, which is constructed from a homogeneous and isotropic correlation function. We have also changed the variance–covariance matrix D so that the correlations of the empirical model observations are taken care of (no longer diagonal). These refinements to the error model imply that NKG2016LU_abs is now non-optimal. The reason for not making a new LSC solution based on the above updates is that the NKG2016LU_abs model has already been officially released (in June 2016). The next subsection examines how the uncertainties of such a new optimal LSC collocation solutions differ from the uncertainty estimates in Fig. 17.

3.4 Heterogeneous least square collocation in one step

As mentioned in Sects. 1 and 2, calculating the combined semi-empirical uplift model in one step is an alternative to the two-step approach used for the officially released NKG2016LU model. A one-step approach is more straightforward and makes the uncertainty calculation easier and more optimal in a statistical sense. This is a promising future development, which is presented here mainly to check the quality of the above uncertainty grid for NKG2016LU illustrated in Fig. 17.

In the one-step approach, we replace the 5th degree polynomial surface used as trend in the purely empirical model of Sect. 2.1 with the GIA model, and calculate the difference between GIA model and observation as signals. Just as in Sect. 2.3, we assume no significant offset between ITRF2008 and GIA model. We use the above-mentioned variance and covariance as statistical characteristics of these signals, meaning that we use the GIA model uncertainty grid in Fig. 15 and the analytical correlation function in Fig. 16. I.e. the GIA uncertainty grid is \( \varGamma_{\text{GIA}} = \sqrt {\varGamma_{\text{STD}}^{2} + \left( {0.34\,{\text{mm/a}}} \right)^{2} } \) and the correlation length is 100 km.

The fundamental difference concerning uncertainty estimation is that we now use the given observation uncertainties directly when forming the D matrix, later used to calculate the normal equations.

Since we no longer estimate a trend surface, we must modify the observation equation given in Subsection 1.2 in Supplementary Material. We have for the levelling observations,

$$ l_{j} = \, - 1{\mathbf{h}}_{{\mathbf{a}}} + \, 1{\mathbf{h}}_{{\mathbf{b}}} + \, t_{j} \left\{ {{\text{GIA}}_{b - a} } \right\} \, {-}t_{j} {\mathbf{s}}_{{\mathbf{a}}} + \, t_{j} {\mathbf{s}}_{{\mathbf{b}}} + \, n_{j} , $$
(15)

and for the CORS data,

$$ l_{j} = {\text{GIA}}_{i} + \, 1{\mathbf{s}}_{{\mathbf{i}}} + \, n_{j} , $$
(16)

where GIAb−a is the difference in GIA model value between point a and point b, and GIAi is the GIA model value in CORS i. Furthermore, ha and hb are the unknown height in point a and b, respectively, and sa and sb are the unknown signals in the same points. For the CORS, si is the unknown signal in CORS i. For both observations types, lj is the observed value and nj the error in observation j. Finally, tj is the year of observation of levelling line j. Since the rise of the geoid affects the two types of observations differently, we correct the levelling observations from levelled to absolute uplift using the model in Fig. 5 before the adjustment.

In matrix form, we have

$$ {\mathbf{l}} \, = \, {\mathbf{Ax}} \, + \, {\mathbf{Bs}} \, + \, {\mathbf{n}} $$
(17)

where l is an observation vector, A is the design matrix for the unknowns, B is the design matrix for the signals, n is the noise vector, x is the unknown vector, i.e. unknown heights at all levelling benchmarks, and s is the unknown signal vector.

All GIA model terms in Eqs. (15) and (16) can be moved to the left-hand side of the equality sign and become a part of the observation vector l. As parametric part of the system of equations, only the unknown heights remain.

As explained in Vestøl (2006) and mentioned in the Supplementary Material, we here do not calculate the signals in the normal way described in, e.g. Moritz (1980), that is used in Sect. 2.3. Instead, we follow the equivalent method suggested by Schwarz (1976) and consider also the stochastic signals as unknown parameters. This is allowed if we add the inverted variance–covariance matrix for the signals, i.e. the C matrix, to the lower right corner of the normal matrix N (cf. Figure 4 in Vestøl 2006). By doing so, the variance–covariance matrix of the calculated semi-empirical land uplift model in the observation points is just the variance–covariance of the estimated signals, easily found as m 2o N−1, where N−1 is the inverted normal matrix and m 2o is the variance of the unit weight.

We also want to find the variance and covariance for the predicted signals, i.e. in points with no observation (e.g. grid points). This is done, similar to ordinary least squares adjustment, by expanding B in Eq. (17) with a null matrix and, in addition, expanding the C matrix to also include variances for predicted signals and covariance between signals in points with and without observations,

$$ B_{\rm new} = \left[\begin{array}{*{20}c} B & {\mathop {\rm null}\limits^{{k} \times {m}} } \\ \end{array} \right]\;{\rm and}\;{C}_{\rm new} = \left[ {\begin{array}{*{20}c} {\mathop {{C}_{\rm obs} }\limits^{n \times n} } & {\mathop {C_{{\rm obs}/{\rm pred}}^{T} }\limits^{n \times m} } \\ {\mathop {C_{{\rm obs}/{\rm pred}} }\limits^{m \times n} } & {\mathop {C_{{\rm pred}} }\limits^{m \times m} } \\ \end{array} } \right], $$
(18)

where n is the number of points with observation(s), m is the number of prediction points, and k is the number of observations. The D matrix, i.e. the variance–covariance matrix for the observations, is in use when calculating the normal matrix N and forming the right side of the normal equations.

The uncertainty model obtained by following this procedure is shown in Fig. 18. When we compare this with Fig. 17, we see very similar values. It is also interesting to compare with the uncertainty of the purely empirical model in Fig. 13. We note that we now have realistic uncertainty estimates also outside the observation area.

Fig. 18
figure 18

a Alternative NKG2016LU uncertainty grid using the one-step approach. b Difference to the official uncertainty grid in Fig. 17. The uncertainty grid in a is not part of the official NKG2016LU release

Full size image

Using this one-step approach results in an uplift model that differs slightly from the officially released one; see Fig. 19 in Sect. 4.

Fig. 19
figure 19

Difference between NKG2016_abs and a model created using the one-step approach

Full size image

4 Discussion and conclusions

We have shown how we have calculated a purely empirical land uplift model from geodetic observations and a GIA model based on a geophysical approach and further how we have merged these two models to obtain a more accurate and reliable semi-empirical land uplift model, NKG2016LU. We have then developed the corresponding uncertainty grids, where we have tried to take all error sources into account. The two-step procedure has been used to compute the final NKG2016LU model and its uncertainty grid. Since the official release of NKG2016LU, we have developed a new one-step method, which is more correct from the statistical point of view. Here, the one-step method has been used mainly to check and verify the NKG2016LU uncertainty grid.

4.1 The role of the tide gauges

In the previous model from 2005, one important dataset was the tide gauge rates of Ekman (1996), consisting of 58 stations in the Baltic Sea and adjacent waters, including coastal Norway. We have omitted this dataset and rejected the whole idea of using tide gauge rates to estimate the land uplift.

There are several reasons behind this decision. One is that there are spatial variations in the mean sea-level rise that makes it difficult to separate this effect from the land uplift. There is also a temporal variation in the mean sea-level rise, making the separation even more difficult. The apparent uplift computed at the tide gauges will always refer to a certain time interval, and it is difficult to combine tide gauge rates referring to different periods.

By omitting the tide gauge rates, the postglacial land uplift model becomes independent of tide gauge and sea-level-related information. This is of principal importance when correcting for vertical land motion in climate studies when mean sea-level variation is involved. We can then use an independent model to correct relative sea-level rise for vertical land motion. However, an omitting/including test shows that the effect of omitting the tide gauge observations is generally less than 0.1 mm/a (after conversion to the same type of uplift). One exception is the northernmost part of Norway where the effect reaches 0.4 mm/a. This could be due to the geographical variation of sea-level rise, but further studies are necessary to find the reason for the difference in that region.

4.2 Reference frame

It is important to be aware of the relationship between a physical phenomenon and the reference frame in which we express this phenomenon. The NKG2016LU land uplift model refers to ITRF2008. This is not necessarily the right frame to describe a geophysical process. ITRF2008 is a global reference frame with its theoretical origin at the Earth’s centre of mass (including oceans and atmosphere) (Altamimi et al. 2011), while GIA is mainly related to solid earth processes and might be better described in a centre of Earth reference frame where the atmosphere is omitted.

Using another reference frame will give slightly different values. Uncertainties, differences and drift in the reference frame might bias the results (see e.g. Argus 2012, for a thorough discussion of the uncertainty of ITRF2008). A possible drift of − 0.6 ± 0.2 mm/a along the z-axis in ITRF2008, as estimated by Wu et al. (2011), will have a systematic effect on the GNSS rates and the final uplift model. Kierulf et al. (2014) showed that the scale drift, causing the difference between ITRF2000 and ITRF2008, resulted in different best fit GIA models.

To avoid the problem with possible drift in the reference frame, and to ensure consistency between geophysical processes (e.g. GIA) and the frame, Kierulf et al. (2014) transformed the GPS velocity field to a reference frame determined by the GIA model subject to validation. The method was named the GIA-frame approach. However, for GIS and positioning applications, a model consistent with a conventional reference frame, like ITRF2008, is preferable.

4.3 One- or two-step approach

The one-step procedure is a recent development of our method, which will be used to compute future semi-empirical models. The reason for including it already here for uncertainty is to check and verify the final standard uncertainty grid of NKG2016LU, which was based on the two-step approach described in Sect. 3.3 and illustrated in Fig. 17.

The two alternative uncertainty grids for the NKG2016LU_abs uplift model differ very little (see Fig. 18b). For the observation area, the uncertainty varies from 0.1 to 0.3 mm/a in both alternatives. Outside this area, the uncertainty of the GIA model influences both alternatives in almost the same way with extremes in the Norwegian Sea outside the north-western coast of Norway, and at the northern coast of the Kola Peninsula in Russia.

Both alternatives have up and down sides. One is more optimal in the statistical sense, but is not strict in accordance with the way the released NKG2016LU model was actually calculated. The other represents the calculation procedure, but takes a few extra shortcuts concerning variance and covariance for signal and noise, mainly due to the two-step division of the problem.

To demonstrate the effect on the uplift model itself, we have tentatively calculated an alternative grid model using the one-step approach, see Fig. 19. The effect is, with a few exceptions, within the calculated standard error. Comparing this with Fig. 17, showing the uncertainty, we see that the differences are mostly within the level of the uncertainty. However, in some areas there seem to be systematic differences larger than the standard error. We have most confidence in the alternative model since the calculation is more optimal in a statistical sense.

4.4 Future developments

We have fitted our GIA model only to RSL data and the uplift component of the GNSS data. Future steps will therefore include additional fits to satellite and terrestrial gravity data as well as GNSS-derived horizontal velocities. The latter may require the inclusion of an extra sub-lithospheric layer in the Earth model as suggested by Peltier and Drummond (2008) to get a good fit. This will at the same time increase the number of model combinations and may also affect the weighing ratios in the fitting test. This also applies when gravity observations are used. Further constraints that could reduce the number of reliable models may be obtained from a comparison of stress changes at the location of prominent glacially induced faults evident in northern and central Europe (e.g. Lagerbäck and Sundh 2008; Brandes et al. 2015; Lund 2015).

Concerning the calculation of semi-empirical models, the described two-step combination approach will be replaced by the more direct one-step counterpart, where we just improve the GIA model using LSC applied to geodetic observations. Nevertheless, we still think that the calculation of purely empirical models that can be compared with GIA models will continue to be very valuable. For example, when we compare the empirical and GIA models in the present study, we notice remarkable differences in some areas such as the Norwegian Atlantic coast. In Kierulf (2017), a denser network of both campaign-based observations and permanent GNSS stations at the west coast of Northern Norway was analysed in detail. He found that the GIA models agree with the observations at 0.5 mm/a level for most stations except for the outermost coastal ones. Here, the GIA model was not able to reproduce the steep gradient and gave more than 1 mm/a higher uplift. Kierulf (2017) points at neo-tectonic processes and erroneous GIA models, e.g. the omission of lateral heterogeneities in the Earth model, as possible explanations. There is much to learn by comparing strictly empirical and GIA models.

As we have not tried to separate the long wavelength uplift caused by postglacial rebound from uplift having other geophysical explanations, there is an interesting and likely difficult issue left for further studies to identify the sources of this difference.