Identifying the sensitivity of GPS to non-tidal loadings at various time resolutions: examining vertical displacements from continental Eurasia

We examine the sensitivity of the Global Positioning System (GPS) to non-tidal loading for a set of continental Eurasia permanent stations. We utilized daily vertical displacements available from the Nevada Geodetic Laboratory (NGL) at stations located at least 100 km away from the coast. Loading-induced predictions of displacements of earth’s crust are provided by the Earth-System-Modeling Group of the GFZ (ESMGFZ). We demonstrate that the hydrological loading, supported by barystatic sea-level changes to close the global mass budget (HYDL + SLEL), contributes to GPS displacements only in the seasonal band. Non-tidal atmospheric loading, supported by non-tidal oceanic loading (NTAL + NTOL), correlates positively with GPS displacements for almost all time resolutions, including non-seasonal changes from 2 days to 5 months, which are often considered as noise, intra-seasonal and seasonal changes with periods between 4 months and 1.4 years, and, also, inter-annual signals between 1.1 and 3.0 years. Correcting the GPS vertical displacements by NTAL leads to a reduction in the time series variances, evoking a whitening of the GPS stochastic character and a decrease in the standard deviation of noise. Both lead, on average, to an improvement in the uncertainty of the GPS vertical velocity by a factor of 2. To reduce its impact on the GPS displacement time series, we recommend that NTAL is applied at the observation level during the processing of GPS observations. HYDL might be corrected at the observation level or remain in the data and be applied at the stage of time series analysis.


Introduction
Vertical displacements of permanently mounted receivers of Global Navigation Satellite Systems (GNSSs) are induced by a range of processes occurring at various temporal and spatial resolutions. Temporal resolutions from weekly up to interannual or secular changes are well covered using dailysampled GNSS displacements, from which the longest time series available so far span even 25 years. Spatial resolutions from site-dependent phenomena up to regional-or continental-scale changes are resolved, respectively, from single GNSS stations or from distributed networks of stations (Dixon et al. 1990;Park et al. 2004;Kierulf et al. 2014). The plurality of geophysical effects, systematic errors or GNSSrelated signals, and the still unknown GNSS sensitivity itself may raise problems when separating various contributors in GNSS displacements (Chanard et al. 2020). While contributors to the long-term and seasonal signals are widely acknowledged (Argus et al. 2014;Kierulf et al. 2014), the short-term part of the spectra of the GNSS displacements remains regularly unexplained. Much research related these residuals to common-mode errors, GNSS-related errors, or site-specific effects, a combination of which was simply explored under the term "noise" (Bogusz and Kontny 2011;Tian and Shen 2016) and characterized by the power-law noise using spectral indices and amplitudes (He et al. 2017). The majority of stations around the globe are affected by power-law noise close to flicker noise (spectral index equal to − 1), while others experiencing local or random effects are primarily affected by random-walk noise (spectral index equal to − 2). In addition to these two, white noise (spectral index equal to 0) is often used to represent the highest frequencies reliably.
The variance of non-tidal mass variability in geophysical fluids of the earth surface is generated mainly by the annual signal for hydrological loading and by changes outside the seasonal band for both non-tidal atmospheric and oceanic loading models (Gruszczynska et al. 2018). As a consequence, loading-predicted displacements explain a variance of GNSS vertical displacements not only in the seasonal frequency band (Dong et al. 2002;Dill and Dobslaw 2013;Chanard et al. 2018), but also contribute to changes frequently considered as noise (Klos et al. 2020;Mémin et al. 2020). Regarding geographical location, hydrological loading is prominent for large river basins (Moreira et al. 2016;Karegar et al. 2018) or for areas excited by snow accumulation and those close to artificial reservoirs (Neelmeijer et al. 2018). Non-tidal oceanic loading appears to be significant only for the coastal regions and islands, while nontidal atmospheric loading seems to dominate over Eastern Europe and Asia (Brondeel and Willems 2003;Williams and Penna 2011). Williams and Penna (2011) and Bian et al. (2020) showed that removing loading-predicted displacements from the Global Positioning System (GPS) displacement time series leads to an improvement in their standard deviations. The standard deviations improve even more when the same loading effects are applied at the observation level (Dach et al. 2011;Männel et al. 2019). Since standard deviations indicate only a broad-band signal reduction, their improvement cannot be directly related to any frequency band. Mémin et al. (2020) attempted to overcome this problem by filtering the estimated loading-predicted displacements and the vertical GPS displacements into several frequency bands and calculating band-dependent root-mean-square (RMS) reductions once the loading time series is removed. The three following broad frequency bands were considered: the longperiod band above 2 months, containing long-term, seasonal and non-seasonal changes, the intermediate-period band with periods between 1 week and 2 months, and the shortperiod band with periods of less than a week, both containing non-seasonal changes and noise. The authors emphasized that continental hydrological loading agrees with the GPS displacements only in the long-period frequency band, while non-tidal atmospheric and non-tidal oceanic loadings agree with the GPS displacements in both the long-and intermediate-period bands. Short periods remained unexplained. This study takes a step further and provides a division of GPS observed and loading-predicted displacements into narrower frequency bands, clearly distinguishing them by periods. We present a set of correlation coefficients for different period bands and discuss their significance. Our analyses of hydrological loading are supported by barystatic sea-level changes to close the global mass budget, unprovided before. Moreover, we discuss the impact of various non-tidal loading models on the uncertainty of GPS velocity and provide practical implications of our study.

