Effect of Surface Mass Loading on Geodetic GPS Observations

We investigated the effect of mass loading (atmospheric, oceanic and hydrological loading (AOH)) on Global Positioning System (GPS) height time series from 30 GPS stations in the Eurasian plate. Wavelet coherence (WTC) was employed to inspect the correlation and the time-variable relative phase between the two signals in the time–frequency domain. The results of the WTC-based semblance analysis indicated that the annual fluctuations in the two signals for most sites are physically related. The phase asynchrony at the annual time scale between GPS heights and AOH displacements indicated that the annual oscillation in GPS heights is due to a combination of mass loading signals and systematic errors (AOH modelling errors, geophysical effects and/or GPS system errors). Moreover, we discuss the impacts of AOH corrections on GPS noise estimation. The results showed that not all sites have an improved velocity uncertainty due to the increased amplitude of noise and/or the decreased spectral index after AOH corrections. Therefore, the posterior mass loading model correction is potentially feasible but not sufficient.


Introduction
The strong annual oscillations in Global Positioning System (GPS) observations are known to be driven by the redistribution of environmental mass loading [1,2].The atmospheric, oceanic and hydrological loading (AOH) is a load that is likely to be important [3][4][5].Surface deformation is usually large enough to be monitored using the space geodetic techniques [6][7][8] and can cause a significant reduction in the accuracy of GPS observations [9][10][11].Globally distributed surface mass variabilities can be estimated to a certain extent using the method of a mass loading model and Gravity Recovery and Climate Experiment (GRACE) technology [12,13].Van Dam et al. [14] corrected 51 GPS sites in Europe with AOH derived from GRACE and found that only 10 GPS stations have a reduction in their root mean square (RMS) values.Compared to GRACE, the mass loading model can estimate the small-scale variability [15].Recently, Jiang et al. [16] adjusted 233 globally distributed GPS sites by AOH derived from the mass loading model, provided by the Global Geophysical Fluid Center (GGFC).The results indicated that 64% of the stations underwent a decrease in their RMS values, although the reasons for the inconsistency between the mass loading displacements and GPS observations were not discussed in detail.Considering that the RMS analysis is only used in the time domain, wavelet coherence (WTC) is a useful exploratory tool to expose regions with high correlations and to further reveal the relative phases between two time series in the time and frequency domain.In this study, we adopted the mass loading model to correct GPS observations using RMS analysis and WTC.Moreover, we will discuss the effect of mass loading correction on site velocity uncertainty and noise estimation.

GPS Data
We used daily GPS data provided by the International Global Navigation Satellite System Service (IGS), which combines the GAMIT (SPOCA) and GIPSY (JPL) solutions.The GAMIT and GIPSY solutions contain the unconstrained geodetic station coordinates, which have been corrected for some geophysical effects (Table 1).The combined solutions are used so that random processing noise introduced by any individual analysis will be compressed to enhance the signal-to-noise ratio.Fritsche et al. [2] confirmed that the scatter in the vertical component for the combined solutions is smaller by 5-15 mm compared to the scatter in any other analyses.We removed offsets, outliers and the linear trend as described in Nikolaidis in GPS observations [17].Due to the removal of offsets and outliers as well as other reasons (e.g., the equipment changes or failures), missing data is an inevitable problem in GPS observations.At present, it is difficult to use interpolation methods to calculate the true value due to the complexity of GPS, although the filling of missing data with interpolation will introduce an artificial error.Therefore, GPS site selection is performed using the following criteria.First, the time span of GPS time series should be longer than 2.5 years [18].Additionally, to minimize the interpolation effect, the number of continuous data gaps should not exceed 15 and the ratio of missing data to entire observations must not be more than 1% for each site [19].Finally, the selected stations should be uniformly distributed.In this analysis, we filled the data gaps using linear interpolation.It should be noted that the mass loading signal still remains in the GPS observation.Considering that the Eurasian plate is a tectonically stable region and the tectonic noise is minimal, we selected 30 GPS stations located in the Eurasian plate for this analysis from January 2010 to January 2014 (Figure 1).

