New 3D velocity model of Estonia from GNSS measurements

. The aim of this study was to create a 3D crustal deformation model for Estonia, based on dense Global Navigation Satellite System (GNSS) data (geodetic points with velocities) and validate the existing models of horizontal and vertical crustal deformations with velocities from Estonian GNSS measurements. The observations performed for at least eight years at Estonian GNSS permanent stations and during the GNSS campaign measurements of 1997, 2008 and 2017 on the Estonian 1st-order geodetic reference network were used as input data. Coordinates of the geodetic points were calculated in the ITRF2008 reference frame using the Precise Point Positioning method. Horizontal and vertical velocities (in the North, East and Up directions) were calculated for a total of 22 GNSS points. Models for horizontal and vertical velocities were calculated using the remove–compute–restore method. The model of glacial isostatic adjustment (GIA) of the Nordic Geodetic Commission NKG2016GIA was used as a reference model. Residual velocities of GNSS points showed a good fit with respect to the reference model. The residual velocities were analysed by geostatistical methods and the prediction surfaces of the residual velocities were modelled. After adding the surface of the residual velocities back to the reference model NKG2016GIA, the modelled surface EST2020VEL was obtained. The obtained model was compared with the up-to-date intraplate deformation model NKG_RF17VEL. It was found that recent Fennoscandian intraplate deformation models NKG2016LU and NKG_RF17VEL fitted well with the Estonian GNSS data. However, both models are systematically shifted with respect to the Estonian GNSS data. For applications in Estonia, it is better to use the fitted model EST2020VEL. The uncertainty of the model is estimated to be lower than ±0.5 mm/a.


INTRODUCTION
In terrestrial reference systems such as the Earth-centred Earth-fixed (e.g.ITRF2014) or Eurasian plate-fixed (e.g.ETRF2014) systems, the measured crustal deformations on the Estonian territory are primarily caused by two processes: (i) plate tectonics, which predominantly cause horizontal movements and (ii) postglacial rebound, which primarily causes vertical deformations but also, to a lesser extent, horizontal deformations.Postglacial land uplift has been studied in Estonia for decades, mostly using repeated levelling data (Zhelnin 1966;Vallner et al. 1988;Kall et al. 2014), sea-level observations (Jevrejeva et al. 2002;Suursaar & Kall 2018), as well as Global Navigation Satellite System (GNSS) data obtained from continuously operating GNSS stations (CORS) (Kall et al. 2019).Until now, the determination of horizontal deformations was complicated, primarily due to the scarcity and low precision of the repeated measurements (mostly triangulation measurements) (Rüdja 2004).Thanks to the densification of the CORS network, studies on horizontal deformations have also gained increasing attention in Estonia.Several studies have already used Estonian CORS data to create horizontal and vertical crustal deformation models (Lidberg et al. 2007;Kierulf et al. 2014;Häkli et al. 2019;Lahtinen et al. 2019;Vestøl et al. 2019).Models of vertical crustal deformations have also been used in geology, environmental sciences, oceanography, etc. (Hulisz et al. 2016;Rosentau et al. 2017;Vilumaa et al. 2017).Repeated GNSS campaign measurements on the 1st-order points of the Estonian Geodetic Reference Network (GRN) in 1997, 2008 and 2017 have been accrued as new valuable data for the determination of the horizontal and vertical crustal deformations in Estonia (Metsar et al. 2018(Metsar et al. , 2019)).
Recently, 3D velocity models came to light in connection with the decision of several countries to change their geodetic reference frame from static to semidynamic or dynamic (Bitharis et al. 2017;Poutanen & Häkli 2018;Kierulf et al. 2019).The motivation for the current study was to validate the existing models of horizontal and vertical crustal deformations with velocities from precise GNSS measurements.The objective was to find out how the velocities based on episodic GNSS measurements on GRN points fit into the velocities based on CORS time series and to create a 3D crustal deformation model for Estonia.Accordingly, consistent 3D coordinates of selected CORS and GRN points were calculated by the Precise Point Positioning (PPP) method using the software package Gipsy-Oasis II v. 6.4 (GOA) (Zumberge et al. 1997) in the ITRF2008 reference frame (Altamimi et al. 2011); North, East and Up (NEU) topocentric coordinate time series were formed.From time series analysis, velocities of the points in ITRF2008 were found.Eurasian plate velocities were removed from the N and E velocities using the ITRF2008-PMM model (Altamimi et al. 2012) to represent intraplate velocities, e.g.postglacial rebound velocities in Estonia.For the 3D intraplate deformation model, the remove-compute-restore (RCR) method (Sjöberg 2005) was used.In the first stage, residuals from the glacial isostatic adjustment (GIA) model were calculated.The GIA model NKG2016GIA (Vestøl et al. 2019) was used as the reference model.In the second stage, horizontal and vertical components (2D + 1D) of the model of the residuals were calculated separately using a least squares fit of the planar surface.In the third stage, velocities of the NKG2016GIA model were added back to the models obtained in the second stage.The obtained final deformation model was compared with a horizontal and vertical intraplate deformation model for the Nordic and Baltic countries NKG_RF17VEL (Häkli et al. 2019).
This paper is organized as follows.In the 'Materials and methods' section, an overview of the modern Estonian geodetic reference network (passive GRN and active CORS) is given.In addition, a methodology of coordinate calculation and time series analysis based on episodic GNSS measurements on GRN points and CORS is presented.In the 'Results' section, the calculated velocities based on time series of GRN points and CORS are presented.Different methods for the calculation of a 3D velocity model are discussed.The results of the comparison of the calculated model, called EST2020VEL, with a similar model for the Nordic and Baltic Countries, NKG_RF17VEL, are introduced.Finally, the most important findings are discussed and summed up in 'Discussions and conclusions'.