GPS data
GPS vertical displacements were obtained from the Nevada Geodetic Laboratory (NGL, Blewitt et al. 2018). We selected 100 stations located in continental Eurasia with observations spanning 5 years and longer. We restricted ourselves to stations situated at least 100 km away from the coast. The time series were screened to remove outliers and offsets. For the outliers, three times the interquartile range rule was used. To remove offsets, we employed epochs reported within the NGL database but augmented those by manual inspection. GPS displacements are complete to 50-99%. Gaps are not interpolated; they have no impact on our analyses.

Non-tidal loading models
We use globally gridded crustal displacements in the vertical direction in the center-of-mass reference frame provided by the Earth-System-Modeling Group of Deutsches Geo-ForschungsZentrum Sect. 1.3 (ESMGFZ; Dill and Dobslaw 2013). Elastic surface deformations are calculated in the spatial domain by convolving the load Green's function for the elastic earth model "ak135" (Kennett et al. 1995) with mass loads. Station time series are extracted from the gridded loading displacements provided on a regular 0.5° × 0.5° global grid using bicubic interpolation to the GPS location and avoiding interpolations over the coastline.
Non-tidal atmospheric mass loads (NTAL) are derived from the atmospheric surface pressure given by 3-hourly ECMWF (European Centre for Medium-Range Weather Forecasts) reanalysis, ERA-40, ERA-Interim, and operational ECMWF data on a global regular 0.5° grid. To isolate the non-tidal variability of the atmosphere, we removed atmospheric tides considering 12 major tidal constituents: main solar tides S1, S2, and S3, and the main semi-diurnal lunar tide M2, each with two additional side-bands added.
Non-tidal oceanic mass loads (NTOL) are derived from the 3-hourly ocean bottom pressure from the Max Planck Institute Ocean Model MPIOM (Jungclaus et al. 2013) given on a global 1.0° grid. Oceanic tides excited by the timevarying atmospheric surface pressure were considered in a similar way as in NTAL.
Hydrological mass loads (HYDL) are based on terrestrial water storage simulated 24-hourly on a global regular 0.5° grid with the Land Surface Discharge Model (LSDM, Dill 2008). Mass loads are given as soil moisture, snow, surface water, and water in rivers and lakes. Water stored in rivers and lakes was relocated from the 0.5° model grid to a 0.125° georeferenced river map  to capture the exceptional high local loading signals in the proximity of larger rivers. To validate the HYDL-predicted displacements, we used GRACE (Gravity Recovery and Climate Experiment) observations in the form of spherical harmonic coefficients up to degree and order 90, available from GFZ as RL06 solution de-noised with DDK3 filter. Spherical harmonic coefficients were transformed into vertical displacements using load Love numbers.
We also utilize a loading model representing barystatic sea-level changes (SLEL) calculated by solving the sea-level equation for the net-mass deficit of the sum of all masses in the atmosphere, oceans, and the continental hydrosphere . SLEL, therefore, represents a term required to close the global mass budget in line with the reasoning of Chen (2005).
For comparison, we also utilize loading models provided for individual stations by the EOST (Ecole et Observatoire des Sciences de la Terre, Strasbourg, France; Mémin et al. 2020) loading service, where 3-D displacements induced by non-tidal loadings are available for a set of global stations. We use three non-tidal atmospheric loading models, which we further refer to as "ERAIN", "ATMIB", and "ATMMO" and one continental hydrological loading being referred to as "MERRA2", in accordance with the file extensions provided by EOST. The ERAIN model is estimated from the surface pressure from the ERA-Interim (ECMWF Reanalysis) model, assuming an inverted barometer ocean response to pressure loading. The ATMIB and ATMMO models are estimated from the surface pressure from the ECMWF operational model. In the ATMIB model, an inverted barometer ocean response to pressure force is assumed. In the ATMMO model, a dynamic ocean response to pressure and wind is assumed. For this purpose, the TUGO-m barotropic model is employed (Carrere and Lyard 20032003). MERRA2 hydrological loading is estimated from the MERRA2 model (Gelaro et al. 2017). It includes soil moisture and snow hydrology components, excluding areas permanently covered with ice.

