Introduction

GNSS velocity fields have been derived from position time series shortly after the beginning of the Global Navigation Satellite Systems (GNSS) era, especially after continuous recordings were started at permanent stations. Every year has brought new insights into the technical capabilities of GNSS and geophysical phenomena interpreted from the data. Reprocessing of the historic data has been carried out in many projects to take advantage of improved orbit and clock parameters, models, and methods. One of the major efforts has been carried out by the International GNSS Service (IGS, Rebischung et al. 2016).

The glacial isostatic adjustment (GIA) and its effects on the national reference frames have been the main driving force for dense GNSS position and velocity solutions in Fennoscandia. The first GNSS-derived velocities have been estimated under the Baseline Inferences for Fennoscandian Rebound, Sea-level, and Tectonics (BIFROST) project (Scherneck et al. 1998) and regularly updated every few years. Kierulf et al. (2014) published the latest and partly densified solutions compared to the original BIFROST network (Lidberg et al. 2010). All these solutions include only a few stations from the Baltic countries, but an equal coverage of stations in the Baltics as in the Nordic countries would be important for Fennoscandian GIA modeling (Wu et al. 2010). These studies have also been focused on GIA modeling, and they did not provide a station position solution. To produce a homogenous and high-quality GNSS solution covering all the Nordic and Baltic countries, the Nordic Geodetic Commission (NKG) launched its GNSS Analysis Centre in 2012. It builds on distributed and highly homogeneous processing and analysis of national reference stations (Lahtinen et al. 2018). Each national analysis center (AC) processes its subnet, and the subnet solutions are combined into daily solutions.

Realistic uncertainties of the estimated velocities are necessary when using the data as a constraint for the GIA modeling. It is well known that the GNSS position time series include temporal correlated noise, and the assumption of white noise produces overly optimistic uncertainties (Agnew 1992). There are several methods for estimating the noise in a GNSS time series and constant discussion on the type and amount of the noise. Most of the studies have focused on power-law (PL) forms of noise: \(p \sim f^{n}\), where the \(f\) is the frequency and \(n\) is the spectral index. White noise (WN, \(n = 0\)), flicker noise (FN, \(n = - 1\)), and random walk (RW, \(n = - 2\)) are special cases of the power-law noise. In particular, the FN + WN and PL + WN combinations, where spectral index is close to flicker noise, are supported by many studies (Williams et al. 2004; King and Williams 2009; Santamaría-Gómez et al. 2011, Amiri-Simkooei 2016; Klos et al. 2016a). The origin of the flicker noise has been addressed in the GNSS system and related methods (equipment, modeling deficiencies, site-specific issues, etc.), although it is difficult to distinguish the sources.

The amount of data used in the noise analysis studies have gradually increased during the years. Santamaría-Gómez et al. (2011) demonstrated that the noise content is dependent on the time period of the analysis and that it is probably related to the amount and quality of the recorded data. The old data were noisier. The presence of random walk noise due to monument motion was suggested in the early years (Johnson and Agnew 1995), although there are several studies (Beavan 2005; King and Williams 2009; Klos et al. 2016b) that do not support it. However, Langbein (2012) demonstrated the differences between the PL and FN + RW models, and how the ignorance of RW would significantly underestimate the uncertainties. Later, Dmitrieva et al. (2015) pointed out the difficulty of estimating small parts of RW and the resulting significant changes in uncertainties. They suggested a network method, where all stations were analyzed together instead of independently, if the individual station uncertainties are not the goal. Further, Dmitrieva et al. (2016) also showed how the noise results vary significantly depending on whether or not the trend is estimated simultaneously with the noise parameters. The FN + WN and PL + WN, where spectral index \(- 2 < n < - 1\), were relatively insensitive to de-trending, whereas the presence of RW was difficult to either confirm or reject if FN and WN exist.

In this study, we present a new ITRF2014 densification for the Nordic and Baltic countries. We describe the methodology of the combined position and velocity solution including time-series analysis. We estimate uncertainties using the most cited combinations of power-law noise models. We aim to find a model that quantifies realistically the station velocity uncertainty. Finally, we discuss the quality of the presented solution.

