Geophysical interpretation of GPS loading deformation over western Europe using GRACE measurements

We analyze more than 10 years of Global Positioning System (GPS) height residuals and vertical displacements predicted from surface mass loading observed by the Gravity Recovery and Climate Experiment (GRACE) for 36 International GNSS Service (IGS) stations over Europe. Seasonal surface displacements, mostly due to atmospheric and hydrological loading, are significant in both GPS and GRACE measurements. With an extended time period, our new analysis based on release 05 GRACE data from Center for Space Research (CSR) shows considerably improved agreement between GPS and GRACE than that from previous studies, for not only annual but also interannual signals. The GPS height residual series at most stations exhibit reduced weighted root-meansquares (WRMS) after removing GRACE-derived vertical displacements, which is attributed to improved accuracy of both GPS and GRACE data products. Furthermore, we demonstrate the necessity of reducing leakage bias in GRACE estimates for the study of surface loading deformation using GRACE satellite gravity observations.


Introduction
Mass transport and redistribution within the Earth system, including the atmosphere, hydrosphere, ocean, cryosphere, mantle and core, are driving forces of temporal variations of geodetic observables, such as gravity change, Earth's rotation, geocenter motion and elastic surface displacements [Mangiarotti et al. 2001].These geodetic observables can be measured by modern space geodetic techniques with increasing and unprecedented accuracy.Surface displacements can be determined fairly well by GPS observations and other remote sensing techniques (e.g., InSAR), owing to continuous improvements in data quality and data processing methods.
Time series generated from continuous GPS observations reveal significant non-linear variations that contain contributions from various sources.One of the main sources is surface mass loading that includes atmospheric, oceanic and hydrological effects.The surface loading deformation can be estimated using predictions from numerical models for different components (e.g., atmosphere, ocean and hydrosphere) of the climate system.Atmospheric surface pressure estimates from advanced atmospheric general circulation models are reasonably accurate owing to the large amount of available in situ meteorological observations that can be assimilated into the atmospheric models [Kalnay et al. 1996].However, model-estimated ocean bottom pressure (OBP) and terrestrial water storage (TWS) changes, especially the latter, are believed to be with large uncertainty, due to the lack of adequate in situ observations and immaturity of the models [Döll et al. 2003, Rodell et al. 2009].
The Gravity Recovery and Climate Experiment (GRACE) satellite gravity mission, whose primary objective is to measure global mass redistribution through its integrated gravitational effect, provides a unique means for measuring mass redistribution on global basis with unprecedented accuracy [Tapley et al. 2004].GRACE has accumulated over 14 years of monthly measurements of the global surface mass change, including TWS change, oceanic mass change, and ice mass change over polar ice sheets and mountain glaciers [Cazenave and Chen 2009].Although the spatial resolution of GRACE observations is not comparable with those of model estimates, the global coverage, extended record (over 14 years now), and increasing accuracy have made GRACE a useful data resource for studying global mass variations and their various applications.
Some previous studies have compared GPS dis-placements with GRACE measurements, in order to interpret the loading deformations observed in GPS data, especially for those in vertical direction.On one hand, some researchers have found good agreement.Davis et al. [2004] found good consistency between GPS and GRACE-derived height time series for the sites in the Amazon river basin, and demonstrated that GRACE can be used to monitor deformation of the Earth's surface.Tregoning et al. [2009] also reported good agreement in different regions (for seasonal deformations with both large and small amplitudes) globally, and attributed these to reductions of spurious signals in GPS analysis.Fu and Freymueller [2012] found close correlations in seasonal variations between GPS and GRACE in the Nepal Himalaya region.Valty et al. [2013] showed reasonable agreement between GPS and GRACE results at interannual time scales in southern Europe.On the other hand, some researchers have found discrepancies between GPS and GRACE.For instance, van Dam et al. [2007] showed that the annual amplitudes of vertical deformation signals varied considerably in GPS heights from site to site over western Europe, and deviated from GRACE-derived vertical displacements.They suggested that the discrepancy was mainly due to technique errors in GPS data processing.Similarly, King et al. [2006] also found poor agreement between GPS and GRACE observations, and indicated that the discrepancy might be due to errors in GPS data and un-modeled non-geophysical effects in GPS time series.
In this paper, we reanalyze the vertical displacements in western Europe using GPS observations and GRACE surface loading measurements.The signal amplitudes in this region are relatively small, where previous studies [e.g., van Dam et al. 2007] showed large discrepancy between GPS and GRACE observations.We focus on the seasonal (mainly annual) and interannual variations in vertical displacements using updated GPS data, and compare them with the new version (release-5, RL05) of GRACE data (previous studies in the region were based on GRACE RL04 or earlier solutions), with the consideration of degree-1 geocenter terms provided by Swenson et al. [2008].We also apply a scaling factor method to correct the amplitude attenuation effect in GRACE-derived vertical loading displacements, in order to find out whether better agreement in amplitudes could be obtained between GPS and GRACEderived height time series.