Wavelet filtering
Wavelet decomposition enables signals to be captured at various predefined time resolutions. We used a Meyer mother wavelet (Meyer 1990) to divide the loading-predicted and the GPS displacements into 9 decomposition levels ( Table 1). This choice was motivated by previous analyses, where we used the Meyer mother wavelet to effectively model both short-and long-term periods in geodetic observations (Bogusz et al. 2013;Klos et al. 2018a). To overcome the problem of edge effects, we extrapolated the time series of five years, both before the actual observations started and after they ended, thereby extending each time series by an additional ten years. For this extrapolation, we used a leastsquares model derived from the GPS time series, including trend, annual and semi-annual cycles, to which we added a white noise of variance equal to the variance of GPS residuals. Once the decomposition was performed, those ten years were removed again from the time series, leaving just the original length of data. Wavelet decomposition was applied individually to the GPS displacements and the displacements predicted by the geophysical loading models. We note that wavelet decomposition does not allow for a separation between signals and noise, meaning that all changes in the predefined frequency bands are retrieved, i.e., both, signals of interest and noise (Klos et al. 2018a). This may reduce the application of wavelet decomposition for retrieving signals of interest but does not hamper the comparisons of certain decomposition levels, as performed in this study.
We present GPS displacements and loading-predicted displacements for randomly selected GPS station ARTU (Russia), which is situated close to the city of Jekaterynburg (Fig. 1). The GPS displacements are characterized by high data quality only affected by a small period of missing data in 2019. The signal magnitudes are quite comparable between GPS and loading models. A total environmental loading (NTAL + NTOL + HYDL + SLEL) almost perfectly fits the original GPS displacements. NTAL + NTOLpredicted displacements clearly contribute to the high frequency changes present in GPS displacements, while HYDL + SLEL-predicted displacements project seasonal and long-term changes also seen in GPS, with clear nonlinear variations in 2011. Vertical GPS displacements plotted versus loading-predicted displacements point to a linear relationship between both (Fig. 1). We see clear agreement between the variability of GPS displacements and the one estimated for NTAL + NTOL + HYDL + SLEL at individual decomposition levels (Fig. 2). Two exceptions are noticed here: D1 and D8 decomposition levels, for which NTAL + NTOL + HYDL + SLEL is characterized by smaller and larger RMS than GPS, respectively. When looking at the contribution of NTAL + NTOL to different levels of decomposition of NTAL + NTOL + HYDL + SLEL, we note that the entire variability of total loading is generated by NTAL + NTOL for D1-D6 levels of detail (2 days-5 months). HYDL + SLEL contributes to lower frequencies, i.e., details D7-D9 (4 months-3 years) with significant contribution also noticed for inter-annual changes, level A9 (see also Fig. 1). For the smallest levels of detail (D1-D5), almost no signal is induced by HYDL + SLEL. The above HYDL + SLEL contribution is confirmed by GRACE-predicted displacements, which match almost perfectly the HYDL + SLEL RMS values. Since GRACE-predicted displacements are available with monthly sampling intervals, these can only be used to validate decomposition levels from D6 onwards. We further note that the contribution of NTOL and SLEL is hardly observed for inland stations we utilized; discussed further. Therefore, we present the sum of HYDL + SLEL and NTAL + NTOL.