Methods

The estimation of a GNSS velocity field consists of several phases. First, we reprocess the historic GNSS data into daily position solutions. Then we analyze the time series of individual stations to find out discontinuities and to exclude periods of bad data. The final velocity estimation can be done either by combining all daily solutions using least-squares methods or by station-wise analysis. We describe the input data, models, and methods used to produce the densified velocity field and the related uncertainties.

Data

The time span of the reprocessing was set to 1997–2017 corresponding to GPS weeks 887–1933. In practice, each national AC started processing from the beginning of their data history. The Nordic countries started between 1997 and 2001, and the Baltic countries between 2006 and 2008. Figure 1 shows the total coverage of the network and the station density through the years. The network includes some stations with a remote location in respect to the Nordic and Baltic area, as they belong to the operative NKG solution of the EUREF Permanent GNSS Network (EPN) used as a backbone solution. Most of the Nordic stations have been established on bedrock, whereas many of the Baltic ones are on roof tops due to geological features, e.g., lack of bedrock.

Fig. 1
figure 1

An overview of station density and coverage. Black circles represent national stations, green circles represent the station in the original ITRF2014 solution, and orange circles represent the used datum points

Full size image

The historic GNSS data were reprocessed to produce consistent solutions for the full time span. We mainly used the same models and methods as described for the operational work (Lahtinen et al. 2018) following the EPN guidelines for its ACs. First, each national AC processed daily solutions of their own subnet, including common stations outside the country. Additionally, we processed the daily backbone solutions covering the whole area. All daily subnet solutions were double-difference solutions computed using Bernese GNSS Software (Dach et al. 2015). We used satellite orbits, earth orientation parameters, and ionosphere models from the Center for Orbit Determination in Europe (CODE). The REPRO_2013 products were used until 2014.0 and the CODE’s final products thereafter (Dach et al. 2018). We applied individual antenna calibrations if available, and if not, the igs08.atx antenna model. The subnet solutions were aligned with the IGb08, as the ITRF2014 was not available then, using minimum constraints over no-net-translation (NNT) condition on the IGS station coordinates at the observation epoch. The NNT method is recommended for regional networks by the Bernese GNSS Software, as the fixed satellite orbits define the orientation of the network. We used only GPS observations to avoid possible effects of the changing GLONASS constellation during the years. Each AC processed two solutions: the main solution with 3-degree cutoff angle that is recommended by EPN, and an alternative 10-degree solution. It was motivated by the fact that the actual elevation angle may have varied due to changes at the site surroundings over the long time span mainly due to tree growth, and many sites do not have visibility lower than the 10-degree angle.

Next, the subnet solutions were combined into daily solutions of the full network. The different datum definition possibilities were analyzed by Lahtinen et al. (2018), and based on that, we applied minimum constraints over all seven parameters (translations, rotations, and scale) using CATREF software (Altamimi et al. 2018). The need of estimating a scale parameter can be questioned, as there should not be a significant scale difference between the frame and our solutions, but our tests have shown very small effects of it. Misfitting station solutions were searched and rejected stepwise down to the criteria of (5, 5, 10) mm in (north, east, up; NEU) and having a normalized residual larger than six. The limit for normalized residuals prevents too low rejections in the case of larger discrepancies in station solutions. The misfitting datum points were rejected down to (15, 15, 30) mm in NEU. The limits were set so that the amount of rejected data was kept reasonable.

Pre-analysis of the time series

GNSS position time series are traditionally modeled with a linear trend plus seasonal harmonics. However, the time series typically contains some periods of outliers, discontinuities, and/or nonlinearities, which possibly affect the estimated positions and velocities. In particular, in the Nordic countries, snow typically accumulates on the antenna radomes during wintertime, causing systematic series of outliers in the time series. The nonlinearities may be originating from geophysical processes, such as earthquakes, post-seismic relaxation, contemporary ice-melting-induced GIA, or from site-specific issues like monument stability. The discontinuities of the time series have the largest effect on the final velocity estimates. In particular, if the time series does not strictly behave linearly, the effect of a single discontinuity may be larger than the realistic uncertainty. Much effort has been put into automatic discontinuity detection (Borghi et al. 2012; Vitti 2012; Bruni et al. 2014; Kowalczyk and Rapinski 2018). However, Gazeaux et al. (2013) provided simulated data to be tested with automatic approaches, but the manual discontinuity detection produced significantly higher performance.

