Open-ocean tides simulated by ICON-O, version icon-2.6.6

This paper evaluates barotropic tides simulated by a newly developed multi-layer ocean general circulation, ICON-O, and assesses processes and model configurations that can impact the quality of the simulated tides. Such an investigation is crucial for applications addressing internal tides that are much more difficult to evaluate than the barotropic tides. Although not specially tuned for tides and not constrained by any observations, ICON-O is capable of producing the main features of the open-ocean barotropic tides as described by the geographical distributions of amplitude, phase, and amphidromic points. An error analysis shows, however, that the open-ocean tides simulated by ICON-O are less accurate than those simulated by two other ocean general circulation models (OGCMs), especially when not properly adjusting the time step and the parameters used in the time-stepping scheme. Based on a suite of tidal experiments, we show that an increase in horizontal resolution only improves tides in shallow waters. Relevant for using ICON-O with its telescoping grid capacity, we show that spatial inhomogeneity does not deteriorate the quality of the simulated tides. We further show that implementing a parameterization of topographic wave drag improves the quality of the simulated tides in deep ocean independent of the model configuration used, whereas the implementation of a self-attraction and loading (SAL) parameterization in a low-resolution (40 km) version of ICON-O degrades the quality of tides in shallow ocean. Finally, we show that the quality of tides simulated by ICON-O with low resolution (40 km) can be significantly improved by adjusting the time step or the parameters in the time-stepping scheme used for obtaining the model solution.

1 Introduction

The development of a new climate and Earth system model, ICON, carried out at the Max Planck Institute for Meteorology in cooperation with the German meteorological service (DWD), has been advancing quickly in the last few years. After having documented the atmospheric and oceanic components, ICON-A and ICON-O (Giorgetta et al., 2018; Korn et al., 2022), a standard application of the coupled ICON model in form of a set of CMIP6 simulations at a resolution typical of models participating in the Coupled Model Intercomparison Project (CMIP) has been assessed by Jungclaus et al. (2022). The ICON model is now used to push the frontiers of simulating the climate based on first principles by resolving the major energetic motions, including those at small scales such as the deep convection in the tropics and eddies and waves in the ocean. The activity takes the form of running storm-resolving atmosphere-only and storm-resolving coupled atmosphere–ocean simulations (Stevens et al., 2004; Hohenegger et al., 2022), a development referred to as ICON-Sapphire. Here “storm-resolving” refers to resolutions of a few kilometers (well below 10 km). This new development, however, requires further efforts. Some of them are needed to better represent potentially important processes that still elude a physical description, such as the mixing processes that occur at hectometer and finer scales. Others are needed in order to better represent phenomena resolved by kilometer-scale models. One such phenomenon is internal tides that are generated as barotropic tides flow over topographic features, such as underwater ridges and sea mounts. To ensure a realistic representation of internal tides, one needs to make sure that the model used is capable of realistically simulating the barotropic tides. This paper evaluates barotropic tides simulated by ICON-O and assesses processes and model configurations that can impact the quality of the simulated tides.

Internal tides are internal waves at the tidal frequencies. As one of the most energetic components of the ocean internal wave field, internal tides play an important role for the ocean general circulation. It is thought that the breaking of internal waves provides the power needed to maintain the oceanic overturning circulation (Munk and Wunsch, 1998; Wunsch and Ferrari, 2004; Ferrari and Wunsch, 2009). The storm-resolving models provide a new tool for studying internal tides. In fact, the collaborative research center TRR181 founded by the German Research Foundation made a proposal to advance the investigation of internal tides and their interactions with mesoscale eddies using a storm-resolving ICON-O. To further enhance the horizontal resolution in specific regions, a telescoping configuration has been developed (Korn et al., 2022). The target configuration has a fine resolution below 1 km, down to about 600 m, within a focus domain of about several thousand kilometers. With such a fine resolution, the proposal aims to simulate a larger range of vertical modes associated with internal tides.

A prerequisite for a realistic simulation of internal tides is a realistic simulation of barotropic tides. Modeling barotropic tides has a long history. As reviewed by Stammer et al. (2014), the modern tidal models, which are barotropic models tightly constrained by satellite altimeter data, produce much more accurate tides than unconstrained hydrodynamic tidal models. Typically, the differences between the observed bottom pressure recorders and the unconstrained models are about a factor of 10 larger than corresponding differences for the modern tidal models that assimilate altimetry. This large difference reflects the overall difficulty in accurately simulating tides without constraining the model with observations, as is usually the case when developing models for climate purposes. Generally, the practice of modeling tides using global ocean general circulation models (OGCMs) is much less mature than that of modeling tides using barotropic tidal models. To our knowledge, modeling tides using multi-layer OGCMs (for the purpose of studying internal tides) has only been carried out in HYCOM simulations using the HYbrid Coordinate Ocean Model (Arbic et al., 2010, 2012) and in STORMTIDE/STORMTIDE2 using the Max Planck Ocean Model (MPIOM) (Mueller et al., 2012; Li and von Storch, 2020). For the new ocean model ICON-O, only the coastal tides have been evaluated (Logemann et al., 2021). No evaluation has been performed concerning the open-ocean tides. It is also not clear whether and how various model configurations impact the quantity of the open-ocean tides simulated by ICON-O. Addressing these issues is an important step towards a storm-resolving OGCM capable of simulating both mesoscale and sub-mesoscale eddies and a large range of modes of internal tides in general and towards a successful implementation of the numerical experiments planned by TRR181 in particular.

We note that accurately simulating tides can be beneficial for a large range of applications. Apart from studying internal tides, an important application is to provide tidal corrections to various measurements so that small non-tidal signals can be identified and investigated (Knudsen and Andersen, 2002). It is conceivable that different applications require different degrees of accuracy in the simulated tides. In particular, tides for the purpose of studying overall behavior of internal tides may not need to be simulated as accurately as tides used for providing tidal corrections. Adequately addressing this difference is beyond the scope of this paper. We will only briefly discuss this issue at the end of paper.

This paper provides a first systematic evaluation of the tides simulated by ICON-O. Such an evaluation becomes possible after the implementation of the lunisolar tidal potential into ICON-O by Logemann et al. (2021). Different from the study by Logemann et al. (2021), which focuses on coastal tides simulated by the ICON-O with highly irregular meshes suitable to coastal studies, this paper focus on open-ocean tides for studying internal tides in the ocean interior.