Atmospheric loading signals in GPS
We now present the Pearson correlation coefficients between the vertical GPS displacements and NTAL loading-predicted displacements for different levels of decomposition (Fig. 3). We note (weak) positive correlations for almost all of the 100 stations considered in our study, even for the shortest periods between 2 and 10 days (D1-D2). The correlations are close to zero for particular stations with small NTAL magnitudes but increase substantially for stations with NTAL RMS of 1 mm or more. No station at all shows a significant anticorrelation between GPS and NTAL even for sub-weekly periodicities, implying that GPS is indeed recording highfrequency variations in NTAL.
For D3, we still find only positive correlations but note an additional distinct dependency with latitude. The higher surface pressure variability can explain this in the extra-tropical latitudes, where frontal systems advected with the mean flow are causing surface pressure variations at periods of around a week. In the tropics, surface pressure variations are instead characterized by distinct periodic variations induced by the atmospheric tides, which are not present in NTAL and which should also not be present in GPS displacements since it is assumed that all tidal signatures are removed during GPS processing in order to avoid temporal aliasing into the daily estimates (Ray et al. 2013). We also note that no obvious systematics were found for station height.
Very similar conclusions are also drawn for the following levels of detail D4-D6 (12 days-5 months). It is important to emphasize that no stations with negative correlations were found even for those long periods. The results are completely different for the inter-annual band containing periods between 1.1 and 3 years. The correlations range between − 0.6 and almost 1 with no obvious systematics regarding latitude, station height, or magnitude of the NTAL signal. The same was also observed for the A9 approximation, which covers long-term trends. This is justified by the fact that GPS displacements are dominated in those longperiods by other geophysical processes than atmospheric surface loading. These include tectonic processes combined with groundwater extraction for the Nepal Himalaya (Fu and Freymueller 2012), sediment compaction or anthropogenic influences for Ganges delta (Brown and Nicholls 2015), groundwater withdrawal for areas close to the Caspian Sea (Joodaki et al. 2014), or hydrological changes in endorheic basins of central Asia (Wang et al. 2018).
We assess how much of the signal in the GPS displacements can be explained by NTAL by calculating the RMS reductions for individual levels of detail (Fig. 4). Decreasing     by the NTAL loading model. NTOL signals far away from the coast are very low, and the correlation between GPS displacements and NTOL-predicted displacements is barely significant (Supplementary Material). Therefore, adding NTOL to NTAL does not change the results significantly for the inland stations specifically selected for this study, demonstrated before by Petrov and Boy (2004). We conclude that continental Eurasian GPS displacements are sensitive to the NTAL loading model from the shortest to the longer periods (D1-D6). Its contribution to the seasonal period band remains dependent on the phase agreement between the GPS-observed and loading-predicted displacements. The NTAL models are not able to explain the long-term trends present in the GPS displacements, confirmed by an almost 0% reduction in the RMS for the A9 approximation. The above results are in line with the study of Memin et al. (2020), who demonstrated that the RMS reduction in GPS displacements reaches more than 5% for periods larger than 2 months (details from D5 onwards) for stations situated in the northern hemisphere at latitudes higher than 50°. They found a much smaller RMS reduction which the different Fig. 3 Pearson correlation coefficients estimated between GPS vertical displacements and vertical displacements predicted by NTAL, computed for corresponding levels of decomposition. Significance was estimated using t-Student test. GPS vertical displacements were corrected for NTOL + HYDL + SLEL before correlation was computed. Correlation is presented for frequency bands defined outside the seasonal band. The area marked with a black rectangle is enlarged in the right plot. Relationships between correlation coefficients and the station's latitude and height are also presented at right, where the dots are color-coded according to the RMS of the loading-predicted vertical displacements estimated for certain frequency bands. Black encircled dots plotted in the figure with detail D1 indicate the 50 stations for which we notice a clear reduction in power spectral density once NTOL + NTAL-loading predicted displacements are removed from the GPS displacements, see Fig. 7. Black encircled dots plotted in the figure with detail D3 indicate the 71 stations for which GPS velocity uncertainty reduces more once the HYDL + SLEL-loading predicted displacements are removed from the GPS displacements instead of the NTOL + NTAL displacements. Yellow encircled dots indicate the 15 stations for which removing the NTOL + NTAL-loading predicted displacements from the GPS displacements brings the largest improvement over HYDL + SLEL choice of period bands might explain, i.e., levels of detail D9 and A9 were considered along with the shorter periods.
Applying alternative NTAL loading models from EOST, i.e., ERAIN, ATMIB, and ATMMO, does not significantly alter the conclusions, thereby supporting the statement that NTAL should be removed from the GPS displacement time series by means of any of the available NTAL models without the need to give a general preference to a particular model provider.

