Tide simulation in a global eddy-resolving ocean model

1. Introduction

Diapycnal mixing is of great significance for global ocean stratification and driving meridional overturning circulation (Marshall and Speer, 2012), which can significantly affect local and global climate change through the diapycnal transport of matter and heat (Kostov et al., 2014). This process can have a profound impact on regional and global climate systems. The wind-driven upwelling of the Southern Ocean and the diapycnal mixing of the Pacific and Indian Oceans are considered two crucial driving processes (Talley, 2013), while the tidal process in the open ocean is the main source of mixing energy (Munk and Wunsch, 1998). Munk and Wunsch (1998) analyzed the tidal energy budget based on the altimeter results of Egbert (1997). The results showed that the Sun-Moon-Earth system delivered a total of 3.5 TW of energy to the ocean, of which 2.6 TW dissipated in the shallow water in the form of barotropic tides, and the remaining 0.9 TW was dissipated in the deep ocean. Deep ocean dissipation could be divided into local dissipation and remote dissipation. Nevertheless, local and remote dissipation could provide considerable energy for maintaining stratification and driving thermohaline circulation in the deep sea. Egbert and Ray (2000), Jayne and St. Laurent (2001), and Niwa and Hibiya (2014) estimated that the magnitude of dissipation due to tides in the deep ocean was approximately 1 TW. Therefore, studying tides and their dissipation is crucial for understanding ocean circulation.

Due to limitations in computing power and theoretical understanding, early tidal and global ocean circulation simulations were conducted separately. Tidal numerical models are mainly based on a two-dimensional numerical model of the tidal wave dynamic equation. In contrast, ocean general circulation models (OGCMs) are based on simplified Navier‒Stokes equations without considering tidal processes. However, with advancements in the development of satellite altimeter technology, two-dimensional tidal models have undergone significant improvements, resulting in relatively mature prediction and simulation capabilities (Stammer et al., 2014). Tidal data that closely match observations can be obtained by integrating two-dimensional tidal models with data assimilation (Shum et al., 1997). The two-dimensional tidal model can be refined by incorporating altimeter data, resulting in an increasingly precise tidal map. Nevertheless, several issues remain unresolved, including the correlation between tides and climate change, the interplay between tides and large-scale ocean circulation, the mechanism behind internal tidal dissipation in the open ocean, and the prediction of internal tides. Addressing these questions necessitates the concurrent simulation of global ocean circulation and tides (Arbic et al., 2018).

There are usually two methods for incorporating tides into OGCMs. The first approach is implicit, where the impacts of tides on the ocean, such as mixing induced by internal tides, are parameterized (St. Laurent et al., 2002). This parameterization enhances mixing in the open ocean and makes the coefficient of diapycnal mixing not uniform in the horizontal direction. Simmons et al. (2004) added the parameterization scheme of St. Laurent et al. (2002) to a low resolution OGCM. The scholars found that adding the parameterization of tidal mixing reduces the deviation of ocean simulation. Yu et al. (2017) showed that this parameterization significantly enhances the Atlantic meridional overturning circulation (AMOC). Although several works have shown that the parameterization of tidal mixing has effects on improving ocean circulation simulation (Jayne, 2009; Melet et al., 2013; Saenko and Merryfield, 2005; Simmons et al., 2004; Song et al., 2023; Yu et al., 2017), uncertainty remains a problem in the parameterization scheme. Changes in some coefficients in this scheme can significantly impact the results of ocean model simulations (Jayne, 2009). The second approach is explicit, where the tidal potential is directly introduced into the sea surface height formula of the ocean model. Thomas et al. (2001) and Schiller (2004) incorporated the tide explicitly into an OGCM. Nevertheless, due to some numerical problems with the model, Schiller (2004) could only strictly limit the tidal potential in the Indonesian sea. Schiller and Fiedler (2007) successfully incorporated the tidal potential of the eight constituents into an OGCM, and the results showed that the inclusion of tides enhanced the simulation of temperature and salinity in the Indonesian throughflow. Müller et al. (2010) discussed the impact of incorporating tidal potential into a climate model and demonstrated that the temperature simulation in the North Atlantic was improved. It is worth mentioning that these studies were conducted with low resolution OGCMs (Müller et al., 2010; Schiller, 2004; Schiller and Fiedler, 2007; Thomas et al., 2001).