As the number of stations in each country was reasonable, we chose to do the pre-analysis manually. Each AC did the pre-analysis of their national stations, and the remaining non-Nordic/Baltic stations were distributed among the ACs. The removal of the snow and other short and deviating periods were marked with the help of the Tsview tool (Herring 2003) and rejected from the daily solution before estimating velocity.

The station log files or similar information was the primary input for the discontinuity analysis. Additionally, other sources such as multipath plots and site pictures were utilized in case of a questionable discontinuity. To harmonize and validate the pre-analysis, the results were cross-checked before and after with a subset of stations. At some sites, gradual tree growth above 10 degree and a sudden cutting of trees have been seen as a jump in the time series, especially in the up component. Figure 2 shows an example time series with roughly a 10-mm jump in height after the tree-cut. In those cases, we tried to reject data from all affected years. The start and end of the period were approximated case-by-case based on the existing information, such as nonlinearities in the time series, data quality information, tree-cut dates, and finally the effect of the rejection on the trend.

Fig. 2
figure 2

An example of a tree-cut effect on the up component at VAAS station. The tree-cut date is shown in red-dashed line

Full size image

Velocity estimation

The station velocities can be estimated by either combining the daily solutions with full covariances using least-squares methods or estimating trends from the station position time series individually for each component (NEU). The analysis of single station time series is fast, and there are several tools for this. The linear trend, sizes of shifts at discontinuities, and magnitudes of annual, semi-annual, and other seasonal signals can be analyzed. However, there are some important advantages in the least-squares combination. By combining the daily solutions with full covariance matrices, the spatial correlations between the stations can be taken into account, although it is computationally heavy. The velocities of the stations located at the same site can be constrained. In addition to the velocities, positions for each station are estimated. Finally, the full network can be aligned to the reference frame using minimum constraints on both station positions and velocities.

We selected the least-squares combination approach to produce both positions and velocities and to estimate site velocities for the co-located stations. We used CATREF software. The combination model is based on seven-parameter similarity transformation, where each of the daily solutions is transformed into the combined or stacked solution at the fixed epoch (Altamimi et al. 2016, 2018). We aligned our solution to the IERS release of the ITRF2014 (Altamimi et al. 2016). We used 2010.0 as the combination epoch, that is, the same as in ITRF2014. Each discontinuity introduces a new position solution for the station. Therefore, one station may have several cumulative position solutions, which we refer to as positions solutions. We used all available ITRF2014 stations as datum points of which position solutions were consistent with our discontinuity analysis (Fig. 1). The remaining outliers were screened based on the size of the residual and/or standardized residual. We did not add seasonal terms into the model. As the combination model includes the transformation parameters for each daily solution with respect to the combined frame, each frequency would have introduced seven singularities in the normal equation system (Altamimi et al. 2016), and treatment of them in the regional network is not obvious. The seasonal terms would also have increased the need for computational resources, and in the long time series, their effect on trends is minimized.

Uncertainties

We estimated the velocity uncertainties using the Maximum Likelihood Estimation (MLE) method implemented by Bos et al. (2013) in Hector software. We used the residuals of the station positions of the cumulative solution as an input instead of original position time series, as the residuals describe the uncertainty of our final velocity estimates. We modeled a linear trend, annual and semi-annual signals, and noise model parameters. The estimated linear trends were insignificant, as we used the residual time series as the input. There were six different noise models implemented in Hector, and any combination of those can be estimated using the MLE technique. We analyzed the following three most cited combinations of power-law-formed alternatives to find out the generally best-fitting model to our time series: Power-law + White noise (PL + WN), Flicker + White noise (FN + WN), and Random Walk + Flicker + White noise (RW + FN + WN).