After describing the suite of tidal experiments performed and the analysis methods used in Sect. 2, we assess the quality of the tides simulated by a standard configuration of ICON-O in Sect. 3. We then evaluate the effect of the horizontal resolution, the effect of the horizontal in-homogeneity of the grid, and the effect of the vertical coordinate on the quality of the tides in Sect. 4. Section 5 describes the effects of parameterizations of topographic wave drag and self-attraction and loading on the quality of the simulated tides. A summary and discussions are given in Sects. 6 and 7.

2 Experiments and error analysis methods

2.1 The suite of tidal experiments

The tidal experiments are performed with ICON-O, an ocean circulation model based on primitive equations. The momentum equation includes as the lunisolar tidal forcing the horizontal component of the difference between the gravitational force due to the celestial body (moon and sun) and the centrifugal force caused by the Earth's rotation around the center of mass of the Earth–celestial body system (Logemann et al., 2021). The novelty of the ICON-O numerics is the use of specific reconstructions that are required to combine scalar and velocity fields for flux calculation (Korn, 2017). These reconstructions are compatible with a discrete scalar product for the state space of ICON-O, consisting of the horizontal velocity, potential temperature, salinity, and surface elevation. This allows for writing the discrete primitive equations in a discrete weak or variational form, a necessary condition to derive discrete conservation principles. The design of ICON-O has been centered on these discrete conservation principles.

The horizontal grid of ICON is created by recursively dividing the original 20 triangles of the icosahedron via bisecting the edges, which results in a so-called R2Bn grid with n indicating the number of subdivisions used (Giorgetta et al., 2018). The horizontal resolution is estimated as the square root of the average triangular area, which may overestimate the true effective resolution (Danilov, 2022). The bottom topography is interpolated from SRTM30 (Farr et al., 2007). The vertical coordinate axis of ICON-O is given either in the z coordinate or in the z* coordinate. When using the z* coordinate, the primitive equations are solved following Adcroft et al. (2004) on a virtual grid, with the vertical coordinate being scaled in proportion with the sea surface elevation. The formulation of the z* coordinate ensures that the structure-preserving machinery developed for the z coordinate can be adopted so that tracer content, variance, and energy are conserved. Both the z and z* coordinate employ a total of 128 layers, with a vertical spacing that is smaller than 10 m (apart from the first layer in the z coordinate) in the upper 140 m, smaller than 100 m in the upper 3450 m, and identical to 200 m below 4544 m.

In all numerical experiments shown here we employ a grid-size-dependent biharmonic viscosity as horizontal friction. The vertical mixing is parameterized using the TKE (turbulent kinetic energy) scheme following Gaspar et al. (1990). The effect of mesoscale eddies is parameterized following Griffies (1998), which is switched on for the simulation at a resolution of about 40 km (R2B6, defined below) but switched off otherwise.

The suite of tidal experiments (Table 1) are designed to address the following three questions. First, how accurate are the tides simulated by a standard configuration of ICON-O? We consider ICON-O configured with a R2B6 grid with a horizontal resolution of about 40 km in the z coordinate as our standard configuration. This choice is made because ICON-O at R2B6 resolution is normally used for climate simulations (Jungclaus et al., 2022). More importantly, a resolution of 40 km generally allows an adequate representation of topographic features characterized by underwater ridges, sea mounts, passages, and straits. This could be (at least partially) the reason that the modern data-constrained tidal models can already reach a high accuracy at a resolution of about 0.5^∘ (Stammer et al., 2014).

Secondly, we address the question of whether and how different model configurations set by horizontal resolution, spatial inhomogeneity of the grid, and different vertical coordinates (z coordinate versus z* coordinate) affect the quality or accuracy of the simulated tides. We do so by performing tidal experiments in which one of the three aspects – horizontal resolution, horizontal inhomogeneity, and vertical coordinate – is modified, while keeping all other aspects of ICON-O almost unchanged.

Generally, it is expected that increasing horizontal resolution could improve the accuracy of the simulated tides. Does this improvement concern only the tides in the coastal region, or do we see also improvement for the open-ocean tides? We assess the effect of horizontal resolution by increasing the horizontal resolution from the standard R2B6 configuration (about 40 km) to R2B8 (about 10 km).

A strong spatial inhomogeneity is introduced by a telescoping grid, which can affect the quality of the simulated tides. When wobbling from a region with higher resolution to a region with lower resolution, tides can be deformed and deteriorated. We assess the effect of spatial inhomogeneity by comparing the tides simulated by ICON-O on a telescoping grid with the tides simulated by ICON-O on the standard R2B6 grid. We use a “base camp telescope” (BCT) grid shown in Fig. 1, whose focus domain lies in the South Atlantic over the Walvis Ridge and has a resolution of about 8 km. The focus domain was dictated by the observational campaign carried out for TRR181. Outside the focus domain of about several thousand kilometers, the resolution decreases gradually to the coarsest value of about 80 km. Even though it is much coarser than the telescoping grid that will be used in the final runs proposed by TRR181, the BCT grid has roughly the same degree of inhomogeneity, i.e., roughly the same ratio for the smallest grid size to the largest grid size.

In addition to horizontal resolution and spatial inhomogeneity, we also consider the effect of vertical coordinates. The experiments proposed by TRR181 will be done with ICON-O in the z* coordinate with a thin surface layer to better resolve sub-mesoscales near the surface. When using the standard z coordinate, the first layer must be sufficiently thick to “absorb” both a time-varying sea ice growth and a time-varying sea surface elevation, which could be large locally in the presence of tides. In fact, even with a surface layer of 11 m, the model becomes unstable when switching on the tides. In order to complete the simulations with tides, we had to further increase the thickness of the first layer from 11 to 14 m (indicated in the fifth column of Table 1). When switching to the z* coordinate, we employ a surface layer of 2 m, which is needed for a better representation of near-surface sub-mesoscale motions. Employing a thin surface layer leads to a change in the distribution of the 128 vertical layers. This change, only noticeable in the upper ocean, is illustrated in Fig. 2. Generally, we do not expect big changes in the quality of the simulated tides when transforming from z to z* coordinates. We instead consider the result obtained with a R2B6 ICON-O in z* coordinates as a further technical check.

Figure 1 Resolution of the telescope grid BCT in kilometers.

Figure 2 Layer thickness as a function of layer number used in the z- (red) and z*-coordinate (blue) configuration.

Finally, we also address the question of whether and to what extent parameterizations of two secondary tidal processes affect the quality of the simulated tides. One of these process is self-attraction and loading (SAL) and the other is the topographic wave drag (TWD), which describes the transfer of the tidal energy to the internal tide energy. While the SAL effects have been mostly included in the barotropic tidal models (Ray, 1998; Accad et al., 1978), an attempt is also made to include the SAL effects in an OGCM (Shihora et al., 2022).