Arbic et al. (2010) were the earliest researchers to incorporate explicit tides into an eddy-resolving global OGCM. The scholars incorporated the tidal gravitational potential of eight major tidal constituents in the hybrid coordinated ocean model (HYCOM) with realistic ocean stratification. The parameterized topographic wave drag was considered to simulate the tides accurately. The simulated barotropic tides were found to be as accurate as previous tidal models without data assimilation (Stammer et al., 2014). Since then, global eddy-resolving OGCMs that incorporate atmospheric and tidal forcing have emerged successively, such as STORMTIDE (Müller et al., 2012), Massachusetts Institute of Technology General Circulation Model (MITgcm) (Rocha et al., 2016a, b), and Icosahedral Nonhydrostatic Weather and Climate Model (ICON-O) (Logemann et al., 2021; von Storch et al., 2023). In addition, some researchers have considered tides using boundary conditions in regional high resolution ocean models (Peng et al., 2021; Siyanbola et al., 2023; Wang et al., 2016). In future studies of eddy-resolving ocean models, concurrent atmospheric and tidal forcings (wind plus tide) should be increasingly utilized (Arbic et al., 2018).

Global eddy-resolving OGCMs prioritize high-frequency movements in the ocean, the reciprocal impacts of internal tides and ocean circulation, and the prediction of internal tides (Arbic et al., 2012, 2018; Li and von Storch, 2020; Müller et al., 2012, 2014, 2015; Ansong et al., 2015, 2017; Buijsman et al., 2016; Savage et al., 2017a, b). However, evaluating the simulation of barotropic tides in a global eddy-resolving OGCM is essential before studying related physical processes. Arbic et al. (2010) reported an error of 8.26 cm for the M₂ constituent, which they found through a comparison of the 102 open-ocean tide observations made by Shum et al. (1997); this value is equivalent to the tides simulated by Arbic et al. (2004). Müller et al. (2014) conducted the same evaluation, and the error of the M₂ constituent reached 8.19 cm. Additionally, seasonal variability in the amplitude of the M₂ constituent was found. von Storch et al. (2023) found slightly larger errors in simulated barotropic tides than Arbic et al. (2010) and Müller et al. (2014). The researchers noted that improving the horizontal resolution of OGCM did not necessarily reduce open-ocean tide simulation errors.

We have performed some related studies with the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics/Institute of Atmospheric Physics (LASG/IAP) Climate System Ocean Model (LICOM). Yu et al. (2016) incorporated tidal forcing into LICOM2.0 with a 1° horizontal resolution, analyzed its impact on ocean circulation simulation, and investigated its effect on climate via coupled models (Yu et al., 2020). Jin et al. (2022) adopted a new formulation of an explicit tidal scheme and successfully tested it in LICOM2.0 with the same resolution. Furthermore, the global eddy-resolving version of LICOM3.0 has undergone preliminary evaluation without tidal forcing (Li et al., 2020). In addition, incorporation and long-term tests of tidal forcing have been recently completed. In the present study, we aim to incorporate tidal forcing into the eddy-resolving model LICOM3.0 and systematically evaluate the simulated barotropic tides, laying a foundation for future research. Section 2 details the model used, the model configuration, the tidal forcing incorporation method, and the error analysis method. Evaluations of the errors of barotropic tides are described in Sections 3 and 4, providing the conclusion of this article and the future outlook.

2. Model, experiment, and methods

2.1 Astronomical tidal forcing

Based on the theory of equilibrium tides (Pugh, 1996), Yu et al. (2016) incorporated the tidal forcing of eight major constituents into LICOM2.0, including four major diurnal constituents (K₁, P₁, Q₁, and O₁) and four major semidiurnal constituents (M₂, S₂, N₂, and K₂). In this article, based on Yu et al. (2016) and Arbic et al. (2018), we incorporated the tidal forcing of eight constituents into LICOM3.0. Ocean tides were generated by the gravitational forces of the moon and the sun on seawater. The diurnal tide $ {\eta }_{\mathrm{EQ}1} $ and the semidiurnal tide $ {\eta }_{\mathrm{EQ}2} $ could be calculated by Eqs (1) and (2), respectively, as follows:

