Optimal Strategy of a GPS Position Time Series Analysis for Post-Glacial Rebound Investigation in Europe

We describe a comprehensive analysis of the 469 European Global Positioning System (GPS) vertical position time series. The assumptions we present should be employed to perform the post-glacial rebound (PGR)-oriented comparison. We prove that the proper treatment of either deterministic or stochastic components of the time series is indispensable to obtain reliable vertical velocities along with their uncertainties. The statistical significance of the vertical velocities is examined; due to their small vertical rates, 172 velocities from central and western Europe are found to fall below their uncertainties and excluded from analyses. The GPS vertical velocities reach the maximum values for Scandinavia with the maximal uplift equal to 11.0 mm/yr. Moreover, a comparison between the GPS-derived rates and the present-day motion predicted by the newest Glacial Isostatic Adjustment (GIA) ICE-6G_C (VM5a) model is provided. We prove that these rates agree at a 0.5 mm/yr level on average; the Sweden area with the most significant uplift observed agrees within 0.2 mm/yr. The largest discrepancies between GIA-predicted uplift and the GPS vertical rates are found for Svalbard; the difference is equal to 6.7 mm/yr and arises mainly from the present-day ice melting. The GPS-derived vertical rates estimated for the southern coast of the Baltic Sea are systematically underestimated by the GIA prediction by up to 2 mm/yr. The northern British Isles vertical rates are overestimated by the GIA model by about 0.5 mm/yr. The area of the Netherlands and the coastal area of Belgium are both subsiding faster than it is predicted by the GIA model of around 1 mm/yr. The inland part of Belgium, Luxemburg and the western part of Germany show strong positive velocities when compared to the GIA model. Most of these stations uplift of more than 1 mm/yr. It may be caused by present-day elastic deformation due to terrestrial hydrology, especially for Rhein basin, or non-tidal atmospheric loading, for Belgium and Luxembourg.


Introduction
The global positioning system (GPS) position time series have been employed for almost three decades to analyze and interpret numerous dynamic phenomena that cause deformations of the Earth's surface. Plate motion [1], earthquake-related processes [2], hydrospheric effects [3,4], and the vertical land motion (VLM) employed to correct tide-gauge records and to study local subsidence [5][6][7] are, among others, the principal applications of these position time series. To analyze any of the geodynamic phenomena, the velocity of the permanent station, i.e., a long-term trend estimated directly from the position time series, is of a great interest. However, its uncertainty is indispensable to decide on the significance of the phenomena sensed by GPS.  We use a time series of a minimum 5-year time span. The time series ends in mid-2018, meaning that the longest observation span is 23 years, i.e., the years from 1995 to 2018. We pre-process the vertical position time series, removing the outliers using three times the interquartile range rule (IQR, i.e., 3IQR). The epochs of offsets are employed from the NGL database (http://geodesy.unr.edu/ NGLStationPages/steps.txt), and supported by manual inspection of the series. In total, 1240 epochs of offsets are used. Since, in the consecutive part of the research, we use the modified maximum likelihood estimation (MLE) [48], no interpolation is performed to fill in the missing observations. The average number of gaps within the series is equal to 10%.
To compare the GPS-derived uplift with the expected present-day motion arising from the PGR, we employ the newest version of the GIA predictions, i.e., ICE-6G_C (VM5a) [43], referred to as ICE-6G herein ( Figure 2). We make use of 1 • per 1 • global grid estimates at the time 0.0 Ka (kiloannus, 10 3 years). This model was refined in comparison to the previous versions, i.e., ICE-5G (VM2) [49], by adding the GPS measurements of the vertical motion as an additional input, and by specifying the exact reference frame of these observations. Also, additional constraints have been applied to optimally fit, but only the Antarctic component to the GPS observations [50]. The present-day uplift predicted from the ICE-6G GIA model is also employed to construct the hinge line.