Water storage loading signals in GPS
We subtract NTAL + NTOL from the GPS displacements and compare the residuals with the HYDL + SLEL loadingpredicted displacements. Barely any importance is found for HYDL for periods shorter than a month, i.e., D1-D4 levels of detail. For level D5, we find consistently positive correlations in India and other regions characterized by a humid Monsoon climate (Fig. 5). This even extends to high mountain regions of the Himalayas, where rainfall and snowfall are largely governed by rather moist air masses arriving from the Pacific and Indian oceans, leading to above-average correlations for stations at altitudes between 2000 and 4000 m above sea level. Instead, correlations are particularly poor at latitudes around 40°N that are characterized by a typically arid or semi-arid climate. LSDM is known to have deficits in representing high evaporation rates in arid regions as its evapotranspiration parameterization is based only on the 2 m temperature and ignores humidity, radiation, and wind. The correspondence between GPS displacements and HYDL + SLEL-predicted displacements gradually increases for longer periods and culminates at level D9 which contains the inter-annual changes (Fig. 5). The vast majority of stations exhibit positive correlations, but results for even neighboring stations sometimes vary substantially. We do not identify dependencies with station height and/or latitude.
Instead, it appears that the local conditions of the stations and the proximity to surface water bodies such as rivers, lakes, and artificial reservoirs affect the correspondence of the model and data.
An even more pronounced spatial variability in the results is found for the remaining low-pass filtered signal component A9, where correlations from the 100 stations are almost equally distributed over the interval -1 to 1. No systematic dependencies with height, latitude, or HYDL + SLEL magnitude are identified. We do, however, find a substantial difference between correlation coefficients estimated for NTAL and these estimated for HYDL + SLEL for A9 approximation. GPS displacements for a number of stations situated in central Asia are anti-correlated with HYDL + SLEL, while they are positively correlated with NTAL, showing a sensitivity of GPS stations to NTAL for these longest periods. For stations situated in Nepal, a positive correlation is observed between GPS and HYDL + SLEL, which may arise from snow accumulation and snow melting captured in LSDM. Long-term below average precipitation might also lead to a long-term water storage decrease in LSDM correlated with GPS observations.
On this basis, we state that a part of GPS long-term trends is explained by non-tidal loading, but this contribution is directly related to the location of the GPS station. Since other contributors to the long-term trends, such as tectonic deformation, also influence GPS displacements, recognition of non-tidal loading impact might be biased depending on the magnitude of the loading signal.
We assessed the ability of HYDL + SLEL to explain the residual signal of GPS displacements obtained after subtracting NTAL + NTOL (Fig. 6). We find no improvement at all for the levels D1-D6 that represent periods shorter than a few months. For level D7 (4-9 months), we note positive improvements for more than 70% of the stations, with 40% of them experiencing RMS reductions of more than 40%. This indicates that several stations indeed experience vertical deformations that are representative for regional  Fig. 4 Reduction in RMS of GPS vertical displacements (GPS-NTOL-HYDL-SLEL) once the NTAL loading-predicted vertical displacements were removed. The reduction is shown for NTAL, ERAIN, ATMIB, and ATMMO loading models. Various frequency bands are presented. To demonstrate the impact of removing the NTOL model on the GPS displacements, we also include the RMS reduction once the NTAL + NTOL models are removed from the GPS vertical displacements (GPS-HYDL-SLEL) water storage variations. For the level D8 that contains the annual period, even more stations (> 60%) have RMS reductions larger than 40% after subtracting HYDL + SLEL, meaning that the sensitivity of the GPS displacements to HYDL + SLEL is prominent in the seasonal band. Moreover, the vast majority of the seasonal signals present in the GPS displacements is well explained by hydrology variations.
Very similar conclusions can be drawn by looking at the second hydrological loading data set called MERRA2, provided by the EOST loading service. This supports our Fig. 5 The same as Fig. 3, but for GPS vertical displacements (GPS-NTAL-NTOL) and vertical displacements predicted by HYDL + SLEL findings that hydrology loading is the dominant contributor to the GPS displacements in the seasonal time resolution. In the non-seasonal, inter-annual, or long-term frequency bands, hydrological loading is caused by transient events of heavy rain or flooding that are not mapped into multiyear correlations and RMS analyses. Based on our study, we again cannot suggest favoring any of the loading model providers. The above results are slightly different from the RMS reduction in GPS displacements once GRACE-predicted displacements are removed. This is especially noted for D6 and D7 levels of detail, for which RMS histograms are skewed to the left rather than to the right. This may be caused by the monthly sampling interval of GRACEpredicted displacements, which is really not representative of daily sampling interval provided by the GPS. Another reason might be related to the spatial resolution of GRACE, which is larger than 300 km, and still not representative for GPS displacements reflecting local to regional changes. GRACE-predicted displacements are, however, very helpful in explaining long-term trends, covered by A9 approximation, for which they match almost perfectly the histograms presented for HYDL + SLEL.

