A Novel Method to Improve the Estimation of Ocean Tide Loading Displacements for K 1 and K 2 Components with GPS Observations

: The accurate estimation of ocean tide loading displacements is essential and necessary for geodesy, oceanic and geophysical studies. It is common knowledge that K 1 and K 2 tidal constituents estimated from Global Positioning System (GPS) observations are unsatisfactory because their tidal periods are nearly same to the revisit cycle or orbital period of GPS constellation. To date, this troublesome problem is not fully solved. In this paper, we revisit this important issue and develop a novel method based on the unique characteristic of tidal waves to separate GPS-system errors from astronomical K 1 /K 2 tides. The well-known credo of smoothness indicates that tidal admittances of astronomical constituents in a narrow band can be expressed as smooth functions of tidal frequencies, while the interference of GPS-system errors seriously damages the smooth nature of observed tidal admittances. Via quadratic ﬁtting, smooth functions of tidal frequencies for tidal admittances can be determined, thus, astronomical K 1 and K 2 tides can be interpolated using ﬁtted quadratic functions. Three GPS stations are selected to demonstrate our method because of their typicality in terms of poor estimates of K 1 /K 2 tidal parameters related to GPS-system errors. After removing GPS-systematical contributions based on our method, corrected K 1 /K 2 tides at three GPS stations are much closer to the modeled K 1 /K 2 tides from FES2014, which is one of the most accurate tide models. Furthermore, the proposed method can be easily applied to other areas to correct GPS-system errors because their smooth nature is valid for global tidal signals.


Introduction
Ocean tides are periodical sea water movements (mainly diurnal and semi-diurnal) driven by the astronomical forcing from the Moon and the Sun [1][2][3]. As a result of ocean mass redistribution caused by ocean tides, the periodical deformation of the earth's surface is called ocean tide loading (OTL) displacement [4]. The spatial features of OTL are mainly decided by nearby ocean tides and the amplitudes of OTL can reach as high as 10 cm in some coastal regions [5,6]. Precise estimation of OTL displacement for major tidal constituents is fundamental and vital for multidisciplinary fields such as the global navigation satellite system (GNSS), as well as oceanic, geophysical and seismological studies [7][8][9].
As a powerful and widely used tool to observe OTL displacement, GPS is superior to other geodesic tools (e.g., very long baseline interferometry) in terms of high accuracy, wide coverage and low cost [5]. However, GPS-observed OTL displacements have some defects. It is well known that the revisit cycle of the GPS constellation is nearly same as K 1 tidal period (23.934 h), while the orbital period of the GPS constellation is nearly the same as K 2 tidal period (11.967 h) [6]. Consequently, K 1 and K 2 tidal amplitudes/phases estimated from GPS observations are significantly poorer than other major tidal constituents [10]. 2 of 16 To solve this knotty problem, numerous efforts have been made. Abbaszadeh et al. (2020) [11] try to use the combination of GLONASS and GPS observations because the orbital period of GLONASS constellation is 11.25 h, which is different from K 2 tidal period (11.967 h). However, it is found that the estimations for K 1 and K 2 tides are still unsatisfactory although there are some improvements [11]. Wei et al. (2022b) [6] go a step further by using the combination of multi-GNSS including GPS, GLONASS, Galileo and BeiDou. Every GNSS has unique orbital parameters, for example, the orbital period of Galileo satellites is~14 h [6]. Thus, it is obvious that tidal estimations from multi-GNSS observations are less influenced by the systematic errors of GPS compared to GPS-only observations. Although multi-GNSS observations improve the estimation of OTL displacement, the method proposed by Wei et al. (2022b) [6] does not completely solve the problem because GPS-system errors still exist as long as GPS observations are used.
The purpose of this paper is to revisit this troublesome but vital issue and to propose a novel yet simple method based on the unique ability of tidal signals to fully address the problem. Our paper is structured as follows. Data and methodology are shown in Section 2. Section 3 shows the results as well as discussions, followed by the summary in Section 4.