GRACE data
We use GRACE RL05 GSM time-variable gravity solutions provided by the University of Texas Center for Space Research (UTCSR) for the period from August 2002 to October 2014 [Bettadpur 2012].Each monthly solution consists of fully normalized spherical harmonic Stokes coefficients to degree and order 60.We have replaced the GRACE C 20 coefficients by independent estimates from the satellite laser ranging (SLR) provided by UTCSR [Cheng et al. 2013].As atmospheric and oceanic effects have been removed during GRACE data processing through the dealiasing process [Bettadpur 2012], GRACE GSM spherical harmonic coefficients mostly represent contributions from terrestrial water mass on land (e.g., soil moisture, snow and ice, and groundwater) and residual mass change over the ocean.
Vertical displacements due to hydrological loading are calculated from GRACE coefficients as follows [van Dam et al. 2007] where, dr(i,z) is the vertical displacement of the Earth's surface, R is the radius of Earth, P ~l,m are normalized Legendre functions of degree l and order m, and C lm and S lm are the Stokes coefficients.h´l and k´l are the load Love numbers of degree l, which are used in this study with the values defined in a Gutenberg-Bullen (G-B) Earth model provided by Farrell [1972].W l is the weighting function of Gaussian smoothing [Jekeli 1981, Wahr et al. 1998], with a filtering radius of 300 km in this study.Note that the l=1 terms that represent geocenter motion are not included in the GRACE solutions.