The suite of tidal experiments used to address the above questions are listed in Table 1. They can be classified into three groups. The first one (italic) comprises spinup runs (Hertwig et al., 2021). Since we aim to simulate tides together with the oceanic circulation, and since a change in the model configuration, e.g., a change in the horizontal resolution, can make the model drift to a somewhat different state, all tidal experiments start from a state at the end of a respective spinup run. A spinup run starts from the temperature and salinity fields interpolated using the Polar Science Center Hydrographic Climatology PHC2 (Steele et al., 2001) from rest and is forced by the Ocean Model Intercomparison Project (OMIP) forcing (Roeske, 2006) without tides. After 100 years, the climate has spun up reasonably. Figure 3 shows, for example, that the Atlantic meridional overturning circulation (AMOC) at 26^∘ N is by and large stabilized in all spinup runs in the last few decades.

Figure 3 Time series of Atlantic meridional overturning circulation defined as the maximum of zonally averaged streamfunction at 26^∘ N and derived from the four spinup experiments listed in Table 1.

Table 1 Names and further specifications of the suite of simulations considered in this paper, consisting of spinup simulations (italic); simulations used to quantify the effect of resolution, spatial inhomogeneity, and vertical coordinates (bold); and simulations used to assess the effects of self-attraction and loading (SAL) and topographic wave drag (TWD) on the quality of tides in different configurations (roman). R2B6 denotes the use of a resolution of about 40 km, R2B8 denotes the use of a resolution of about 10 km, and BCT stands for base camp telescope. L128 means 128 unequally spaced vertical levels. PHC3 refers to Polar Science Center Hydrographic Climatology (Steele et al., 2001). Only the height of the first levels varies between the runs and is given together with the vertical coordinate (z or z*).

The second (bold) and third (roman) groups in Table 1 are the actual tide experiments (Hertwig et al., 2021, 2022). The experiments in bold are used to assess the effects of horizontal resolution, vertical coordinates, and spatial inhomogeneity on the simulated tides. The other experiments (roman) are used to assess the effect of SAL and TWD in different model configurations. All tidal experiments start from year 97 of the respective spinup run and are forced by both the OMIP forcing and the tidal potential for 3 years. For the last of the 3 years, we output the sea surface elevation at a hourly frequency (Hertwig et al., 2021). The analysis will be carried out for the 1-year hourly output of each experiment.

2.2 Error analysis

To evaluate the quality of the simulated tides, we mainly use the tidal model product TPXO9, the most recent version of a global barotropic model of ocean tides, which best fits (in a least-squares sense) the Laplace tidal equations and altimetry data and will hereafter be considered our proxy observation. More details about TPXO can be found in Egbert et al. (1994) and Egbert and Erofeeva (2002). The tides are harmonically analyzed by the least-squares-fitting method of Foreman et al. (2009). From this harmonic analysis the amplitudes and phases of the eight major diurnal and semi-diurnal tidal constituents (M₂, S₂, N₂, K₂, K₁, O₁, P₁, and Q₁) are derived. The amplitudes and phases of the simulated tides are compared with those obtained from our proxy observation, TPXO9. For the comparison, the model grid is interpolated into the TPXO9 grid.

Following Arbic et al. (2004) and Stammer et al. (2014), we quantify the fidelity of the simulated kth tidal constituent at a grid point in terms of the total error defined as the squared difference $d_k^{\mathrm{2}}$ between the simulated and observed tides at that grid point averaged over a tidal period:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{d}_k^{\mathrm{2}}& =<\left(A_{\mathrm{o}}\cos \right(\mathit{\omega }t-{\mathit{\varphi }}_{\mathrm{o}})-A_{\mathrm{m}}\cos (\mathit{\omega }t-{\mathit{\varphi }}_{\mathrm{m}}){)}^{\mathrm{2}}>\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(A_{\mathrm{o}}^{\mathrm{2}}+A_{\mathrm{m}}^{\mathrm{2}})-A_{\mathrm{o}}{A}_{\mathrm{m}}\cos ({\mathit{\varphi }}_{\mathrm{o}}-{\mathit{\varphi }}_{\mathrm{m}})\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(A_{\mathrm{o}}-A_{\mathrm{m}}{)}^{\mathrm{2}}+A_{\mathrm{o}}{A}_{\mathrm{m}}(\mathrm{1}-\cos ({\mathit{\varphi }}_{\mathrm{o}}-{\mathit{\varphi }}_{\mathrm{m}})),\end{array}\end{array}$$

where A_o and A_m are the amplitudes; ϕ_o and ϕ_m are the phases of the kth tidal constituent identified from TPXO9 (denoted by the subscript o for “observations”) and from ICON-O simulation (denoted by m for “model”), respectively; and $<\cdot >$ denotes an average over the respective tidal period.

It is obvious from Eq. (1) that the total error as measured by $d_k^{\mathrm{2}}$ can be decomposed into two parts:

$$\begin{array}{}\text{(2)}& d_k^{\mathrm{2}}=d_{k,\mathrm{AM}}^{\mathrm{2}}+d_{k,\mathrm{PH}}^{\mathrm{2}},\end{array}$$

with