Impact of adding barystatic sea-level to hydrology loading
The impact of adding SLEL to HYDL is barely noticeable for the Eurasian inland stations in the seasonal frequency Reduction in RMS of GPS vertical displacements once vertical displacements predicted by HYDL + SLEL, MERRA2 and GRACE were removed. The reduction is presented for various frequency bands band (D7 and D8); correlation coefficients estimated between GPS displacements reduced for NTAL + NTOLpredicted displacements and hydrological loading HYDL agree with those estimated once SLEL is added to HYDL. There are only a few exceptions. The YIBL station (Oman) is located close to the sea (but still > 100 km), for which adding SLEL to HYDL changes the correlation coefficient from positive to slightly negative. For the KKN4 station (Nepal), the correlation coefficient estimates between GPS displacements and hydrology-predicted displacements are highly positive once HYDL is considered alone. Adding SLEL to HYDL results in a large reduction in the correlation coefficient. The correlation coefficients also decrease slightly for the SG31 (China) and MANM (Tajikistan) stations once SLEL is added to HYDL. The correlation coefficients estimated between the GPS and HYDL differ significantly from HYDL + SLEL for nonseasonal changes, however. They are positive for almost all GPS stations once their displacements are correlated with HYDL alone, while they decrease after adding SLEL to HYDL. This is clearly noticeable for stations situated outside major river basins. For stations included in the major river basins, adding SLEL to HYDL has an infinitesimal impact on the GPS-hydrology correlations. Two possible reasons for seasonal agreement and non-seasonal disagreement between HYDL and HYDL + SLEL can be mentioned. The first is that the variance of the HYDL-predicted displacements is very small for non-seasonal frequency bands (Figs. 1, 2). The second is that the SLEL signal occurs due to the mass balance of continental water storage variations in the ocean, and therefore, it is anti-correlated with the HYDL signal. For both reasons, adding any other contributor to HYDL increases the time series variance for the sub-weekly and sub-monthly frequency bands, leading to the change of correlation coefficients.