where the equilibrium tide $ {\eta }_{\mathrm{E}\mathrm{Q}} $ is the sum of $ {\eta }_{\mathrm{E}\mathrm{Q}1} $ and $ {\eta }_{\mathrm{E}\mathrm{Q}2} $, $ \phi $ is the latitude, $ \lambda $ is the longitude, $ t $ is the time, $ {t}_{\mathrm{r}\mathrm{e}\mathrm{f}} $ is the reference time of the model, $ {V}_{0}({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}) $ is the initial phase of the constituents at the reference time, $ A $ and $ \omega $ are the amplitude and frequency of the equilibrium tide for each constituent, respectively, $ f({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}) $ and $ \mu ({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}) $ are the nodal factor and nodal angle, respectively, representing the modulations of the amplitude and delay angle of the tides over 18.61 years (Pugh, 1996), and $ {h}_{2} $ and $ {k}_{2} $ are Love numbers, which explain the deformation of the solid Earth caused by astronomical forcing and the influence of the redistribution of Earth’s mass on the gravitational potential, respectively. In this study, $ 1+{k}_{2}-{h}_{2} $ is taken as a constant value of 0.7.

The nodal factor $ f\left({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}\right) $ and the phase-dependent values $ {V}_{0}({t}_{\mathrm{r}\mathrm{e}\mathrm{f}})+\mu ({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}) $ are only related to astronomical variables. The formula for $ {V}_{0}({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}) $ is

where n₀−n₆ are Doodson numbers, found in Appendix Table 4 in Zuo et al. (2018). $ {\tau }_{0} $, $ {s}_{0} $, $ {h}_{0} $, $ {p}_{0} $, $ {N}_{0} $, and $ {p}_{0}' $ are mean lunar time, mean longitude of moon, mean longitude of sun, mean longitude of lunar perigee, mean longitude of lunar ascending node, and mean longitude of solar perigee, respectively. The formula for $ f\left({t}_{\mathrm{ref}}\right) $ and $ \mu ({t}_{\mathrm{ref}}) $ are shown in Table 1. The astronomical parameters are calculated by

Table  1.  Calculation formulas for the nodal factor f and the nodal angle μ

Constituent	f	μ/(°)
K1	1.0060 + 0.1150 cos$ {N}_{0} $ – 0.0088 cos(2$ {N}_{0} $) + 0.0006 cos(3$ {N}_{0} $)	–8.86° sin$ {N}_{0} $ + 0.68° sin(2$ {N}_{0} $) − 0.07° sin(3$ {N}_{0} $)
O1	1.0089 + 0.1871 cos$ {N}_{0} $ − 0.0147 cos(2$ {N}_{0} $) + 0.0014 cos(3$ {N}_{0} $)	10.80° sin$ {N}_{0} $ − 1.34° sin(2$ {N}_{0} $) + 0.19° sin(3$ {N}_{0} $)
P1	1	0
Q1	same as O1	same as O1
M2	1.0004 − 0.0373 cos$ {N}_{0} $ + 0.0003 cos(2$ {N}_{0} $)	–2.14° sin$ {N}_{0} $
S2	1	0
N2	same as M2	same as M2
K2	1.0241 + 0.2863 cos$ {N}_{0} $ + 0.0083 cos(2$ {N}_{0} $) − 0.0015 cos(3$ {N}_{0} $)	–17.74° sin$ {{N}}_{0} $ + 0.68° sin(2$ {N}_{0} $) − 0.04° sin(3$ {N}_{0} $)

 | Show Table

DownLoad: CSV

where T is Julian century since January 1, 1900, T = $ [365(Y\ - 1\,900)+(D-1)+i]/{36\,525} $, $ Y $ is the year of model time, $ D $ is the day of the model time in year $ Y $, and $ i $ is the number of leap years from 1900 to year $ Y $. In the subsequent analysis, we compared the nodal factor $ f\left({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}\right) $ and the phase-dependent $ {V}_{0}({t}_{\mathrm{r}\mathrm{e}\mathrm{f}})+ \mu ({t}_{\mathrm{r}\mathrm{e}\mathrm{f}}) $ values of the theoretical and LICOM3.0 model results; moreover, the equilibrium tide $ {\eta }_{\mathrm{E}\mathrm{Q}} $ is shown.

By incorporating the sea surface height caused by tides into the sea surface height $ \eta $ of the barotropic equation in the model, the gradient scheme for the sea surface height became the following formula:

where the coefficient $ \alpha $ represents the self-attraction and loading term (Hendershott, 1972). In this article, $ \alpha $ is taken as 0.948.

The tide can change with the relative position of the Sun, Moon, and Earth. The tide is the strongest when the three celestial bodies are in a straight line, where it is called a spring tide. When the angle between the straight line connecting the Sun and the Moon to the Earth is a right angle, the tide is relatively weak, and it is called a neap tide.

To ensure model stability, we added the parameterized topographic wave drag term of Jayne and St. Laurent (2001), while in the work of Yu et al. (2016), only the bottom friction term was implemented. The formula for the wave drag term used was as follows:

