A Loading Correction Model for GPS Measurements Derived from Multiple-Data Combined Monthly Gravity

: Time-dependent loading deformations of the Earth’s surface, due to nontidal changes in the atmosphere, ocean, land water/ice, etc., contribute signiﬁcantly to the seasonal and secular Global Positioning System (GPS) site displacements, especially for the up component. While loading deformations derived from general circulation model (GCM) outputs are usually used to correct loading signals in the GPS site displacements, this study aims to provide a loading correction model based on the multiple-data combined monthly gravity products LDCmgm90. We have adopted GPS measurements from 249 IGS reference frame stations and 3 different GCM-based loading models to test the reliability of the LDCmgm90 model. Compared to the GCM-based models, the LDCmgm90 loading correction is more effective in attenuating seasonal (especially annual) loading signals and can bring more signiﬁcant improvements to most stations for both the data-trend-removed and the data-trend-retained cases. Thus, we have validated the LDCmgm90 model from the loading aspect and proved it to be a reliable loading-correction model for GPS displacements. The relatively better secular loading signals provided by the LDCmgm90 loading model may provide us a chance to study the long-term, nonloading signals in GPS data.


Introduction
Despite successful applications in various fields, Global Positioning System (GPS) measurements suffer from scatters in GPS position time series, which reduce GPS measurement accuracy significantly. Ray [1] listed many sources of measurement error and other systematic effects, such as draconitic errors, local site errors, temperature cycle errors and errors in IERS models [2] (also see Ray et al. [3,4]; more details can be found in Appendix B). Besides these causes, time-dependent loading deformations, due to mass fluctuations in the Earth's surficial geophysical fluids, i.e., atmosphere, ocean and land water/ice, lead to variable loads on the crust and thus variable crustal deformations, thus also contributing greatly to GPS position data scatters, especially at seasonal frequency bands [5][6][7][8][9][10][11][12][13][14][15][16]. Therefore, accurate loading deformation corrections are urgently needed in GPS data analyses, or else the explorations of the geophysical causes of the remaining crustal deformations will be affected. Some studies on loading are reviewed below.
Van Dam and Wahr [17] studied the displacement of the Earth's surface due to atmospheric loading. Mangiarotti and Cazenave [18] compared loading results derived from 16 Doppler orbitography and radiopositioning integrated by satellite (DORIS) stations and those from annual surface mass redistributions (atmosphere and ocean mass, soil moisture and snow loadings). Van Dam et al. [5] compared vertical displacements modeled by general circulation models (GCMs) and GPS heights observed by 147 globally distributed sites (∆r M and ∆r O ) and found that the root mean square (RMS) of seasonal, especially annual, ∆r O adjusted by ∆r M was significantly reduced. Dong et al. [8] conducted a comprehensive study on the effectiveness of loading corrections and showed that ∼40% of the power of the annual vertical variations in continuous GPS time series can be explained by the joint contribution of seasonal surficial mass redistributions. Blewitt [19] discussed the self-consistency in reference frames, geocenter definition and surface loading of the solid Earth. Since these studies, GCM-derived loadings were widely used to correct the GPS measurements. However, with longer and more accurate time series of GPS measures and GCM outputs available, it became obvious that GCM-derived loading corrections do not always improve the GPS data, but sometimes lead to increased noises, partly due to limitations in the GCMs, such as hydrostatic approximation in atmospheric and oceanic models, and violation of the conservation of global mass when summing atmospheric, oceanic and hydrological masses (see Chen et al. [20,21] for more details). Therefore, the GCM-based loadings may have biases in both phase and amplitude and sometimes lead to even more excess seasonal signals (worse results) in GPS time series (in Section 4, one can see that seasonal noises in GPS data are sometimes amplified rather than attenuated after correction by GCM-derived loadings).
The 2002 launch of the Gravity Recovery and Climate Experiment (GRACE) twin satellites has made available several time series of geopotential coefficients [22], which can provide information for surficial mass redistributions and therefore loading deformations meeting the conservation of global mass. Notable discrepancies are found between GCMbased and GRACE-based loading deformation models. Davis et al. [23] compared the annual deformations measured by GPS stations in the Amazon River Basin with predictions calculated from GRACE measurements and found high correlation. They also confirmed the superiority of GRACE-measured surface water variations over some hydrology models, indicating that the adoption of GRACE measurements can be a better approach to study loading deformation compared with GCMs. Tregoning et al. [11] computed elastic deformation using continental water loading estimates derived from GRACE, which agreed well with GPS observations. Yan et al. [24] proved that using either GRACE or GCM 36 degree/order (d/o) spherical harmonic coefficients can explain 90% of the variance of total crustal vertical deformation. Fu and Freymueller [25] and Chanard et al. [26] confirmed that surficial mass redistributions observed by GRACE, combined with an elastic spherical and layered Earth model, can be used to provide first-order corrections for loading deformation observed in both horizontal and vertical components of GNSS station position time series.
Despite the advantages of GRACE-derived loading, there are also some significant defects in the GRACE-based loading corrections. Due to the longitudinal striped errors and coarse spatial resolution of traditional GRACE spherical harmonic products, additional corrections (destriping and smoothing) are required to reveal valuable signals [27,28]. Consequently, with different GRACE data processing strategies adopted, the loading deformation at the same GPS station derived by different researchers can differ from each other. Further, destriping would cause distortions of signals of interest while smoothing would attenuate them (signal leakages). These problems, together with limited spatial and temporal resolutions, make the GRACE-based results not quite reliable. The newly released Mascon products are free of the longitudinal striped errors but also contain some problems. For example, atmospheric, oceanic and hydrological model outputs (used to generate the RL06 AOD1B) violate the conservation of global mass, so a sea-level mass term is introduced to conserve the mass. However, only the ocean model outputs are used to produce the RL06 GAB (namely, parts of oceanic contributions are not modeled). In addition, like traditional GRACE spherical harmonic products, there are notable differences among Mascon products released by different institutes (more details can be found in Chen et al. [29]). Thus, further improvements are needed for GRACE products.
By inverting GPS displacement series measured at roughly 450 continuously tracking sites and ocean bottom pressure estimates of a data assimilated ocean circulation model, Wu et al. [30] obtained monthly global surface mass distributions in the spherical harmonic domain with a complete spectrum up to degree and order 50 and then combined them with GRACE gravity data to provide enhanced spectral and geographic coverage. Chen et al. [29] developed a multiple-data-based monthly geopotential model set LDCmgm90 (available at https://doi.org/10.6084/m9.figshare.7874384.v4 (accessed on 20 September 2021)) by assimilating various versions of GRACE (including Mascon products); satellite laser ranging (SLR) monthly gravity data [31,32]; and outputs from various global atmospheric, oceanic and hydrological circulation models by using the least difference combination method [20,21]. The LDCmgm90 dataset provides more reliable data trends (see [21,29] for some validations) and has no stripes as seen in most GRACE data, and it was shown to have better performance in explaining atmospheric, oceanic and hydrological long-term and seasonal time-scale mass redistributions compared with most meteorological models and time-varying gravity products, including Mascon solutions.
This study aims to provide a particular choice of improved surficial loading models by using the LDCmgm90 dataset and, in the meantime, to validate the LDCmgm90 from the loading aspect. We test our LDCmgm90 loading model at 249 IGb14 GPS reference frame stations and compare them with three GCM-based loading models respectively provided by the International Mass Loading Service (IMLS, http://massloading.net/, accessed on 27 June 2020), the School and Observatory of Earth Sciences (EOST) Loading Service (http://loading.u-strasbg.fr/index.php, accessed on 8 December 2020) and the Earth System Modeling group at the Deutsches GeoForschungsZentrum (ESMGFZ) (http://rz-vm115.gfz-potsdam.de:8080/repository, accessed on 6 January 2021). The results show that the LDCmgm90-based loading model explains more vertical crustal deformation caused by surficial loading at most GPS stations than other loading models involved in this research. This paper is arranged as follows: Datasets and methods used in this research are described in Section 2. LDCmgm90-based load deformation is derived in Section 3 and compared with GCM-based load deformations at 249 IGS GPS stations in Section 4. Conclusions are presented in Section 5.

GPS Data
The IGb14 reference frame is an update to the previous IGS14 reference frame in which a large fraction of the IGS14 GPS station coordinates have been "repaired" after previously being rendered unusable due to the discontinuities caused by unpredictable events such as earthquakes or equipment changes. The transformation parameters between IGb14 and IGS14/ITRF2014 are therefore all zero because IGb14 is aligned in origin, scale and orientation to IGS14, hence ITRF2014 (please refer to ftp://igs-rf.ign.fr/pub/IGb14/IGb14. ssc (accessed on 20 September 2021) for more details about IGb14 reference frame and the 261 GPS stations involved). However, only 259 of 261 IGb14 GPS stations' 24-h-sampled final solutions are released on the website of Nevada Geodetic Laboratory (NGL).
NGL-processed GPS data contain time series of east-, north-and up-components of nonlinear displacements in the IGS14 reference frame, where multiple measurement models and corrections have been applied [33][34][35][36][37][38][39][40], including ocean tidal loading correction. However, nontidal loading correction of atmosphere, ocean and hydrology is not applied. The NGL modeling of tropospheric delays is particularly notable and effective at retaining the full nontidal atmospheric load signal in the observed station coordinate time series [41]. Please refer to NGL GPS Data Analysis Strategy and Products Summary (http://geodesy. unr.edu/gps/ngl.acn.txt, accessed on 6 December 2020) for more details.
We first abandoned 10 GPS stations whose durations are less than 2 years (730 samples) within the LDCmgm90 timespan due to the unreliability of their annual fluctuations. Then, we preprocessed the 249 remaining time series to remove their gross errors (larger than 3 times standard deviation) and used a first-order difference method to correct jumps in GPS time series, where iteration and decreasing thresholds are required. The distribution of the 249 selected GPS stations is shown in Figure 1.

GCM-Based Loading Data
The IMLS [42,43] provides 3-h-sampled nontidal loading deformation models for 1272 stations, containing GEOS-FPIT atmospheric and hydrological models developed by Global Modeling and Assimilation Office at NASA Goddard Space Flight Center and MPIOM06 oceanic model developed at ESMGFZ. This loading model set is termed as the IMLS loading for short.
The EOST Loading Service [42] provides 1-h-sampled loading deformation datasets (in the form of Stokes coefficients up to degree and order 120; available at http://loading. u-strasbg.fr/listdata.php?dirn=st, accessed on 8 December 2020) based on the European Centre for Medium-Range Weather Forecasts (ECMWF) fifth reanalysis (ERA5) atmospheric surface pressure fields, induced TUGO-m ocean model [44] and GLDAS/Noah continental hydrology model [45]. However, unlike other classical ocean general circulation models, TUGO-m model is forced not only by atmospheric pressure but also by surface winds (please notice that both atmospheric and oceanic loadings are included in their ERA5_TUGO product while ERA5_IB contains only atmospheric loading). This loading model set is termed as the EOST loading for short.
Both the IMLS and EOST loading data are provided in the center of mass (CM) frame without corresponding data for geocenter motion corrections. We converted them to the center of figure frame to match the NGL GPS data, using geocenter motion derived from the LDCmgm90 model as described in Section 2.3 and Appendix A.
The ESMGFZ provides grid data for nontidal elastic surface loading deformation in the center of figure (CF) frame with 0.5 • × 0.5 • spatial resolution [46]. ESMGFZ nontidal atmospheric loading is 3-h-sampled and based on the operational ECMWF model which is typically updated about twice a year to incorporate advanced algorithms and physical models into the system, which affects its long-term stability but improves near real-time and forecasting quality, leading to substantially higher spatial resolution compared with alternatively available atmospheric reanalysis models. ESMGFZ nontidal oceanic loading is 3-h-sampled and based on Max Planck Institute for Meteorology Ocean Model (MPIOM) ocean model and forced by the ECMWF atmospheric pressure. ESMGFZ hydrological loading is 24-h-sampled and based on the version-2 Land Surface Discharge Model (LSDM) hydrological model, in which soil moisture, snow, surface water and water in rivers and lakes are all considered with mean mass during 2013~2014 subtracted. High-resolution geographic information system (GIS)-based river network is used in this model for higher accuracy of mass. Unlike other loading datasets, ESMGFZ also provides a special 24-hsampled sea-level loading for the conservation of global mass. Therefore, when ESMGFZ loading is used, the sum of the above four components represents the complete surface loading deformation. This loading model set is termed as the ESMGFZ loading for short.
The above-mentioned model sets are briefly summarized in Table 1. In this study, all the GCM-based loading data were downsampled to monthly time series (with signals with frequencies not exceeding 6 cycles per year (cpy)) to match the monthly LDCmgm90 loading. Table 1. Atmospheric, oceanic and hydrological nontidal loading based on different meteorological models or products involved in four loading datasets used in this research.

Dataset
Atmosphere Loading

Ocean Loading
Hydrology Loading

LDCmgm90-Based Loading Data and Geocenter Motion Correction
In this research, LDCmgm90 time-variable gravity [29] was adopted to derive LD-Cmgm90 global deformation model. The LDCmgm90 GSM dataset provides Stokes coefficients for degree/order 2~90, so we needed to determine the corresponding degree-1 coefficients ∆C 10 , ∆C 11 and ∆S 11 in this study.
The GRACE monthly gravity products are usually released together with the GRACE AOD1B products, which provide a model-based dataset (including GAA, GAB, GAC and GAD) that describes the time variations of the gravity potential at satellite altitudes that are caused by nontidal mass variability in the atmosphere and oceans [47]. The GAA product describes the monthly nontidal atmospheric mass anomalies simulated by the operational run of the atmosphere model ECMWF. GAB refers to monthly nontidal oceanic mass anomalies simulated by the operational run of the (unconstrained) ocean model MPIOM. GAC is the sum of GAA and GAB. GAD can be regarded as a revised version of GAC with nontidal atmospheric and oceanic mass anomalies only over ocean areas. GSM is just the gravity residual after GAA and GAB are removed from the GRACE observations. Like the monthly GRACE gravity models, the LDCmgm90 dataset also has five components, namely the GSM, GAA, GAB, GAC and GAD products [29]. However, unlike common GRACE solutions, destripe correction is not used in processing LDCmgm90 solutions, which do not contain the stripe errors.
Recognizing the limited accuracies of glacial isostatic adjustment (GIA) models, Chen et al. [29] did not include the GIA correction into the LDCmgm90 dataset. The GIA effects reflect the viscoelastic behavior of the Earth's mantle, while loading deformations are only relevant with the Earth's elastic responses. Thus, we needed to remove the GIA contributions from the LDCmgm90 monthly gravity, and in this research, we chose the ICE6Gavg-GPS model (available at https://www.mdpi.com/2072-4292/12/7/1209/s1, accessed on 2 February 2020) [48], which combines multiple GIA models and global GPS observations, to remove GIA signals in the LDCmgm90 dataset. The loading model set derived from the LDCmgm90 dataset, after GIA correction, is termed as the LDCmgm90 loading for short.
The 3D NGL GPS displacements are provided in the CF frame, while the LDCmgm90 and all the relevant GRACE products are given in the CM frame. Therefore, we needed to convert the LDCmgm90-derived loading results from the CM frame to the CF frame. The CM and CF move relative to each other due to the degree-1 surficial mass loading. Let ∆X(t), ∆Y(t) and ∆Z(t) be three components of geocenter variations, reflecting differences between the Earth's CM and CF. We have [30,49] In Equations (1) and (2), ρ ave = 5517 kg/m 3 is the average density of the Earth, a = 6378 km is the semimajor axis of the Earth, ∆σ is the incremental surface density and can be derived from Equation (14) of Wahr et al. [27] and k l is the degree-l elastic loading Love number for geopotential variation. In this study, we used degree-1 loading number k 1 = 0.027 derived by Han and Wahr [50] since the origin of the coordinate system is the center of figure (also see Wahr et al. [27] for more details).
To derive ∆C 10 , ∆C 11 and ∆S 11 corresponding to a global time-dependent gravity product, one possible approach is using an independent ocean bottom pressure model that is the ocean function that equals 1 in ocean areas and 0 in land areas. Another way, which was utilized in this research, is using the corresponding GRACE GAD product (please notice that GRACE GAA, GAB, GAC and GAD products are derived from meteorological model instead of GRACE observation). According to the definition, degree-1 Stokes coefficients derived from the GAD product are exactly those defined by Equation (2). Please notice that in Equation (3), GAC and GSM products should be both considered and transformed into mass spherical harmonic coefficients to present the complete global effect.

Method to Determine Loading Deformations
Loading deformations at colatitude θ and longitude φ caused by mass redistribution on the surface of the Earth presented by fully normalized Stokes coefficients can be derived as [51] where the ∆R(θ, φ, t), ∆E(θ, φ, t) and ∆N(θ, φ, t) are respectively the radial, eastward and northward loading deformations; g ave = 9.826057 m/s 2 is the average gravitational acceleration on the surface of the Earth; and h l and L l are the degree-l vertical and horizontal elastic loading Love numbers, respectively. For degree > 2, we adopt k l , h l and L l derived from the MATLAB code developed by Chen et al. [52]. V l (θ, φ, t) is degree-l gravitational potential on the surface of the Earth and can be expressed as where G = 6.67259 × 10 −11 N· m 2 /kg 2 is the gravitational constant and M = 5.975562 × 10 24 kg is the total mass of the Earth.

Results
By using the derived degree-1 Stokes coefficients (please refer to Appendix A for more details) and degree > 1 Stokes coefficients of the LDCmgm90 dataset, we derive the corresponding horizontal and radial loading deformation separately after removing GIA effects. Obvious annual fluctuation is observed in degree-1 radial loading deformation in 2010 (Figure 2), which is similar in other years. Global distribution of degree > 1 loading deformation is significantly more complex due to the contributions of time-dependent precipitation, discharge of hydrology and atmospheric and oceanic mass redistribution.
Radial loading deformation is significantly larger than horizontal loading deformation. Globally, in 2010, absolute values of radial loading deformation in middle and northern high latitudes remain less than 15 mm, while Antarctica and the tropical continents of South America, central Africa, southern Asia and northern Australia have significant annual loading deformation variations with amplitudes larger than 15 mm (see Figure 3, line 1, GAC + GSM).
To separate the contributions of different sources, we derive the global radial loading deformation of LDCmgm90 GSM and GAC products separately. We find that continental anomalies along the equator are only observed in GSM loading deformation ( Figure 3, Line 2, GSM). These prominent annual variations reflect the changes of terrestrial water storage in tropical continents of both hemispheres since precipitation brings a large amount of water in their respective rainy seasons, leading to downward loading deformation, while high temperature and drought during dry seasons vaporize water, leading to upward loading deformation. In middle-and high-latitude areas, annual variations of GSM loading deformation are not so obvious due to less precipitation and lower temperature compared with the tropical area. Seasonal melting and freezing in high latitudes also contribute a little to GSM loading deformation, among which the Greenland coastline is the most prominent. However, annual loading deformation variations in middle and high latitudes are mainly caused by changes in atmospheric pressure and ocean bottom pressure, especially in Antarctica.  (Figures 4 and 5). On the whole, loading deformation diverges from a positive anomaly and converges to a negative anomaly. The direction of horizontal loading deformation is the same as the gradient of the vertical loading deformation field.

Discussion
We applied our LDCmgm90 GAC-and GSM-derived loading products to 249 IGb14 GPS time series processed by NGL to correct nontidal loading effects. Likewise, EOST and ESMGFZ atmospheric, oceanic and hydrological loading products were applied for comparison. Since these stations are scattered globally with diverse data characteristics and qualities, the results of loading corrections will reflect the reliability of load deformation models adopted. We also selected 12 of 249 GPS stations as examples, in which IMLS loading products were also applied.

Comparison of Annual Residual
The up-component displacements of 12 IGb14 GPS stations are compared in Figure 6, in which loading effects have been removed by the four different loading products mentioned above (only in Figure 6 the trends are removed for easier comparisons of the seasonal signals). In most cases in the time domain, according to the remaining residuals, we find that more loading effects are removed by EOST, ESMGFZ and LDCmgm90 loadings compared with IMLS loading, especially at NRIL and WSRT. According to the corresponding residuals' power spectrum density (PSD) shown in Figure 7 and annual amplitudes of time series listed in Table 2, we find that, among these 12 stations, significantly more annual loading is removed by LDCmgm90 datasets at 6 stations: ADE1, BREW, CHPI, MOBS, TNML and WSRT. The EOST loading is better at three stations: AMC2, NKLG and TIXI. ESMGFZ loading is better at three stations: HOLM, LHAZ and NRIL. No advantage of IMLS loading is found among these stations. In fact, applying IMLS loading often increases the annual amplitudes of the 12 residuals.  We also notice that at several stations, applying loading corrections leads to excess annual amplitudes instead of attenuating them, which indicates that worse results are induced by loading models. According to Table 2, compared with original GPS time series, annual loading amplitudes are successfully reduced at 6 stations with IMLS loading applied, at all 12 stations with EOST loading applied, at 10 stations with ESMGFZ loading applied and at 11 stations with LDCmgm90 loading models applied. Numerically, among these 12 stations, the most prominent improvement is brought by ESMGFZ loading at HOLM, where 91% up-component annual fluctuation in the original GPS signal is removed (only 0.21 mm remains).
The accuracy of hydrological loading is not as reliable as that of atmospheric and oceanic loadings due to the uneven distribution of continent water storage and much less accurate monitoring generally. Different hydrological models diverge from each other and may enlarge the residuals instead of attenuating them further. In Table 2, we notice that at eight stations (ADE1, AMC2, HOLM, MOBS, NKLG, NRIL, TIXI and WSRT), applying IMLS hydrological loading leads to larger annual residuals than applying IMLS atmospheric and oceanic loadings only, which is significantly improved by EOST, ESMGFZ and LDCmgm90 loading corrections. Therefore, we infer that EOST, ESMGFZ and LDCmgm90 models contain better hydrological loading information compared with IMLS loading.
Divergences between different loading models are smaller in other frequency bands at these 12 stations (such as semiannual loading, see frequency = 2 cpy in Figure 7 and Table 3). However, since loading effects are mainly annual fluctuations [5,8], discrepancies among loading models are mainly introduced by their annual signals, which should represent the quality of loading models. Therefore, due to the smaller time-and frequency-domain annual residual and better hydrological loading information involved, we believe that EOST, ESMGFZ and LDCmgm90 loadings are generally better than IMLS loading. Table 2. Effects of applying loading corrections: the case of annual signals (unit: mm).  Table 3. Effects of applying loading corrections: the case of semiannual signals (unit: mm).

Loading Due to Seasonal Mass Redistributions
We removed the means and trends of all the relevant data in order to estimate the improvement of seasonal loading correction brought by our LDCmgm90 loading models compared with EOST and ESMGFZ datasets. First of all, we calculated the mean PSD of 249 IGb14 GPS time series with different loading correction models applied (see Figure 8).
One can see the LDCmgm90 is the most efficient in attenuating the annual signal in GPS data, especially in the up-component. Besides, long-period variations of up-component GPS time series are also improved by LDCmgm90 compared with EOST and ESMGFZ. However, there are also strong peaks of draconitic year error (with period 351.2 day or frequency 1.04 cpy; also see Ray et al. [3]) and its higher harmonics (peaks near the 2, 3, 4 and 5 cpy) as shown by Figure 8. These draconitic year errors in GPS data have nothing to do with loading but also contribute greatly to the scatter of GPS position residual, which can be well corrected by least-square fitting (see three bottom subplots in Figure 8). On the other hand, both the EOST and ESMGFZ have better performances than LDCmgm90 for up-component variations between draconitic peaks, the cause of which is not quite certain but could be a limitation of the frequency and spatial resolutions (~6 cpy and~2 degrees, respectively) of the LDCmgm90 data (it is rather common that there is no data point at the frequency of interest and thus the plots are impacted by its neighboring frequency points, due to the limited frequency domain sampling rate caused by the finite length of the time series of interest). Root mean square (RMS) of detrended GPS residual is a vital index to estimate the quality of a loading model (the mean value of each time series is removed before time-domain analyses). Globally, compared with applying EOST and ESMGFZ loadings, applying LDCmgm90 loading improves the residual RMS at more of the 249 IGb14 GPS stations involved in this research. Up-component RMS variation (∆RMS) brought by LDCmgm90 loading compared with applying either EOST loading or ESMGFZ loading is shown in Figures 9 and 10, respectively. Besides, statistical information about RMS improvement ratio is listed in Tables 4 and 5.
Compared with original GPS signal (Table 4), applying LDCmgm90 loading successfully reduces up-component RMS at 218 (equaling 87.55%) GPS stations, which is slightly more than applying EOST or ESMGFZ loading. When significantly improved (RMS improvement ratio > 5%) stations are considered only, LDCmgm90 loading is about 4% better than EOST and matches ESMGFZ loading.
In horizontal components, LDCmgm90 roughly matches EOST and ESMGFZ; differences between RMS with either loading model applied at most stations are within 5%. Numerically, according to    1 Positive RMS improvement ratio means that RMS of LDCmgm90 residual is smaller than that of EOST or ESMGFZ residual, indicating that LDCmgm90 is better, while a positive value means that EOST or ESMGFZ is better. 2 "Station Improved" column shows the total number and ratio of GPS stations whose loading effects are better removed by LDCmgm90 compared with ERA5 or ESMGFZ, equaling ">5%" column plus "0~5%" column. For each component, the total number of improved stations is listed outside the brackets while the corresponding ratio is given in the brackets. 3 RMS improvement ratio brought by LDCmgm90 compared with EOST is calculated by 100%-RMS (GPS-LDCmgm90)/RMS (GPS-EOST), which is similar to that compared with ESMGFZ. Decrease in RMS means that applying LDCmgm90 loading reduces the residual's RMS compared with applying EOST loading, indicating that LDCmgm90 loading is better than EOST at these GPS stations, while increase in RMS means that EOST is better. Besides, locations of 249 stations are marked with "+" (blue "+": >5% RMS is removed; cyan "+": 0~5% RMS is removed; magenta "+": 0~5% RMS is induced; red "+": >5% RMS is induced).
To better understand these three loading models, we also compared EOST, ESMGFZ and LDCmgm90 together, finding the best loading model for each station involved in this research and calculating the RMS improvement ratio of GPS residuals brought by LDCmgm90 compared with the better result of applying EOST and ESMGFZ loadings ( Figure 11). To reach the smallest RMS of GPS residuals, applying EOST is the best choice for 54 (21.69%) stations, applying ESMGFZ is the best choice for 85 (34.14%) stations and applying LDCmgm90 is the best choice for 110 (44.18%) stations. Regionally, EOST-best stations lie mainly in the middle-latitude area of North America and Asia. ESMGFZ-best stations are mainly in western Europe. LDCmgm90-best stations are widely distributed around the world, especially along the coastline of Greenland.
Therefore, according to the comparison of 249 globally distributed IGb14 GPS stations, especially the improvement in removing up-component effects, we conclude that LD-Cmgm90 loading is reliable and generally better than EOST loading and ESMGFZ loading. Decrease in RMS means that applying LDCmgm90 loading reduces the residual's RMS compared with applying ESMGFZ loading, indicating that LDCmgm90 loading is better than EOST at these GPS stations, while increase in RMS means that ESMGFZ is better. Besides, locations of 249 stations are marked with "+" (blue "+": >5% RMS is removed; cyan "+": 0~5% RMS is removed; magenta "+": 0~5% RMS is induced; red "+": >5% RMS is induced).  Figures 9 and 10, for each station, red "+" means EOST loading is the best, magenta "+" means ESMGFZ is the best and blue "+" means LDCmgm90 is the best.

Contributions of Secular or Long-Term Mass Redistributions
Loading affects not only seasonal variations but also trends in GPS measurements, which should be appropriately considered by a complete loading model (note that many long-period geophysical signals, e.g., the~18.6-year tide, may look like trends for the relatively short GPS time series; thus, "trends" here may contain contributions from these long-period geophysical signals). Therefore, to test the secular or long-term signals of LDCmgm90 loading, we only removed the means but retained the trends of all data, and we applied GIA corrections (ICE6Gavg-GPS model) to GPS data and LDCmgm90 loading (termed "trend-retained" for short hereafter). The analyses are similar to those in Section 4.2. Comparisons are listed in Tables 6 and 7.
According to the comparison of 249 GPS stations in Table 6, when trends are considered, the number of stations significantly improved (RMS improvement ratio > 5%) by applying EOST and ESMGFZ loadings decreases by 15~20%, while the number of stations significantly improved by applying LDCmgm90 loading varies only 2%, which is 16% better than EOST loading and 18% better than ESMGFZ loading.
Applying LDCmgm90 loading successfully reduces up-component RMS at 201 (equaling 80.72%) GPS stations, which is similar to applying EOST loading and 12% more than applying ESMGFZ loading. On the other hand, the~14% larger improvement of the EOST model with respect to the ESMGFZ model is not found in trend-removed comparison. This may be due to the facts that the former adopts the ECMWF reanalysis (ERA5) of atmospheric data and the associated oceanic and hydrological models are forced by both surface pressure and winds (more practical), which induces better trend signals, while the ESMGFZ only uses ECMWF operational data and its associated ocean model is forced only by surface pressure (less practical; operational data usually have poorer quality than reanalysis data). Improvement brought by the trend of LDCmgm90 loading can also be observed in Table 7, where significantly improved (RMS improvement ratio > 5%) stations rise 11% and 15% when compared with EOST and ESMGFZ results, respectively. Total numbers of improved stations also slightly increase. We also repeat the three-model comparison similar to Figure 11, the result of which is given in Figure 12. We find that the largest RMS improvement ratio is about 80%, which is nearly double that in Figure 11. Besides, when trends are considered, applying LDCmgm90 is the best choice for 116 (46.59%) stations, slightly more than trend-removed condition. In order to show the advantages of long-term and annual signals in LDCmgm90 loading more clearly, we calculate the improvement of LDCmgm90 loading in long-term (with periods > 5 years) and annual frequency bands (with periods between 338 and 410 days) compared with EOST and ESMGFZ loadings (Table 8). We find that 49% EOST and 49.8% ESMGFZ long-term residuals and 57.83% EOST and 55.02% ESMGFZ annual residuals are significantly improved, larger than 43.37% and 44.58% in Table 7 where all <6 cpy signals (period > 2 months) are considered, supporting the contributions of LDCmgm90 loading in long-term and annual bands. In the meanwhile, total numbers of improved stations in Table 8 remain similar compared with Table 7, indicating that (1) general advantages of LDCmgm90 loading still exist and (2) EOST and ESMGFZ loading models also perform well in other seasonal variations, which corresponds with PSD shown in Figure 8.
Therefore, according to the trend-retained comparison in this section and band comparison in Table 8, we confirm that not only seasonal variations but also long-term signals in LDCmgm90 loading are better than those in EOST and ESMGFZ loadings.   Figure 11, but for trendretained comparison.

Conclusions
By using LDCmgm90 multiple-data-based monthly gravity dataset and estimating its degree-1 Stokes coefficients, we have obtained the complete LDCmgm90 atmospheric, oceanic and hydrological loading models. Globally, we compared LDCmgm90 loading model with GCM-based loading models provided by the EOST Loading Service and ESMGFZ by applying them to 249 globally distributed IGb14 GPS stations processed by NGL to remove their loading effects. Loading models provided by IMLS were also applied to 12 GPS stations.
Our research shows that the LDCmgm90 loading model can reliably remove loading effects in GPS time series and is likely to achieve similar or even better results compared with other latest loading models. In general, according to our numerical results, we infer that the EOST and ESMGFZ loading models are better than the IMLS loading, but the LDCmgm90 loading is even better, not only for GPS sites on the continents but also for those on the islands surrounded by the oceans. The better annual, trend and longperiod information, especially that of continent water, contribute to the superiority of LDCmgm90 loading compared with other loading products and therefore lead to significant improvements in the LDCmgm90 model with respect to the EOST and ESMGFZ models.
When all the data trends are removed, applying LDCmgm90 loading successfully attenuates up-component loading effects at 218 (87.55%) GPS stations. In addition, timeand frequency-domain analyses show that LDCmgm90 has better performance than EOST at 155 (62.25%) GPS stations and ESMGFZ at 142 (57.03%) GPS stations. These values increase to 160 (64.26%) for EOST and 152 (61.04%) for ESMGFZ if all the data trends are retained, implying that LDCmgm90 can provide better secular or long-term global mass redistributions. Therefore, the LDCmgm90 loading model is the best choice to remove not only seasonal but also secular or long-term loading signals in GPS site displacement data, providing a chance to explore other geophysical causes buried in the residual secular or long-term signals in the GPS data. The LDCmgm90 loading model is available upon request.
It should be mentioned that there are also some limitations in the LDCmgm90 loading model, such as limited spatial (~2-degree) and frequency (~6-cpy) resolutions and availability limited to the GRACE mission period (~2002-2017, but the LDCmgm90 will be extended to include the GRACE Follow-On data in the future). While the~6-cpy frequency resolution should be enough for global mass loading, the~2-degree spatial resolution could be insufficient for some small regions with large local hydrological signals, such as California. Thus, no load model is perfect; each contains certain defects and even errors. Users should beware of applying load corrections to GPS time series data blindly assuming that no new errors will be introduced.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Estimations of Degree-1 Stokes Coefficients and Geocenter Variation of LDCmgm90
We estimate the LDCmgm90 degree-1 Stokes coefficients (∆C ) and compare them with degree-1 coefficients of CSR and JPL RL06 GRACE Mascon solutions ( Figure A1). To recover the corresponding total global effects of the original Mascon solution, we remove GAD products, reconsider the ICE-6G-D model and restore GAC products (total Mascon = original Mascon product-GAD + GIA + GAC). Our results, especially ∆S 11 , have good agreement with other degree-1 coefficients. The annual amplitudes and phases of LDCmgm90 degree-1 Stokes coefficients are compared with values from RL06 Mascon solutions from CSR and JPL in Table A1. Ycomponent annual fluctuation of LDCmgm90 total geocenter variation has good agreement with CSR and JPL Mascon solutions. However, phases of X-and Z-components of LD-Cmgm90 total geocenter variation are slightly larger than others. Our X-component trend also deviates from CSR and JPL Mascon solutions, while X-component and Z-component are similar to CSR and JPL Mascon solutions separately. Thermal cycle of the Earth will cause expansions of antenna structure and surrounding bedrock [13].

• Errors in IERS Models
The recommended IERS Conventions (2010) models suffer weaknesses in various respects, including deficiencies and errors in some models as well as a lack of any standard model for some effects. Models with known errors include ocean tidal loading [16], diurnal and semidiurnal EOP tides that cause aliases into longer-period signals [53] though an alternative model has been put forward outside the formal conventions (https://ivscc.gsfc. nasa.gov/hfeop_wg/, accessed on 27 October 2021), possible lateral heterogeneity effects in the solid Earth tidal displacement [54] and the atmospheric pressure loading (e.g., S1, S2). Time variation of the low-degree geopotential terms has been recognized, but no standard model is recommended for all techniques. In fact, there have been no major updates to any IERS models in recent years apart from the return to a secular mean pole model in 2018 that replaced the ill-considered quadratic model recommended for several years. As a consequence, the techniques and analysis centers have adopted divergent approaches in some cases.

Antenna-Related Effects
Local multipath errors and antenna calibration errors may add aliasing signals to observation records. Snow, ice cover and rain may also affect antennas by adding annual signals or changing sky visibility, especially in high latitudes [4] Clear differences between products for the same stations released by ACs sometimes exist [4].