Significance of correlation coefficients
We employed the t-Student test with a probability of observing the null hypothesis being below 0.05 and found that almost all Pearson correlation coefficients computed between original GPS displacements and displacements predicted by NTAL + NTOL + HYDL + SLEL are significant, starting from the lowest level of decomposition and ending on the final signal approximation (Supplementary Material). This means that the impact of environmental loading on GPS displacements is obvious. Pearson correlation coefficients computed between GPS displacements reduced by NTAL + NTOL and displacements predicted by HYDL + SLEL are insignificant for a majority of GPS stations for the D1 decomposition level. However, with an increasing decomposition level and increasing periods, more correlation coefficients are becoming significant. For decomposition levels D7-A9, all correlation coefficients are significant. This proves a resemblance between GPS displacements and HYDL + SLEL for the seasonal band and long-term changes, as expected. Almost half of the correlation coefficients computed between GPS displacements reduced by NTAL + HYDL + SLEL and displacements predicted by NTOL are insignificant for the lowest levels of decomposition D1-D3 (Supplementary Material). This is a direct result of choosing a continental set of stations and is expected to change once coastal stations are considered. Starting from D4 and on, we notice more and more correlation coefficients became significant, to end with D8, for which only 4 out of 100 stations are characterized by insignificant correlation with GPS displacements. Pearson correlation coefficients computed between GPS displacements reduced by NTOL + HYDL + SLEL and displacements predicted by NTAL prove that there exists a clear correlation between those, starting from the lowest level of decomposition, i.e., D1. This is observed for the northern part of Eurasia, for which the impact of NTAL is expected to be the greatest. For the remaining part of the continent, correlation coefficients are significant for all GPS stations starting from D3 decomposition level. This appears to be true for decomposition levels higher than D3, up to the final approximation A9, with only a few exceptions noticed.