$$\begin{array}{}\text{(3)}& d_{k,\mathrm{AM}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(A_{\mathrm{o}}-A_{\mathrm{m}}{)}^{\mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(4)}& d_{k,\mathrm{PH}}^{\mathrm{2}}=A_{\mathrm{m}}{A}_{\mathrm{o}}(\mathrm{1}-\cos ({\mathit{\varphi }}_{\mathrm{o}}-{\mathit{\varphi }}_{\mathrm{m}}\left)\right).\end{array}$$

Here, $d_{d,\mathrm{AM}}^{\mathrm{2}}$ measures the errors arising solely from model's inability to correctly simulate the amplitude of a tide and can be obtained by letting ϕ_m=ϕ_o in Eq. (1), and $d_{d,\mathrm{PH}}^{\mathrm{2}}$ measures the errors arising solely from model's inability to correctly simulate the phase of a tide and can be obtained by letting A_m=A_o in Eq. (1).

The errors described by the local measures $d_k^{\mathrm{2}}$, $d_{k,\mathrm{AM}}^{\mathrm{2}}$, and $d_{k,\mathrm{PH}}^{\mathrm{2}}$ can be further summarized by averaging over an area (or over a set of collected grid points). The results are denoted by

$$\begin{array}{}\text{(5)}& D_k^{\mathrm{2}}=\stackrel{\mathrm{‾}}{\left(d_k^{\mathrm{2}}\right)},D_{k,\mathrm{AM}}^{\mathrm{2}}=\stackrel{\mathrm{‾}}{\left(d_{k,\mathrm{AM}}^{\mathrm{2}}\right)},D_{k,\mathrm{PH}}^{\mathrm{2}}=\stackrel{\mathrm{‾}}{\left(d_{k,\mathrm{PH}}^{\mathrm{2}}\right)},\end{array}$$

where $\stackrel{\mathrm{‾}}{(\cdot )}$ indicates an average over grid points.

To quantify the overall performance of the M=8 tidal constituents listed above, we follow Arbic et al. (2004) and consider two additional measures. One measures the combined error for all eight tidal constituents, referred to as the RSS (root square sum) of the error,

$$\begin{array}{}\text{(6)}& {\mathrm{RSS}}^{\mathrm{2}}=\sum _{k=\mathrm{1}}^M{D}_k^{\mathrm{2}}.\end{array}$$

The other is the skill score

$$\begin{array}{}\text{(7)}& E=\mathrm{100}\mathit{\%}\times (\mathrm{1}-{\displaystyle \frac{{\mathrm{RSS}}^{\mathrm{2}}}{\sum _{k=\mathrm{1}}^M{S}_k^{\mathrm{2}}}}),\end{array}$$

defined using the total error measured by RSS² relative to the sum over the observed signals S_k

$$\begin{array}{}\text{(8)}& S_k^{\mathrm{2}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{A_{\mathrm{o}}^{\mathrm{2}}}}{\mathrm{2}}},\end{array}$$

where A_o denotes the amplitude of the kth tidal constituent obtained from TPXO. E equals 100 % if all eight tidal constituents are perfectly simulated by the model. Both RSS defined in Eq. (6) and the skill score E defined in Eq. (7) can also be calculated for errors with respect to amplitude or phase by replacing d_k in RSS with d_k,AM and d_k,PH, respectively (i.e., RSS_AM, RSS_PH, E_AM, and E_PH).

3 Tides simulated by ICON-O in the standard R2B6 configuration

Generally, there is an overall good agreement between tides simulated by ICON-O in its standard R2B6 configuration and the tides in TPXO9. For the M₂ tide, which is the strongest tidal constituent in the ocean, the geographical distributions of amplitude (colors) and phase (black lines), including the amphidromic points, obtained from ICON-O (top panel of Fig. 4) are by and large comparable to those obtained from TPXO9 (bottom panel of Fig. 4). Similar comparability between R2B6 and TPXO9 is found for the K₁ tides – the strongest diurnal tide (Fig. 5). The most obvious difference is the too strong amplitude in R2B6, especially for the M₂ tide in the Atlantic. Some noticeable phase differences are found for the K₁ tide in the Southern Ocean. Another difference concerns the tiny “wiggles” of the lines of constant phase, which are found in R2B6 (especially for the K₁ tide in the Pacific in Fig. 5) but are clearly absent in TPXO9. Arbic et al. (2012) found similar wiggles in the phase lines of M₂ tide simulated by HYCOM and attributed them to the simulated M₂ internal tides. The wavelength of the K₁ internal tide, about 200–250 km equatorward of the critical latitudes (Li et al., 2017), is longer than that of M₂ internal tide, which is about 150 km at maximal extent (Li et al., 2015). Since waves with longer wavelengths can be simulated more easily, the wiggles are more pronounced for the phase lines of the K₁ tide than those of the M₂ tide.

Figure 4 Amplitude (in cm) and phase of the M₂ tide simulated by ICON-O in the standard R2B6 configuration (a) and derived from TPXO9 (b).

Figure 5 The same as Fig. 4 but for the K₁ tide.

To further quantify the errors against TPXO9, we show maps of the differences d_k, d_k,AM, and d_k,PH for k=M₂ in R2B6 in the top panels of Figs. 6–8. The largest errors in M₂ tide are found off the coast of Iberian Peninsula, in the Labrador Sea, and west of central America (top panel in Fig. 6), where the M₂ tide has its largest amplitude in TPXO9 (dark shading in the bottom panel of Fig. 4). It seems that there is a correlation between the magnitude of the errors and the strength of the signal. The top panels in Figs. 7 and 8 show that these errors arise from the inaccuracy in simulating both the amplitude and the phase of the tide, with the contribution from amplitude error being larger than that from phase error in the North Atlantic but vice versa for the errors west of central America.

Large errors are also found in the Southern Ocean along the cost of West Antarctica west of the Drake Passage, where M₂ tide has small amplitudes in TPXO9 (color shadings in the bottom panel of Fig. 4). The study by Pal et al. (2023) suggests that these errors can be induced by ignoring Antarctic ice shelf cavities in ICON-O. The top panels in Figs. 7 and 8 further show that these errors result mainly from amplitude errors.

Table 2 summarizes the signals; the errors measured by D_k, D_k,AM, and D_k,PH for each constituent; and the overall skill E averaged for all constituents (last column). To focus on the open-ocean tides, all numbers are obtained by averaging over grid points deeper than 1000 m. We see that the signals of the R2B6 tides (second row of numbers) are stronger than those of the TPXO9 tides (first row). For the four semi-diurnal tides, the largest value of D_k, which amounts to 14.25 cm, is found for the M₂ tide (third row). This error is related to a strong overestimation of amplitude: while the TPXO9 M₂ signal is about 24.5 cm, the R2B6 M₂ signal is 31.67 cm. This overestimation in amplitude is reflected in the larger value of $D_{M_{\mathrm{2}},\mathrm{AM}}$ (fourth row) than the value of $D_{M_{\mathrm{2}},\mathrm{PH}}$ (fifth row). For the four diurnal tides, the largest value of D_k of 5.19 cm is found for the O₁ tide, the second-strongest diurnal tide. As for the M₂ tide, the amplitude seems to contribute more to the total error than the phase.

For the M₂ tide, the averaged error D_k of 14.25 cm is clearly larger than those produced by all the unconstrained barotropic and baroclinic hydrodynamic models discussed in Stammer et al. (2014). Their Table 12 shows averaged errors mostly smaller than 10 cm. When averaged over all eight constituents, the overall skill E is only about 63 %. Ignoring the errors in amplitude or in phase increases the skill to about 81 %.

