Non-tidal loading of the Baltic Sea in Latvian GNSS time series


 The objective of this study is to investigate the effect of the Baltic Sea non-tidal loading in the territory of Latvia using observations of the GNSS continuously operating reference stations (CORS) of LatPos, EUPOS®-Riga, EPN and EstPos networks. The GNSS station daily coordinate time series obtained in a double-difference (DD) mode were used. For representation of the sea level dynamics, the Latvian tide gauge records were used. Performed correlation analysis is based on yearly data sets of these observations for the period from 2012 up to 2020. The approach discloses how the non-tidal loading can induce variations in the time series of the regional GNSS station network. This paper increases understanding of the Earth’s surface displacements occurring due to the non-tidal loading effect in Latvia, and is intended to raise the importance and necessity of improved Latvian GNSS time series by removing loading effects.


Introduction
The thorough analysis of GNSS time series provides the velocity field, which can indicate the vertical and horizontal crustal motions. The detailed models and algorithms for eliminating the influence of different geodynamic effects on geodetic station positions are given in the IERS (International Earth Rotation and Reference Systems Service) Conventions [1] and supported by the post-processing software packages. The correction for non-tidal ocean loading (NTOL) is usually not applied, but is an important condition for estimating highly reliable rates of Earth's crust movements in the coastal regions.
Post-glacial rebound is the most notable geodynamic process in Fennoscandia. The effect of post-glacial land uplift in the territory of Latvia is relatively small: according to the semi-empirical deformation model NKG2016LU_abs for absolute land uplift [2,3], the range of velocities is from 0 to +2 mm/yr [4]. In the case of land uplift relative to the mean sea level [5], another scenario can be observed: velocities have also negative values. An example of the future development of the Baltic Sea within the thousands of years based on geodynamic model was given [6]. The authors marked that in contrast to the land emergence in the northern parts of the Baltic Sea, land at the southern coastline drowns due to sea level rise.
The biggest uncertainty for the future scenarios is related to the sea level rise, which attracts the scientific interest today. No less important factor is to obtain velocity field eliminating the effect of the Baltic Sea non-tidal loading in GNSS time series, especially in the case of Latvia, where it is critically important to understand current situation of the geodynamic processes and future scenarios of the shoreline changes.
In [7] the authors first looked at the effect of NTOL on geodetic station height coordinate using output from an Ocean General Circulation Model (OGCM). More recently [8,9], there was demonstrated stronger impact of the NTOL on GPS height time series, thus justifying modern reanalysis with the best "global" data available at that time. Using height coordinate time series from a network of 344 globally distributed GPS reference stations in [8], there was examined the spatial variability of NTOL and found that the largest displacements of the Earth's surface due to this effect are expected to be found near semi-enclosed basins, where the redistribution of water due to wind or atmospheric pressure forcing is hindered by the local coastal geometry and bathymetry.
Several authors reported on the effect of NTOL on regional GPS networks [10,11,12,13]. In [10], there was investigated the effect of NTOL at 4 sites of tropical Pacific islands. Using sea-level data, the authors found an inverse correlation between sea-level and the GPS height coordinates at these sites. The effect of NTOL on height coordinates for 17 GPS sites around the southern North Sea region was analyzed in [11]. The authors reduced the scatter on all the height time series by 20-30 %, and highlighted the usefulness of a high-resolution model compared to a low-resolution global model. In [12] the tide gauge data were used to estimate NTOL effect on the vertical coordinates of 7 GPS stations on the Baltic Sea. The authors found that removing the computed loading from the GPS time series reduces the standard deviation, but neither for all series nor for all combinations of loads; the maximum reduction they found: 23 % at Metsähovi, 43 % at Mårtsbo and 34 % at Visby.
The capability of GRACE gravity solutions to recover the mass variation in the Baltic was investigated in [13]. The authors concluded that due to the small areal extent and complex shape of the studied geographical region, the Baltic Sea is probably close to the limit of detection for GRACE. They found that mass variations on monthly scales can be well described by surfaces fitted to tide gauge measurements alone.
All these studies are focused on the influence of NTOL on vertical component. The exceptions are three studies [14,15,16]. In [14] and [15] the authors detected the 3-D deformation during a storm surge events in the North Sea with sub-daily resolution for the one-month and twomonth time series. The 3-D effect of the Baltic Sea loading in GPS daily time series using data of 1.5 years was studied in [16]. In order to attenuate the common-mode type effects of other load signals in the GPS time series in [16], the authors formed differences of both coordinates and of modelled load for station pairs. The results have shown that due to correcting for the loading, the variance of time series reduces by up to 56 % in the East component and up to 30 % and 16 % in North and Up; the reduction depends on the station pair. For the pairs, where the Latvian GNSS station RIGA presents, using the PPP time series provided by the Jet Propulsion Laboratory, the authors found: 7 %, 15 % and −1 % in North, East and Up respectively (MAR6-RIGA); 12 %, 39 % and 14 % in North, East and Up respectively (VIS0-RIGA).
This study shows the impact of NTOL in three coordinate components and for longer observation period (9 years) using a relatively dense network of GNSS stations (in scope 38 stations) in the territory of Latvia. The approach discloses more obviously how the effect of NTOL can induce seasonal variations in the time series of the regional GNSS station network, and how regionally related results are affected from year to year.