Methods
We employ the vertical (Up) position changes and construct the time series model using both deterministic and stochastic parts, similarly to [51]: where U0 is the initial value, the long-term trend in the form of a linear or non-linear signal is modelled as a j-degree polynomial within the j j i a t ⋅ sum, the seasonal signatures from k=1 to O are characterized by their amplitudes A and phase shifts ϕ and assumed as in [10]. The number of M=1240 offsets are represented by their epochs ti and magnitudes p. All of the above constitutes the deterministic part of the time series. Finally, ε stands for the residuals which remained after the mathematical model was fitted; the so-called stochastic part. The entire analysis is performed using the maximum likelihood estimation algorithm (MLE, [52]) in the Hector software [48]. Within the analysis, we examine polynomial degrees from 1 to 6. Our choice of the preferred degree is based on the bayesian information criterion (BIC, [53]), Figure 3 presents an exemplary vertical time series from station Lochgilphead located in Scotland.
As was previously shown, the assumption on white-noise-only leads to an underestimation of velocity uncertainty of several times [32,54], as the covariance matrix of observations is created based on the amplitude of white-noise-only a and the identity matrix I, as in: Adding a power-law noise component, the covariance matrix is being rebuilt to include both noise models inside, as in:

Methods
We employ the vertical (Up) position changes and construct the time series model using both deterministic and stochastic parts, similarly to [51]: where U 0 is the initial value, the long-term trend in the form of a linear or non-linear signal is modelled as a j-degree polynomial within the a j · t j i sum, the seasonal signatures from k = 1 to O are characterized by their amplitudes A and phase shifts φ and assumed as in [10]. The number of M = 1240 offsets are represented by their epochs t i and magnitudes p. All of the above constitutes the deterministic part of the time series. Finally, ε stands for the residuals which remained after the mathematical model was fitted; the so-called stochastic part. The entire analysis is performed using the maximum likelihood estimation algorithm (MLE, [52]) in the Hector software [48]. Within the analysis, we examine polynomial degrees from 1 to 6. Our choice of the preferred degree is based on the bayesian information criterion (BIC, [53]), Figure 3 presents an exemplary vertical time series from station Lochgilphead located in Scotland.
As was previously shown, the assumption on white-noise-only leads to an underestimation of velocity uncertainty of several times [32,54], as the covariance matrix of observations is created based on the amplitude of white-noise-only a and the identity matrix I, as in: Adding a power-law noise component, the covariance matrix is being rebuilt to include both noise models inside, as in: where b κ is the amplitude of the power-law noise characterized by the spectral index κ, and J κ is the power-law noise matrix, which includes the dependencies between individual observations. compare this approach with the previously used white-noise-only assumption, showing how much of the velocity uncertainty can be underestimated with improper assumptions. In our research, we do not account for the elastic uplift arising from the melting of present-day glaciers and ice sheets, as well as continental hydrology [56]. Instead, we analyze the differences between the predicted contemporary GIA uplift and the GPS-derived vertical velocities and present the areas, for which the GIA model may under-or over-estimate the evident large-scale crustal movements. Figure 3. The GPS vertical position time series (in mm) for the LOCG station (Lochgilphead, UK) along with the time series model plotted, respectively, in red and black. The 4 th polynomial degree is employed to model the non-linearity within the series. The LOCG station is estimated to subside by 0.05±1.23 mm/yr. The large uncertainty of this vertical rate is caused by the significant non-linearity that the time series is characterized by. Due to its uncertainty, the vertical rate has been found as insignificant and removed from the GPS-GIA comparison, see text for details.

Results
In the following section, we present the results of the deterministic part modelling, as well as the noise analysis we perform with the MLE algorithm. Then, we focus on the statistical significance of the vertical velocities, showing which of the GPS-derived vertical rates can be used to model the along with the time series model plotted, respectively, in red and black. The 4th polynomial degree is employed to model the non-linearity within the series. The LOCG station is estimated to subside by 0.05±1.23 mm/yr. The large uncertainty of this vertical rate is caused by the significant non-linearity that the time series is characterized by. Due to its uncertainty, the vertical rate has been found as insignificant and removed from the GPS-GIA comparison, see text for details.
Beyond the spectral index, the amplitude of noise is another factor which contributes to the velocity uncertainty estimates. The higher the amplitude of noise is, the greater impact it has on the uncertainty. When combined together, the velocity variance may be estimated from [55]: where N is the number of observations, Γ is the Gamma function and ∆T is the sampling rate.
In this research, we examine the position time series using the MLE algorithm. We find the preferred time series model using the BIC criterion assuming a combination of white and power-law noises, which is optimal for the GPS position time series, Equation (1). In this way, we derive the velocity uncertainties to assess their statistical significance. We assume that the velocity is significant and can be interpreted in terms of the PGR phenomena, only if the velocity value exceeds its uncertainty estimated using the model that is preferred for the particular time series. Finally, we compare this approach with the previously used white-noise-only assumption, showing how much of the velocity uncertainty can be underestimated with improper assumptions.
In our research, we do not account for the elastic uplift arising from the melting of present-day glaciers and ice sheets, as well as continental hydrology [56]. Instead, we analyze the differences between the predicted contemporary GIA uplift and the GPS-derived vertical velocities and present the areas, for which the GIA model may under-or over-estimate the evident large-scale crustal movements.

