# Comparative analysis of different atmospheric surface pressure models and their impacts on daily ITRF2014 GNSS residual time series

Zhao Li, Wu Chen, Tonie van Dam, Paul Rebischung, Zuheir Altamimi

► To cite this version:

Zhao Li, Wu Chen, Tonie van Dam, Paul Rebischung, Zuheir Altamimi. Comparative analysis of different atmospheric surface pressure models and their impacts on daily ITRF2014 GNSS residual time series. Journal of Geodesy, 2020, 94 (4), pp.42. 10.1007/s00190-020-01370-y. hal-03936628

HAL Id: hal-03936628

https://hal.science/hal-03936628v1

Submitted on 8 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comparative analysis of different atmospheric surface pressure models and their impacts on daily ITRF2014 GNSS residual time series

Zhao Li¹, Wu Chen¹, Tonie van Dam², Paul Rebischung³, Zuheir Altamimi³

1 Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University, 181 Chatham Road South, Hung Hom, Kowloon 999077, Hong Kong

2 Faculté des Sciences, de la Technologie et de la Communication, University of Luxembourg, L-4365 Esch-sur-Alzette, Luxembourg

3 IGN LAREG, Univ. Paris Diderot, Sorbonne Paris Cité, Paris, France

## Abstract

To remove atmospheric pressure loading (ATML) effect from GNSS coordinate time series, surface pressure (SP) models are required to predict the displacements. In this paper, we modeled the 3-D ATML surface displacements using the latest MERRA-2 SP grids, together with four other products (NCEP-R-1, NCEP-R-2, ERA-Interim and MERRA) for 596 globally distributed GNSS stations, and compared them with ITRF2014 residual time series. The five sets of ATML displacements are highly consistent with each other, particularly for those stations far away from coasts, of which the lowest correlations in the Up component for all the four models w.r.t. MERRA-2 become larger than 0.91.

ERA-Interim-derived ATML displacement performs best in reducing scatter of the GNSS height for 90.3 % of the stations (89.3 % for NCEP-R-1, 89.1 % for NCEP-R-2, 86.4 % for MERRA, and 85.1 % for MERRA-2). We think that this may possibly due to the 4-D variational data assimilation method applied. Considering inland stations only, more than 96 % exhibit WRMS reduction in the Up direction for all five models, with an average improvement of 3~4 % compared with the original ITRF2014 residual time series before ATML correction. Most stations (>67 %) also exhibit horizontal WRMS reductions based on the five models, but of small magnitudes, with most improvements (>76 %) less than 5 %. In particular, most stations in South America, South Africa, Oceania and the Southern Oceans show larger WRMS reductions with MERRA-2, while all other four SP datasets lead to larger WRMS reduction for the Up component than MERRA-2 in Europe.

Through comparison of the daily pressure variation from the five SP models, we conclude that the bigger model differences in the SP-induced surface displacements and their impacts on the ITRF2014 residuals for coastal/island stations are mainly due to the IB correction based on the different land-sea masks. A unique high spatial resolution land-sea mask should be applied in the future, so that model differences would come from only SP grids. Further research is also required to compare the ATML effect in ice-covered and high mountainous regions, for example, the Qinghai-Tibet Plateau in China, the Andes in South America, etc, where larger pressure differences between models tend to occur.

## Keywords

Atmospheric pressure loading (ATML), MERRA-2, NCEP Reanalysis, ERA-Interim, MERRA, GNSS position time series

## 1 Introduction

Previous research has already given evidence for correlations between atmospheric pressure loading (ATML) induced surface displacements and global GNSS height time series (van Dam et al., 1994; Brondeel and Willems, 2003; Tregoning and van Dam, 2005; van Dam et al., 2010; Dach et al., 2011; Wijaya et al., 2013). These environmentally driven surface displacements add non-linear variations to the GNSS time series being used for investigating Earth processes
associated with sea level rise, postglacial rebound, etc. To remove this environmental signal from the GNSS data, surface pressure (SP) models are used to predict the surface displacements. A commonly used dataset for modeling ATML effects include the 6-hourly, 2.5°x2.5° SP grids from the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) global reanalysis datasets (NCEP-R-1) (Kальнай et al., 1996). The 6-hourly NCEP-R-1 SP induced displacements were usually averaged into weekly ATML displacements and compared with weekly GNSS station heights (van Dam et al., 1994; Tregoning and Watson, 2009; Jiang et al., 2013). The correlations between ATML displacements and GNSS station heights at sub-weekly periods were however not studied in sufficient detail yet (Dach et al., 2011). Several new SP datasets have recently become available, that may potentially predict ATML displacements better than NCEP-R-1.

During and after the main production phase of the NCEP-R-1, several processing errors were discovered. To correct the known errors and also update the parameterizations of physical processes in NCEP-R-1, the NCEP-Department of Energy (NCEP-DOE) global reanalysis products (NCEP-R-2) was developed. Besides fixes of the human processing errors, NCEP-R-2 also incorporated upgrades to the forecast model and a diagnostic package that had been developed since the time the NCEP-R-1 system was frozen (Kanamitsu et al, 2002). Therefore, NCEP-R-2 should not be considered as a next-generation reanalysis, but should be regarded as an updated and corrected version of NCEP-R-1.

In addition to the NCEP products, there are other publicly available atmospheric reanalysis products from the European Centre for Medium-Range Weather Forecasts (ECMWF), i.e. the ERA-Interim (Berrisford et al., 2011). The model runs from 1979 to the present and is normally updated once per month (with a delay of two months) to allow for quality assurance and to correct technical problems with the production if any. The spatial resolution of ERA-Interim is 0.75°x0.75° in latitude and longitude.

In 2009, the Modern-Era Retrospective Analysis for Research and Applications (MERRA) dataset was released. The MERRA project focuses on historical climate analyses at a broad range of weather and climate time scales. It uses a state-of-the-art data assimilation system that includes a variety of modern observing systems in a climate framework to produce a temporally and spatially consistent synthesis of observations and analyses of variables not easily observed (Rienecker et al., 2011). Compared to other updated reanalysis, MERRA has achieved its goals by including improvements in precipitation and water vapor climatology. Like NCEP-R-1 NCEP-R-2 and ERA-Interim, MERRA also provides 6-hourly SP grids. However, the spatial resolution of MERRA SP grids is 2/3°x1/2° in longitude and latitude.

Finally, without a substantial investment to update assimilation routines, the MERRA data assimilation scheme was frozen in 2008 and, thus, lacked the capability to assimilate the most recent observations. As a result, a new full reanalysis was undertaken, and the second Modern-Era Retrospective analysis for Research and Applications (MERRA-2) was introduced in early 2016 to replace the original MERRA reanalysis. For this reason, The MERRA reanalysis was completed on February 29, 2016, and MERRA-2 is currently in production. The output of MERRA-2 are on a regular 0.625°x0.5° longitude by latitude grid (Bosilovich et al., 2016). MERRA-2 is the first long-term global reanalysis to assimilate space-based observations of aerosols and to represent their interactions with other physical processes in the climate system.

In this paper, we model the 3-D loading displacement induced by the latest MERRA-2 SP grids and another four products, NCEP-R-1 and NCEP-R-2, ERA-Interim and MERRA, at the same 6-hourly temporal resolution. We assess the quality of the five SP models by inter-comparing them with one another and by comparing the displacements predicted by each model with the ITRF2014 daily residual time series at 596 globally distributed GNSS stations (Altamimi et al., 2016). In particular, we determine whether SP products with a higher spatial resolution are better at reducing the scatter of daily GNSS time series in both the Up and horizontal components.

# 2 Data description

## 2.1 SP datasets for ATML modeling

We considered five SP models in this study. They all have a temporal resolution of 6-hours and
were compared over the period from 2000.01.01 until 2016.02.29, the last day of the MERRA dataset. General information about the five SP datasets is given in Table 1. The 6-hourly SP grids from NCEP-R-1 can be accessed directly via http<fnref n="1" />, with a latency of about 3-5 days. Like NCEP-R-1, the 6-hourly NCEP-R-2 SP grids can be downloaded via http<fnref n="2" />. They have the same spatial resolution as NCEP-R-1, but the latency is about 4 months due to the error fixing and parameterization update process. The ERA-Interim SP grids can be retrieved from the ECMWF server<fnref n="3" />. They have a higher spatial resolution than the NCEP-R-1 and NCEP-R-2 grids.