The large averaged errors and the low skill suggest that the tides in ICON-O are less accurate than those found in the other unconstrained models. To further confirm this impression, the tides simulated by the different models must be evaluated against the same observational evidence. For the two OGCMs considered in Stammer et al. (2014), HYCOM with a horizontal resolution of $\mathrm{1}/\mathrm{12.5}{}^{\circ }$ and MPIOM with a horizontal resolution of $\mathrm{1}/\mathrm{10}{}^{\circ }$, error analyses against the 102 pelagic tide measurements (Shum et al., 1997) are available. We hence repeated the analysis for R2B6 but now against the pelagic data. We obtained for the R2B6 tides the E value of 77.4 %, which is clearly lower that 89.02 % in STORMTIDE2 (Li and von Storch, 2020), 92.8 % in STORMTIDE (Mueller et al., 2012), and 92.6 % in HYCOM (Arbic et al., 2010). Although the value of 77.4 % averaged over pelagic stations is somewhat larger than the 63.38 % averaged over grid points deeper than 1000 m (last column in Table 2), the result clearly shows that the ICON-O in its standard R2B6 configuration is less skillful in simulating tides than the $\mathrm{1}/\mathrm{12.5}{}^{\circ }$ HYCOM and $\mathrm{1}/\mathrm{10}{}^{\circ }$ MPIOM. It is tempting to attribute the low skill to the lower horizontal resolution of about 40 km relative to the resolution of about 10 and 8 km in STORMTIDE (and STORMTIDE2) and the HYCOM simulation, respectively. As will be discussed in the next section, increasing horizontal resolution by using a R2B8 grid improves some aspects but not the skill in the deep ocean.

Table 2 Values of the signals D_k, D_k,AM, and D_k,PH (first eight columns of numbers) in centimeters and the skills E, E_AM, and E_PH (last column of numbers) as derived from the ICON-O in the standard R2B6 configuration. All values are averaged over grid points with water deeper than 1000 m.

Figure 6 Difference d_k resulting from the total error in the M₂ tide in ICON-O in the R2B6 (a), R2B8 (b), and BCT (c) configurations; d_k is defined in Eq. (1) and evaluated against TPXO9 (values are given in centimeters).

Figure 7 The same as Fig. 4 but for the difference d_k,AM resulting from errors in amplitude only.

Figure 8 The same as Fig. 4 but for the difference d_k,PH resulting from errors in phase only.

4 Effects of the horizontal resolution, spatial inhomogeneity, and vertical coordinates

When increasing the horizontal resolution from about 40 km (R2B6) to about 10 km (R2B8), the error d_k is not reduced for all tidal constituents at all grid points. For the M₂ tide, the middle panel in Fig. 6 shows that in many places the error seems to be enhanced rather than reduced. When decomposing $d_k^{\mathrm{2}}$ into $d_{k,\mathrm{AM}}^{\mathrm{2}}$ and $d_{k,\mathrm{PH}}^{\mathrm{2}}$, the middle panels of Figs. 7 and 8 show that increasing horizontal resolution leads to a general reduction in the error in amplitude but worsens the error in phase, especially in the Pacific.

Figure 9 Effects of increasing horizontal resolution as quantified by the globally averaged errors $D_{M_{\mathrm{2}}}$, $D_{M_{\mathrm{2}},AM}$, and D_M2,PH (a, b) for the M₂ tide and by the skill scores E, E_AM, and E_PH averaged over eight constituents (c, d) simulated by ICON-O with a resolution of about 40 km (blue) and 10 km (red) and with the BCT grid (orange). The bars in (a) and (c) are obtained by averaging over grid points deeper than 1000 m, whereas the bars in (b) and (d) are obtained by averaging over grid points shallower than 1000 m.

That increasing horizontal resolution does not improve the accuracy of the simulated M₂ tide in deep oceans is further quantified in Fig. 9, which compares the errors (top) and the skills (bottom) averaged over grid points shallower than 1000 m (right) with those averaged over grid points deeper than 1000 m (left) for the M₂ tide. Consider first the situation in the deep ocean. We see that the errors in simulating the amplitude of the M₂ tide ($D_{M_{\mathrm{2}},\mathrm{AM}}$, Fig. 9a) is reduced in R2B8 (red) relative to that in R2B6 (blue) as already suggested by Fig. 7. This slight improvement is however outweighed by a larger phase error $D_{M_{\mathrm{2}},\mathrm{PH}}$, making D_M2 larger in the R2B8 run than in the R2B6 run. When considering all eight tidal constituents (Fig. 9c and d), the overall skill E is lower in R2B8 (red, Fig. 9c) than in R2B6 (blue). This somewhat lower skill in R2B8 is due to the lower skill in simulating the phase as measured by E_PH, while the skill in simulating the amplitude as measured by E_AM is slightly enhanced.

A different picture is found for tides in shallow seas (Fig. 9b and d). For the M₂ tide, all three error measures $D_{M_{\mathrm{2}}}$, $D_{M_{\mathrm{2}},\mathrm{AM}}$, and $D_{M_{\mathrm{2}},\mathrm{PH}}$ are smaller in R2B8 (red bars in Fig. 9b) than in R2B6 (blue bars in Fig. 9b). When considering all constituents, all three skill scores E, E_AM, and E_PH are higher in R2B8 (red bars in the Fig. 9d) than in R2B6 (blue bars in Fig. 9d). We hence conclude that the major effect of increasing horizontal resolution is to improve the quality of the simulated tides in shallow seas, without having significant impact on the tides in the deep ocean.

Using a spatially inhomogeneous telescoping grid does not seriously deteriorate the quality of the simulated tides. The spatial distributions of $d_{M_{\mathrm{2}}}$, $d_{M_{\mathrm{2}},\mathrm{AM}}$, and $d_{M_{\mathrm{2}},\mathrm{PH}}$ in BCT (bottom panels in Figs. 6–8) are comparable to those in R2B6 (top panels). Figure 9 further shows that the error related to the M₂ tide, $D_{M_{\mathrm{2}}}$, is slightly reduced in BCT (orange) in deep ocean and that the skill E is slightly enhanced relative to the respective values in R2B6 (blue) in the deep ocean. In shallow seas, the quality of the tides in BCT is also improved rather than deteriorated relative to R2B6, albeit to a less extent than in R2B8. In terms of skills, the values of E, E_AM, and E_PH are higher than those in R2B6 but lower than those in R2B8. The result indicates that a spatially inhomogeneous telescope grid does not strongly deteriorate the quality of the simulated tides.