Mass Loading
For the non-tidal atmospheric loading deformation, we adopted the provisional product provided by the International Earth Rotation and Reference Systems Service (IERS) Global

Mass Loading
For the non-tidal atmospheric loading deformation, we adopted the provisional product provided by the International Earth Rotation and Reference Systems Service (IERS) Global Geophysical Fluids Center (GGFC)-German Research Centre for Geosciences (GFZ), which is available at http://www.gfzpotsdam.de/en/esmdata.Non-tidal atmospheric mass loads with 0.5 • × 0.5 • spatial resolution and 6-hourly temporal resolution were estimated by using the surface pressure from the European Centre for Medium-Range Weather Forecasts (ECMWF) operational model [15].For non-tidal ocean loading deformation, we also adopted the provisional product provided by the GFZ.Non-tidal ocean loads with 0.5 • × 0.5 • spatial resolution and 6-hourly temporal resolution were derived from non-tidal ocean bottom pressure from the Ocean Model for Circulation and Tides (OMCT) [15].The GFZ processes the consistent deformation fields from the daily means of the ECMWF atmospheric surface pressure and the OMCT non-tidal ocean bottom pressure to derive the total non-tidal deformation signal.For hydrological loading (HYDL) deformation, we used the provisional product provided by the GFZ.The HYDL with 0.5 • × 0.5 • spatial resolution and the daily temporal resolution was predicted from the Global Land Surface Discharge Model (LSDM) [20].The LSDM was found to be composed of the masses of soil moisture, snow and runoff from glaciers, especially water in rivers and lakes.To maintain consistency with GPS heights, the mass loading displacements were all calculated in the center of the figure frame (CF).

Wavelet Coherence
Firstly, we will briefly outline the continuous wavelet transform (CWT).The CWT of a time series (x n , n = 1, 2, ..., N) is defined as follows [21]: where W X n (s) is the wavelet coefficient, δt is the time step, s and ϕ 0 are the wavelet scale and the mother function, respectively, and n and √ δt/s represent the reversed time and the normalization factor, respectively.The idea behind the wavelet transform is to apply the wavelet as a band-pass filter to the time series where the wavelet scale is linearly related to the characteristic period of the filter.Based on the CWT, the WTC of the two-time series X n and Y n (n = 1, 2, ..., N) is defined as follows [21]: where S is the smoothing operator, W X n (s) is the wavelet coefficient of the time series X n , W Y n (s) is the wavelet coefficient of the time series Y n , W XY n (s) is the cross wavelet transform of the two time series X n and Y n , the complex argument arg(W XY n (s)) is the local relative phase between X n and Y n in time frequency space, s is the wavelet scale, R n (s) is the localized correlation coefficient and n is the number of stations.The advantage of WTC is its ability to locate regions in the time-frequency domain where the two time series are significantly correlated, and further reveal the relative phase between the two time series in the time-frequency space.In our analysis, the relative phasing of the two time series can be quantified by using the circular mean of the phase or the WTC-based semblance.The definition of the circular mean of a set of angles (α i , i = 1, 2, ..., N) is expressed as follows [21]: Meanwhile, the WTC-based semblance is defined as follows [22]: The value of ρ ranges from −1 (inversely correlated) to 0 (uncorrelated) to 1 (fully correlated).In this analysis, we used the WTC-based semblance as the key parameter to compare the two time series.

