# Analysis of GNSS Displacements in Europe and Their Comparison with Hydrological Loading Models

Alexandre Michel, Alvaro Santamaría-Gómez, Jean-Paul Boy, Félix Perosanz, Sylvain Loyer

► To cite this version:

Alexandre Michel, Alvaro Santamaría-Gómez, Jean-Paul Boy, Félix Perosanz, Sylvain Loyer. Analysis of GNSS Displacements in Europe and Their Comparison with Hydrological Loading Models. Remote Sensing, 2021, 13 (22), pp.4523. 10.3390/rs13224523. hal-03466945

HAL Id: hal-03466945

https://hal.science/hal-03466945

Submitted on 9 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Analysis of GNSS Displacements in Europe and Their Comparison with Hydrological Loading Models

Alexandre Michel <fnref n="1" />*, Alvaro Santamaría-Gómez <fnref n="2" /> ID, Jean-Paul Boy <fnref n="1" /> ID, Félix Perosanz <fnref n="2" /> ID and Sylvain Loyer <fnref n="3" />

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Michel, A.; Santamaría-Gómez, A.; Boy, J.-P.; Perosanz, F.; Loyer, S. Analysis of GNSS Displacements in Europe and Their Comparison with Hydrological Loading Models. *Remote Sens.* **2021**, *13*, 4523. https://doi.org/10.3390/rs13224523

Academic Editor: Chung-yen Kuo

Received: 11 October 2021

Accepted: 8 November 2021

Published: 10 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Thanks to the increasing number of permanent GNSS stations in Europe and their long records, we computed position solutions for more than 1000 stations over the last two decades using the REPRO3 orbit and clock products from the IGS CNES-CLS (GRGS) Analysis Center. The velocities, which are mainly due to tectonics and glacial isostatic adjustment (GIA), and the annual solar cycle have been estimated using weighted least squares. The interannual variations have been accounted for in the stochastic model or in the deterministic model. We demonstrated that the velocity and annual cycle, in addition to their uncertainties, depend on the estimation method we used and that the estimation of GPS draconitic oscillations minimises biases in the estimation of annual solar cycle displacements. The annual solar cycle extracted from GPS has been compared with that from loading estimates of several hydrological models. If the annual amplitudes between GPS and hydrological models match, the phases of the loading models were typically in advance of about 1 month compared to GPS. Predictions of displacements modelled from GRACE observations did not show this phase shift. We also found important discrepancies at the interannual frequency band between GNSS, loading estimates derived from GRACE, and hydrological models using principal component analysis (PCA) decomposition. These discrepancies revealed that GNSS position variations in the interannual band cannot be systematically interpreted as a geophysical signal and should instead be interpreted in terms of autocorrelated noise.

**Keywords:** GNSS; IGS REPRO3; tectonic velocity; hydrological loading; principal component analysis

## 1. Introduction

Time series of station coordinates derived from global navigation satellite systems (GNSS) have been used for decades to investigate various geophysical phenomena. Horizontal components have mostly been used to estimate tectonic activity (plate motion [1,2] or seismic events [3,4]), and vertical components have been used for seasonal signals due to mass redistribution (hydrological loading [5–7], atmospheric loading [8], nontidal ocean and atmospheric loading [9], and ice–snow loading [10,11]). All components are also used to determine reference frames, such as the International Terrestrial Reference Frame (ITRF) [12]. Aside from these signals, GNSS position time series also contain large broadband variations of unknown origin, typically represented by the combination of power-law (PL) and white noise (WH) models [13], which impacts the determination and interpretation of other parameters, especially the velocity and its uncertainty [14–16], but also possibly the seasonal signal [17]. The origin of PL noise in GNSS position time series may have several sources: orbit mismodelling [18], tropospheric delay [19], instability of the station monumentation [20,21], and multipath [22], etc. Moreover, [23] shows that the assessment

<footnote n="1">¹ Institut Terre et Environnement Strasbourg (ITES, UMR7063: Université de Strasbourg, CNRS, ENGEES), 5 Rue René Descartes, 67084 Strasbourg, France; jeanpaul.boy@unistra.fr</footnote>

<footnote n="2">² Geosciences Environnement Toulouse (GET, UMR5563: Université de Toulouse, CNES, CNRS, IRD, UPS), 14 Avenue Edouard Belin, 31400 Toulouse, France; alvaro.santamaria@get.omp.eu (A.S.-G.); felix.perosanz@cnes.fr (F.P.)</footnote>

<footnote n="3">³ Collecte Localisation Satellites, 11 Rue Hermès, Parc Technologique du Canal, 31520 Ramonville Saint-Agne, France; sloyer@groupcls.com</footnote>

<footnote>* Correspondence: alexandre.michel@unistra.fr; Tel.: +33-(0)3-68-85-01-09</footnote>
of noise content in GNSS position series is heavily impacted by the estimation of position discontinuities, which make the interpretation even more difficult and compromise the validity of long station records. Taking into account that the broadband noise variations are also correlated spatially [24,25], the proper interpretation of the interannual variations in GNSS position time series and their common modes across a region remains a challenging task. The deformations that are observed in GNSS position time series may depend not only on the geological nature of the crust (karst aquifers [26]) and its thermal deformation [27,28] but also on human activities (groundwater pumping [29] and mining [30,31]). We can suspect some apparent deformation related to the stability of the station monument. Finally, in [32], the authors show that the power spectral density (PSD) from GNSS position time series contains harmonics of 1.04 cpy (cycle per year), known as the GPS draconitic oscillation. The origin of the draconitic oscillation and its harmonics in GNSS time series is still not well established but is most likely due to satellite orbit mismodelling [32,33]. The first draconitic harmonic is very close to the annual solar cycle (1 cpy), so the estimated amplitudes of these two waves are strongly correlated. Thus, the interpretation of the seasonal signal, and especially the annual signal, is highly dependent on the draconitic signal estimation.