Figure 10 The overall skill score E averaged over eight constituents in six different areas, with area 1 containing grid boxes with the smallest sizes (of less than 8 km) and area 6 containing grid boxes with the largest sizes (of larger than 40 km).

We further examine the dependence of the skill in BCT on the inhomogeneous grid size by calculating E in BCT averaged over grid points in six selected areas classified according to the grid sizes in the areas. The grid size increases from smaller than 8 km in area 1 to grid size larger than 40 km in area 6. The areas are indicated by the black lines in the bottom panel of Fig. 6. We calculate the skill E for each area and compare the result with the skill E obtained from R2B6 in the respective area (Fig. 10). For areas 1–3, which have grid sizes clearly smaller than the R2B6 grid size, the BCT tides are more accurate than the R2B6 tides. The picture is less clear in areas 4 and 5, where BCT has a grid size close to the R2B6 resolution. For some unknown reason, the BCT skill is higher than the R2B6 skill in area 6, where the BCT grid sizes are larger than the R2B6 resolution. Figure 10 provides some additional confidence that the BCT grid does not seriously deteriorate the quality of the simulated tides.

We now consider the impact arising from switching from the z coordinate to the z* coordinate in the R2B6 configuration. We find hardly any change when switching from the z coordinate with a first-layer thickness of 14 m to the z* coordinate with a first-layer thickness of 2 m. The overall skill as measured by E is 63.13 %, almost identical to the 63.38 % obtained using the z coordinate.

5 Effect of topographic wave drag and SAL

Previous studies have shown that improving the accuracy of tides simulated by a barotropic tidal model requires the consideration of both the self-attraction and loading (SAL) effect (Hendershott, 1972) and the energy loss from barotropic tides to internal tides over rough topography (Jayne and St. Laurent, 2001; Arbic et al., 2004). The latter is often represented by a topographic wave drag (TWD). It is, however, not completely clear whether and how the parameterizations of SAL and TWD used in a barotropic tidal model that simulates only tides should be modified to get an accurate simulation of barotropic tides in an OGCM that simulates both tides and non-tidal motions. With respect to the TWD, a general circulation model (GCM) with a horizontal resolution of about 10 km can simulate the conversion of barotropic tides to low-mode internal tides (Arbic et al., 2012; Li et al., 2015). With such a GCM, part of the energy loss, namely the part that is transformed to the energy of the low-mode internal tides, is represented by the model and does not need to be parameterized separately. This seems to suggest that we still need to parameterize the energy loss to high-mode internal tides. When using the 102 pelagic tidal gauges as reference, the skill E of the simulated tides is 92.8 % for a MPIOM simulation that does not include a TWD (Mueller et al., 2012) and 93 % for an HIM (Hallberg isopycnal model) simulation that includes a TWD (Arbic et al., 2004). These results do not suggest that the inclusion of the energy loss to high modes plays a decisive role in improving the accuracy of the simulated tide.

In Sect. 5.1 and 5.2 we describe the two parameterizations implemented in ICON-O: the topographic wave drag (TWD) that parameterizes the energy loss over rough topography from barotropic tides to internal tides and the parameterization of SAL in form of a simple scalar approximation by Ray (1998). The impacts of these parameterizations on the quality of the simulated tides are quantified in Sect. 5.3.

5.1 Topographic wave drag

The energy conversion from barotropic tides to internal tides, which causes energy loss of barotropic tides, can be parameterized as an additional bottom drag. In addition to the standard quadratic bottom drag, $c_{\mathrm{B}}\left|{\mathit{u}}_{\mathrm{B}}\right|{\mathit{u}}_{\mathrm{B}}$, with c_B being a bottom drag coefficient and u_B the bottom velocity, we add the linear topographic wave drag:

$$\begin{array}{}\text{(9)}& {\mathit{D}}_{\mathrm{TWD}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\kappa }{h}^{\mathrm{2}}{N}_{\mathrm{B}}{\mathit{u}}_{\mathrm{B}}.\end{array}$$

The form of this drag is motivated by the scaling relation of the energy conversion from barotropic tides to internal tides (Jayne and St. Laurent, 2001). Here, h² represents the bottom roughness, which we compute as the topography variance within a 100 km radius around each grid point, and N_B is the bottom stratification.

Zarroug et al. (2010) argue that using the stratification averaged over some vertical distance above the ocean floor instead of using the bottom stratification N_B provides a better estimate of the internal tide generation and therefore also a more realistic topographic wave drag. In this scenario, the respective vertical distance should be linked to the vertical scale of the locally dominant internal tide mode. However, in this first attempt to model tides in ICON-O, we use the most widely adopted version based on the bottom stratification only.

κ represents a typical horizontal wavelength of the generated internal tides and is in reality dependent on the dominant topographic length scale. Within the framework of the TWD parameterization, κ is used as a geographically homogeneous tuning parameter that allows us to minimize the error in the modeled barotropic tide. It turns out that in ICON-O the best result can be obtained when choosing κ=50 km. This is in the range of what we expect as typical wavelength for internal tides and also what was used previously (Jayne and St. Laurent, 2001; Exarchou et al., 2012).

u_B in Eq. (9) stands for the bottom flow. Ideally, D_TWD should remove only energy from the tidal flow and not from the non-tidal bottom flow, implying that u_B should contain only the tidal bottom flow. For this purpose, Arbic et al. (2010) apply a 24 h running-mean filter to the bottom flow at each time step that approximately separates the tidal from the non-tidal flow. They then apply the drag only to the tidal flow.

The importance of such a data-intensive online filtering can be estimated by comparing the energy that is removed by the TWD when first using as u_B the full bottom flow (including tidal and non-tidal flow) and second using only tidal flow (obtained from the reconstruction based on the harmonic analysis). According to this analysis, which we performed for the R2B6 simulation, the non-tidal bottom flow accounts for only about 5 % of the energy removed by the TWD, indicating that the tidal flow is by far the dominant flow in the abyssal oceans. We conclude that using the full bottom flow as u_B in TWD parameterization will not lead to an unreasonably large error. This allows us to not apply the online time filtering and thereby keep the parameterization in its simplest form.

5.2 Self-attraction and loading (SAL)

The parameterization of SAL is supposed to cover the effects of deformation of the ocean seafloor due to the weight of the water column, the associated mass redistribution and the according changes in the gravitational field, and the gravitational attraction of the waterbody on itself. Generally, a more complete description of SAL should utilize a decomposition of mass anomalies into their spherical harmonic constituents, as was done in Gordeev et al. (1977), Brus et al. (2012), Barton et al. (2022), and Shihora et al. (2022). As a first attempt, however, here we use the scalar approximation of the SAL effect (Ray, 1998), which is less accurate than the consideration based on spherical harmonics.