GNSS campaign measurements on 1st-order points of the national geodetic reference network
The modern geodetic reference network of Estonia was established in 1996-1997 in order to more precisely utilize the new European Terrestrial Reference System 1989 (ETRS89) in the Estonian territory.The first attempt to realize ETRS89 in Estonia was made in 1992 during the EUREF-BAL92 GPS measurement campaign, where five of 43 1st-order network points were connected to the EUREF network under the guidance of the Nordic Geodetic Commission (NKG).Unfortunately, after the EUREF-BAL92 verification campaign in 1994 it was discovered that this realization had systematically shifted around 45 mm in the E direction, primarily due to the shortcomings in the computational strategy.Therefore, it was decided at the beginning of 1996 to replace the old 1st-order reference network with a new geodetic reference network divided into 1st, 2nd and 3rd orders (Rüdja 2004).To connect the new network with the old coordinate systems used in Estonia, as many old network points as possible were integrated to the 3rd order.For the 1st-and 2nd-order networks, brand new points, 212 in total, were installed.
The new 1st-order network consists of 13 points, including the first CORS Suurupi (SUR4), which started operation in 1996.The average distance between 1st-order points is ~100 km (Fig. 1).The first GNSS campaign of the new 1st-order network, known as EUREF-Estonia-1997, was carried out on 19-28 July 1997.Ten GPS receivers Ashtech Z-12 with Choke Ring antennas were used.Due to the low number of GPS receivers, measurements were performed using a complicated scheme: some receivers were stationary during the whole campaign and the rest moved from point to point, according to the session table.The measurement sessions ranged from 60 to 204 h (Table 1).To calculate the coordinates of the measured points, eight surrounding International GNSS Service (IGS) CORS stations in the neighbouring countries were selected for fiducial points.Coordinates were calculated in the ITRF96 reference frame in the central epoch of the measurements 1997.56 and were transformed into the ETRS89 system (reference frame ETRF96).These coordinates of the 1st-order points of the GRN belong to EUREF B-class (precision of the coordinates ±1 cm in the measurement epoch) according to the resolution of the EUREF symposium in Prague on 2-4 June 1999 (Rüdja 1999).
The GNSS measurements on the 1st-order points of the GRN were repeated in 2008, in order to (i) determine precise coordinates for the Estonian network of CORS, called ESTPOS, in the national geodetic system and (ii) detect spatio-temporal deformations of the GRN.A detailed overview of the ESTPOS network is given in Metsar et al. (2018).The GNSS campaign was carried out from 28 July to 08 August 2008August (mean epoch 2008.59.59) and the same antennas and receivers were used as in 1997 (12 Ashtech Z-12 GPS receivers and GPS Choke Ring antennas).The increased number of receivers and antennas allowed simultaneous measurements of all 12 points during a single measuring session.The length of the sessions ranged from 103 to 161 h (Table 1) depending on the location of the receiver.In addition to 12 GRN points, Suurupi (SUR4), Kuressaare (KURE), Tõravere (TOR2), Toila (TOIL) and Audru (AUDR) CORS were included in the campaign (Fig. 1) (Ellmann et al. 2008).
The third GNSS campaign was carried out on 7-15 August 2017 (mean epoch 2017.61);lengths of the sessions varied from 73 to 168 h (Table 1).Measurements were divided into two sessions: the duration of the first session (7-11 August 2017) was 84 h and of the second session (11-15 August 2017), 91 h.Leica antennas LEIAT504GG (Choke Ring), LEIAR25.R4 (3D Choke Ring) and Leica GNSS receivers GRX1200 and GR25 were used.The GPS and GLONASS observables were recorded with a 1-s interval (Kollo et al. 2017).

