The Magnitude of Diurnal/Semidiurnal Atmospheric Tides (S1/S2) and Their Impacts on the Continuous GPS Coordinate Time Series

: The site displacement due to diurnal and semidiurnal atmospheric tides (S1/S2) is often neglected in the routine precise GPS data processing. We recall the S1/S2 modeling method and show the magnitude of the S1/S2 tides under the Center of Mass (CM) frame. The results show that the annual amplitudes caused by both S1 and S2 tides exceed 1 mm for stations near the equator. The impact of S1/S2 on the 24-h Global Positioning System (GPS) solution is typically at the sub-mm level, and the scatter of the S1/S2 caused displacement difference increases with the decreasing latitude for northern hemisphere stations, among which the maximum Standard Deviation (STD) reaches up to 1.5 mm, 1 mm and 0.7 mm for the Up, East and North components, respectively, at low-latitude stations. We also ﬁnd that 65% of the stations’ vertical velocity change caused by S1/S2 is larger than 1%, among which the maximum velocity variation rate reaches more than 40% for some coastal/island stations, while stations with the weighted root mean square reduced account for 65%, 39%, and 38% for the up, east, and north components respectively, in particular for most coastal/island stations. Furthermore, the S1/S2 correction could partly reduce the annual and the semi-annual signals of the global stacked vertical component together with the semi-annual amplitude of the east component. The power of some anomalous harmonics of 1.04 cycle per year also decreased a lot. These results further prove the beneﬁts of S1/S2 correction in the precise GPS data processing.


Introduction
The diurnal heating of the atmosphere causes the oscillation of the surface pressure at the diurnal (S1), semi-diurnal (S2) periods, as well as other higher frequency harmonics [1]. These "atmospheric tides" induce the periodic deformation of the Earth's surface [2]. The amplitude of the vertical deformation caused by S1 and S2 atmospheric tides (hereafter, S1/S2) reached several millimeters, which is equal to the magnitude of some ocean tidal loading components, thus, they must be precisely modeled for geodynamic studies, for example, sea level change or gravity [3][4][5]. The International Earth Rotation and Reference Systems Service (IERS) Conventions 2010 also suggests that the resulting deformation caused by S1/S2 should be considered in the station motion model [6].
In the most recent studies, Tregoning and Watson found that the neglect of the S1/S2 would introduce anomalous propagated signals with periods that closely match the Global Positioning System (GPS) draconitic annual (~351.4 days) and semiannual period (~175.7 days) [7,8]. Being close to the orbital period of the GPS satellites, modeling of the S2 effect is especially important in order to minimize aliasing [7]. Jiang et al. proposed that the S1/S2 could be one of the main reasons that cause the annual motion of IGS stations inside China, especially for the vertical annual motion of stations in the central and southern regions [9]. In the geodetic field, however, in particular, for real-time and precise point positioning (PPP) applications, the S1/S2 effects did not come into consideration for routine data reduction [10]. Additionally, the most recent second International GNSS Service (IGS) data reprocessing campaign did not include the S1/S2 correction [11]. For regional studies, until most recently, the routine GNSS data processing for the Crustal Movement Observation Network of China (CMONOC) still did not consider the S1/S2 effects [12]. What is the effect of S1/S2 on the spectrum, especially the annual and semiannual variation of the global IGS stations' position time series? What is the magnitude of the S1/S2 on the 24-h GPS solution, long-term velocity together with the non-linear variations of global IGS stations? Furthermore, without using the existing precise GPS data processing software, for example, GAMIT [13], how do we implement the S1/S2 correction during our own self-developed GPS data processing? These are the focuses of this paper.
In this paper, we first recall the modeling method of the S1/S2 correction during the precise GPS data processing, determining the magnitude and spatial distribution of global IGS station's annual amplitudes caused by S1/S2 under the center of mass (CM) frame. We then reprocess the GPS data of 109 globally distributed IGS stations spanning from January 1999 to December 2003 with state of the art models according to IERS Conventions 2010 based on GAMIT software. The reprocessing includes two runs, one is with the S1/S2 correction, while the other is without S1/S2 correction. Finally, quantitative analyses have been done on two sets of GPS coordinate time series in both time and frequency domains to evaluate the contributions of S1/S2 to global GPS coordinate time series. The results of this paper may give readers a thorough understanding of the S1/S2 effects during GPS data processing. It also provides numerical support for interpreting crustal movement and other geophysical signals.