Root Mean Square Analysis and Correlation Analysis
Given that the mass loading behaves in the vertical direction, we only considered the vertical component in the following analysis.
Figure 2 shows the time series of GPS heights of four stations along with AOH-derived vertical displacements.It is clear from Figure 2 that GPS heights and AOH displacements prominently have an annual variation pattern.AOH displacements track GPS heights well at the ANKR, BOR1 and LHAZ sites, indicating that AOH may be responsible for the non-linear vertical motion of GPS stations.However, for the KUNM site, the AOH deformation signal is not in agreement with GPS height signal due to the phase asynchrony between these two time series.
Meanwhile, the WTC-based semblance is defined as follows [22]: The value of  ranges from −1 (inversely correlated) to 0 (uncorrelated) to 1 (fully correlated).In this analysis, we used the WTC-based semblance as the key parameter to compare the two time series.

Root Mean Square Analysis and Correlation Analysis
Given that the mass loading behaves in the vertical direction, we only considered the vertical component in the following analysis.
Figure 2 shows the time series of GPS heights of four stations along with AOH-derived vertical displacements.It is clear from Figure 2 that GPS heights and AOH displacements prominently have an annual variation pattern.AOH displacements track GPS heights well at the ANKR, BOR1 and LHAZ sites, indicating that AOH may be responsible for the non-linear vertical motion of GPS stations.However, for the KUNM site, the AOH deformation signal is not in agreement with GPS height signal due to the phase asynchrony between these two time series.From Table 2, we can see that the RMS value of GPS height time series varies from 5.8 to 14.5 mm, whereas the RMS value of AOH deformation time series ranges from 3.1 to 7.7 mm.Therefore, AOH cannot fully explain GPS height fluctuations.The Pearson correlation analysis was employed to investigate the correlation between GPS heights and AOH displacements.As shown in Figure 3a, GPS height and AOH deformation are positively correlated for all the selected sites, while GPS height and AOH deformation have a strong correlation for most sites (correlation coefficient is greater than 0.6).The strong correlation between GPS heights and AOH displacements may be attributed to the significantly related annual signal.Apart from the annual component, the other components in AOH deformation may act as the potential candidates for determining the strong correlation between GPS heights and AOH displacements.Therefore, by adopting the least squares method, the annual term was removed from GPS heights and AOH displacements to generate GPS height residuals and AOH deformation residuals, respectively.Figure 3b presents the correlation coefficients between GPS height residuals and AOH deformation residuals.All the selected sites are positively correlated, and most sites have a medium correlation (correlation coefficient is between 0.4 and 0.6).Therefore, the correlation between GPS height and AOH deformation is reduced after removing the annual signal.As shown in Figure 3a, GPS height and AOH deformation are positively correlated for all the selected sites, while GPS height and AOH deformation have a strong correlation for most sites (correlation coefficient is greater than 0.6).The strong correlation between GPS heights and AOH displacements may be attributed to the significantly related annual signal.Apart from the annual component, the other components in AOH deformation may act as the potential candidates for determining the strong correlation between GPS heights and AOH displacements.Therefore, by adopting the least squares method, the annual term was removed from GPS heights and AOH displacements to generate GPS height residuals and AOH deformation residuals, respectively.Figure 3b presents the correlation coefficients between GPS height residuals and AOH deformation residuals.All the selected sites are positively correlated, and most sites have a medium correlation (correlation coefficient is between 0.4 and 0.6).Therefore, the correlation between GPS height and AOH deformation is reduced after removing the annual signal.

Wavelet Coherence Analysis
Considering that the Pearson correlation analysis is aimed to simply find a correlation in the time domain, there are several issues in determining the specific relationship between GPS heights and AOH displacements with this analysis.Firstly, since the annual component plays an essential role, a more detailed insight aimed at the annual scale (frequency) should be considered.Secondly, in addition to the annual component, other components may potentially be relevant.Lastly, the relationship (e.g., correlation coefficient and relative phase) may change over time.Here, we propose to employ the WTC to overcome all the above-mentioned issues.