In this study, we used specific IGS REPRO3 GNSS satellite orbit and clock products computed by the CNES-CLS team on behalf of the Groupe de Recherche en Géodésie Spatiale (GRGS, https://grgs obs-mip fr/ (accessed on 7 November 2021)) in France. They were used to calculate station positions in precise point positioning (PPP) mode using the GINS software developed by CNES [34]. We computed the daily position solutions of more than one thousand stations mainly located in Western Europe and Scandinavia. We extracted the linear velocity and seasonal signals together with their uncertainties to validate the accuracy of the estimation of the new products using different methods. We also assessed the spectral contents of the time series, especially in terms of the noise level at the interannual band.

Among the current techniques used for determining the parameters of the model, we can cite the weighted least squares method (WLS), Kalman filtering, (multi-) singular spectrum analysis (M-SSA), and the Wiener filter, all reported and analysed in [35]. We demonstrated the importance of the contribution of the interannual variation in GNSS time series by using two different methods of estimation.

The first method is based on a WLS fit using an optimized covariance matrix obtained from an MLE (maximum likelihood estimation) of a PL noise model [36]. The interannual signal is considered as time-correlated noise by the estimator and is directly propagated in the parameter uncertainties given by the WLS.

The second method is a WLS method where only white noise is accounted for in the covariance matrix, but in which we modified the classic deterministic model (linear velocity and seasonal signal) to account for interannual variation as polynomials. Adding polynomials in the estimated model should reduce the parameter uncertainties given by WLS as it should better fit the data. However, given these polynomials, we can compute the instantaneous velocity and deduce a statistical uncertainty on the linear velocity, which reflects the variability of the instantaneous velocity around its average value. This method should provide a more realistic empirical velocity uncertainty than that directly given by the WLS in the case which is known to be far too small when only WH noise is accounted for [15,23].

Then, the fundamental difference between these two approaches (respectively, WLS + MLE and WLS + polynomials + WH) is the presumed nature of the interannual variation (respectively, stochastic and deterministic) and, consequently, the estimation of the parameter uncertainties. Moreover, we highlight the need to estimate both draconitic and solar periods to prevent biasing estimates of the annual cycle. Finally, in order to search for the true geophysical signals in interannual variations, the extraction of spatial common modes from the signal is usually carried out [37]. This can be conducted using different techniques: principal component analysis (PCA) [38–40], independent component analysis (ICA) [41], or robust statistical methods [42]. Since the sources are not necessarily independent, it is reasonable to choose PCA in order to compare the interannual signal from GNSS with the
signals derived from hydrological models, GRACE (Gravity Recovery and Climate Experiment), and the GRACE Follow-On time variable gravity field [43]. Although recent studies on interannual variations show correlations between GNSS and some global hydrology models [44], we show here that models can also differ significantly from each other and that the interannual signal from GNSS should be interpreted very carefully.

## 2. Materials and Methods

### 2.1. CNES-CLS Multi-GNSS Orbits and Clock Products

The GNSS satellite orbit and clock files used in this study were generated as part of the GRGS participation in the third International GNSS Service (IGS) reprocessing campaign effort (REPRO3) to contribute to the realisation of ITRF2020. These products (referenced here as MG3) are the result of a homogeneous reprocessing of multi-GNSS data (GPS, GLONASS, and Galileo) between 2000 and 2020 using the zero difference with integer ambiguity fixing method described in [34,45]. We used up-to-date models listed in Table 1, following the IGS recommendations (http://acc.igs.org/repro3/repro3.html (accessed on 7 November 2021)). The MG3 reference frame solutions were first evaluated and then preliminarily combined by the IGS [46]. The station network includes a set of about 300 selected IGS stations that are distributed over the globe (Figure 1). The number of available satellites for each constellation and the daily number of available stations varied over the observing period as shown in Figure 2. The CNES-CLS MG3 products are available in the CDDIS (Crustal Dynamic Data Information System) archive at https://cddis.nasa.gov/archive/gnss/products/www/repro3/GRG6RE3FIN*.gz (accessed on 7 November 2021).

<table><caption>**Table 1.** Dynamical and loading models used for MG3 products.</caption><thead><tr><th></th><th>MG3 Products</th></tr></thead><tbody><tr><td>Gravity field</td><td>EIGEN-GRGS.RL04.MEAN-FIELD [47]</td></tr><tr><td>Ocean tides (gravity)</td><td>FES2014b (Finite Element Solution) [48]</td></tr><tr><td>Planet ephemerides</td><td>de421bdlf.ad [49]</td></tr><tr><td>Relativistic acceleration</td><td>Schwarzschild and geodetic precession and Lense-Thirring</td></tr><tr><td>Antex</td><td>IGSR3.atx [50]</td></tr><tr><td>Mean pole (C21/S21)</td><td>IERS conventions (from geopotential model) [51]</td></tr><tr><td>Subdaily EOP model</td><td>[52]</td></tr><tr><td>Atmospheric tides (S1/S2)</td><td>[53]</td></tr><tr><td>Ocean tide loading</td><td>FES2014b [48]</td></tr><tr><td>Centre of mass correction</td><td>FES2014b [48]</td></tr><tr><td>Solid tides (station)</td><td>IERS conventions [51]</td></tr><tr><td>Reference frame</td><td>IGS_R3 [IGSMAIL-8026]</td></tr><tr><td>Galileo ponderation</td><td>3.5 cm/1 m</td></tr><tr><td>GPS ponderation</td><td>3.5 mm/60 cm</td></tr><tr><td>GLONASS ponderation</td><td>3.5 cm/2 m</td></tr></tbody></table>

### 2.2. Time Series Analysis

#### 2.2.1. Selection of Stations

We computed the GNSS position time series of 1077 stations over Europe using the MG3 products and the GINS software in PPP mode. We used the same models as the ones listed in Table 1 in order to maintain consistency between the products used and the individual station processing. For the tropospheric delays, we use the global mapping function (GMF) [54] tropospheric model and the global pressure and temperature empirical function GPT2 [55]. The numerous agencies and organisations providing the raw GNSS data in RINEX format are listed in the acknowledgements section.

The selected stations have a minimum time span of five years and a completeness of 50%. Some exceptions of a time span between three and four years have been considered for recent stations (ending after 2019) but with a more selective completeness criterion (minimum 70%). The distribution of stations and the statistics of the time series are shown
in Figure 3. The analysed network is homogeneous and dense in Great Britain, France, Spain, and Italy, while in the remaining European countries, the station distribution is not as dense but still spatially homogeneous.

<figure><img src="image_5.png" /><figcaption>Figure 1. GNSS station network used for MG3 product computation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2. (a) Number of satellites and (b) number of stations used in MG3 product computation as a function of time.</figcaption></figure>

### 2.2.2. Parameter Estimation

For each time series, we determined different parameters corresponding to the classical model used in GNSS position adjustment:

$$ 
\begin{aligned}
y(t) = & y_0 + v(t - t_0) \\ 
& + \sum_{m=1}^{2} a_m \cos(m\omega_a t) + b_m \sin(m\omega_a t) \\ 
& + \sum_{n=1}^{8} c_n \cos(n\omega_d t) + d_n \sin(n\omega_d t) \\ 
& + \sum_{k=1}^{n} s_k \mathcal{H}(t - t_k) + \varepsilon(t),
\end{aligned} 
\quad (1) 
$$
where $y_0$ and $v$ are, respectively, the intercept (at epoch $t_0$) and the linear velocity of the station; $\omega_a/2\pi = 1$ cpy is the solar frequency; and $\omega_d/2\pi = 1.04$ cpy is the draconitic frequency. The offsets in the time series, which can be linked to known events (earthquakes and antenna or receiver changes) but are also mostly of an unknown origin [56], were modelled with a Heaviside function $H$. Although automatic methods exist for offset detection [57], we chose to do it manually for all the stations using an initial seismic database and monumentation change data provided by sitelog files. We calculated the residuals of the time series by removing the linear trend and offsets from this preliminary database. We first removed major outliers from these residual time series using automatic detection of the largest outliers with a criterion of $5\sigma$, where $\sigma$ is the standard deviation. Then, we visually checked the precleaned residual time series in order to validate the outliers/offsets removed during the first step and refine the database if needed (adding or removing offset positions and outliers). The manual detection of offsets and outliers is essential even after applying any automatic procedure. Indeed, regardless of its robustness, the automatic process could eliminate data, especially those corresponding to unmodelled geophysical signals. However, manual checking remains a subjective method that depends on the criteria set by each analyst. The term $\varepsilon(t)$ is the stochastic part of the model that has to be modelled and then estimated.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** (a) GNSS station network and some characteristics in terms of (b) station time span, (c) station availability over time, and (d) daily solution in each station.</figcaption></figure>

We used the CATS software developed by [36], which uses WLS estimation of the deterministic parameters of Equation (1) along with MLE determination of the stochastic component. We worked with week-averaged time series in order to obtain a good compromise between the computation time and time resolution of the series. We modelled the stochastic part using white noise (WH) and power-law noise (PL), where the spectral index was also estimated using CATS. Nevertheless, for some stations, we used only PL noise without WH noise in the stochastic part since the downsampling to weekly time series removes a large part of the WH noise contribution that can no longer be estimated properly
by the software. We computed one solution with the complete model of Equation (1), called *cats_d*, and one without the estimation of the draconitic periods, simply called *cats*.

### 2.2.3. Interannual Polynomial Model

Since interannual variations play an important role in the determination of the model parameters and their uncertainties [58], especially in GNSS, instead of capturing these variations in the stochastic model, we alternatively chose to directly model most of the interannual variation within the deterministic model function. In addition to the different terms in Equation (1), we added degree 2 and degree 3 polynomial terms. This model is called *tiasd* (for trend interannual (semi-)annual steps draconitic). The low degree polynomial function was used to fit the long-term (relative to the length of the time series) interannual signal. We adjusted the model with the WLS method using the least squares algorithm of [59] implemented in the Python programming language. The stochastic term $ε(t)$ was then reduced to a simple WH noise estimation. We first obtained the optimal parameters of the model and calculated the residuals of the time series. These residuals could still have contained a faster interannual signal (between a year and a decade, depending on the length of the time series) and were then modelled by a polynomial of degree 4 to 11. The degree of the adjusting polynomial was chosen such that it is always inferior to the time series length (in year). For time series with a low completeness (<70%) and high standard deviation (beyond 7 mm in the vertical component), we determined the polynomial degree by analysing each time series individually while verifying that the polynomial did not overfit the series. The polynomial fitting procedure in these two steps (low degrees first and then high degrees on residuals) was decided upon to avoid the strong correlation that could have occurred between high-degree polynomials and the sinusoidal and offset terms. In order to reconstruct the full interannual variations, we added the slow and fast polynomial contributions. We calculated the instantaneous velocity by taking the first derivative of the resulting polynomial function and deducing a value of linear velocity by taking the mean of instantaneous velocities. The uncertainty on this definition of linear velocity should reflect the magnitude of the interannual variation. Then, in addition to considering only the dispersion of the positions in the calculation of the velocity uncertainty, as in the case of WH-only WLS estimation without polynomial adjustment, we also took into account the dispersion of the instantaneous velocity, which should give a more realistic uncertainty on the linear velocity. Thus, we define the uncertainty as the standard deviation of the instantaneous velocity divided by the square root of the time series length [60].

### 2.3. *Hydrological Loading Computations*

We computed surface displacements due to continental water storage variations using two state-of-the-art global hydrological models, GLDAS-2.1 (Global Land Data Assimilation System)/Noah [61] and MERRA-2 (Modern-Era Retrospective Analysis for Research and Applications) [62] land component, and estimates derived from the latest (RL06v1.0) GRACE and GRACE Follow-On iterated global mascons from the NASA Goddard Space Flight Center [63]. Among other differences, the GLDAS-2.1 model includes soil moisture and snow and canopy water, whereas MERRA2 only includes soil moisture and snow.

We used the classical Green's function approach [64], assuming a spherically symmetric nonrotating elastic isotropic (SNREI) Earth model, using PREM [65] rheological parameters. More details of the loading computations can be found in [9,66]. In particular, we ensured the total water mass conservation of the hydrological models by adding/removing a uniform ocean layer to compensate for any lack/excess of water over land.

All the loading time series are available at the EOST loading service (http://loading.u-strasbg.fr (accessed on 7 November 2021)).

## 3. Results

In this section, we present the parameters that we obtained for three different estimated solutions described in the previous section (i.e., *cats_d*, *cats*, and *tiasd*) along with further
analysis of the results. The properties of these three models are summarized in Table 2. Then, we can compare the effect of draconitic adjustment by comparing *cats_d* and *cats*, and we can compare the modelisation of interannual variations and stochastic parts by comparing *cats_d* and *tiasd*.

<table><caption>**Table 2.** Estimation properties of the three models considered in this article.</caption><thead><tr><th></th><th><em>cats_d</em></th><th><em>cats</em></th><th><em>tiasd</em></th></tr></thead><tbody><tr><td>Draconitic frequencies adjustment</td><td>Yes</td><td>No</td><td>Yes</td></tr><tr><td>Interannual variation</td><td>—</td><td>—</td><td>polynomials</td></tr><tr><td>Stochastic model</td><td>WH + PL noise</td><td>WH + PL noise</td><td>WH noise</td></tr></tbody></table>

### 3.1. Tectonic Velocity

The velocity maps of *cats_d* are presented in Figure 4. We used two scales (red and green arrows) in the horizontal map to distinguish the largest velocities in Greece and Turkey, while the yellow dots represent the stations with horizontal velocities lower than 1 mm/y. In order to provide the horizontal velocity field relative to the Eurasian plate, we removed the rotation of the Eurasian plate with the Euler pole coordinates estimated in [67] (lon<sub>p</sub> = -99.094(7)°E, lat<sub>p</sub> = 55.070(4)°N, ω = 0.261(1) × 10<sup>-6</sup> deg/y). This enabled us to compare our solution with previous regional studies [1] or the EPOS (European Plate Observing System) solution available at http://doi.osug.fr/data/public/GNSS_products/ Europe/ (accessed on 7 November 2021). The vertical velocity map is represented together with the ICE-6G_D GIA model developed by [68,69]. In general, there is a good match between the two, although some differences are noticeable, especially in Southern Italy, where other non-GIA geophysical signals may occur.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Velocity field of *cats_d* in (a) horizontal and (b) vertical direction, along with ICE-6G_D GIA model of [68]. The yellow dots in (a) are stations for which the horizontal velocity is lower than 1 mm/y.</figcaption></figure>

In Figure 5, we provide the velocity differences between *cats_d* and *cats*, *cats_d* and *tiasd*, and *cats_d* and *epos*. We observe that there are no significant differences between the velocity values of the two CATS estimations, while *tiasd* slightly underestimated the vertical velocities compared to *cats_d*. The distribution shown by the histogram remained Gaussian for both vertical velocity differences. Moreover, we do not observe any preferential direction
in the horizontal maps, such that we can make two assumptions. The first is that none of the estimation methods seem to be biased with respect to the others. The second is that we can consider the velocity differences as a random field. In order to compare our GNSS solution to other published GNSS solutions, we chose to compute the difference between *cats_d* and the EPOS solution for the common stations, as both networks are slightly different. We see that there are some systematic effects on both components, especially for the vertical component, whose distribution is shifted near 0.5 mm/y. These effects, being visible at large spatial scales, probably indicate a difference in reference frame realisation. In fact, the reference frames of the two solutions are obviously different since the IGS_R3 reference frame is specifically used in REPRO3 products. Nevertheless, when we removed this mean shift, the distribution seemed to follow the same behaviour as *cats_d*–*tiasd* (*cats_d* slightly overestimates the velocity compared to *tiasd* and *epos*).

<figure><img src="image_9.png" /><figcaption>**Figure 5.** (a,b) Difference in velocity field between *cats_d* and *cats*. The panel (c) represents the number of stations where the vertical velocity was, respectively, overestimated (red) or underestimated (blue) by *cats_d* compared to *cats*. (d–f) are the same for the difference between *cats_d* and *tiasd*, and (g–i), between *cats_d* and the EPOS solution.</figcaption></figure>
The uncertainties of the three estimated solutions are plotted in Figure 6. They are greater for *cats* since the remaining draconitic signal in the time series contributes to a greater dispersion of the series, increasing the uncertainty estimated by CATS. In addition, *tiasd* underestimated both horizontal and vertical uncertainties by a factor varying from 3 to 4 compared to *cats_d*. Consequently, evaluating interannual variations with a complete stochastic model (WH + PL noise) rather than with WH noise only + the polynomials model provides more realistic uncertainties on velocity. Indeed, the strong correlation between the estimated coefficients of the polynomials used to model the interannual signal could lead to a significant underestimation of the instantaneous velocity dispersion and, therefore, to an underestimation of the linear velocity uncertainty.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Uncertainty of (a–c) horizontal and (d–f) vertical velocity field for *cats_d*, *cats*, and *tiasd*.</figcaption></figure>

## 3.2. Annual Signal

We now focus on the annual solar cycle of the time series derived from the three estimated solutions. Among the known sources of seasonal variations in GNSS position time series, we can cite the hydrological loading of amplitude ~4 mm in Europe [70], thermoelastic deformation of the crust of amplitude ~1 mm [27,28], nontidal loading of amplitude < 1 mm in Europe [9], and thermal dilation of the GNSS monumentation with variable amplitudes. In addition to the velocity field, the annual cycle recovery also depends on the choice of the estimation method. Figure 7 represents the phase and amplitude of *cats_d* for the vertical component with their uncertainties. The horizontal components are not shown here because their amplitudes were at the level of the resolution of the technique (the uncertainties had the same order of magnitude as the signal) and could have been dominated by GNSS monument motion, which makes the interpretation very difficult.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** (a) Phase and (b) amplitude of *cats_d* along with their respective uncertainties (c,d) estimated using CATS software. The colour bar of the phase map indicates the month of maximum displacement towards the up direction.</figcaption></figure>

In spite of spatial variations in both amplitude and phase from Figure 7, the results are coherent with the expectations of finding maximum displacement in the mid-summer or beginning of fall, with 3 to 4 mm of mean amplitude. We can observe some typical patterns: for example, a gradient of amplitude over Great Britain and larger amplitudes in Eastern Europe and Scandinavia due to important snow covering (compared to Western Europe) and atmospheric loading signals caused by Siberian anticyclones, both in winter. In addition, the uncertainties are mostly below 1 month for the phase and 1 mm for the amplitude, which make the results suitable for proper interpretation.

As for velocities, we compared the results of *cats* and *tiasd* with respect to *cats_d* in Figure 8. The phase differences are coloured in red when the tested model is delayed compared to *cats_d* and in blue when it is in phase advance. Note that we computed the difference in amplitudes (resp. phase) and not the amplitude (resp. phase) of the
differences. We will evaluate both the effects of the interannual estimation method and draconitic adjustment in the annual solar cycle recovery.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** (a–d) Difference in annual solar cycle between *cats_d* and *cats*. Panels (b,d) represent, respectively, the number of stations where *cats* is in phase advance (blue) or delay (red) with respect to *cats_d* and the number of stations with larger (blue) or smaller (red) amplitude than *cats_d*. (e–h) are the same for the difference between *cats_d* and *tiasd*.</figcaption></figure>

For the difference between *cats_d* and *cats*, there are two distinct effects resulting from the estimation of the draconic frequencies which affect the amplitude and phase of the annual solar cycle. The first is the small phase delay of *cats*, which has also a smaller amplitude compared to *cats_d*. This is consistent with the fact that the draconitic cycles can disrupt the solar cycle and that the modulation between the draconitic first harmonic and the annual solar oscillation can bias the estimation of the latter. We also remark that the distribution of the differences in the histograms reflects the spatial variability of the draconitic terms. For the difference between *cats_d* and *tiasd*, there are also two distinct effects resulting from the estimation method of interannual variations which affect the amplitude and phase of the annual solar cycle. We observe that estimating the interannual variations with polynomials seems to increase the amplitude of the annual cycle and to create a slight advance of phase. The distribution of the differences in the histograms is tighter than in the previous case since the spatial variability of the difference is less important.

The estimation of the interannual polynomials has opposite consequences on the solar cycle amplitude/phase determination than the estimation of draconitic frequencies, but the differences remain small and near the uncertainty level observed in Figure 7: the majority of stations show differences shorter than a month for the phase and smaller than 1 mm for the amplitude. There are also some isolated and randomly distributed stations with large anomalies (phase or amplitude), meaning that the choice of a particular solution
significantly impacts the results very locally. These stations should be considered as outliers
or studied individually to determine the cause of these differences.

3.3. *Comparison with Hydrological Models and GRACE*

Taking into account that a large part of the solar annual cycle in GNSS time series
in Europe is likely due to hydrological loading, we compared it with loading estimates
computed from hydrological models and GRACE/GRACE Follow-On-derived continental
water storage variations. The site displacements that we computed with the GLDAS2.1,
MERRA2, and GRACE models are expressed in the centre of figure (CF) reference frame for
SNREI Earth. The time series of loading models were adjusted according to Equation (1),
excluding offsets and draconitic harmonics but adding the interannual polynomials. Since
the hydrological models and EOST loading service do not provide uncertainties, we
adjusted the models with the nonweighted least squares method (LS). We show in Figure 9
the differences between the amplitude and phase of the cats_d annual solar cycle with
those of each loading model. Important differences compared with centred Gaussian
distribution are shown by the histograms, especially in the phase shift, where we see that
the centre of the distribution is systematically shifted towards negative phase shifts. In other
words, hydrological loading models are systematically in advance of phase compared to
GNSS. Concerning the amplitudes, the GLDAS2.1 and MERRA2 amplitudes are slightly
larger than those for GNSS, even if the distribution is well centred on zero. For GRACE,
the amplitude seems to be quite underestimated compared to GNSS. Figure 9 highlights
the presence of a coherent spatial pattern, suggesting the existence of disparate common
modes between the GNSS solution and the models. Nevertheless, some stations with high
signal intensity have the same values for every model so that they simply indicate the
difference existing locally between the GNSS solution and the models.

3.4. *Principal Component Analysis of the Interannual Signal*

The interannual signals derived from GNSS and derived from loading models are
compared using PCA decomposition. All the time series were previously detrended (using
the trend estimate from *cats_d* for GNSS) and resampled to nearly every 10 days by taking
the mean of the successive intervals of the year (01/01-10/01, 11/01-20/01, 21/01-31/01,
01/02-10/02, 11/02-20/02, 21/02-28/02, ...). The intervals that contained no data were
left empty. Then, we completed the time series using linear interpolation. We chose to treat
GRACE time series by keeping the original sampling of one month before interpolation. We
finally removed the mean seasonal signal in the same way as [44]. We obtained this cycle by
calculating the mean of the collection of the same dates within years: for example, taking
the mean of every 05/01, the mean of every 15/01, and so on. This is arguably the best way
to filter the mean seasonal cycle compared to a sinusoidal fit. We choose to select stations
only based on the following completeness criterion. The selected stations have available
data between 2010 and 2020, which have accumulated gaps no longer than 60 days, the
largest of which is less than 30 days, and have 90% completeness such that the interpolation
should not distort the signal too much. Implementing these criteria left 268 stations from
the initial network presented in Figure 3. We call EOFs (empirical orthogonal functions)
the spatial function and PCs (principal components) the temporal associated time series of
the PCA output. The first three principal components (also called modes) of the residual
time series are presented in Figure 10, where we report the variance fraction of the total
variance associated with each mode in each panel. We see that the variance fraction of
the GNSS is equal to or lower than the models for each mode and that the GNSS’s first
PCs are much noisier than the ones from the model, even after resampling the data (which
removes a part of the GNSS white noise). The first mode is spatially quite homogeneous.
The associated PCs are very different between GNSS and the loading models. Among the
three models, we can find some common temporal patterns, but also some differences that
can be important, especially with GRACE. The second and third modes are under 15% of
the variance fraction but show common spatial patterns between the solutions. It seems
that the two modes have to be interpreted simultaneously since PCs and EOFs belonging to both appear to be very similar. For example, the PCs and EOFs of panels (i), (f), (g), and (l) seem to be consistent when they belong to two different modes (the same is the case for (j), (k), and (h)).

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Differences in phase and amplitude of annual solar cycle and the associated histograms as described in Figure 8 between *cats_d* and GLDAS2.1 (a–d), MERRA2 (e–h), and GRACE (i–l).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** PCA (EOFs and associated PCs) of GNSS and the three loading models' (including GRACE) residual time series between 2010 and 2020, where the trend and the seasonal signal have been removed: (a–d) first component, (e–h) second component, and (i–l) third component. The percentage of the total variance corresponding to each mode is given in each plot. The PCs are scaled to unit variance and plotted for the period 2010 to 2020, while the corresponding EOFs are given in terms of the correlation between the initial time series and the PCs.</figcaption></figure>
### 3.5. Frequency Content and Interannual Variations

The frequency content of the detrended and 10-day-resampled loading models, along with the GNSS MG3 solution, is provided in Figure 11. The mean Lomb–Scargle [71,72] periodograms showing the amount of variance per frequency band were computed with Python using the algorithm described in [73]. Before stacking, the individual periodograms were un-normalized from the length of the time series in order to compare all the periodograms in a consistent manner. Moreover, the frequency range we used was the same for every time series (GNSS and models). The GNSS has a higher noise level than the loading models (the GNSS periodogram shows a larger amount of variance than the loading models), which is quite understandable if we consider the multiple sources of noise for the GNSS techniques that were listed in the Introduction section of this article. The estimation of the spectral index delivered by CATS for *cats_d* and *cats* is reported in Table 3. In *cats_d*, the spectral index seems to be overestimated when WH noise is taken into account and estimated. Moreover, since for PL noise only, *cats_d* and *cats* provide quite different results, we conclude that the index estimation is also greatly affected by the choice of the estimation model (namely including or not the draconitic cycle). As the spectral indices estimated by CATS impact the uncertainties calculation, it is important to confirm the order of magnitude of these values with the profiles of periodograms of Figure 11. First, the GNSS periodogram is well described by simple PL noise at any frequencies (a slope of constant value). This means that the PL noise contribution dominates the WH noise contribution in the plotted range of frequencies. This can be the reason why CATS fails to estimate a WH noise component in the weekly series while it is robust in estimating only PL noise. The periodograms of loading models seem to be closer to a Gauss–Markov process with a very strong annual signal. They also seem to have the same global behaviour as GNSS (especially the slope) in the interannual band, even if they have a lower variance signal. However, for frequencies greater than 1 cpy, the slopes of periodograms correspond to spectral indexes around −2 for hydrological models, while this remains around −0.7 for GNSS. The origin of this change of slope remains unknown, but we can argue that if the models do not contain WH noise, then we only see the coloured noise dominating at high frequencies. As MERRA2 and GLDAS2.1 had initial samples shorter than a day (1 and 3 hours, respectively), we further investigated to search for WH noise in the model time series. We computed the stacked periodograms of the raw model time series in order to reach frequencies around 300 cpy, where there was still no WH noise, which corroborated our previous assumption.

<table><caption>**Table 3.** PL noise mean spectral index for each of the two CATS-estimated solutions with or without WH noise estimation along PL noise.</caption><thead><tr><th></th><th>East</th><th>North</th><th>Up</th></tr></thead><tbody><tr><td><i>cats_d</i> (WH + PL)</td><td>-0.89</td><td>-1.03</td><td>-0.77</td></tr><tr><td><i>cats_d</i> (PL)</td><td>-0.69</td><td>-0.84</td><td>-0.63</td></tr><tr><td><i>cats</i> (PL)</td><td>-0.86</td><td>-0.99</td><td>-0.75</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Lomb–Scargle periodograms for the three loading models and the detrended MG3 GNSS solution. The periodograms of all stations were stacked and then divided by the number of stations in order to obtain mean periodograms for each solution. The vertical black lines indicate the solar cycle harmonics, and the vertical green lines indicate the draconitic harmonics.</figcaption></figure>

# 4. Discussion

## 4.1. *Interannual Signal in GNSS Time Series*

Concerning the difference between the deterministic (*tiasd*) and stochastic (*cats_d*) accounting of interannual variations, we can point out some key elements. Even if our statistical method using the dispersion of instantaneous velocity provides more realistic uncertainties (more than 10 times larger) than that directly given by the WLS method, we also observe in Figure 6 that the uncertainties for *tiasd* are three to four times smaller than those for the WSL + MLE method, which contain a complete stochastic part (WH + PL noise). Moreover, these small uncertainties are accompanied by important differences in
the velocity values compared to that estimated by CATS: the differences are nested in an interval of ±1.5 mm/y, which is quite important considering the precision requirement on the terrestrial reference frame (such as ITRF) of 0.1 mm/y [12].

Even if it is mathematically correct to fit the time series with polynomials, this is clearly not suitable for geophysical interpretation. Interannual signals in GNSS time series are due to geophysical deformations in only few special cases: regions with melting ice caps or tectonic activity such as slow slip events [4,74–76]. If some GNSS interannual variations have been related to a geophysical origin [44], several previous studies demonstrate that there is a limitation in interpreting this type of signal in terms of geophysics at large spatial scales [77]. Interannual variations in GNSS time series are properly captured by a spatially consistent PL noise model with a relatively large amplitude, regardless of the GNSS station location or the geophysical phenomena affecting them. Therefore, to avoid a misleading interpretation of the interannual signal as a deterministic signal, we do not recommend the use of polynomials in GNSS time series adjustment models. It still remains very difficult to distinguish the true geophysical signal from the noise contribution in the residual time series.

## 4.2. Importance of Draconitic Adjustment

Although it is very common to see GNSS time series adjusted with only the trend and solar cycle, we emphasise here the importance of additionally fitting draconitic harmonics. The differences can be larger than the uncertainty level. The horizontal velocity uncertainties for a large number of stations almost doubled when the draconitic harmonics are not adjusted (Figure 6 and [15]). Even if the numerical estimations of velocity are close in *cats_d* and *cats*, estimating the draconitic periods, or not, could statistically affect the estimated velocity. The vertical component is less affected since the relative power of the draconitic oscillations is lower than that for the horizontal component, as can be seen in the periodograms in Figure 11. If we look at the annual cycle determination, the draconitic signals influence not only the uncertainties but also the values of the phase and amplitude. Even though these differences and the parameter uncertainties are of the same order of magnitude for the majority of stations, there are several stations for which the annual cycle is strongly affected by the draconitic adjustment. In fact, the separation between the solar annual and first draconitic frequencies in the estimation process is conditioned by the length of the time series. The minimum theoretical length for good separation is 25 years. Even though some of the GNSS time series are close to this duration criterion, we are still not able to properly separate the two components for the majority of stations. Moreover, there is actually no evidence that both solar and draconitic signals are stationary. There are potential amplitude variations over time, especially because of their relation to environmental changes, orbit calculation, and the contribution of local effects such as multipath, which is closely related to the time-variable antenna environment and observation geometry. A clean separation of the two signals could then be even more difficult, even for the longest records. To evaluate the correlation between these two terms, we plot the correlation coefficient values depending on the length of the time series, shown in Figure 12. We represent only the correlation coefficient for the annual/draconitic cosine and sine terms of *tiasd* that we extracted from PYTHON's WLS estimation function. We chose to only represent the correlation coefficient for the vertical component of GNSS since we previously verified that it was extremely similar for the East and North components.

The correlation between the solar and draconitic cycles increases when the time series length increases from 4 to 10 y. For longer time spans, the correlation decreases slowly and linearly.

The correlation coefficients of the low time span stations (<10 y) are impacted by a large dispersion. The low values of these coefficients are then relatively not significant and should not be misinterpreted. As we suggested before, the parameters of the longest time span stations (19 y) are still correlated (around 0.4). Extrapolating a decorrelation slope (−0.02 y⁻¹) for longer time spans, we found that the total decorrelation between the two
cycles should happen beyond a 25-year time span (the theoretical limit), which confirms
the difficulty of separating these signals, as they can change with time. In conclusion,
even though the correlation between the solar annual and the first draconitic term is quite
important, it is strongly recommended to include draconitic frequencies in the fit model
in order to reduce the uncertainties, to avoid a joint beat frequency (which can create
the illusion of no stationary annual amplitude) being unmodelled, and to try not to mix
geophysical signals (mostly at solar frequency) and signals coming from orbital errors.

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Coefficient of correlation between annual solar wave and draconitic first harmonic for cosine and sine terms of *tiasd* as a function of the length of the time series.</figcaption></figure>

4.3. *Model Phase Advance over GNSS Seasonal Signal*

The phase advance of the loading models over GNSS displacement could be linked to
a single shortcoming of most global hydrology models: the misrepresentation of horizontal
fluxes. In general, vertical fluxes are well modelled within each individual cell of the
model, but any horizontal runoff often immediately disappears in the oceans. In reality,
this water is still stored for a certain time over land and flows through rivers. Examples
of the importance of the surface water runoff when computing hydrological loading can
be found in [78,79]. Figure 9 shows that the spatial distribution of this phase advance is
not compatible with such a hypothesis because there are no long rivers in Great Britain,
but it can play a role in other regions of the world (the Amazon basin, for example, [78])
and has to be carefully investigated. Globally, there is a good match between the models
and the GNSS, as the differences between them for the majority of stations are of the same
order of magnitude as the GNSS uncertainties (Figure 7). Nevertheless, the differences in
phase values confirm that the GNSS annual cycle actually contains signals other than just
those from hydrology. If we consider that the agreement between GNSS and the models
is good when the histogram of the difference is close to a centred Gaussian distribution
(which should indicate random errors), the agreement between GNSS and GLDAS2.1
seems worse than the agreement between GNSS and MERRA2, which seems worse than
the agreement between GNSS and GRACE. This result shows the importance of dedicated
gravity missions, such as GRACE and GRACE-FO, in helping to improve the modelling of
continental water storage variations. The stations where the difference is similar for each
loading model are stations where the GNSS annual term contains phenomena other than
just the hydrological signature. For example, we did not take into account the nontidal
ocean and atmospheric loading [9,80,81] in the loading models that we used here. Indeed,
the associated annual cycle is quite small in Europe (~0.5 mm) but could explain some
local differences that we see in Figure 9, mostly along the coast. Given the amplitude of the differences, they could also be due to the nature of the ground (karst aquifers [26] and mining [30,31]), but also to the thermal deformation of the surface and antenna monuments [27,28]. A meticulous study of each station time series should provide initial intuition into the sources of these differences, but this goes far beyond the goal of this study. In any case, we can observe that groups of stations at large spatial scales that have similar differences are expected to be more likely affected by a geophysical signal, while anomalies on isolated stations are expected to be more likely due to monumentation deformation or a multipath effect. Unlike the thin peaks corresponding to the annual solar frequency in the loading model periodograms of Figure 11, the large peaks centred on the annual solar cycle in the GNSS periodograms are another important piece of evidence for the multiplicity of annual signal content in GNSS.

### 4.4. Common Mode Estimation in GNSS

Since Europe has had weak tectonic activity and a stable climate over the years without great meteorological events such as El Nino or ice melting, we expected a quite low interannual signal. The PCA decomposition of GNSS in Figure 10a confirms this hypothesis since the PCs seem to be dominated by noise. Loading models have a much lower noise level than the GNSS observations and exhibit significant differences due to different estimates of continental water storage variations. By analysing the second and third modes together in Figure 10, we can find similarities between the EOFs and PCs of (i), (f), (g), and (l), on the one hand, and (j), (k), and (h), on the other hand. However, considering the discrepancies between the models themselves, it seems inadequate to indicate the superiority of any model over another. The choice of one model rather than another in the comparison with GNSS should then be particularly justified, especially if the interpretation is only based on PCA results. The difference in the total variance fraction of the first mode between GNSS and the loading models is most likely due to the difference in noise content (Figure 11).

We qualitatively compared our PCA results with those of [44] who performed PCA of the UNR/NGL (University of Nevada Reno/Nevada Geodetic Laboratory) GNSS solutions [82] over the same region and the same time span as in our study and removed the seasonal signal with the same method as the one presented in Section 3.4. We note that the second mode of MG3 GNSS (Figure 10e) does not appear in [44], and we can speculate about the origin of such a uniform signal over Great Britain associated with this singular time signature. This mode could be associated with the correlated noise of the MG3 time series. As this noise component could differ from one GNSS analysis centre to another, it could be very different for NGL products.

Moreover, we would like to emphasise the importance of the choice of the GNSS network for performing PCA [40]. The choice of a homogeneous network is wise in order to equally distribute the signal across the entire region. However, the stations selected for the PCA also need to meet a completeness criterion during a time interval, which, most of the time, results in an inhomogeneous network. There are then two options. The first is to keep this inhomogeneous network unchanged, knowing that PCA could over-represent the regions with a denser station distribution. This could be the reason why the modes that had the largest variance fraction in our study corresponded to an important signal over the most dense regions (north of Spain, Great Britain, and France). The second option is to select specific stations in order to produce a homogeneous network, asking for subjective criterion selection. For example, what would be the criterion for choosing between two stations, both reaching the completeness criterion, being 50 km away but showing different time series? The choice of particular stations for producing the homogeneous network should thus also impact the PCA results. In conclusion, as the choice of the network has such an impact on the PCA results, and as the GNSS noise could be specific to each GNSS solution, the discrepancies or similarities of the interannual signal given by PCA between models and GNSSs should then be considered and interpreted very carefully.
**Author Contributions:** This study was designed by all the authors. F.P. and S.L. were responsible for the GNSS MG3 products and GINS software. J.-P.B. provided the loading models computed from hydrology models and GRACE. A.M. computed the GNSS solutions. A.S.-G. ran the CATS treatments. A.M. and A.S.-G. analysed the results. A.M. wrote the preliminary draft, and A.M. and S.L. created the figures. All the authors contributed to the original draft preparation. All authors have read and agreed to the published version of the manuscript.

**Funding:** Partial funding for this study was obtained from CNES (Centre National d'Etudes Spatiales).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the
corresponding author.

**Acknowledgments:** We thank all the different agencies, networks, etc., for providing the daily RINEX GNSS data. All the GPS measurements used in this study were obtained from IGN/RGP (Institut Géographique National/Réseau GNSS permanent, http://rgp.ign.fr, (accessed on 7 November 2021)), RENAG (Réseau National GNSS, http://webrenag.unice.fr, (accessed on 7 November 2021)), SONEL (Système d'Observation des variations du Niveau de la mÉr à Long terme, https://www.sonel.org, (accessed on 7 November 2021)), EUREF (EUREF Permanent GNSS Network, https://www.epncb.oma.be, (accessed on 7 November 2021)), RING (Rete Integrata. Nazionale GPS, http://ring.gm.ingv.it, (accessed on 7 November 2021)), BIGF (British Isles continuous GNSS Facility, http://www.bigf.ac.uk, (accessed on 7 November 2021)), CDDIS (Crustal Dynamics Data Information System, https://cddis.nasa.gov, (accessed on 7 November 2021)), NOA (National Observatory of Athens, http://www.gein.noa.gr, (accessed on 7 November 2021)), NIEP (National Institute for Earth Physics, http://gps.infp.ro, (accessed on 7 November 2021)), DGPA (Dutch Permanent GNSS Array, http://gnss1.tudelft.nl/dpga/, (accessed on 7 November 2021)), ERGNSS (Red Geodésica Nacional de Estaciones de Referencia GNSS, https://www.ign.es/, (accessed on 7 November 2021)), and ITACYL (Red de estaciones GNSS de Castilla y León, http://gnss.itacyl.es, (accessed on 7 November 2021)). We also thank anonymous reviewers whose comments and suggestions helped improve and clarify this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Nocquet, J.M. Present-day kinematics of the Mediterranean: A comprehensive overview of GPS results. *Tectonophysics* **2012**, *579*, 220–242. [CrossRef]

2. Masson, C.; Mazzotti, S.; Vernant, P.; Doerflinger, E. Extracting small deformation beyond individual station precision from dense Global Navigation Satellite System (GNSS) networks in France and western Europe. *Solid Earth* **2019**, *10*, 1905–1920. [CrossRef]

3. Kouba, J. Measuring seismic waves induced by large earthquakes with GPS. *Stud. Geophys. Geod.* **2003**, *47*, 741–755. [CrossRef]

4. Klein, E.; Fleitout, L.; Vigny, C.; Garaud, J. Afterslip and viscoelastic relaxation model inferred from the large-scale post-seismic deformation following the 2010 M w 8.8 Maule earthquake (Chile). *Geophys. J. Int.* **2016**, *205*, 1455–1472. [CrossRef]

5. van Dam, T.; Wahr, J.; Milly, P.; Shmakin, A.; Blewitt, G.; Lavallée, D.; Larson, K. Crustal displacements due to continental water loading. *Geophys. Res. Lett.* **2001**, *28*, 651–654. [CrossRef]

6. Argus, D.F.; Fu, Y.; Landerer, F.W. Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion. *Geophys. Res. Lett.* **2014**, *41*, 1971–1980. [CrossRef]

7. Chanard, K.; Fleitout, L.; Calais, E.; Rebischung, P.; Avouac, J.P. Toward a global horizontal and vertical elastic load deformation model derived from GRACE and GNSS station position time series. *J. Geophys. Res. Solid Earth* **2018**, *123*, 3225–3237. [CrossRef]

8. Martens, H.R.; Argus, D.F.; Norberg, C.; Blewitt, G.; Herring, T.A.; Moore, A.W.; Hammond, W.C.; Kreemer, C. Atmospheric pressure loading in GPS positions: Dependency on GPS processing methods and effect on assessment of seasonal deformation in the contiguous USA and Alaska. *J. Geod.* **2020**, *94*, 1–22. [CrossRef]

9. Mémin, A.; Boy, J.P.; Santamaria-Gomez, A. Correcting GPS measurements for non-tidal loading. *GPS Solut.* **2020**, *24*, 45. [CrossRef]

10. Johansson, J.M.; Davis, J.L.; Scherneck, H.G.; Milne, G.A.; Vermeer, M.; Mitrovica, J.X.; Bennett, R.A.; Jonsson, B.; Elgered, G.; Elósegui, P.; et al. Continuous GPS measurements of postglacial adjustment in Fennoscandia 1. Geodetic results. *J. Geophys. Res. Solid Earth* **2002**, *107*, ETG-3. [CrossRef]

11. Grapenthin, R.; Sigmundsson, F.; Geirsson, H.; Arnadóttir, T.; Pinel, V. Icelandic rhythmics: Annual modulation of land elevation and plate spreading by snow load. *Geophys. Res. Lett.* **2006**, *33*. [CrossRef]

12. Altamimi, Z.; Rebischung, P.; Métivier, L.; Collilieux, X. ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions. *J. Geophys. Res. Solid Earth* **2016**, *121*, 6109–6131. [CrossRef]
13. Mao, A.; Harrison, C.G.; Dixon, T.H. Noise in GPS coordinate time series. *J. Geophys. Res. Solid Earth* **1999**, *104*, 2797–2816.
[CrossRef]

14. Williams, S. The effect of coloured noise on the uncertainties of rates estimated from geodetic time series. *J. Geod.* **2003**, *76*, 483–494.
[CrossRef]

15. Santamaría-Gómez, A.; Bouin, M.N.; Collilieux, X.; Wöppelmann, G. Correlated errors in GPS position time series: Implications for velocity estimates. *J. Geophys. Res. Solid Earth* **2011**, *116*. [CrossRef]

16. Klos, A.; Olivares, G.; Teferle, F.N.; Hunegnaw, A.; Bogusz, J. On the combined effect of periodic signals and colored noise on velocity uncertainties. *GPS Solut.* **2018**, *22*, 1. [CrossRef]

17. Dong, D.; Fang, P.; Bock, Y.; Cheng, M.; Miyazaki, S. Anatomy of apparent seasonal variations from GPS-derived site position time series. *J. Geophys. Res. Solid Earth* **2002**, *107*, ETG-9. [CrossRef]

18. Griffiths, J.; Ray, J.R. On the precision and accuracy of IGS orbits. *J. Geod.* **2009**, *83*, 277–287. [CrossRef]

19. Tregoning, P.; Watson, C. Atmospheric effects and spurious signals in GPS analyses. *J. Geophys. Res. Solid Earth* **2009**, *114*.
[CrossRef]

20. Johnson, H.O.; Agnew, D.C. Monument motion and measurements of crustal velocities. *Geophys. Res. Lett.* **1995**, *22*, 2905–2908.
[CrossRef]

21. King, M.A.; Williams, S.D. Apparent stability of GPS monumentation from short-baseline time series. *J. Geophys. Res. Solid Earth* **2009**, *114*. [CrossRef]

22. King, M.A.; Watson, C.S. Long GPS coordinate time series: Multipath and geometry effects. *J. Geophys. Res. Solid Earth* **2010**, *115*.
[CrossRef]

23. Santamaría-Gómez, A.; Ray, J. Chameleonic noise in GPS position time series. *J. Geophys. Res. Solid Earth* **2021**, *126*, e2020JB019541.
[CrossRef]

24. Williams, S.D.; Bock, Y.; Fang, P.; Jamason, P.; Nikolaidis, R.M.; Prawirodirdjo, L.; Miller, M.; Johnson, D.J. Error analysis of continuous GPS position time series. *J. Geophys. Res. Solid Earth* **2004**, *109*. [CrossRef]

25. Amiri-Simkooei, A. On the nature of GPS draconitic year periodic pattern in multivariate position time series. *J. Geophys. Res. Solid Earth* **2013**, *118*, 2500–2511. [CrossRef]

26. Silverii, F.; d'Agostino, N.; Métois, M.; Fiorillo, F.; Ventafridda, G. Transient deformation of karst aquifers due to seasonal and multiyear groundwater variations observed by GPS in southern Apennines (Italy). *J. Geophys. Res. Solid Earth* **2016**, *121*, 8315–8337.
[CrossRef]

27. Fang, M.; Dong, D.; Hager, B.H. Displacements due to surface temperature variation on a uniform elastic sphere with its centre of mass stationary. *Geophys. J. Int.* **2014**, *196*, 194–203. [CrossRef]

28. Xu, X.; Dong, D.; Fang, M.; Zhou, Y.; Wei, N.; Zhou, F. Contributions of thermoelastic deformation to seasonal variations in GPS station position. *GPS Solut.* **2017**, *21*, 1265–1274. [CrossRef]

29. Argus, D.F.; Landerer, F.W.; Wiese, D.N.; Martens, H.R.; Fu, Y.; Famiglietti, J.S.; Thomas, B.F.; Farr, T.G.; Moore, A.W.; Watkins, M.M. Sustained water loss in California's mountain ranges during severe drought from 2012 to 2015 inferred from GPS. *J. Geophys. Res. Solid Earth* **2017**, *122*, 10–559. [CrossRef]

30. Gourmelen, N.; Amelung, F.; Casu, F.; Manzo, M.; Lanari, R. Mining-related ground deformation in Crescent Valley, Nevada: Implications for sparse GPS networks. *Geophys. Res. Lett.* **2007**, *34*. [CrossRef]

31. Muntean, A.; Mocanu, V.; Ambrosius, B. A GPS study of land subsidence in the Petrosani (Romania) coal mining area. *Nat. Hazards* **2016**, *80*, 797–810. [CrossRef]

32. Ray, J.; Altamimi, Z.; Collilieux, X.; van Dam, T. Anomalous harmonics in the spectra of GPS position estimates. *GPS Solut.* **2008**, *12*, 55–64. [CrossRef]

33. Guo, S.; Shi, C.; Wei, N.; Li, M.; Fan, L.; Zhang, D. Effect of ambiguity resolution on the draconitic errors in sub-daily GPS position estimates. *GPS Solut.* **2021**, *25*, 101. [CrossRef]

34. Loyer, S.; Perosanz, F.; Mercier, F.; Capdeville, H.; Marty, J.C. Zero-difference GPS ambiguity resolution at CNES-CLS IGS Analysis Center. *J. Geod.* **2012**, *86*, 991–1003. [CrossRef]

35. Klos, A.; Bogusz, J.; Bos, M.S.; Gruszczynska, M. Modelling the GNSS Time Series: Different Approaches to Extract Seasonal. In *Geodetic Time Series Analysis in Earth Sciences*; Springer: Cham, Switzerland, 2019; p. 211.

36. Williams, S.D. CATS: GPS coordinate time series analysis software. *GPS Solut.* **2008**, *12*, 147–153. [CrossRef]

37. He, X.; Montillet, J.P.; Fernandes, R.; Bos, M.; Yu, K.; Hua, X.; Jiang, W. Review of current GPS methodologies for producing accurate time series and their error sources. *J. Geodyn.* **2017**, *106*, 12–29. [CrossRef]

38. Serpelloni, E.; Faccenna, C.; Spada, G.; Dong, D.; Williams, S.D. Vertical GPS ground motion rates in the Euro-Mediterranean region: New evidence of velocity gradients at different spatial scales along the Nubia-Eurasia plate boundary. *J. Geophys. Res. Solid Earth* **2013**, *118*, 6003–6024. [CrossRef]

39. Shen, Y.; Li, W.; Xu, G.; Li, B. Spatiotemporal filtering of regional GNSS network's position time series with missing data using principle component analysis. *J. Geod.* **2014**, *88*, 1–12. [CrossRef]

40. Wu, S.; Nie, G.; Liu, J.; Wang, K.; Xue, C.; Wang, J.; Li, H.; Peng, F.; Ren, X. A sub-regional extraction method of common mode components from IGS and CMONOC stations in China. *Remote Sens.* **2019**, *11*, 1389. [CrossRef]

41. Liu, B.; Dai, W.; Peng, W.; Meng, X. Spatiotemporal analysis of GPS time series in vertical direction using independent component analysis. *Earth Planets Space* **2015**, *67*, 189. [CrossRef]
42. Kreemer, C.; Blewitt, G. Robust estimation of spatially varying common-mode components in GPS time-series. *J. Geod.* **2021**, *95*, 13. [CrossRef]

43. Tapley, B.D.; Watkins, M.M.; Flechtner, F.; Reigber, C.; Bettadpur, S.; Rodell, M.; Sasgen, I.; Famiglietti, J.S.; Landerer, F.W.; Chambers, D.P.; et al. Contributions of GRACE to understanding climate change. *Nat. Clim. Chang.* **2019**, *9*, 358–369. [CrossRef] [PubMed]

44. Elia, L.; Zerbini, S.; Raicich, F. Interannual Variability of GPS Heights and Environmental Parameters over Europe and the Mediterranean Area. *Remote Sens.* **2021**, *13*, 1554. [CrossRef]

45. Katsigianni, G.; Perosanz, F.; Loyer, S.; Gupta, M. Galileo millimeter-level kinematic precise point positioning with ambiguity resolution. *Earth Planets Space* **2019**, *71*, 76. [CrossRef]

46. Rebischung, P. *Terrestrial Frame Solutions from the IGS Third Reprocessing*. In *Proceedings of the EGU General Assembly Conference Abstracts*, online, 19–30 April 2021; p. EGU21-2144.

47. Lemoine, J.M.; Biancale, R.; Reinquin, F.; Bourgogne, S.; Gégout, P. *CNES/GRGS RL04 Earth Gravity Field Models, from GRACE and SLR Data*; GFZ Data Services. 2019. Available online: https://doi.org/10.5880/ICGEM.2019.010 (accessed on 1 October 2021).

48. Lyard, F.H.; Allain, D.J.; Cancet, M.; Carrère, L.; Picot, N. FES2014 global ocean tide atlas: Design and performance. *Ocean Sci.* **2021**, *17*, 615–649. [CrossRef]

49. Folkner, W.M.; Williams, J.G.; Boggs, D.H. The planetary and lunar ephemeris DE 421. *IPN Prog. Rep.* **2009**, *42*, 1–34.

50. Villiger, A.; Dach, R. *International GNSS Service: Technical Report 2019*; Technical Report; IGS Central Bureau and University of Bern: Bern, Switzerland, 2020.

51. Petit, G.; Luzum, B. *IERS Conventions (2010)*; Technical Report; Bureau International des Poids et Mesures : Sèvres, France, 2010.

52. Desai, S.D.; Sibois, A.E. Evaluating predicted diurnal and semidiurnal tidal variations in polar motion with GPS-based observations. *J. Geophys. Res. Solid Earth* **2016**, *121*, 5237–5256. [CrossRef]

53. Ray, R.D.; Ponte, R.M. Barometric tides from ECMWF operational analyses. *Ann. Geophys.* **2003**, *21*, 1897–1910. [CrossRef]

54. Boehm, J.; Niell, A.; Tregoning, P.; Schuh, H. Global Mapping Function (GMF): A new empirical mapping function based on numerical weather model data. *Geophys. Res. Lett.* **2006**, *33*. [CrossRef]

55. Lagler, K.; Schindelegger, M.; Böhm, J.; Krásná, H.; Nilsson, T. GPT2: Empirical slant delay model for radio space geodetic techniques. *Geophys. Res. Lett.* **2013**, *40*, 1069–1073. [CrossRef]

56. Williams, S.D. Offsets in global positioning system time series. *J. Geophys. Res. Solid Earth* **2003**, *108*. B002156. [CrossRef]

57. Bruni, S.; Zerbini, S.; Raicich, F.; Errico, M.; Santi, E. Detecting discontinuities in GNSS coordinate time series with STARS: case study, the Bologna and Medicina GPS sites. *J. Geod.* **2014**, *88*, 1203–1214. [CrossRef]

58. Santamaría-Gómez, A.; Mémin, A. Geodetic secular velocity errors due to interannual surface loading deformation. *Geophys. J. Int.* **2015**, *202*, 763–767. [CrossRef]

59. Moré, J.J. The Levenberg-Marquardt algorithm: Implementation and theory. In *Numerical Analysis*; Springer: Berlin/Heidelberg, Germany, 1978; pp. 105–116.

60. Collectif Jolidon. *Physique Éxpérimentale*; EDP Sciences: Les Ulis, France, 2021.

61. Rodell, M.; Houser, P.; Jambor, U.; Gottschalck, J.; Mitchell, K.; Meng, C.J.; Arsenault, K.; Cosgrove, B.; Radakovich, J.; Bosilovich, M.; et al. The global land data assimilation system. *Bull. Am. Meteorol. Soc.* **2004**, *85*, 381–394. [CrossRef]

62. Gelaro, R.; McCarty, W.; Suárez, M.J.; Todling, R.; Molod, A.; Takacs, L.; Randles, C.A.; Darmenov, A.; Bosilovich, M.G.; Reichle, R.; et al. The modern-era retrospective analysis for research and applications, version 2 (MERRA-2). *J. Clim.* **2017**, *30*, 5419–5454. [CrossRef]

63. Loomis, B.; Luthcke, S.; Sabaka, T. Regularization and error characterization of GRACE mascons. *J. Geod.* **2019**, *93*, 1381–1398. [CrossRef]

64. Farrell, W. Deformation of the Earth by surface loads. *Rev. Geophys.* **1972**, *10*, 761–797. [CrossRef]

65. Dziewonski, A.M.; Anderson, D.L. Preliminary reference Earth model. *Phys. Earth Planet. Inter.* **1981**, *25*, 297–356. [CrossRef]

66. Petrov, L.; Boy, J.P. Study of the atmospheric pressure loading signal in very long baseline interferometry observations. *J. Geophys. Res. Solid Earth* **2004**, *109*. [CrossRef]

67. Altamimi, Z.; Métivier, L.; Rebischung, P.; Rouby, H.; Collilieux, X. ITRF2014 plate motion model. *Geophys. J. Int.* **2017**, *209*, 1906–1912. [CrossRef]

68. Peltier, W.R.; Argus, D.; Drummond, R. Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C (VM5a) model. *J. Geophys. Res. Solid Earth* **2015**, *120*, 450–487. [CrossRef]

69. Peltier, R.W.; Argus, D.F.; Drummond, R. Comment on “An assessment of the ICE-6G_C (VM5a) glacial isostatic adjustment model” by Purcell et al. *J. Geophys. Res. Solid Earth* **2018**, *123*, 2019–2028. [CrossRef]

70. van Dam, T.; Wahr, J.; Lavallée, D. A comparison of annual vertical crustal displacements from GPS and Gravity Recovery and Climate Experiment (GRACE) over Europe. *J. Geophys. Res. Solid Earth* **2007**, *112*. [CrossRef]

71. Lomb, N.R. Least-squares frequency analysis of unequally spaced data. *Astrophys. Space Sci.* **1976**, *39*, 447–462. [CrossRef]

72. Scargle, J.D. Studies in astronomical time series analysis. II-Statistical aspects of spectral analysis of unevenly spaced data. *Astrophys. J.* **1982**, *263*, 835–853. [CrossRef]

73. Townsend, R. Fast calculation of the Lomb-Scargle periodogram using graphics processing units. *Astrophys. J. Suppl. Ser.* **2010**, *191*, 247. [CrossRef]
74. Laxon, S.; Peacock, N.; Smith, D. High interannual variability of sea ice thickness in the Arctic region. *Nature* **2003**, *425*, 947–950.
[CrossRef]

75. Klein, E.; Duputel, Z.; Zigone, D.; Vigny, C.; Boy, J.P.; Doubre, C.; Meneses, G. Deep transient slow slip detected by survey GPS in the region of Atacama, Chile. *Geophys. Res. Lett.* **2018**, *45*, 12–263. [CrossRef]

76. Tobita, M. Combined logarithmic and exponential function model for fitting postseismic GNSS time series after 2011 Tohoku-Oki earthquake. *Earth Planets Space* **2016**, *68*, 41. [CrossRef]

77. Rosat, S.; Gillet, N.; Boy, J.P.; Couhert, A.; Dumberry, M. Interannual variations of degree 2 from geodetic observations and surface processes. *Geophys. J. Int.* **2021**, *225*, 200–221.

78. Nicolas, J.; Verdun, J.; Boy, J.P.; Bonhomme, L.; Asri, A.; Corbeau, A.; Berthier, A.; Durand, F.; Clarke, P. Improved Hydrological Loading Models in South America: Analysis of GPS Displacements Using M-SSA. *Remote Sens.* **2021**, *13*, 1605. [CrossRef]

79. Nahmani, S.; Bock, O.; Bouin, M.N.; Santamaría-Gómez, A.; Boy, J.P.; Collilieux, X.; Métivier, L.; Panet, I.; Genthon, P.; de Linage, C.; et al. Hydrological deformation induced by the West African Monsoon: Comparison of GPS, GRACE and loading models. *J. Geophys. Res. Solid Earth* **2012**, *117*. [CrossRef]

80. Williams, S.; Penna, N. Non-tidal ocean loading effects on geodetic GPS heights. *Geophys. Res. Lett.* **2011**, *38*. [CrossRef]

81. van Dam, T.; Collilieux, X.; Wuite, J.; Altamimi, Z.; Ray, J. Nontidal ocean loading: Amplitudes and potential effects in GPS height time series. *J. Geod.* **2012**, *86*, 1043–1057. [CrossRef]

82. Blewitt, G.; Hammond, W.C.; Kreemer, C. Harnessing the GPS data explosion for interdisciplinary science. *EOS* **2018**, *99*, 485.
[CrossRef]