where c is a coefficient and can be adjusted, $ {N}_{0{\mathrm{b}}} $ represents the buoyancy frequency at the bottom boundary, $ {U}_{0} $ denotes the barotropic velocity in the direction of the reference drag layer, $ {N}_{0{\mathrm{b}}} $ and $ {U}_{0} $ are computed in LICOM3.0 with model density and velocities, $ \kappa $ denotes the topographic wavenumber and is set to 2π/10 km (Jayne and St. Laurent, 2001), and $ {h}^{2} $ represents the topographic roughness and is computed using the (1/30)° Earth Topography v2 (ETOPO2; National Centers for Environmental Information, 2006) grid as the average squared depth anomaly over a 1° square for a plane fit. Arbic et al. (2010) and Buijsman et al. (2020) included a parameterized topographic wave drag term and found that the wave drag term had a certain impact on the accuracy and energy distribution of the simulated tides.

2.2 Model and experiment

The model used in this study was LICOM (LASG/IAP Climate System Ocean Model). This model is an OGCM developed by the National Key Laboratory for Numerical Simulation of Atmospheric Science and Earth Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences, China. The version used in this paper was the eddy-resolving LICOM3.0 (Li et al., 2020). The horizontal grid of LICOM3.0 was a tripolar grid (Murray, 1996), which could effectively address the singularity issues of latitude and longitude grids. The eddy-resolving version had a resolution of (1/10)°, with 55 layers of $ \eta $ coordinates vertically. The sea ice model used was Community Ice Code version 4 (CICE4).

To study the influence of wave drag coefficient c on the accuracy of barotropic tides simulation, c is set to 0.5, 1.0, and 1.5 in the wave drag scheme; the three experiments are called wd0.5, wd1, and wd1.5. The astronomical tidal forcing is added to these three experiments, starting from January 1, 2015, to March 31, 2015, using the results of the phase 2 of the Ocean Model Intercomparison Project (OMIP2) experiment (Li et al., 2020) as initial values. The tide has gone through dozens of cycles in three months, so the barotropic tides are approximately stable. The sea surface height of the three experiments is conducted for harmonic analysis of tides. The results can be seen in Table 2. As c increases, the amplitude and phase errors increase of the four major constituents, as well as the total error. Due to the high cost of conducting high-resolution global experiments with hourly output, c = 0.5 is considered the relatively optimal parameter so far. So, in the following experiment, c is set to 0.5.

Table  2.  Total (d), amplitude (d_a), and phase (d_p) errors of the four major constituents of LICOM3.0 relative to TPXO

Experiment	Error/cm
	M2				S2				K1				O1
	d	da	dp		d	da	dp		d	da	dp		d	da	dp
wd0.5	9.86	7.21	5.37		5.12	4.00	2.47		3.58	2.36	2.25		2.59	1.26	1.99
wd1	10.70	8.16	5.45		5.38	4.25	2.57		3.77	2.49	2.36		2.57	1.25	1.94
wd1.5	11.39	8.88	5.59		5.59	4.43	2.67		3.92	2.59	2.45		2.57	1.26	1.91

 | Show Table

DownLoad: CSV

Based on the OMIP2 experiment (Li et al., 2020) using JRA55-do (Tsujino et al., 2018), the astronomical tidal forcing in the form of tidal potential was further incorporated, and the wind plus tide experiment, as defined by Arbic et al. (2018), was conducted. The standard OMIP2 experiment covered the period from 1958 to 2018, while the tidal experiment was integrated starting from January 1, 2015, to December 31, 2016, using the results of the OMIP2 experiment as initial values. The data from 2016 were analyzed in this paper. For tidal harmonic analysis, the data output interval of the tidal experiment was one hour.

2.3 Harmonic analysis and error analysis

To evaluate the simulated tides, we used the products of Topex Poseidon Cross-Overs v9 (TPXO9; Egbert and Erofeeva, 2002), which is the latest product of the global barotropic ocean tidal model TPXO, in this paper. We interpolated these results to a grid with a horizontal resolution of (1/10)°. In addition, we used 102 open-ocean tide observations (hereafter referred to as st102) from Shum et al. (1997) and data from two tide gauge stations, Yap (9.508 3°N, 138.128 3°E) and Kodiak (57.731 7°N, 152.511 7°W), selected from World Ocean Circulation Experiment (WOCE; Ponchaut et al., 2001). Moreover, we used the Python version of UTIDE (Codiga, 2011) to perform the harmonic analysis of the sea surface height of the tide experiment and obtained the average amplitude and phase values of the eight constituents.