Generally, SAL can arise from all barotropic motions, including non-tidal ones, especially those at high frequencies. For low-frequency motions, it is quite safe to assume that there is a quasi-instantaneous adaptation of the seawater to the time-variable external gravity field of the Earth (Shihora et al., 2022). This timescale dependence of the effect of SAL can lead to complications when implementing SAL in an OGCM that simulates both non-tidal and tidal flows. Being expressed as an additional body force, a SAL parameterization impacts simulated motions on all timescales. One cannot expect the same parameterization that works well for a barotropic tidal model to also work well for an OGCM. In the study by Shihora et al. (2022), implementing a SAL parameterization into an OGCM does not improve the accuracy of the simulated tides to the same degree as it does when implementing the same SAL parameterization into a barotropic tidal model, suggesting the complexity of the problem and the value of first concentrating on the simple scalar approach.

Following the scalar approach (Ray, 1998), we add to the astronomical tide potential an additional SAL potential that is proportional to tide-induced sea surface height (SSH) η_tidal,

$$\begin{array}{}\text{(10)}& {\mathit{\eta }}_{\mathrm{SAL}}=\mathit{\beta }{\mathit{\eta }}_{\mathrm{tidal}},\end{array}$$

with β=0.06. While the model SSH η is equivalent to η_tidal in pure tidal models (e.g., Jayne and St. Laurent, 2001), this is no longer true in a three-dimensional OGCM that has a heterogeneous stratification and non-tidal flow components (Arbic et al., 2010). To keep things simple, we consider an approximation of the tide-induced SSH η_tidal by correcting the full model SSH η for the two largest non-tidal SSH contributions:

$$\begin{array}{}\text{(11)}& {\mathit{\eta }}_{\mathrm{tidal}}\approx \mathit{\eta }-{\mathit{\eta }}_{\mathrm{tm}}-{\mathit{\eta }}_{\mathrm{steric}},\end{array}$$

where η_tm is the contribution due to the large-scale, time-mean circulation and η_steric is the so-called steric contribution. We take η_tm as the time-mean model SSH of the last 2 years of the spinup runs. Removing η_tm is consistent with the expectation that the SAL effect does not arise from low-frequency motions. The steric contribution is computed as