We introduced the same discontinuities for the residual time series that were used in the cumulative combination in CATREF. A station with many discontinuities should get a larger uncertainty compared to continuous observations. Otherwise, the uncertainties of those stations would become too optimistic, as if we had discontinuity-free time series.

The uncertainties of the co-located and constrained stations were finalized separately, as Hector can process only one station at a time. To derive a formula for the total uncertainty of the constrained stations, we approximate the final velocity as a weighted mean of the two velocities. Then the final uncertainty can be computed using the error propagation law. The weighting could be done based on the length of the time series or the estimated single stations’ velocity uncertainties. The weighting approach did not affect much on the total uncertainty. The rms of the differences between the two alternative weighting approaches was 0.05 mm/year, so we used the time-series lengths as weights. Thus, by using the error propagation law, the total uncertainty was computed by

$$\sigma_{\text{tot}} = \sqrt {\left( {\frac{{t_{i} }}{{t_{i} + t_{j} }}} \right)^{2} \sigma_{i}^{2} + \left( {\frac{{t_{j} }}{{t_{i} + t_{j} }}} \right)^{2} \sigma_{j}^{2} }$$

where \(\sigma\) is the uncertainty and \(t\) is the number of epochs for the station pair \(i\) and \(j\).

Results

We analyze the results of the combined cumulative 3-degree solution and present the final velocity field for the Nordic and Baltic countries. We describe the results of the noise analysis of residual time series and present the final uncertainties. We also compare the 3-degree solution to our alternative 10-degree solution and to the other published velocity solutions.

Cumulative position and velocity solution

Our cumulative solution included 289 stations. We excluded 37 stations from these results, as they had less than 1095 epochs in the final combined solution corresponding to 3 years of data. A few stations with a shorter time span were kept if they were constrained to a co-located station. The final products included station positions and velocities for 252 stations.

Figure 3 shows the weighted RMS (WRMS) of the daily solutions, based on the station-wise residuals weighted by the corresponding sigmas. We reached a 1-mm level for horizontal components and 3-mm level for the up component, verifying good internal precision and long-term stability of the daily solutions. The scatter is slightly larger for the first years as expected.

Fig. 3
figure 3

The weighted RMS of the daily solutions for north (top), east (middle), and up (bottom) components shows the internal stability of the solutions through the years

Full size image

We used 22 stations in total for the datum definition. We achieved a good agreement with the ITRF2014 reference frame: The WRMS were 0.4, 0.7, 1.6 mm for the positions and 0.1, 0.1, 0.2 mm/year for the velocities in NEU, respectively.

The time-series analysis showed significant velocity differences at some of the Swedish co-located stations. The old stations on concrete pillar monuments have been established in the 1990s and the newer ones on steel grid masts in 2011. Both stations are still operational, and they are located only a few meters apart. The difference was significant, especially in horizontal components, where typically 5–6 years of data should give converging estimates from our experience. Constraining velocity solutions of those stations would have been a trade-off between the old and new one and that was visible as trends in residual time series. Therefore, we removed the constraints of all 20 Swedish co-located station pairs, although the difference was not significant for all station pairs. The final horizontal velocity differences ranged from 0.11 to 0.36 mm/year, and the first and third quartiles ranged from 0.22 to 0.29 mm/year. We did not observe a similar significant difference at the other co-located stations, but the overlapping time span was typically shorter in those cases.

Figures 4 and 5 show the horizontal and vertical velocities. The Eurasian and North American plate tectonics dominate the horizontal velocities, and the Fennoscandian land uplift dominates the vertical velocities. Iceland stands out with large local velocity differences due to its location on both plates. The numerical dataset of the NEU velocities is given in electronical supplement.

Fig. 4
figure 4

The densified ITRF2014 horizontal velocities mostly reflect the plate tectonics. The color of the circles shows the magnitude of the velocity, and the vectors show the direction. The vectors are shown only for a part of the stations for clarity. The plate border between the Eurasian and North American plates is shown with a gray line

Full size image
Fig. 5
figure 5