Results
In the following section, we present the results of the deterministic part modelling, as well as the noise analysis we perform with the MLE algorithm. Then, we focus on the statistical significance of the vertical velocities, showing which of the GPS-derived vertical rates can be used to model the crustal movements. Finally, we compare the present-day uplift predicted by the GIA model with the one derived from the GPS observations, indicating the differences between both rates.

GPS-Derived Vertical Velocities
The vertical velocities computed using the MLE algorithm for stations situated in Europe range between −2.3 mm/yr for the LIL2 (Lille, France) station and 11.0 mm/yr for the SKE0 (Frostkåge, Sweden) station; see Figure 4. The largest uplift, as expected, is found in Scandinavia, affected the most by the PGR [36]. The vertical rates derived for Sweden are greater than those derived for the remaining areas. All 27 Swedish stations are uplifting more than 1 mm/yr, with 18 stations exceeding the value of 5 mm/yr. Moving closer to the south coast of the Baltic Sea-the Danish, German and Polish borders-the vertical rates change either values or their signs from strongly positive to slightly negative. The area of the United Kingdom is divided into two halves; the northern part is characterized by uplift, while the southern part is subsiding. Many local movements are noticed for the areas of France, Belgium and the Netherlands with vertical rates of −2 to 2 mm/yr, which we describe further in this paper.
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 17 crustal movements. Finally, we compare the present-day uplift predicted by the GIA model with the one derived from the GPS observations, indicating the differences between both rates.