Based on the definition of Arbic et al. (2004), the square of the total error of tides $ {d}^{2} $ between LICOM3.0 and the observation data at a certain grid point could be defined as the sum of the square of the amplitude error $ {d}_{\mathrm{a}}^{2} $ and the square of the phase error $ {d}_{\mathrm{p}}^{2} $:

where $ {A}_{\text{model}} $ and $ {A}_{\text{obs}} $ represent the amplitude of the model and observation of a certain constituent, respectively, and $ {\phi }_{\text{model}} $ and $ {\phi }_{\mathrm{o}\mathrm{b}\mathrm{s}} $ represent the phase of the model and observation of a certain constituent, respectively. Obviously, the corresponding square of the regional or grid point average errors are as follows:

where $ \overline{(\cdot )} $ indicates an average over grid points or area.

To evaluate the simulations of all constituents, the Root Square Sum (RSS) values of all constituent errors were considered. The percentage of sea surface height variance captured was defined as follows, where RSS² is the sum of the square of average error D² of all constituents and S is the signal of tides:

when the percentage of sea surface height variance captured (V_capt) was 100%, the model perfectly simulated the barotropic tides. Therefore, the higher this value was, the more accurate the model simulation of the tides.

3. Results

3.1 Tidal forcing

First, we verified whether the phase-dependent $ {V}_{0}+\mu $ and nodal factor f in LICOM3.0 are correct. We output these two items from the module of tidal forcing and compare them with the results of Fisheries and Oceans Canada (2017, 2018) (DFO) and the tidal harmonic analysis program UTIDE. We recorded the results on January 1 and July 1, as shown in Table 3. We found few differences between phase-dependent $ {V}_{0}+\mu $ and nodal factor f in the DFO and UTIDE. Excluding the Q₁ constituent, the phase-dependent $ {V}_{0}+\mu $ simulated by LICOM3.0 is within 0.31° of DFO and UTIDE, and the difference in the Q₁ constituent is slightly larger at approximately 1.632°. The relative error of nodal factor f simulated by LICOM3.0 compared to DFO and UTIDE is less than 1.22%. Therefore, the phase-dependent $ {V}_{0}+\mu $ and nodal factor f calculated by the tidal module in LICOM3.0 are basically correct when compared with the references, with no significant differences.

Table  3.  Phases-dependent $ {V}_{0}+\mu $ and nodal factor f values of each constituent calculated by DFO, tidal harmonic analysis program UTIDE and LICOM3.0 tidal module in 2016, corresponding to the dates of January 1 and July 1 of that year

Constituent	($ {V}_{0}+\mu $) (Jan. 1)/(°)				f (Jul. 1)
	DFO	UTIDE	LICOM3.0		DFO	UTIDE	LICOM3.0
M2	210.745	210.729	210.857		1.037	1.037	1.037
K1	9.242	9.227	9.299		0.886	0.886	0.886
O1	202.050	202.106	202.000		0.809	0.809	0.813
S2	0.001	0	0		0.998	0.998	1.000
P1	349.847	349.852	349.903		1.011	1.011	1.000
Q1	43.050	43.384	41.752		0.821	0.823	0.813
N2	50.690	50.690	50.610		1.037	1.037	1.037
K2	198.451	198.412	198.722		0.752	0.752	0.755

 | Show Table

DownLoad: CSV

The patterns of the daily changes of the spring and neap tides calculated by LICOM3.0 are shown in Fig. 1, where four time points on January 10 and January 18, 2016, are selected, respectively, at 04:00, 10:00, 16:00, 22:00, and 00:00, 06:00, 12:00, 18:00 (UTC+8). There are obvious differences in the amplitude of sea surface height between spring tides and neap tides. The maximum tide height during a spring tide is above 0.4 m, while the maximum tide height during a neap tide is only half of that during a spring tide, i.e., approximately 0.2 m. Both tides exhibit distinct characteristics of semidiurnal and diurnal cycles, with a phase shift from east to west. The change in tides in high-latitude areas is minimal. From west to east, the sequence of spring tides and neap tides follows a positive-negative-positive-negative pattern in latitude, which is consistent with the findings of the EXP1 experiment by Jin et al. (2022) and is aligned with the actual situation regarding ocean tides.

Figure  1.  Changes in the sea surface caused by the equilibrium tide calculated by the LICOM3.0 tidal module: the results during a spring tide at 04:00, 10:00, 16:00, and 22:00 (a–d) and the results during a neap tide at 00:00, 06:00, 12:00, and 18:00 (e–h). Time in UTC+0.