The densified ITRF2014 vertical velocities show the land uplift pattern in the Fennoscandia. The plate tectonic-related effects cause large uplift rates at some of the stations in Iceland. Smaller circles have been used for the newer Swedish co-located stations that were not constrained

Full size image

Uncertainty analysis of the velocities

The analysis of the different noise models showed noise characteristics between white and flicker noise. Table 1 summarizes the statistics of the flicker noise fractions of the FN + WN solution and the spectral indices of the PL + WN solution. The flicker noise fractions describe the percentage of the flicker noise of the total noise. The mean flicker noise fractions of all components ranged between 43% and 78% (first and third quartiles, respectively), and the spectral indices between − 0.7 and − 0.9 (first and third quartiles). The east component was closest to the flicker noise, whereas the north and up components are somewhere between flicker and white noise. The RW + FN + WN combination produced mostly the same results as the FN + WN solution, as the estimated fractions of the RW model were less than one percent for roughly 90 percent of the stations. The median increase was 0.01, 0.01, 0.04 mm/year in NEU. Typically, the larger fractions or larger amplitudes for RW were estimated for the stations, of which time series were either close to the minimum 3 years or had non-standard behavior such as nonlinear, short pieces, or unknown shifts.

Table 1 Statistics of the FN fractions (FN + WN solution) and spectral indices (PL + WN solution)
Full size table

There appeared to be visually some spatial correlation in the spectral indices of horizontal components in the Fennoscandian area, most visible in the east component (Fig. 6, top panel). Stations in the eastern part, covering eastern Sweden, Finland, and Baltic, are close to flicker noise, whereas the indices are on a − 0.5 level in the western part, covering western Sweden, Norway, and Denmark, independent of the time-series length. The spectral indices of the up component were more random, without visual spatial correlation (bottom panel). We do not see any obvious cause for spatially different spectral indices in horizontal components.

Fig. 6
figure 6

The visualized estimated spectral indices show some spatial correlation, especially in east component (top). The pattern in up component (bottom) is more random. The size of the circle changes with respect to the length of the time series

Full size image

Table 2 summarizes the statistics of the estimated FN + WN and PL + WN uncertainties. The results are very similar, that is, in line with the spectral analysis. The quartiles are on the same level within ± 0.02 mm/year both horizontally and vertically. Mostly the Akaike and Bayesian information criteria (AIC/BIC), describing the goodness of the fit of the model, were very close for both noise model alternatives not giving a clear support for either one of them. However, the vertical uncertainties are more sensitive to the choice of the noise model, especially in the case of shorter time series. In a few cases, the PL + WN seems to overestimate the uncertainties. For example, power-law gives vertical uncertainties of 0.7–0.8 mm/year for some stations with roughly 4–5 years of data without any discontinuities or issues, whereas the flicker noise estimates were only 0.2 mm/year. In these cases, the AIC/BIC estimators support the FN + WN alternative. The results suggest that the FN + WN choice is more robust to be used for the whole dataset, especially when the length of the time series varies.

Table 2 Statistics of estimated uncertainties in mm/year for FN + WN and PL + WN solutions
Full size table

Figure 7 shows the final FN + WN uncertainty estimates. The uncertainties reflect well the length of the time series and the overall quality of the time series in the area. The uncertainties are less than 0.05 mm/y horizontally and 0.10 mm/year vertically for the best stations with almost 20 years of solutions. The introduction of the breaks also in the uncertainty estimation was a key method to achieve more realistic uncertainties, e.g., for the Icelandic stations with geodynamical issues. On the other hand, it did not affect significantly the uncertainty estimates of stable stations. The uncertainty estimates describe only the noise of the time series excluding, e.g., the error sources originating from the reference frame itself.

Fig. 7
figure 7

The visualized final (FN + WN) horizontal (top), and vertical (bottom) uncertainties for the area of main interest. Note the different scales in the top and bottom panels

Full size image

Comparisons of the velocities and uncertainties

The comparison of the velocity field to other solutions gives insight into the quality and stability of our solution and the existing solutions. We show the differences between our 3-degree main solution and the alternative 10-degree solution. We also analyze the velocity differences at the ITRF2014 stations that were excluded from the datum points. Finally, we show the differences to the latest published Nordic–Baltic velocity field.