Methodology and Data
Our method is established on the basis of the credo of smoothness proposed by Munk and Cartwright (1966) [12]. The credo of smoothness indicates that tidal admittances are smooth functions of tidal frequency in a narrow band [13,14]. Such smooth functions can be simply expressed as quadratic functions [12][13][14], namely, Equation (1): (1) A (w) represents tidal admittances while w is tidal frequency. Model parameters a 1 , a 2 , a 3 can be determined via least squares. Tidal admittances for amplitudes (i.e., normalized amplitudes) are defined by the ratios of observed tidal amplitudes to theoretical equilibrium tidal amplitudes, while tidal admittances for phases are defined by the phase lags between observed tidal phases and theoretical equilibrium tidal phases. Note that equilibrium tidal amplitudes can be calculated based on the equilibrium tide theory via the s_equilibrium_tide function in the S_TIDE toolbox developed by the authors [15][16][17], while phase lags can be directly estimated from observations using tidal harmonic analysis. Tidal admittances represent actual tidal responses of ocean environment to the astronomical forcing, which mainly depend on the frequencies of tidal waves, water depth and coastlines. For tides in specific sea areas, tidal admittances only depend on tidal frequencies because the ocean environment is the same for all kinds of tidal waves. Semi-diurnal tides have similar tidal periods (close to 0.5 days), which means that their wave lengths are also similar. As a result, tidal responses for semi-diurnal tides are highly consistent; namely, their co-tidal charts are closely similar in terms of spatial patterns (see Figure A1 in the Appendix A). Therefore, tidal admittances for semi-diurnal tides can be approximated by the smooth functions of tidal frequencies. Diurnal tides also have similar features. Figure A2 in Appendix A displays tidal admittances for diurnal ocean tides at Cuxhaven (Germany), located at the mouth of the Elbe river, as a typical example. Normalized amplitudes are quadratic functions of diurnal frequencies while phase lags are nearly linear functions of diurnal frequencies over these frequency ranges. As a result of the interference of non-astronomical factors (i.e., river flow), observed tidal admittances at Cuxhaven slightly deviate from the fitted curves, but their smoothness is generally robust and reliable. Since OTL originated from ocean tides, it also inherits the nature of smooth tidal admittances. Consequently, co-tidal charts of OTL for tides with similar periods also have similar spatial patterns (see in Hart-Davis et al. (2021) [9] for examples). Figure 1 displays diurnal tidal admittances for the up component based on the FES2014 OTL model [18] at Aira (31.8241 • N, 130.5996 • E) (location shown in Figure A3 in Appendix A), as a typical example. Both normalized amplitudes and phase lags of diurnal OTL are quadratic functions of diurnal frequencies at Aira. The existence of GPS-system errors severely destroys the smoothness of tidal admittances but also provides opportunities to remove non-astronomical contributions to K 1 and K 2 tidal estimations (see Section 3 for details).
reliable. Since OTL originated from ocean tides, it also inherits the nature of smooth tidal admittances. Consequently, co-tidal charts of OTL for tides with similar periods also have similar spatial patterns (see in Hart-Davis et al. (2021) [9] for examples). Figure 1 displays diurnal tidal admittances for the up component based on the FES2014 OTL model [18] at Aira (31.8241°N, 130.5996°E) (location shown in Figure A3 in Appendix A), as a typical example. Both normalized amplitudes and phase lags of diurnal OTL are quadratic functions of diurnal frequencies at Aira. The existence of GPS-system errors severely destroys the smoothness of tidal admittances but also provides opportunities to remove non-astronomical contributions to K1 and K2 tidal estimations (see Section 3 for details). To demonstrate our method, we use OTL displacement estimates at three GPS stations (Aira, Bill and Helj) in the CM frame provided by Yuan et al. (2013) [5]. Tidal amplitudes (denoted by H) and phase lags (denoted by G) of eight major constituents are obtained by classical harmonic analysis (Equation (2). Z(t) are GPS-observed OTL displacements while Z0 is a constant. fn, un, and vn are nodal factor, nodal angle and astronomical argument of n-th tidal constituent, which are all known: To demonstrate our method, we use OTL displacement estimates at three GPS stations (Aira, Bill and Helj) in the CM frame provided by Yuan et al. (2013) [5]. Tidal amplitudes (denoted by H) and phase lags (denoted by G) of eight major constituents are obtained by classical harmonic analysis (Equation (2). Z(t) are GPS-observed OTL displacements while Z 0 is a constant. f n , u n , and v n are nodal factor, nodal angle and astronomical argument of n-th tidal constituent, which are all known: Aira is discussed in detail due to its typicality in terms of poor estimates of K 1 /K 2 tidal parameters related to GPS-system errors. Semi-diurnal tides are generally larger than diurnal tides at Aira. As shown in Table 1, M 2 phase lag (144.0 • ) for the north component is close to N 2 phase lag (134.3 • ) and their phase difference is only 9.7 • . The frequency difference between S 2 and K 2 tides is much smaller than that between M 2 and N 2 , therefore, the phase difference between S 2 and K 2 phase lags should be smaller than 9.7 • . However, the GPS-observed S 2 and K 2 phase difference is as high as 85.4 • , which clearly indicates the interference of GPS-system errors. Similarly, the phase difference between K 1 and P 1 tides should be much smaller than that between Q 1 and O 1 tides. However, unfortunately, the observed K 1 -P 1 phase difference (20.5 • ) for the north component is noticeably larger than the observed Q 1 -O 1 phase difference (4.3 • ). Table 1. Amplitudes (mm) and Greenwich phase lags ( • ) at the up (U), north (N), east (E) components of eight main tidal constituents at Aira GPS station. Equilibrium (Eq) tidal amplitudes (mm) are also shown.

Results
Figure 2 displays semi-diurnal tidal admittances for the up component at Aira. Both normalized amplitudes and phase lags are quadratic functions of semi-diurnal tidal frequencies. Note that fitted quadratic curves are determined via ordinary least squares using N 2 , M 2 and S 2 tidal admittances. K 2 tidal admittances are not used because they are contaminated by GPS-system errors. As shown in Figure 2, K 2 tidal admittances significantly deviate from fitted quadratic curves. According to the tidal theory, S 2 and K 2 tides have highly close periods, thus, their admittances should also be close as long as they have astronomical origins. However, the interference of the systematic errors of GPS observations induces abnormal changes in K 2 admittances. Via the interpolation, astronomical K 2 admittances can be calculated (red dots in Figure 2), which are close to the observed S 2 admittances. Since equilibrium K 2 tidal amplitude is known, astronomical K 2 amplitude and phase lag at the up component at Aira can be obtained as 2.38 mm and 149.4 • , respectively. Similar to semi-diurnal tidal admittances, diurnal tidal admittances for the up component at Aira are also quadratic functions of diurnal frequencies ( Figure 3). K 1 admittances unnaturally deviate from the fitted curves, especially for phase lags, due to GPS-system errors. Astronomical K 1 admittances can be estimated using the interpolation (red dots in Figure 3), which are close to observed P 1 admittances. Based on known equilibrium K 1 tidal parameters, astronomical K 1 amplitude and phase lag at the up component can be obtained as 10.17 mm and 243.43 • , respectively.
The observed K 1 tides by GPS have two main energy sources: one is the astronomical K 1 tide driven by astronomical forcing and the other is the non-astronomical K 1 tide mainly induced by GPS-system errors. Hence, as displayed in Figure 4a, the observed K 1 tide (black arrow) is the vectorial composition of the astronomical K 1 tide (red arrow) and non-astronomical K 1 tide (blue arrow). The non-astronomical K 1 amplitude and phase lag are 2.64 mm and 139.0 • , respectively. Although non-astronomical and astronomical K 1 tides have same tidal frequency, their features (namely, amplitudes and phase lags) are distinct because they are originated from different processes. Similar to semi-diurnal tidal admittances, diurnal tidal admittances for the up component at Aira are also quadratic functions of diurnal frequencies ( Figure 3). K1 admittances unnaturally deviate from the fitted curves, especially for phase lags, due to GPSsystem errors. Astronomical K1 admittances can be estimated using the interpolation (red dots in Figure 3), which are close to observed P1 admittances. Based on known equilibrium K1 tidal parameters, astronomical K1 amplitude and phase lag at the up component can be obtained as 10.17 mm and 243.43°, respectively.  The observed K1 tides by GPS have two main energy sources: one is the astronomical K1 tide driven by astronomical forcing and the other is the non-astronomical K1 tide mainly induced by GPS-system errors. Hence, as displayed in Figure 4a, the observed K1 tide (black arrow) is the vectorial composition of the astronomical K1 tide (red arrow) and non-astronomical K1 tide (blue arrow). The non-astronomical K1 amplitude and phase lag are 2.64 mm and 139.0°, respectively. Although non-astronomical and astronomical K1 tides have same tidal frequency, their features (namely, amplitudes and phase lags) are distinct because they are originated from different processes.  The modeled K1 amplitude and phase lag for the up component at Aira are 9.11 mm and 242.3°, respectively, based on FES2014 tide model [18]. Figure 4a clearly indicates that astronomical K1 tide obtained from our method is highly consistent with modeled K1 tide (green arrow). Note that FES2014 tide model is one of the most accurate tide models [9], thus, it is used as a reference. The vector difference (VD) between modeled K1 tidal vector and uncorrected K1 tidal vector (namely, GPS-observed K1 tide) is 2.41 mm. After correction using our method (namely, removing non-astronomical contributions), VD sharply decreases to 1.08 mm (Figure 4b), which indicates the effectiveness and reliability of our method. Similarly, corrected K2 tide is also much closer to modeled K2 tide compared to uncorrected K2 tide (Figure 4b). Figures 5 and 6 display semi-diurnal and diurnal tidal admittances for the north component at Aira. Similar to the up component, K1 and K2 tidal admittances at the north component anomalously deviate from the fitting curves, especially for phase lags. The interpolated K1 and K2 tidal admittances (red dots) are generally close to observed P1 and S2 tidal admittances, which is consistent with tidal theory. Results for the east component are depicted by Figures 7 and 8 and, in addition, indicate the effectiveness of our method. Note that nearly all global tide models only provide OTL for the up component, thus, we do not compare corrected results at the north and east components with model results due to lack of model data. The vectorial composition of the astronomical K 1 tide (red arrow) and non-astronomical K 1 tide (blue arrow) furnishes the observed K 1 tide (black arrow) at Aira. Green arrow represents the modeled K 1 tidal vector (using FES2014 tidal model). (b) The vector differences (VDs) between modeled K 1 /K 2 tidal vectors and uncorrected/corrected K 1 /K 2 tidal vectors (corrections are performed on the basis of the credo of smoothness).
The modeled K 1 amplitude and phase lag for the up component at Aira are 9.11 mm and 242.3 • , respectively, based on FES2014 tide model [18]. Figure 4a clearly indicates that astronomical K 1 tide obtained from our method is highly consistent with modeled K 1 tide (green arrow). Note that FES2014 tide model is one of the most accurate tide models [9], thus, it is used as a reference. The vector difference (VD) between modeled K 1 tidal vector and uncorrected K 1 tidal vector (namely, GPS-observed K 1 tide) is 2.41 mm. After correction using our method (namely, removing non-astronomical contributions), VD sharply decreases to 1.08 mm (Figure 4b), which indicates the effectiveness and reliability of our method. Similarly, corrected K 2 tide is also much closer to modeled K 2 tide compared to uncorrected K 2 tide (Figure 4b).
Figures 5 and 6 display semi-diurnal and diurnal tidal admittances for the north component at Aira. Similar to the up component, K 1 and K 2 tidal admittances at the north component anomalously deviate from the fitting curves, especially for phase lags. The interpolated K 1 and K 2 tidal admittances (red dots) are generally close to observed P 1 and S 2 tidal admittances, which is consistent with tidal theory. Results for the east component are depicted by Figures 7 and 8 and, in addition, indicate the effectiveness of our method. Note that nearly all global tide models only provide OTL for the up component, thus, we do not compare corrected results at the north and east components with model results due to lack of model data.         To further validate our method, we tested it using two other typical GPS stations called Bill (33.5782 • N, 242.9354 • E) and Helj (54.1745 • N, 7.8931 • E). As shown in Figures 9 and 10, corrected K 1 /K 2 tides at Helj and Bill are much closer to modeled K 1 /K 2 tides compared to uncorrected K 1 /K 2 tides. The details on smooth tidal admittances at Bill and Helj can be found in the Appendix A (Figures A4-A7) Figures A4-A7) to avoid repetitions. In summary, the results of three GPS stations all show the validity of proposed method.  The vectorial composition of the astronomical K 2 tide (red arrow) and non-astronomical K 2 tide (blue arrow) furnishes the observed K 2 tide (black arrow) at Helj. Green arrow represents the modeled K 2 tidal vector (using FES2014 tidal model). (b) The vector differences (VDs) between modeled K 1 /K 2 tidal vectors and uncorrected/corrected K 1 /K 2 tidal vectors (corrections are performed on the basis of the credo of smoothness).

Discussions: The Advantages and Limitations of the Proposed Methodology
The proposed method has many advantages. Its principle is simple and easy to understand. The realization of this method is uncomplicated, without numerous computation resources. Although this method can be applied to the whole globe, in theory, cautions should be taken in some coastal resonant areas where significant nonlinear interactions between main tidal constituents exist. Taking K2 tide, for example, with the exception of astronomical forcing, K2 tide can also be generated by the nonlinear interaction between M2, O1 and K1 tides (i.e., M2 − O1 + K1 or K1 + K1). Furthermore, S2 seasonality can also contribute to K2 tide. In most sea areas, ocean tides are semi-diurnal, thus, astronomical contributions are dominant in the generation of K2 tide. However, for some sea areas where diurnal tidal resonance exists (such as the Gulf of Tonkin in the South China Sea [19]), local diurnal tides are abnormally larger than semi-diurnal tides, therefore, the nonlinear origin of K2 tides is important and not negligible. The proposed method not only removes the contribution of GPS-system errors to K2 tides, but it also eliminates the contributions of nonlinear interactions and S2 seasonality. Namely, the obtained K2 tide via our method is mainly astronomical. Figure 10. (a) The vectorial composition of the astronomical K 1 tide (red arrow) and non-astronomical K 1 tide (blue arrow) furnishes the observed K 1 tide (black arrow) at Bill. Green arrow represents the modeled K 1 tidal vector (using FES2014 tidal model). (b) The vector differences (VDs) between modeled K 1 /K 2 tidal vectors and uncorrected/corrected K 1 /K 2 tidal vectors (corrections are performed on the basis of the credo of smoothness).

Discussions: The Advantages and Limitations of the Proposed Methodology
The proposed method has many advantages. Its principle is simple and easy to understand. The realization of this method is uncomplicated, without numerous computation resources. Although this method can be applied to the whole globe, in theory, cautions should be taken in some coastal resonant areas where significant nonlinear interactions between main tidal constituents exist. Taking K 2 tide, for example, with the exception of astronomical forcing, K 2 tide can also be generated by the nonlinear interaction between M 2 , O 1 and K 1 tides (i.e., M 2 − O 1 + K 1 or K 1 + K 1 ). Furthermore, S 2 seasonality can also contribute to K 2 tide. In most sea areas, ocean tides are semi-diurnal, thus, astronomical contributions are dominant in the generation of K 2 tide. However, for some sea areas where diurnal tidal resonance exists (such as the Gulf of Tonkin in the South China Sea [19]), local diurnal tides are abnormally larger than semi-diurnal tides, therefore, the nonlinear origin of K 2 tides is important and not negligible. The proposed method not only removes the contribution of GPS-system errors to K 2 tides, but it also eliminates the contributions of nonlinear interactions and S 2 seasonality. Namely, the obtained K 2 tide via our method is mainly astronomical.
In the construction of smooth functions of tidal admittances, we only use eight main tidal constituents but do not use minor tidal constituents such as 2Q 1 , J 1 , 2N 2 , NU 2 , LAD 2 and L 2 . These minor diurnal and semi-diurnal tides are much smaller than major constituents in terms of equilibrium tidal amplitudes, which means that they are more easily influenced by background noises and nonlinear interactions between main tidal constituents. The frequency of 2N 2 tide is the same as 2NM 2 tide, which can be generated by nonlinear interactions between N 2 and M 2 tides (i.e., N 2 + N 2 − M 2 ). The frequency of NU 2 tide is the same as MLS 2 tide, which can be generated by nonlinear interactions between M 2 , L 2 and S 2 tides (i.e., M 2 + L 2 − S 2 ). The frequency of LAD 2 tide is the same as SNM 2 tide, which can be generated by nonlinear interactions between S 2 , N 2 and M 2 tides (i.e., S 2 + N 2 − M 2 ). The frequency of L 2 tide is the same as 2MN 2 tide, which can be generated by nonlinear interactions between M 2 and N 2 tides (i.e., M 2 + M 2 − N 2 ). In some coastal areas, the contributions of nonlinear interactions to these minor constituents may be comparable or even larger than direct astronomical contributions. Therefore, including these minor constituents may significantly distort the smoothness of tidal admittances.
In fact, major constituents can also obtain energy from nonlinear interactions. For example, M 2 tidal frequency is the same as KO 2 tide, which is generated from the nonlinear interaction between K 1 and O 1 tides (i.e., K 1 + O 1 ). M 2 tide can also obtain energy from the nonlinear interaction of N 2 , S 2 and LAD 2 tides (i.e., N 2 + S 2 − LAD 2 ). S 2 tidal frequency is the same as KP 2 tide, which is derived from the nonlinear interaction between K 1 and P 1 tides (i.e., K 1 + P 1 ). Furthermore, S 2 tide can obtain energy from the nonlinear interaction of M 2 and MU 2 tides (i.e., M 2 + M 2 − MU 2 ). However, compared to minor constituents, the amplitudes of major constituents contributed to by astronomical forcing are large enough to resist the interference of nonlinear contributions in most sea areas. In some resonant sea areas, like the Gulf of Tonkin, the smoothness of semi-diurnal tidal admittances may be distorted to some extent. In fact, with the exception of nonlinear interactions, radiational S 2 tides driven by solar radiation also contribute to observed S 2 tides. Observed S 2 admittances are altered in phases by adding 5.9 • and in amplitudes by multiplying by 0.97 due to radiational S 2 tides [20]. Overall, radiational S 2 tides only slightly distort S 2 tidal admittances.

Conclusions
Originating from ocean tides, OTL displacements have complex and irregular spatial variations [21]. Accurate knowledge of OTL displacement is beneficial for numerous geophysical and geodetic studies [22,23]. Although GPS is widely used for observing OTL displacements due to its various advantages, it also has some disadvantages which can severely influence the accuracy of K 1 and K 2 tidal estimations. Tidal constants of major constituents are estimated from GPS observations via classical harmonic analysis, which is a frequency-dependent method. Therefore, tidal harmonic analysis can only obtain tidal parameters of given tidal frequencies and cannot distinguish astronomical contributions from GPS-system errors.
To fully solve this question, we propose a novel but straightforward method based on the distinct physical properties of non-astronomical and astronomical tides. For astronomical tides, their admittances in a narrow band are similar and can be expressed as smooth functions of tidal frequencies [24]. Via simple interpolation, astronomical K 1 and K 2 tides can be separated from GPS-system errors. Three representative GPS stations have demonstrated the validity of our method. Moreover, as a universal method, it can be easily applied to other areas to correct GPS-system errors because the credo of smoothness of tidal admittances is effective for all kinds of tidal signals. However, care should be taken in some coastal resonant sea areas where the contribution of nonlinear interactions may make up a considerable part of observed main tidal constituents. It should be noted that the proposed method always tends to eliminate all non-astronomical contributions. Furthermore, in addition to our method, artificial intelligence algorithms [25,26] may be helpful in solving this problem, which deserves further exploration.
Funding: This study is jointly supported by National Natural Science Foundation of China (42206022, 42076024, 41821004), China Postdoctoral Science Foundation (2022M713677) and the Qingdao postdoctoral application research project (QDBSH202108).

Conflicts of Interest:
The authors declare no conflict of interest. Appendix A Figure A1. Co-tidal charts for K 1 , O 1 , M 2 and S 2 tides in the Gulf of Tonkin. Tidal constants are obtained from EOT20 model, which is developed based on FES2014 [9]. Co-tidal charts are drawn using the s_draw_tidalchart function in the S_TIDE toolbox. Co-tidal charts for M 2 and S 2 tides are similar, while those for K 1 and O 1 tides are highly consistent.