DownLoad: Full-Size Img PowerPoint

3.2 Tides simulated by LICOM3.0 and TPXO

We mainly analyzed the largest semidiurnal constituent, M₂, and the largest diurnal constituent, K₁. Figure 2 shows the amplitudes and phase of M₂ and K₁ constituents from LICOM3.0 and TPXO. Overall, the geographical distribution of the amplitude extreme values is approximately the same, and the amphidromic point locations can be mapped. For M₂, the amplitude of the tide simulated by LICOM3.0 is larger than that of TPXO. The global average M₂ amplitudes of LICOM3.0 and TPXO are 37.60 cm and 30.70 cm, respectively. The large amplitude values for TPXO displays are mostly distributed on the continent or the periphery of islands, such as the Bay of Alaska, the Weddell Sea, and the North Atlantic. In addition, in these regions, the amplitude of the M₂ constituent simulated by LICOM3.0 is relatively large, which accounts for the high global average amplitude. The phase of the M₂ constituent differs from TPXO in some areas, such as the Tasman Sea and the Southern Ocean. The gap is obvious in the Southern Ocean, which may be related to the absence of LICOM3.0 data at high latitudes, resulting in a northward shift of the amphidromic points. In contrast to M₂, large amplitude values for K₁ are mainly distributed along the entire coast of Antarctica and around Indonesia.

Figure  2.  Spatial patterns of the amplitude (colored, cm) and phase (contour, °) characteristics of the M₂ constituent for TPXO (a) and LICOM3.0 (b) and of the K₁ constituent for TPXO (c) and LICOM3.0 (d). The black solid lines are separated by 30°.

DownLoad: Full-Size Img PowerPoint

Similarly, the amplitude of the K₁ constituent simulated by LICOM3.0 is larger than that of TPXO, with global average amplitudes of 11.11 cm and 10.2 cm, respectively. From a spatial distribution perspective, the simulated tides surrounding the Antarctic mainland are particularly strong in LICOM3.0. The position of the amphidromic points simulated by LICOM3.0 is consistent with those simulated by TPXO, with significant differences in the North Atlantic and Indian oceans.

Notably, in both the M₂ and K₁ constituents, the phase contour lines simulated by LICOM3.0 display obvious perturbations in regions such as the North Pacific and the tropical Pacific. These perturbations are more obvious in areas where the internal tides are active. However, there is no such phenomenon in TPXO. These perturbations represent baroclinic tides, which signal internal tides. As TPXO is a two-dimensional barotropic model without internal tide processes, no such signals exist. The horizontal resolution of the LICOM3.0 eddy-resolving OGCM is approximately 10 km, which is much smaller than the shortest wavelength of mode-1 M₂ and K₁ internal tides (Carrere et al., 2021). LICOM3.0 is sufficient to partially resolve the M₂ and K₁ internal tides, resulting in such perturbations in the phase. A similar phenomenon appeared in the results reported by Arbic et al. (2012).

Figure 3 displays the amplitude error $ {d}_{\mathrm{a}} $, phase error $ {d}_{\mathrm{p}} $, and total error $ d $ values of the M₂ and K₁ constituents. Table 4 presents the global average results. The results in Table 4 only comprise an ocean between 66°S and 66°N and a depth greater than 1 000 m. For the M₂ constituent, the global average total error of LICOM3.0 is 10.19 cm relative to TPXO9, with an amplitude error of 6.83 cm and a phase error of 6.09 cm. The amplitude error accounts for a relatively large proportion. From a spatial distribution perspective, the maximum value of the M₂ constituent total error appears near Greenland in the North Atlantic and the Ross Sea, and the error is large near other semi-enclosed bays. In areas where the total error is large, the amplitude and phase errors are both large, but the amplitude error is even larger. In such areas, the amplitude of the M₂ constituent is the largest in both TPXO9 and LICOM3.0, and the error in LICOM3.0 is the largest. For the K₁ constituent, the global average total error of LICOM3.0 is 2.85 cm relative to TPXO9, with an amplitude error of 1.66 cm and a phase error of 1.92 cm. The phase error accounts for a relatively large proportion. From a spatial distribution perspective, the total error is large near the Southern Ocean. However, the phase error contributes more to the total error in the K₁ constituent than in the M₂ constituent. ICON-O (von Storch et al., 2023) is a recently developed OGCM with tidal forcing, and we compared the results of LICOM3.0 and ICON-O. Table 4 shows the results of ICON-O (von Storch et al., 2023). A comparison shows that the results of LICOM3.0 are significantly better than those of ICON-O, both in terms of the total error and the errors of the two components. The total error of the two constituents is more than 50% lower in LICOM3.0. In addition, it is worth noting that, unlike LICOM3.0, the amplitude error is larger than the phase error of the M₂ constituent in ICON-O results, while amplitude and phase errors contribute almost equally in the K₁ constituent.