3- and 10-degree solutions

We combined the 10-degree solution parallel to the 3-degree main solution. Figure 8 shows the distribution of the differences. Most of the differences (75%) were smaller than 0.04, 0.02, 0.15 mm/year in NEU, respectively. For the subset of stations having more than 10 years of data, the corresponding figures were 0.03, 0.02, 0.11 mm/year. The horizontal differences can be considered negligible, and the vertical differences correspond well to the 1-sigma of the FN + WN uncertainties. Figure 9 demonstrates that the largest differences were linked to the short time series, and there are no geographical areas with systematically larger differences. However, even for few stations having more than 10 years of data, a vertical difference of 0.25 mm/year was reached, suggesting some station-specific issues. Those were EPN stations outside Nordic and Baltic countries.

Fig. 8
figure 8

Histogram of the velocity differences between 3- and 10-degree solutions

Full size image
Fig. 9
figure 9

The velocity differences between 3- and 10-degree solutions were very small and non-systematical. Station-specific issues, such as obstacles at low elevations, may cause the observed differences. The colored circles represent the vertical differences, and the vectors represent the horizontal differences. The size of the circle changes with respect to the length of the time series

Full size image

When comparing the velocity differences to the uncertainties, we found five stations with a larger difference than the 3-sigma uncertainty (FN + WN) in some component. These cases mostly reflect the overly optimistic uncertainties of those stations. The horizontal uncertainties were as low as 0.01 mm/year and vertical uncertainties 0.03 mm/year. A few cases can be considered as bad-quality stations due to the strong dependence of the velocity on the elevation cutoff setting and to be excluded from the geodynamical modeling.

The comparison shows the high stability of the solutions. As we did not observe any deficiencies or systematics, we consider the 3-degree main solution as the final solution.

Original ITRF2014 solution

Some of the ITRF2014 stations were rejected from the datum as we had modeled the stations differently, either by rejecting more data or by adding/removing discontinuities based on the local knowledge of the actions at the stations. Table 3 summarizes the differences between the estimated velocities. For TRO1 and the up component of VAAS, the differences are clearly larger than the 3-sigma of the FN + WN uncertainty estimates. The ITRF2014 GNSS time-series cover years 1996.0–2015.1, and the time span difference may somewhat affect the differences. The differences emphasize the importance of the time-series analysis before the final velocity estimation.

Table 3 Differences in velocities between ITRF2014 and our results
Full size table

Latest Nordic–Baltic velocity field

We had 80 points in common with the latest published Nordic–Baltic BIFROST velocity field by Kierulf et al. (2014). The BIFROST solution was aligned to our solution using a similarity transformation with translation and rotation parameters before comparison to reduce the differences in datum alignment. The amount of data is approximately 4 years less in the BIFROST solution. The exact time span has not been reported. Figure 10 shows the distribution of the horizontal and vertical velocity differences, and Fig. 11 visualizes them. The largest horizontal differences were clearly linked to the stations with short time series, reaching a level of 0.5 mm/year and exceeding significantly the uncertainty estimates of the BIFROST solution. Half of the vertical differences were within ± 0.15 mm/year that corresponds to the 1-sigma level of our uncertainty estimates. The figure shows large vertical differences at single stations, e.g., in Northern Finland, but the BIFROST uncertainties of those stations are also large, roughly 0.5 mm/year, compared to our results. We achieved an uncertainty level of 0.10–0.15 mm/year for the same stations. Both uncertainty estimates are FN + WN-based estimates and thus comparable. The differences demonstrate the need of an up-to-date solution for the Nordic and Baltic countries.

Fig. 10
figure 10

Histogram of the differences between our solution and the aligned BIFROST solution

Full size image
Fig. 11
figure 11

The station velocity differences between our solution and the aligned BIFROST solution show partly good agreement of the solutions, and partly differences larger than 0.5 mm/year. The symbols are the same as in Fig. 9. The black and white circles represent values exceeding ± 0.8 mm/year

Full size image

Discussion