S1/S2 Modeling
The tidal signals in global pressure data can be problematic for several reasons, including the potential aliased sampling of the S2 tide, as well as the presence of various modeling or timing errors [14]. In 2003, Ray and Ponte presented an effective S1 and S2 atmospheric tidal model to handle the S1/S2 tides for the oceanographic and geodetic applications (hereafter referred to as the RP03 model) based on Reference [15], which has become a conventional recommendation to calculate the periodic motion of the stations caused by S1 and S2 atmospheric tides [3,6].
The S1 and S2 RP03 tidal model is derived from the European Centre for Medium-Range Weather Forecasts (ECMWF) operational global surface pressure fields. It provides a global grid of the 4 annual mean tidal pressure coefficients (cos S1, sin S1, cos S2, sin S2), with a spatial resolution of 1.125 • × 1.125 • . The surface displacement caused by S1/S2 could then be obtained through convolution between Farrell's elastic Green's functions under the specific reference frame and the RP03 model [16][17][18][19]. Please note that during the generation of the RP03 model, no ocean response to pressure is assumed [20]. This means that they did not consider an inverted barometer (IB) response, and just used the original value from global surface pressure grid for the ocean area to calculate the mean amplitude and phase of S1 and S2. We agree that this is a disadvantage for the RP03 model since, from our experience, IB response affects a lot to the coastal stations for the non-tidal atmospheric loading, and we think that it would have a similar effect on the air tides too. At the current stage, as geodetic users, we just use the RP03 model, but could not assess the errors of this model. In the future, we will try to generate our own S1/S2 model based on the IB response, and then make a detailed comparison with the existing RP03 model to find out the validity.
From October 2010, the Global Geophysical Fluid Center (GGFC), established by the International Earth Rotation and Reference Systems Service (IERS), began to provide the online S1/S2 tidal calculation service [20]. Users could obtain the 3-D amplitudes and phases, or the amplitudes of the sines and cosines for the S1 and S2 tides under both the center of mass (CM) and the center of the solid Earth (CE) frame by either specifying the stations' longitude and latitude together with the target reference frame to let the calculator do the convolution, or download the global gridded 3-D deformations and a Fortran program to interpolate the gridded deformation to any point of interest. The 3-D displacements at the sites for any epoch are determined using d S1 (u, n, e) = A S1 (u, n, e) × cos(ω 1 * t) + B S1 (u, n, e) × sin(ω 1 * t) d S2 (u, n, e) = A S2 (u, n, e) × cos(ω 2 * t) + B S2 (u, n, e) × sin(ω 2 * t) d T (u, n, e) = d S1 (u, n, e) + d S2 (u, n, e) where d S1 (u, n, e) and d S2 (u, n, e) represent the S1, S2 induced 3-D displacements for the Up (U), North (N) and East (E) components, d T (u, n, e) represents the total displacements of the S1 and S2 tides, A S1 , B S1 , A S2 , B S2 are the in-phase and out-of-phase amplitudes of the S1 and S2 tides in unit of millimeter (mm), t is in fractions of a universal time (UT1) day, ω 1 and ω 2 denote the frequencies of the S1 and S2 tides respectively, with ω 1 = 2πradians/day and ω 2 = 4πradians/day.
To have a global view of the magnitude and spatial characteristics of the S1 and S2 effects, we calculate the annual amplitudes and phases of the S1/S2 induced surface displacement under the CM frame for 109 evenly distributed global IGS stations using the above online tidal loading calculator (Supplementary, Figure S1). The spatial distribution of the selected IGS stations is shown in Figure 1 [21]. Figure 2 shows the displacement amplitudes caused by S1 and S2 tides under the CM frame for each station with respect to the latitude. Here, we should keep in mind that the amplitude and phase obtained has an uncertainty since the RP03 model itself definitely has an uncertainty in it. However, at the current stage, we could not assess the uncertainty as ECMWF only provides surface pressure data, and there is no precision information along with it. Correspondingly, the RP03 model, which comes from the operational ECMWF analysis, also provides no error information. Therefore, we could not evaluate or validate the precision of the calculated S1/S2 amplitude and phase. This is a disadvantage of using the RP03 model to model atmospheric tides. Nevertheless, it was the only choice until now and has been suggested as the convention [6]. In the future, with the development of environmental loading service (e.g., there exists a new surface pressure data with precision information, or a more operational atmospheric tidal model comes into existence), we may find an effective way to evaluate the uncertainty of the S1/S2 induced surface displacement.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 15 S1/S2 tidal calculation service [20]. Users could obtain the 3-D amplitudes and phases, or the amplitudes of the sines and cosines for the S1 and S2 tides under both the center of mass (CM) and the center of the solid Earth (CE) frame by either specifying the stations' longitude and latitude together with the target reference frame to let the calculator do the convolution, or download the global gridded 3-D deformations and a Fortran program to interpolate the gridded deformation to any point of interest. The 3-D displacements at the sites for any epoch are determined using d S1 (u,n,e) = A S1 (u, n,e) × cos(ω 1 * t) + B S1 (u,n,e) × sin(ω 1 * t) d S 2 (u,n,e) = A S 2 (u,n,e) × cos(ω 2 * t) + B S 2 (u, n,e) × sin(ω 2 * t) d T (u, n,e) = d S1 (u,n,e) + d S 2 (u,n,e) T d u n e represents the total displacements of the S1 and S2 tides, B are the in-phase and out-of-phase amplitudes of the S1 and S2 tides in unit of millimeter (mm), t is in fractions of a universal time (UT1) day, 1 ω and 2 ω denote the frequencies of the S1 and S2 tides respectively, with To have a global view of the magnitude and spatial characteristics of the S1 and S2 effects, we calculate the annual amplitudes and phases of the S1/S2 induced surface displacement under the CM frame for 109 evenly distributed global IGS stations using the above online tidal loading calculator (Supplementary, Figure S1). The spatial distribution of the selected IGS stations is shown in Figure 1 [21]. Figure 2 shows the displacement amplitudes caused by S1 and S2 tides under the CM frame for each station with respect to the latitude. Here, we should keep in mind that the amplitude and phase obtained has an uncertainty since the RP03 model itself definitely has an uncertainty in it. However, at the current stage, we could not assess the uncertainty as ECMWF only provides surface pressure data, and there is no precision information along with it. Correspondingly, the RP03 model, which comes from the operational ECMWF analysis, also provides no error information. Therefore, we could not evaluate or validate the precision of the calculated S1/S2 amplitude and phase. This is a disadvantage of using the RP03 model to model atmospheric tides. Nevertheless, it was the only choice until now and has been suggested as the convention [6]. In the future, with the development of environmental loading service (e.g., there exists a new surface pressure data with precision information, or a more operational atmospheric tidal model comes into existence), we may find an effective way to evaluate the uncertainty of the S1/S2 induced surface displacement.   From Figure 2, we observe that both the S1 and S2 tides could result in big vertical annual movements, and the predominant effect is latitude dependent, among which the amplitudes of stations near the equator are larger than 1 mm. This is easy to understand, as both the diurnal and semidiurnal surface pressure variation are decreased with the increasing latitude. For the horizontal components, the impacts are much smaller, in particular for the S2 tide, whose global annual amplitudes are less than 0.2 mm for both the North (N) and East (E) components. The magnitude of the S1 tide is a little bigger compared with the S2 tide. Specifically, the annual amplitude for the East (E) component of S1 is within the range of 0.5-0.7 mm globally without the latitudinal dependence, while the amplitude for the North (N) component becomes larger with the increasing latitude, and the biggest variation is concentrated above the latitude of ±60° at around 0.6 mm. This is because the S1 and S2 tides here are induced by surface heating, and the horizontal thermally driven force is predominantly diurnal. Additionally, there is a strong North-South temperature gradient existing in the daily heating system, which caused the latitudinal dependence of S1 amplitude for the North (N) component.