$$\begin{array}{}\text{(12)}& {\mathit{\eta }}_{\mathrm{steric}}=\underset{\mathrm{0}}{\overset{H}{\int }}{\displaystyle \frac{<\mathit{\rho }>-\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{d}z,\end{array}$$

where H is the ocean depth, $<\mathit{\rho }>$ is the time-mean of density ρ, and ρ₀=1025 kg m⁻³ is a reference density. Correcting for the steric SSH contribution efficiently filters out the SSH imprint of time-varying baroclinic motions such as those related to mesoscale eddies. In this way, we can make sure that our additional SAL potential does not introduce too many spurious effects due to the time-mean circulation and the non-tidal baroclinic motion in ICON-O.

5.3 Results

Further tidal simulations (indicated by the roman letters in Table 1) are carried out, with either one or both parameterizations switched on. We consider the three configurations described in Sect. 3, namely the standard R2B6 configuration, the higher-resolution configuration with the R2B8 grid, and the configuration with the non-homogenous BCT grid. These simulations are then evaluated against the TPXO9 by calculating the skill defined in Eq. (7); both E values related to the total errors; and E_AM and E_PH values related to the amplitude and phase errors, respectively. The results are shown in Fig. 11.

Figure 11 Effects of parameterizations SAL and TWD as quantified by the skill scores E (a, b), E_AM (c, d), and E_PH (e, f) for ICON-O in R2B6, R2B8, and BCT configurations. Blue bars indicate the case without employing SAL and TWD, red bars indicate the case with SAL only, orange bars indicate the case with TWD only, and magenta bars indicate the case with both SAL and TWD. The values are obtained by averaging over eight constituents and grid points deeper (a, c, e) and shallower (b, d, f) than 1000 m.

We found that the two parameterizations (Fig. 11) have different effects on tides in deep and shallow oceans. When considering the tides in deep oceans (left column of Fig. 11) with respect to the overall skill E (top left panel), the largest difference within each group of bars is found for the R2B8 configuration (middle), suggesting that the two parameterizations have a larger impact on the simulated tides in R2B8 than in R2B6 and in the BCT configuration. For R2B8, E increases from barely 60 % (blue) without SAL and TWD to about 71 % when including SAL (red) and to about 74 % when including both parameterizations (magenta). For R2B6, E increases from 63 % (blue) without SAL and TWD to 69.5 % (orange) when employing TWD, while employing SAL (red) slightly degrades the simulated tides. Both SAL and TWD together (magenta) still give an improvement for R2B6. For BCT, the skill is only slightly enhanced when employing both parameterizations (individually and together).

The improvement achieved with the R2B8 configuration arises mainly from including the SAL parameterization. As suggested by E_AM and E_PH for R2B8, the SAL-induced improvement is related to a reduction in both amplitude and phase errors, leading to larger values of E_PH and E_AM, with the increase in E_PH being more prominent than that in E_AM. The result is consistent with the previous consideration that the SAL can affect both the amplitudes and the phases (Ray, 1998).

With the R2B6 configuration, the improvement achieved arises mainly from including the TWD parameterization, which damps the simulated tides. This alleviates the overestimation of the tidal signals, leading to a higher skill E_AM with TWD than without TWD. When using the BCT configuration, the SAL parameterization slightly reduces the phase error and thus slightly increases the skill E_PH. The TWD parameterization slightly reduces the amplitude error and thus slightly increases the skill E_AM.

The situation is quite different in shallow oceans (right column of Fig. 11). In these instances, the SAL parameterization does not work properly in a low-relation configuration. We find a strong reduction in the skill E_AM and E_PH in the R2B6 and BCT configurations. Recall that relative to R2B8, BCT has a coarser resolution over a major part of the world ocean. The total skill E is even negative in R2B6. The TWD parameterization, when implemented in low-resolution configuration, also degrades the skills, although the degradation is less severe compared to that induced by SAL. Such strong degradation is not found in the R2B8 configuration. The SAL parameterization in R2B8 is able to reduce the phase errors, thereby enhancing E_PH in both deep and shallow oceans.

6 Discussions

This paper suggests that for ICON-O with 128 vertical levels, the quality of open-ocean tides cannot be satisfactorily improved by including parameterizations of SAL and TWD. Generally, when implementing parameterizations of SAL and TWD developed for a single-layer model (Egbert et al., 2004; Gordeev et al., 1977; Pal et al., 2023) in a more realistic baroclinic OGCM, some additional procedures are required. For TWD, which is a drag acting on tidal velocity only, we need a procedure that separates tidal velocity from the full velocity. For SAL, we need a procedure that separates sea level height variations arising from high-frequency barotropic motions from those arising from the large-scale time-mean circulation and low-frequency baroclinic features. Arbic et al. (2010) addressed this problem. Shihora et al. (2022) showed that implementing a parameterization of SAL in a GCM is less successful compared to the case for a barotropic tidal model. Here we try out the simplest approach. For TWD, we compare the energy removed by TWD equipped with full bottom velocity with the energy removed by TWD equipped with tidal velocity only and conclude that the error arising by using full (instead of tidal) velocity is negligible. For SAL, we use only the sea level height obtained by removing the time-mean sea level height and a steric contribution computed as a vertical integral of normalized density anomaly. Further investigation is needed to figure out whether the quality of tides simulated by ICON-O can be further improved by optimizing the separation procedures.

Apart from the SAL and TWD parameterizations, the quality of tides can also be affected by other parameters in the model. We hence performed a suite of additional experiments to address the sensitivity of the quality of the simulated tides to some of these parameters. Given that the tidal signals in our standard R2B6 configuration are too strong (first two rows in Table 2), we paid special attention to the following parameters that have the potential to affect the strength of barotropic motions: the strength of the bottom drag, the biharmonic viscosity coefficient, and the time step. To this end, we point out that the time-stepping scheme in ICON-O (Korn, 2017) has a damping effect, which strengthens with increasing time step for fixed values of parameters (i.e., for fixed values of β and γ in Eqs. 33 and 34 in Korn, 2017). To exclude the time step sensitivity, in all experiments summarized in Table 1 we used the same values of β and γ and the same time step of 300 s, which is required for a stable R2B8 simulation.

Figure 12 Sensitivity of skill score E to bottom drag, biharmonic viscosity, and time step obtained from ICON-O in R2B6 configuration. Shown are nine simulations performed with R2B6 ICON-O on the new high-performance computer Levante at DKRZ, which has recently replaced Mistral. One corresponds to the R2B6 run described before (blue); two correspond with the strength of bottom drag being doubled (green) and halved (orange), respectively; two correspond with biharmonic viscosity being increased by 14 % (red) and by 40 % (purple), respectively; four correspond with the time step being reduced from 300 to 180 s (brown), increased from 300 s to 600 (pink), increased from to 1200 s (grey), and increased from to 1800 s (yellow), respectively.

Figure 12 shows the sensitivity of the skill E to bottom drag, biharmonic viscosity, and time step obtained from this new suite of experiments. We found that doubling and halving the bottom drag coefficient both have only a small effect on the quality of the simulated tides. The same is true for increasing the biharmonic viscosity. In contrast, we find larger changes in skill E when we alter the model time step. With a time step of 1800 s, which would be normally used in R2B6 configuration, the overall skill E increases from about 65 % (blue bullet) to almost 80 % (yellow bullet). When considering only the errors in amplitudes, the skill E_AM is found to be well above 90 % (yellow cross). Similar skill values are obtained with a time step of 300 s when the values of β and γ are properly adjusted (not shown).

7 Summary

This paper assesses the ability of a newly developed OGCM, ICON-O, to realistically simulate open-ocean tides. In the standard R2B6 configuration with a horizontal resolution of 40 km, ICON-O is able to simulate the main features of the tides, as found in the geographical distributions of amplitude, phase, and amphidromic points obtained from TPXO9. When keeping the time step and the parameters in the time-stepping scheme the same as in R2B8 with finer resolution, the overall skill E (averaged over eight constituents) is about 63 % against TPXO9 and about 77 % against pelagic station data. The skill E increases to about 80 % when adjusting the time step or the parameters in the time-stepping scheme. Using the R2B8 configuration, where the horizontal resolution is increased to 10 km, leads to an overall improvement of tides in coastal regions but not in the deep oceans. Using a spatially highly non-homogenous grid does not reduce the quality of the simulated tides. The transformation from ICON-O based on z coordinates to ICON-O based on z* coordinates is successful. The latter allows simulating tides with the same accuracy without employing a thick surface layer.

The effect of SAL and TWD parameterizations on the simulated tides depends on the configurations considered. The SAL parameterization, when implemented in the R2B8 resolution, noticeably reduces the phase error and enhances the accuracy of the simulated tides in both the deep and shallow oceans. The situation is different when implementing SAL in a lower-resolution configuration. In R2B6, no improvement can be achieved for tides in deep oceans, while strong degradation is found for tides in shallow oceans. In BCT, there is a slight improvement of the skill in deep oceans but a clear reduction in skill in shallow oceans. The TWD parameterization damps the overestimated tidal amplitudes, reducing the amplitude error and enhancing the accuracy of the simulated deep-ocean tides in all three configurations (R2B6, R2B8, and BCT). In shallow oceans, TWD has little effect in the R2B8 configuration but degrades the tides in both the R2B6 and BCT configurations. The result suggests that when interested only in deep-ocean tides, the TWD parameterization should be included independent of the model configuration. It is not meaningful to study shallow-ocean tides using a low-resolution model with the SAL parameterization included.

The error and skill analyses were carried out for grid points deeper than the threshold d_* of 1000 m. This value of d_* is chosen quite arbitrarily. Varying d_* does not change the overall picture. In particular, when d_* is slightly reduced, the strength of the observed tidal signal increases and the root sum square of errors also increases, leaving the skill E as defined in Eq. (7) largely unchanged.

The purpose of assessing the ICON-O's ability to realistically simulate the open-ocean tides is to ensure that ICON-O at kilometer-scale resolution is able to realistically simulate internal tides, which are much more difficult to evaluate than the barotropic tides. The work reported here presents only the first efforts toward this goal. Even though ICON-O has a somewhat lower skill when simulating barotropic tides than MPIOM and HYCOM, we believe that the barotropic tides simulated by ICON-O are accurate enough for studying some statistics of internal tides, such as the energy fluxes related to internal tides. For such studies, the phase errors are not crucial. When taking out the phase errors and properly choosing the parameters in the time-stepping scheme, the skill, as measured by E_AM, can be well above 90 %. Although encouraging, there is plenty of room for further improvements.