Wavelet Coherence Analysis
Considering that the Pearson correlation analysis is aimed to simply find a correlation in the time domain, there are several issues in determining the specific relationship between GPS heights and AOH displacements with this analysis.Firstly, since the annual component plays an essential role, a more detailed insight aimed at the annual scale (frequency) should be considered.Secondly, in addition to the annual component, other components may potentially be relevant.Lastly, the relationship (e.g., correlation coefficient and relative phase) may change over time.Here, we propose to employ the WTC to overcome all the above-mentioned issues.
The WTC spectrum of GPS height time series and AOH deformation time series for the BOR1 site and KUNM site are shown in Figure 4.The reason that we chose these two stations is that Figure 1 displays the approximate relative phase between GPS heights and AOH displacements for the BOR1 site and KUNM site, which can be used to validate the results of WTC.The thick contour indicates the 95% confidence interval against the red noise, while the phase of the area within the cone of influence was unaccounted for due to the edge effects.The relative phasing relationship is denoted by arrows (with in-phase pointing right, anti-phase pointing left and GPS height leading AOH deformation to 90 • pointing straight down).Figure 4 shows that the correlation between these two time series is considerably high at a frequency that is closest to 1 cycle per year (cpy).From Figure 4a, we can see that GPS height time series and AOH deformation time series are in phase at the period closest to 1 year for the BOR1 site, indicating that AOH is responsible for the annual fluctuation in GPS heights.From Figure 4b, we can see that GPS height time series and AOH deformation time series are not in phase at the period closest to 1 year for the KUNM site, indicating that AOH is not responsible for the annual fluctuation in GPS heights.For the high-frequency components (in the period band of fewer than 256 days), although a locally significant correlation was detected, the relative phase angles of the high-frequency components fluctuated sharply, indicating that the GPS height signal and AOH deformation signal are not in phase.Consequently, it is demonstrated that the geophysical effects cannot explain the sub-seasonal oscillations in GPS observations [23].
4a, we can see that GPS height time series and AOH deformation time series are in phase at the period closest to 1 year for the BOR1 site, indicating that AOH is responsible for the annual fluctuation in GPS heights.From Figure 4b, we can see that GPS height time series and AOH deformation time series are not in phase at the period closest to 1 year for the KUNM site, indicating that AOH is not responsible for the annual fluctuation in GPS heights.For the high-frequency components (in the period band of fewer than 256 days), although a locally significant correlation was detected, the relative phase angles of the high-frequency components fluctuated sharply, indicating that the GPS height signal and AOH deformation signal are not in phase.Consequently, it is demonstrated that the geophysical effects cannot explain the sub-seasonal oscillations in GPS observations [23].Since the fluctuation range of the time-variable relative phasing at the period closest to 1 year is small, the mean WTC-based semblance was computed to quantitatively express the relative phase between these two signals.Following the same method as depicted in Figure 4, the WTC was used to analyze the relative phase between GPS height time series and AOH deformation time series for the other sites, with the mean WTC-based semblance at the frequency closest to 1 cpy.
As shown in Figure 5, the average value of the WTC-based semblance for most, but not all, sites are close to 1, indicating that AOH is not a complete representation of the driving force for the nonlinear GPS site motion.From Figure 5, we can see that for the SHAO site, the mean semblance is only 0.4, while the other 10 sites (NILI, BJFS, CHAN, KUNM, ONSA, POTS, SOFI, URUM, VILL and WUHN) have a mean semblance value that is close to 0.7-0.9.Other geophysical effects (e.g., groundwater), GPS systematic errors (e.g., inaccurate local multi-path and tropospheric modelling) and AOH modelling errors (e.g., interpolation method) can lead to the phase asynchrony between GPS observations and mass loading displacements.A measure of RMS reduction ratio expressed in percent (ΔRMS%) is employed to quantify the correction performance of AOH models: Since the fluctuation range of the time-variable relative phasing at the period closest to 1 year is small, the mean WTC-based semblance was computed to quantitatively express the relative phase between these two signals.Following the same method as depicted in Figure 4, the WTC was used to analyze the relative phase between GPS height time series and AOH deformation time series for the other sites, with the mean WTC-based semblance at the frequency closest to 1 cpy.
As shown in Figure 5, the average value of the WTC-based semblance for most, but not all, sites are close to 1, indicating that AOH is not a complete representation of the driving force for the non-linear GPS site motion.From Figure 5, we can see that for the SHAO site, the mean semblance is only 0.4, while the other 10 sites (NILI, BJFS, CHAN, KUNM, ONSA, POTS, SOFI, URUM, VILL and WUHN) have a mean semblance value that is close to 0.7-0.9.Other geophysical effects (e.g., groundwater), GPS systematic errors (e.g., inaccurate local multi-path and tropospheric modelling) and AOH modelling errors (e.g., interpolation method) can lead to the phase asynchrony between GPS observations and mass loading displacements.A measure of RMS reduction ratio expressed in percent (∆RMS%) is employed to quantify the correction performance of AOH models: Table 3 indicates that all sites experience a reduction in their RMS values when corrected by Table 3 indicates that all sites experience a reduction in their RMS values when corrected by AOH, revealing the feasibility of correcting GPS height with the mass loading model.From Table 3, we can see that a relatively poor improvement occurs at the SHAO, VILL and WUHN sites, where there is a significant phase asynchrony between GPS heights and AOH displacements.