Coordinate and velocity calculation of the 1st-order geodetic reference network
The velocities of the GRN from repeated GNSS campaigns found by Kruusla (2019) were used as a first input for the calculation of the 3D velocity model.A brief description of the calculation method used is presented below.During the reprocessing of the GNSS campaigns of 1997, 2008 and 2017, Earth-centred Earth-fixed (ECEF) coordinates of the 12 GRN 1st-order points (Fig. 1) were found in the ITRF2008 reference frame using the PPP method and software program GOA.The Suurupi CORS (SUR4 with ID 6392) was excluded from calculations, since its 3D velocity can be more reliably calculated using long time series analysis (see the section 'Network of continuously operating GNSS stations ESTPOS and its velocities').
Due to the specifics of GOA PPP data processing, where GNSS observations at one point are processed by a maximum of 24-h-long measurement sessions, shorter sessions which were performed at the same point, on the same day, with the same receiver, antenna and antenna height, were merged into one 24-h RINEX file using the software package TECQ.Sessions shorter than 4 h were not used for coordinate calculations.Measurement sessions of 2008 were 12 h long, and only shorter sessions were merged.Measurements of 1997 and 2008 were recorded with a 30-s interval and measurements of 2017 with a 1-s interval.In the interest of measurement unification, data from 2017 were also decimated to 30-s observations.As a result of this pre-processing, RINEX observation files of all 12 GRN points on three GNSS campaigns were obtained.All measurement sessions were 12-24 h long with a 30-s interval.Although it is advantageous to use 24-h sessions, as most of the possible sub-daily periodic unmodelled signals could be smoothed out, studies have shown that similar repeatability of the coordinates can be achieved with over 3-h-long sessions (Gandolfi et al. 2017).The default GOA settings were used to estimate troposphere parameters (zenith delay, troposphere gradients).The loading effects due to the atmospheric and hydrological mass variations were not modelled.The IGS (Johnston et al. 2017) final orbit and clock data products in ITRF2008 (IGb08) calculated by the Jet Propulsion Laboratory (JPL) analysis centre were used.The main processing parameters are presented in Table 2.
Initially, several coordinates (from every measurement session) were obtained for one point from PPP processing.They were merged into a single weighted average coordinate with full covariance matrix of a given point of a specific GNSS campaign using GOA's utility stamrg and covariance matrix of each merged solution.Before merging the coordinates, differences of coordinates of one point from separate sessions were studied.Repeatability of the NEU coordinates was calculated for each point (Table 3).As can be seen, average repeatability for the 1997, 2008 and 2017 campaigns was relatively similar.The smallest repeatability was obtained for the 2017 campaign.
Next, time series of NEU coordinates from different measurement sessions of one GNSS campaign for the GRN points were formed (e.g.Fig. 2).Time series were inspected and outliers were removed using software TSAnalyzer (Wu et al. 2017).In the final step, time series of weighted average coordinates of different GNSS campaigns were created for every GRN point (Fig. 3).On the basis of these time series, velocities of the 1st-order points of the GRN were found.
The velocities of the 1st-order GRN points were estimated from the weighted average coordinates of the three GNSS campaigns (central epochs 1997.56; 2008.59; 2017.61) using the software package TSAnalyzer.The program uses the linear weighted least squares method, where weights were found based on the standard deviations of coordinates.In addition, velocities based on the two   , 1997.56 and 2008.59, 2008.59 and 2017.61, and 1997.56 and 2017.61 were found, in order to detect change in velocities over time.The significance of differences between velocities was analysed using a z-test.
The precision and reliability of high-rate GNSS data are heavily influenced by systematic errors.The main systematic errors include the common mode error (CME) caused by the spatial correlations between the GNSS stations (Li et al. 2019).Several methods (e.g.regional stacking filtering (Wdowinski et al. 1997;Teferle et al. 2008), principal component analysis (Jackson & Chen 2004) and Karhunen-Loeve expansion (Dong et al. 2006)) are proposed to reduce the effect of the CME.However, there is no general CME filtering method and since the physical origin and spatial distribution of the CME are unclear or undetermined, it is difficult to extract the correct CME accurately and reliably.For example, Tarayoun et al. (2018) observed increased root mean square (RMS) values of daily positions of campaign stations after CME correction compared to CORS.Further studies are required on this topic (He et al. 2017;Montillet & Bos 2019).As the effect of removing the CME from campaign data remains controversial, no CME has been filtered out from the GRN time series.