GPS-Derived Vertical Velocities
The vertical velocities computed using the MLE algorithm for stations situated in Europe range between -2.3 mm/yr for the LIL2 (Lille, France) station and 11.0 mm/yr for the SKE0 (Frostkåge, Sweden) station; see Figure 4. The largest uplift, as expected, is found in Scandinavia, affected the most by the PGR [36]. The vertical rates derived for Sweden are greater than those derived for the remaining areas. All 27 Swedish stations are uplifting more than 1 mm/yr, with 18 stations exceeding the value of 5 mm/yr. Moving closer to the south coast of the Baltic Sea-the Danish, German and Polish borders-the vertical rates change either values or their signs from strongly positive to slightly negative. The area of the United Kingdom is divided into two halves; the northern part is characterized by uplift, while the southern part is subsiding. Many local movements are noticed for the areas of France, Belgium and the Netherlands with vertical rates of -2 to 2 mm/yr, which we describe further in this paper. We find significant non-linearities for 62 stations, which we model with polynomial degrees from 2 to 6 (Figs. 1 and 3). No spatial dependence of those non-linearities is noticed, meaning that they only arise from the local phenomena affecting the GPS stations in the vertical direction. For each of the polynomials, their coefficients are computed. The largest coefficients are found for the second degree and vary between -1.6 and 2.2 mm/yr 2 for, respectively, the ALUK (Alūksne, Latvia) and MOBK (Obninsk, Russia) stations. The coefficients decrease for higher degrees to reach a maximum of 0.3 mm/yr 6 for the 6 th degree. The preferred polynomial degree is chosen on a station-by-station basis using the BIC values.
The seasonal signals are modelled as described by [10]. The annual (tropical) amplitudes, i.e. the period of 365.25 days, reaches a maximum of 11.2 mm for the KIR8 (Kiruna, Sweden) station. The remaining set of amplitudes fall below 10 mm. The semi-annual period is characterized by a maximum amplitude of 6.2 mm estimated for the MDVO (Mendeleyevo, Russia) station, while the remaining amplitudes do not exceed 4 mm. The draconitic year, i.e. the period of 351.6 days, is also assumed in Equation (1). Its amplitudes reach a maximum of 6.6 mm. This value is estimated for the SLD1 (  We find significant non-linearities for 62 stations, which we model with polynomial degrees from 2 to 6 ( Figures 1 and 3). No spatial dependence of those non-linearities is noticed, meaning that they only arise from the local phenomena affecting the GPS stations in the vertical direction. For each of the polynomials, their coefficients are computed. The largest coefficients are found for the second degree and vary between −1.6 and 2.2 mm/yr 2 for, respectively, the ALUK (Alūksne, Latvia) and MOBK (Obninsk, Russia) stations. The coefficients decrease for higher degrees to reach a maximum of 0.3 mm/yr 6 for the 6th degree. The preferred polynomial degree is chosen on a station-by-station basis using the BIC values.
The seasonal signals are modelled as described by [10]. The annual (tropical) amplitudes, i.e., the period of 365.25 days, reaches a maximum of 11.2 mm for the KIR8 (Kiruna, Sweden) station. The remaining set of amplitudes fall below 10 mm. The semi-annual period is characterized by a maximum amplitude of 6.2 mm estimated for the MDVO (Mendeleyevo, Russia) station, while the remaining amplitudes do not exceed 4 mm. The draconitic year, i.e., the period of 351.6 days, is also assumed in Equation (1). Its amplitudes reach a maximum of 6.6 mm. This value is estimated for the SLD1 (Saldus, Latvia) station and is an outlier, as the rest of the values do not exceed 4.0 mm. The second harmonic of the draconitic year, i.e., 175.8 days, reaches a maximum of 3.0 mm as derived for the MDVO station. The remaining values do not exceed 2 mm.

Noise Analysis
We characterize the stochastic part of the GPS position time series providing the spectral indices and the amplitudes of the power-law noise for a set of 469 European stations. Based on those values, the velocity uncertainties are computed, as indicated in Equation (4). The spectral indices vary between −1.2 and −0.3 for the set of stations we analyze, with the majority of sites falling into fractional Gaussian noise models ( Figure 5). The most extreme values are found for the VIL6 (Lomsjökullen, Sweden) and AMBE (Ambleside, UK) stations, respectively. The former has also been classified by EPN (EUREF Permanent GNSS Network) to class B stations, recognized as very noisy, see the individual EPN analysis centers solutions. In this case some spatial dependencies are noticed; stations situated in eastern Europe are characterized by lower spectral indices, meaning that they are closer to flicker noise, while the highest values of spectral index are found for the UK area and the North Sea European coast, i.e., the Netherlands and Belgium. A significant impact of the Baltic Sea is observed, which also confirms the results obtained by Klos and Bogusz [23] and Bos et al. [57]; position time series from stations situated in the coastal areas are affected by significant temporal correlation. It may arise from large-scale phenomena which were insufficiently modelled at the processing stage of the observations. Similar spatial patterns were noticed by Gruszczynski et al. [28] from the analysis based on the empirical orthogonal functions (EOFs) computed for the non-tidal environmental loading models. They found an evident impact of the Baltic Sea, concluding that this area is mostly affected by large-scale phenomena; this impact will occur as the common mode error when computed for the European sites.

Noise analysis
We characterize the stochastic part of the GPS position time series providing the spectral indices and the amplitudes of the power-law noise for a set of 469 European stations. Based on those values, the velocity uncertainties are computed, as indicated in Equation (4). The spectral indices vary between -1.2 and -0.3 for the set of stations we analyze, with the majority of sites falling into fractional Gaussian noise models ( Figure 5). The most extreme values are found for the VIL6 (Lomsjökullen, Sweden) and AMBE (Ambleside, UK) stations, respectively. The former has also been classified by EPN (EUREF Permanent GNSS Network) to class B stations, recognized as very noisy, see the individual EPN analysis centers solutions. In this case some spatial dependencies are noticed; stations situated in eastern Europe are characterized by lower spectral indices, meaning that they are closer to flicker noise, while the highest values of spectral index are found for the UK area and the North Sea European coast, i.e. the Netherlands and Belgium. A significant impact of the Baltic Sea is observed, which also confirms the results obtained by Klos and Bogusz [23] and Bos et al. [57]; position time series from stations situated in the coastal areas are affected by significant temporal correlation. It may arise from large-scale phenomena which were insufficiently modelled at the processing stage of the observations. Similar spatial patterns were noticed by Gruszczynski et al. [28] from the analysis based on the empirical orthogonal functions (EOFs) computed for the non-tidal environmental loading models. They found an evident impact of the Baltic Sea, concluding that this area is mostly affected by large-scale phenomena; this impact will occur as the common mode error when computed for the European sites.  Similarly to the spectral indices, the amplitudes of the power-law noise are also spatially dependent, with extreme values equal to 28.0 and 7.7 mm/yr −κ/4 for, respectively, VIL6 and RAL1 (Chilton, UK) permanent stations ( Figure 6). The time series from stations located in the Baltic countries are characterized by the greatest amplitudes, exceeding 15 mm/yr −κ/4 with the maximum values found for Sweden, Finland and Latvia. The smallest values are found, again, for the North Sea coast, especially for the southern part of the UK, the Netherlands and Belgium. Germany is the area with average values of power-law noise amplitudes, dividing Europe into two halves with large (East) and small (West) amplitudes of noise. Also, the UK area is divided; the northern part is characterized by amplitudes closer to 20 mm/yr −κ/4 , and larger than the estimates derived for the southern part. Similarly to the spectral indices, the amplitudes of the power-law noise are also spatially dependent, with extreme values equal to 28.0 and 7.7 mm/yr -κ/4 for, respectively, VIL6 and RAL1 (Chilton, UK) permanent stations ( Figure 6). The time series from stations located in the Baltic countries are characterized by the greatest amplitudes, exceeding 15 mm/yr -κ/4 with the maximum values found for Sweden, Finland and Latvia. The smallest values are found, again, for the North Sea coast, especially for the southern part of the UK, the Netherlands and Belgium. Germany is the area with average values of power-law noise amplitudes, dividing Europe into two halves with large (East) and small (West) amplitudes of noise. Also, the UK area is divided; the northern part is characterized by amplitudes closer to 20 mm/yr -κ/4 , and larger than the estimates derived for the southern part.

Statistical Significance of Vertical Rates
The significance of the vertical rates we derive from the GPS position time series processed by the NGL laboratory is examined. Figure 4 presents the vertical velocities estimated with the MLE algorithm for the entire European area, and only those which are recognized here as statistically significant, meaning that the value of velocity is higher than its uncertainty: its standard error computed with Equation (4), Figure 7. All of the Scandinavian, Danish and Estonian vertical rates are found here as significant (see Supplementary Materials 1) and are left for further analyses. Due to their small vertical rates, 172 velocities from central and western Europe are found to fall below their uncertainties. These stations are therefore excluded from further analyses. The uncertainties of vertical velocities range here between 0.1 and 3.9 mm/yr. The largest uncertainties are found for sites affected by the lowest spectral indices and the greatest amplitudes of power-law noise; i.e. for the Baltic Sea area. These two parameters, when combined together in Equation (4), result in a high velocity variance. Naturally enough, the stations situated within the North Sea areai.e. the UK,

Statistical Significance of Vertical Rates
The significance of the vertical rates we derive from the GPS position time series processed by the NGL laboratory is examined. Figure 4 presents the vertical velocities estimated with the MLE algorithm for the entire European area, and only those which are recognized here as statistically significant, meaning that the value of velocity is higher than its uncertainty: its standard error computed with Equation (4), Figure 7. All of the Scandinavian, Danish and Estonian vertical rates are found here as significant (see Supplementary Materials Table S1) and are left for further analyses. Due to their small vertical rates, 172 velocities from central and western Europe are found to fall below their uncertainties. These stations are therefore excluded from further analyses. The uncertainties of vertical velocities range here between 0.1 and 3.9 mm/yr. The largest uncertainties are found for sites affected by the lowest spectral indices and the greatest amplitudes of power-law noise; i.e., for the Baltic Sea area. These two parameters, when combined together in Equation (4), result in a high velocity variance. Naturally enough, the stations situated within the North Sea area-i.e., the UK, Belgium, the Netherlands and Denmark-are all characterized by small velocity uncertainties, which resulted in the fact that a large number of the vertical rates are recognized as significant (Figure 4a). Few of the English and Belgian sites are found to have large velocity uncertainties, although both the spectral index and the amplitude of power-law noise are small enough to result in a low velocity uncertainty. These sites may have been affected by local phenomena which occurred as unmodelled events in the stochastic part and therefore leads to overestimation of the velocity uncertainties. All other western European vertical rates are characterized by uncertainties lower than 0.6 mm/yr, as computed with the power-law and white noise combination. The Swedish stations, despite having large velocity uncertainties, are still recognized as statistically significant since the values of vertical rates are large enough.
despite having large velocity uncertainties, are still recognized as statistically significant since the values of vertical rates are large enough.
Next, the uncertainties of velocity estimates with considering a combination of power-law and white noise are compared with the white-noise-only assumption. The latter imposes that a covariance matrix of observations includes only a white noise term, as in Equation (2). Therefore, as shown before, the velocity uncertainties are underestimated. In our research, the general dilution of precision (GDP, [32]) of velocity uncertainty is equal to eight times on average. The largest ratios are found, as might be expected, for the Baltic countries, for which the spectral indices are far from 0 (white noise).
Taking all of the above into consideration, we conclude that any of the time series has to be reliably modelled, i.e. all of the time series components should be included to decide on the statistical significance of vertical rates and to recognize those which can be interpreted in terms of crustal deformations.  Next, the uncertainties of velocity estimates with considering a combination of power-law and white noise are compared with the white-noise-only assumption. The latter imposes that a covariance matrix of observations includes only a white noise term, as in Equation (2). Therefore, as shown before, the velocity uncertainties are underestimated. In our research, the general dilution of precision (GDP, [32]) of velocity uncertainty is equal to eight times on average. The largest ratios are found, as might be expected, for the Baltic countries, for which the spectral indices are far from 0 (white noise).

Comparison to the GIA model
Taking all of the above into consideration, we conclude that any of the time series has to be reliably modelled, i.e., all of the time series components should be included to decide on the statistical significance of vertical rates and to recognize those which can be interpreted in terms of crustal deformations.

Comparison to the GIA Model
The GPS-derived vertical rates are compared to the present-day uplift predicted by the ICE-6G model (Figure 8). The differences between both reach 0.5 mm/yr on average. Sweden, for which the vertical rates agree within 0.2 mm/yr with the predicted uplift, is the area of best agreement. The largest difference between GPS and PGR rates of 6.7 mm/yr is found for the NYAL (Ny-Ålesund, Norway) station, which has been also classified by EPN as noisy class B station. This site is extremely affected by the elastic uplift of present-day ice mass loss which is not represented by the GIA model. As shown by Rajner [58], this difference is reduced when the present-day ice melting and the deglaciation of the little ice age are accounted for. The maximum of the GPS vertical rates of 11.0 mm/yr is found for the SKE0 station and is situated between the present-day centers of uplifts as predicted by the ICE-6G and ICE-5G GIA models, respectively. The ice footprints at the last glacial maximum (LGM) agree for both GIA models, whereas the hinge lines computed for both of them are shifted with respect to each other. Worth noting is the fact that the hinge line computed for the newest ICE-6G model is much closer to what the GPS vertical rates provide.
maximum. This ice footprint almost overlaps the hinge line. It is only somehow shifted to the south for central Europe-i.e. Poland and Germany-and for the British Isles ( Figure 9).
Then, the GPS-derived vertical velocities are interpolated to the continuous vertical velocity field (Figure 9). We focus on the areas surrounding the hinge line, as the agreement between Scandinavian GPS uplifts and those predicted by the GIA model is quite good due to significant vertical rates. The hinge line area is the most interesting one, as the vertical rates estimated from observations may help to assess the sensitivity of the GIA model. Unlike in central Europe, the eastern areas agree well with the hinge line. If this line was plotted for the Polish and German territory based on the GPS-derived vertical rates, it would be shifted to the south, compared to what we obtain from the GIA model predictions. The central Europe area however stays in a good agreement with the ice footprint, meaning that the change in the sign of vertical rate from positive to negative occurs at the edge of the area entirely covered by ice during the LGM. The British Isles uplift of 1-2 mm/yr in the north while slightly subsiding in the south by 1-2 mm/yr is in good agreement with the hinge line in England, but disagreeing in Ireland. The area of the Netherlands subsides significantly with the average vertical rate of -2 mm/yr. The area of Belgium is affected by many local phenomena, showing locally strong subsidence and strong uplift of a few mm/yr. The present-day uplift predicted by the ICE-6G model is employed to compute the hinge line. Looking at eastern Europe, this hinge line goes through Russia and transects the area of Belarus, Poland, Germany and Denmark ( Figure 9). Then, it divides the area of the British Isles into two halves: the northern and southern parts. The areas situated to the north of this line are affected by uplift due to the PGR. The southern area might subside slightly, but the vertical rate is close to zero in most cases, however being statistically significant (see also Figure 2). Furthermore, we consider the ice footprint at the time of 21 Ka, meaning the area entirely covered with ice at the last glacial maximum. This ice footprint almost overlaps the hinge line. It is only somehow shifted to the south for central Europe-i.e., Poland and Germany-and for the British Isles ( Figure 9).
Then, the GPS-derived vertical velocities are interpolated to the continuous vertical velocity field (Figure 9). We focus on the areas surrounding the hinge line, as the agreement between Scandinavian GPS uplifts and those predicted by the GIA model is quite good due to significant vertical rates. The hinge line area is the most interesting one, as the vertical rates estimated from observations may help to assess the sensitivity of the GIA model. Unlike in central Europe, the eastern areas agree well with the hinge line. If this line was plotted for the Polish and German territory based on the GPS-derived vertical rates, it would be shifted to the south, compared to what we obtain from the GIA model predictions. The central Europe area however stays in a good agreement with the ice footprint, meaning that the change in the sign of vertical rate from positive to negative occurs at the edge of the area entirely covered by ice during the LGM. The British Isles uplift of 1-2 mm/yr in the north while slightly subsiding in the south by 1-2 mm/yr is in good agreement with the hinge line in England, but disagreeing in Ireland. The area of the Netherlands subsides significantly with the average vertical rate of −2 mm/yr. The area of Belgium is affected by many local phenomena, showing locally strong subsidence and strong uplift of a few mm/yr. Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 17 Figure 8. The vertical rates derived from the NGL position time series, plotted in mm/yr as colored dots. The present-day vertical rates predicted by the ICE-6G GIA model are given as the background colored map. The present-day center of uplift of the ICE-6G GIA model is plotted as a black square with a vertical rate of 9.7 mm/yr. The present-day center of uplift predicted by the previous GIA model, i.e. ICE-5G, is plotted as a green square with a vertical rate of 9.3 mm/yr. The GPS-derived vertical rate estimated for the southern coast of the Baltic Sea-i.e. Estonia, Latvia, Lithuania, Poland, Germany, Denmark and Sweden-is systematically underestimated by the GIA prediction by up to 2 mm/yr ( Figure 10). This means that the GIA uplift should probably not be employed alone and adapted in those areas as the vertical land motion to correct the tide gauge records and receive the absolute sea-level values. It should be definitely supported by other observations. This underestimation reaches central Germany, for which the GIA model predicts a vertical rate close to zero, while the GPS observations point to evident uplift of up to 1.5 mm/yr. This The GPS-derived vertical rate estimated for the southern coast of the Baltic Sea-i.e., Estonia, Latvia, Lithuania, Poland, Germany, Denmark and Sweden-is systematically underestimated by the GIA prediction by up to 2 mm/yr ( Figure 10). This means that the GIA uplift should probably not be employed alone and adapted in those areas as the vertical land motion to correct the tide gauge records and receive the absolute sea-level values. It should be definitely supported by other observations. This underestimation reaches central Germany, for which the GIA model predicts a vertical rate close to zero, while the GPS observations point to evident uplift of up to 1.5 mm/yr. This is probably due to the ice load history in that region. The GIA model overestimates the uplift for the northern UK-the GPS rates are systematically lower by about 0.5 mm/yr. The central part of the UK is underestimated by the GIA predictions by up to 0.5 mm/yr, while the southern part does not show any systematic spatial pattern: the individual locations are under-or over-estimated, which evidences many local phenomena. The areas close to London are systematically overestimated, meaning that the ground is subsiding more than is expected from the GIA predictions by more than 1 mm/yr. According to the GIA predictions, southern Ireland should not be affected by any crustal movement, while the GPS observations prove that it is uplifting by around 0.5 mm/yr. The area of the Netherlands and the coastal area of Belgium are both subsiding more than is predicted by the GIA model of around 1 mm/yr. The inland part of Belgium, Luxembourg and the western part of Germany show strong positive velocities when compared to the GIA model. Most of these stations' uplift of more than 1 mm/yr is probably caused by present-day elastic deformation due to hydrology, especially in Rhine basin, or non-tidal atmospheric loading. is probably due to the ice load history in that region. The GIA model overestimates the uplift for the northern UK-the GPS rates are systematically lower by about 0.5 mm/yr. The central part of the UK is underestimated by the GIA predictions by up to 0.5 mm/yr, while the southern part does not show any systematic spatial pattern: the individual locations are under-or over-estimated, which evidences many local phenomena. The areas close to London are systematically overestimated, meaning that the ground is subsiding more than is expected from the GIA predictions by more than 1 mm/yr. According to the GIA predictions, southern Ireland should not be affected by any crustal movement, while the GPS observations prove that it is uplifting by around 0.5 mm/yr. The area of the Netherlands and the coastal area of Belgium are both subsiding more than is predicted by the GIA model of around 1 mm/yr. The inland part of Belgium, Luxembourg and the western part of Germany show strong positive velocities when compared to the GIA model. Most of these stations' uplift of

Summary and conclusions
We present a comprehensive analysis of the GPS vertical position time series which should be employed to analyze the geodynamic phenomena known as the post-glacial rebound. Our methodology is validated for the set of 469 European stations which were processed by the NGL Figure 10. The differences between the GPS-derived vertical rates and the uplift predicted by the ICE-6G GIA model at the GPS station coordinates. The hinge line, constructed based on the present-day GIA uplift, and the ice footprint at the time 21 Ka are plotted with the dashed and green lines, respectively. Individual maps show eastern Europe (a), the British Isles (b) and western Europe (c).

Summary and Conclusions
We present a comprehensive analysis of the GPS vertical position time series which should be employed to analyze the geodynamic phenomena known as the post-glacial rebound. Our methodology is validated for the set of 469 European stations which were processed by the NGL laboratory. We prove that the proper treatment of either deterministic or stochastic components of the time series is indispensable to obtain reliable vertical velocities along with their uncertainties. The non-linearity of trend and seasonal signatures are critically important, as they bring significant correlation to the stochastic part, if unmodelled. Also, the assumptions of the noise character will influence the velocity uncertainties and lead to their underestimation of up to 5 mm/yr [59], up to 6 times [54], or 8 times on average (this study).
Similarly to Klos and Bogusz [23], and Gruszczynski et al. [28], we also observed a distinguished spatial pattern in the distribution of the parameters describing the stochastic part, i.e., the spectral index and the amplitudes of the power-law noise with an impact for the Baltic countries. These parameters influence the velocity uncertainty, as in Equation (4). The largest uncertainties are observed in this study for Sweden, Latvia, Estonia and Lithuania, for which the greatest parameters of the stochastic part are delivered. Based on these uncertainties, we decide about the statistical significance of the GPS vertical rates. As shown by [28] with empirical orthogonal functions, areas of central and northern Europe, i.e., especially those surrounding the Baltic coasts, are significantly impacted by non-tidal atmospheric loading, which is for now not included within the GPS processing. We argue that this loading may mostly affect the character of the position time series, i.e., the amplitudes of power-law noises, and also, the character of the noise, and the velocity uncertainties, consequently. The impact of non-tidal atmospheric loading is removed these days on the stage of data processing using different modelled values.
The GPS vertical velocities delivered in this research reach the maximum values for Scandinavia, as expected. The maximal uplift is found for the SKE0 station and is equal to 11.0 mm/yr. This station has been also found by EPN to significantly uplift by 10 mm/yr. The differences between both rates may arise from different assumptions made in Equation (1) during the vertical rate estimates and from different lengths of time series as well. This station is situated between two previously estimated centers of the ICE-5G and ICE-6G GIA models-the cities of Umea and Furuogrund, for which with the vertical rate exceeds 9 mm/yr [36]. Moving to the south, the vertical velocities decrease to reach values close to zero at the GIA hinge line.
Finally, we provide an extensive comparison between the GPS-derived rates and the newest GIA ICE-6G model for a dense European network of GPS stations, which had not been examined before. We prove that these rates agree at a 0.5 mm/yr level on average with some areas being clearly overand under-estimated by the GIA predictions. For these, the GIA uplift should not be employed as the representative land motion to e.g., transfer the relative into the absolute sea-level rise estimated from the tide gauge records.
The methodology we present is also useful to show areas affected by local phenomena for which the elastic uplift exceeds the post-glacial rebound arising from the glacial isostasy. Such an elastic uplift due to present-day ice melting or hydrology has already been modelled by van der Wal et al. [60], Simon et al. [56], Nocquet et al. [61] and Rajner [58]. It is expected the present-day melting of glaciers and hydrology contribute to the elastic uplift estimated from the vertical rates from 0.5 to 2 mm/yr.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/10/1209/s1, Table S1: A set of 469 European GPS permanent stations with latitudes higher than 50 • N which observations are processed by the Nevada Geodetic Laboratory in a PPP mode.