Comparative Analysis of The Correction Effect of Different Environmental Loading Products On Global GNSS Coordinate Time Series

The global navigation satellite system (GNSS) coordinate time series is affected by the environmental loading (including atmospheric loading (ATML), hydrological loading (HYDL), non-tidal oceanic loading (NTOL), etc.) and many organizations now provide grid products of these loadings. The temporal and spatial resolutions of these products, the loading models and data sources used are not the same, so the effect of correcting the nonlinear deformation of the GNSS coordinate time series is obviously different. This study mainly selects the three agencies, namely, School and Observatory of Earth Sciences (EOST) in France, German Research Center for Geosciences (GFZ) in Germany, and International Mass Loading Service (IMLS) in the United States, including 6 types of ATML models, 7 types of HYDL models and 5 NTOL models. The classication of these 18 environmental loading models was discussed, and the root mean square (RMS) reduction rate of the GNSS coordinate time series after environmental loading corrections (ELCs) was used to evaluate the performance differences of various models. Our results show that both the different models provided by the same organization and the same model provided by different organizations have different correction effects. Regardless of the models, it has a signicant impact on the vertical coordinate time series. In order to correct the nonlinear deformation of the GNSS stations to the greatest extent, based on the above analysis, this study selects the optimal model combination of three environmental loadings as ECMWF_IB+MERRA2+ECCO1, and then explores its inuence on the periodic signals in the GNSS coordinate time series. Research suggests that environmental loadings have a signicant impact on the amplitude and phase of GNSS time series. Especially in the vertical direction, the largest RMS value can reach 8.42 mm. Before and after ELCs, the maximal difference of the annual amplitude and the half-annual amplitude at global 631 stations can reach 8.96 mm and 1.51 mm, respectively. Among them, 84.60% of the stations were corrected by the optimal environmental loading combination model, thus the nonlinear deformation was weakened.


