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 daily-sampled 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 GNSS-related 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 non-tidal 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 long-period 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 short-period 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 GeoForschungsZentrum 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 time-varying 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 (Dill et al. 2018) 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 (Dill and Dobslaw 2019). 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 least-squares 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.

Table 1 Frequency bands retrieved for vertical displacements by 9 decomposition levels
Full size table

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 + NTOL-predicted 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 non-linear variations in 2011. Vertical GPS displacements plotted versus loading-predicted displacements point to a linear relationship between both (Fig. 1).

Fig. 1
figure 1

Vertical displacements observed by GPS and predicted by loading models for GPS station ARTU (Russia). Left: GPS displacements plotted versus loading models. GPS displacements were reduced by respective loading models, as indicated on the plots. Right: Respective scatter plots. A 1:1 relationship is also included for a clarity, pointing to a linear relationship between displacements observed by GPS and predicted by loading models

Full size image
Fig. 2
figure 2

GPS displacements and displacements predicted by the NTAL + NTOL + HYDL + SLEL, HYDL + SLEL and NTAL + NTOL loading models for ARTU GPS station (Russia), divided into different decomposition levels, from D1 to A9. The anomalies in 2019 are an artifact of modeling a gap present in the GPS displacements with a few outlying observations in the middle of the gap. For the purpose of HYDL + SLEL validation, GRACE-predicted displacements are also plotted

Full size image

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 anti-correlation between GPS and NTAL even for sub-weekly periodicities, implying that GPS is indeed recording high-frequency variations in NTAL.

Fig. 3
figure 3figure 3figure 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

Full size image

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 long-periods 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 RMS indicates the ability of the NTAL model to explain the GPS displacements. Thus, distribution functions heavily skewed to the right would document NTAL model performance suitable for geodetic applications. For levels of detail between D1 and D5, a positive impact of NTAL was found for more than 60% of the GPS stations. For the seasonal frequency band, i.e., levels of detail D7 and D8, the RMS of the GPS displacements reduced by NTAL decreased for 50% of the stations. Stations for which significant increases in RMS values were noticed are characterized by a phase difference between the signals generating the D7 and D8 levels of detail for GPS displacements and those predicted 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 choice of period bands might explain, i.e., levels of detail D9 and A9 were considered along with the shorter periods.

Fig. 4
figure 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)

Full size image

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

Fig. 5
figure 5figure 5figure 5

The same as Fig. 3, but for GPS vertical displacements (GPS-NTAL-NTOL) and vertical displacements predicted by HYDL + SLEL

Full size image

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

Fig. 6
figure 6

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

Full size image

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 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 multi-year 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 GRACE-predicted 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 band (D7 and D8); correlation coefficients estimated between GPS displacements reduced for NTAL + NTOL-predicted 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 non-seasonal 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.

Fig. 7
figure 7

Power spectral densities plotted for four randomly selected GPS stations. GPS displacements are plotted in black. GPS displacements reduced by HYDL + SLEL-loading predicted displacements are plotted in green. GPS displacements reduced by NTOL + NTAL-loading predicted displacements are plotted in red

Full size image

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.

Acronyms

ATMIB (atmospheric loading), ATMMO (atmospheric and induced oceanic loading), ECMWF (European Centre for Medium-Range Weather Forecasts), EOST (Ecole et Observatoire des Sciences de la Terre, Strasbourg, France), ERAIN (ERA-Interim), ESMGFZ (Earth-System-Modeling Group of Deutsches GeoForschungsZentrum), GNSS (Global Navigation Satellite System), GPS (Global Positioning System), GRACE (Gravity Recovery and Climate Experiment), HYDL (hydrological loading), LSDM (Land Surface Discharge Model), MERRA (Modern Era Retrospective-Analysis), NGL (Nevada Geodetic Laboratory), NTAL (non-tidal atmospheric loading), NTOL (non-tidal oceanic loading), RMS (root-mean-square), SLEL (barystatic sea-level changes), TUGO-m (Toulouse Unstructured Grid Ocean model).