The MERRA SP grids are available at the Modeling and Assimilation Data and Information Services Center (MDISC), managed by the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC), as one grid file for each day<fnref n="4" />, while the MERRA-2 SP grids are also available online through GES DISC<fnref n="5" />, from January 1, 1980 to present, with a latency of about one and a half months. The spatial resolution of MERRA-2 SP grids is slightly higher as compared to MERRA.

To be consistent with our selected ITRF2014 residuals, we do not remove the diurnal (S1) and the semi-diurnal (S2) frequencies from all the five SP grids, so that the modeled ATML displacement includes the impacts of both tidal and non-tidal ATML effects (see second paragraph of Section 2.3 for details).

<table><caption>**Table 1** General information of the five SP models</caption><thead><tr><th></th><th>Number of<br>grid points in<br>longitude</th><th>Number of<br>grid points in<br>latitude</th><th>Temporal<br>resolution</th><th>Spatial<br>resolution<br>(lon x lat)</th><th>Latency</th></tr></thead><tbody><tr><td>MERRA-2</td><td>576</td><td>361</td><td>6-hourly</td><td>5/8° x 1/2°</td><td>1.5 months</td></tr><tr><td>MERRA</td><td>540</td><td>361</td><td>6-hourly</td><td>2/3°x 1/2°</td><td>Completed on<br>29 Feb., 2016</td></tr><tr><td>NCEP-R1</td><td>144</td><td>73</td><td>6-hourly</td><td>2.5° x 2.5°</td><td>3-5 days</td></tr><tr><td>NCEP-R2</td><td>144</td><td>73</td><td>6-hourly</td><td>2.5° x 2.5°</td><td>4 months</td></tr><tr><td>ERA-Interim</td><td>480</td><td>241</td><td>6-hourly</td><td>0.75° x 0.75°</td><td>3 months</td></tr></tbody></table>

## 2.2 Integer land-sea masks for ATML modeling

In addition to the global surface pressure grids, an integer land-sea mask with the same spatial resolution as the surface pressure grids is also needed to define the land-sea boundary to model the ATML effects for the land and ocean areas separately. Within the integer land-sea mask, if the grid locates in the ocean, then the grid value is assigned as 0. Otherwise, it is equal to 1. For both the NCEP-R-1 and NCEP-R-2 data, an integer land-sea mask is provided at the same spacing as the SP grids. An integer land-sea mask is also provided together with the ERA-Interim SP grids<fnref n="6" />. These three land-sea masks can be used directly during ATML modeling.

The fraction of land-covered surface in each cell of the MERRA grids can be found in the MERRA constants file collection<fnref n="7" />. Unlike NCEP and ERA-Interim, the MERRA surface parameterization includes fractional coverage of land, ocean, lake, and land ice (Lucchesi, 2012). For example, at coastlines, if the grid cell includes a mixture of land and sea coverage, thus the grid value may assign as 0.75, or 0.125, depending on the distance from the ocean. The smaller the distance, the closer the grid value to 0. Under this circumstance, if we simply use the provided file as fractional land-sea mask, then all the grid cells with non-zero value would be considered as land, which definitely results in an overestimation of the land area. Correspondingly, quite different ATML loading displacement would be obtained for stations at coastlines due to the IB correction. Therefore, in order to obtain more realistic result than using the original provided fractional land-sea mask, the MERRA constants file that describes surface coverage should be tailored to an integer land-sea mask before modeling ATML displacements. The MERRA-2

<footnote n="1">¹https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html</footnote>

<footnote n="2">²https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2 surface.html</footnote>

<footnote n="3">³http://apps.ecmwf.int/datasets/data/Interim-full-daily/levtype=sfc/</footnote>

<footnote n="4">⁴https://disc.gsfc.nasa.gov/datasets/MAI6NVANA_V5.2.0/summary?keywords=merra</footnote>

<footnote n="5">⁵https://disc.gsfc.nasa.gov/datasets/M2I6NVANA_V5.12.4/summary?keywords=MERRA-2</footnote>

<footnote n="6">⁶http://apps.ecmwf.int/datasets/data/Interim-full-invariant/</footnote>

<footnote n="7">⁷https://disc.gsfc.nasa.gov/datasets/MAC0NXASM_V5.2.0/summary?keywords=merra</footnote>
constants file<fnref n="8" /> has exactly the same data structure as the MERRA constants file, and was tailored
to a binary integer land-sea mask in the same way. For the generation of MERRA and MERRA-2
integer land-sea mask, please refer to the Appendix.

2.3 GNSS position time series

ATML driven surface displacements can be observed in GNSS vertical position time series (van Dam et al., 1994; Tregoning and van Dam, 2005; van Dam et al., 2010). Comparison between ATML and GNSS time series can be used to test the different ATML modeling methods. However, in addition to ATML signals and other known geophysical effects, e.g., continental water storage and non-tidal ocean tidal loading induced displacements (van Dam et al., 1997, 2001, 2012; Tregoning et al., 2009; Fritsche et al., 2012; Williams and Penna, 2011; Yuan et al., 2018), GNSS-related systematic errors also contribute to the observed non-linear variations in GNSS time series (Ray et al., 2008, 2013). Previous results have shown that imperfect GNSS data processing models and strategies could result in spurious non-linear displacements of GNSS stations, which mask real crustal displacements, thus affecting the comparison between GNSS and predicted loading signals (Deng et al., 2017; Li et al., 2014, 2018; King and Watson, 2010; Penna and Stewart, 2003; Penna et al., 2007; Tregoning and Watson, 2009).