Introduction
In the past 20 years, many scholars have used GNSS coordinate time series to carry out different aspects of geoscience research work, mainly including periodic signal analysis in GNSS coordinate time series  (Farrell, 1972;. Existing studies have shown that environmental loading factors including ATML, NTOL and HYDL are the main cause of nonlinear changes in GNSS time series (Jiang et al., 2010;Sun et al., 1995). By establishing the surface loading model and calculating the elastic deformation displacement sequence of GNSS stations, the geophysical in uence mechanism at the station position can be effectively studied. At the same time, it compares and analyzes the characteristic changes of nonlinear deformation in the GNSS coordinate time series before and after the surface loading correction, and detects the effects of various surface loadings on the GNSS coordinate time series. This is essential to correctly understand the geophysical sources of nonlinear changes in GNSS coordinate time series.
For ATML model, van Dam et al. (1994) analyzed the elevation time series of 19 global positioning system (GPS) stations in the northern hemisphere. They concluded that when the atmospheric loading correction was applied, the variance of the GPS elevation time series could be reduced by 24%, and the atmospheric loading displacement series uctuates greatly in high latitude regions. van Dam et al. (2010) also found that when the in uence of terrain factors on ATML was taken into account, the nonlinear motion of GPS elevation time series could be better explained. Li  The increase of HYDL on land will cause the ground to sink, and the vertical position of GPS stations will move downward. The weakening of HYDL on land will cause surface to rebound, which will cause the vertical position of GPS stations to move upward.
For NTOL model, Williams and Penna (2011)  (2013b) studied the causes of nonlinear changes in 11 IGS stations in China, and believed that ELCs could effectively weaken the annual amplitude of GPS elevation time series, but the correction effect on horizontal component of GPS stations was not ideal. The displacement of the GNSS stations caused by the environmental loadings can reach several centimeters and produce a seasonal signal with an amplitude of millimeters, which not only has a signi cant impact on the estimated speed of these stations (Johnson et al., 2017), but also deviates the establishment of the earth reference frame (Freymueller, 2009;Zou et al., 2014). By correcting environmental loadings, the nonlinear deformation in GPS coordinate time series can be reduced to a certain extent. Taking into account the differences of the surface loading data sources used by different agencies, the current surface loading models can explain up to 50% of the nonlinear deformation of GNSS stations in the ideal case, and less than 20% in the horizontal direction (Yan et al., 2009;Xu et al., 2017). Therefore, the ELCs cannot completely eliminate the nonlinear deformation in GNSS coordinate time series.
At present, many research institutions around the world provide environmental loading products, including GFZ (http://rz-vm115.gfz-potsdam.de:8080/repository) (Dill and Dobslaw., 2013), EOST (http://loading.u-strasbg.fr/) in Strasbourg, France (Mémin et al., 2020) and IMLS (http://massloading.net/#Download). In addition, softwares such as QOCA (Dong et al., 2002) and OMD (Jiang et al., 2013a) can also calculate the ground displacement caused by environmental loadings. Due to different measured data (such as global atmospheric pressure, ocean bottom pressure, land water reserves), different geophysical models and methods, the environmental loading displacement at the same stations obtained by the above institutions and softwares may be quite different. namely ABOA, and compared the environmental loading products provided by GFZ, EOST and IMLS. The study found that the non-tidal oceanic loadings given by GFZ and EOST were signi cantly different.
Therefore, this study selects 21-year coordinate time series of 631 International GNSS Service (IGS) reference stations around the world, and mainly compares 6 kinds of ATML products, 7 kinds of HYDL products and 5 kinds of NTOL products provided by EOST, GFZ and IMSL. We analyze the effects of various models on GNSS coordinate time series, then select the optimal three environmental loading models, and further explore their effects on periodic signals in GNSS coordinate time series.

GNSS coordinate time series
The GNSS coordinate time series data is acquired from Scripps Orbit and Permanent Array Center (SOPAC), a data analysis center organized by IGS. This organization provides not only raw data with outliers and offsets in GNSS coordinate time series, but also clean data removing outliers. In this study, we mainly use the detrended mean GNSS time series (CleanDetrend) after removing offsets (caused by coseismic or noncoseismic reasons) and outliers, and retain the seasonal period term. The data is downloaded from ftp://garner.ucsd .edu/archive/garner/timeseries/measures/ats/.
Before data analysis, this study needs to select the acquired data from 2000 to 2021. The selection criteria are as follows. First, within recent past 21 years, the length of GNSS time series shall not be less than 3 years to ensure data integrity and the reliability of following processing results. Then, delete sites with obvious abnormal non-linear motions, including obvious post seismic deformations, local surface deformations, or abnormal motions caused by other unknown causes. Through the elimination of the above steps, 631 sites around the world are nally selected for data processing. The distribution of GNSS stations is shown in Fig. 1.

ATML model
The ATML model used in this paper includes 6 kinds of loading displacement data relative to the Earth's center of gure, among which three are provided by EOST. The rst ECMWF_IB model estimated from European Centre for Medium Range Weather Forecasts (ECMWF) operational model, assuming an inverted barometer ocean response to pressure forcing. The second is the ECMWF model, which estimated from surface pressure from ECMWF operational model, assuming a dynamic ocean response to pressure & winds from TUGO-m barotropic model (Carrère and Lyard., 2003). The third type is the era interim (ERA-interim) model. It is estimated from surface pressure from ERA interim (ECMWF Reanalysis) model, assuming an inverted barometer ocean response to pressure forcing. In addition, there is the ECMWF model provided by GFZ. The data is the elastic surface loading deformation value calculated by the method of repairing the Green function. The other two models are the Global Earth Observing System Forward Processing Instrumental Team (GEOSFPIT) model and the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA2) model provided by IMLS of the National Aeronautics and Space Administration (NASA). The spatial, time resolution and time span of the six ATML models provided by the above three agencies are shown in Table 1.  Dill, 2008). The physics and parameterization of LSDM is based on the research of Hagemann and Dümenil (1998). LSDM includes soil moisture, shallow groundwater, snow cover, and surface water stored in rivers and lakes (Hagemann and Dümenil,1998). The last two HYDL models are the Global Earth Observing System Forward Processing Instrumental Team (GEOSFPIT) and MERRA2 models provided by IMLS. The time and spatial resolution and time span of the seven HYDL models used are shown in Table   2.  Table 3. f k refers to the corresponding frequency; g i and T gi are the offset and the corresponding epoch respectively; r is the residual time series. H represents the Heaviside step function, which is expressed as follows: (2) In this article, we use correlation coe cients and RMS reduction rates to evaluate the correction effects of different environmental loading products. Before calculating these indicators, we need to deduct the linear trend term from the GNSS time series. The RMS reduction rate is expressed as follows: (3) Among them, RMS original and RMS corrected respectively represent RMS value of the GNSS coordinate time series before and after ELCs, and RMS reduction rate can be used as an index to measure the contribution of the loading models to GNSS coordinate time series. The larger of RMS reduction value, the better of the correction effect of environmental loading models. This paper also introduces amplitude and phase of annual and half-annual signals to quantitatively evaluate the contribution of environmental loading effects to nonlinear changes of GNSS coordinate time series. The difference of annual/half-annual amplitude and phase before and after ELCs is as follows. (4) Among them, A original , A corrected represents the annual amplitude or phase before and after ELCs, and the half-annual amplitude or half-annual phase. When A difference is a positive value, it indicates that the non- First, obtain the coordinate time series of 631 GNSS stations around the world from 2000 to 2021, and then use the six ATML models mentioned above, including ECMWF_IB, ECMWF, and ERA interim provided by EOST, ECMWF provided by GFZ, GEOSFPIT and MERRA2 provided by IMSL. RMS reduction rate of GNSS coordinate time series before and after the ATML correction is used to compare and analyze the correction effects of the six ATML models.
As shown in Fig. 2 to Fig. 4, this paper calculates RMS reduction rate of three-dimensional coordinate components of GNSS stations after the correction of six ATML models. Figure 2 shows the distribution of RMS reduction rates after six ATML corrections, and the percentage of GPS stations whose reduction rate is positive is marked at the upper right of each panel. It can be seen from Fig. 2 that in N direction, the model with the largest percentage of positive reduction rate (64.98%) is ECMWF_IB provided by EOST. The model with the smallest percentage of positive reduction rate (50.40%) is GEOSFPIT provided by IMSL. It can be seen from Fig. 3 that in E direction, the model with the largest percentage of positive reduction rate (66.40%) is still the ECMWF_IB model provided by EOST. The model with the smallest reduction rate is still GEOSFPIT model provided by IMSL, with a ratio of 52.14%. It can be seen from Fig. 4 that in U direction, the model with the largest percentage of positive reduction rate (80.98%) is ECMWF model provided by EOST. The reduction rate with the smallest positive value is GEOSFPIT model provided by IMSL, with a ratio of up to 60.70%.
At the same time, this article also counts the maximum reduction rate of each model in N, E and U directions, as shown in Table 4. The maximum reduction rate in each direction is the ECMWF model provided by EOST, which can reach 81.68%, 75.26% and 82.57% respectively. It indicates that the ATML model, ECMWF from EOST, has played a very active role in some GNSS stations. From Fig. 2 to Fig. 4, we can also see that even if the same ATML model provided by different agencies, the RMS reduction rate is not exactly the same. For example, the ECMWF model provided by EOST has a ratio of 60.22%, 61.97%, and 80.98% in the three directions of N, E, and U, respectively. However, the ECMWF model provided by GFZ has a positive RMS reduction rate of 51.03%, 63.07%, and 77.97% in these three directions. This may be related to the method of model establishment. The ATML model provided by EOST is estimated by the surface pressure of the ECMWF operating model, and the GFZ data is the elastic surface loading deformation value calculated by the method of repairing the Green's function. In addition, comparing Fig.   4 with Fig. 2 and Fig. 3, it can be seen that in the vertical component, the RMS reduction rate of the six ATML models is generally higher than that of the horizontal component. This shows that no matter what kind of ATML model, the in uence on the vertical direction is signi cantly greater than that on the horizontal direction.  Fig. 7 that in U direction, the model with the largest percentage of positive reduction rate (81.93%) is MERRA2 provided by EOST, and the smallest percentage of positive reduction rate (60.06) is still LSDM provided by GFZ. In general, after comparing these seven hydrological loading models, it is found that the relatively optimal performance is the MERRA2 model provided by EOST.  It can be seen from Fig. 8 to Fig. 10 that the RMS reduction rate of the NTOL is different in different dimensions. In general, the reduction rate in the vertical component is greater than the horizontal component. Moreover, the ve NTOL models corresponding to the respective components have different percentage of positive reduction rates. As shown in Fig. 8, in N direction, the model with the largest percentage of positive reduction rate (62.60%) is ECCO1 provided by EOST, and the smallest percentage of positive reduction rate (49.29%) is EMPIOM provided by GFZ. It can be seen from Fig. 9 that in E direction, the model with the largest percentage of positive reduction rate (64.18%) is still ECCO1 provided by EOST, and the smallest percentage of positive reduction rate (50.24%) is still EMPIOM model provided by GFZ. From Fig. 10, in U direction, the model with the largest percentage of positive reduction rate (80.35%) is EMPIOM provided by GFZ, and the smallest percentage of positive reduction rate (61.97%) is GLORYS2v3 provided by EOST. For some NTOL models, the impact on the horizontal component may be small, but when compared on the vertical component, the impact is more signi cant, such as EMPIOM model provided by GFZ. Therefore, when analyzing the impact of non-tidal oceanic loading on GNSS coordinate time series, the components must be discussed separately.
Finally, this article counts the maximum value of the RMS reduction rate, as shown in

The optimal combination of environmental loading models
In this paper, statistics are made on 6 atmospheric loading models, 7 hydrological loading models and 5 non-tidal oceanic loading models provided by EOST, GFZ and IMSL. The percentage of GNSS stations with positive RMS reduction rates before and after the correction of these loading models are calculated, and the results are discussed and explained above one by one. In this section, the positive RMS reduction rates of above 18 models in N, E and U directions are averaged. Then compare these averages to nd the optimal combination of environmental loading models on a global scale. The result is shown in Fig. 11.
It can be seen from Fig. 11 that among the six ATML models, the positive RMS reduction rate with the largest average percentage is ECMWF_IB provided by EOST. As indicated by the purple histogram in Fig.  11(a), the percentage can reach 69.41%. Therefore, this article selects ECMWF_IB as the optimal ATML model. At the same time, among the seven HYDL models, the MERRA2 provided by EOST has the largest average percentage of positive RMS reduction rate, which is represented by the brown histogram in Fig.  11(b), and the percentage can reach 67.14%. Therefore, the MERRA2 provided by EOST is selected as the optimal HYDL model. By analogy, the optimal model among the ve NTOL models is the ECCO1 provided by EOST, which is represented by the red histogram in Fig. 11(c). The average percentage of positive RMS reduction rate can reach 67.67%. This article next focuses on the combination of these three optimal environmental loading models as ECMWF_IB + MERRA2 + ECCO1 to analyze their impact on nonlinear changes in GNSS coordinate time series.

The in uence of the optimal ELCs on GNSS coordinate time series
The optimal atmospheric loading model ECMWF_IB, hydrological loading model MERRA2 and non-tidal oceanic loading model ECCO1 provided by EOST are used as the optimal environmental loading model combination to be discussed. In this study, the displacement sequences of these three environmental load models are cumulatively summed to further explore the nonlinear impact of the ELCs on global GNSS coordinate time series. First of all, this study analyzes and explains the magnitude of its impact. As shown in Fig. 12, in N direction, the maximum RMS value is 1.13 mm, in E direction, the maximum RMS value is 1.25 mm, and in U direction, the maximum RMS value can reach 8.42 mm. And from the RMS distribution diagrams in the three directions, it can be seen that the magnitude of the environmental loading impact in horizontal direction is still smaller than vertical direction. This shows that the environmental loading has a greater impact on GNSS coordinate time series in vertical direction, and it can be seen from Fig. 12 that the impact value is signi cantly greater in the middle and high latitude regions than in the low latitude regions.
From above analysis, it can be seen that the environmental loading has a non-negligible impact on the nonlinear deformation of GNSS coordinate time series, so this study uses Hector software to obtain the annual and half-annual amplitude and phase of GNSS coordinate time series (Bos et al., 2013(Bos et al., , 2016.
Discuss the in uence magnitude of annual and half-annual amplitude and phase before and after the ELCs. The result is shown in Fig. 13.
When only white noise + icker noise is considered, the GNSS coordinate time series have obvious annual amplitude and phase changes before and after ELCs. In E direction, the maximum annual amplitude difference can reach 1.68 mm. In 98.57% of the stations, the absolute values of the annual amplitude differences are less than 1 mm. In N direction, the maximum annual amplitude difference is 1.69 mm. In 98.09% of the stations, the absolute values of the annual amplitude differences are less than 1 mm. In U direction, the maximum annual amplitude difference can reach 8.96 mm, and 97.94% of the stations have an absolute value of the annual amplitude difference between 0 and 6 mm, and the magnitude of the impact is signi cantly greater than in the horizontal direction. At the same time, in E, N and U directions, 46.19%, 43.02% and 84.60% of the stations were corrected for environmental loading, and their non-linear amplitude was weakened. It can be further seen that non-linear in uence of environmental loading on GNSS coordinate time series is mainly re ected in the vertical direction.
In addition to annual amplitude difference, this article also statistically explains annual phase difference, as shown in Fig. 14. In E, N and U directions, the maximum absolute value of annual phase difference before and after ELCs are 178.55°, 179.15° and 174.62°, respectively. More than half of the stations have absolute annual phase differences between 0° and 60°, accounting for 80.47%, 73.01% and 72.22%, respectively. Compared with the in uence of annual amplitude, the ELCs on the annual phase is not signi cantly different in the vertical and horizontal directions, but the magnitude is more obvious in the three directions. Therefore, when considering the impact of environmental loading on annual phase, the three directions of N, E and U need to be taken into consideration.
In addition to annual amplitude and phase, the environmental loading also has a more obvious impact on half-annual amplitude and half-annual phase of GNSS coordinate time series. As shown in Fig. 15, compared with the half-annual amplitude change in horizontal direction, the environmental loading has a more obvious in uence in vertical direction, and the maximum difference can reach 1.51 mm. The maximum value of the half-annual amplitude difference in horizontal direction is less than 1 mm. Compared with the in uence of annual amplitude, the impact of environmental loading on half-annual amplitude is smaller. When considering the impact of environmental loading on half-annual amplitude, the vertical direction is the top priority. In addition to half-annual amplitude, this study also provides a statistical description of half-annual phase difference, as shown in Fig. 16. In the three directions of E, N and U, the maximum values of half-annual phase differences are 179.68°, 179.73° and 178.39°, respectively. There are 89.21%, 81.90% and 85.71% of the stations, and the absolute value of half-annual phase difference is less than 60°. Therefore, compared with annual phase, the ratio is larger, which can show that the impact on annual phase is greater than half-annual phase before and after ELCs.

Summary
This study selects 631 GNSS coordinate time series that span from 2000 to 2021 globally. A total of 18 environmental loading products, including 6 atmospheric loading, 7 hydrological loading and 5 non-tidal oceanic loading provided by EOST, GFZ and IMLS, were compared and analyzed. Among them, for the ve ATML products, in three directions of N, E and U, the maximum average percentage of positive RMS reduction is ECMWF_IB provided by EOST, and its rate can reach 69.41%. For the seven HYDL products, the maximum average percentage of positive RMS reduction is MERRA2 provided by EOST, and its ratio can reach 67.14%. For the ve NTOL products, the maximum average percentage of positive RMS reduction is ECCO1 provided by EOST, and the rate can reach 67.67%. Therefore, these three environmental loading models ECMWF_IB + MERRA2 + ECCO1 are selected as the optimal model for cumulative analysis to obtain their in uence on the nonlinear factors of GNSS coordinate time series.
Using the RMS reduction rate to classify and analyze the 18 environmental loading models, the results show that whether it is an ATML product, a HYDL product or a NTOL product, the in uence on the vertical roles. Therefore, in the later discussion and analysis, the data source and model name must be indicated.
Using the optimal environmental loading model combination obtained from the above analysis, the cumulative effect of ELCs on the periodic factors in GNSS coordinate time series is obtained. The results show that annual amplitude and half-annual amplitude mainly have a non-negligible effect on vertical direction. The annual amplitude difference and half-annual amplitude difference can reach 8.96 mm and 1.51 mm, respectively. After the ELCs, 84.60% of the non-linear deformation of the stations have been corrected. In horizontal direction, the impact is small, both less than 2mm. As for the results of annual phase and half-annual phase, it can be seen that whether it is horizontal component or vertical component, the environmental loading will have a non-negligible impact on it. Especially for annual phase, there are 19.67%, 26.98% and 27.95% of the stations in E, N and U directions, and the absolute value of phase differences are greater than 60°. Therefore, when studying the in uence of environmental loading on the phase in GNSS coordinate time series, the horizontal and vertical component need to be discussed separately.

Competing interests
The authors declare that they have no competing interests.   The distribution of the RMS reduction rate in N direction after ATML correction at different latitudes Page 21/33 Figure 3 The distribution of the RMS reduction rate in E direction after ATML correction at different latitudes  The distribution of the RMS reduction rate in U direction after ATML correction at different latitudes Page 23/33 Figure 5 The distribution of the RMS reduction rate in N direction after HYDL correction at different latitudes Figure 6 The distribution of the RMS reduction rate in E direction after HYDL correction at different latitudes Figure 7 The distribution of the RMS reduction rate in U direction after HYDL correction at different latitudes