Network of continuously operating GNSS stations ESTPOS and its velocities
A detailed overview of the calculations of coordinates and velocities of CORS is given in Kall et al. (2019).To summarize, the ECEF coordinates in ITRF2008 were calculated for ten CORS stations (Fig. 1), which had an operating period of at least eight years and contained the least interruptions in time series due to the antenna or site change.Eight stations belonging to the ESTPOS were equipped with the Leica GNSS receivers GR25 and AR25 GNSS Choke Ring antennas with radome LEIS.Two stations (KARG and MISS) belong to the private network Trimble VRS Now and are equipped with the Trimble NetR5 and NetR9 receivers together with the TRM55971.00antennas.Daily coordinates for the period of 2008-2016 were calculated using the software program GOA with the same processing parameters which were used for coordinate calculation of the GRN points (Table 2).The software package Hector 1.6 (Bos et al. 2013) was used for the coordinate time series analysis and velocity calculations.An example of the time series for one station is presented in Fig. 4. For realistic velocity uncertainty estimations, six different models for accounting for temporal correlated noise were tested.Final velocities and their uncertainties (Table 4) were calculated using the Flicker noise + White noise model for the N and E components and a Generalized Gauss-Markov model for the U component.Decisions were made based on the Akaike and Bayesian information criteria estimated by Hector.

Velocities of the GRN points and CORS
The velocities of the GRN points and their standard deviations are given in Table 5.The intraplate velocities, which were obtained after velocities of the points due to plate tectonics (ν ITRF 2008-PMM ) according to the model ITRF2008-PMM were subtracted from the observed velocities (ν obs ).These are presented in Fig. 5.In addition, Table 5 also includes residual velocities (ν res ), which were obtained after the velocities of the points due to the GIA (ν GIA ) according to the model NKG2016GIA (Vestøl et al. 2019) were subtracted from the intraplate velocities (ν obs -ν ITRF 2008-PMM ), according to the formula (1) Similar to the GRN points (Table 5), the velocities of the CORS and their standard deviations are given in Table 4.The intraplate velocities of CORS (ν ITRF 2008-PMM ) are presented in Fig. 6.Table 4 also includes residual velocities (ν res ), ac cording to Eq. (1).
Tables 4 and 5 indicate that the velocities of the CORS and GRN points present relatively similar features related to the GIA model.The largest shift is in the N direction, the least in the E direction.Moreover, standard deviations of the residual velocities are similar in the N and E directions (for CORS ±0.31 and ±0.30, respectively).However, the standard deviation of residual velocities in the U direction is twice as large (±0.65 for CORS), also indicating larger uncertainty of the reference GIA model in the U direction.

Change in the velocities of the GRN over time
The statistical significance of the velocity differences of GRN points were tested using the z-test in MS Excel.The z-test is used to determine whether two population means are different when the variances are known.Zero hypothesis was that the difference of velocities between the periods of 1997-2008 and 2008-2017 is zero.
The velocities ν for point i in N, E and U directions were calculated using the formula (2) ( , , ) ( , , ) and the variances of velocities (σ ν i 2 ) for point i needed for the z-test were calculated using the formula (Holdahl 1978) (3) where x 1 and x 2 are the N, E or U coordinates and σ 1 and σ 2 mark the standard deviations of the N, E or U coordinates for epochs 1 and 2 (1997 and 2008 or 2008 and 2017, respectively), and t is the length of the period in years (1997.56-2008.59 = 11.03 or 2008.59-2017.61 = 9.02, respectively).The standard deviations of the velocities were multiplied by a factor of 3.29 to increase the confidence level to 99.9%.The results of the test are presented in Table 6.
Velocity differences between the periods of 1997-2008 and 2008-2017 are mostly statistically insignificant, as Table 4.The velocities (ν) of the CORS in ITRF2008 and standard deviations of the velocities (σ at 68% confidence level) in North (N), East (E) and Up (U) directions.Residual velocities were obtained according to Eq. ( 1) by removing the motion due to plate tectonics (model ITRF2008-PPM) and GIA (model NKG2016GIA).Unit is mm/a Table 5.The velocities (ν) of the 1st-order points of the GRN in ITRF2008 and standard deviations of the velocities (σ at 68% confidence level) in North (N), East (E) and Up (U) directions.Residual velocities were obtained according to Eq. ( 1) by removing the motion due to plate tectonics (model ITRF2008-PPM) and GIA (model NKG2016GIA).Unit is mm/a B (°) L (°) Name of the station can be seen from Table 6.The velocity difference in the E direction is statistically significant for three points (MÄEBE97, LONDI97 and SUURSOO97).It is known that defective tribrach was used in point SUURSOO97 during the 2008 campaign (Ellmann et al. 2008), which could cause such a velocity difference.However, when the N and E directional velocities were recalculated to one horizontal velocity vector, the horizontal velocity differences were also statistically insignificant for those three points.Therefore, it can be concluded that the rates of horizontal and vertical velocities of 1st-order points of the GRN between two separate periods are equal.