GPS data
In this analysis, we use the JPL GIPSY solutions of daily GPS heights for 36 sites in Europe (Bock and Webb [2012]; ftp://sopac-ftp.ucsd.edu/pub/timeseries/measures/ats/).The locations of used GPS sites are shown in Figure 1.The GPS data sets cover the same period as GRACE (2002GRACE ( .08 -2014.10.10).JPL's reanalysis orbit and clock products are determined using a set of models in a free-network reference frame, including gravity from Earth, Sun, Moon, and other planets, DE421 planetary ephemeris, GSPM10 satellite solar pressure model, IAU06 model for precession and nutation, IERS2010 tides model, FES2004 ocean loading model, and IGS satellite and receiver antenna phase center models.Outliers, mean, trend and coseismic jumps are removed from the time series.There are more than 1100 global GPS sites from the JPL GIPSY products, and about 100 sites located in Europe.For comparison purposes, we select the 36 sites in western Europe that have been included in previous studies [e.g., van Dam et al. 2007].

Seasonal and interannual vertical displacements from GPS and GRACE
In order to compare GPS results with GRACE estimated hydrological loading deformation (based on GRACE GSM gravity solutions as described in 2.1), the GPS data have to be corrected for atmospheric and non-tidal oceanic loading.We derive the atmospheric and non-tidal oceanic loading vertical displacements using the GRACE AOD1B dealiasing products [Flechtner et al. 2014].In particular, when using Equation (1) here, the degree-1 components of vertical displacement should be calculated using corresponding Stokes coefficients and Love numbers that are converted to values in the center of the figure (CF) reference frame.The conversion of Stokes coefficients in different reference frames has been discussed in Swenson et al. [2008], and the different values of degree-1 load Love numbers in different reference frames have been summarized in Blewitt [2003].
For the sake of comparison, the GPS and GRACE data should be represented in a consistent reference frame.We calculate the vertical displacement using degree-1 "GSM-like" Stokes coefficients provided by GRACE Tellus website (http://grace.jpl.nasa.gov/;Swenson et al. [2008]) and remove them from the GPS heights so that both GPS and GRACE data sets are now in the CM (center of mass) reference frame.Note that the degree-1 Stokes coefficients are defined in the CF frame and thus the corresponding load Love number should also be used in this frame.
We analyze the periodic components of the displacement series by estimating their amplitudes and phases from least squares fitting.Using a weighted least squares fit, first we extract the annual and semiannual amplitudes and phases of GPS height residuals (after removing the atmospheric and oceanic loading effects), then we derive those from the GRACE-predicted vertical displacements.The results for annual components are listed in Table 1.The semiannual component shows considerably weaker amplitudes than annual ones, thus are not included for brevity sake of Table 1. Figure 2 shows the GPS height residuals and GRACE-derived vertical displacements due to hydrological loading in 8 selected sites.Annual variations dominate in both GPS and GRACE-derived time series, and the two results agree quite well in both amplitudes and phases for the annual terms at most of the selected sites.
To compare the results at interannual time scales, we apply a 13-month moving average (with half weight assigned to the first and last values and full weight to the rest) to both GPS and GRACE monthly time series, and then compute correlation coefficients and weighted root-mean-square (WRMS) reductions for the 8 selected stations (with the coefficient values and percentages shown at the bottom of each subplot in Figure 3).The correlation coefficients are computed at zero phase lags, reflecting the correlation between two observations affected by the same loading effects (i.e., hydrological loading in this study).The WRMS of the GPS height or GRACE-derived vertical displacement of a given site is computed as Then the percentage of WRMS reduction reflecting the degree of agreement between GPS and GRACE reads where x ij is the vertical displacement, x − i is the weighted mean value of it and P ij is the weight at site i in month j.
We show the phasor diagram in Figure 4, with arrows representing the amplitudes and phases of annual vertical displacements from GPS or GRACE data at each of the 36 selected GPS sites.The estimated errors WRMS GPS or GRA RACE P i (3)   of annual amplitudes are also presented in Figure 4 by circles at the arrow heads.The error of annual GPS height amplitude at each site is estimated from the GPS formal error provided by the JPL product, through error propagation in least squares.We estimate the error of GRACE derived displacements by analyzing the residual Stokes coefficients following Wahr et al.
[2006] and applying the error propagation law in the loading deformation calculation with spherical harmonic series as well as in the least squares fitting of the seasonal signals.
The amplitudes and phases of GPS observed annual vertical displacements are less spatially coherent than GRACE results, especially in coastal regions (Figure 4).However, the mean amplitudes and phases (shown by red arrows in the upper left corners in the two subfigures of Figure 4) of GPS and GRACE estimates agree well, i.e. on the order of 2 mm in amplitude and peakvalue month between June and July in phase.
The weighted root-mean-square (WRMS) reduction is commonly used to verify whether the GPS time series agree well with GRACE hydrological loading contributions [van Dam et al. 2007].Figure 5 shows the WRMS reductions when GRACE-derived vertical displacements have been removed from GPS height residuals.For the sites where the WRMS reduction values are positive (marked with warm or red color in Figure 5), the agreement between GPS and GRACE displacements is good.Larger positive WRMS reduction values indicate better agreement.Whereas for sites where the WRMS reduction values are negative (marked with cool or blue color in Figure 5), the agreement between GPS and GRACE displacements is bad.Larger negative WRMS reduction values indicate worse agreement.The distribution of WRMS reductions is shown in Figure 6 as a histogram.For the 36 GPS sites over the studied European region, we find 22 sites showing positive WRMS reduction val-  ues.For the positive WRMS reductions, many of them (17 of 22 sites) are in the range of 0-20%.While for the negative numbers, they are mainly in the range of −10-0% (10 of 14 sites), which is smaller than the positive.
While van Dam et al. [2007] found only 10 out of 51 sites with reduced WRMS, our results show better agreement between GPS vertical displacements and GRACE loading contribution.Moreover, the inland sites exhibit larger positive WRMS reduction values, which suggests that the areas with more evident hydrological loading effect (see the background color map in Figure 5) have better agreement between GPS and GRACE displacements.This conclusion could further be derived from Figure 7.In Figure 7, we separate 36 sites into 3 groups according to their WRMS reductions, which are larger than 15%, between 0 and 15%, and smaller than 0, and show their GPS and GRACE annual amplitudes and mean amplitude of each group.Evidently, the sites with larger amplitudes of loading signal have better agreement between GPS and GRACE data (larger WRMS reductions).

Restoring amplitudes of GRACE estimates
At high degrees and orders, GRACE Stokes coefficients are dominated by noise, including longitudinal stripes, and other errors.In the present study, a 300 km Gaussian smoothing is applied to GRACE data, in order to suppress the spatial noise in GRACE data.However, applying 300 km Gaussian smoothing will not only reduce most of the short wavelength spatial noise, but also attenuate the amplitude of real signals and meanwhile cause leakage effects [e.g., Chen et al. 2007].Therefore, the GRACE-derived vertical displacements probably underestimate the real deformation.Thus the amplitude of real signal needs to be restored.Landerer and Swenson [2012] elaborated the restoring method of using three scaling factors systematically when estimating surface mass changes using GRACE data.We apply the method here in order to investigate the impact of correcting leakage effects in the study of loading deformation.
First, we adopt the Global Land Data Assimilation System (GLDAS) -NOAH terrestrial water storage (TWS) model data [Rodell et al. 2004], and convert them into Stokes coefficients up to degree and order 180, which approximately correspond to spatial scales of ~110 km (roughly 1 degree at the equator).Next, we calculate the vertical displacements using GLDAS Stokes coefficients with truncation at degree 180, and show the corresponding annual amplitudes in the stud-ied region in Figure 8b.Then, in order to be consistent with GRACE data processing, we truncate the GLDAS Stokes coefficients at degree 60 and apply the 300km Gaussian smoothing.The annual amplitudes of the truncated and smoothed results are shown in Figure 8a.We assume that the amplitude reduction from Figure 8b to Figure 8a represents the attenuation effect.At each grid point, the ratio (i.e., a scaling factor) between Figure 8a and Figure 8b can be computed.Then the scaling factors for the entire region can be applied to restore the true magnitude of the signals.
Assuming that the GRACE data are subject to an equivalent attenuation effect, we can apply the GLDAS model derived scaling factors to vertical displacements predicted by GRACE GSM products (Figure 8c) to obtain the "restored" GRACE results (Figure 8d).After restoration by these scaling factors, we find vertical displacements for the areas in the northwestern Iberian Peninsula, midwestern Balkan Peninsula and around the Black Sea, are significantly recovered in the signal amplitudes.
The annual amplitudes after the restoring for all the 36 GPS sites are listed in Table 1 for comparisons.The differences between the two GRACE amplitudes (with or without restoring), which can easily be computed (column 7 -column 4), are not shown for brevity sake of Table 1.24 of 36 differences are positive, reflecting that the amplitudes have increased after restoring.To some extent, the restoring leads to improved agreement in annual amplitude.We compute the differences between the annual amplitudes of GRACE and GPS estimates before and after restoring (shown in column 8 and 9 of Table 1), and find that 15 of 36 sites have reduced differences after restoring (|column 9|< |column 8|), marked with "y" in column 10.We also show the differences |column 9|−|column 8| of these "y-mark" sites in column 10, which could be used to quantify the effectiveness of the signal restoring method.For example, sites ebre, hers, mate, nya1, nyal, vill have values larger than 10%, which could be regarded as being relatively well restored.Note that the value at pado (also >10%) is unreliable because it is far larger than the others.The reason may be that the leakage effect of GRACE signal from nearby field causes the huge discrepancy between GPS and GRACE data, and also the large percentage.

Conclusions and discussion
Using an extended record of GPS and GRACE data, we reanalyze seasonal vertical displacements due to hydrological loading of 36 GPS sites over Europe, and show improved agreement between GPS and GRACE results than those from previous studies.We find 22 GPS sites show notably reduced WRMS, when GRACE-derived hydrological loading deformations are removed from GPS observations.The mean amplitudes and phases of GPS observations and GRACE-derived loading deformations for the 36 sites agree remarkably well.In addition, GPS observations and GRACE results show good agreement at interannual time scales at a few GPS sites (which could reach up to >0.8 correlation coefficients and >35% WRMS reduction).
The improved agreement can be attributed to a number of reasons.First, we use a longer record of GPS time series with improved data quality.The phasor diagrams (Figure 4a) show a good degree of spatial coherence, especially at the inland GPS sites, which reveals considerable improvement than those in previous studies [van Dam et al. 2007], demonstrating the improved quality of GPS data.The current formal errors of GPS height measurements are about half of those several years ago.The newer version (RL05) GRACE data and the extended record of GRACE time series also help improve the estimated hydrological loading deformations.
It is noticeable that the GRACE results show more spatial coherence than GPS data do, and the inland sites in Europe show better spatial coherence than the coastal sites do.GPS data are in situ point-wise measurements, which tend to be affected by local disturbances, such as artificial effects and small-scale weather impacts.For coastal GPS sites or those located on islands, hydrological loading effects are relatively small.Therefore, GPS observed seasonal vertical deformation signals are considerably smaller in the coastal regions (see Figure 4).Thus the GPS estimates of annual amplitudes and phases in those regions are expected to be with relatively larger uncertainty.The inland sites have higher signal-to-noise-ratio and show better spatial coherence because larger loading effects are expected.
In addition, the annual amplitude of geocenter correction we applied ranges from 0.4 mm to 0.6 mm over the studied region, which is relatively small compared to the seasonal loading signal (a few mm).Given the relatively small magnitude of the geocenter correction, the uncertainty in estimated seasonal geocenter variations [Swenson et al. 2008] is unlikely to be the main error source responsible for the discrepancy between GPS and GRACE results, even though the exact uncertainty of the geocenter correction is presently not clear.
As the GRACE-derived loading deformations are expected to be affected by attenuation due to the truncation and spatial filtering applied to GRACE data, we estimate potential attenuation effects on annual amplitudes of GRACE-derived vertical displacements using the scaling factor method [Landerer and Swenson 2012] based on the TWS changes estimated from GLDAS models.After restoring the GRACE annual amplitudes, only less than half of the GPS sites (i.e., 15 over 36) show improved agreement between GPS and GRACE estimates.Possible reasons are as follows: (a) since the scaling factors are derived from GLDAS-predicted TWS changes, the method is only valid when the model-predicted TWS could perfectly resemble the true TWS change signal for the studied region; Thus any uncertainty in GLDAS TWS predictions will contaminate the estimated scaling factors and affect restored GRACE estimates; (b) the loading deformation signals in western Europe are relatively small (<5 mm in annual amplitude in vertical direction, as shown by the background map of Figure 5), and more likely to be affected by observing noise in the GPS data; (c) a number of studied GPS sites are located along the coastline or nearby the coastal regions, so for these sites the GLDAS-derived scaling factors are more subject to large uncertainty than for those sites of inland regions (due to the complicated leakage effect between land and ocean); (d) the scaling factor method might be less reliable when applied at gridded points (as applied in this study) than at regional scales (for total or average changes of the re-gion); and (e) the JPL GPS products may underestimate the true loading deformation.Despite of the above factors, the restoring experiment at least suggests that non-negligible attenuation effects exist in the study of surface loading deformation using GRACE satellite gravity observations.Furthermore, we have done additional tests with similar GPS products provided by other institutions, such as the Scripps Orbit and Permanent Array Center (SOPAC) GPS solution, and the MEaSUREs combination solution which combines the JPL and SOPAC solutions (also on ftp://sopac-ftp.ucsd.edu/pub/timeseries/measures/ats/).We have found significant discrepancies in the estimated seasonal vertical deformations from different GPS products, which is likely due to the different data calibration strategies and inconsistent definitions of reference frame used in the different solutions.
For the present study we adopt the JPL GPS height solutions in the present study to be consistent with previous studies for comparison sake.
Another factor that possibly affects GPS-estimated annual amplitudes of vertical deformations is the systematic error associated with the draconitic period of 351.2 days in GPS data [Ray et al. 2007, Fu et al. 2015].Accurate separation of this systematic from the true annual deformation signal (~365.25 days) is challenging due to the close periods of the two terms.Therefore, in the present study we do not take into account the potential systematic error, whereas this might affect the comparisons between GPS and GRACE results.
The GRACE mission is entering its 15th year, and the GRACE Follow-On mission is scheduled to be launched in 2017.With the extended data record of GRACE time series (now over 14 years), and improvement of data quality and data processing methods, GRACE timevariable gravity measurements will continue offering great potential to improve understanding of the global mass redistribution and vertical deformation of the Earth surface at a broad band of frequencies.

Figure 2 .
Figure2.Comparisons of monthly GPS height residuals and GRACE-derived vertical displacements at 8 selected GPS sites.Atmospheric and oceanic loading effects have been removed from both GPS height residuals (black curves) and GRACE-derived vertical displacements (red curves), so the GPS and GRACE results here reflect hydrological loading deformation.

Figure 3 .
Figure 3. Comparisons of GPS height residuals and GRACE-derived vertical displacements after 13-month moving average filtering.WRMS reduction in GPS heights is computed after subtracting GRACE-derived displacements.Correlation coefficients are computed at zero phase lags.

Figure 4 .
Figure 4. Annual amplitudes and phases of GPS height residuals (left) and GRACE-derived vertical displacements (right).The annual phase is defined as { in sin [2r(t−t 0 )+{], where t 0 refers to the start of the year (e.g., 2002.0).The length of the arrow represents the annual amplitude of the site, and the direction of the arrow is the peak-value month counted counterclockwise from the east ( January to December).The red arrow on the upper left represents the average annual amplitudes and phases of the 36 selected sites in Europe.The scale of the arrow is marked on the lower left.Background color maps show gridded estimates of annual amplitudes of vertical displacements predicted from GRACE GSM products.A 300 km Gaussian smoothing has been applied to the GRACE data.

Figure 5 .
Figure 5. WRMS reductions when GRACE signals have been removed from GPS heights.Colors of the circles and scales on right represent the percentages of WRMS reduction.Background color map and scale (at the bottom) are the same as in Figure 4.

Figure 8 .
Figure 8. Gridded estimates of annual amplitudes of vertical displacements before and after restoring.The estimates are derived from: (a) GLDAS-NOAH with Stokes coefficients truncation at degree 60 and 300 km Gaussian smoothing; (b) GLDAS-NOAH with Stokes coefficients truncation at degree 180; (c) GRACE GSM products with Stokes coefficients truncation at degree 60 and 300 km Gaussian smoothing; (d) restored annual amplitudes of GRACE-derived vertical displacements.