Figure  3.  Errors of the simulated M₂ constituent (a–c) and K₁ constituent (d–f) relative to TPXO in LICOM3.0: amplitude error (a, d), phase error (b, e), and total error (c, f).

DownLoad: Full-Size Img PowerPoint

Table  4.  Errors relative to TPXO in LICOM3.0 and ICON-O (von Storch et al., 2023), including total error $ d $, amplitude error $ {d}_{\mathrm{a}} $, and phase error $ {d}_{\mathrm{p}} $

Constituent	Error/cm
	LICOM3.0				ICON-O
	$ d $	$ {d}_{\mathrm{a}} $	$ {d}_{\mathrm{p}} $		$ d $	$ {d}_{\mathrm{a}} $	$ {d}_{\mathrm{p}} $
M2	10.19	6.83	6.09		14.25	10.63	9.39
K1	2.85	1.66	1.92		4.60	3.26	3.25
O1	4.30	2.55	2.98		5.19	3.89	3.43
S2	5.37	3.12	3.84		7.90	4.71	1.50
P1	1.12	0.52	0.86		1.49	1.08	1.02
Q1	0.77	0.47	0.56		0.86	0.48	0.72
N2	1.51	0.87	1.06		2.30	1.50	1.74
K2	1.19	0.78	0.73		2.30	1.19	1.96

 | Show Table

DownLoad: CSV

We have also investigated other six constituents (not shown). We found that the results for other constituents are similar to the most dominant semidiurnal (M₂) and diurnal (K₁) constituents.

In summary, we found that the errors of simulated M₂ and K₁ constituents vary in different sea areas in LICOM3.0, and the amplitude and phase contribute differently to the total error. The results from different OGCMs are different, and it is impossible to identify any apparent commonalities. Further analysis of results from additional models is needed.

3.3 LICOM3.0 simulated tides and observed tides

We further utilized the st102 data from Shum et al. (1997) and the data from two tidal gauge stations of World Ocean Circulation Experiment (WOCE) for validation. These observational data are often used for tide simulation verification. Table 5 shows the bias of the eight major constituents in the LICOM3.0 and st102 (Shum et al., 1997) dataset. Unlike TPXO, the data from these 102 stations are observation data, which, although they do not cover as wide a range as TPXO, are representative of the true tides. The comparison shows that the simulation error is still the largest for the M₂ constituent, with a total error of 13.85 cm and a relative error of 41.7%, where the amplitude contributes more to the total error. The total error for the K₁ constituent is 3.44 cm, and the relative error is 30.6%, with the phase contributing more to the total error. Overall, the relative sizes of the total errors for the M₂ and K₁ constituents and the contributions of amplitude and phase to the total errors are the same as those found in the previous section by LICOM3.0 and TPXO9. For comparison, we evaluated the results of other models, with the total errors for the M₂ and K₁ constituents being 8.26 cm and 2.48 cm for HYCOM (Arbic et al., 2010) and 8.91 cm and 2.52 cm for STORMTIDE (Müller et al., 2014), respectively, which are smaller than those from LICOM3.0.

Table  5.  Signal $ S $, total error $ d $, amplitude error $ {d}_{\mathrm{a}} $, and phase error $ {d}_{\mathrm{p}} $ of the eight major constituents of LICOM3.0 relative to st102 data, the percentages of captured sea surface height variance $ {V}_{\mathrm{c}\mathrm{a}\mathrm{p}\mathrm{t}} $ for each constituent of LICOM3.0, STORMTIDE (Müller et al., 2014), and HYCOM (Arbic et al., 2010), and the root square sum (RSS) values of the eight major constituents

Constituent and RSS	Error/cm					$ {V}_{\mathrm{capt}} $/%
	$ S $	$ d $	$ {d}_{\mathrm{a}} $	$ {d}_{\mathrm{p}} $		LICOM3.0	STORMTIDE	HYCOM