Calculation of the 3D deformation model
The RCR method (Sjöberg 2005) was used for the calculation of the deformation model.In brief, to employ the RCR method, GIA velocities based on NKG2016GIA were subtracted from calculated 3D velocities (in the N, E and U directions) of the GRN points and CORS.From the horizontal velocities, the Eurasian plate velocities based on ITRF2008-PPM were subtracted (remove phase, Eq. ( 1)).The obtained residual velocities are presented in Tables 4 and 5.The autocovariance or semi-variance of the residual signal was estimated on the basis of the residual velocities at observation points (dN, dE, dU) and prediction surfaces of dN, dE, dU were modelled (compute phase).By combining the obtained prediction surfaces with the NKG2016GIA, e.g.restoring the GIA signal (restore phase), the RCR solution is obtained.
For studying the spatial variability and correlation of the residual velocities, empirical isotropic variograms of the N, E and U directional residual velocities were calculated using the software program Vesper 1.6.3(Whelan et al. 2002) (Fig. 7).As expected, variation in the residual velocities is much less in the N and E directions compared to U. It is obvious from the models that horizontal movements caused by GIA and plate tectonics are similar in direction and size for a small territory such as Estonia and therefore are well estimable from the models.It is clear from the variograms that the semi-variance in the N direction quickly turns to constant (Fig. 7A), i.e. correlation between residual  Table 6.Velocities (ν) between periods 1997Velocities (ν) between periods -2008Velocities (ν) between periods and 2008Velocities (ν) between periods -2017 in North (N), East (E) and Up (U) directions, standard deviations of velocities (σ at 99.9% confidence level) and probability of significance of velocity differences (p) of 1st-order points of GRN.Probabilities marked with grey background are statistically significant at the level α = 0.001.Unit is mm/a Name of the point velocities is insignificant.This suggests that, for fitting N-directional velocities with the GIA model, it is sufficient to use a constant shift between them.On the other hand, semi-variances of E and U residual velocities rise almost linearly until the distance 5° (Fig. 7B, C), suggesting that there is a trend in the data.Therefore, for fitting with the GIA model, constant shift is not sufficient here, as was the case with the N-directional residual velocities.
Vesper software was also used for the modelling of empirical variograms.The best fitting model was selected based on the least sum of squared error and Akaike Information Criterion.The lowest Akaike Information Criterion pertains to the best model (Akaike 1974).The Gaussian model (Minasny et al. 2005) was identified as the best model for all, i.e.N-, E-and U-directional empirical variograms of residual velocities.Modelled variograms are also presented in Fig. 7.The variogram models were used for spatial interpolation of the residual velocities using the ordinary Kriging method where the number of pairs was used as the weight.As a result, N-, E-and U-directional residual velocities were obtained in grid form with a spacing of 0.1666667° in the E direction and 0.0833333° in the N direction.
For further analysis of spatial correlation between the pointwise data, autocovariance plots of the residual velocities were created (Fig. 8).These plots show that, in the case of N and E, the covariance of the residual velocities is small and quickly decreases to zero with increasing distance between points.This means that the employed GIA model predicts velocities correctly and the remaining residual velocities can be characterized as random noise plus a constant or linear (described with the slope plane) shift.This also means that the modelling of the residual velocities with Kriging or least squares collocation (LSC) is probably not beneficial and it could be sufficient to fit a constant shift or slope plane to the residual velocities.However, the semi-variance of U is increasing linearly and the curve of autocovariance is going almost linearly to negative, which indicates a trend in the data.When the linear trend is removed using a least squares fit, an autocovariance curve similar to the N and

