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

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.

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,6,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,15,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 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\mathrm{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\mathrm{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,39,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. Dynamical and loading models used for MG3 products.

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

, 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/wwww/repro3/GRG6RE3FIN*.gz (accessed on 7 November 2021).

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.

2.2.2. Parameter Estimation

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

$$\begin{array}{cc}\hfill y\left(t\right)& =y_0+v(t-t_0)\hfill \\ \hfill & +\sum _{m=1}^2{a}_{m}cos\left(m{\omega }_{\mathrm{a}}t\right)+b_{m}sin\left(m{\omega }_{\mathrm{a}}t\right)\hfill \\ \hfill & +\sum _{n=1}^8{c}_{n}cos\left(n{\omega }_{\mathrm{d}}t\right)+d_{n}sin\left(n{\omega }_{\mathrm{d}}t\right)\hfill \\ \hfill & +\sum _{k=1}^n{s}_k\mathcal{H}(t-t_k)+\epsilon \left(t\right),\hfill \end{array}$$

(1)

where $y_0$ and v are, respectively, the intercept (at epoch $t_0$) and the linear velocity of the station; ${\omega }_{\mathrm{a}}/2\pi =1\mathrm{cpy}$ is the solar frequency; and ${\omega }_{\mathrm{d}}/2\pi =1.04\mathrm{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 $\mathcal{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 $\epsilon \left(t\right)$ is the stochastic part of the model that has to be modelled and then estimated.

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 $\epsilon \left(t\right)$ 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\mathrm{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. Estimation properties of the three models considered in this article.

	cats_d	cats	tiasd
Draconitic frequencies adjustment	Yes	No	Yes
Interannual variation	–	–	polynomials
Stochastic model	WH + PL noise	WH + PL noise	WH noise

. 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.

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\mathrm{mm}/\mathrm{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] (${\mathrm{lon}}_{\mathrm{p}}=-99.094\left(7\right)°$ E, ${\mathrm{lat}}_{\mathrm{p}}=+55.070\left(4\right)°$ N, $\omega =0.261\left(1\right)\times {10}^{-6}\mathrm{deg}/\mathrm{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.

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\mathrm{mm}/\mathrm{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).

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.

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.

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\mathrm{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.

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\mathrm{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)).

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. PL noise mean spectral index for the each of the two CATS-estimated solutions with or without WH noise estimation along PL noise.

	East	North	Up
cats_d (WH + PL)	$-0.89$	$-1.03$	$-0.77$
cats_d (PL)	$-0.69$	$-0.84$	$-0.63$
cats (PL)	$-0.86$	$-0.99$	$-0.75$

. 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\mathrm{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\mathrm{cpy}$, where there was still no WH noise, which corroborated our previous assumption.

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 $\pm 1.5\mathrm{mm}/\mathrm{y}$, which is quite important considering the precision requirement on the terrestrial reference frame (such as ITRF) of $0.1\mathrm{mm}/\mathrm{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,75,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\mathrm{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${}^{-1}$) 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.

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.