M2	33.22	13.85	9.56	8.12		82.62	93.9	93.8
K1	11.26	3.44	1.89	2.39		90.67	95.0	95.1
O1	7.76	6.21	3.31	4.62		35.89	83.2	89.7
S2	12.62	6.35	4.33	4.03		74.66	86.9	83.2
P1	3.62	1.61	0.73	1.26		80.22	94.7	95.2
Q1	1.62	1.06	0.66	0.73		57.19	64.7	82.1
N2	6.86	2.02	1.36	1.21		91.07	96.0	95.9
K2	3.43	1.58	1.23	0.74		78.78	89.7	76.9
RSS	39.04	17.11	11.36	10.64		80.76	92.8	92.6

 | Show Table

DownLoad: CSV

Regarding the percentage of sea surface height variance that we captured, the simulated M₂ and K₁ constituents of LICOM3.0 were found to reach 82.62% and 90.67%, respectively. Among the eight constituents, the accuracy of the simulated N₂ constituent is the best, with a variance of 91.07%. The poorest accuracy constituent was O₁, with a variance of 35.89%. The RSS of the eight constituents in LICOM3.0 is approximately 80.76% relative to the st102 data, which indicates that LICOM3.0 can generally simulate the eight major constituents well. In addition, this RSS was compared with the results of Müller et al. (2014) and Arbic et al. (2010), whose variances reached 92.8% and 92.6%, respectively. In addition, von Storch et al. (2023) increased the variance explained by adjusting the model time step and altering the coefficients of the bottom friction term and the topographic wave drag term; however, this variance was still below 80%. The tide simulation of LICOM3.0 is slightly inferior to that of Müller et al. (2014) and Arbic et al. (2010) but superior to that of von Storch et al. (2023).

Figure 4 displays the spectral analysis of the sea level data from the Yap and Kodiak tide gauge stations (Ponchaut et al., 2001) and the sea surface height data from LICOM3.0. The results from LICOM3.0 at the frequency bands of the eight major constituents are similar to the tide gauge station results; however, the amplitude at the frequency band of the M₂ constituent is significantly lower in LICOM3.0 than in the tide gauge station results. The amplitudes of the major constituents are similar in both the LICOM3.0 and tide gauge station results. Furthermore, there are some errors in low-frequency bands between LICOM3.0 and the tide gauge station results, but tides do not cause these errors and are outside the scope of this discussion. These low-frequency errors have been observed in the results of other articles, such as Jin et al. (2022).

Figure  4.  Spectral analysis of sea levels observed (obs) by Yap (9.5083°N, 138.1283°E) and Kodiak (57.7317°N, 152.511 7°W) in WOCE (Ponchaut et al., 2001) and simulated by LICOM3.0.

DownLoad: Full-Size Img PowerPoint

4. Conclusions and discussion

In this study, we incorporated the astronomical tidal potentials of eight major constituents into a global eddy-resolving OGCM, LICOM3.0. Herein, we comprehensively evaluated the tidal simulation ability of the model. The resulting sea surface height was analyzed through tidal harmonic analysis to obtain the harmonic constants. The harmonic constants were compared with the data from TPXO9 and st102 and with previous results. In addition, the sea surface height was compared with WOCE tide gauge stations through spectral analysis. The main conclusions are as follows:

(1) LICOM3.0 can effectively simulate barotropic tides. The amplitude and phase of the M₂ and K₁ constituents simulated by LICOM3.0 are approximately the same as the TPXO data in spatial distribution, and the position of the amphidromic points can also be correspondingly mapped. The simulation error is smaller in the open ocean and larger in the nearshore regions.

(2) Relative to the TPXO data, the total errors for the M₂ and K₁ constituents are 10.19 cm and 2.85 cm, respectively. Specifically, the amplitude error is more significant in the M₂ constituent, whereas the K₁ constituent exhibits a larger phase error. This situation is different from those in other models.

(3) In LICOM3.0, the simulated tidal amplitude is greater than that of TPXO, and the errors in these high-amplitude regions are relatively great.

(4) LICOM performs comparably to other models. Relative to the st102 data, the tidal error in LICOM3.0 is smaller than that of ICON-O (von Storch et al., 2023) but slightly larger than those of STORMTIDE (Müller et al., 2014) and HYCOM (Arbic et al., 2010).

In the discussion of this article, the tidal results simulated by LICOM3.0 were found to be suboptimal, indicating a need for further adjustment of the coefficients of the topographic wave drag term. Although this paper featured only an initial evaluation of the tidal simulation in the eddy-resolving version of LICOM3.0, the simulated tides were sufficiently accurate to capture the signal of internal tides, allowing for the study of some statistical characteristics. To investigate the influences of simulated tides on ocean energy distribution and circulation, we should further explore the synergistic impact of internal tides and circulation, building on existing barotropic tide simulation results.