Sea level variations at the Latvian coast
According to [17], the long-term variations of the Baltic Sea (within months and years) are caused by persistent winds redistributing water between the Atlantic Ocean (the North Sea) and the Baltic Sea, and affecting the Baltic as a whole. The short-term sea level variations (less than one month) are mainly internally driven variations, with maximum amplitudes in the far north and the far south, and a nodal line close to Stockholm in the middle. The range of tidal variations of the Baltic Sea is in the order of centimetres, but due to winds and complex shape of gulfs and islands, the range of non-tidal variations can reach 3 m on the coasts of gulfs [18,19,20,21].
The Latvian coast is also located close to the nodal line of the Baltic Sea level variations as in the case of Stockholm, where in spite of strong impact of the open sea flows due to wind, short-term variations are nearly eliminated [17]. Latvia has 500 km long coastline. The west coast of Latvia is washed by the open sea up to north (Cape Kolka), where meets waters of the Gulf of Riga.
The Gulf of Riga is a shallow semi-enclosed basin separated by Estonia's Saaremaa Island and connected to the Baltic Sea by two relatively narrow and shallow straits towards north and west. According to [22], the highest water levels in the interior of the Gulf of Riga are developed under the joint impact of three major drivers: water volume of the entire Baltic Sea, water occasionally pushed by a sequence of cyclones into the Gulf of Riga for 1-2 days and local storm surges with a typical duration of a few hours. Each of these drivers may add about 1 m to the resulting water level.
The Latvian Environment, Geology and Meteorology Centre (LEGMC) is maintaining totally 9 tide gauges: 2 of them are located on the coast washed by the open sea, but the others are located on the coast of the Gulf of Riga. Records of 7 tide gauges were used for this study: Liepāja, Ventspils, Kolka, Mērsrags, Lielupes grīva, Daugavgrīva and Salacgrīva (locations are shown in Figure 4). Sea level data are daily mean values of hourly records. Figure 1 presents daily sea level variations at the selected tide gauges from 2012 to 2020. The time series have an increase in the amplitude during the winter season. This effect was pointed out by [23]. He found a significant secular increase of the seasonal sea level variation in the Baltic Sea in December-January and a significant secular decrease in February-March, explaining this by changes of wind conditions over the transition area between the North Sea and the Baltic Sea. However, there is some inconsistency through the years in the Latvian tide gauge data for the observation period: the time of increase and decrease of sea level is changing.
According to the daily mean values of sea level records for the period of 9 years, the range of sea level variations at the tide gauges Liepāja, Ventspils and Kolka is close to 1.5 m; at Mērsrags, it is 1.8 m; but at the tide gauges Lielupes grīva, Daugavgrīva and Salacgrīva, this range runs up