Noise and Velocity Analysis
Errors in the mass loading model are currently difficult to quantify.Nonetheless, they will certainly enter into our estimates of surface deformation [24].In this case, aside from investigating the RMS reduction of GPS position time series, discussing the effect of mass loading corrections on the velocity and noise estimates is equally important.In this study, we adopted the Hector software to determine the velocity and its uncertainty with pre-assumed power-law noise and white noise [25].The software can sufficiently describe the stochastic properties of most daily GPS position time series.The noise parameters, including the spectral indices and noise amplitudes, were also estimated before and after the AOH corrections by using the maximum likelihood estimation to analyze the change in noise properties.The accuracy of velocity is calculated based on the following formula [26]: where m v is the accuracy of velocity, A PL is the amplitude of noise, K is the spectral index, Γ is the gamma function, N is the number of data in the time series and ∆T is the sampling interval.The middle histogram in Figure 6 demonstrates that the power-law noise amplitude is mostly decreased when AOH corrections are applied.The small amplitude of noise after AOH corrections indicates the small scatter of residuals, and consequently a considerably stable solution for the determination of site velocity.Meanwhile, the bottom histogram in Figure 6 indicates that the spectral index shifts from a flicker noise (k = −1) to white noise (k = 0) in most sites.Hence, the loading signals may be superimposed on the stochastic process realization.The removal of AOH deformation for GPS stations causes a reduction in median velocity error from 1.17 to 0.87 mm/year in median power-law noise amplitude from 14.99 to 11.84 mm/year-k/4 and in the median spectral index from −0.84 to −0.81 mm/year-k/4.The results indicate that AOH corrections altered the noise properties of GPS height time series.Based on Equation (6), changes in the noise properties will certainly influence the estimates of velocity uncertainties.Accordingly, AOH corrections will potentially affect the velocity-field determination, and thus, its geodynamical interpretation.From the top histogram in Figure 6, we can see that for 24 out of the 30 sites the velocity uncertainty estimates decrease after Appl.Sci.2018, 8, 1851 8 of 10 AOH corrections.However, for 6 out of the 30 sites (KUNM, PADO, SOFI, TLSE, VILL and YEBE) the velocity uncertainty increases due to the increased amplitude of noise and/or the decreased spectral index after AOH corrections.This indicates that errors or noises are introduced when implementing the mass loading correction, suggesting that the correction of GPS heights with the mass loading model is potentially feasible but not sufficient.
GPS stations causes a reduction in median velocity error from 1.17 to 0.87 mm/year in median powerlaw noise amplitude from 14.99 to 11.84 mm/year-k/4 and in the median spectral index from −0.84 to −0.81 mm/year-k/4.The results indicate that AOH corrections altered the noise properties of GPS height time series.Based on Equation ( 6), changes in the noise properties will certainly influence the estimates of velocity uncertainties.Accordingly, AOH corrections will potentially affect the velocityfield determination, and thus, its geodynamical interpretation.From the top histogram in Figure 6, we can see that for 24 out of the 30 sites the velocity uncertainty estimates decrease after AOH corrections.However, for 6 out of the 30 sites (KUNM, PADO, SOFI, TLSE, VILL and YEBE) the velocity uncertainty increases due to the increased amplitude of noise and/or the decreased spectral index after AOH corrections.This indicates that errors or noises are introduced when implementing the mass loading correction, suggesting that the correction of GPS heights with the mass loading model is potentially feasible but not sufficient.

