# Analysis of GNSS Displacements in Europe and Their Comparison with Hydrological Loading Models
## 2 Materials and Methods
### Time Series Analysis
#### 2.2.2 Parameter Estimation  
For each time series, we determined different parameters corresponding to the classical model used in GNSS position adjustment:  
\[\begin{split} y(t)&=y_{0}+v(t-t_{0})\\ &+\sum_{m=1}^{2}a_{m}\cos(m\omega_{\text{a}}t)+b_{m}\sin(m\omega_{ \text{a}}t)\\ &+\sum_{n=1}^{8}c_{n}\cos(n\omega_{\text{d}}t)+d_{n}\sin(n\omega_{ \text{d}}t)\\ &+\sum_{k=1}^{n}s_{k}\mathcal{H}(t-t_{k})+\varepsilon(t),\end{split} \tag{1}\]  
Figure 1: GNSS station network used for MG3 product computation.  
Figure 2: (**a**) Number of satellites and (**b**) number of stations used in MG3 product computation as a function of time.  
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\,\mathrm{c}\mathrm{p}\mathrm{y}\) is the solar frequency; and \(\omega_{d}/2\pi=1.04\,\mathrm{c}\mathrm{p}\mathrm{y}\) 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 \(\varepsilon(t)\) 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  
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.  
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_.