GNSS station time series
The time series of 38 permanent GNSS stations were used for this study. Stations belong mainly to the LatPos and EUPOS ® -Riga networks, as well as to the EUREF Permanent GNSS Network (EPN) (stations IRBE and RIGA; RIGA is also station of the International GNSS Service (IGS)) and EstPos network (stations RUHN and IKLA).
Daily coordinate time series were computed for the period of 9 years (from 2012 up to 2020) in a double-difference (DD) mode using Bernese GNSS software v5.2 [24]. The final CODE (Center for Orbit Determination in Europe) precise orbits, Earth orientation and clock products, along with the CODE final ionosphere product, were used for the GNSS data processing. The dry Global Mapping Function (GMF) was used as the a priori troposphere model, while zenith path delay parameters were estimated using the wet GMF. In the processing the Phase-Based Widelane (L5) and Quasi-Ionosphere-Free (QIF) ambiguity resolution strategies were applied. A cut-off elevation angle of 3°was selected. Station positions were corrected for the effect of solid Earth tides and the ocean tidal loading (FES2004 ocean tide model was used). Obtained positions are without corrections of the atmospheric tidal loading.
Since 2015, the combined processing of GPS and GLONASS observations was applied (only GPS observations were used before). Since 2020, the Galileo obser- Coordinates obtained in IGb08 and IGS14 were transformed to ETRF2000 using one-step transformation with 14 transformation parameters [25].
Atmospheric loading is out of focus of this study. As CORS network is relatively dense and stations distributed over the plain territory of about 65000 km², the atmospheric surface pressure anomaly can be considered as near-uniform, it produces a near-uniform vertical deformation and near-zero horizontal deformations.
During whole observation period the CORS network is not consistent [26]. From 2012 to 2015, the observations of 28 stations were used on a regular basis. Since 2016, the observations of 2 Estonian stations have been introduced together with the data of new Latvian stations. Station receivers and antennas have been changed, as well as few stations have stopped their operation. The usage of yearly correlations allows to neglect the unavailability of some stations for the whole observation period, and to cover the territory of Latvia with maximum number of stations, observing the changing NTOL effect.
Initially, the influence of the Baltic Sea loading on the Latvian GNSS station coordinate time series was observed analyzing cumulative horizontal velocities, and then representing annual horizontal velocities. Some of stations had clearly visible common orientation of annual displacement [27]. This was the reason to perform correlation analysis using GNSS station time series and tide gauge observation data.

Results of the correlation analysis
The correlation coefficients were obtained using GNSS station position residuals in North, East and Up components and sea level data for each year and GNSS station separately. The maximum values of correlation coefficients (among 7 tide gauges) were used as output. Figures 2 and 3 present statistics of the correlation analysis for East and Up components at locations of the GNSS stations for each observation year. The results have a high probability of being statistically significant, as the number of daily data pairs during one year is higher than 300, with few exceptions. The number of GNSS CORS observation days and the number of stations per year are summarized in the Table 1. Names of the GNSS stations with less than 300 observation days per year are also listed.
The correlation greater than 0.20 for 300 pairs of data sets has less than 0.05 % probability of occurring randomly, and the correlation greater than 0.20 for 365 pairs of data sets has less than 0.01 % probability of occurring randomly. However, if the correlation is 0.10 with 300 pairs, the probability of occurring randomly becomes around  8 % [28]. For example, in Table 2       For the East component, sites display negative correlations (with two exceptions) or very weak positive correlations in few cases. The inexplicable fact is the inverse correlation at the station SLD1 during the whole observation period. Possibly, this can be explained by some geological effect. The other case with the positive correlation is observed in 2015 at station SALP. This station is located close to the Riga hydroelectric power plant and its position could be additionally affected by changing water masses in the river Daugava and reservoir of the plant.
In the case of North component (see Figure 5C), correlation coefficients are changing from very weak to positive moderate correlations, with outstanding values for some cases. For example, in 2020 the observations are strongly correlated at the EPN/IGS station RIGA (0.68) (with the sea level residuals of the tide gauges Lielupes grīva and Daugavgrīva).
In this context, the geographical location of the territory of Latvia should be considered. The main water masses are concentrated in the west for the territory of Latvia with shallow connections from the Baltic Proper to the Gulf of Riga [29]. Therefore, uniform water rise in the whole of the Baltic Sea will induce maximum downward vertical deformation (subsidence) in the west coast of Latvia [16]. Obtained results, i. e. negative correlations in the Up component at the western sites, have confirmed this fact. As well as the positive correlations at the eastern sites indicate the uplift in the eastern part of Latvia. In turn, in the East component this leads to the negative correlations. Periodically changing correlation in the North component can be admitted in this sense as well, because it can be less affected by the uniform change of water mass in the whole of the Baltic Sea, but more affected by water variations in the north and the south, and in the Gulf of Riga itself.

