The Impact of Di ﬀ erent Ocean Tide Loading Models on GNSS Estimated Zenith Tropospheric Delay Using Precise Point Positioning Technique

: Global navigation satellite systems (GNSSs) have become an important tool to derive atmospheric products, such as the total zenith tropospheric delay (ZTD) and precipitable water vapor (PWV) for weather and climate studies. The ocean tide loading (OTL) e ﬀ ect is one of the primary errors that a ﬀ ects the accuracy of GNSS-derived ZTD / PWV, which means the study and choice of the OTL model is an important issue for high-accuracy ZTD estimation. In this study, GNSS data from 1 January 2019 to 31 January 2019 are processed using precise point positioning (PPP) at globally distributed stations. The performance of seven widely used global OTL models is assessed and their impact on the GNSS-derived ZTD is investigated by comparing them against the ZTD calculated from co-located radiosonde observations. The results indicate that the inclusion or exclusion of the OTL e ﬀ ect will lead to a di ﬀ erence in ZTD of up to 3–15 mm for island stations, and up to 1–2 mm for inland stations. The di ﬀ erence of the ZTD determined with di ﬀ erent OTL models is quite small, with a root-mean-square (RMS) value below 1.5 mm at most stations. The comparison between the GNSS-derived ZTD and the radiosonde-derived ZTD indicates that the adoption of OTL models can improve the accuracy of GNSS-derived ZTD. The results also indicate that the adoption of a smaller cuto ﬀ elevation, e.g., 3 ◦ or 7 ◦ , can signiﬁcantly reduce the di ﬀ erence between the ZTDs determined by GNSS and radiosonde, when compared against a 15 ◦ cuto ﬀ elevation. Compared to the radiosonde-derived ZTD, the RMS error of GNSS-derived ZTD is approximately 25–35 mm at a cuto ﬀ elevation of 15 ◦ , and 15–25 mm when the cuto ﬀ elevation is set to 3 ◦ . about 3–15 mm in the GNSS-derived ZTD with the PPP at the coastal and island stations. Therefore, an optimal OTL model must be used in the estimation of ZTD with the GNSS PPP technique, especially for the coastal and island stations. The comparison between di ﬀ erent OTL models indicated that the vertical displacement calculated using di ﬀ erent OTL models can reach 10 mm at island stations, while for inland stations the di ﬀ erence is only about 1–2 mm. The study also indicates that the di ﬀ erence between the ZTDs derived from di ﬀ erent OTL models is quite small and the inclusion of an OTL model in the GNSS data processing will improve the accuracy of the ZTD estimation. With the radiosonde-derived ZTD as a reference, the accuracy of GNSS-derived ZTDs was studied with di ﬀ erent cuto ﬀ elevations and the results show that the RMS error is approximately 25–35 mm at a cuto ﬀ elevation of 15 ◦ , while the RMS error was only 15–25 mm with a cuto ﬀ elevation of 3 ◦ at all stations. Therefore, a low cuto ﬀ elevation (e.g., 3 ◦ or 7 ◦ ) is helpful in reducing the RMS error of the ZTD estimates and thus is recommended when estimating ZTD with GNSS.


