An Analysis of Vertical Crustal Movements along the European Coast from Satellite Altimetry, Tide Gauge, GNSS and Radar Interferometry

: The main aim of the article was to analyse the actual accuracy of determining the vertical movements of the Earth’s crust (VMEC) based on time series made of four measurement techniques: satellite altimetry (SA), tide gauges (TG), ﬁxed GNSS stations and radar interferometry. A relatively new issue is the use of the persistent scatterer InSAR (PSInSAR) time series to determine VMEC. To compare the PSInSAR results with GNSS, an innovative procedure was developed: the workﬂow of determining the value of VMEC velocities in GNSS stations based on InSAR data. In our article, we have compiled 110 interferograms for ascending satellites and 111 interferograms for descending satellites along the European coast for each of the selected 27 GNSS stations, which is over 5000 interferograms. This allowed us to create time series of unprecedented time, very similar to the time resolution of time series from GNSS stations. As a result, we found that the obtained accuracies of the VMEC determined from the PSInSAR are similar to those obtained from the GNSS time series. We have shown that the VMEC around GNSS stations determined by other techniques are not the same.


Introduction
Vertical movements of the Earth's crust are widely studied along the European coast, which provides a reference for comparing the obtained results. The movements of the Earth's crust are significantly related to the maintenance and updating of coordinate systems, human activity, monitoring of flood hazards and changes in the mean sea level [1][2][3][4][5][6]. These movements are manifested by the horizontal and vertical displacement of tectonic plates [7]. Vertical plate tectonics are more difficult to determine, and they have both natural and sometimes anthropogenic causes [8]. Vertical movements of the earth's crust are classified as relative or absolute (geocentric) movements, direct or indirect movements and point or surface movements. They can be determined based on measurements involving various techniques, including geometric levelling, global navigation satellite systems (GNSS), Doppler orbitography and radiopositioning integrated by satellite (DORIS), satellite laser ranging (SLR), synthetic aperture radar (SAR) and very long baseline interferometry (VLBI). Vertical movements determined with the use of different techniques in the same area may not be identical [9][10][11]. They may also differ for the same measuring technique [12]. These differences can be attributed to environmental and anthropogenic factors, the types of applied data [10] and the data compilation method.
Currently, vertical movements of the earth's crust are most often determined as absolute (geocentric) movements directly in permanent GNSS stations [13]. The density of GNSS stations may vary for the different continents. Some stations are located on the coast near tide gauge stations (TG) and colocated with them. Absolute vertical movements of the Earth's crust on the coast can be determined indirectly from satellite altimetry (SA) and tide gauge (TG) data [13,14].
Tide gauges are widely used to detect sea level changes along European coasts. Changes in sea level on the regional and global scale have been monitored with high accuracy based on satellite altimetry (SA) data since 1993 [15]. Satellite altimetry observations provide information about absolute sea level within a geocentric reference frame. In turn, TGs measure sea level relative to a land benchmark.
The difference between SA and TG observations is the geocentric vertical crust motion at the TG site; therefore, both sea level measurements are combined to assess vertical displacement at TG sites. In recent research, a combination of these measurements has been used to estimate vertical crust movements [10,16,17].
Collaborative elaboration of results from different measurement techniques can be done by interpolation or simultaneous data alignment. Both for the joint elaboration of the results and their mutual verification, it is necessary to provide information on the accuracy of the obtained vertical movements of the earth's crust from individual measurement techniques. This fact is often neglected in the literature, which may lead to misinterpretation of the final results [18]. To obtain the accuracy at the assumed level for each of the techniques, it is sought to determine the optimal conditions that should be met by the time series. Therefore, for mareographic data, the series should be from several dozen to several hundred years long with a monthly resolution. Similarly, requirements are placed on the time series composed of altimetric data [19]. The areas were selected so that the mareographic data was at least 50 years old and, in places, even up to several hundred years. In the absence of such extensive series, the stations available with a shorter time interval were adopted, as was the SA data in the longest possible time interval. This allows the influence of the circulation of the moon's nodes around the sun to be taken into account, which many researchers have overlooked.
In this study, vertical crust movement was investigated based on the differences between the time series of daily/monthly sea levels generated with the use of SA and TG data in each selected TG station on the European sea coast. The change trends in mean sea level were determined from TG observations and SA data in the European coastal zone.
Vertical crustal movements in nearby GNSS stations were also determined. GNSS time series should contain data from at least 3 years of the measurement period [20] (optimally it is 5 years) with daily resolution. GNSS data is at least 5 years old, and PPP compiled is to be unaffected by GNSS network alignment.
An InSAR analysis was carried out in GNSS sites based on the data obtained from a SAR-C sensor mounted on Sentinel-1A/B satellites for one frame of ascending orbits and one frame of descending orbits. The persistent scatterer InSAR (PSI) method supported the determination of deformation at coherent points on the ground over a long time interval [21]. Displacement values were estimated by reducing error sources related to temporal and geometrical decorrelation and atmospheric phase delay [22]. The following issues were taken into account to reliably assess the identified persistent scatterer (PS) points:

1.
The locations at which PS and GNSS data are measured do not coincide; therefore, spatial interpolation is required [9]; 2.
Sentinel-1 data are not synchronised spatially, which means that their start and end times differ within each orbit.
Due to the different technical approach and the time-consuming nature of the calculations, few studies indicate the optimal nature of the time series generated from the PSInSAR data together with the actual accuracy assessment. Typically, these are series with low time resolution (one interferogram per year, quarter or month) in a 1-, 2-or 3-year time frame. Recently, the literature has shown a resolution of 290 interferograms over 5 years [18].
Our article uses an average of 110 interferograms for ascending satellites and 111 interferograms for descending satellites along the European coast for a selected 27 GNSS stations, which gives over 5000 interferograms in a 3-year time interval, providing a time resolution of nearly 1 week.

