1. 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.
2. 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 (
,
,
,
), 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
where
and
represent the S1, S2 induced 3-D displacements for the Up (U), North (N) and East (E) components,
represents the total displacements of the S1 and S2 tides,
,
,
,
are the in-phase and out-of-phase amplitudes of the S1 and S2 tides in unit of millimeter (mm),
is in fractions of a universal time (UT1) day,
and
denote the frequencies of the S1 and S2 tides respectively, with
and
.
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.
3. 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.
5. 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.