Introduction
Atmospheric products derived from the global navigation satellite system (GNSS), such as zenith tropospheric delay (ZTD) and precipitable water vapor (PWV), have been widely used in weather and climate studies [1][2][3][4][5]. Since there is a strong correlation between ZTD, height coordinates and ocean tide loading (OTL) effects during GNSS data processing, it is important to take the OTL effect into account when estimating ZTD, especially for stations located near the ocean. Otherwise, the OTL effect will map into the ZTD and station height estimates [6,7]. As suggested by previous studies, neglecting the OTL effect could lead to an error of up to 50 mm in station height estimates. Dach and Dietrich [8] pointed out that the difference between the ZTD estimated with inclusion or exclusion of the OTL model may reach up to 10 mm with a double-difference processing strategy. As suggested by the World Meteorology Organization (WMO) the threshold for using ZTD in the weather forecasting should be 15 mm [9]. Therefore, to meet this requirement and to further improve the accuracy of the GNSS-derived ZTD, it is very important to carefully investigate the impact of different data processing strategies on the estimation of the ZTD, including the OTL effects.
As the elastic response of the earth's crust to ocean tides, the OTL effect is dominated by semi-monthly periods, with amplitudes of several centimeters. Therefore, an accurate OTL model is needed in GNSS data processing to achieve a high-accuracy result for a variety of glaciological and oceanographic applications [10,11]. Several OTL models have been proposed thanks to the rapid development of satellite altimetry measuring techniques from several missions including TOPEX, Jason1, Jason2 and Grace. The main difference between different OTL models is the model construction method and the dataset adopted. For example, GOT4.7 is a global ocean tide model, with a resolution of 0.5 • × 0.5 • , that used the European Remote Sensing (ERS) altimetry data and the Geosat Follow-On (GFO) data in shallow and polar seas [12,13]. The Finite Element Solutions (FES) series was produced by the French tidal group. It has many versions, such as FES2014, FES2012, FES2008, FES2004, and FES99. Although FES2014, FES2012 and FES2008 have a higher accuracy and resolution than FES2004, FES2004 is the conventional model recommended by International Earth Rotation Service (IERS) conventions [14]. The FES2004 used a refined mesh with a resolution of 0.125 • × 0.125 • . This ocean tide model assimilated the tide gauge data, Topex/Poseidon and ERS altimetry data [15].
Because these OTL models were constructed with different data, methods and spatial resolution, their accuracy and impact on the ZTD estimation might be different in different regions. Therefore, in this study, the performance of these OTL models and their impact on the ZTD estimation at different regions was studied in detail. Although some studies have focused on the difference between GNSS-ZTDs estimated with and without OTL models [8,16,17], few have been conducted on the effect of different OTL models on the ZTD estimation with GNSS [18,19]. In addition, the aforementioned studies were all conducted using double-difference data processing strategies, and few studies have been conducted with precise point positioning (PPP) [20]. As it moves forward in the high-accuracy orbit and clock products, PPP has become an important method in ZTD estimation, since it is much more efficient and flexible compared to double-difference data processing when estimating real-time and high-rate ZTD for a large network. The International GNSS Service (IGS) final troposphere products were generated with PPP [21]. In this study, the difference of using different OTL models in GNSS-derived ZTD was studied using the PPP strategy; the difference between OTL models in different areas has also been studied.
The structure of the paper as follows. Section 2 describes the data and the adopted methodology. Section 3 provides the details on the results. Section 4 provides a discussion of this study. The main findings and conclusions are shown in Section 5.