GPS Data Reprocessing
In the previous section, we confirm that the S1/S2 tides could cause significant annual movement of the global IGS stations under the CM frame. To investigate the impact of S1/S2 on the position of global IGS stations, a contrast experiment was designed using GAMIT 10.4 by reprocessing the GPS data of the above 109 IGS stations. In experiment A, we do not consider the S1/S2 effects, while in experiment B, we calculate the S1/S2 induced surface displacement [7,8]. In this way, the contribution of S1/S2 to daily GPS solution and continuous GPS time series of the global IGS stations could be determined by differencing between experiment A and B. The global grid of the amplitudes and phases for the S1 and S2 tides is obtained from the above-mentioned RP03 model (ftp://everest.mit.edu/pub/GRIDS/atl.grid). Only GPS data on Wednesdays during the period from January 1999 to December 2003 are processed. In the beginning, we tried to process the daily GPS data for one month and found that station movement during one week was very small, which could almost be neglected. Additionally, for the ITRF establishment, although the IGS analysis center did the daily GPS data processing, these daily time series were first stacked into weekly solutions and then combined with other technologies to obtain the ITRF [6]. Due to these reasons, we finally decided to choose only one day (Wednesdays here) from one week to represent  Figure 2. The annual amplitudes of the surface displacements caused by the S1 and S2 tides with respect to the latitude. The red and blue circle represent the S1 and S2 amplitudes, respectively. From Figure 2, we observe that both the S1 and S2 tides could result in big vertical annual movements, and the predominant effect is latitude dependent, among which the amplitudes of stations near the equator are larger than 1 mm. This is easy to understand, as both the diurnal and semidiurnal surface pressure variation are decreased with the increasing latitude. For the horizontal components, the impacts are much smaller, in particular for the S2 tide, whose global annual amplitudes are less than 0.2 mm for both the North (N) and East (E) components. The magnitude of the S1 tide is a little bigger compared with the S2 tide. Specifically, the annual amplitude for the East (E) component of S1 is within the range of 0.5-0.7 mm globally without the latitudinal dependence, while the amplitude for the North (N) component becomes larger with the increasing latitude, and the biggest variation is concentrated above the latitude of ±60 • at around 0.6 mm. This is because the S1 and S2 tides here are induced by surface heating, and the horizontal thermally driven force is predominantly diurnal. Additionally, there is a strong North-South temperature gradient existing in the daily heating system, which caused the latitudinal dependence of S1 amplitude for the North (N) component.

GPS Data Reprocessing
In the previous section, we confirm that the S1/S2 tides could cause significant annual movement of the global IGS stations under the CM frame. To investigate the impact of S1/S2 on the position of global IGS stations, a contrast experiment was designed using GAMIT 10.4 by reprocessing the GPS data of the above 109 IGS stations. In experiment A, we do not consider the S1/S2 effects, while in experiment B, we calculate the S1/S2 induced surface displacement [7,8]. In this way, the contribution of S1/S2 to daily GPS solution and continuous GPS time series of the global IGS stations could be determined by differencing between experiment A and B. The global grid of the amplitudes and phases for the S1 and S2 tides is obtained from the above-mentioned RP03 model (ftp://everest.mit.edu/pub/GRIDS/atl.grid). Only GPS data on Wednesdays during the period from January 1999 to December 2003 are processed. In the beginning, we tried to process the daily GPS data for one month and found that station movement during one week was very small, which could almost be neglected. Additionally, for the ITRF establishment, although the IGS analysis center did the daily GPS data processing, these daily time series were first stacked into weekly solutions and then combined with other technologies to obtain the ITRF [6]. Due to these reasons, we finally decided to choose only one day (Wednesdays here) from one week to represent the station's weekly movement. This could be treated as an approximate weekly solution, not rigorous, but reasonable from a practical point of view.
Besides, the RINEX GPS observation files are on a daily basis, e.g., one file per day, and the obtained daily time series is indeed a 24-h average value. Therefore, the "continuous" in GNSS Geodesy usually refers to continuous observations, not the solutions, as both the daily or weekly solutions are a kind of average. As long as the GPS observation data during the selected time span is continuous, and the time span is long enough, either solution (daily, weekly, monthly, etc.) could be used to study the characteristics of the stations' movement. Their difference only comes from the data selection of different sampling rates. Here, we calculate an approximate weekly solution. Our whole data span is five consecutive years, and the shortest time period is about 1.6 years (Supplementary File, Table S1). As the behavior of the S1/S2 mainly shows an annual and semi-annual pattern, this continuous 5-year time span should be enough to describe the seasonal characteristics caused by S1/S2.
During the data reprocessing, satellite orbits, earth orientation parameters (EOPs), site coordinates together with the tropospheric delay and the horizontal gradient parameters are resolved simultaneously. Loose constraints are implemented on the stations, among which the constraints of the IGS core stations are set as 5 cm, and non-core stations are 1 dm. Where possible, ambiguities are fixed to integers [22]. The satellite cut-off elevation angle is chosen to be 10 degrees and the site-specific, elevation-dependent weighting of the observations based on an assessment of the post phase residuals was applied [7,8]. Corrections of solid earth tide, ocean tide, as well as pole tide have been implemented [23][24][25]. Please note that at the S2 frequency, the ocean responds to both gravitational potential and air pressure (called the ocean radiational tides). In the FES2004 ocean tide model, the ocean's response to these atmospheric tides is already modelled separately through the site displacement [6]. The Vienna Mapping Function 1 (VMF1) troposphere mapping functions are used to calculate the tropospheric delay [26,27]. Where possible, receiver-independent exchange (RINEX) meteorological files (.m) are used to provide pressure and temperature for the a priori zenith hydrostatic delay; otherwise, values from VMF1 are used [28]. The absolute antenna phase center offsets and variations are used (igs08_1636.atx) [29]. Impacts of the second and third-order ionospheric delay are considered, during which the International Geomagnetic Reference Field 11 (IGRF 11) is selected to calculate the second-order ionospheric delay [22]. Non-tidal loading corrections are not implemented according to the IERS Conventions 2010 [6]. Finally, the gross error rejection and datum transformation are implemented to the daily baseline resolutions using GLOBK to obtain the stations' coordinate time series and velocities under the ITRF08 [30]. During the datum transformation, only 6 parameters, that is 3 rotation and 3 translation parameters, are estimated to reduce the aliasing effects of the unmodeled surface mass loadings, e.g., the non-tidal atmospheric loading [31].
For the data selection, we agree that using old data from 1999 to 2003 to do the experiment is a disadvantage since IGS data during that period was with low accuracy orbits and models. However, because we only focus on the difference caused by S1/S2, this disadvantage would not impact the relative coordinate difference much. Moreover, from our experience, a more recent and longer time series would not impact much on the final result, as we did a comparative analysis on the minor ocean tide using both a short time (1999-2003) and a longer period (1998-2010) [32], and the conclusions were almost the same. Therefore, we think that the experiment here is reasonable and effective. In the near future, we will use the most recent GPS data (for example, 2011-2017) to repeat this experiment and see whether higher-precision orbit and models would affect the S1/S2 induced surface displacement or not.
Another important thing that should be addressed here is the tropospheric delay correction. In fact, during our GPS data processing, we estimate the Zenith Tropospheric Delay (ZTD) every hour together with all the other parameters. However, since the topic of this paper only focuses on the GPS coordinate time series, we did not perform the analysis related to the ZTD here. The tropospheric delay is highly related to the GPS height and it is indeed possible to absorb part of the signal in the height component. Nevertheless, until now there has been no effective way to precisely separate these two items on a large scale. Because of its high correlation, the impact of S1/S2 should also propagate into ZTD. From this point of view, our analysis here is not complete. In our next step, we will make a more detailed analysis of the impact of S1/S2 on both the GPS derived coordinate and ZTD using longer and more recent GPS observations.

Results and Discussion
In this section, we will discuss the impact of S1/S2 on the GPS solutions from four aspects, including its influence on the 24-h GPS daily position difference and global long-term velocity, together with the weighted root mean square (WRMS) and the spectral characteristics of the coordinate time series.

The 24-h GPS Daily Position Difference Caused by S1/S2 Tides
Through reprocessing the GPS data using our contrast experiment, we obtain the S1/S2 effects on the 24-h daily GPS solutions. Figure 3 shows an example of the GPS time series for station BRAZ (Brasilia, Brazil) and REYK (Reykjavik, Iceland) due to S1/S2 effects. We observe that the maximum daily position difference caused by S1/S2 for the low-latitude station BRAZ reaches more than 6 mm for the horizontal components, and up to 3 mm for the vertical component ( Figure 3, top panels). With respect to the high-latitude station REYK, the magnitude is much smaller, with a peak-to-peak variation of up to 3 mm and 2 mm for the vertical and horizontal components, respectively (see Figure 3, bottom panels). This is reasonable, although the maximum position difference is bigger than the maximum annual amplitude of S1 and S2 since both the S1 and S2 would interact with all the other periodic signals (e.g., tidal signal), and sometimes the total signal becomes extremely strong when all the signals exhibit resonance behavior.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 15 separate these two items on a large scale. Because of its high correlation, the impact of S1/S2 should also propagate into ZTD. From this point of view, our analysis here is not complete. In our next step, we will make a more detailed analysis of the impact of S1/S2 on both the GPS derived coordinate and ZTD using longer and more recent GPS observations.

Results and Discussion
In this section, we will discuss the impact of S1/S2 on the GPS solutions from four aspects, including its influence on the 24-h GPS daily position difference and global long-term velocity, together with the weighted root mean square (WRMS) and the spectral characteristics of the coordinate time series.

The 24-h GPS Daily Position Difference Caused by S1/S2 Tides
Through reprocessing the GPS data using our contrast experiment, we obtain the S1/S2 effects on the 24-h daily GPS solutions. Figure 3 shows an example of the GPS time series for station BRAZ (Brasilia, Brazil) and REYK (Reykjavik, Iceland) due to S1/S2 effects. We observe that the maximum daily position difference caused by S1/S2 for the low-latitude station BRAZ reaches more than 6 mm for the horizontal components, and up to 3 mm for the vertical component ( Figure 3, top panels). With respect to the high-latitude station REYK, the magnitude is much smaller, with a peak-to-peak variation of up to 3 mm and 2 mm for the vertical and horizontal components, respectively (see Figure 3, bottom panels). This is reasonable, although the maximum position difference is bigger than the maximum annual amplitude of S1 and S2 since both the S1 and S2 would interact with all the other periodic signals (e.g., tidal signal), and sometimes the total signal becomes extremely strong when all the signals exhibit resonance behavior.   Figure 4 gives the Standard Deviation (STD) of the S1/S2 induced GPS time series for each station with respect to the latitude. Note that the average position difference of each station is 0 in our case, that is, the bias value is 0; therefore, we only show the STD value here to describe the magnitude of S1/S2. This is also expected as the site displacement caused by S1/S2 is calculated from the mean annual amplitude and phase of S1/S2, while our data spans the whole year around. We observe that the S1/S2 loading could cause typically sub-mm variations in the 24-h GPS daily position difference in general. As expected, the vertical component exhibits the biggest impact. Most stations (80%) have an induced STD of more than 0.5 mm, including all the stations within the latitude of ±40 degree, and the magnitude increases with the decreasing latitude for stations located at the northern hemisphere ( Figure 4, bottom panel), among which the low-latitude station GUAM (Guam, West Pacific) shows a scatter of more than 1.4 mm. The impacts of the S1/S2 correction on the horizontal components are slightly smaller, and the magnitude also shows a strong latitude dependent pattern for the northern hemisphere stations, which is quite similar to the vertical component, with a maximum scatter of up to 1 mm and 0.7 mm for the East and North components, respectively. Compared with the signal itself (Figure 2), the vertical annual amplitude of both of S1 and S2 could reach more than 1.2 at a low latitude. If we add the S1 and S2 displacement together, the maximum displacement could reach more than 2.4 mm. This is consistent with Figures 3 and 4. The maximum annual amplitude of S1 and S2 for the horizontal components reach more than 0.6 mm and 0.1 mm, respectively, while for most stations, the STD of the S1/S2 induced horizontal displacement is below 0.7 mm. This also shows a consistency between the horizontal signal and displacement.
To be more complete, Table S1 in the supplement lists the average uncertainty of the GPS coordinate time series for each station together with the time span. We can see that the mean error of the horizontal and vertical displacement caused by S1/S2 are 3 mm and 9 mm respectively, while the average time span for the selected stations is 3.6 years. This information could prove that our result is valid within the noise level, and the S1/S2 induced surface displacement is significant compared with the noise level.  Figure 3. The GPS coordinate time series for station BRAZ (top) and REYK (bottom) due to the S1/S2 effects. Figure 4 gives the Standard Deviation (STD) of the S1/S2 induced GPS time series for each station with respect to the latitude. Note that the average position difference of each station is 0 in our case, that is, the bias value is 0; therefore, we only show the STD value here to describe the magnitude of S1/S2. This is also expected as the site displacement caused by S1/S2 is calculated from the mean annual amplitude and phase of S1/S2, while our data spans the whole year around. We observe that the S1/S2 loading could cause typically sub-mm variations in the 24-h GPS daily position difference in general. As expected, the vertical component exhibits the biggest impact. Most stations (80%) have an induced STD of more than 0.5 mm, including all the stations within the latitude of ±40 degree, and the magnitude increases with the decreasing latitude for stations located at the northern hemisphere ( Figure 4, bottom panel), among which the low-latitude station GUAM (Guam, West Pacific) shows a scatter of more than 1.4 mm. The impacts of the S1/S2 correction on the horizontal components are slightly smaller, and the magnitude also shows a strong latitude dependent pattern for the northern hemisphere stations, which is quite similar to the vertical component, with a maximum scatter of up to 1 mm and 0.7 mm for the East and North components, respectively. Compared with the signal itself (Figure 2), the vertical annual amplitude of both of S1 and S2 could reach more than 1.2 at a low latitude. If we add the S1 and S2 displacement together, the maximum displacement could reach more than 2.4 mm. This is consistent with Figures 3 and 4. The maximum annual amplitude of S1 and S2 for the horizontal components reach more than 0.6 mm and 0.1 mm, respectively, while for most stations, the STD of the S1/S2 induced horizontal displacement is below 0.7 mm. This also shows a consistency between the horizontal signal and displacement.
To be more complete, Table S1 in the supplement lists the average uncertainty of the GPS coordinate time series for each station together with the time span. We can see that the mean error of the horizontal and vertical displacement caused by S1/S2 are 3 mm and 9 mm respectively, while the average time span for the selected stations is 3.6 years. This information could prove that our result is valid within the noise level, and the S1/S2 induced surface displacement is significant compared with the noise level.

The S1/S2 Effects on the Long-Term Velocity of Global IGS Stations
In the above section, we confirm the contribution of S1/S2 to the 24-h GPS daily solutions. Now, we focus on its impact on the long-term GPS solutions, as the 24-h mean of S1/S2 caused displacement is non-zero at most of the time (Figure 3) and would interact with the other tidal signals, which may cause a linear trend. Figure 5 shows the stations' velocity variation rate with respect to the latitude. Here, we only estimate the velocity of the IGS stations with a continuous observation span of more than 1.5 years, that is, 78 GPS weeks, so as to ensure the reliability of the statistics. Finally, the effective station number remains 104. The basic equation can be written as where neu V Δ denotes the velocity variation rate for each component, original neu v represents the station's long-term velocity without S1/S2 correction, neu v Δ is the velocity difference before and after considering the S1/S2 effects. Negative means the velocity decrease after implementing the S1/S2 corrections. Figure 6 illustrates the spatial distribution of the stations' vertical velocity variation rate caused by the S1/S2 effects. From Figures 5 and 6, we observe that the S1/S2 effects on the global long-term velocity have no latitude dependency, which is quite different from that of the 24-h GPS solution (Figure 4). It could only cause a very small horizontal velocity change, with 98% and 95% of the selected stations' velocity variations smaller than 1% for the N and E components, respectively ( Figure 5, the top and middle panels). With respect to the very small vertical velocity itself, the vertical velocity change caused by S1/S2 is much larger ( Figure 5, bottom panel, Figure 6). A total of 65% of the stations' variation is larger than 1% and those stations with a vertical velocity change larger than 5% account for 25%, including many coastal/island stations, among which the maximum variation rate reaches more than 40%, for example, station TCMS (Taiwan, China) and SUTH (Sutherland, South Africa).

The S1/S2 Effects on the Long-Term Velocity of Global IGS Stations
In the above section, we confirm the contribution of S1/S2 to the 24-h GPS daily solutions. Now, we focus on its impact on the long-term GPS solutions, as the 24-h mean of S1/S2 caused displacement is non-zero at most of the time (Figure 3) and would interact with the other tidal signals, which may cause a linear trend. Figure 5 shows the stations' velocity variation rate with respect to the latitude. Here, we only estimate the velocity of the IGS stations with a continuous observation span of more than 1.5 years, that is, 78 GPS weeks, so as to ensure the reliability of the statistics. Finally, the effective station number remains 104. The basic equation can be written as where ∆V neu denotes the velocity variation rate for each component, v neu origina l represents the station's long-term velocity without S1/S2 correction, ∆v neu is the velocity difference before and after considering the S1/S2 effects. Negative means the velocity decrease after implementing the S1/S2 corrections. Figure 6 illustrates the spatial distribution of the stations' vertical velocity variation rate caused by the S1/S2 effects. From Figures 5 and 6, we observe that the S1/S2 effects on the global long-term velocity have no latitude dependency, which is quite different from that of the 24-h GPS solution (Figure 4). It could only cause a very small horizontal velocity change, with 98% and 95% of the selected stations' velocity variations smaller than 1% for the N and E components, respectively ( Figure 5, the top and middle panels). With respect to the very small vertical velocity itself, the vertical velocity change caused by S1/S2 is much larger ( Figure 5, bottom panel, Figure 6). A total of 65% of the stations' variation is larger than 1% and those stations with a vertical velocity change larger than 5% account for 25%, including many coastal/island stations, among which the maximum variation rate reaches more than 40%, for example, station TCMS (Taiwan, China) and SUTH (Sutherland, South Africa).  . The spatial distribution of the S1/S2 induced long-term vertical velocity variation rate. The unit of the variation rate is in percentage (%). The black and white dots in the figure indicate that the velocity variation rate for the station is smaller than the minimum and larger than the maximum value on the scale.

Weighted Root Mean Square (WRMS) Analysis
To investigate the impact of S1/S2 on the non-linear characteristics of global IGS stations, we calculate the WRMS of the coordinate time series for the 109 stations before and after considering the S1/S2 effects (Supplementary, Table S2). Here, we also only estimate the WRMS of IGS stations with a continuous observation span of more than 1.5 years. Figure 7 illustrates the spatial distribution of the S1/S2 induced WRMS variation rate. Negative means that the WRMS reduced after implementing the S1/S2 correction. For the definition of the WRMS and its variation rate, please refer to Reference [33].   . The spatial distribution of the S1/S2 induced long-term vertical velocity variation rate. The unit of the variation rate is in percentage (%). The black and white dots in the figure indicate that the velocity variation rate for the station is smaller than the minimum and larger than the maximum value on the scale.

Weighted Root Mean Square (WRMS) Analysis
To investigate the impact of S1/S2 on the non-linear characteristics of global IGS stations, we calculate the WRMS of the coordinate time series for the 109 stations before and after considering the S1/S2 effects (Supplementary, Table S2). Here, we also only estimate the WRMS of IGS stations with a continuous observation span of more than 1.5 years. Figure 7 illustrates the spatial distribution of the S1/S2 induced WRMS variation rate. Negative means that the WRMS reduced after implementing the S1/S2 correction. For the definition of the WRMS and its variation rate, please refer to Reference [33].  Figure 6. The spatial distribution of the S1/S2 induced long-term vertical velocity variation rate. The unit of the variation rate is in percentage (%). The black and white dots in the figure indicate that the velocity variation rate for the station is smaller than the minimum and larger than the maximum value on the scale.

Weighted Root Mean Square (WRMS) Analysis
To investigate the impact of S1/S2 on the non-linear characteristics of global IGS stations, we calculate the WRMS of the coordinate time series for the 109 stations before and after considering the S1/S2 effects (Supplementary, Table S2). Here, we also only estimate the WRMS of IGS stations with a continuous observation span of more than 1.5 years. Figure 7 illustrates the spatial distribution of the S1/S2 induced WRMS variation rate. Negative means that the WRMS reduced after implementing the S1/S2 correction. For the definition of the WRMS and its variation rate, please refer to Reference [33]. The spatial distribution of the S1/S2 induced WRMS variation rate. The black and white dots in the figure indicate that the WRMS variation rate for the station is smaller than the minimum and larger than the maximum value on the scale. From top to bottom are the U, E, and N components. The unit of the variation rate is in percentage (%).
From Figure 7, we can see that S1/S2 could reduce the vertical WRMS of 65% of the stations, concentrating mainly in the continental regions except for South America and South Africa, but the magnitude is quite small. Stations with a WRMS reduction larger than 1% only account for 13%, among which most of them are coastal stations, e.g., OHIG (Antarctica), whose WRMS reduction reaches 4%. For the horizontal components, stations with a reduced WRMS account for less than 40%. The maximum reduction in the E component is located in most coastal and island stations, such as station STJO (St. John's, NL, Canada), RBAY (Richardsbay, South Africa), with the WRMS reduction reaching more than 3%. With respect to the N component, stations with a WRMS reduction rate larger than 1% account for 13%, but there is no obvious spatial pattern. Even today there is still a big discrepancy between the global IGS stations' non-linear position time series and the environmental loading induced displacement, in particular, for the horizontal components [9,[33][34][35][36][37][38]. The S1/S2 modeling at the observation level may indicate a possible source that would narrow Figure 7. The spatial distribution of the S1/S2 induced WRMS variation rate. The black and white dots in the figure indicate that the WRMS variation rate for the station is smaller than the minimum and larger than the maximum value on the scale. From top to bottom are the U, E, and N components. The unit of the variation rate is in percentage (%).
From Figure 7, we can see that S1/S2 could reduce the vertical WRMS of 65% of the stations, concentrating mainly in the continental regions except for South America and South Africa, but the magnitude is quite small. Stations with a WRMS reduction larger than 1% only account for 13%, among which most of them are coastal stations, e.g., OHIG (Antarctica), whose WRMS reduction reaches 4%. For the horizontal components, stations with a reduced WRMS account for less than 40%. The maximum reduction in the E component is located in most coastal and island stations, such as station STJO (St. John's, NL, Canada), RBAY (Richardsbay, South Africa), with the WRMS reduction reaching more than 3%. With respect to the N component, stations with a WRMS reduction rate larger than 1% account for 13%, but there is no obvious spatial pattern. Even today there is still a big discrepancy between the global IGS stations' non-linear position time series and the environmental loading induced displacement, in particular, for the horizontal components [9,[33][34][35][36][37][38]. The S1/S2 modeling at the observation level may indicate a possible source that would narrow the gap. Further work still needs to be done on the contribution of S1/S2 to the inconsistency between the GPS coordinate time series and the surface displacement driven by environmental loading.

Spectral Analysis
To investigate the impact of S1/S2 on the periodic motion of global IGS reference stations, we calculate the Power Spectrum Density (PSD) of the 109 stations in both Experiment A and B using the CATS software [39]. Since there is no big difference in the PSD characteristics between coastal/island and inland stations, or stations in different latitudinal regions, here, we only show the global stacked PSD results. Figure 8 shows the global filtered PSD stacking results for each component before and after S1/S2 correction (top panels), together with their differences (bottom panels). In the bottom panels of Figure 8, the positive value means that the PSD amplitude is decreased after correcting the S1/S2 effects. For the implementation of the PSD stacking and filtering, please refer to Reference [35] for details.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 15 the gap. Further work still needs to be done on the contribution of S1/S2 to the inconsistency between the GPS coordinate time series and the surface displacement driven by environmental loading.

Spectral Analysis
To investigate the impact of S1/S2 on the periodic motion of global IGS reference stations, we calculate the Power Spectrum Density (PSD) of the 109 stations in both Experiment A and B using the CATS software [39]. Since there is no big difference in the PSD characteristics between coastal/island and inland stations, or stations in different latitudinal regions, here, we only show the global stacked PSD results. Figure 8 shows the global filtered PSD stacking results for each component before and after S1/S2 correction (top panels), together with their differences (bottom panels). In the bottom panels of Figure 8, the positive value means that the PSD amplitude is decreased after correcting the S1/S2 effects. For the implementation of the PSD stacking and filtering, please refer to Reference [35] for details.  From Figure 8, we observe that no matter whether S1/S2 correction is applied or not, the annual signal and the anomalous harmonics of 1.04 cpy ± 0.008 cpy (known as the GPS draconitic year, [35]) until the 4th draconitic signals are the dominant periodic signals in both the vertical and horizontal components of the global IGS stations, and the magnitude of the vertical component is about 10-times larger than the horizontal components ( Figure 8, top panels). We find that the annual and the semi-annual amplitude of the U component, together with the semi-annual amplitude of the E component, decreased to some extent after considering the S1/S2 effects, which further extends the previous conclusion on the contribution of S1/S2 in the IGS stations inside China [9]. It could also help reduce the gap between the seasonal variations detected from the GNSS sites, loading models and GRACE data [40]. However, the annual and the semi-annual amplitudes of the N component exhibit a significant increase. We also find that the amplitude of the 4th draconitic signal in the U component and the 2nd draconitic signal in the N component decreased a lot, but the other anomalous harmonics exhibit an amplitude increase. This result partly agrees with Tregoning and Watson's findings that the neglect of S1/S2 would introduce anomalous propagated signals with periods that closely match the GPS draconitic semiannual period (~175.7 days, [7,8]); thus, we conclude that S1/S2 may also be a possible source that caused the anomalous harmonics of 1.04 cpy. This result may also help interpret part of the draconitic errors found by Reference [12] that results in the large perturbations in residuals over mainland China. However, our stacking result does not show an amplitude decrease of the GPS draconitic annual signal. This may due to the differences in the IGS station selection, time period, and the applied GPS processing strategies in these two papers. For example, we use solutions between 1999-2003, while they have solutions from 2000-2008. We choose more stations within Europe, but they choose much more southern hemisphere stations. We considered the effects of higher-order ionospheric delay, absolute antenna phase center offsets and variations correction, but it was not clear whether they were considered or not from References [7,8].
Another interesting result from Figure 8 is that most long-term periodic signals with a frequency equal to or smaller than 1 cpy in the U component also decreased to some degree. This phenomenon confirms the recent finding that the unmodeled diurnal and sub-diurnal periodic signals could indeed propagate and produce spurious long-term signals in the GPS coordinate time series [41].

Conclusions
The effects of S1/S2 are often neglected in routine precise long-baseline GPS data processing. In this paper, we first recall the S1/S2 modeling method and then give the magnitude and spatial distribution of the global IGS station's annual amplitude caused by S1/S2 itself under the CM frame. We find that both S1 and S2 tides could result in big vertical annual movements, in particular, for those near the equator, whose annual amplitudes are larger than 1 mm. The S1 tide could also cause in-negligible horizontal annual movement, with a magnitude within the range of 0.5-0.7 mm globally for the East component, and around 0.6 mm for the North component above the latitude of ±60 • .
Secondly, we reprocess the GPS data of 109 globally distributed IGS stations based on the GAMIT software through two runs, with and without S1/S2 correction. We find that the S1/S2 loading causes typically sub-mm variations in the 24-h GPS daily displacement in general, and the vertical component exhibits the biggest impact. A total of 80% of the stations have a vertical scatter of more than 0.5 mm. The magnitude of the scatter increases with the decreasing latitude for northern hemisphere stations, among which the maximum STD reaches up to 1.5 mm, 1 mm and 0.7 mm for the Up, East and North components respectively at low-latitude stations, e.g., station GUAM.
We then study the impact of S1/S2 on the long-term velocity of global IGS stations. We find that 65% of the stations' vertical velocity change caused by S1/S2 is larger than 1%; 25% of the selected stations, including many coastal/island stations' variation, are larger than 5%; and the maximum variation rate reaches more than 40% (e.g., station TCMS, SUTH). Thus, we suggest that the S1/S2 effects should be better corrected during high-precision GPS data processing for coastal and ocean regions, so as to obtain a more reliable vertical motion for the applications of GPS in geodynamic studies, e.g., sea level change.
Moreover, the non-linear characteristics of global IGS stations caused by S1/S2 are discussed. We find that, in general, S1/S2 could only introduce small WRMS variations in the global IGS stations, even for the coastal/island stations, where the maximum variation rate is smaller than 4%. 65%, 39%, and 38% of the stations' WRMS reduced for the U, E and N components, respectively, after considering the S1/S2 effects, among which the maximum reduction has a magnitude of 1-4% for the U and E components gathered mostly along coasts or on the island (e.g., station OHIG, STJO). Further work still needs to be done on whether the S1/S2 modeling at the observation level of GPS data processing would narrow the gap between the global IGS stations' non-linear position and the environmental loading induced displacement time series.
At last, the spectral characteristics of the global IGS stations caused by S1/S2 is analyzed. Our results show that the S1/S2 correction could partly reduce the annual and semi-annual amplitudes of the U component together with the semi-annual amplitude of the E component, but the annual and the semi-annual amplitudes of the N component exhibit a significant increase. Furthermore, most long-term periodic signals with a frequency equal to or smaller than 1 cpy in the U component also decreased to some extent after correcting the S1/S2 effects. We also find that adding the S1/S2 correction could significantly reduce the power of the 4th draconitic signal of the anomalous harmonics of 1.04 cpy ± 0.008 cpy in the U component and the 2nd draconitic signal in the N component, but the other anomalous harmonics exhibit an amplitude increase, which partly confirms Tregoning and Watson's finding that S1/S2 may also be a possible source that caused the anomalous harmonics of 1.04 cpy.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/10/7/1125/ s1. Figure S1: Annual amplitudes and phases of the S1/S2 induced surface displacement, Table S1: Average uncertainty of GPS time series for each station and its time span, Table S2: WRMS for the GPS time series before (left) and after S1/S2 correction (right).