Conclusions
Considering the limitations of the RMS analysis-i.e., that it can only be used in the time domain-in this paper, WTC was applied to identify the correlation and relative phases between GPS height time series and AOH deformation time series in the time-frequency space.The results of the WTC-based semblance analysis indicate that AOH is not a complete representation of the driving force for the non-linear GPS site motions.AOH modelling errors (e.g., interpolation method), other geophysical effects (e.g., groundwater) and GPS system errors (e.g., local multi-path and tropospheric mis-modelling) can cause the phase asynchrony between GPS observations and AOH displacements.

Conclusions
Considering the limitations of the RMS analysis-i.e., that it can only be used in the time domain-in this paper, WTC was applied to identify the correlation and relative phases between GPS height time series and AOH deformation time series in the time-frequency space.The results of the WTC-based semblance analysis indicate that AOH is not a complete representation of the driving force for the non-linear GPS site motions.AOH modelling errors (e.g., interpolation method), other geophysical effects (e.g., groundwater) and GPS system errors (e.g., local multi-path and tropospheric mis-modelling) can cause the phase asynchrony between GPS observations and AOH displacements.Furthermore, the impact of AOH corrections on noise estimates were investigated.The results showed that, for most sites, the velocity uncertainty is improved after AOH corrections, whereas the velocity uncertainty increases for some sites due to the large amplitude of noise and/or the decreased spectral index after AOH corrections.Overall, this indicates that the posterior mass loading model correction is potentially feasible but not sufficient.In conclusion, the wavelet coherence analysis is a useful exploratory tool to identify correlated behaviors and the time-variable relative phases between GPS observations and mass loading displacements.

10 Figure 1 .
Figure 1.Spatial distribution of the 30 GPS stations used in this study.

Figure 1 .
Figure 1.Spatial distribution of the 30 GPS stations used in this study.

Figure 3 .
Figure 3. (a) Correlation between the GPS height and the AOH deformation; (b) Correlation between the GPS height residuals and AOH deformation residuals.

Figure 3 .
Figure 3. (a) Correlation between the GPS height and the AOH deformation; (b) Correlation between the GPS height residuals and AOH deformation residuals.

Figure 4 .
Figure 4.The wavelet coherence (WTC) spectrum of GPS heights and AOH displacements for (a) BOR1 site and (b) KUNM site.

Figure 4 .
Figure 4.The wavelet coherence (WTC) spectrum of GPS heights and AOH displacements for (a) BOR1 site and (b) KUNM site.

10 Figure 5 .
Figure 5. Mean value of the WTC-based semblance at the period closest to 1 year for 30 sites.

Figure 5 .
Figure 5. Mean value of the WTC-based semblance at the period closest to 1 year for 30 sites.

Figure 6 .
Figure 6.Comparison of the velocity uncertainty and noise parameters before and after AOH corrections.

Figure 6 .
Figure 6.Comparison of the velocity uncertainty and noise parameters before and after AOH corrections.

Table 1 .
Strategies and geophysical models applied when processing GPS time series with GAMIT and GIPSY.

Table 2 .
The RMS value of GPS height time series and AOH deformation time series respectively.

Table 3 .
The RMS reduction ratio in percent when correcting GPS heights with AOH displacements.