Figure 1 shows the bathymetry based on the gridded bathymetric dataset from GEBCO (General Bathymetric Chart of the Oceans). The GEBCO gridded bathymetric dataset is a global terrain model of oceans and lands, providing elevation data in meters at 15 arc-second intervals. These grids can be downloaded or accessed through Web Map Services (https://www.gebco.net).

In this study, the tidal harmonic constants were obtained from 31 GPS observation data, which was provided by Yuan et al. (2013). The amplitudes and phases of M₂, S₂, N₂, K₂, K₁, O₁, P₁ and Q₁ constituents in the Center of Mass (CM) reference frame are listed in Table 1. Furthermore, the corresponding values in the Center of Figure (CF) reference frame are accurately listed in Table 2. Due to the proximity of some stations, their locations have been rounded to two decimal places, resulting in identical latitudes and longitudes.

The global VDL tide models used in our study are classified into two categories: one is the empirical model mainly relying on results obtained from satellites (e.g., GOT4.10, GOT4.8, and EOT11a); the other is the assimilation model (e.g., FES2014 and NAO.99b) which are constrained by empirical observations through different assimilation approaches. Table 3 lists the countries, resolutions, major constituents, and types of the five VDL tide models adopted. The details of these models for accuracy assessment are as follows.

The GOT4.10 and GOT4.8 models are developed at Goddard Space Flight Center (GSFC) in the United States and follow a long series of similar efforts starting with Schrama and Ray (Ray, 2013), which are the result of an empirical harmonic analysis of satellite altimetry relative to an adopted prior model. Of the modern models examined here GOT4.10 and GOT4.8 have a spatial resolution of 0.5° × 0.5° with a grid size of 720 × 361, covering latitudes from 90°S to 90°N and longitudes from 0°E–180°–0.5°W. These models provide information on the vertical displacement load tides of 10 constituents in the CM reference frame.

EOT11a is a global empirical ocean tide model obtained through residual analysis of multi-mission satellite altimeter data in 2011 (Savcenko and Bosch, 2012), which is developed at Deutsches Geodtisches Forschungs Institut (DGFI) in Germany. EOT11a is provided on the regular (1/8)° × (1/8)° grid with a grid size of 2 881 × 1 441, covering latitudes from 90°S to 90°N and longitudes from 0°E–180°–0°W. This model provides information on the vertical displacement load tides of 13 constituents in the CF reference frame.

FES2014 is the latest finite elements hydrodynamic model, which absorbs tide gauge observations and multi-mission altimeter data which is developed at French Tidal Group (Lyard et al., 2021). FES2014 is based on the resolution of spectral-configured shallow water hydrodynamic equations, and uses global finite element grids with increasing resolution in coastal and shallow water regions, with a spatial resolution of (1/16)° × (1/16)° and a grid size of 5760 × 2881, covering latitudes from 90°S to 90°N and longitudes from 0°E–180°–0.0625°W. This model provides information on the vertical displacement load tides of 33 constituents in the CM reference frame.

NAO.99b model is an assimilation model developed by the National Astronomical Observatory (NAO) in Japan (Matsumoto et al., 2000). It has a resolution of 0.5° × 0.5° with a grid size of 720 × 360, covering latitudes from 89.75°S to 89.75°N and longitudes from 0.25°E–180°–0.25°W. This model provides information on the vertical displacement load tides of 23 constituents in the CM reference frame.

Munk and Cartwright (1966) introduced the concept of tidal admittance when proposing the response methodology for tidal analysis. The tidal admittance $ {A_n} $ for a specific constituent n is defined as follows:

where $ A_n^* = {H_n}/{C_n} $, with $ {H_n} $ is amplitude and $ {C_n} $ is the coefficient of tidal generating force, and $ {g_n} $ is phase lag.

For the same constituent, the relative admittance of the equilibrium tide with respect to the constituent obtained from GPS analysis can be defined as

where,

where $ R_{\text{GPS/equ}}^* $ represents amplitude ratios, $ {r_{\text{GPS/equ}}} $ represents phase lag.

For any observation point $ \left(\lambda ,\varphi \right) $，the corresponding equilibrium tide $ {\zeta _{{\text{EQ}}}} $ can be expressed as

where $ m $ is the eight tidal constituents; $ {f_i} $ and $ {u_i} $ are nodal correction factors of the i-th component; $ {V_0} $ may be expressed by astronomy phase in the corresponding expansion term of tidal potential at $ t = 0 $ in Greenwich system of the i-th component; $ {\omega _i} $ is the frequency of the i-th tidal constituent; $ p $ is tide group number, $ p = 1 $ represents diurnal tide, $ p = 2 $ represents semi-diurnal tide; $ \lambda $ is longitude; $ \varphi $ is latitude; $ S = 0 $ represents Greenwich universal time, $ {H_{\text{equ}}}_i $ is the theoretical equilibrium amplitude of the i-th constituent corrected by the earth tide, $ {H_{\text{equ}}} $ can be expressed as

Our methodology is based on Pan et al. (2023) which is built on the principle of smoothness about tidal admittance. The method (Pan et al., 2023) is based on the different physical characteristics between non-astronomical tides and astronomical tides. The admittance in a narrow band is similar to astronomical tides, moreover, it can be represented as a smoothing function of the tidal frequency (Zetler, 1971). Through the quadratic interpolation, astronomical K₂ tide can be separated from GPS system errors. Due to the VDL tide is generated by the displacement deformation of OTL in the vertical direction, the tidal responses for the VDL tide inherit the nature of smooth tidal admittance. namely, the tidal admittance for the VDL tide can be represented by the approximate algorithms of smoothing functions of tidal frequencies.

The tidal constituents obtained from GPS gauge measured data were compared with the corresponding VDL tide model results.

To quantify the error of each tidal constituent between GPS observations and the corresponding VDL tide models, the mean absolute difference values were calculated using the following formula:

where H and g are amplitude and phase, the subscripts mod and obs represent the simulated values of the VDL tide models and the observation values of the GPS gauges, respectively. j represents the serial number of the GPS gauges, N is the number of GPS gauges used, j = 1, 2, ···, N.

Further, we used the Root-Mean-Square (RMS) value known as the standard deviation, to represent the overall deviation of the VDL tide models and GPS observation values, the RMS value was calculated using the following formula:

where a is the cosine component of the harmonic constant and b is the sine component of the harmonic constant, which are calculated as follows:

The RMS is the distance between the model simulated value and the actual observed value, representing the degree of deviation between the model simulated value and the actual observed value, while the relative degree of deviation between the RMS and the actual observed value can be expressed as the relative deviation $ \text{δ} $,

where,

where ${\overline{a}_{{\text{obs}}}}$ represents the average value over all the GPS stations of $ {a_{{\text{obs}},j}} $， ${\overline{b}_{{\text{obs}}}}$ represents the average value over all the GPS stations of $ {b_{{\text{obs}},j}} $.

We can also use the degree of fit ( $ {r^2} $) between the model values and the actual observed values,

where $ r $ is equivalent to the correlation coefficient in linear regression.

In addition, the Root-Square-Sum (RSS) value was used to quantify the accuracy of VDL tide models, the RSS value of the m major constituents was calculated by the following formula:

where $ m $ is the eight tidal constituents mentioned above.

According to the tidal theory, the tidal responses for the VDL tide inherit the nature of smooth tidal admittance, the tidal admittance for the VDL tide can be represented by the approximate algorithms of smoothing functions of tidal frequencies, which are here approximated by quadratic functions of tidal frequencies for the normalized amplitudes and phase lags. For the semi-diurnal tidal group, we use the tidal admittances of N₂, M₂, and S₂ constituents and the quadratic interpolation method to determine the interpolated quadratic curves. It is worth noting that K₂ tidal admittances which are contaminated by GPS-system errors are not used. Due to typicality in the poor estimation of K₂ tidal parameters associated with GPS-system errors, Mald (4.19°N, 73.53°E) and Pre2 (25.75°S, 28.22°E) are discussed in detail as the typical examples. As shown in Figs 2 and 3 at Mald GPS station, we calculate the astronomical K₂ admittances (red dots in Figs 2 and 3) by using the quadratic interpolation method for interpolation. Due to the known equilibrium amplitude of the K₂ constituent (as shown in Section 3.1), the amplitude and phase lag of the astronomical K₂ constituent at Mald GPS station can be figured out as 2.66 mm and 104.4° in the CM frame, as 2.63 mm and 103.2° in the CF frame, respectively.

There are two main energy sources for the K₂ tide obtained from GPS observations. One is a consequence of the astronomical K₂ tide generated by astronomical factors, and the other is a consequence of the non-astronomical K₂ tide majorly induced by GPS-system errors. Therefore, Fig. 4 shows the Mald GPS station as an example, the observed K₂ tide (represented by the black arrow) is composed of the vector sum of astronomical K₂ tide (represented by the red arrow) and non-astronomical K₂ tide (represented by the blue arrow). The observed K₂ amplitude and phase lag at the Mald GPS station are 2.41 mm and 52.0° in the CM frame, and 2.35 mm and 49.1° in the CF frame, respectively. The amplitude and phase lag of the astronomical K₂ constituent are 2.66 mm and 104.4° in the CM frame, and 2.63 mm and 103.2° in the CF frame, respectively. Figure 4 clearly indicates that the astronomical K₂ tide is highly consistent with the modeled K₂ tide (represented by the green arrow). While the amplitude and phase lag of non-astronomical K₂ constituent are 2.25 mm and 342.4° in the CM frame, and 2.28 mm and 339.9° in the CF frame, respectively. It is shown that although the tidal frequencies of non-astronomical K₂ tide and astronomical K₂ tide are the same, the amplitudes and phase lags characteristics are significantly different due to their different origins and processes.

Similarly, Figs 5 and 6 display the tidal admittances of the semi-diurnal VDL tide at the Pre2 GPS station in the CM frame and in the CF frame. The calculation principle is the same with the Mald GPS station. Via the interpolation, the amplitude and phase lag of the astronomical K₂ constituent at the Pre2 GPS station can be obtained as 2.37 mm and 252.9° in the CM frame, as 2.34 mm and 255.5° in the CF frame, respectively. Figure 7 displays the observed K₂ amplitude and phase lag at the Pre2 GPS station are 1.78 mm and 232.9° in the CM frame, and 1.85 mm and 236.6° in the CF frame, respectively. The amplitude and phase lag of the astronomical K₂ constituent are 2.37 mm and 252.9° in the CM frame, and 2.34 mm and 255.5° in the CF frame, respectively. Figure 7 clearly indicates that the astronomical K₂ tide is highly consistent with the modeled K₂ tide (represented by the green arrow). The amplitude and phase lag of non-astronomical K₂ constituent are 0.92 mm and 114.2° in the CM frame, and 0.84 mm and 120.9° in the CF frame, respectively.

Further, we correct K₂ tide for OTL vertical displacement estimates at 31 GPS stations (Fig. 1) in the CM and CF frame which was provided by Yuan et al. (2013). The corrected harmonic constants of K₂ at 31 GPS stations for the equatorial and Indian Ocean are shown in Table 4.

The amplitude and phase of the GPS gauge considered in this section were obtained from 31 stations (Fig. 1). The harmonic constants of 31 stations were extracted from each tide model using the spline interpolation method. Then, the harmonic constants for M₂, S₂, N₂, K₂, K₁, O₁, P₁, and Q₁ tidal constituents of tide models (FES2014, EOT11a, GOT4.10c, GOT4.8, and NAO.99b) were compared with corresponding values of 31 GPS gauge stations (Section 2.2).

The comparison results of the five VDL tide models are listed in Table 5. As is described in Table 5, the discrepancies between the VDL tide models and GPS measurements were large for the positions located along the continental coastlines, with the RMS values ranging from 1.30–3.51 mm for the S₂ and K₁ tidal constituents and 0.39–1.79 mm for the M₂ and O₁ tidal constituents. The remaining constituents (i.e., N₂, K₂, Q₁, and P₁) had comparatively lower RMS values of 0.13–1.22 mm. In detail, the FES2014 model had the best results for M₂, N₂, and P₁ constituents for the equatorial and Indian Ocean with the RMS are 0.61 mm, 0.21 mm, and 0.67 mm, respectively; the EOT11a model had the best results for S₂, K₂, O₁, and Q₁ constituents for the equatorial and Indian Ocean with the RMS are 1.33 mm, 0.39 mm, 0.39 mm, and 0.13 mm, respectively; the GOT4.10c model had the best results for K₁ constituent for the equatorial and Indian Ocean with the RMS is 1.30 mm. Overall, EOT11a and FES2014 models exhibited the smaller RSS value (2.29 mm and 2.34 mm) of the five VDL tide models, indicating a slightly superior performance.

The VDL cotidal charts for the eight major tidal constituents derived from the EOT11a tide model, which accuracy was relatively high in the equatorial and Indian Ocean regions, are presented in Fig. 8. Notably, the M₂ VDL tide exhibits a peak amplitude exceeding 45 mm in the Mozambique Channel, followed by a secondary maximum of over 40 mm along the northwest coast of Australia, and a third maximum exceeding 35 mm in the central Indian Ocean. As expected, the VDL amplitudes on land are comparatively lower due to the absence of load effects. The distribution patterns of the S₂ VDL tide closely resemble those of the M₂ VDL tide. However, given that the amplitudes of the S₂ ocean tide are inherently smaller than those of the M₂ ocean tide, the S₂ VDL amplitudes are correspondingly lower. The first two maxima amplitudes of the S₂ VDL tide, occurring in the Mozambique Channel and along the northwest coast of Australia, have magnitudes exceeding 23 mm, which are nonetheless smaller than those of the M₂ VDL tide. Similarly, a third maximum with amplitudes over 21 mm is observed in the middle of the Indian Ocean. Phase distributions of both the S₂ and M₂ VDL tides exhibit similarities, with a notable difference being the presence of an amphidromic point for the S₂ VDL tide on the west coast of Australia, whereas the M₂ VDL tide displays a degenerated amphidromic point near the southern Australian landmass. The distribution characteristics of the N₂ VDL tide mirror those of the M₂ VDL tide, with similarly reduced amplitudes due to the smaller amplitudes of the N₂ ocean tide. Furthermore, the K₂ VDL tide exhibits a distribution pattern similar to that of the S₂ VDL tide. Notably, the K₁ VDL tide reaches its peak amplitude in the central Arabian Sea, exceeding 18 mm, while a secondary maximum, surpassing 16 mm, is observed along the northwest coast of Australia. Similarly, the O₁ VDL tide exhibits distribution characteristics akin to the K₁ VDL tide, with the two highest amplitudes occurring in the same regions, albeit with magnitudes of over 9 mm and 10 mm, respectively. However, the amplitudes of the P₁ VDL tide are relatively smaller than those of the O₁ VDL tide, reflecting the corresponding reduction in amplitude compared to the O₁ ocean tide. Likewise, the Q₁ VDL tide shares similar distribution characteristics with the O₁ VDL tide.

Due to various factors such as inherent errors in the GPS system, errors in satellite signal transmission, the complexity of tidal phenomena, and limitations in data processing and analysis methods, tidal constants derived from GPS observations inevitably contain relatively large errors (Table 1). As a consequence, RMSs between GPS results and tidal models have large uncertainties (Table 5), which may exert potential influences on the results of accuracy evaluation. Furthermore, when evaluating the accuracy of VDL models, differences in the number of stations may have potential impacts on the evaluation results, which are all factors that need to be taken into account. It is believed that more effective methods to reduce uncertainties in GPS results will be proposed in the future with the advancement of observation technology.