of 24
To compare the PSInSAR results against GNSS, an innovative procedure was developed: the workflow of determining the value of vertical movement velocities in GNSS stations based on InSAR data.

Materials and Methods
In this study, vertical movements of the earth's crust along the European coast were estimated using tide gauge (TG), satellite altimetry (SA), GNSS and satellite interferometric synthetic aperture radar (InSAR) data (from Sentinel mission). Vertical movements were determined in the locations shown in Figure 1.
year time frame. Recently, the literature has shown a resolution of 290 interferograms over 5 years [18].
Our article uses an average of 110 interferograms for ascending satellites and 111 interferograms for descending satellites along the European coast for a selected 27 GNSS stations, which gives over 5000 interferograms in a 3-year time interval, providing a time resolution of nearly 1 week.
To compare the PSInSAR results against GNSS, an innovative procedure was developed: the workflow of determining the value of vertical movement velocities in GNSS stations based on InSAR data.

Materials and Methods
In this study, vertical movements of the earth's crust along the European coast were estimated using tide gauge (TG), satellite altimetry (SA), GNSS and satellite interferometric synthetic aperture radar (InSAR) data (from Sentinel mission). Vertical movements were determined in the locations shown in Figure 1. The obtained or created time series were decomposed. A linear trend and the trend standard error were determined. The linear regression method and Fourier analysis were used. Data were obtained from:  The absolute vertical crustal movement is calculated with the use of the following formula: where v SA -absolute change in sea level determined from altimetric data (referred to an ellipsoid); v h -absolute vertical crustal movement (referred to an ellipsoid) and v TGrelative change in mean sea level determined from TG observations. The standard error in the determination of the absolute vertical movement of the Earth's crust is given by the below formula: where σ v h -the standard error of determination absolute vertical crustal movement; σ v SA -the standard error of determination change in sea level determined from altimetric data and σ v TG -the standard error of determination relative change in mean sea level determined from TG observations.
The calculated vertical movements of the Earth's crust are verified against the values determined in neighbouring points with the use of other techniques, and the values are varied [23,24]. The basic prerequisite is the closest proximity of available TG [25,26]. The relative movements v R between a GNSS station and a TG should be taken into account [27]. GNSS stations located on the coast provide the best neighbourhood. If GNSS stations are not available, vertical movements are determined at the nearest alternative points using SAR data. Deformations and deformation velocities are calculated based on SAR data [9,28,29], and they should closely correspond to the velocities determined based on GNSS station data.
Data from GNSS stations and TG form time series. SA and SAR data require additional processing to form a time series. The time series for each data set differ in time resolution (daily, several days, weekly, monthly), white and colour noise, a number of outliers and discontinuities (gaps and jumps). The quality of velocity calculations depends on the above factors, as well as the applied methods for eliminating or reducing the impact of undesirable factors.
Several approaches are used to detect and determine the jump value (vertical shifts in the data series) [30][31][32][33][34][35]. Jumping is caused by technical issues, human errors and environmental factors [33,36]. The determination of the vertical interval, which is regarded as a jump, poses the main technical problem. The vertical interval is difficult to determine because the ranks are affected by environmental conditions, as well as outliers and data gaps.
Several solutions for eliminating outliers have been proposed, and most of them are based on the assumptions presented in the literature [37,38]. Outliers are also eliminated by filling the gaps in the time series. Interpolation methods are used for this purpose [39]. Statistical and spectral methods are most commonly used to detect white and coloured noise [40,41]. In most cases, only the main annual and semiannual periods of seasonality are used, but considerably longer or shorter periods can also be applied due to the influence of other factors [27,[42][43][44][45].

Analysis of Vertical Crustal Movement Velocities Based on SA and TG Data
Vertical crust movements in TG sites were estimated along the European coast based on different time series for SA minus TG (Equation (1)). The time series for SA and TG data describing changes in sea level were characterised by nearly identical behaviour. All the time series had seasonal components (annual, semiannual and 18.61-year cycles related to the relative movements of the Moon) and trend, which can be expressed as follows: f F (t) = a + bt + A a cos (ω a t − ϕ a ) + A sa cos (ω sa t − ϕ sa ) + A 18.61 cos (ω 18.61 t − ϕ 18.61 ) where f F (t) is a Fourier function; a is the bias; b is the trend; t is time; A a and A sa are the annual and semiannual amplitudes, respectively; ϕ a and ϕ sa are the annual and semian-  [47]. The least squares method was used to fit the time series of sea level variations for every station, and annual and semiannual amplitudes and the long-term trend in seasonal sea level variations were estimated.
TG measurements constitute long-time series of mean sea level data; they are obtained from the Permanent Service for Mean Sea Level (PSMSL) [48] and provide primary evidence for the rise in the globally averaged sea level. For this study, 27 TG stations on the European coast ( Figure 1 the time series had seasonal components (annual, semiannual and 18.61-year cycles related to the relative movements of the Moon) and trend, which can be expressed as follows: where fF(t) is a Fourier function; a is the bias; b is the trend; t is time; Aa and Asa are the annual and semiannual amplitudes, respectively; φa and φsa are the annual and semiannual phase, respectively; ωa and ωsa are the annual and semiannual angular frequency, respectively [46]; A18.61 is the 18.61-year amplitude; φ18.61 is the 18.61-year phase and ω18.61 is the 18.61-year angular frequency. The 18.61-year cycle is a lunar nodal cycle caused by the relative movements of the Moon. This important precession of the Moon, namely the 18.61-year lunar nodal cycle, causes tidal modulations over a range of interannual time scales. These modulations affect the interpretation of TG data spanning several years, particularly when dealing with water level extremes [47]. The least squares method was used to fit the time series of sea level variations for every station, and annual and semiannual amplitudes and the longterm trend in seasonal sea level variations were estimated.
TG measurements constitute long-time series of mean sea level data; they are obtained from the Permanent Service for Mean Sea Level (PSMSL) [48] and provide primary evidence for the rise in the globally averaged sea level. For this study, 27 TG stations on the European coast ( Figure 1   Gridded daily sea level anomalies with a resolution of 0.25 × 0.25 degrees and the time series for January 1993 to December 2017 from the Copernicus Marine and Environment Monitoring Service (CMEMS) were used. The time series of altimetric data (one series for each altimetric observation point) were obtained. This set combines data from several altimetry missions. Altimetric measurements were corrected for atmospheric effects (ionospheric delay and dry/wet tropospheric effects) and geophysical processes (solid, ocean, pole tides, loading effect of ocean tides, sea state bias and the inverted barometer response of the ocean). Detailed information on the introduced corrections can be found in AVISO [50] and CMEMS [51]. The differences in SA and TG data sets were investigated by correlation analysis. An SA grid point that was highly correlated with TG data was selected. It is simultaneously the closest SA grid to the TG data. The daily SA data were averaged to correspond with monthly TG data. The mean value of the correlation coefficient was 0.93.
The two main goals of the time series analysis are to describe the character of the analysed phenomenon based on a series of observations and to forecast future values. These goals can be achieved if time series elements are identified and described. The elements of a time series include a systematic component which can be described as a trend (linear or not) and a seasonal component where the duration of seasonal fluctuations can vary. Fluctuations can be regarded as a seasonality when their duration does not exceed 1 year, but when the corresponding period is longer, an economic cycle and a random component (noise) appear. These factors have to be identified in a formalised forecasting method. Periodic phenomena are identified in a harmonic analysis where a priori assumptions are not made. The harmonic analysis aims to decompose a time series with the use of cyclic factors on sine and cosine functions related to a given wavelength.
Harmonic analyses are performed to determine the average value of the studied phenomena. A trend was identified in the analysed data series; therefore, oscillations were determined with the use of the following model [46]: where f H (t) is a harmonic function; a is the bias; b is the trend; t is time and n is the number of observations. The amplitudes and phase shift were calculated for different harmonics. The annual, semiannual and 18.61-year harmonics were calculated. The annual, semiannual and 18.61year harmonic functions for Sassnitz stations based on TG and SA time series are presented in Figure 3.    The correlation coefficient for annual variations from two independent observation techniques was determined at 0.92. The mean annual amplitude was ±5.59 cm with an estimation error of ±1.32 cm for TG data and ±6.17 cm for SA data with and estimation error of ±1.44 cm. The mean semiannual amplitude was ±2.57 cm with an estimation error of ±0.63 cm for TG data and ±1.79 cm for SA data with estimation error of ±0.42 cm. The correlation coefficient for semiannual variations was 0.72, which is a satisfactory result, but it was determined at only 0.34 for the 18.61-year cycle. The correlation coefficient for annual variations from two independent observation techniques was determined at 0.92. The mean annual amplitude was ±5.59 cm with an estimation error of ±1.32 cm for TG data and ±6.17 cm for SA data with and estimation error of ±1.44 cm. The mean semiannual amplitude was ±2.57 cm with an estimation error of ±0.63 cm for TG data and ±1.79 cm for SA data with estimation error of ±0.42 cm. The correlation coefficient for semiannual variations was 0.72, which is a satisfactory result, but it was determined at only 0.34 for the 18.61-year cycle.

Analysis of Vertical Crustal Movement Velocities Based on GNSS Data
Vertical crustal movements can be classified as relative v r or absolute v h . Relative vertical crustal movements refer to any point (or surface) that is constant over time. Absolute vertical crustal movements are referred to an ellipsoid [52], and they comprise relative vertical crustal movements, average changes in sea level, eustatic movements and geoid changes over time. Absolute crustal movements are geocentric crustal movements [2] that are determined based on direct GNSS measurements as well as combined TG and SA measurements [13]. The location of a GNSS station relative to the TG station and the location of the point where changes in sea level were measured based on altimetric data have to be determined to compare the results generated by different methods.
The GNSS is commonly used for geodetic measurements because it monitors land movements with high precision [53,54]. This study analysed the measurements from a total of 27 GNSS stations that are colocated with or positioned in the proximity of a TG station on the European coast. The time series of vertical coordinates of 27 GNSS stations on the European coast ( Figure 1) were obtained from the Nevada Geodetic Laboratory (NGL) [55] and SONEL (www.sonel.org, accessed on 12 December 2019). Time series from 3 selected GNSS stations from a total of 27 stations are shown in Figure 5.
The time series covered a period of 4 to 17 years with daily time resolution. The distance from the nearest TG ranged from 0.002 km to 16.069 km ( Figure 6).
The GNSS is commonly used for geodetic measurements because it monitors land movements with high precision [53,54]. This study analysed the measurements from a total of 27 GNSS stations that are colocated with or positioned in the proximity of a TG station on the European coast. The time series of vertical coordinates of 27 GNSS stations on the European coast ( Figure 1) were obtained from the Nevada Geodetic Laboratory (NGL) [55] and SONEL (www.sonel.org, accessed on 12 December 2019). Time series from 3 selected GNSS stations from a total of 27 stations are shown in Figure 5.  The data from the GNSS station form a daily time series covering several to more than 10 years. All the time series feature data discontinuities, jumps and noise. According to [13], the uncertainty of velocity calculations based on GNSS data ranges from 0.1 mm/year to 1.3 mm/year with a median of 0.21 mm/year. The time series of the coordinates from the GNSS station were decomposed. Missing values were filled by interpolation. Vertical shifts were defined with an algorithm developed for the needs of Figure 6. Distance between GNSS stations and the nearest TG (the longest distance between a TG station and a GNSS station was 16.069 km, and the shortest distance between a TG station and a GNSS station was 0.002 km).
The data from the GNSS station form a daily time series covering several to more than 10 years. All the time series feature data discontinuities, jumps and noise. According to [13], the uncertainty of velocity calculations based on GNSS data ranges from 0.1 mm/year to Remote Sens. 2021, 13, 2173 9 of 24 1.3 mm/year with a median of 0.21 mm/year. The time series of the coordinates from the GNSS station were decomposed. Missing values were filled by interpolation. Vertical shifts were defined with an algorithm developed for the needs of this study [32]. The VSED algorithm filters the outlier series with the Grubbs method. It defines and eliminates jumps and fits a trend line. The operations performed by the algorithm are presented in Figure 7. Figure 6. Distance between GNSS stations and the nearest TG (the longest distance between a TG station and a GNSS station was 16.069 km, and the shortest distance between a TG station and a GNSS station was 0.002 km).
The data from the GNSS station form a daily time series covering several to more than 10 years. All the time series feature data discontinuities, jumps and noise. According to [13], the uncertainty of velocity calculations based on GNSS data ranges from 0.1 mm/year to 1.3 mm/year with a median of 0.21 mm/year. The time series of the coordinates from the GNSS station were decomposed. Missing values were filled by interpolation. Vertical shifts were defined with an algorithm developed for the needs of this study [32]. The VSED algorithm filters the outlier series with the Grubbs method. It defines and eliminates jumps and fits a trend line. The operations performed by the algorithm are presented in Figure 7.  The mathematical model used in this algorithm is a straight line with steps in different epochs: where t-epoch, h(t)-height difference in epoch t; v-velocity between stations; h 0 -height difference at epoch 0; c 1 , c 1 , . . . , c m -elements of the C matrix and s 1 , s 2 , . . . , s mmagnitude of "jumps". The interval defining the "jump" was based on the standard deviation, which ranged from ±0.2 mm (values centred around the mean) to ±0.7 mm (values scattered around the mean). The number of jumps for the analysed GNSS stations, depending on the received height determination error, is shown in Figure 8.
The number of jumps does not affect the accuracy of height determination for individual epochs in the time series of coordinates from the GNSS station.
The interval defining the "jump" was based on the standard deviation, which ranged from ±0.2 mm (values centred around the mean) to ±0.7 mm (values scattered around the mean). The number of jumps for the analysed GNSS stations, depending on the received height determination error, is shown in Figure 8. The number of jumps does not affect the accuracy of height determination for individual epochs in the time series of coordinates from the GNSS station.

Analysis of the Vertical Movement Model Based on InSAR Data
The InSAR technique relies on interferometric comparison of SAR phase images to determine relative surface motions from the millimetre to centimetre [56][57][58]. In contrast to the pointwise information provided by GNSS, InSAR can provide a spatially dense image of surface displacements. The InSAR technique based on time series analysis (multitemporal InSAR, MTI) was used to achieve such a high level of accuracy and to reduce the basic error related to temporal and geometrical decorrelation and atmospheric phase delay. The course of deformation in the selected area was reconstructed with the use of point methods based on the selection of pixels that maintain coherence in time [21,59,60].
The interferometric phase in Equation (6) is the sum of contributions from several factors, including the components to be extracted, and it specifies ground deformation in the LOS (line of sight) direction in the time interval of the SAR image pair. The remaining

Analysis of the Vertical Movement Model Based on InSAR Data
The InSAR technique relies on interferometric comparison of SAR phase images to determine relative surface motions from the millimetre to centimetre [56][57][58]. In contrast to the pointwise information provided by GNSS, InSAR can provide a spatially dense image of surface displacements. The InSAR technique based on time series analysis (multitemporal InSAR, MTI) was used to achieve such a high level of accuracy and to reduce the basic error related to temporal and geometrical decorrelation and atmospheric phase delay. The course of deformation in the selected area was reconstructed with the use of point methods based on the selection of pixels that maintain coherence in time [21,59,60].
The interferometric phase in Equation (6) is the sum of contributions from several factors, including the components to be extracted, and it specifies ground deformation in the LOS (line of sight) direction in the time interval of the SAR image pair. The remaining factors relating to topography, atmosphere and orbit should be taken into account and eliminated from the interferometric phase.
where ϕ defo -LOS deformation, ϕ atm -atmospheric delay, ϕ orbit -orbit error, ϕ dopo -DEM error and ϕ noise -noise The analysed GNSS stations are located in built-up areas along the European coast ( Figure 1). The method StaMPS (Stanford Method for Persistent Scatterers), based on the PSI algorithm, was selected, and StaMPS software was used in this analysis. The PSInSAR algorithm detects stable scatterers and extracts phase information for these points. The StaMPS framework is a collection of spatial and temporal filtering routines that can be used to estimate the phase components in Equation (6) by assuming a spectral structure. More details are contained in [60]. The processed series of radar images (for 2015-2017) provided mean values of displacement velocities for persistent scatterers in the LOS direction ( Figure 9). On average, a calculated 110 interferograms for ascending and 111 interferograms for descending orbital passes were determined at each point of the GNSS station.  Different objects make good permanent scatterers (PS), including buildings, lanterns and fragments of various structures (bridges, fences, etc.). In areas without infrastructure, rock outcrops, hard unvegetated earth surfaces and boulders can be used as PS points. However, the locations of PS points and GNSS stations do not coincide; the distances between selected PS points and the GNSS station are shown in Figure 10. The longest distance (of more than 300 m) was noted with the WLAD station, which is located on port breakwaters. The distance between PS points and ROTG, SCOA, ACOR and TARI stations, which are also located in the direct proximity of hydraulic structures, did not exceed 64 m. The distance between the PS point and the SABL was determined to be 227.34 m because it is the only GNSS station surrounded by a forest with no permanent structures in the vicinity. The scheme in Figure 11 shows the workflow of determining the vertical movement model based on InSAR data. In the scheme, in Step 4, groups of scatterers that exhibit similar behaviour are analysed. Statistical procedures are used to evaluate their spatial relationships and the consistency of estimated parameters in a given environment. During processing, many points are rejected based on different quality criteria [70].
The spatial distribution of ascending and descending PS point targets, in addition to the effect of temporal decorrelation, is highly related to the orientation of slopes and correlates well with the terrain aspect [61,64,71]. The proposed approaches included terrain aspect and interpolation performed in each range of time, both ascending and descending targets to be the same range time. The nearest neighbour vector (NNV) solution was used (according to the procedure http://gmt.soest.hawaii.edu, accessed on 10 January 2020). The PS points are not regularly spaced in space; the distances from their nearest neighbour with opposite geometry are different. Therefore, the search radius was determined based on the slope terrain analysis and the point density.
The locations of PS points and GNSS stations do not coincide; therefore, approaches were adopted according to the diagram shown in Figure 11. In the first approach (Step 5), the PS point was selected on the following assumptions: 1. The PS point is located on the same type of infrastructure or facility as the GNSS station. Sufficient PS points should be available to detect outliers. Next, the height of the PS points should be estimated to check whether the PS comes from the same technical infrastructure object and not from the surface level; 2. Points should be selected from the area with the same slope. The point and slope heights were derived from data based on Shuttle Radar Topography Mission Global 1 arc second (SRTMGL1); 3. PS point targets from ascending or descending tracks are located in the proximity of a GNSS station, ~10 m. The three-dimensional motion was only partly captured in the calculated rate of change in the LOS direction [61]. Therefore, the components of vertical movement had to be determined from SAR observations in the LOS. Radar satellites acquire image data of the same area on ascending (south to north) and descending (north to south) orbital passes. An LOS velocity V los is composed of the 3D velocity components V EW (East-West), V NS (North-South) and V up (Up-Down) . The LOS displacements can be mathematically transformed to V EW and V up components if available at the same location for all analysed tracks and within the same time period. The V NS displacements are small due to the polar orbit of all SAR satellites and the side-looking image geometry [62]. Various approaches to determining 3D surface motion from InSAR data are discussed in publications [63][64][65][66][67].
The PSI deformation value should be predicted at the location of the GNSS stations. Every location is different, and complex deformation phenomena may occur; however, an assessment of each station is possible [68]. In this article, vertical motion in GNSS stations was calculated based on various view geometries and viewing angles. During data postprocessing, ascending and descending PSI measurements were paired to calculate the vertical component for each PS point using Equation (7 where V los -deformation along the LOS, V up -vertical deformation, V h ald -projection of horizontal deformation in descending azimuth look direction, θ-incident angle and ∆αsatellite heading difference between ascending and descending mode [22]. The scheme in Figure 11 shows the workflow of determining the vertical movement model based on InSAR data. In the scheme, in Step 4, groups of scatterers that exhibit similar behaviour are analysed. Statistical procedures are used to evaluate their spatial relationships and the consistency of estimated parameters in a given environment. During processing, many points are rejected based on different quality criteria [70]. The PS points that met the criteria listed above were 25% of the points (6 stations) ( Figure 10). In the second approach (Step 6), the following was adopted:

•
There are no PS points directly (condition 3 from the first attempt); • Points were selected in a surrounding area with a radius of 500 m; • Eliminate outliers; • Displacements behave linearly in time within a radius of 500 m of the GNSS station.
Then, the NNV was calculated, and the spatial interpolation was performed, taking into account the geostatistical properties of displacements in each measurement epoch. The PSI results were spatially interpolated by the Kriging method. PS point targets from ascending and descending tracks were interpolated separately. Based on a review of the literature [9,72], the ordinary Kriging method with a spherical, exponential and Gaussian semivariogram model were applied to the PSI vector data [61]. The spatial distribution of ascending and descending PS point targets, in addition to the effect of temporal decorrelation, is highly related to the orientation of slopes and correlates well with the terrain aspect [61,64,71]. The proposed approaches included terrain aspect and interpolation performed in each range of time, both ascending and descending targets to be the same range time. The nearest neighbour vector (NNV) solution was used (according to the procedure http://gmt.soest.hawaii.edu, accessed on 10 January 2020). The PS points are not regularly spaced in space; the distances from their nearest neighbour with opposite geometry are different. Therefore, the search radius was determined based on the slope terrain analysis and the point density.

Results
The locations of PS points and GNSS stations do not coincide; therefore, approaches were adopted according to the diagram shown in Figure 11. In the first approach (Step 5), the PS point was selected on the following assumptions:

1.
The PS point is located on the same type of infrastructure or facility as the GNSS station. Sufficient PS points should be available to detect outliers. Next, the height of the PS points should be estimated to check whether the PS comes from the same technical infrastructure object and not from the surface level; 2.
Points should be selected from the area with the same slope. The point and slope heights were derived from data based on Shuttle Radar Topography Mission Global 1 arc second (SRTMGL1); 3.
PS point targets from ascending or descending tracks are located in the proximity of a GNSS station,~10 m.
The PS points that met the criteria listed above were 25% of the points (6 stations) ( Figure 10). In the second approach (Step 6), the following was adopted:

•
There are no PS points directly (condition 3 from the first attempt); • Points were selected in a surrounding area with a radius of 500 m; • Eliminate outliers; • Displacements behave linearly in time within a radius of 500 m of the GNSS station.
Then, the NNV was calculated, and the spatial interpolation was performed, taking into account the geostatistical properties of displacements in each measurement epoch. The PSI results were spatially interpolated by the Kriging method. PS point targets from ascending and descending tracks were interpolated separately. Based on a review of the literature [9,72], the ordinary Kriging method with a spherical, exponential and Gaussian semivariogram model were applied to the PSI vector data [61].

Results
Vertical crust movements in TG stations were estimated from different time series of TG data ( Figure 12a) and SA data (Figure 12b) using robust linear regression [73]. Each difference in the time series of SA and coastal TG data was analysed by considering a linear trend and an annual, semiannual and 18.61-year cycle.  The results of the analysis were used to calculate linear trends in vertical crust movement on the TG sites along the European coast based on SA data minus TG data and at GNSS time series (GNSS stations closest to TG stations were considered). The velocity in coastal stations, estimated from the different time series of SA data minus TG data from the beginning of data to 2018, is shown in Figure 13a. The velocity in coastal stations, estimated from the different time series of SA data minus TG data from 1993-2018, is The results of the analysis were used to calculate linear trends in vertical crust movement on the TG sites along the European coast based on SA data minus TG data and at GNSS time series (GNSS stations closest to TG stations were considered). The velocity in coastal stations, estimated from the different time series of SA data minus TG data from the beginning of data to 2018, is shown in Figure 13a. The velocity in coastal stations, estimated from the different time series of SA data minus TG data from 1993-2018, is shown in Figure 13b. The results of the analysis were used to calculate linear trends in vertical crust movement on the TG sites along the European coast based on SA data minus TG data and at GNSS time series (GNSS stations closest to TG stations were considered). The velocity in coastal stations, estimated from the different time series of SA data minus TG data from the beginning of data to 2018, is shown in Figure 13a. The velocity in coastal stations, estimated from the different time series of SA data minus TG data from 1993-2018, is shown in Figure 13b.  In the TG dataset, trend changes in the European coastal zone ranged from −1.75 ± 1.59 mm/year (Tarifa2 station in Spain) to +5.26 ± 0.19 mm/year (Borkum station in Germany).
In the SA dataset, trend changes in the European coastal zone ranged from +0.87 ± 0.08 mm/year (Ibiza station in Spain) to +4.48 ± 0.08 mm/year (Gdansk station in Poland). The mean trend for all stations was determined at +2.92 ± 0.08 mm/year based on SA data in a different period (1993 to 2017). A comparison of both data sets revealed the highest trend values in the northern European coastal zone (Germany and Poland) and the lowest trend values in the southern European coastal zone (Spain).
The differences in the trend values calculated based on TG and SA data were greatest in the Tarifa2 station (+4.24 ± 1.59 mm/year) and smallest in the Roscoff station (+0.44 ± 0.08 mm/year).
Vertical crust movement values were estimated based on multitime SAR data series, and they ranged from −2.24 ± 0.39 mm/year (IJMU) to +2.88 ± 0.30 mm/year (TGBF). The mean minimum standard error was ±0.20 mm/year (GESR), and the maximum standard error was ± 0.49 mm/year (SAS2). The greatest changes in velocity were observed in the coastal zone of Northwestern Europe. The mean velocity values in the stations in France, Belgium, the Netherlands and Germany were negative in the range of −2.24 ± 0.39 mm/year (IJMU) to −0.14 ± 0.46 mm/year (SMTG). The highest positive values were noted in TGBF (Germany) at +2.88 ± 0.30 mm/year and TERS (Netherlands) at +0.27 ± 0.40 mm/year. The stations in the Baltic Sea zone in Northern Europe moved with mean positive velocities: +0.65 ± 0.45 mm/year in SWIN, +0.19 ± 0.48 mm/year in WLAD, +0.21 ± 0.20 mm/year in GESR and +0.07 ± 0.40 mm/year in WARN. Stations on Rügen island in Germany moved with mean negative velocities: +0.66 ± 0.49 mm/year in SAS2 and +0.53 ± 0.33 mm/year in SASS. The stations in Southern Europe (Spain) were characterised by positive velocities in the range of +0.06 ± 0.22 mm/year (CANT) to +0.47 ± 0.39 mm/year (SCOA). Negative velocity was observed only in the ACOR station (−0.35 ± 0.35 mm/year). A map of vertical crust movements in coastal stations estimated from different time series of GNSS and PSInSAR data is presented in Figure 14. and they ranged from −2.24 ± 0.39 mm/year (IJMU) to +2.88 ± 0.30 mm/year (TGBF). The mean minimum standard error was ±0.20 mm/year (GESR), and the maximum standard error was ± 0.49 mm/year (SAS2). The greatest changes in velocity were observed in the coastal zone of Northwestern Europe. The mean velocity values in the stations in France, Belgium, the Netherlands and Germany were negative in the range of −2.24 ± 0.39 mm/year (IJMU) to −0.14 ± 0.46 mm/year (SMTG). The highest positive values were noted in TGBF (Germany) at +2.88 ± 0.30 mm/year and TERS (Netherlands) at +0.27 ± 0.40 mm/year. The stations in the Baltic Sea zone in Northern Europe moved with mean positive velocities: +0.65 ± 0.45 mm/year in SWIN, +0.19 ± 0.48 mm/year in WLAD, +0.21 ± 0.20 mm/year in GESR and +0.07 ± 0.40 mm/year in WARN. Stations on Rügen island in Germany moved with mean negative velocities: +0.66 ± 0.49 mm/year in SAS2 and +0.53 ± 0.33 mm/year in SASS. The stations in Southern Europe (Spain) were characterised by positive velocities in the range of +0.06 ± 0.22 mm/year (CANT) to +0.47 ± 0.39 mm/year (SCOA). Negative velocity was observed only in the ACOR station (−0.35 ± 0.35 mm/year). A map of vertical crust movements in coastal stations estimated from different time series of GNSS and PSInSAR data is presented in Figure 14.  The estimates of vertical crust movements obtained from GNSS stations and the combination of SA data minus TG data were compared in 27 sites along the European coast where both types of measurements were available.
The movement of the Earth's crust determined from the combination of TG and SA data ranged from −2.31 ± 0.21 mm/year (Borkum, Fischerbalje) to + 2.54 ± 0.06 mm/year (Kolobrzeg, Gedser). The mean standard error was ±0.13 mm/year (TG data from the beginning of data to 2018) and ±0.47 mm/year (TG data for 1993-2017). The greatest standard error values were noted in Les Sables D'Olonne (±0.29 mm/year-TG to 2018) and TARIFA 2 (±1.59 mm/year-TG since 1993) ( Figure 15). The greatest differences in vertical movement calculated based on PSI data and GNSS data were noted in TGBF, SCOA, COUD, SAS2 and SASS stations. Differences were also observed in ACOR, IJMU, SABL and TERS stations. The vertical component of interpolated PS points was calculated according to Equation (7). The values obtained in GNSS stations and interpolated PS points are shown in Figure 16.
coast where both types of measurements were available.
The movement of the Earth's crust determined from the combination of TG and SA data ranged from −2.31 ± 0.21 mm/year (Borkum, Fischerbalje) to + 2.54 ± 0.06 mm/year (Kolobrzeg, Gedser). The mean standard error was ±0.13 mm/year (TG data from the beginning of data to 2018) and ±0.47 mm/year (TG data for 1993-2017). The greatest standard error values were noted in Les Sables D'Olonne (±0.29 mm/year-TG to 2018) and TARIFA 2 (±1.59 mm/year-TG since 1993) ( Figure 15). The greatest differences in vertical movement calculated based on PSI data and GNSS data were noted in TGBF, SCOA, COUD, SAS2 and SASS stations. Differences were also observed in ACOR, IJMU, SABL and TERS stations. The vertical component of interpolated PS points was calculated according to Equation (7). The values obtained in GNSS stations and interpolated PS points are shown in Figure 16.

Discussion
The obtained results present vertical deformations in the proximity of selected GNSS stations located on the European coast along with an assessment of the actual accuracy.
Differences were observed in the rate of changes in average sea level determined based on TG data (long-term observations) along the European coast. The greatest changes were observed in tide gauges in the Mediterranean area. The rate of changes in the mean sea level increased to the south (positive values), and a decrease was observed in only two tide gauge stations (Tarifa2 station in Spain and Concarneau station in France).

Discussion
The obtained results present vertical deformations in the proximity of selected GNSS stations located on the European coast along with an assessment of the actual accuracy.
Differences were observed in the rate of changes in average sea level determined based on TG data (long-term observations) along the European coast. The greatest changes were observed in tide gauges in the Mediterranean area. The rate of changes in the mean sea level increased to the south (positive values), and a decrease was observed in only two tide gauge stations (Tarifa2 station in Spain and Concarneau station in France). The mean changes in sea level along the European coast, determined based on satellite altimetry data, increased towards the north. All changes in sea level were positive.
To validate the presented calculations, the velocity determined on selected stations were compared with SONEL measurements based on TG data. Based on TG data and SA data for 1993-2011, 1993-2014 and 1993-2015 [74]. The results are shown in Figure 17. The comparison produced relatively satisfactory results. The greatest difference between the values calculated in this study and those determined in the SONEL network based on the TG time series was noted in the Borkum station. This discrepancy could be attributed to differences in the analysed time periods (1993-2011 in the SONEL network, 1993-2017 in this study).
Differences were observed in the vertical movements of the Earth's crust determined from GNSS data (Figure 14a). The greatest negative vertical movements were observed in ACOR and SCOA stations. The greatest standard errors in velocity determination were noted in SAS2, ILDX and SCOA stations.
In real data time series, an accuracy of 0.1 mm per year can be achieved over a 100year observation period. In time series without seasonal variations, an accuracy of 0.1 mm per year requires an observation period of 49 years [19]. To achieve an accuracy of 0.5 mm/year, a minimum observation period of 3 years is needed [20]. The above assumptions do not apply to AUT1, ILDX and SCOA stations. The calculated standard error values differ from the predicted values, which indicates that the time series is disrupted (the changes in coordinate values proceed in a disordered manner, and their cause is unknown). According to [13], the uncertainty of velocity determinations based on GNSS data ranges from 0.1 mm/year to 1.3 mm/year with a median of 0.21 mm/year. Both higher and lower accuracy values were obtained in this study. The velocity calculated in this study based on GNSS data differs from that calculated by the analytical centres in ULR, NGL and JPL ( Table 1). The greatest differences were noted in IJMU, COUD, SMTG, SCOA, TARI, IBIZ and AUT1. Considerable differences were also observed in standard error values. These discrepancies indicate that the method of calculating and interpreting vertical shifts in time series (jumps) has a significant impact on the final result.
The results obtained during our research of vertical movement velocity on the coasts of Poland, Denmark and Germany was verified on the basis of the article [75]. Differences in values may be connected with the different length of the time series and chosen time The comparison produced relatively satisfactory results. The greatest difference between the values calculated in this study and those determined in the SONEL network based on the TG time series was noted in the Borkum station. This discrepancy could be attributed to differences in the analysed time periods (1993-2011 in the SONEL network, 1993-2017 in this study).
Differences were observed in the vertical movements of the Earth's crust determined from GNSS data (Figure 14a). The greatest negative vertical movements were observed in ACOR and SCOA stations. The greatest standard errors in velocity determination were noted in SAS2, ILDX and SCOA stations.
In real data time series, an accuracy of 0.1 mm per year can be achieved over a 100year observation period. In time series without seasonal variations, an accuracy of 0.1 mm per year requires an observation period of 49 years [19]. To achieve an accuracy of 0.5 mm/year, a minimum observation period of 3 years is needed [20]. The above assumptions do not apply to AUT1, ILDX and SCOA stations. The calculated standard error values differ from the predicted values, which indicates that the time series is disrupted (the changes in coordinate values proceed in a disordered manner, and their cause is unknown). According to [13], the uncertainty of velocity determinations based on GNSS data ranges from 0.1 mm/year to 1.3 mm/year with a median of 0.21 mm/year. Both higher and lower accuracy values were obtained in this study. The velocity calculated in this study based on GNSS data differs from that calculated by the analytical centres in ULR, NGL and JPL ( Table 1). The greatest differences were noted in IJMU, COUD, SMTG, SCOA, TARI, IBIZ and AUT1. Considerable differences were also observed in standard error values. These discrepancies indicate that the method of calculating and interpreting vertical shifts in time series (jumps) has a significant impact on the final result. The results obtained during our research of vertical movement velocity on the coasts of Poland, Denmark and Germany was verified on the basis of the article [75]. Differences in values may be connected with the different length of the time series and chosen time periods.
The trends in mean sea level were determined by [76] based on TG data for the coastal areas in France (Dunkerque, Roscoff, Saint-Malo, Concarneau, Les Sables D'Olonne, Saint Jean De Luz, Sete and other locations that were not analysed in this study). Relative sea level trends and vertical land movements were determined for 1993-2018. The vertical land movements estimated in the cited study are consistent with the presented findings for the Saint-Malo station (−0.63 ± 0.55 mm/year in the cited study; −0.68 ± 0.30 mm/year in this study), Concarneau station (−0.46 ± 0.42 mm/year and −0.46 ± 0.24 mm/year, respectively), Les Sables D'Olonne station (−0.05 ± 0.37 mm/year and −0.08 ± 0.07 mm/year, respectively), Sete station (−0.87 ± 0.43 mm/year and −0.83 ± 0.13 mm/year, respectively) and Roscoff station (−1.28 ± 0.43 mm/year and −1.67 ± 0.26 mm/year, respectively). Differences were observed only in the Dunkerque station (−0.18 ± 0.71 in the cited study, +0.45 ± 0.26 mm/year in this study).
The values of crust movement determined with the use of Equation (1) are shown in Figure 13a,b. Crust movements were well above zero in 12 tide gauges (based on data for 1993-2017). Tide gauges were below zero in Lampedusa and Sassnitz, which could be attributed to a greater decrease in average water levels calculated based on TG data rather than SA data. According to [13], the uncertainty of velocity determinations based on GNSS data is three to four times lower on average in comparison with combined satellite radar altimetry minus TG data which ranges from 0.3 mm/year to 3.0 mm/year, with a median of 0.80 mm/year. The uncertainty of trends for time series of daily SA data ranges from 0.04 mm/year to 0.09 mm/year, and it is higher when monthly data are used (around 0.40 mm/year). The uncertainty of trends based on the entire TG time series (monthly data) ranges from 0.03 mm/year (time series covering approximately 200 years) to 1.59 mm/year (time series covering approximately 20 years). The adopted method for verifying the velocity of GNSS stations was based on radar interferometry data generated with the use of the PSI approach. The results indicate that the velocities calculated in selected stations (TGBF, COUD, SCOA, TERS, SAS2 and SASS) based on SAR data differ in sign relative to the values calculated based on GNSS data. As shown in Figure 14a,b, the velocities calculated based on GNSS and SAR data differ in stations located on the North Sea and the Bay of Biscay. A considerable difference was observed in the SCOA station (−3.67 mm/year) where the standard error of velocity determination reached ±0.39 mm/year in the PSI method and ±1.92 mm/year in the calculations based on GNSS data. A large discrepancy was also noted in the ACOR station (−2.97 mm/year), where the standard error reached ±0.35 mm/year in the PSI approach and ±0.23 mm/year in the calculations based on GNSS data. In the TGBF station, the difference in the velocity was +3.84 mm/year with a standard error of ±0.30 mm/year in the PSI method and ±0.30 mm/year in the calculations based on GNSS data. On the other hand, the velocities calculated based on GNSS and SAR data are more similar in stations located on the Baltic Sea and the Mediterranean. The smallest difference was observed in the SWIN and LAMP stations. In the SWIN station, the difference in the velocity was +0.07 mm/year with a standard error of ±0.45 mm/year in the PSI method and ±0.18 mm/year in the calculations based on GNSS data. In the LAMP station, the difference in the velocity was −0.10 mm/year with a standard error of ±0.35 mm/year in the PSI approach and ±0.05 mm/year in the calculations based on GNSS data. The remaining velocity values resulted from the calculation of relative displacements in the PSI method as well as from differences in the analysed data sets and time periods.

Conclusions
TG and SA data were used to determine variations in sea level in the coastal zone. The results point to considerable differences in sea level trends, and they indicate that local and regional processes influence TG time series The results presented in the study point to land subsidence in the south and southwest along the European coast. These observations were confirmed by the data from permanent GNSS stations and InSAR systems. The InSAR technology helps with the determination of the vertical movements of the Earth's crust (deformation) in the areas surrounding GNSS stations.
Different measurement techniques have different real accuracies in determining the vertical movements of the earth's crust. The standard error from the combination of TG and SA data ranged from ±0.06 mm/year to ±0.21 mm/year at a long-time resolution and from ±0.31 mm/year to ±1.59 mm/year at a short-time resolution. The use of an unprecedented number of PSInSAR datasets allowed the estimation of a new level of accuracy in the range of ± 0.20 mm/year to ±0.49 mm/yy. The accuracy of the vertical movements of the Earth's crust from GNSS data ranged from ±0.04 mm/year to ±1.92 mm/year. Comparing the PSInSAR results required the development of a new procedure: the workflow of determining the value of vertical movement velocities in GNSS stations based on InSAR data.