Our solution is a robust position and velocity estimate for the whole area of the Nordic and Baltic countries. The solution also includes estimates for some northern Europe and Arctic stations that have been involved in the backbone solution, but those have been analyzed only based on the site log information. The stations in Iceland suffer from geodynamical issues, e.g., plate tectonics, earthquakes, volcanoes, and the velocity estimates are more uncertain there. Moreover, the secular velocities may not be sufficient for accurate reference frame maintenance there.

The most significant improvements in our solution with respect to previous estimate (Kierulf et al., 2014) are the better reliability and the improved coverage, especially in all Baltic countries. For part of the stations, the agreement is on the 0.1, 0.2 mm/year level (horizontal and vertical), but for many stations, we achieved significantly smaller uncertainties of the velocity estimates. Our solution benefits from longer time series, but the thorough pre-analysis of the station position time series has been the main key to high quality. In particular, acknowledging the biasing effect of gradual tree growth on the trends, heavy removal of data has been carried out in some cases. Most of our analysis centers have also contributed to the EPN’s dense velocity field solution (Kenyeres 2018), but the same level of homogeneity and accuracy cannot be achieved without local knowledge of the stations, as demonstrated in the comparison to some of the ITRF2014 velocities.

The spectral analysis showed a more robust estimate when using the flicker noise model instead of the general power-law noise model. Klos et al. (2017) linked very small spectral indices (smaller than − 0.5) to the short time series in their analysis of the vertical velocity uncertainties, using a global set of IGS reprocessed time series (Rebischung et al. 2016). We found those also for some of the longest time series in southern Sweden. Klos et al. also found spectral indices close to flicker noise at stations near the Baltic Sea, but we could not observe any similarity there with a much denser network. This points out the different stochastical properties of global and regional solutions. We should have long enough time series from the oldest stations for detecting potential random walk (Klos et al. 2017), but it had a negligible role in our results. However, we recognize the difficulty in estimating small parts of RW in the presence of FN noise pointed out by Dmitrieva et al. (2015). The removal of the common mode error in the network may have revealed better RW characteristics, but we did not expect to get significant differences in the final uncertainty estimates in that way.

The estimated uncertainties for the Finnish and Swedish stations, having almost 20 years of data, were lowest on the 0.01, 0.03 mm/year level (horizontal and vertical), raising some speculations about still too optimistic uncertainties. We removed only annual and semi-annual signals in the uncertainty analysis, and the potentially remaining periodical signals, like GPS draconitic year, should evenly increase the uncertainty. However, the difference between the 3- and 10-degree solutions was less or equal to the 3-sigma noise estimate, showing very robust velocity estimates. In practice, the large amount of data drastically decreases the uncertainty despite the statistically correct noise modeling.

Our uncertainty estimates describe only the noise component of the time series. The horizontal velocity differences found at Swedish co-located stations suggest an unknown uncertainty that is not observable without co-located stations. A potential reason for systematic differences could have been the differences in the time span combined with the changes in the satellite system (satellite types and constellation) over time. However, in our case study in the pre-analysis phase, we observed a similar size of differences though the same time span was used. Thus, the reason for the systematic differences remains still unknown and needs to be re-analyzed in the next reprocessing.

Conclusions

We have produced an ITRF2014 position and velocity solution including uncertainties for 252 stations mainly in the Nordic and Baltic countries. We combined up to 20 years of daily solutions with full covariance matrices instead of a station-wise and component-wise analysis. The velocity field is more accurate and denser than any previous estimates through thorough and individual pre-analysis of the station position time series. We found significant velocity differences at the co-located stations, pointing out biases not included in the formal uncertainties. We concluded the flicker plus white noise uncertainty estimates were more robust compared to the general power-law estimates, though they may be still somewhat too optimistic for single stations.

This ITRF2014 densification solution is a long-awaited update of the Nordic–Baltic velocity solution and now includes a position solution as well. High-quality station velocities and positions are crucial for maintaining the Nordic–Baltic reference frames that are affected by the GIA. In the future, we aim to produce regular cumulative updates on this solution to take full benefit of the newly established stations in many countries and to consider, for example, equipment changes at the stations.