Variations in the GNSS time series
In GNSS time series about 40 % of the observed variations can be explained by the joint contribution of seasonal surface mass redistributions due to varying atmospheric, hydrological and ocean masses [30]. This leads to the derived rate uncertainties larger than expected [31,32]. In Figure 6A the time series with the negative correlation (−0.71) in the Up component at station IRBE are  shown (year 2020). Red dashed line presents the GNSS station coordinate residuals (with trend line) and grey line with an area depicts the sea level residuals (at tide gauge Ventspils) divided by factor 35. In this case, the annual rate of the IRBE station reaches 10 mm/yr. However, in 2013 (see Figure 6B) the rate of this station is −1 mm/yr. As it can be seen from the graphs, such difference in the annual rates is caused by the loading effect mainly from the changing water masses in the Baltic Sea.
Additionally, according to the graphs in Figure 6, more evident outliers in the GNSS time series are observed at time spans with the maximum sea level during the year. Therefore, in this context, the peaks in the GNSS time series usually occurred during the winter time and previously assumed to be due to snow coverage of GNSS antenna [33], need further investigation.
In . This proves that for the relatively dense GNSS network of the territory of Latvia, the sites are affected similarly in the Up component within small local parts, and in a different way within whole territory, showing inconsistency in rates for the same observation period. Figure 8 displays the time series with the strong negative correlation (−0.78) in the East component at station VAL1 (year 2018). In this case, the sea level residuals (at tide gauge Ventspils) divided by factor 250. The strong correlation between two data sets is obvious.
According to [16], one-meter layer of water in the whole of the Baltic Sea can induce vertical deformation about −16 mm for the territory of west coast of Latvia, and about 4 mm for the horizontal displacement. The time series of the VAL1 station show approximately the same range of caused displacements in the East component, when the range of sea level residuals at the Latvian tide gauges is from 0.98 m (Ventspils) up to 1.38 m (Salacgrīva). Concerning vertical displacement, the range is higher: for the IRBE station, it is about 25 mm, and for the eastern station BALV, it is 16 mm (the range of sea level residuals at the Latvian tide gauges in 2020 is from 0.99 m (Ventspils) up to 1.42 m (Daugavgrīva)). Furthermore, in contrast to the 3D deformation model for the whole of the Baltic Sea presented in [16], the increase of water masses results in positive correlations at sites of the eastern part of Latvia, i. e. the obtained results show the uplift in the eastern part of the territory under the loading of water masses.

Conclusions
This paper raises the importance and necessity of im- For the East component, sites display negative correlations with few exceptions. Correlation in the East component is most prominent for stations located close to the eastern coast of the Gulf of Riga. The observations at 10 sites are strongly correlated in 2018, and at 6 sites in 2020. The maximum correlation coefficients in the East component are −0.78 (VAL1), −0.73 (LIMB) and −0.80 (IKLA). In the North component correlation coefficients are changing with outstanding values for some cases. The RIGA station has shown the maximum correlation in 2020, which is 0.68.
The statistics of the correlation analysis show that despite the weak and very weak correlations, the correlation results are regionally related and consistently changing from year to year. The GNSS time series reveal the loading effect of the changing water masses, which leads to the inconsistency in the annual rates of GNSS stations within the dense Latvian network. The study also highlights the importance in the planning of geodetic measurement campaigns considering the Baltic Sea non-tidal loading effect in Latvia.
The evident outliers in the GNSS time series during the winter time need further investigation. Previously, they were assumed to be due to snow coverage of GNSS antenna. As the maximum sea levels of the Baltic Sea are observed during the winter season, the peaks in the time series can be caused by the changing water masses.
According to the time series of station IRBE for the year 2020, the range of caused vertical displacements can reach 25 mm. Stations IRBE and RIGA are EPN stations, RIGA belongs also to the IGS; therefore, they are sites of the highest importance. Observations at these sites have shown moderate and strong correlations. This should be noted, because the data of these stations are often used for GNSS data processing in the Baltic Sea region.
Atmospheric loading is out of focus of this study. As the GNSS network is relatively dense, the atmospheric surface pressure anomaly can be considered as near-uniform, it produces a near-uniform vertical deformation and nearzero horizontal deformations. The research will be continued using modelled loading time series, considering atmospheric loading. Hopefully, further research will allow to reduce the variance of the Latvian GNSS coordinate time series correcting for the Baltic Sea non-tidal loading, and ensure computation of reliable horizontal and vertical velocity fields.