Accurate Evaluation of Vertical Tidal Displacement Determined by GPS Kinematic Precise Point Positioning: A Case Study of Hong Kong

Global Positioning System (GPS) kinematic precise point positioning (KPPP) is an effective approach for estimating the Earth’s tidal deformation. The accuracy of KPPP is usually evaluated by comparing results with tidal models. However, because of the uncertainties of the tidal models, the accuracy of KPPP-estimated tidal displacement is difficult to accurately determine. In this paper, systematic vector differences between GPS estimates and tidal models were estimated by least squares methods in complex domain to analyze the uncertainties of tidal models and determine the accuracy of KPPP-estimated tidal displacements. Through the use of GPS data for 12 GPS reference stations in Hong Kong from 2008 to 2017, vertical ocean tide loading displacements (after removing the body tide effect) for eight semidiurnal and diurnal tidal constituents were obtained by GPS KPPP. By an in-depth analysis of the systematic and residual differences between the GPS estimates and nine tidal models, we demonstrate that the uncertainty of the tidal displacement determined by GPS KPPP for the M2, N2, O1, and Q1 tidal constituents is 0.2 mm, and for the S2 constituent it is 0.5 mm, while the accuracy of the GPS-estimated K1, P1, and K2 tidal constituents is weak because these three tidal constituents are affected by significant common-mode errors. These results suggest that GPS KPPP can be used to precisely constrain the Earth’s vertical tidal displacement in the M2, N2, O1, and Q1 tidal frequencies.


Introduction
The tidal force of the Sun and Moon can cause deformation of the solid Earth, which is named "body tide", and it also causes periodic rise and fall of the ocean mass, which is called "ocean tide". The response of the solid Earth to this mass distribution is known as "ocean tide loading" (OTL). Modern high-precision geodetic technology can determine large-scale deformation at the sub-millimeter level, but the effect of tidal displacement must be considered [1][2][3][4]. Body tide can now be modeled with an accuracy of usually less than 1% [5]. OTL displacement can be obtained by convolution integral of an ocean tide model and the Green's function calculated from an Earth model [6], which we call an OTL model. Penna et al. [7] reported that the height error of M2 OTL displacement can reach around 20% (~8 mm), depending on the ocean tide model adopted and the handling strategy for the grid cells in coastal areas. The change of the vertical OTL displacement in some areas when using various Earth models can reach 2-3 mm [8]. Therefore, OTL is still the main uncertain factor in tidal displacement analysis.
Since 2000, the Global Positioning System (GPS) has been widely used to estimate OTL displacement, including eight semidiurnal and diurnal tidal constituents. The methods for estimating

Method of KPPP-estimated OTL Displacement
The influence of OTL in the horizontal direction is about one-third of that in the vertical direction [8]. Therefore, in this paper, we mainly analyzed the OTL displacement in the vertical direction. According to the International Earth Rotation and Reference Systems Service (IERS) Convention 2010 [26], the vertical OTL displacement of a station can be represented by four semidiurnal (M2, S2, N2 and K2) and four diurnal (K1, O1, P1 and Q1) tidal constituents, depending on the frequencies: ∆c = 8 j=1 f j A cj cos(ω j t + χ j (t 0 ) + u j − φ cj ) (1) where A cj and φ cj respectively represent the amplitude and phase of the j-th tidal constituent at the station; f j and u j denote the nodal modulation corrections; ω j is the angular frequency; and χ j (t 0 ) is the astronomical argument at reference time t 0 . Table 1 lists the periods of the eight semidiurnal and diurnal tidal constituents. For a given location, except that the amplitude and phase of the j-th tidal constituent need to determined, the remaining parameters can be obtained according to the expression of these arguments in the IERS Convention 2010 [26]. GPS kinematic precise point positioning (KPPP) is an absolute positioning method, which can directly determine the position of a receiver based on satellite constellations [27]. Therefore, an absolute displacement series in vertical direction with a high sampling rate can be obtained by the GPS KPPP solution. OTL displacements can usually be effectively estimated by using a 4-years of GPS observation [25], but simultaneous nodal correction is also needed. The OTL amplitude and phase of the semidiurnal and diurnal tidal constituents can be estimated from the KPPP coordinate series by using the classical harmonic analysis method [28].