To minimize the impact of the GNSS-related errors, we used here daily GNSS residual time
series that obtained in the frame of the ITRF2014 computation (Altamimi et al., 2016). The GNSS
contribution to ITRF2014 consists of daily solutions provided by the International GNSS Service
(IGS; Johnston et al., 2017) in the frame of its second reprocessing campaign (repro2; Rebischung
et al., 2016), and covering the period 1994.01.02-2015.02.14. Compared to previous ITRF
realizations, the IGS repro2 solutions have improved due to the implementation of better
background models and methodology (http://acc.igs.org/reprocess2.html). In addition, the IGS
repro2 solutions represent a sum of contributions from nine independent Analysis Centers (ACs),
and are as such potentially more precise than any one individual AC contribution. During Repro2
implementation, the IGS analysis centers (ACs) were asked not to apply model corrections for the
load displacements caused by large-scale non-tidal atmosphere, ocean and hydrological surface
fluid motions, thus the daily ITRF2014 residuals contain the non-tidal ATML signals. In addition,
only two ACs (GFZ and GRG) applied the S1/S2 atmospheric tidal corrections, while the other 7
ACs did not (Griffiths, 2019). Therefore, we finally assume that the combined ITRF2014
residuals also contain the tidal ATML effect. This is actually an inconsistency here. However,
since the impact of S1/S2 is rather small, we think that the assumption will not affect our final
results much.

In the first step of the ITRF2014 computation, the IGS repro2 daily solutions were stacked
and functional models, including piecewise linear functions, annual and semi-annual sine waves,
as well as logarithmic and/or exponential functions in case of post-seismic displacements
(Altamimi et al., 2016), were adjusted to the IGS repro2 station position time series. The GNSS
displacement time series we use here are the residuals from this stacking, to which the estimated
annual and semi-annual signals have been restored. They are free of linear trends, offsets and
post-seismic signals, and can be directly compared with ATML displacements predicted from the
different SP models. Note that our obtained ITRF2014 daily residuals (directly from IGN) contain
both position information and the formal error at each epoch, so that in the later data processing
stage we could compute the weighted root mean square (WRMS) of the GNSS coordinate time
series.

To illustrate the advantage of our used GNSS time series, Figure 1 plots the ITRF2014
residual time series together with the other three products from the Making Earth Science Data
Records for Use in Research Environment (MEaSUREs) program<fnref n="9" />, including the cleaned residual
time series from SOPAC, JPL, as well as the combined solution from SOPAC and JPL, for station
ARTU as an example. To make the figure more readable, time series of only 5 years between
2010 and 2015 are shown here. It can be clearly seen that the ITRF2014 residuals exhibit the
smallest scatter for all the three components among the four products without any sudden jumps.

<footnote n="8">⁸https://disc.gsfc.nasa.gov/datasets/M2C0NXASM_V5.12.4/summary?keywords=MERRA-2</footnote>

<footnote n="9">⁹https://cddis.nasa.gov/Data_and_Derived_Products/GNSS/SESES_time_series_products.html</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1 Cleaned residual time series for station ARTU from different institutions. The blue curve represents the product from SOPAC, red indicates the JPL solution, orange illustrates the combined solution from SOPAC and JPL, while purple demonstrates the ITRF2014 residuals.</figcaption></figure>

# 3 Data processing

## 3.1 Modeling SP induced surface displacement

We computed the 3-D surface displacements at the ITRF2014 GNSS stations induced by SP variations over a spherical grid using Farrell’s Green’s function approach (Farrell, 1972; van Dam and Wahr, 1987; van Dam et al., 2010). They can be expressed as:

$$d_n(\theta, \varphi) = \sum_{i=1}^{nlon} \sum_{j=1}^{nlat} \Delta P_{i,j} G_{i,j}^n A_{i,j}$$

$$d_e(\theta, \varphi) = \sum_{i=1}^{nlon} \sum_{j=1}^{nlat} \Delta P_{i,j} G_{i,j}^e A_{i,j}$$

$$d_u(\theta, \varphi) = \sum_{i=1}^{nlon} \sum_{j=1}^{nlat} \Delta P_{i,j} G_{i,j}^u A_{i,j}$$

where *i* and *j* denote a unique loading grid cell, *A*<sub>*i*,j*</sub> is the area of the grid cell, *nlon* and *nlat* are the numbers of grid cells in longitude and latitude respectively. *θ* and *φ* represent the colatitude and longitude of a given station where the effect of the load is being determined. *d*<sub>*n*,θ*</sub>, *d*<sub>*e*,θ*</sub>, *d*<sub>*u*,θ*</sub> are the loading displacements of the station in the North (N), East (E) and Up (U) directions. *G*<sub>*i*,j*</sub><sup>*n*</sup>, *G*<sub>*i*,j*</sub><sup>*e*</sup> and *G*<sub>*i*,j*</sub><sup>*u*</sup> denote the Green's functions for each component (Farrell, 1972), which depend on the angular distance between the loading grid cell and the station. The ITRF2014 residual displacement time series were expressed in the Center of Figure (CF) reference frame (Dong et al., 1997; Dong et al., 2003; Blewitt, 2003). Therefore, we used Green's functions in the CF frame to model ATML displacements, so as to maintain consistency between ATML and GNSS displacements. To avoid overestimating the loading displacement, we applied different data processing strategies for stations far away and close to the load grid cell. If the angular distance between the loading grid cell and the station is smaller than 5 degree, we divide the grid unit into 100 smaller units and recalculate the Green's functions for the new refined grid cells. If the angular distance between the new grid cell and the station is less than 0.5 degree, we then divide the new grid into another 100 smaller units and recalculate the Green's function for them separately. In this way, if 2.5 degree SP grid is used, for the station very close to the load grid cell, the highest spatial resolution is actually 0.025 degree (2.5/10θ/10θ). Therefore, no
matter the original load grid is sparse or dense, our code works well in calculating the modeled
ATML displacement.

Finally, $\Delta P_{i,j}$ is the pressure variation at the grid cell. Since we are only interested in time variable effects, a 10-year mean surface pressure from 2000 to 2009 (hereafter called reference pressure) is removed, so that $\Delta P_{i,j}$ corresponds to the original surface pressure from each dataset minus its corresponding reference pressure. Note that the reference pressure for each grid cell is independent of the selected continuous period, which refers to the uninterrupted SP data span from the first epoch on January 1st of the starting year until the last epoch on December 31st of the ending year. For example, if the data span is from 2000 to 2009, then the reference pressure at each grid cell equals to the summation of the gridded pressure value from the first epoch of January 1st 2000 until the last epoch of December 31st 2009 divided by the total epochs, that is, 4 epochs per day, 365 or 366 days per year, depending on the leap year or not.

During the calculation, an integer land-sea mask with the same spatial resolution as the surface pressure grids was used to define the land-sea boundary, and the pressure loading over the oceans was determined using a modified inverted barometer (IB). For a detailed discussion of the modified IB, readers are referred to van Dam and Wahr (1987). Concerning this, we should keep in mind that any observed differences in the displacement for coastal sites would depend not only on the pressure grid itself, but also on the IB correction based on different land-sea mask that used for each SP data-set.

## 3.2 Comparison of ATML displacements and ITRF2014 residuals

After having obtained predicted ATML displacement time series from the above five SP datasets, we compared them to the ITRF2014 residual time series of a selected set of GNSS stations. Due to the relatively lower quality of IGS products before 2000 (Griffiths, 2019), we restricted this comparison to the period after 2000, and first selected the ITRF2014 GNSS stations with at least 5.5 years of observations after 2000, so as to obtain a more reliable global statistical result, and also to ensure that there are as many stations as possible that evenly distributed around the globe. We additionally removed stations with abnormally scattered residual time series, i.e. with standard deviations larger than 10 mm and 5 mm for the Up and horizontal components, respectively. This selection procedure left us with 596 globally distributed GNSS stations shown in Figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2 Spatial distribution of the 596 selected GNSS stations</figcaption></figure>

Our ATML displacement time series all have four epochs per day at 00:00, 06:00, 12:00 and 18:00. In order to make them comparable with our daily GNSS time series, we first averaged them into daily estimates (no weighting has been applied here for each epoch), and then applied them as a posteriori corrections to remove the ATML effect from the daily ITRF2014 residuals. This procedure is easy to implement and also reasonable, as Tregoning and van Dam (2005) have already proved that the majority of the ATML effect could be modeled by applying a daily-averaged correction. Moreover, it has advantages for the rigid processing scheme applied within the IGS, since in this way the ATML effect can be corrected after data analysis and before combination in a consistent way using identical values for solutions based on the same or different space-geodetic techniques (Dach et al., 2011). A greater improvement would be expected for the
performance of each SP model in correcting the GNSS scatter if the ATML correction is applied at the observation level (Tregoning and van Dam, 2005; Dach et al., 2011). Nevertheless, no matter whether the ATML correction are applied on the observation level or as a posteriori daily averages, we think that the SP model difference would not be affected much, since we use the same GNSS data during comparison.

Finally, according to our previous research on the continental water storage (CWS) models, we found that the MERRA-Land product with the highest spatial resolution at the time performed the best in reducing the weekly scatter of the GPS height (Li et al., 2015). Therefore, we speculate that the ATML model would exhibit the same characteristics. Since MERRA-2 grid is currently the most recent product with the highest spatial resolution among the five SP datasets, we first described the ATML displacements derived from MERRA-2 and compared them with the selected ITRF2014 GNSS residual time series. We then used the MERRA-2 results as the reference to compare the differences between ATML models.

# 4 Comparison of MERRA-2 derived ATML displacements with ITRF2014 GNSS residual time series

## 4.1 Global variability of MERRA-2-derived ATML displacements

To have a global view of the MERRA-2-derived ATML displacement, Figure 3 depicts the time series of some example stations with maximum ATML displacement for the Up, North and East component, while Figure 4 depicts the scatter (top panels) and the absolute maximum displacement (bottom panels) of the 6-hourly ATML time series for the Up component of the selected 596 stations during the period 2000.01.01-2016.02.29. For clear comparison, the ITRF2014 residuals time series for the example stations are also presented in Figure 3. Note that in the middle panel of Figure 3, there's a sudden large deviation in the North component of CAS1 on 2008-12-17. Then from 2008-12-18 until 2008-12-31, there is no data in the ITRF2014 residuals. We checked the log file for station CAS1 from IGS website<fnref n="10" /> and only found that the receiver type had been changed on 2008-12-3 from ASHTECH to LEICA. Since we could not find the exact reason that caused the sudden deviation, this epoch was excluded in the later comparison. Also, because coastal regions cannot be represented accurately enough by the relatively sparse spatial resolution of our selected SP grids, here we defined a special group of inland stations, which was model-dependent based on the corresponding land-sea mask, so as to separate the impact of the IB correction. Specifically, inland stations indicated all stations 50 km, 54 km, 60 km, and 200 km away from nearest coast for MERRA-2, MERRA, ERA-Interim and the NCEP grids, while the number of the inland stations was 174, 168, 166 and 113, respectively.

We notice that the spatial patterns in Figure 4 are globally consistent with previous research focused on ATML displacements (van Dam et al., 1994; Dach et al., 2011). Due to the small pressure variation and the assumed IB response, the smallest ATML displacements in the Up component predicted from MERRA-2 SP mostly locate in the low and mid-latitude oceanic and coastal regions, with scatters of less than 2 mm and maximum displacements of 5 mm (Figure 4). Most continental stations in the Northern Hemisphere, especially the inland stations, have scatters of greater than 3 mm. In particular, the absolute maximum ATML displacements of some stations located at high latitudes reach more than 18 mm due to the larger storms found there (bottom panel of Figure 4), for example, station ARTU (Russia, 24.2 mm, top panel of Figure 3), etc. Because of the irregular distribution of land and ocean on the Earth's surface, there are significantly more ocean and much less land in the Southern Hemisphere, so that the Southern Hemisphere climates tend to be slightly milder than those at similar latitudes in the Northern Hemisphere. Correspondingly, the pressure variation and its induced ATML displacement for the Up component of continental stations in the Southern Hemisphere are also about 1/3~1/2 of the magnitude for those in the Northern Hemisphere (Figure 4).

The MERRA-2 SP variations could only induce very small displacements in the East and North components (horizontal hereafter), with scatters of less than 1 mm and maximum displacements of 4 mm, respectively. To save space, in the following analyses, most of the figures

<footnote n="10">¹⁰ ftp://ftp.igs.org/pub/station/log/cas1_20190412.log</footnote>
that illustrate the characteristics of the horizontal components are presented in the supplementary material. The same as the Up component, stations at low and mid-latitudes exhibit the smallest horizontal displacements due to the small pressure variation, particularly for those located in the oceans and along the coasts (Figure S1, S2), where the IB response exhibits the biggest impact. The largest horizontal ATML displacements are also found at high latitudes where big storms occur, especially in Antarctica, Northern Europe and Russia. For example, station CAS1 (Antarctica, middle panel of Figure 3) and SCOR (Greenland, bottom panel of Figure 3) exhibit the biggest horizontal scatter of 0.99 mm and 0.94 mm for the North and East component, respectively.

<figure><img src="image_4.png" /><figcaption>Figure 3 Time series of stations with maximum MERRA-2-derived ATML displacement (blue curve) together with their corresponding ITRF2014 residuals (red curve)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4 RMS (top) and maximum amplitudes (bottom) of the 6-hourly MERRA-2-derived ATML displacement time series for the Up component. Units are mm.</figcaption></figure>

## 4.2 Correlation between MERRA-2-derived ATML displacements and ITRF2014 residuals

To compare MERRA-2-derived ATML displacements with the selected ITRF2014 GNSS residual time series, we first calculated the correlation coefficients between both data sets (Figure 5; Figure 6). The correlations are generally low in the horizontal components, with only 9.6 % and 5.7 % of all stations showing correlations of higher than 0.4, while 11.9 % and 7.4 % even exhibiting negative correlation in the East and North components, respectively, among which most belong to stations near the coasts.

The correlations for the Up component are larger, with a mean value of 0.33, and 34.4 % of all stations exhibit correlations of higher than 0.4, mostly in the middle- and high-latitude continental regions (top panel, Figure 5). For some Eurasian continental stations, the correlations even reach more than 0.7, for example at station IKRM (Russia, 0.78) and BJFS (China, 0.74). Stations with poor correlations are those located in the low and mid-latitude oceanic and coastal regions, some of which also having negative correlations, e.g. station XMIS (Christmas Island, Australia) with the largest negative correlation of -0.22. Among the 174 inland stations (Figure 6), there exists no negative correlations in the Up direction, and the mean correlation is 0.43, which is much larger than that for all stations (0.33).
<figure><img src="image_6.png" /></figure>

Figure 5 Correlation coefficients between MERRA-2-derived ATML displacements and the ITRF2014 residual time series for all 596 stations.

<figure><img src="image_7.png" /></figure>

Figure 6 Histogram of correlation coefficients between MERRA-2-derived ATML displacements and the ITRF2014 residual time series. The X label describes the correlation coefficient. The Y label illustrates the number of stations with correlation inside different ranges. Blue color represents the results for all 596 stations, while
## 4.3 WRMS reduction of ITRF2014 residuals

We calculated the relative WRMS reductions (see detailed definition in Jiang et al., 2013) obtained when removing the MERRA-2-derived ATML displacements from the ITRF2014 residual time series. The relative WRMS reductions of individual stations are shown in Figure 7, while Figure 8 depicts the number of stations with relative WRMS reductions inside different ranges. We find that 85.1 %, 73.3 % and 71.8 % of all stations have their WRMS reduced by ATML corrections in the Up, East and North components respectively. In the Up component, this number is larger than those obtained in previous studies (e.g., van Dam et al., 1994, 2010; Brondeel and Willems, 2003; Wijaya et al., 2013). However, here we cannot yet determine whether this improvement stems from using improved GNSS coordinate time series or the MERRA-2 SP grids. We will explore this question further in the next section.

<figure><img src="image_8.png" /><figcaption>Figure 7 Relative WRMS reductions obtained when correcting ITRF2014 residual time series from MERRA-2-derived ATML displacements for all 596 stations. Unit is %</figcaption></figure>

The relative horizontal WRMS reductions are very small, usually less than 5 %, with mean reduction of 3.1 % for the North component and 3.5 % for the East component (Figure 8; Table 2, Section 5.2). For inland stations, the mean horizontal WRMS reduction is only slightly larger, at around 3.6 %. The largest improvements are obtained for station POL2 in Kyrgyzstan and station SCOR in Greenland, with relative WRMS reductions of 15.0 % in the North and 22.2 % in the East.
The relative WRMS reductions for the Up component are globally larger, with a mean reduction of 8.3 %, and their spatial distribution is similar to that of the correlation coefficients in the Up direction (Figure 5). Due to the small pressure variation and the assumed IB response, stations in the low and mid-latitude oceanic and coastal regions mainly show small relative WRMS reductions (< 5 %), while those stations in the middle- and high-latitude continental regions, where larger storms occur, mostly exhibit large relative WRMS reductions (> 10 %), particularly for the Northern Hemisphere because of the irregular distribution of the land and ocean (see Section 4.1). The largest improvement reaches more than 30 % at some stations, for example, station ULAB in Mongolia (34.3 %), stations BJFS and CHAN in China (33.9 %). If only considering inland stations, 97.1 % (174 in total) show WRMS reduction in the Up component, with mean reduction of 11.4 %, and the best improvement concentrates mainly in the Eurasia continent (Figure 7, 8).

<figure><img src="image_9.png" /><figcaption>Figure 8 Histogram of relative WRMS reduction. The X label describes the WRMS reduction (Unit: %). The Y label illustrates the number of stations with WRMS reduction inside different ranges. Blue color represents the results for all 596 stations, while orange denotes the 174 inland stations</figcaption></figure>

# 5 Discussion

From the results presented in Section 4, we conclude that the scatter of the daily ITRF2014 GNSS residual time series can be globally reduced using the latest, high-resolution MERRA-2 SP grids, particularly in the Up component. In this section, we compared the ATML displacements derived from different SP datasets (Section 5.1) and their abilities to reduce the scatter of GNSS time series (Section 5.2 and Section 5.3). Then we compared the daily pressure variation from the five SP dataset to further interpret the ATML displacement differences (Section 5.4).

## 5.1 Correlation between ATML displacements derived from MERRA-2 and other SP datasets

To investigate the differences of the ATML displacements derived from different SP datasets, we calculated the correlation coefficients between MERRA-2-derived displacement time series and time series derived from the other SP dataset at the 596 selected GNSS stations. The correlation results for the Up component are shown in Figure 9, while those for the horizontal components are plotted in Figure S3a and Figure S4a. Figure 10 depicts the histogram of the correlation coefficients for both all stations (blue color) and inland stations (orange color). We find that in general all the five sets of ATML displacement time series are very consistent with each other.

Among the four other datasets, MERRA exhibits the highest correlations with MERRA-2 (bottom panel of Figure 9, S3a, S4a; the last three panels of Figure 10), with mean correlations of more than 0.995 in all three components. The correlations between NCEP-R-2 and MERRA-2 are slightly lower (top panel of Figure 9, S3a, S4a; the first three panels of Figure 10), with mean
correlations of 0.986, 0.985 and 0.989 in the North, East, and Up components respectively. Further, only a few island and coastal stations (less than 3 % of all stations) have correlations less than 0.90. NCEP-R-1 correlates closely with MERRA-2 as NCEP-R-2 (the second panel of Figure 9, S3a, S4a; the three panels in the second row of Figure 10), while ERA-Interim shows the lowest correlations with MERRA-2 (the third panel of Figure 9, S3a, S4a; the three panels in the third row of Figure 10), but with mean correlations still larger than 0.96 in all three components. Again, the lowest correlations ($< 0.8$) are found in the low and mid-latitude oceanic and coastal regions, particularly for the island stations. There are two possible reasons that may cause this low correlation: the land-sea masks used for MERRA-2 and ERA-Interim may be quite different, or there may be quite different SP distribution over the oceans for these two products.

<figure><img src="image_10.png" /><figcaption>Figure 9 Correlation coefficients for the Up component between ATML displacements derived from MERRA-2 and the other 4 SP datasets for all 596 stations. Black dots indicate the correlation coefficients are smaller than the minimum value of the scale</figcaption></figure>

With respect to inland stations, we observe that the correlation between ATML displacement derived from MERRA-2 and the other four SP datasets are larger than those near the coasts. In particular, the lowest correlations in the Up direction are larger than 0.91 for all the four dataset (Figure 10), while the horizontal correlation is slightly smaller, with the lowest correlation locates at Station HARB (0.61, South Africa, 480 km from coast) and RIOP (0.83, South America, 146
km from coast) in the North and East component for the ERA-Interim product. The difference between all stations and inland stations further confirms our deficiency during ATML modeling for coastal regions. However, at this stage, we cannot decide whether the SP grid itself, or the land-sea mask is the main cause for the bigger differences near the coast. We will answer this question later.

<figure><img src="image_11.png" /><figcaption>Figure 10 Histogram of correlation coefficient between MERRA-2 and other model derived ATML displacement for the Up, East and North components. The X label describes the correlation coefficient. The Y label illustrates the number of stations with correlation coefficient inside different ranges. Blue color represents the results for all 596 stations, while orange denotes the results for inland stations</figcaption></figure>

## 5.2 Relative WRMS reductions of ITRF2014 residuals based on different SP datasets

In this section, we removed the predicted ATML displacement based on various SP models from the GNSS coordinate time series. We first calculated the relative WRMS reductions of the ITRF2014 residual time series for each SP dataset, as we have done in Section 4.3. Table 2 gives the percentages of stations with WRMS reduced after correcting the ATML effects using different models together with their average and maximum WRMS reduction for both all stations and inland stations. In general, all five SP datasets are able to reduce the non-linear variations in the Up component of GNSS time series. ERA-Interim performs the best, NCEP-R-1 and NCEP-R-2 perform slightly worse, MERRA ranks fourth, and MERRA-2 shows the worst performance, although MERRA-2 is currently the most advanced SP dataset with the highest spatial resolution. This is a quite different situation as with continental water storage, for which the higher spatial resolution MERRA dataset was found to perform best (Li et al., 2015). We think that the differences mainly come from the different assimilation methods that used for each model, since similar raw observations have been assimilated for all the SP products. The advantage of ERA-Interim over the other datasets may be due to the four-dimensional variational data assimilation method (4D-Var) applied. In addition, its data assimilation benefits from more extensive use of radiances with an improved fast radiative transfer model. Both NCEP-R-1 and NCEP-R-2 only adopt a three-dimensional variational data assimilation system (3D-Var), thus their performances are slightly worse. The MERRA and MERRA-2 also use 3D-Var as the assimilation framework, and the incremental analysis updates (IAU) procedure is implemented to slowly adjust the model state toward the observed state. However, the focus of the assimilation process is to simulate the hydrological cycle correctly (Lindsay et al., 2014). The water cycle has been improved (Li et al., 2015), but the quality of the surface pressure seems to be slightly degraded.

Depending on the SP dataset used, between 85.1 % and 90.3 % of the selected GNSS stations have their WRMS reduced by ATML corrections in the Up component for all 596 stations, with average WRMS reduction of around 8.3 % (Table 2). These values are even 11 %~16 % larger than the weekly comparison result reported after considering the topographic effects (van Dam et
al., 2010). Since our numbers are similar among the different tested SP datasets, we think that this improvement over previous studies is mainly due to our use of latest-generation, reprocessed GNSS time series. In addition, most of the stations (> 67 %) also have their horizontal WRMS reduced by ATML corrections using any of the five SP datasets, although the magnitude of the horizontal WRMS reductions is quite small, with average WRMS reduction at around 3.1 % (Table 2), and most improvements (>76 %) are less than 5 %. An important thing that should be noticed is that although ERA-Interim does the best performance in reducing the WRMS for the Up and East components among the five SP models, it performs the worst in reducing the WRMS for the North component (Table 2). One possible reason may be the interference between ERA-Interim-derived ATML displacement and the impact of the temperature gradient in the North-South direction. Further research is still required to investigate the exact reasons.

With respect to inland stations only, we observe that more than 96 % of the remaining stations exhibit WRMS reduction in the Up component for all the five SP datasets, and the average WRMS reduction is also 3~4 % larger compared to the result of all 596 stations (Table 2). Nevertheless, with respect to the horizontal component, there is only a slight improvement on the average WRMS reduction (~3.5 %) and the percentage of stations with WRMS reduced (> 67 %). In particular, the maximum horizontal WRMS reduction exhibits mostly near the coasts using any of the products, especially for the E component, which is quite different from the Up component.

<table><caption>Table 2 Statistics of the relative WRMS reductions obtained with different SP datasets for all 596 stations and inland stations</caption><thead><tr><th rowspan="2">Model name</th><th rowspan="2">Number of stations</th><th colspan="2">U component</th><th colspan="2">E component</th><th colspan="2">N component</th></tr><tr><th>Percentage of stations with WRMS reduced (%)</th><th>Average/Maximum WRMS reduction (%)</th><th>Percentage of stations with WRMS reduced (%)</th><th>Average/Maximum WRMS reduction (%)</th><th>Percentage of stations with WRMS reduced (%)</th><th>Average/Maximum WRMS reduction (%)</th></tr></thead><tbody><tr><td rowspan="2">MERRA-2</td><td>596</td><td>85.1</td><td>8.3/34.3</td><td>73.3</td><td>3.5/22.2</td><td>71.8</td><td>3.1/15.0</td></tr><tr><td>174</td><td>97.1</td><td>11.4/34.3</td><td>75.3</td><td>3.6/21.7</td><td>69.0</td><td>3.6/15.0</td></tr><tr><td rowspan="2">MERRA</td><td>596</td><td>86.4</td><td>8.3/33.8</td><td>74.0</td><td>3.5/22.0</td><td>72.5</td><td>3.1/15.1</td></tr><tr><td>168</td><td>96.4</td><td>11.5/33.8</td><td>76.2</td><td>3.5/21.6</td><td>67.9</td><td>3.6/15.1</td></tr><tr><td rowspan="2">ERA-Interim</td><td>596</td><td>90.3</td><td>8.5/35.2</td><td>75.0</td><td>3.3/22.0</td><td>67.8</td><td>2.7/13.8</td></tr><tr><td>166</td><td>98.8</td><td>11.9/35.2</td><td>77.1</td><td>3.4/20.2</td><td>69.9</td><td>3.1/13.8</td></tr><tr><td rowspan="2">NCEP-R-1</td><td>596</td><td>89.3</td><td>8.0/33.9</td><td>74.0</td><td>3.4/22.5</td><td>72.5</td><td>2.9/13.5</td></tr><tr><td>113</td><td>96.5</td><td>12.3/33.9</td><td>75.2</td><td>3.5/15.0</td><td>73.4</td><td>3.5/12.8</td></tr><tr><td rowspan="2">NCEP-R-2</td><td>596</td><td>89.1</td><td>8.1/34.5</td><td>74.7</td><td>3.4/22.1</td><td>72.0</td><td>2.9/13.8</td></tr><tr><td>113</td><td>96.5</td><td>12.4/34.5</td><td>77.0</td><td>3.5/15.0</td><td>72.6</td><td>3.5/12.3</td></tr></tbody></table>

## 5.3 Relative WRMS reduction difference of the SP models w.r.t MERRA-2

We then compared the relative WRMS reductions obtained with the four other SP datasets to the relative WRMS reductions obtained with MERRA-2. The differences in relative WRMS reductions for the Up component are shown in Figure 11a, while the results for the E and N components are plotted in Figure S5 and S6, respectively. Figure 11b illustrates the histogram of relative WRMS reduction differences for the U, E and N components. Negative values mean that MERRA-2 reduces the scatter of GNSS time series more than the other SP datasets. We observe that there is very little difference between MERRA and MERRA-2 in terms of relative WRMS reductions (bottom panel of Figure 11a; S5; S6). This was expected from the very high
correlations between MERRA-derived and MERRA-2-derived ATML displacement time series
(see Section 5.1).

<figure><img src="image_12.png" /><figcaption>Figure 11a Relative WRMS reduction differences (%) of different models w.r.t MERRA-2 for the Up component of all 596 stations.</figcaption></figure>

The relative WRMS reduction differences between ERA-Interim and MERRA-2 are larger (the third panel of Figure 11a; S5; S6). This is also expected given the lower correlations between ERA-Interim-derived and MERRA-2-derived ATML displacement time series (see Section 5.1). ERA-Interim is better at reducing the scatter of GNSS height time series than MERRA-2 for most stations (65.1 %), especially in North America, Eurasia and Antarctica. However, the majority of stations in South America, South Africa, Australia, together with some island stations in the Indian Ocean and the Pacific Ocean have their WRMS in the Up component more reduced with MERRA-2. On the other hand, MERRA-2 performs better at reducing the scatter of horizontal GNSS time series than ERA-Interim for more than half stations, especially in coastal regions. Nevertheless, most stations in Europe have their horizontal WRMS more reduced with ERA-Interim.
Using NCEP-R-1, 41.9 %, 49.0 % and 47.1 % of the stations show improved WRMS reductions for the North, East, and Up components, respectively, as compared to using MERRA-2 (see the second panel of Figure 11a). The biggest improvements are located mainly in Europe, together with several coastal stations in North America, Antarctica, and New Zealand. The spatial distribution of the relative WRMS reduction differences between NCEP-R-2 and MERRA-2 is quite similar as between NCEP-R-1 and MERRA-2 (top panel of Figure 11a). NCEP-R-2 is slightly better at reducing the scatter of GNSS time series for most stations in Europe than NCEP-R-1 and MERRA-2. However, for most island and coastal stations, MERRA-2 is superior to NCEP-R-2.

<figure><img src="image_13.png" /><figcaption>Figure 11b Histogram of relative WRMS reduction differences (%) of different models w.r.t MERRA-2 for the Up (U), East (E) and North (N) components. The X label describes the WRMS reduction difference. The Y label illustrates the number of stations with WRMS reduction difference inside different ranges. Blue color represents the results for all 596 stations, while orange denotes the results for inland stations</figcaption></figure>

The WRMS reduction difference between model and MERRA-2 become remarkably smaller for the inland stations, particularly for MERRA, of which 98.8 % of the stations show differences of less than 1 % for the Up component, while all the horizontal differences are within 1 % (bottom panel of Figure 11b). NCEP-R-1 and NCEP-R-2 show slightly bigger difference (top and the second panel of Figure 11b). Stations with WRMS reduction differences of smaller than 1 % account for more than 82 %. ERA-Interim exhibits the biggest difference with respect to MERRA-2, but still more than 93 % of the stations showing WRMS reduction difference of smaller than 2 % for all the three component (the third panel of Figure 11b). This comparison result further confirms the finding from Section 5.1 that stations near the coasts exhibit bigger difference than those far away from coasts in the SP-derived ATML displacement.

## 5.4 Daily surface pressure variation derived from different SP datasets

To investigate the differences between modeled ATML displacements from the SP grid itself, Figure 12 illustrates the daily pressure variation difference of the four SP products w.r.t MERRA-2 on 2014.07.15 as an example. Positive (reddish color) means the MERRA-2 pressure variation is smaller than the selected model, while negative (bluish) means MERRA-2 is larger. The Generic Mapping Tools (GMT) software was used here to interpolate all the grids into the same spatial resolution as MERRA-2 for differencing (Wessel and Smith, 2018). Note that here the date 2014.07.15 was arbitrarily selected, as the daily pressure variation differences were almost the same for different days in different seasons.

We observe that in general MERRA shows the smallest pressure variation difference w.r.t MERRA-2 (bottom panel of Figure 12). Grid points with difference of smaller than 20 millibar (mbar) account for 93.8 %, including the entire ocean. ERA-Interim exhibits slightly bigger difference in the mountain regions compared with MERRA (the third panel of Figure 12), while NCEP-R-2 and NCEP-R-1 (top and the second panel of Figure 12) present much bigger pressure variation difference w.r.t MERRA-2 in ice-covered areas, like Greenland and Antarctica, as well
as all the other mountain regions, especially for the NCEP-R-2 product. In particular, MERRA-2
exhibits much bigger pressure variation (>360mbar) in Mount Roraima (South America) than
MERRA and ERA-Interim, while both the NCEP products show the maximum positive difference
of more than 700mbar w.r.t MERRA-2 in the New Guinea Highlands (Oceania). Our observed
pressure difference between NCEP and ERA-Interim is more than 10-times bigger than that from
Salstein et al. (2008). The big discrepancy should be caused by the different surface topographic
height assumed for the two models. If referenced all the five SP models into the same surface
topography, the pressure differences should be at least 10-times smaller. However, the observed
global trend of the pressure difference between models here is correct. Due to the less
well-constrained analyses, larger pressure differences between models tend to occur in the
ice-covered area and the high mountainous regions (Salstein et al., 2008).

<figure><img src="image_14.png" /><figcaption>Figure 12 Daily pressure variation differences of the four SP models w.r.t MERRA-2 on 2014.07.15. Unit of the pressure variation difference is mbar. From top to bottom panels represent the NCEP-R-2, NCEP-R-1, ERA-Interim and MERRA product.</figcaption></figure>

Compared with the spatial pattern of the model differences described in the previous sections,
we can therefore answer the question in Section 5.1 that the bigger model differences for
coastal/island stations are mainly due to IB correction based on the different sparse land-sea
masks. This is one of the deficiencies in this research. Currently, the Special Bureau for the
Atmosphere (SBA) under the Global Geophysical Fluid Center (GGFC) within the International
Earth Rotation and Reference Systems Service (IERS) provides operational ATML displacement
for GNSS sites worldwide from NCEP and ECMWF<fnref n="11" />. They use a much more sophisticated

<footnote n="11">11 https://www.aer.com/science-research/earth/earth-mass-and-rotation/special-bureau-atmosphere/</footnote>
treatment of the coastal zone based on the IB assumption, e.g. a land-sea mask with 0.25-degree for NCEP<fnref n="12" /> or 0.1-degree for ECMWF<fnref n="13" />, thus could represent more realistic displacement for coastal/island sites. To avoid the shortcoming of our ATML calculation, readers could download the SBA provided ATML displacement for stations near the coasts from different SP models and compare their differences instead. Moreover, with respect to those stations far away from coasts, Figure 12 confirms that slightly lower correlations exhibit at higher altitude with bigger pressure variations. Since most of our selected GNSS sites are located at very low altitude, further research is still required to focus on the ATML effects in the high mountainous regions, e.g. the Qinghai-Tibet Plateau in China, the Andes in South America, etc.

# 6 Summary

ATML displacements contribute to the non-linear variations in the GNSS coordinate time series. To remove this environmental signal from the GNSS data, SP models are required to predict ATML displacements. Based on Farrell's Green's Function approach, we modeled the 3-D surface displacements induced by the latest MERRA-2 SP grids, together with the other four products (NCEP-R-1, NCEP-R-2, ERA-Interim and MERRA) over the period 2000.01.01-2016.02.29 at 596 GNSS stations from the ITRF2014 network. Our results confirm that SP variations induce only small displacements in the low- and mid-latitude oceanic and coastal regions, while most continental stations far away from coasts in the Northern Hemisphere exhibit ATML displacement scatters of larger than 3 mm for the Up component. The five sets of ATML displacement are highly consistent with each other, among which MERRA exhibits the highest correlation with MERRA-2, NCEP-R-1 and NCEP-R-2 rank the second and the third, while ERA-Interim shows the lowest correlations, with only several coastal/island stations having correlation of lower than 0.9.

Comparing the modeled ATML displacements with daily ITRF2014 residual time series, we find that all five SP products perform well in correcting the scatter of GNSS height, and the percentage of stations with WRMS reduced improve by more than 10 % compared with previous results, due to our used latest-generation, reprocessed GNSS time series. Depending on the different assimilation methods, ERA-Interim performs best, reducing the scatter for 90.3 % of the GNSS height, with maximum improvement of 35.2 %. NCEP-R-1 and NCEP-R-2 do slightly worse (89.3 %, 89.1 %), MERRA ranks the fourth (86.4 %), while MERRA-2 exhibits the worst performance (85.1 %), although it is the most recent advanced product. Most stations (> 67 %) also exhibit horizontal WRMS reductions, but the improvement is usually less than 5 %. Further improvement would be expected if correcting the ATML effect at the observation level rather than as daily averaged a posteriori corrections here.

Compared with MERRA-2, more than half of the stations show larger WRMS reductions in the Up direction using MERRA and ERA-Interim. However, MERRA-2 performs better in reducing the horizontal scatter than ERA-Interim, as well as 3-D scatter than both NCEP products for more than half stations. In particular, the majority of stations in South America, South Africa, Oceania and the southern oceans show larger WRMS reductions with MERRA-2 than with any other SP datasets, while all other four SP datasets enable larger WRMS reduction in the Up component than MERRA-2 in Europe.

The model difference w.r.t MERRA-2 become remarkably smaller for inland stations only, with the lowest correlations in the Up component with MERRA-2 of larger than 0.91 for all the four datasets, while more than 96 % of the inland stations exhibit WRMS reduction for all the five SP datasets in the Up direction, and the average WRMS reduction is also improved by 3~4 %. Through comparison of the daily pressure variation between models, we conclude that the bigger model differences for coastal/island stations mainly due to the different sparse land-sea masks used during IB correction. In future work, a unique high spatial resolution land-sea mask should be applied during comparison, so that the model differences would come from only SP grids. Further investigation is also required to evaluate the ATML effect in Greenland, Antarctica and other mountain areas where larger pressure differences occur among different SP models, e.g. the Qinghai-Tibet Plateau in China.

<footnote n="12">12 http://geophy.uni.lu/ggfc_atmosphere/NCEP-loading.html</footnote>

<footnote n="13">13 http://loading.u-strasbg.fr/displ_all.php</footnote>
Finally, only ATML effect is considered here. Taking into account the hydrology loading and the non-tidal ocean loading would help reduce the GNSS scatter more, but not the full remaining part. In-depth research is still expected to investigate the possible reasons that may cause the inconsistency between loading and the GNSS displacement, e.g. the temperature gradient, etc.

## Author Contribution Statement

Z.L. designed and performed the research; P.R. and Z.A. provided the GPS data; T.V.D provided the software;

Z.L. analyzed the data and wrote the manuscript; W.C, T.V.D and P.R. revised the manuscript.

## Data Availability Statement

The ATML time series involved in the research and the scripts to generate the integer land-sea masks for MERRA and MERRA-2 are available upon request.

## Acknowledgement

We thank Mike Bosilovich for helping us generate the MERRA integer land-sea mask. We thank the IGN for providing us the latest ITRF2014 residuals. We also thank Dr. Xavier Collilieux for helping us improve our paper. Figures in this paper are plotted with the GMT and MATLAB software. This research is supported by the National Key Research and Development Program of China (Project 2016YFB0502101), and the European Commission / Research Grants Council (RGC) Collaboration Scheme sponsored by the Research Grants Council of Hong Kong Special Administrative Region, China (Project No. E-PolyU 501/16), the National Science Foundation for Distinguished Young Scholars of China (Grant No. 41525014).

## References

Altamimi Z, Rebischung P, Métivier L., Collilieux X (2016) ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions. J Geophys Res 121(8):6109-6131.
doi:10.1002/2016JB013098

Berrisford P, Dee D, Poli P, Brugge R, Fielding K, Fuentes M, Kačlberg P, Kobayashi S, Uppala S, Simmons A (2011) ERA report series: The ERA-Interim archive Version 2.0, pp1-27. European Centre for Medium Range Weather Forecasts, Shinfield Park, Reading, Berkshire RG2 9AX, United Kingdom

Blewitt G (2003) Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth.
J Geophys Res 108(B22103). doi:10.1029/2002JB002082

Bosilovich MG, Lucchesi R, Suarez M (2016) MERRA-2: File Specification. GMAO Office Note No. 9 (Version 1.1), 73 pp, available from http://gmao.gsfc.nasa.gov/pubs/office_notes

Brondeel M, Willems T (2003) Atmospheric pressure loading in GPS height estimates. Adv Space Res 31(8): 1959-1964. doi: 10.1016/S0273-1177(03)00157-1

Dach R, Böhm J, Lutz S, Steigenberger P, Beutler G (2011) Evaluation of the impact of atmospheric pressure loading modeling on GNSS data analysis J Geod 85: 75–91. doi:10.1007/s00190-010-0417-z.

Dee DP, Uppala SM, Simmons AJ, Berrisford P, Poli P, Kobayashi S, Andrae U, Balmaseda MA, Balsamo G, Bauer P, Bechtold P, Beljaars ACM, van de Berg L, Bidlot J, Bormann N, Delsol C, Dragani R, Fuentes M, Geer AJ, Haimberger L, Healy SB, Hersbach H, Hólm EV, Isaksen L, Kållberg P, Köhler M, Matricardi M, McNally AP, Monge-Sanz BM, Morcrette JJ, Park BK, Peubey C, de Rosnay P, Tavolato C, Thépaut JN, Vitart F (2011) The ERA-Interim reanalysis: configuration and performance of the data assimilation system. Q J R Meteorol Soc 137: 553–597. doi:10.1002/qj.828

Deng LS, Jiang WP, Li Z, Chen H, Wang KH, Ma YF (2017) Assessment of second- and third-order ionospheric effects on regional networks: case study in China with longer CMONOC GPS coordinate time series. J Geod 91:207–227. doi: 10.1007/s00190-016-0957-y

Dong D, Dickey JO, Chao Y, Cheng MK (1997) Geocenter variations caused by atmosphere, ocean and surface ground water. Geophys Res Lett 24(15):1867-1870
Dong D, Yunck T, Heflin M (2003) Origin of the International Terrestrial Reference Frame. J Geophys Res 108(B42200). doi:10.1029/2002JB002035

Farrell WE (1972) Deformation of the Earth by Surface Loads. Rev Geophys Space Phys10(3): 751-797

Fritsche M, Döll P, Dietrich R (2012) Global-scale validation of model-based load deformation of the Earth's crust from continental water mass and atmospheric pressure variations using GPS. J Geodyn 59-60:133-142

Griffiths J (2019) Combined orbits and clocks from IGS second reprocessing. J Geod 93:177-195. https://doi.org/10.1007/s00190-018-1149-8

Jiang WP, Li Z, van Dam T, Ding WW (2013) Comparative Analysis of Different Environmental Loading Methods and Their Impacts on the GPS Height Time Series. J Geod 87(7):687-703. doi: 10.1007/s00190-013-0642-3

Johnston G, Riddell A, Hausler G (2017) The International GNSS Service. In Teunissen, Peter JG, Montenbruck O. (Eds.). Springer Handbook of Global Navigation Satellite Systems (1st ed.): 967-982. Cham, Switzerland: Springer International Publishing. doi: 10.1007/978-3-319-42928-1

Kalnay E, Kanamitsu M, Kistler R, Collins W, Deaven D, Gandin L, Iredell M, Saha S, White G, Woollen J, Zhu Y, Chelliah M, Ebisuzaki W, Higgins W, Janowiak J, Mo K, Ropelewski C, Wang J, Leetmaa A, Reynolds R, Jenne R, Joseph D (1996) The NCEP/NCAR 40-year reanalysis project. Bull Amer Meteor Soc 77(3):437-471. DOI:10.1175/1520-0477

Kanamitsu M, Ebisuzaki W, Woollen J, Yang Sk, Hnilo JJ, Fiorino M, Potter GL (2002) NCEP-DOEAMIP-II reanalysis (R-2). Bull Amer Meteor Soc 83(83):1631-1643. doi: http://dx.doi.org/10.1175/BAMS-83-11-1631

King MA, Watson CS (2010) Long GPS coordinate time series: Multipath and geometry effects. J. Geophys. Res. 115 (B04403). doi:10.1029/2009JB006543

Li Z (2012) Research on the non-linear variation of GPS coordinate time series. Ph.D dissertation. Wuhan University

Li Z, Jiang WP, Ding WW, Deng LS, Peng LF (2014) Estimates of Minor Ocean Tide Loading Displacement and Its Impact on Continuous GPS Coordinate Time Series. Sensors, 14(3):5552-5572. doi:10.3390/s140305552

Li Z, van Dam T, Collilieux X, Altamimi Z, Rebischung P, Nahmani S (2015) Quality evaluation of the weekly vertical loading effects induced from continental water storage models. IAG 150 Years, International Association of Geodesy Symposia 143, Springer International Publishing Switzerland. 45-54. doi:10.1007/1345_2015_174

Li Z, van Dam T (2015) The phase 2 North America land data assimilation system (NLDAS-2) products for modeling water storage displacements for plate boundary observatory GPS stations. International Association of Geodesy Symposia, Springer International Publishing Switzerland. 1-9. doi:10.1007/1345_2015_176

Li Z, Chen W, Jiang WP, Deng LS, Yang RH (2018) The magnitude of diurnal/semidiurnal atmospheric tides (S1/S2) and their impacts on the continuous GPS coordinate time series. Remote Sens 10 1125. doi: 10.3390/rs10071125

Lindsay R, Wensnahan M, Schweiger A, Zhang J (2014) Evaluation of seven different atmospheric reanalysis products in the arctic. J Climate 27: 2588-2606, https://doi.org/10.1175/JCLI-D-13-00014.1

Lucchesi R (2012) File Specification for MERRA Products. GMAO Office Note No. 1 (Version 2.3), available from http://gmao.gsfc.nasa.gov/pubs/office_notes

Penna NT, Stewart MP (2003) Aliased tidal signatures in continuous GPS height time series. Geophys Res Lett 30 (23) 2184. doi:10.1029/2003GL018828

Penna NT, King MA, Stewart MP (2007) GPS height time series: Short-period origins of spurious long-period signals. J Geophys Res 112 (B02402). doi:10.1029/2005JB004047

Ray J, Altamimi Z, Collilieux X, van Dam T (2008) Anomalous harmonics in the spectra of GPS position estimates. GPS Solut 12:55-64

Ray J, Griffiths J, Collilieux X, Rebischung P (2013) Subseasonal GNSS positioning errors. Geophys Res Lett 40(22):5854-60. doi:10.1002/2013GL058160

Rienecker MM, Suarez MJ, Gelaro R, Todling R, Bacmeister J, Liu E, Bosilovich MG, Schubert SD, Takacs L, Kim GK, Bloom S, Chen JY, Collins D, Conaty A, da Silva A, Gu W, Joiner J, Koster RD, Lucchesi R, Molod A, Owens T, Pawson S, Pegion P, Redder CR, Reichle R, Robertson FR, Ruddick AG, Sienkiewicz
M, Woollen J (2011) MERRA-NASA's Modern-Era Retrospective Analysis for Research and Applications. J Climate 24(14): 3624-3648. doi:10.1175/JCLI-D-11-00015.1

Rebischung P, Altamimi Z, Ray J, Garayt B (2016) The IGS contribution to ITRF2014. J Geod 90(7): 611-630.
doi:10.1007/s00190-016-0897-6

Salstein DA, Ponte RM, Cady-Pereira K (2008) Uncertainties in atmospheric surface pressure fields from global analyses. J Geophys Res 113 (D14107). doi:10.1029/2007JD009531

Steigenberger P, Rothacher M, Schmid R, Rülke A, Fritsche M, Dietrich R, Tesmer V (2009) Effects of different antenna phase center models on GPS-derived reference frame. Geodetic Ref Frames Int Assoc Geodesy Symposia 134:83-88

Tregoning P, van Dam T (2005) Atmospheric pressure loading corrections applied to GPS data at the observation level. Geophys Res Lett 32(L22310). doi:10.1029/2005GL024104

Tregoning P, Watson C (2009) Atmospheric effects and spurious signals in GPS analyses. J Geophys Res 114:B09403. doi:10.1029/2009JB006344

Tregoning P, Watson C, Ramillien G, McQueen H, Zhang J (2009) Detecting hydrologic deformation using GRACE and GPS. Geophys Res Lett 36(L15401). doi:10.1029/2009GL038718

van Dam TM, Wahr J (1987) Displacements of the Earth's surface due to atmospheric loading: Effects on gravity and baseline measurements. J Geophys Res 92:1281-1286

van Dam T, Blewitt G, Heflin M (1994) Atmospheric pressure loading effects on global positioning system coordinate determinations. J Geophys Res 99(B12):23939-23950

van Dam TM, Wahr J, Chao Y, Leuliette E (1997) Predictions of crustal deformation and of geoid and sea-level variability caused by oceanic and atmospheric loading. Geophys J Int 129(3):507-517

van Dam T, Wahr J, Milly PCD, Shmakov AB, Blewitt G, Lavallée D, Larson KM (2001) Crustal displacements due to continental water loading. Geophys Res Lett 28(4):651-654

van Dam T, Altamimi Z, Collilieux X, Ray J (2010) Topographically induced height errors in predicted atmospheric loading effects. J Geophys Res 115(B07415). doi:10.1029/2009JB006810

van Dam T, Collilieux X, Wuite J, Altamimi Z, Ray J (2012) Nontidal ocean loading: amplitudes and potential effects in GPS height time series. J Geod 86(11):1043-1057. doi:10.1007/s00190-012-0564-5

Wessel P, Smith WHF (2018) The Generic Mapping Tools. Version 4.5.18 Technical reference and cookbook. Laboratory for satellite altimetry, NOAA/NESDIS/STAR

Wijaya DD, Böhm J, Karbon M, Krásná H, Schuh H (2013) Atmospheric pressure loading, in J. Böhm and H. Schuh (eds): Atmospheric Effects in Space Geodesy. Springer Verlag, pp. 137-157,
doi:10.1007/978-3-642-36932-2_4

Williams SDP, Penna NT (2011) Non-tidal ocean loading effects on geodetic GPS heights. Geophys Res Lett
38(L09314). doi:10.1029/2011GL046940

Yuan P, Li Z, Jiang WP, Ma YF, Chen W, Sneeuw N (2018) Influences of environmental loading corrections on the nonlinear variations and velocity uncertainties for the reprocessed global positioning system height time series of the crustal movement observation network of China. Remote Sens 10(6) 958. doi: 10.3390/rs10060958

Appendix

To help users deal with the non-integer constants files, here we give a simple step of how to generate integer
land/sea mask for MERRA and MERRA-2 datasets:

(1) Download the MERRA constant file with short name as MAC0NXASM.

(2) Create a local binary file that has integer masks for land and ocean in MERRA using gradsdap script. Then, generate a descriptor file that can be used to open the above binary integer mask file in grads.

(3) Write the land and sea variables from the above descriptor file into separate files in netcdf format using grads command sdfwrite<fnref n="14" />. Note that for using the command sdfwrite, users should install a fairly recent version of GrADS, e.g., version 2.0.a3 at least.

(4) Combine the separate land and sea netcdf files to generate the MERRA integer land/sea mask.

(5) Download the MERRA-2 constant file with short name as M2C0NXASM, repeat the above procedures from (1) to (4) to generate the MERRA-2 integer land-sea mask

After implementing the above 5 steps, we can use the obtained integer land/sea mask for modeling the ATML

<footnote n="14">14 http://www.iges.org/grads/gadoc/gradcomdsdfwrite.html</footnote>
effects using both MERRA and MERRA-2 PS grids. When creating the binary file from MERRA and MERRA-2
that containing integer masks for land and ocean, a fractional cutoff value is used to determine that if a point is
classified land or sea. Here, we choose the cutoff value as 0.25 for global study (from Mike Bosilovich's
cookbook). For regional studies where the land and ocean are intended to be discretized, more attempts should be
examined to determine an appropriate cutoff value.