Data and Methodologies
In this study, the GNSS data from 1 January 2019 to 31 January 2019 are processed using PPP at 21 globally distributed GNSS stations. To investigate the impact of the OTL effect on the stations at different regions, 21 stations which are located on islands, coastal areas or on inland areas that are far from the sea were selected. To investigate the impact of different OTL models on the estimation of ZTD, the ZTD estimated without using OTL is compared to that estimated with seven commonly used OTL models, including DUT10 [22], EOT11a [23], FES2004 [14,15,[24][25][26][27], GOT4.7 [12,13], Hamtide11a [28], Nao.99b [29] and TPXO7.2 [30,31]. For the 10 GNSS stations that each have a co-located radiosonde site, their GNSS-derived ZTDs were compared to their radiosonde-derived ZTDs; more details can be found in Section 2.1. In Section 2.2, the OTL models and the method used to calculate OTL are analyzed. The ZTD calculation, as derived from either radiosonde or GNSS, is shown in Sections 2.3 and 2.4.

GNSS and Radiosonde Data
In this study, GNSS observations at 21 globally distributed GNSS stations, including 4 inland stations, 7 island stations, and 10 coastal stations, were taken from 1 January 2019 to 31 January 2019 and processed to analyze the impact of different OTL models on the estimation of ZTD. Figure 1 indicates the 21 globally distributed GNSS stations. Among the 21 stations, 10 stations had a co-located radiosonde site within a distance of 5 km. A comparison of the GNSS and radiosonde derived ZTD was measured at the ten co-located stations; the results are presented in Table 1. The ZTDs at these 10 radiosonde stations were calculated as references to investigate the accuracy of the GNSS-derived ZTD estimates using different OTL models.

GNSS and Radiosonde Data
In this study, GNSS observations at 21 globally distributed GNSS stations, including 4 inland stations, 7 island stations, and 10 coastal stations, were taken from 1 January 2019 to 31 January 2019 and processed to analyze the impact of different OTL models on the estimation of ZTD. Figure 1 indicates the 21 globally distributed GNSS stations. Among the 21 stations, 10 stations had a colocated radiosonde site within a distance of 5 km. A comparison of the GNSS and radiosonde derived ZTD was measured at the ten co-located stations; the results are presented in Table 1. The ZTDs at these 10 radiosonde stations were calculated as references to investigate the accuracy of the GNSSderived ZTD estimates using different OTL models. In this study, the observations files at four stations (ALRT, PBRI, POVE, and NOVM) were downloaded from the site ftp://garner.ucsd.edu/pub/rinex/ and the others were downloaded from the site ftp://cddis.gsfc.nasa.gov/. For these 10 radiosonde stations, there are nine stations located within a distance less than 2 km from their co-located GNSS stations and one station with a distance of 3.682 km. It should be noted that, for some extreme weather conditions, this horizontal distance may result in a difference in ZTD at the GNSS site versus the nearby radiosonde site, due to the spatial variation of water vapor. However, as suggested by previous studies [32][33][34][35], it is acceptable in most cases to neglect this horizontal distance and regard these sites as co-located. The radiosonde data is provided by the Integrated Global Radiosonde Archive (IGRA) [36][37][38]. More details regarding the GNSS data processing and ZTD calculation with radiosonde data can be found in the following subsections. In this study, the observations files at four stations (ALRT, PBRI, POVE, and NOVM) were downloaded from the site ftp://garner.ucsd.edu/pub/rinex/ and the others were downloaded from the site ftp://cddis.gsfc.nasa.gov/. For these 10 radiosonde stations, there are nine stations located within a distance less than 2 km from their co-located GNSS stations and one station with a distance of 3.682 km. It should be noted that, for some extreme weather conditions, this horizontal distance may result in a difference in ZTD at the GNSS site versus the nearby radiosonde site, due to the spatial variation of water vapor. However, as suggested by previous studies [32][33][34][35], it is acceptable in most cases to neglect this horizontal distance and regard these sites as co-located. The radiosonde data is provided by the Integrated Global Radiosonde Archive (IGRA) [36][37][38]. More details regarding the GNSS data processing and ZTD calculation with radiosonde data can be found in the following subsections.

OTL Models
Most global ocean tide models provide 11 main tidal constituents to reflect the variations in sea surface height [16,39]. The accuracy of the OTL value at a GNSS site calculated with an OTL model is mainly dependent on the accuracy of the adopted OTL and its spatial resolution. As shown in Table 2, a comparison between the different OTL models is presented. The accuracy of an OTL model is normally lower in a coastal area than in an inland area because the complex environment in the coastal area makes it very difficult to obtain high-accuracy satellite altimetry data for model construction [40]. Equation (1) provides a method to compute the OTL by Farrell [42]. It uses an ocean tide model and a Green's function, written as: where ϕ and λ are the latitude and longitude, respectively; ρ and Ω represent the mean density of ocean water and the domain of integration, respectively; h and dS are the elevation and surface area of sea surface, respectively; (ϕ, λ) and (ϕ , λ ) are the estimation point and loading point, respectively; α represents the angular distance and L is the displacement vector which can be calculated with SPOTL software [43]. G denotes the Green's function, which can be calculated using Equation (2): where m E and g(a) are the mass and gravitational acceleration of the Earth, respectively; h n and k n represent the Love numbers of loading which were determined by Farrell; and P n denotes the Legendre function. Regarding the Green's function, there is a small difference in the use of different Green's functions to calculate the OTL. Generally, the difference between the OTL values calculated with different Green's functions is less than 2% [44]. Equation (1) can be rewritten as a displacement vector of vertical, north-south, and east-west components at a grid point and at the epoch t, as [45]: where ω i and ϕ i represent the frequency and astronomical argument of the different tidal constituents, represent the Greenwich phase lags; and n is the number of main tidal constituents.

ZTD Calculation with Radiosonde
The radiosonde observations were used to calculate the ZTD as a reference for assessing the accuracy of the GNSS-estimated ZTD and the effect of different OTL models on the ZTD estimation with GNSS. As shown in Table 1, all of the GNSS stations and their co-located radiosonde sites were within 5 km of each other and had a height difference under 50 m. The ZTD was estimated as the integral of refractivity through the atmosphere in the zenith direction and Equation (4) was used to determine the refractivity with observations at each pressure level. Refractivity is written as [32]: where N is the refractivity index; K 1 = 77.6890K/hPa, K 2 = 22.13K/hPa, and K 3 = 3.739 × 10 5 K/hPa are empirically determined constants; e expresses the pressure of water vapor; and T and P d are the temperature and the dry air pressure, respectively. With the refractivity determined using Equation (4), the ZTD is calculated as: where ∆h is the height difference between the two pressure levels, and i is the pressure level. When the GNSS site is below the radiosonde site, a vertical extrapolation procedure to extrapolate the meteorological parameters from the radiosonde height to the GNSS height is conducted to make sure that the radiosonde-derived ZTD is at the same height as the GNSS. If the GNSS site is above the radiosonde site, then the meteorology parameter at the GNSS site is interpolated with the radiosonde data at two nearest pressure levels. After this extrapolation or interpolation procedure, the ZTD at the GNSS site was calculated with Equations (4) and (5), using radiosonde observations [5].

ZTD Calculation with GNSS
Compared to the double-difference solution, PPP is a faster and more efficient means of estimating geodetic parameters such as the coordinate and tropospheric delay. The RTKLIB software is used in this study for the GNSS data processing with PPP. The latest version (V2.4.3) of RTKLIB is adopted in this study, which can be downloaded at https://github.com/tomojitakasu/RTKLIB/tree/rtklib_2.4.3 [46]. In RTKLIB, the tropospheric zenith hydro-static delay is given by the Saastamoinen model. In regards to the mapping function, both the Neill mapping function (NMF) and global mapping function (GMF) [47,48] are available. As suggested by previous studies, although the Vienna mapping functions (VMF) have a better performance than the others, the difference among these commonly used mapping functions is quite small [49]. The gradient model parameters for the north and east are estimated together with the ZTD with the form proposed by Macmillan [50]. RTKLIB employs an extended Kalman filter (EKF) to estimate coordinates, ZTD and other parameters. The ZTD results estimated during the convergence period were removed in the following study. In this study, the IGS final products, including the satellite orbits and clock, and the earth rotation and differential code bias (DCB) [51], were used for the estimation of ZTD. Since RTKLIB does not support the inclusion of Remote Sens. 2020, 12, 3080 6 of 18 an atmospheric pressure loading model, that effect was not considered in this study. As suggested by previous studies [52], a random walk stochastic model is used; the constraint of ZTD is set to 20 mm/( √ h) and the standard deviation of the carrier phase bias in the process noise is 10 -4 cycle/( √ s). The processing strategies in Table 3 were adopted in this study.

Results
The effects of different OTL models and elevation cutoffs on the ZTD estimates were studied in this section. In Section 3.1, the magnitudes and variations of the OTL effects are studied at a global scale. In Section 3.2, the GNSS-derived ZTD using the seven OTL models and calculated without OTL models in different regions are compared. In Section 3.3, the difference between the GNSS-derived ZTD and the co-located radiosonde data calculated ZTD are investigated and the impact of the elevation cutoff on the estimation of ZTD is also studied.

Temporal and Spatial Characteristics of the OTL Models
Before investigating the impact of different OTL models on the estimation of ZTD, the temporal and spatial characteristics of the OTL models were investigated. In this study, SPOTL software was used to compute the magnitudes of the site displacements due to OTL effects on model grids, and then bilinear interpolation was adopted to calculate the magnitude of the OTL effect at any specific site. ZTD estimation is mainly influenced by the OTL effects in the vertical direction; those effects have a magnitude and variation that are much larger than the corresponding values in the horizontal direction. As shown in Figure 2, the magnitudes of the OTL effect in the vertical direction exhibit an obvious variation with location and time during January 2019, with a maximum value of around 60 mm. The OTL effect in the coastal region is larger than in the inland area, which only has a magnitude smaller than 10 mm. It can be seen in Figure 2 that there is an obvious semi-monthly variation of OTL effects.
Remote Sens. 2020, 12, 3080 7 of 18 have a magnitude and variation that are much larger than the corresponding values in the horizontal direction. As shown in Figure 2, the magnitudes of the OTL effect in the vertical direction exhibit an obvious variation with location and time during January 2019, with a maximum value of around 60 mm. The OTL effect in the coastal region is larger than in the inland area, which only has a magnitude smaller than 10 mm. It can be seen in Figure 2 that there is an obvious semi-monthly variation of OTL effects.   Figure 3 is an inland station located in Novosibirsk, Russia; the BRMU station in Figure 4 is an island station located on Bermuda island, United Kingdom; and the NYA1 is a coastal station located on Svalbard, Norway. As shown in Figures 3-5, it can be seen that the magnitudes of the OTL effects are much greater at the coastal and island stations than they are at the inland stations. The variation of the vertical displacements   Figure 3 is an inland station located in Novosibirsk, Russia; the BRMU station in Figure 4 is an island station located on Bermuda island, United Kingdom; and the NYA1 is a coastal station located on Svalbard, Norway. As shown in        To investigate the differences among the OTL models, the widely used FES2004 model was used as the reference to calculate the difference between FES2004 and the other OTL models in the calculation of OTL effect corrections at the aforementioned three stations. As shown in Figures 6 and 7, the biggest difference between each OTL model was found at the island station BRMU, with a root-mean-square (RMS) error of around 4-5 mm and a maximum value of around 10 mm. The OTL models show good agreement with each other at the inland station NOVM, with an RMS error of about 1 mm. The RMS error at the coastal station NYA1 was about 3 mm and an obvious temporal variation can be seen in the differences between OTL models at this station. In addition, the EOT11a model is in better agreement with the FES2004 model than the other models. mean-square (RMS) error of around 4-5 mm and a maximum value of around 10 mm. The OTL models show good agreement with each other at the inland station NOVM, with an RMS error of about 1 mm. The RMS error at the coastal station NYA1 was about 3 mm and an obvious temporal variation can be seen in the differences between OTL models at this station. In addition, the EOT11a model is in better agreement with the FES2004 model than the other models.   about 1 mm. The RMS error at the coastal station NYA1 was about 3 mm and an obvious temporal variation can be seen in the differences between OTL models at this station. In addition, the EOT11a model is in better agreement with the FES2004 model than the other models.

GNSS-Derived ZTD
GNSS-derived ZTDs were calculated at 21 stations using the PPP data processing strategy with or without the OTL models. Figures 8 and 9 show the difference between the ZTDs calculated with different OTL models and without OTL models. As shown in Figure 8, the RMS of differences between ZTDs calculated with and without OTL models for the island and coastal stations is between 3-15 mm at most stations. The largest difference was found at the island station FAA1, with an RMS error of about 30 mm. For the inland stations, the differences between the ZTD determined with and without the OTL models were relatively small, with an RMS error of 1-2 mm. or without the OTL models. Figures 8 and 9 show the difference between the ZTDs calculated with different OTL models and without OTL models. As shown in Figure 8, the RMS of differences between ZTDs calculated with and without OTL models for the island and coastal stations is between 3-15 mm at most stations. The largest difference was found at the island station FAA1, with an RMS error of about 30 mm. For the inland stations, the differences between the ZTD determined with and without the OTL models were relatively small, with an RMS error of 1-2 mm.   different OTL models and without OTL models. As shown in Figure 8, the RMS of differences between ZTDs calculated with and without OTL models for the island and coastal stations is between 3-15 mm at most stations. The largest difference was found at the island station FAA1, with an RMS error of about 30 mm. For the inland stations, the differences between the ZTD determined with and without the OTL models were relatively small, with an RMS error of 1-2 mm.   As shown in Figure 9, the ZTDs calculated with different OTL models show slight differences. To further investigate the impact of the different OTL models on the estimation of ZTD, a comparison between the ZTD calculated with the commonly used FES2004 model and other models was conducted. Figures 10 and 11 show that the largest difference between each OTL model was found at coastal stations, with an RMS error of about 0.5-1.5 mm. The OTL models show good agreement with each other at the inland stations, with an RMS error of about 0.1-0.2 mm. The RMS error at the island stations was about 0.5-1 mm. In addition, the difference shows a twice-monthly variation at the BRMU station.
between the ZTD calculated with the commonly used FES2004 model and other models was conducted. Figures 10 and 11 show that the largest difference between each OTL model was found at coastal stations, with an RMS error of about 0.5-1.5 mm. The OTL models show good agreement with each other at the inland stations, with an RMS error of about 0.1-0.2 mm. The RMS error at the island stations was about 0.5-1 mm. In addition, the difference shows a twice-monthly variation at the BRMU station.   conducted. Figures 10 and 11 show that the largest difference between each OTL model was found at coastal stations, with an RMS error of about 0.5-1.5 mm. The OTL models show good agreement with each other at the inland stations, with an RMS error of about 0.1-0.2 mm. The RMS error at the island stations was about 0.5-1 mm. In addition, the difference shows a twice-monthly variation at the BRMU station.

Radiosonde
To further investigate the impact of the OTL models on the accuracy of the GNSS-derived ZTDs, the radiosonde data at 10 stations co-located with GNSS stations were used to calculate the ZTDs, which were used as a reference in this study. More details regarding the 10 radiosonde and GNSS stations are listed in Table 1. The ZTDs derived from radiosonde data were used as the reference to evaluate the accuracy of the GNSS-derived ZTDs. As presented in Table 4, the ZTDs derived from different OTL models show a good agreement with the radiosonde-derived ZTD, and all have a better agreement than the ZTD derived without the OTL models. For example, at the BRMU station, the RMS error of the GNSS-ZTD derived without the OTL model was 25.68 mm and it decreased to about 19 mm when an OTL model was used during data processing. Table 4. The statistic presented in this table is the difference between the GNSS-derived ZTD and the ZTD derived from the co-located radiosonde station (as a reference), for GNSS stations with a co-located radiosonde station, with comparisons across OTL models. It should be noted that the cutoff elevation adopted in the above data processing was 10 • . As there is a large correlation between height, ZTD and OTL displacement, it may be useful to lower the elevation cutoff to reduce their correlation and increase the accuracy of the GNSS-derived ZTD. In this study, the GNSS-derived ZTD was calculated with 6 cut-off angles, i.e., 3 • , 5 • , 7 • , 10 • , 12 • and 15 • . Results for each were compared and the ZTDs derived from them were compared with radiosonde data. Figure 12 shows the RMS error between the GNSS-derived ZTD with different cut-off elevations and radiosonde-derived ZTDs. The results show that the RMS error of the GNSS-derived ZTDs decreases when the cutoff elevation is lower. For example, the GNSS-derived ZTDs had the smallest RMS error at the cutoff elevation of 7 • , and the largest RMS error at 15 • at the NOVM station. For the BRMU station and NYA1 station, the best agreement between the ZTD calculated with GNSS and that calculated with radiosonde is found at the elevation of 3 • or 7 • . The RMS error of GNSS-ZTD was approximately 28 mm at a cutoff elevation of 15 • and only 20 mm at the cutoff elevation of 3 • at the BRMU station. This improvement can be found for the experiments with all of the studied OTL models.  Figure 13 shows the RMS differences between the GNSS-derived ZTD and the radiosondederived ZTD at 10 stations with six cutoff angles using the FES2004 model. The results indicated that the adoption of a smaller cutoff elevation, e.g., 3° or 7°, significantly reduces the difference between the ZTDs determined with GNSS versus radiosonde, compared with using a 15° cutoff elevation. Compared to the radiosonde-derived ZTD, the RMS error of the GNSS-derived ZTD was approximately 25-35 mm at a cutoff elevation of 15°, and only 15-25 mm with a cutoff elevation of  Figure 13 shows the RMS differences between the GNSS-derived ZTD and the radiosonde-derived ZTD at 10 stations with six cutoff angles using the FES2004 model. The results indicated that the adoption of a smaller cutoff elevation, e.g., 3 • or 7 • , significantly reduces the difference between the ZTDs determined with GNSS versus radiosonde, compared with using a 15 • cutoff elevation. Compared to the radiosonde-derived ZTD, the RMS error of the GNSS-derived ZTD was approximately 25-35 mm at a cutoff elevation of 15 • , and only 15-25 mm with a cutoff elevation of 3 • at all stations. At the inland station, POVE, the RMS error was reduced from 25 mm to 10 mm when the cutoff elevation changed from 15 • to 7 • . Figure 12. The RMS differences between the GNSS-derived ZTD and radiosonde stations with six cutoff angles and different OTL models at the (a) NOVM station, (b) BRMU station, and (c) NYA1 station. Figure 13 shows the RMS differences between the GNSS-derived ZTD and the radiosondederived ZTD at 10 stations with six cutoff angles using the FES2004 model. The results indicated that the adoption of a smaller cutoff elevation, e.g., 3° or 7°, significantly reduces the difference between the ZTDs determined with GNSS versus radiosonde, compared with using a 15° cutoff elevation. Compared to the radiosonde-derived ZTD, the RMS error of the GNSS-derived ZTD was approximately 25-35 mm at a cutoff elevation of 15°, and only 15-25 mm with a cutoff elevation of 3° at all stations. At the inland station, POVE, the RMS error was reduced from 25 mm to 10 mm when the cutoff elevation changed from 15° to 7°.

Discussion
GNSS has become an established technique for atmospheric remote sensing in the last decade. With the development of high-accuracy IGS products, including the orbits and clocks, PPP has been widely used to retrieve the ZTD because it is much more efficient and faster than the double-difference method. Since there is a strong correlation among ZTD, height coordinates and OTL effects during GNSS data processing, it is important to take the OTL effect into account when estimating ZTD with PPP, especially for stations located near the ocean. Otherwise, the OTL effect will map into the ZTD and station height estimation results.
In this study, GNSS data from 1 January 2019 to 31 January 2019 are processed using PPP at globally distributed stations. The performance of seven ocean tide models was assessed and their impact on the GNSS-derived ZTD was investigated by comparing them against the ZTD calculated from co-located radiosonde observations. As shown in Figures 8 and 9, the RMS error of the differences between ZTD calculated with and without OTL models is about 3-15 mm at the most island and coastal stations. The largest difference was found at the island station FAA1, which had an RMS error of about 30 mm. For inland stations, the differences between the ZTD determined with and without OTL models were relatively small, with an RMS of 1-2 mm. This means that in the estimation of ZTD with GNSS, the ocean tide models must be considered when the GNSS station is located on an island or in a coastal area. The differences of the ZTDs determined with different OTL models are shown in Figures 10 and 11, and they are seen to be quite small. The largest difference between each OTL model was found at coastal stations, with an RMS of about 0.5-1.5 mm. The OTL models show a good agreement with each other at the inland station, with an RMS error of the differences of about 0.1-0.2 mm and 0.5-1 mm at most island stations. The OTL has semi-annual periods, so the OTL effect on the ZTD estimation might be different in different seasons. These seasonal impacts of the OTL effects on the ZTD estimation need to be investigated with a longer period of data.
To further investigate the impact of OTL models on the accuracy of the estimated ZTD with GNSS, radiosonde data at 10 co-located stations were used to calculate ZTDs, which were used as a reference. The comparison results presented in Table 4 show that the RMS error of the differences between ZTD determined with GNSS versus radiosonde is 10-30 mm at all stations. The ZTDs derived from different OTL models show a good agreement and all have a better agreement with the radiosonde-derived ZTD than the ZTD derived without OTL models. The comparison between the ZTD determined with different cutoff elevations shows that the adoption of a lower cutoff elevation (e.g., 3 • or 7 • ) significantly improves the accuracy of the ZTD compared with a 15 • cutoff elevation, as shown in Figures 12 and 13. The RMS error of GNSS-derived ZTD was approximately 25-35 mm at a cutoff elevation of 15 • , and only 15-25 mm with a cutoff elevation of 3 • , at all stations. For example at the inland station POVE, the RMS error was reduced from 25 mm to 10 mm when the cutoff elevation changed from 15 • to 7 • ; at the island station BRMU, the RMS error was reduced from 35 mm to 20 mm when the cutoff elevation changed from 15 • to 3 • .

Conclusions
In this study, the impact of the OTL models on the estimation of GNSS-ZTD with PPP has been comprehensively studied. The results indicate that the vertical OTL displacement is very large in the coastal and island regions and the omission of the OTL effect will lead to a difference of about 3-15 mm in the GNSS-derived ZTD with the PPP at the coastal and island stations. Therefore, an optimal OTL model must be used in the estimation of ZTD with the GNSS PPP technique, especially for the coastal and island stations. The comparison between different OTL models indicated that the vertical displacement calculated using different OTL models can reach 10 mm at island stations, while for inland stations the difference is only about 1-2 mm. The study also indicates that the difference between the ZTDs derived from different OTL models is quite small and the inclusion of an OTL model in the GNSS data processing will improve the accuracy of the ZTD estimation. With the radiosonde-derived ZTD as a reference, the accuracy of GNSS-derived ZTDs was studied with different cutoff elevations and the results show that the RMS error is approximately 25-35 mm at a cutoff elevation of 15 • , while the RMS error was only 15-25 mm with a cutoff elevation of 3 • at all stations. Therefore, a low cutoff elevation (e.g., 3 • or 7 • ) is helpful in reducing the RMS error of the ZTD estimates and thus is recommended when estimating ZTD with GNSS.