Impact of non-tidal loading on the GPS noise character
The impact of non-tidal loading on the GPS noise character is examined based on power spectral densities, representing the time series power dependent on the period. Figure 7 presents the power spectral densities plotted for four randomly selected Eurasian stations. We note that the contribution of NTOL + NTAL to the GPS displacements differs from the contribution of HYDL + SLEL. For both, a reduction in power in the seasonal frequency band is noticed (1 cycle per year, i.e., cpy, represents an annual signal). The contribution of the HYDL + SLEL-and NTOL + NTAL-predicted displacements is, however, station specific. It shows that beyond HYDL + SLEL and NTOL + NTAL, the contribution of other geophysical processes may be prominent for some stations or the phase difference between the GPS displacements and the loading-predicted displacements precludes the reduction in the annual signal. HYDL + SLEL also contributes to semi-annual and tri-annual peaks, 2 and 3 cpy, which are reduced once the GPS displacements are reduced for HYDL + SLEL. The most interesting impact of non-tidal loading is, however, noticed for GPS displacements of 50% of stations reduced for NTOL + NTAL (Fig. 3). A clear reduction in time series variance is observed for periods between 5 and 80 cpy, pointing out an apparent sensitivity of GPS displacements to the atmospheric loading modeled by NTAL.
A change in the GPS time series variance evokes a change in the noise parameters. To examine those, we use the Maximum Likelihood Estimation method (MLE; Langbein 2012) and assume a time series model containing only trend, annual and semi-annual signals. We find that the noise character of the GPS displacements is whitened (closer to white noise), once they are reduced by non-tidal loading models. The ratio of the spectral indices (GPS/GPS-HYDL-SLEL or GPS/GPS-NTOL-NTAL) is always positive. Its maximum value reaches 1.7 for GPS-NTOL-NTAL. The noise amplitudes are reduced for more than 70% of the stations once non-tidal loading models are applied. The above effects have a combined impact on the uncertainties of GPS velocities. We find that uncertainties of the GPS vertical velocities are being reduced for 100% of the stations we considered when GPS displacements are reduced by non-tidal loadings. For a number of 70% of the stations (Fig. 3), the velocity uncertainty reduces more once HYDL + SLEL is removed from the GPS displacements instead of NTOL + NTAL. For these stations, the velocity uncertainty is improved on average by a factor of 1.7. This impact arises obviously from the fact that HYDL + SLEL reduces the seasonal signals present in the GPS displacements more than NTOL + NTAL. Unmodeled seasonal peaks in the GPS displacements have a larger impact on the GPS velocity uncertainty than a change of the stochastic part (Klos et al. 2018b). Since we model annual and semi-annual signatures in the time series model before the noise parameters are estimated, this statement also comprises the time-variable signatures mismodeled in the time series model. For 30% of the stations, however, the change in the stochastic part induced by NTOL + NTAL is more prominent than the mismodeled time-variable seasonal signal covered by HYDL + SLEL. For these stations, the GPS velocity uncertainty improves on average by a factor of 2.

Summary and conclusions
Improper recognition of the loading impact will bias agreements between various geodetic techniques (Männel et al. 2019). Therefore, NTAL should be applied to the observed vertical earth's crust displacements since it consistently reduces the scatter in the vertical coordinate time series at high frequencies at almost all stations considered. For large continental land masses such as Eurasia, variations in the surface pressure can be quite large. Such pressure variations are present in GPS observed displacements (Jiang et al. 2013). The positive correlation coefficients between GPS displacements and NTAL-predicted displacements, which we notice for all decomposition details, are therefore induced by those surface pressure variations. GPS is proved to be sensitive to even sub-weekly changes in surface pressure. HYDL, however, is most relevant for longer periods: Almost no signal is observed for D1-D6 level of details, with a prominent contribution of HYDL to longer periods of total environmental load NTAL + NTOL + HYDL + SLEL. These results were also confirmed by GRACE-predicted displacements.
For practical GNSS processing, we confirm significant contributions of NTAL to all periods of GPS displacements. We recommend including loading-predicted NTAL at the observation level during GNSS processing. Since all tested NTAL loading models perform equally well, no preference for any model provider can be given. HYDL has been proved to impact the GPS displacements on longer time scales, including mainly seasonal bands and long-term changes. Therefore, the impact of HYDL can be corrected at the observation level but could also be successfully considered at the stage of time series analysis by means of reduction in seasonal amplitudes and longer terms.
This study deliberately left out all stations within 100 km of the coast. A detailed follow-up study would be valuable to assess the quality of NTOL and SLEL, preferably by using coastal GPS stations such as the TIGA stations co-located with coastal tide gauges.
Our results are crucial for understanding a variety of signals included in the GPS displacements and their probable impact on the most-often-used GPS rate, i.e., GPS velocity. We observe an impact on the GPS stochastic character toward white noise, once the loading-induced displacements are directly removed from the GPS displacements (Gruszczynska et al. 2018). This change toward white noise leads to a decrease in velocity uncertainty by 2 times on average.