GPS Data Acquisition and Processing
Hong Kong is located on the north coast of the South China Sea, and is mainly composed of the Kowloon Peninsula and many islands. The Hong Kong coastline is complex and is affected by OTL. GPS observation data for 2008-2017 provided by the Surveying and Mapping Office of the Hong Kong Lands Department (https://www.geodetic.gov.hk/smo/gsi/programs/tc/index.htm) were used, in which the HKFN station was replaced by the T430 station after 2014. The distribution of the GPS reference stations in the Hong Kong area are shown in Figure 1. semidiurnal and diurnal tidal constituents. For a given location, except that the amplitude and phase of the j-th tidal constituent need to determined, the remaining parameters can be obtained according to the expression of these arguments in the IERS Convention 2010 [26]. GPS kinematic precise point positioning (KPPP) is an absolute positioning method, which can directly determine the position of a receiver based on satellite constellations [27]. Therefore, an absolute displacement series in vertical direction with a high sampling rate can be obtained by the GPS KPPP solution. OTL displacements can usually be effectively estimated by using a 4-years of GPS observation [25], but simultaneous nodal correction is also needed. The OTL amplitude and phase of the semidiurnal and diurnal tidal constituents can be estimated from the KPPP coordinate series by using the classical harmonic analysis method [28].

GPS Data Acquisition and Processing
Hong Kong is located on the north coast of the South China Sea, and is mainly composed of the Kowloon Peninsula and many islands. The Hong Kong coastline is complex and is affected by OTL. GPS observation data for 2008-2017 provided by the Surveying and Mapping Office of the Hong Kong Lands Department (https://www.geodetic.gov.hk/smo/gsi/programs/tc/index.htm) were used, in which the HKFN station was replaced by the T430 station after 2014. The distribution of the GPS reference stations in the Hong Kong area are shown in Figure 1.  By using 24-h GPS observation data in RINEX format, the kinematic PPP solution for every 600 s was carried out with Bernese GNSS software version 5.2 [29] (http://www.bernese.unibe.ch/ docs/DOCU52.pdf), in which the elevation cutoff angle was set to 3 • . The main parameter files, i.e., satellite orbit, 30 s satellite clock, Earth rotation, ionosphere parameter, and differential code bias (DCB), were provided by the European Center for Orbit Determination in Europe (COD). The tropospheric delay was modeled using the Global Mapping Function (GMF) [30]. The antenna phase center file used PCV_COD.I08. Body tide and pole tide were modeled according to the IERS Convention 2010 [26]. No corrections were applied to the ocean and atmospheric tide loading.
The estimated KPPP coordinates were in the IGS08 reference frame. The origin of IGS08 is defined as the center of mass (CM), including the solid Earth, the ocean, and the atmosphere, which is usually called the CM frame. Coordinates under IGS08 reference frame were transformed into the station coordinate system, including east, north, and vertical directions. The linear trend of the vertical series was removed. A threshold of 200 mm was set to eliminate gross errors. The % of Gross Errors are shown in Table 2. The eight semidiurnal and diurnal tidal components in the vertical direction at 12 stations were estimated using the method described in Section 2.1. The amplitudes and phases of OTL displacements for eight tidal constituents determined by GPS KPPP are listed in Tables A1 and A2.

OTL Model and Tidal Signal Analysis
OTL displacement predicted by a model is usually based on the Green's function [6], which is determined by the adopted Earth model and represents the response characteristics of the solid Earth to mass loading. The OTL displacement u for a location r is as follows: where Z is a vector expression of the tide at r provided by the ocean tide model; Ω represents the global ocean area; ρ is the density of sea water; and G is the Green's function, characterizing the load deformation caused by the unit mass of sea water. In this paper, eight recent ocean tide models were considered: FES2014b, GOT4.10c, TPXO.8, DTU10, EOT11a, HAMTIDE, NAO99b, and OSU12 [31,32]. The OTL displacements predicted by the FES2014b and GOT4.10c models based on the STW105 Earth model were provided by the International Mass Loading Service [33], for which the density of sea water was set to 1025 kg/m 3 . The predictions of the OTL displacements by FES2014b, TPXO.8, DTU10, EOT11a, HAMTIDE, NAO99b, and OSU12 were provided by the Ocean Tide Loading Provider [34], for which the density of sea water was set to 1030 kg/m 3 . The FES2014b and TPXO.8 tide models based on the Preliminary Reference Earth Model (PREM) were calculated with CARGA [7]. The DTU10, EOT11a, HAMTIDE, NAO99b, and OSU12 tide models based on the Gutenberg-Bullen A Earth model were calculated with OLFG/OLMPP [35].
After modeling the body tide in the GPS data processing, the GPS-estimated tidal deformation is mainly made up of the OTL displacement. A tidal constituent can be conveniently expressed in terms of the in-phase and out-of-phase amplitudes of the phasors on the complex plane [13]. The phasor vector relationship between the GPS estimate and OTL model is shown in Figure 2. with OLFG/OLMPP [35]. A total of nine OTL models were obtained, and we refer to them as GOT4.10c + STW, FES2014b + STW, FES2014b + PREM, TPXO.8 + PREM, DTU10 + GB, EOT11a + GB, HAMTIDE + GB, NAO99b + GB, and OSU12 + GB, respectively. After modeling the body tide in the GPS data processing, the GPS-estimated tidal deformation is mainly made up of the OTL displacement. A tidal constituent can be conveniently expressed in terms of the in-phase and out-of-phase amplitudes of the phasors on the complex plane [13]. The phasor vector relationship between the GPS estimate and OTL model is shown in Figure 2. The phasor difference between the GPS estimate and OTL model for the j-th tidal constituent at station k is as follows: is the tidal displacement estimated by GPS, and , , denotes the tidal deformation predicted by the OTL model. The complex form of , can be written as: where i is the imaginary unit; and , and Φ , represent the residual amplitude and phase of the j-th tidal constituent, respectively.
The RMS of the phasor difference for the j-th tidal constituent between the GPS estimate and OTL model for N stations is calculated as follows: In a small region, the phasor difference for the j-th tidal constituent between the GPS estimate and OTL model at station k can be divided into two kinds of components. The first kind of component is related to the specific station, i.e., the multipath effect on GPS observation. The second kind of component is applicable to all stations, i.e., the unmodeled body tide, the uncertainties of the OTL model caused by ocean tide model errors and mantle anelasticity, the GPS orbit errors, and the signal propagation delay. Therefore, the phasor difference , of the j-th tidal constituent at station k can The phasor difference between the GPS estimate and OTL model for the j-th tidal constituent at station k is as follows: where Z GPS,k,j is the tidal displacement estimated by GPS, and Z OTL,k,j denotes the tidal deformation predicted by the OTL model. The complex form of Z k,j can be written as: where i is the imaginary unit; and A k,j and Φ k,j represent the residual amplitude and phase of the j-th tidal constituent, respectively. The RMS of the phasor difference for the j-th tidal constituent between the GPS estimate and OTL model for N stations is calculated as follows: In a small region, the phasor difference for the j-th tidal constituent between the GPS estimate and OTL model at station k can be divided into two kinds of components. The first kind of component is related to the specific station, i.e., the multipath effect on GPS observation. The second kind of component is applicable to all stations, i.e., the unmodeled body tide, the uncertainties of the OTL model caused by ocean tide model errors and mantle anelasticity, the GPS orbit errors, and the signal propagation delay. Therefore, the phasor difference Z k,j of the j-th tidal constituent at station k can be decomposed into Z residual,k,j , which is related to the specific station, and Z system,j , which is applicable to all stations. For the N stations, the phasor difference of the j-th tidal constituent can be expressed as follows: Next, we explain how to obtain Z residual,k,j and Z system,j on the complex plane by the least-squares method. Equation (6) is written in matrix form: Equation (7) can be abbreviated as: where I denotes the unit matrix, and Z system,j is the parameter to be estimated. It is assumed that Z residual,k,j (k = 1, . . . , N) obey a Gaussian distribution and are uncorrelated with each other. The expected value and variance of Z residual,k,j (k = 1, . . . , N) are E |Z residual,k,j | = 0 and D |Z residual,k,j | = σ 2 , respectively. Therefore, the variance-covariance matrix of Z j can be written as: where D() is the variance-covariance matrix. If we let V donate the correction value of Z residual,j and Z system, j represent the estimated value of Z system,j , then formula (8) can be transformed as: From formula (9), we know that Z k,j (k = 1, . . . , N) are equal-weight observations. Therefore, to obtainẐ system,j based on the least-squares principle, we only need to satisfy V T V = min, in which V T is the transpose matrix of V. In order to obtain the minimum value of V T V, the derivative of V T V tô Z system, j should be 0. The detailed process is as follows: According to Equations (10) and (11), the formula after eliminating V can be obtained as: Therefore, the least-squares estimate of complex Z system,j is: Finally, Z residual,k,j (k = 1, . . . , N) can be obtained by formula (6).

Comparison with the OTL Models
To analyze the KPPP-estimated tidal constituent stability, according to Equations (6) and (13), we calculated the residual and systematic phasor differences between the GPS estimates (using the observed data for years 4 to 10 since 2008) and the GOT4.10c + STW OTL model, for each of the eight tidal constituents at 12 stations. The corresponding RMSs of the residual and systematic phasor difference of the eight tidal constituents were obtained by formula (5), as shown in Figure 3.

Comparison with the OTL Models
To analyze the KPPP-estimated tidal constituent stability, according to formulas (6) and (13), we calculated the residual and systematic phasor differences between the GPS estimates (using the observed data for years 4 to 10 since 2008) and the GOT4.10c + STW OTL model, for each of the eight tidal constituents at 12 stations. The corresponding RMSs of the residual and systematic phasor difference of the eight tidal constituents were obtained by formula (5), as shown in Figure 3.  Figure 3a,b indicate that, for the M2, N2, O1 and Q1 tidal constituents, after the fourth year, the RMS of the residual and systematic phasors for the four tidal constituents is stabilized. For the S2 tidal constituent, its residual RMS is stabilized after the fourth year; however, it is not until the ninth year that the systematic RMS of the S2 constituent begins to stabilize. Although the residual and systematic RMS of the P1 constituent is stabilized after the fourth year, the magnitude of the residual and systematic RMS for the P1 tidal constituent is significant. The residual RMS of the K1 and K2 constituents is stabilized after the fifth year, but their systematic RMS does not convergence until the tenth year. Figure 4a shows the spatial distributions of the phasor vector differences, between the GPS estimates (using the observed data for 2008-2017) and the GOT4.10c + STW OTL model, of the eight tidal constituents at 12 stations. Figure 4b shows the systematic phasor vector differences estimated from Figure 4a by formula (13) for the eight tidal constituents. The residual vector phasor differences were obtained by differencing Figure 4a,b, as shown in Figure 4c. The expression of the phasor vector difference for a given tidal constituent is given in Figure 2.  Figure 3a,b indicate that, for the M2, N2, O1 and Q1 tidal constituents, after the fourth year, the RMS of the residual and systematic phasors for the four tidal constituents is stabilized. For the S2 tidal constituent, its residual RMS is stabilized after the fourth year; however, it is not until the ninth year that the systematic RMS of the S2 constituent begins to stabilize. Although the residual and systematic RMS of the P1 constituent is stabilized after the fourth year, the magnitude of the residual and systematic RMS for the P1 tidal constituent is significant. The residual RMS of the K1 and K2 constituents is stabilized after the fifth year, but their systematic RMS does not convergence until the tenth year. Figure 4a shows the spatial distributions of the phasor vector differences, between the GPS estimates (using the observed data for 2008-2017) and the GOT4.10c + STW OTL model, of the eight tidal constituents at 12 stations. Figure 4b shows the systematic phasor vector differences estimated from Figure 4a by formula (13) for the eight tidal constituents. The residual vector phasor differences were obtained by differencing Figure 4a,b, as shown in Figure 4c. The expression of the phasor vector difference for a given tidal constituent is given in Figure 2.  From Figure 4a, we can see that the spatial distributions of the phasor vector differences for the M2, K2, K1, O1, and P1 tidal constituents have obvious coherence. The systematic phasor differences could be well derived, as shown in Figure 4b. Figure 4a also shows that the spatial distributions of the phasor vector differences for the S2, N2, and Q1 tidal constituents had no obvious coherence; therefore, there were no distinguishing systematic differences for these three constituents, as shown From Figure 4a, we can see that the spatial distributions of the phasor vector differences for the M2, K2, K1, O1, and P1 tidal constituents have obvious coherence. The systematic phasor differences could be well derived, as shown in Figure 4b. Figure 4a also shows that the spatial distributions of the phasor vector differences for the S2, N2, and Q1 tidal constituents had no obvious coherence; therefore, there were no distinguishing systematic differences for these three constituents, as shown in Figure 4b. Figure 4c indicates that the residual phasor differences of the eight tidal constituents had no spatial correlation, i.e., they are related to the location.

Accuracy Assessment for the KPPP-Estimated OTL
To validate the accuracy of the GPS KPPP estimates of the eight tidal constituents, according to Equations (5), (6) and (13), we calculated the RMS of the total phasor differences, the systematic phasor differences and the corresponding residual phasor differences between the GPS estimates (using the observed data for 2008-2017) and the nine OTL models, for each of the eight tidal constituents at 12 stations. The RMS values of the eight tidal constituents are shown in Figure 5.  Figure 5 also shows that there was obvious systematic RMS of the K1, P1, and K2 tidal constituents between the GPS estimates and the nine OTL models, of around 6 mm, 2 mm, and 9 mm, respectively. This indicates that the K1, K2, and P1 tidal constituents are affected by the significant common-mode errors for all the stations. These errors may include: (1) the unmodeled body tide; (2) the OTL model uncertainty caused by the ocean tide model and mantle inelasticity; and (3) GPSrelated error, such as satellite orbit error, which has almost the same effect for all stations. Firstly, for the three tidal constituents of K1, K2, and P1, the unmodeled body tide error should be at the submillimeter level [12]. Secondly, the OTL displacement differences for the K1, P1, and K2 tidal constituents among the predictions by the nine different OTL models are at the sub-millimeter level, which indicates that that the OTL models are unlikely to have caused such misfits. Therefore, the common-mode errors are likely caused by the GPS-related errors. After removing the systematic differences of the K1, P1, and K2 tidal constituents, the residual RMS was around 0.8 mm, 0.5 mm, and 1 mm, respectively. The residual RMS was still significant, which indicates that these three tidal  Figure 5 shows that the systematic RMS of the N2 and Q1 tidal constituents between the GPS estimates and the different OTL models varies from 0.02 to 0.40 mm and from 0.12 to 0.40 mm, respectively. The systematic RMS varies significantly with the different OTL models, which shows that the systematic of the N2 and Q1 tidal constituents was mainly caused by the uncertainties of the OTL models. The uncertainties of the OTL models result in the systematic differences. After removing the systematic differences, the RMS of the residual phasors of the N2 and Q1 tidal constituents are all around 0.2 mm. The residual RMS does not change with the different OTL models, so the residual RMS reflects the uncertainty of the GPS estimates.
As shown in Figure 5, the systematic RMS change of the M2 and O1 tidal constituents with the various OTL models is from 0.01 to 1.02 mm and 0.20 to 1.01 mm, respectively. These results indicate that the systematic RMS of the M2 and O1 constituents was mainly caused by the uncertainties of the OTL models. After removing the systematic differences, the residual RMS of the O1 tidal constituent was around 0.2 mm, and no longer changed with the various OTL models. The residual RMS values of the M2 tidal constituent between the GPS estimates and the six OTL models of GOT4.10c + STW, FES2014b + STW, FES2014b + PREM, TPXO.8 + PREM, NAO99b + GB, and OSU12 + GB were all 0.2 mm. However, the residual RMS of the M2 tidal constituent between the GPS estimates and the other three OTL models (DTU10 + GB, EOT11a + GB, HAMTIDE + GB) were around 0.35 mm. This indicates that the residual RMS of the M2 tidal constituent changes slightly with the different OTL models. Therefore, the residual RMS should absorb the errors of the three OTL models, which may not exist entirely as systematic differences, because the M2 tidal resolution of the DTU10, EOT11a, and HAMTIDE ocean tide models is 1/8 • × 1/8 • smaller than the 1/2 • × 1/2 • size of Hong Kong. Overall the uncertainty of the GPS-estimated M2 tidal constituent is at the 0.2 mm level. Figure 5 also shows that there was obvious systematic RMS of the K1, P1, and K2 tidal constituents between the GPS estimates and the nine OTL models, of around 6 mm, 2 mm, and 9 mm, respectively. This indicates that the K1, K2, and P1 tidal constituents are affected by the significant common-mode errors for all the stations. These errors may include: (1) the unmodeled body tide; (2) the OTL model uncertainty caused by the ocean tide model and mantle inelasticity; and (3) GPS-related error, such as satellite orbit error, which has almost the same effect for all stations. Firstly, for the three tidal constituents of K1, K2, and P1, the unmodeled body tide error should be at the sub-millimeter level [12]. Secondly, the OTL displacement differences for the K1, P1, and K2 tidal constituents among the predictions by the nine different OTL models are at the sub-millimeter level, which indicates that that the OTL models are unlikely to have caused such misfits. Therefore, the common-mode errors are likely caused by the GPS-related errors. After removing the systematic differences of the K1, P1, and K2 tidal constituents, the residual RMS was around 0.8 mm, 0.5 mm, and 1 mm, respectively. The residual RMS was still significant, which indicates that these three tidal constituents are probably affected by the station-related error. Previous studies [13,14] inferred that K1 and K2 may be affected by the satellite orbit error and multipath effect, because the period of the K1 coincides with that of the satellite constellation, and the period of the K2 is almost the same as the satellite orbital period, while P1 may be disturbed by residual GPS signal propagation errors or non-tidal daily cycles in the atmosphere because its period is very close to solar day.
As shown in Figure 5, the total RMS of the S2 tidal constituent between the GPS estimates and the nine OTL models were between 0.5 to 0.6 mm. There are some not obvious systematic differences of 0.1 to 0.3 mm between the GPS estimates and the nine different OTL models, which may arise from the effect of the atmospheric S2 tide. The effect of the atmospheric S2 tide was analyzed using the GEOSFPIT atmospheric model provided by the International Mass Loading Service [8], with a resolution of 0.625 • × 0.5 • . The results show that the amplitude of the atmospheric tide displacement, as predicted by the GEOSFPIT model, is only 0.22 mm at Hong Kong. Figure 6 shows that, after correcting the atmospheric tide by the GEOSFPIT model, the RMS residual of the S2 tidal constituent between the GPS estimates and some OTL models slightly reduced. Overall there was still an obvious RMS residual of about 0.5 mm between the GPS estimates and the OTL models. Figure 5 shows that the residual RMS of the S2 tidal constituent were all nearly 0.5 mm. This indicates that there are some station-related errors of around 0.5 mm, which arise from GPS-related error. GPS signal propagation error or non-tidal daily cycles in the atmosphere may bias the estimate of S2 because the period of S2 is half of the solar day [13,14]. This leads to the uncertainty of GPS-estimated S2 tidal constituent being at 0.5 mm level.
RMS residual of about 0.5 mm between the GPS estimates and the OTL models. Figure 5 shows that the residual RMS of the S2 tidal constituent were all nearly 0.5 mm. This indicates that there are some station-related errors of around 0.5 mm, which arise from GPS-related error. GPS signal propagation error or non-tidal daily cycles in the atmosphere may bias the estimate of S2 because the period of S2 is half of the solar day [13,14]. This leads to the uncertainty of GPS-estimated S2 tidal constituent being at 0.5 mm level.

Conclusions
In this study, we used GPS kinematic PPP to estimate the vertical OTL displacements of eight semidiurnal and diurnal tidal constituents with 10-year GPS data from 12 GPS reference stations in Hong Kong. The systematic and residual phasor differences between the GPS estimates and nine OTL models were obtained by least squares methods in the complex domain. Our experimental results suggest that: For the M2, N2, O1, and Q1 tidal constituents, the RMS residual between the GPS KPPP estimates and the different OTL models is mainly due to the uncertainties of the OTL models, which exist as a systematic component. When the systematic differences are removed, the residual RMS is basically within 0.2 mm, which indicates that the uncertainty of the four KPPP-estimated tidal constituents is at the 0.2 mm level. The K1, P1, and K2 tidal constituents are affected by significant GPS-related common-mode error. After removing the systematic differences, the residual RMS was still significant, which indicates that the three tidal constituents are also affected by GPS-related residual error. The RMS residual of the S2 tidal constituent between the GPS estimates and the OTL models is mainly caused by GPS-related error. Overall the uncertainty of the GPS-estimated S2 tidal constituent is at the 0.5 mm level. These results indicate that GPS KPPP can provide constraints for the Earth's vertical tidal displacement at the 0.2 mm level for the M2, N2, O1, and Q1 tidal components, and the 0.5 mm level for the S2 tidal component.
Author Contributions: G.W. and Q.W. conceived and designed the experiments; G.W. and W.P. processed data; G.W. and Q.W. carried out the experimental analysis.

Conclusions
In this study, we used GPS kinematic PPP to estimate the vertical OTL displacements of eight semidiurnal and diurnal tidal constituents with 10-year GPS data from 12 GPS reference stations in Hong Kong. The systematic and residual phasor differences between the GPS estimates and nine OTL models were obtained by least squares methods in the complex domain. Our experimental results suggest that: For the M2, N2, O1, and Q1 tidal constituents, the RMS residual between the GPS KPPP estimates and the different OTL models is mainly due to the uncertainties of the OTL models, which exist as a systematic component. When the systematic differences are removed, the residual RMS is basically within 0.2 mm, which indicates that the uncertainty of the four KPPP-estimated tidal constituents is at the 0.2 mm level. The K1, P1, and K2 tidal constituents are affected by significant GPS-related common-mode error. After removing the systematic differences, the residual RMS was still significant, which indicates that the three tidal constituents are also affected by GPS-related residual error. The RMS residual of the S2 tidal constituent between the GPS estimates and the OTL models is mainly caused by GPS-related error. Overall the uncertainty of the GPS-estimated S2 tidal constituent is at the 0.5 mm level. These results indicate that GPS KPPP can provide constraints for the Earth's vertical tidal displacement at the 0.2 mm level for the M2, N2, O1, and Q1 tidal components, and the 0.5 mm level for the S2 tidal component.