A B C
E directions can be seen, i.e. an almost non-existent spatial correlation between residual velocities is obtained.Consequently, it can be assumed that fitting a slope plane could also be sufficient to model the residual velocities in the U direction.A new computation for U after removing the linear trend shows almost no autocovariance for intervals greater than 0.5°.The empirical autocovariance curves in Fig. 8 were fitted with an isotropic Gauss-Markov 2nd-order (GM2) model C(l) = C 0 (1 + l/α)e -l/α (Kasper 1971) with two parameters: variance C 0 and correlation length X 1/2 , where α = 0.595X 1/2 implemented in the software program GRAVSOFT (Forsberg & Tscherning 2008) for the LSC method.The correlation length l = X 1/2 is the value of the argument for which C(l) has decreased to half of its value: C(l = X 1/2 ) = 1/2C 0 (Moritz 1980).These fitted GM2 curves were further employed for modelling the residual velocity prediction surface by the weighted LSC using the same grid space that was used in Kriging.The weighting of data was based on the standard deviations of the estimated velocities from Tables 4 and 5. Since the GRN time series contains only three epochs of data, their noise parameters cannot be accurately estimated (Kierulf 2017).Therefore, standard deviations of the velocities are too optimistic, since not all error sources are taken into account (Alothman et al. 2016).To make standard deviations of GRN points' velocities more consistent with CORS results, standard deviations of GRN points' velocities in Table 5 were rescaled using the formulas (4) where σ(N, E, U) is a standard deviation of the GRN point velocity in the N, E or U direction, respectively (Table 5), 1.84 is a t-statistic at 68.3% confidence level with degrees of freedom = 1 (three velocities minus two parameters of ( , ) ( , ) ( , ) 0.08 1.84 . ., , linear fit), 0.08 and 0.32 mm/a are the average RMS of residual signals for the N and E, and U components, respectively, estimated from CORS time series when the linear trend was removed.These RMS values describe the periodic (annual, seasonal, etc.) and irregular variations of coordinates due to the loading and other effects.
To determine the best/optimal method (least squares fitting of constant shift or slope plane, or more advanced methods like Kriging, or LSC for the modelling of the residual velocities, all the above-mentioned methods were tested.The test computations revealed that the best result in the fitting of the N component was achieved by applying the constant shift (-0.94 mm/a, Tables 4 and 5).The E and U components, however, need the advanced methods mentioned earlier.The variograms and autocorrelation plots (Figs 7,8) also support those results.
Next, the residual velocities were fitted with the simple slope planar surface (z(x, y) = A + Bx +Cy) using the unweighted least squares method.Goodness of the models was estimated by the fit between the model and data points (statistics of the residuals between the modelled surface and data points: mean, standard deviation, RMS difference).The results of the comparison are presented in Table 7.In the calculation of the statistics of the E-directional differences, the GRN point SUURSOO97 was omitted due to the tribrach problem mentioned above.
Although the LSC model showed the lowest standard deviation of residuals between the modelled and observed velocities, it suffered from a bull's-eye effect (Fig. 9), which is presumably due to the low and diminishing correlation between the points' velocities found from the autocovariance analysis (Fig. 8).The plane fit gives the same descriptive statistics (Table 7) as Kriging, therefore the simpler model (slope plane least squares fitting) was preferred for all residual velocity components (N, E and Table 7. Statistics of the differences between modelled surfaces of residual velocities and residual velocities themselves (Tables 4  and 5 in North (N), East (E) and Up (U) directions.Unit is mm/a Fig. 9. Least squares collocation model for residual velocities (red arrows) of the Up component.Note 'bull's-eyes' around data points (triangles), due to the lack of autocorrelation between points.
-0.12 U) as the equivalent solution to more complex surfaces (Kriging and LSC).A good result for the N component was also obtained by applying the constant shift (-0.94 mm/a) to the GIA model.However, statistics of the differences between the model and data points are worse compared to the Kriging and slope plane models.
The fitted slope models of the residual velocities were added back to the predicted velocities of the gridded model NKG2016GIA (restore phase).The obtained horizontal and vertical intraplate velocities of the Earth's crust in the ITRF2008 (vertical movements are referred to ellipsoid GRS80) are presented in Fig. 10.Horizontal velocities are represented as a 2D vector model, where the direction of the vector coincides with the direction of the intraplate horizontal movement and the length and colour of the vector represent the magnitude of the movement.The model of the intraplate vertical movements is presented as the classical model with colours and isolines, where the isoline and colour tone represent the location of the velocities with the same magnitude.The obtained models (Fig. 10) were named EST2020VEL for ease of reference.
To estimate the uncertainty of EST2020VEL, differences between the velocities interpolated from the model and observed velocities at data points were calculated for N, E and U components separately.In addition, cross validation of EST2020VEL was performed using the leave-one-out method (Hastie et al. 2009).Furthermore, even though the Kriging and LSC residual velocity models were not used for the calculation of the final EST2020VEL, their spatial uncertainty estimates will provide additional information about the uncertainties of the final model.The estimates obtained using the different methods and the uncertainty of the reference model NKG2016GIA are summarized in Table 8.
Table 8 shows that all methods for the estimation of EST2020VEL uncertainty will give a more or less similar result of ~0.3 mm/a (except for the RMS of the LSC standard deviation for the U component, which is somewhat larger than other estimates).The average uncertainty of NKG2016GIA is approximately 0.4-0.5 mm/a in Estonia (Vestøl et al. 2019).As a result, the combined uncertainty of EST2020VEL can be estimated to be about 0.30 mm/a for N and E, and 0.45 mm/a for U.

Comparison of the model EST2020VEL
The modelled horizontal and vertical velocities were compared with the recent model of Fennoscandian intraplate deformations NKG_RF17VEL (Häkli et al. 2019), the vertical velocity component of which is represented by the semi-empirical land uplift model NKG2016LU_abs (Vestøl et al. 2019).The latter is compiled using the GIA model NKG2016GIA and the time series of CORS and repeated levellings (including data from Estonia).The horizontal component of NKG_RF17VEL is also based on the NKG GIA model and time series of CORS.The computation of NKG_RF17VEL was performed using the same methodology (RCR) as for EST2020VEL.However, the residual velocities were modelled using the LSC method and the velocities are in a different reference frame (ETRF2014).The vertical velocities are in ITRF2008, as are the EST2020VEL velocities.For proper comparison, the horizontal velocities (N and E components) of the EST2020VEL grid were transformed to ETRF2014.The descriptive statistics of differences between the models are presented in Table 9.
Table 9 shows that standard deviations of differences of the models are small, which shows a good fit between NKG_RF17VEL and EST2020VEL.The RMS differences also remain within the limits of velocity uncertainties of used GNSS velocities.Note that there is a small bias between the N and U components of the models, characterized by the mean difference, which should be taken into account when looking to use NKG_RF17VEL over the Estonian territory.The N component differences are tilted from NW to SE, being largest in SE Estonia.The E direction differences are largest in the central part of Estonia.The largest differences were found in the U component and these are largest in western Estonia (more specifically, in the islands of Saaremaa and Hiiumaa).

DISCUSSION AND CONCLUSIONS
The aim of this study was to improve the 3D velocities of the existing GIA model over Estonia using dense GNSS velocity data.As a result of this study, the 3D (northsouth, east-west, and vertical component) velocity model of intraplate deformations EST2020VEL was obtained.The result was compared with the up-to-date intraplate deformation model NKG_RF17VEL.One important finding of this study was that recent Fennoscandian intraplate deformation models NKG2016LU and NKG_RF17VEL fitted very well with Estonian GNSS data, which was expected since they are correlated, i.e.Estonian GNSS data were used for these models.However, the previously mentioned models have a small systematic bias with respect to the GNSS data used in this study, which should be taken into account when using these models for the Estonian territory.In applications for Estonia, it is recommended to use the fitted model EST2020VEL.
The GIA model NKG2016GIA was used as the EST2020VEL reference model.The differences between the observed GNSS velocities and NKG2016GIA were calculated and so-called residual velocities were obtained.

A B
Next, the residual velocities were analysed by geostatistical methods and the prediction surfaces of the residual velocities were modelled.The modelled surface EST2020VEL was obtained after restoring the surface of the residual velocities to the reference model NKG2016GIA.The systematic shift between GNSS velocities and the GIA model was found for numerous reasons, the most important of which is that GNSS calculations and GIA models are based on different reference frames.Calculations of the ITRF reference frame's Earth centre include the mass of solid Earth as well as the masses of the surrounding hydrosphere and atmosphere, etc. Calculations of the GIA's reference frame only take the mass of the solid Earth into account.However, the standard deviations of the residual velocities are quite small in the N and E directions, considering the uncertainties of the GIA model and GNSS velocities.On the other hand, the standard deviations of the residual velocities in the U direction are somewhat larger because the input data of the GIA model (viscoelastic Earth rheology model, chronological model of glaciers thickness and extent) are more sensitive in the vertical direction (largest movement amplitudes and velocities).The fact that a simple slope plane was sufficient to create model surfaces for the residual velocities (for the N component, almost equivalent would be to apply a constant shift) suggests a very good fit of the reference GIA model to the empirical GNSS data.
The small systematic difference between EST2020VEL and NKG2016LU and NKG_RF17VEL is probably related to the differences in the source data.Although the reference model NKG2016GIA is the same for all compared models, the empirical data used to create the models differ.The model NKG2016LU was calculated using the velocities of five CORS and re-levelling data from Estonia as well as equivalent data from neighbouring countries such as Latvia, Finland and further afield, which may cause the regional trend of the NKG model.The velocities of the CORS for NKG2016LU have been found on the basis of the time series of different lengths from this study.The methods for calculating the coordinates and time series of CORS also differ.For the calculation of NKG_RF17VEL, the observations from fewer Estonian CORS have been used as empirical data compared to this study.In addition, no Estonian GRN time series have been used.The methods for calculating the coordinates and time series of CORS and the length of time series also differ.However, the systematic differences are not significant and the small standard deviations of the differences indicate a good fit between the models.
The developed EST2020VEL model can be used in many different applications and research projects.Postglacial uplift models have been used in numerous studies to eliminate the land uplift effect from coastal tide gauge observations.The developed model EST2020VEL_Up is also suitable for this purpose.In studies of climate-change-driven coastal vulnerability, knowledge about deformations of coastal areas is needed.The coordinates of Estonian geodetic points are based on the GPS measurements of 1997, thus the relative position of the points (the greater the distance between the points, the greater the relative difference) has changed.The coordinates of GNSS permanent stations must also be regularly updated to ensure consistency between pre -  epoch 1997.56).The model EST2020VEL can be used to update the geodetic reference system as well as to convert the coordinates of CORS from the epoch of their coordinates to the epoch of the Estonian geodetic reference system.As a recent trend, several countries (Iceland, Canada, New Zealand, etc.) have switched from a static to a dynamic or semi-dynamic reference frame.Primarily, the need for dynamic reference frames arises from the widespread utilization of precise real-time GNSS applications.The EST2020VEL model can also be used as the underlying deformation model of the dynamic reference frame of Estonia.Currently, the uncertainty of EST2020VEL is estimated at approximately ±0.30 mm/a for the horizontal component and ±0.45 mm/a for the vertical component.The uncertainty of the model warrants further investigation.
However, it should be taken into account that EST2020VEL primarily reflects the intraplate deformations caused by postglacial rebound.To account for local deformations (mining areas, major cities), the model needs to be further improved by including time series of newer CORS, re-levelling and InSAR data and to extend the time series of CORS used in the current study.

109T.
Fig. 1.The 1st-order national geodetic reference network points (blue triangles) and continuously operating GNSS reference stations used in this study (red diamonds).
(in mm) in terms of the standard deviation of the North (N), East (E) and Up (U) coordinates and degree of freedom (DOF) from different measurement sessions of GOA PPP processing of the GRN points 10

Fig. 2 .
Fig. 2. The repeatability of the North, East and Up coordinates of the GRN point OJAKÜLA97 (Fig. 1) from different sessions of the GNSS campaign in 2017.Error bars denote standard deviations of the coordinates based on the PPP method.

Fig. 3 .
Fig. 3.The coordinate time series of the GRN point OJAKÜLA97 (Fig. 1) from three GNSS campaigns in 1997, 2008 and 2017 and velocities of the point from weighted linear regression in North, East and Up directions.

Fig. 4 .
Fig. 4. The North, East and Up coordinate (blue squares) time series of the CORS station TOIL (Fig. 1) from 2008 to 2016.The red line represents the linear trend plus annual and semiannual periodic signals.

Fig. 7 .
Fig. 7.The empirical (blue) and modelled (red) variograms in North (A), East (B) and Up (C) directions estimated for the residual 3D velocity components N, E, U of the 1st-order GRN points and CORS in Estonia.The number of lags is 15 in all variograms.Numbers near points in variogram (A) represent the number of pairs.

Fig. 8 .
Fig. 8.The empirical and modelled autocovariance estimated for the residual 3D velocity components N, E, U of the 1st-order GRN points and CORS in Estonia.Different empirical curves were estimated by using different interval width dl values from 0.25° to 0.5° and software (gstat, cov_func_comp).The isotropic Gauss-Markov 2nd-order model (GM2) Cov = C(l) as a function of interval l between points with two unknown parameters was fitted with empirical curves; see text for more details.

121T.
Fig.10.Models of (A) horizontal and (B) vertical intraplate deformations with respect to the geocentric ITRF2008 reference frame.2D vectors represent horizontal velocities where the colour of the vector represents the magnitude of the movement.Velocities of the Eurasian plate are removed using the ITRF2008-PMM model(Altamimi et al. 2012).All units are mm/a.

Table 1 .
The observation lengths of GRN points in hours

Table 2 .
Parameters used in coordinate calculation by

Table 8 .
Estonian Journal of Earth Sciences, 2021, 70, 2, 107-125 122 Estimates of the uncertainty of EST2020VEL in North (N), East (E) and Up (U) directions.Unit is mm/a

Table 9 .
Statistics of the differences between North, East and Up components of the models EST2020VEL and NKG_RF17VEL over the Estonian territory.Unit is mm/a Up 0.04 cise real-time kinematic GNSS measurements and the Estonian static reference frame (ETRS89 (ETRF96),