Global tide simulations with ICON-O: testing the model performance on highly irregular meshes

The global tide is simulated with the global ocean general circulation model ICON-O using a newly developed tidal module, which computes the full tidal potential. The simulated coastal M2 amplitudes, derived by a discrete Fourier transformation of the output sea level time series, are compared with the according values derived from satellite altimetry (TPXO-8 atlas). The experiments are repeated with four uniform and sixteen irregular triangular grids. The results show that the quality of the coastal tide simulation depends primarily on the coastal resolution and that the ocean interior can be resolved up to twenty times lower without causing considerable reductions in quality. The mesh transition zones between areas of different resolutions are formed by cell bisection and subsequent local spring optimisation tolerating a triangular cell’s maximum angle up to 84°. Numerical problems with these high-grade non-equiangular cells were not encountered. The results emphasise the numerical feasibility and potential efficiency of highly irregular computational meshes used by ICON-O.


Introduction
Numerical modelling of global ocean dynamics using a regular (structured) grid generally leads to two problems: first, the meridional convergence, i.e. the impossibility to uniformly resolve the sphere. Though rotated spherical coordinates are shifting the North Pole onto a land mass (e.g. Marsland et al. 2003) to avoid the most severe computational problems, the fact remains that the resolution varies in a way, which is not exclusively oriented along oceanographic criteria. This also points to the second problem: oceanic fluxes exist on a broad band of scales. Thereby, fluxes caused by smaller scale processes like the flow through a strait or a western boundary current can have crucial consequences up to the basin-scale. However, simulations with sufficiently high-resolution structured grids, that would resolve these small scales, are too costly to be performed effectively.
Hence, unstructured grid ocean modelling became an attractive alternative, because it offers the possibility to circumvent the pole problem and at the same time to freely vary the resolution within the model domain. In the early 1990s, finite element ocean models, being inherently affiliated with unstructured grids, like ADCIRC (Luettich Jr et al. 1992) inspired the ocean modelling community. More recent finite element model developments are FESOM (Danilov et al. 2004), SELFE (Zhang and Baptista 2008) or SLIM (White et al. 2008). However, because of their fundamental conservation property in combination with comparable low computational costs, finite volume models mostly remained the tool of choice. This explains the ground-breaking impact of the UnTRIM model (Casulli and Walters 2000) which uses finite volumes on a so-called unstructured orthogonal grid containing cells which form arbitrary convex polygons within the horizontal plane. A number of similar model developments followed, e.g. ELCIRC (Zhang et al. 2004), FVCOM (Chen et al. 2015) or SCHISM (Zhang et al. 2016), whereas all these models use polygonal cells which share each lateral surface with exactly one neighbour cell. Such a grid is called conformal (e.g. Ives 1982).
Non-conformal grids normally consist of cuboids defined and uniformly oriented in Cartesian coordinates with cells having more than six neighbours if they border a region of higher spatial resolution. This geometrically rather simple structure reduces the computational costs and facilitates the handling of temporally variable refinements, i.e. dynamic mesh adaptions (Khokhlov 1998). Therefore, non-conformal grids are popular within the engineering sciences (e.g. Iousef et al. 2017) whereas oceanographic applications were performed by Popinet and Rickard (2007) (dynamic horizontal mesh refinement), Backhaus (2008) (static, vertical refinement) and Logemann et al. (2013) (static three-dimensional refinement). The latter model (CODE) shows a cross-scale spectrum from Icelandic fjords up the North Atlantic basin scale and was successfully applied in studies of regional oceanography and regional marine biology.
However, in the field of climate science, recent global ocean model developments all use finite volumes on conformal grids. Being the oceanic part of an earth system model MPAS-Ocean (Ringler et al. 2013) uses Voronoi tessellations, i.e. basically hexagonal cells, with an Arakawa C-grid staggering, with scalar variables at the cell centre and normal components of the velocity vector at cell boundaries (Arakawa 1966). FESOM2 (Danilov et al. 2017) uses triangular cells with an Arakawa B-grid staggering, i.e. whereas the scalar variables (e.g. pressure, temperature, salinity) are defined on the triangle vertices the flow vector is located at the cell (triangle) centre. This, at first glance surprising, staggering can be traced back to the avoidance of a numerical problem, which is related to the triangular C-grid staggering. A numerical mode, which leads to checkerboard pattern in the divergence field of the horizontal flow, was first mentioned by Stuhne and Peltier (2009) and clarified theoretically by Danilov (2010) whereas Wolfram and Fringer (2013) discuss different filtering solutions.
Finally, as part of the German ICON climate and weather modelling initiative Crueger et al. 2018), the ocean model ICON-O was developed (Korn 2017), which uses a triangular C-grid. This model overcomes the triangular C-grid dilemma by using a new discretisation technique, which integrates ideas of Finite Elements, Finite Differences and Mimetic Finite Differences and provides a new way to control the divergence noise without violating the discrete conservation conditions. This paper investigates the performance of ICON-O running on different global triangular computational meshes. These meshes will be either uniform with different resolutions or locally refined with a large horizontal resolution range between 148 and 7 km, i.e. with a local resolution increase up to a factor of 21. The mesh refinement is performed by the edge-bisection (quadtree) method wherein the original triangular cell obtains three new vertices at the midpoints of its edges. These three new vertices are connected with three new edges thereby dividing the old cell into four new cells.
Under the term "performance", we understand two aspects: first, the stability, i.e. the question whether ICON-O is able to run on a certain grid and to produce physically meaningful results, and second, the model's ability to use an irregular grid in an efficient way. In this paper, we define efficiency as the quality increase of a simulation (i.e. reduction of the error metric) caused by the replacement of a uniform grid by an irregular one having the same number of grid cells.
Our numerical experiments can be seen as a continuation of the work of Stuhne and Peltier (2009) who had successfully implemented a twofold edge-bisection refinement for the Hudson Bay region and showed that their triangular C-grid model tolerated the related sharp resolution jumps from 74 to 28 to 14 km. In our study, we will apply up to four edgebisection refinements with much smaller spacing between the jumps, i.e. with much higher horizontal resolution gradients. The basic numerical problem with high resolution gradients is caused by the fact that the originally nearly equiangularly shaped triangular cells are distorted which leads to growing numerical errors and, in the case of ICON-O, even to runtime errors if the triangle's maximum angle is greater than 90°. Certainly, there exist other more sophisticated mesh refinement techniques reducing the sharpness of the resolution jumps (e.g. Engwirda 2017). Furthermore, Walko and Avissar (2011) developed a row of algorithms to reduce the sharpness between regions of different edge-bisection refinements. However, the basic idea of this paper is to test the ICON-O performance on computational meshes, which are not geometrically optimal.
The oceanographic background of our experiments is the fact that barotropic tides are the dominant dynamic process along the greatest parts of the world ocean's coasts. Questions related to coastal protection, marine traffic, fishery, coastal ecology and coastal carbon fluxes are strongly connected to coastal tides. Therefore, besides testing the model, we also try to find the direction towards an optimal mesh for the simulation of the global barotropic tide along the coast, whereas this optimisation would be defined with regard to the precision of the coastal tidal wave simulation and the associated computational costs.

The model
ICON-O is a global ocean general circulation model based on finite volume numerics. The surface of the sphere is divided into triangular cells with a C-grid staggered distribution of variables. In the vertical direction, it uses a z-coordinate axis. Besides computing the ocean state, it includes the computation of the dynamics and thermodynamics of sea ice. This is performed by a coupling with the sea ice model FESIM (Danilov et al. 2015), which uses the viscous-plastic rheology after Hibler (1979). For a detailed mathematical and numerical description of ICON-O, the reader is referred to Korn (2017).
In order to parameterise the horizontal turbulent viscosity A H , used for the traditional Laplace formulation of viscosity, the following option was chosen: with Δs being the length of the triangle edge connected to the velocity point and Δx the distance between the two adjacent triangular circumcentres. This linear relation between the grid spacing and the exchange coefficients reminds to Smagorinsky's (1963) approach, though this also considers the velocity shear. The topic of optimizing the velocity closure by biharmonic, or a combination of harmonic and biharmonic operators or flow dependent closures such as the Smagorinsky and the Leith closure is the subject of ongoing research. The coefficient for horizontal tracer diffusion K H is set to zero and only the flux limiter in the flux-corrected transport scheme provides a local diffusive effect. The model uses a constant bottom drag coefficient of 0.003. The vertical turbulent viscosity and diffusivity are parameterised by the modified Richardson number dependent so-called PP scheme (Pacanowski and Philander 1981). Here, the vertical shear is directly computed at the edges (velocity points) whereas the stratification is interpolated from the two ambient cells onto the edge. The water column is resolved by 60 layers, which have a thickness between 15 m within the upper 150 m and 200 m for the depth range below 2000 m. This relatively high vertical resolution near the surface is necessary for a realistic representation of the coastline. The model topography is derived from the 1′ global relief model ETOPO1 (Amante and Eakins 2009). ICON-O uses a static land sea mask; a wetting-drying algorithm is yet not included.

Tidal dynamics
Ocean tides are the results of the gravitational forces of the moon and the sun acting on a unit mass at the earth's surface. Using the geometry given in Fig. 1 and denoting by j either moon M or sun S, we calculate the gravitational force f g at point P with with G being the gravitational constant, m j the moon's or sun's mass, R' the distance between P and the respective celestial body. Furthermore, the geometry gives: with R j being the distance between the earth's centre and the respective celestial body. The angle γ is located between the vector from the earth's centre to the celestial body and the vector from the earth's centre to the reference point P on the earth's surface (Apel 1987).
Hence, it follows This force is mostly balanced by the centrifugal force f c caused by the earth's rotation around the centre of mass of the earth/celestial body system with the angular velocity Ω j . It is With Keppler's third law (complete balance at the earth's centre) if follows Now, we introduce the tidal potential Φ whose negative horizontal gradient is equal to the tide generating horizontal component of the difference between f g and f c .
Based on Eqs. (4) and (7), it follows with Φ 0j = G m j /R j (Apel 1987). Spherical geometry leads to: where φ denotes the geographical latitude of the reference point, α j the declination of the celestial body, and ϕ j its azimuth phase which depends on time and the longitude of the reference point. Finally, α j , ϕ j and R j are determined by the ephemeris of Duffett-Smith and Zwart (2011). The effects of loading and self-attraction (LSA) are neglected in this first version of the module.
We have introduced the horizontal tide-generating acceleration a ! T into the equations of motion and added a new module to the model code which computes the full tidal potential Φ as the sum of the moon's and sun's contribution, Following Love (1909), the horizontal vector a ! T is derived from the tidal potential Φ with: After Thomas (2002), the two Love numbers k 2 and h 2 satisfy: i.e. we assume a 31% reduction of the tide-generating acceleration caused by the tidal deformations of the solid earth.
The reader may wonder why we have implemented the full tidal potential (Eq. 9) and avoided the more traditional and less expensive way of forcing the model by a single or a suite of the most important partial tides (e.g. Stuhne and Peltier 2009). The answer is that the "full approach", apart from being in a sense more straightforward, provides broad frequency tidal dynamics, which include nonlinear interactions between the partial tides. These lead to signals like the M 4 -tide, which can reach substantial amplitudes (Weis et al. 2008). Though these components are not relevant for the present study, we preferred to develop a tidal module, which covers the entire spectrum. The "full approach", the full baroclinic dynamics resolved by the comparatively fine vertical grid resolution with 60 layers and the meteorological forcing (see Section 2.3) are features, which are certainly not necessary for the study of the efficient simulation of barotropic tides. Rather they are development steps of ICON-O preparing future experiments. However, the results of this study should not be influenced by the partly excessive complexity of our computations.
As already stated by Weis et al. (2008) regarding their full tidal potential approach, it should be noted that tidal predictions produced by such a model, even after the inclusion of the LSA terms, will still be far away from the precision of current tidal models using data assimilation techniques like OTIS (Egbert and Erofeeva 2002) or HAMTIDE (Bosch et al. 2009).

Grid generation
The computational mesh's basic geometrical structure is, as the model's name already suggests, the regular icosahedron, which has 12 vertices (one at the north pole, ten evenly distributed at 26.565°N and 26.565°S and one at the south pole), 30 edges and 20 equilateral triangle faces (cells) (Wan 2009). This very coarse global grid is further recursively refined with the edge-bisection method, resulting four triangles for each original one.
If these refinements are performed for all grid cells, the grid conserves the regular icosahedron's uniform, essentially equiangular structure. If only a part of the cells is refined, the grid becomes irregular, i.e. the resolution will vary considerably within the grid. In this case, unrefined cells adjacent to refined cells have to be split-not into four but into two new cells. Thereby, the newly defined vertex of the refined cell, which bisects one of the triangle's edges, a so-called hanging node, is connected to the triangle's opposing vertex. This way, the hanging node is removed and the mesh remains conformal (Walko and Avissar 2011). Figure 2 illustrates these two ways of cell refinement, edge-bisecting (quadtree) and halving.
After such an irregular grid is created one problem remains, which has to be solved before the ICON-O simulation can start: as Fig. 2 shows, the cells in the transition zone between different refinement levels have lost the equiangularity of those cells they were created of. Nevertheless, triangular cells Fig. 1 Geometry for tidal force resolution on a unit mass located at point P on the earth's surface. R j and R' are the vectors between the earth's centre and the celestial body and between P and the celestial body respectively. r denotes the vector from the earth centre to P and γ the angle between r and R j (after Apel 1987) with a maximum angle greater than 90°, i.e. with circumcentres laying outside the triangle, are not suitable for the model numerics and cause immediate runtime errors. These are caused by the distance between the velocity points at the cell boundaries and the tracer points, which are defined at the cell's circumcentre. This distance approaches zero if the maximum angle approaches 90°. Therefore, oriented at the socalled "spring optimisation" by Tomita et al. (2002), an algorithm was developed which contains a stepwise movement of the vertices of those cells with a maximum angle larger than 84°. The accelerating forces of this movement can be thought of as caused by springs, which connect the vertices along the edges. After the maximum angle of all triangles is brought below 84°, a practical value that combines numerical stability with an uncomplicated spring-algorithm, the grid is stored and ready to be loaded by the ocean model.
In the present paper, we discuss the same global tide simulation carried out on different computational meshes. All these meshes are based on a fivefold edge-bisection refinement of the regular icosahedron. Hence, after the fivefold refinement, the number of cells is 20 × 4 5 = 20,480. After interpolating the topography onto the mesh, 4751 cells, either dry cells not neighbouring a wet cell or wet cells not having a connection to the ocean, were removed. Finally, the number of grid cells of our uniform basis mesh, provided with the name uni00 and shown in Fig. 3, is 15,729. The mean horizontal resolution, defined as the spacing of cell centres, of the uni00 grid is 139 km. It varies between 124 and 142 km. The maximum angle of the triangles, a heritage of the original icosahedron, is 72°. On average, the triangles' maximum angle is 64°. Hence, uni00 is an essentially regular grid consisting of essentially equilateral (equiangular) triangular cells. angle of some triangles is greater than 90°. c) Before (red) and after spring-optimisation (black). Now, the maximum angle of all cells is less than 84°. d) Zoom into figure c) The next three uniform grids, uni01, uni02 and uni03, are simply formed by further overall edge-bisection refinements. Hence, the resolution increases to 69 km (uni01), 35 km (uni02) and 17 km (uni03), whereas the grid's triangle shape statistics remain unchanged. The number of cells, i.e. columns consisting of 60 prisms (depth layers), increases to about 60,000 (uni01), 240,000 (uni02) and 950,000 (uni03). Grid uni03 is also shown in Fig. 3.
The remaining 16 grids (experiments) are highly irregular grids. Within these, some areas remain at the 139 km resolution of the uni00 grid, which always is the starting point of the grid generation, whereas other areas experiencing a three-to fourfold edge-bisection refinement, i.e. reaching a resolution of 14 to 7 km. Following Lyard et al. (2006) regarding the computational mesh design for modelling tides, the spatial criteria determining whether a cell is refined are either the coastal, i.e. an oceanic (water depth > 0 m) cell is refined if it contains at least one coastline point defined before in the 1′ topography matrix, or it is the slope criterion, i.e. the oceanic cell is refined if the low-pass filtered absolute value of the horizontal depth gradient exceeds a certain threshold. A third refinement criterion is the neighbourhood of a just refined cell, which leads to a certain broadening of the selected refinement area. In the case of multiple grid refinements, this broadening controls the sharpness of resolution transition zones. Hence, in our grid generation program, a set of three integers: k b , r c and r t defines the structure of the irregular grid, with k b being the number of steps broadening the selected refinement areas, r c the maximal refinement level of coastal cells and r t the maximal refinement level for cell fulfilling the slope criterion. Hence, we decided to name the 16 irregular grid experiments irrk b r c r t .
Although one could guess certain (k b , r c , r t ) settings being favourable for the simulation of coastal tides, e.g. a high value of r c should be crucial, our idea was to test arbitrarily selected settings enabling the experiments to indicate independently the most efficient grid. Hence, some grids focus on coastal refinements (irr241, irr230, irr240) some on topographic slopes (irr313, irr224, irr213), most use sharp resolution gradients but some use rather broad transitions (irr431, irr432). The maximum refinement level is either 3 or 4. The number of cells lays between 205,722 (irr230) and 510,731 (irr242) which are 22% and 54% respectively of the cell number of regular level 3 grid uni03. The reader is referred to Table 1 for more specific information about the different grids. Figure 4 shows histograms of the cell's maximum angle of four different meshes. According to Korn and Linardakis (2018), an equiangular mesh offers the best numerical characteristics. Hence, the uniform grid uni03 shows something of an ideal distribution with the highest frequency in the 63°-64°i nterval and only few cells with a maximum angle above 67°. In contrast, the irregular grid irr241, generated by a fourfold Fig. 3 Selection of computational meshes used in the experiments. uni00: uniform basic low-resolution (124-142 km) mesh. uni03: uniform high-resolution (16-18 km). irr313: irregular with refinements up to 14 km along large topographic gradients. irr241: refinements up to 7 km along the coast. irr431: coastal refinements up to 14 km and smooth transitions. irr230: similar to irr431 but with sharp transitions and complete ignorance of topographic slope refinement coastal and onefold slope refinement, shows a flattened 63-64°maximum and a near-uniform spread of smaller frequencies towards greater angles until it drops abruptly at the 83°-84°interval. The comparison of the histograms of grid irr230 (threefold coastal refinement) and irr240 (fourfold coastal refinement) provides an insight into the mesh geometry's development during the alternating steps of refinement and subsequent spring-optimisation. It shows that the fourth refinement caused a slight frequency reduction in the ideal range between 60°and 67°, a slight increase in the 69°-80°range and another frequency drop between 80°and 83°. Furthermore, the local 80°-81°maximum of the irr230 grid moved to 78°-79°at irr240. Apparently, the additional spring-optimisation operations performed after the fourth refinement step caused a further improvement even of those cells, which had been created during the first three refinement steps. irr240 generally contains more optimisation work.

Experimental setup
The four uniform (uni00-uni03) and the sixteen irregular grids (irrk b r c r t ) are eponyms of the respective ICON-O simulations, which are, apart from the different grids, set up identically. They start at the date 2001-01-01 00:00, use a time step of 60 s and a spin-up period of 30 days for the tidal forcing. The model was operated in the "ocean only" mode, i.e. not coupled with an atmosphere model.
The atmospheric forcing consists of the synthetic cyclic year of Röske (2006), a data set that also provides climatological runoff values. The ocean is assumed initially at rest. The initial temperature and salinity fields were taken from the 0.25°resolution World Ocean Atlas 2013 data set Zweng et al. 2013). Here, winter season fields of period 2005-2012 were selected and interpolated onto the model grid. Thereafter, the temperature data was converted to potential temperature referenced to the sea surface pressure following the algorithm of Fofonoff and Millard Jr (1983), which is based on the work of Bryden (1973). The sea surface salinity (SSS) field derived from the initial data is used by a relaxation scheme preventing the simulated SSS to drift too far away from observations. Soon after the spin-up phase the barotropic tides reach a quasi-stationary state, so that we decided to start a half-hourly output (snapshots) of the sea surface elevation field at 01-05-2001 00:00 and end this at 12-05-2001 00:00. Table 1 Summary of the results of the 20 experiments. k b , r c , and r t are gridding parameters explained in Section 2.2. dx min , dx max denote the minimum and maximum grid spacing (distance between neighbour cell centres). α mean denotes the mean value of the maximum angle of each triangle cell, P 70 , P 75 , P 80 the percentage of cells with a maximum angle above 70°, 75°and 80°respectively. ē denotes the mean difference between the simulated coastal M2-amplitude and those from the TPXO-8 data set. sx denotes the standard deviation of ē. c is the Pearson correlation coefficient between the modelled and the TPXO-8 data. scRMS denotes the "scaled centred root mean square difference" defined in Section 3 Experim. k b r c r t Number of cells dx min [km]

Results
The purpose of our experiments is to test the model's ability to simulate global barotropic tides on different computational meshes. For this, we concentrate on the most important tidal constituentthe M2 tide and compare the respective coastal amplitudes with the corresponding values from the TPXO-8 data set (Egbert and Erofeeva 2002). The TPXO-8 has a spatial resolution of 2′ and is created by a generalised inverse tidal model, which assimilates the TOPEX/Poseidon satellite altimetry data. The precision of the TPXO-8 amplitudes is less than 1 cm for the open ocean and less than 4 cm for high coastal amplitudes. Therefore, it is suitable to serve as a near observational reference when calculating the coastal M2 amplitude error e of our simulations. This is done in the following way: first, the simulated M2 amplitude for each grid cell is determined by a discrete Fourier transformation of the output sea level time series (Fig. 5). In a second step, for each coastal point i of the 2' TPXO-8 matrix (213,575 points), the nearest wet ICON cell k is determined. The two corresponding M2 amplitudes, TPXO i and ICON k , define the error e i : Consequently, the mean error ē, the standard deviation sx, and the Pearson correlation coefficient c are computed with c ¼ Fig. 4 Probability histograms of the triangular cell's maximum angle α max for four different meshes. uni03: uniform mesh after a threefold refinement with a maximal α max of 72°and a mean α max of 64.1°. irr241: irregular mesh after fourfold coastal refinement. irr230 irregular mesh after a threefold coastal refinement. irr240: like irr230 but with one further coastal refinement. The bars' overlap is coloured dark grey whereas N = 213,575 denotes the number of analysed coastal points, and TPXO , ICON the mean values of TPXO i and ICON k . Finally, we define the scaled centred root mean square difference scRMS with Hence, scRMS can be interpreted as the reduced standard deviation of e in the case of adjusted values of ICON k so that ICON = TPXO. The oceanographic meaning of scRMS is that this error metric rather describes the model's ability to capture the structure of the M2 amplitude field neglecting a general over-or underestimation of amplitudes that could be handled quite easily within a model fitting procedure, e.g. through a variation of friction parameters. Fig. 5 Amplitude of the M2 tidal constituent given by the TPXO-8 dataset (Egbert and Erofeeva 2002) and simulated in experiment uni00 (low-resolution uniform mesh), uni03 (highresolution uniform mesh) and irr241 (mesh with coastal refinement). The bottom row shows the difference between the ICON simulations uni03 and irr241 and the TPXO-8 data

Uniform grids
The first simulation was performed on the uni00 grid, i.e. the uniform triangular mesh with a mean resolution of 139 km, which forms the basis of all other grids used within this study. The comparison of the resulting M2 amplitude field with the TPXO-8 data shows that the key structures were correctly reproduced (Fig. 5). For example, even the two amphidromic points north-west and east of Iceland are reproduced. Generally, the M2 amplitude is over-estimated by the model. Along the coasts we find e = 28 cm and scRMS = 61 cm. This over-estimation even increases in the uni01 experiment after the first edge-bisection refinement to a mean resolution of 69 km. Here, e rises to 32 cm. However, the increase of the correlation c from 0.42 to 0.50 and the decrease of scRMS to 58 cm indicate an improved reproduction of the M2 amplitude structure caused by the higher spatial resolution. This trend continues within experiments uni02 (35 km resolution) and uni03 (17 km): scRMS is further reduced to 51 cm and 48 cm respectively, c increases to 0.61 and 0.65 and also the initial over-estimation now develops in the right direction with e = 23 cm in uni02 and 11 cm in uni03 ( Fig. 5 and Table 1).
In order to illustrate the functionality of the entire tide model also comprising all other tidal constituents, we have plotted the uni03 sea surface elevation at the position of Reykjavik, Iceland, and compared it with the astronomical tide prediction of the NOAA (2000). Figure 6 shows a temporally variable amplitude difference between 5 and 50 cm and a rather constant phase error of around − 1.7 h. Figure 7 shows the scRMS of all 20 experiments as a function of the number of grid cells. The uniform grid experiments, uni00 -uni03, are symbolised by red dots connected by a red line. Regarding the irregular grid experiments, irrk b r c r t , symbolised by green and blue dots, the key criterion is whether their points lay above or below the red line. A position above the red line, obtained for irr213 and irr313, means that the applied grid refinement is ineffective because a uniform grid would deliver a smaller scRMS with less computational effort, i.e. with a smaller number of grid cells. Table 1 shows that both grids, irr213 and irr313, were made with the setting r c = 1, i.e. with only one coastal refinement step. Remarkably, the scRMS of both experiments (0.61 m) is above that of uni01 (0.58 m) which has the same coastal resolution but only one fifth of the number grid cells.

Irregular grids
Obviously, and not surprisingly, a reduction of scRMS, i.e. of the coastal M2 amplitude error cannot be achieved without sufficient coastal mesh refinement. However, apart from the two outliers, irr213 and irr313, who have disregarded the coastline with r c = 1, all irregular grid experiments are below the red line, i.e. they show at least some efficiency. Even the grids, which still have more focus on topographic slopes but with a twofold coastal refinement, irr223, irr224, irr323 and irr423, show a scRMS of 0.50 m which is slightly below that of uni02 with 0.51 m. The grids with a threefold coastal refinement, irr230, irr231, irr232, irr233, irr331, irr332 and irr431, lead to a scRMS of 0.48 m, i.e. they reach the precision of the uni03 grid (scRMS = 0.48 m) but using only one third of the number of grid cells. If we define the efficiency as the slope of the line connecting the root point uni00 with the experiment's point then irr230 with d(scRMS)/dn = −6.63 × 10 −7 m is the most efficient grid. Its setting is (k b , r c , r t ) = (2,3,0), i.e. its refinements are performed exclusively along the coast, refinement along topographic slopes is ignored and sharp resolution gradients are accepted. One could describe irr230 as the most spartan of our settings. Finally, the grids with a fourfold coastal refinement, irr240, irr241 and irr242, lead to scRMS = 0.45 m, 0.46 m and 0.45 m respectively. These values are clearly below the uni03 value of 0.48 m. Hence, the higher coastal resolution yielded a more precise simulation while using only around 50% of the grid cells compared to uni03. Within this group, the most efficient setting is that of irr241 with (k b , r c , r t ) = (2,4,1) and d(scRMS)/dn = − 3.45 × 10−7 m closely followed by irr240 with (k b , r c , r t ) = (2,4,0) and d(scRMS)/dn = −3.44 × 10 −7 m. Here, the "luxury" of around 12.000 additional cells at sharp topographic slopes in irr241 and a slightly better mesh geometry (Table 1) is paying off.
In order to examine the influence of grid irregularities on the model solution, we have plotted the sea surface elevation southwest of Iceland at various times simulated on the uniform uni03 and on the irregular irr230 grid (Fig. 8). The comparison shows that the tidal wave arriving from the southeast was reproduced on both grids in a similar way. The deviation between both solutions are below 20 cm. Apart from a weak stationary wavelike structure in the difference fields along the grid resolution transition zone at 63°N between 21°W and 24°W (right column of Fig. 8), disturbances like reflections or sharp stationary anomalies directly linked to the resolution transition zone of irr230 are not visible. Figure 9 compares the according sea surface elevation time series at a location near the transition zone. Again, a clear numerical disturbance is not visible. The temporally variable difference remains below 10 cm and the phases differ by around 30 min.

Summary
Our experiments show that the coastal amplitude error decreases with increasing coastal resolution. In fact, this dependency is nearly independent of the resolution chosen for the inner ocean. The irregular grid irr230 and the uniform grid uni03, both resolving the coast with around 17 km, yield nearly the same coastal amplitude error (scRMS = 0.486 m and 0.482 m respectively) though their resolution of the inner ocean differs sharply (148 km and 17 km respectively). The number of grid cells, i.e. columns consisting of 60 layers, is 205,722 for irr230 and 945,995 for uni03, i.e. the irregular grid irr230 reduced the computational costs of a simulation with a scRMS below 0.49 m by 78%. The grid irr241 contains one further coastal grid refinement and consequently yields a scRMS = 0.454 m which is clearly below the uni03 value but . Green and blue dots denote simulations on irregular grids (irrk b r c r t ) with a maximal resolution of 14 km and 7 km respectively still uses a number of grid cells (472,235) which is 50% lower than that of uni03.
When defining the efficiency as the scRMS reduction divided by the number of additional grid cells related to the uni00 experiment, the irregular grids irr241 from the 148-7 km class, and irr230 from the 148-14 km class have proven to be the two most efficient grids. Both grids are based on economic settings, i.e. high-resolution gradients were allowed and the refinements were performed almost exclusively along the coastline. mesh. The middle column shows the results of experiment irr230 using a 14-km resolution only near-shore and up to 148 km in the basin's interior. The right column shows the difference between both experiments. The first row shows the different grid structures. The two cells marked in red provided the time series shown in Fig. 9 We have developed a module for the newly developed global ocean general circulation model ICON-O (Korn 2017), which computes the full tidal potential and the according horizontal body forces based on the ephemeris of Duffett-Smith and Zwart (2011). The first tests presented here show an essential realistic simulation of global tides. Without any tuning, the coastal M2 amplitudes are reproduced with a scaled rms error between 62 and 45 cm depending on the model's coastal resolution between 140 and 7 km. A further increase of the horizontal resolution will certainly lead to an additional decrease of the amplitude error. However, a better representation of topography by a higher vertical resolution and the tuning of friction parameters (bottom friction, horizontal momentum exchange) should also be effective and should substantially reduce the tidal phase errors (e.g. Davies and Aldridge 1993). These show an error of around 2 h (Figs. 6 and 9), with the simulated wave being ahead of the observed. Future versions of the module will include the effects of loading and self-attraction, which is expected to reduce the phase error (M. Thomas, personal communication, June 11, 2020).
Further tests relate to the performance of ICON-O simulating the global barotropic tide on 20 different triangular grids. To summarize the results of these experiments into one sentence: the application of highly irregular grids within ICON-O is feasible and can be very efficient. The results provide no evidence for severe grid-related numerical errors like reflections at resolution transition zones or sharp stationary, grid-related pressure anomalies. Considering the applied rather straightforward mesh refinement technique with its partly substantial deformation of cells, we conclude that a wide range of other methods, e.g. JIGSAW-GEO (Engwirda 2017) or the multiple transition rows approach of Walko and Avissar (2011) can be used as well. In addition, the span of resolutions within one grid could probably be driven further, i.e. beyond the used 148-7-km range. A local horizontal resolution in the range of 1 km would enable the simulation of internal waves. The effort of simulating baroclinic waves generated by barotropic tides, already proposed by Stuhne and Peltier (2009), could be taken up again by an application of ICON-O with its improved features regarding smallscale noise of the horizontal divergence field (Korn 2017).
Our results show that the simulation of the coastal M2 amplitude, i.e. the error metric we have used in this paper, predominantly depends on the coastal grid size. This becomes clear if one considers a bay with a tide amplifying shape and the expected error of a model, which does not resolve this bay. Of course, a lower resolution of the ocean basin's interior is linked to losses of local precision. This can be seen at the offshore M2 amplitude error of experiment irr241 compared to that of uni03 (Fig. 5, bottom row). Hence, other error metrics of tidal simulations will certainly point to a more important role of the resolution of the inner ocean and of a refinement along topographic slopes like Lyard et al. (2006) propose.
Both aspects, the numerical feasibility and the successful efficiency enhancement should encourage future users of ICON-O to apply mesh refinement techniques even if their  Fig. 8. The red curve shows results from experiment uni03 (uniform high-resolution mesh), the blue curve results from experiment irr230 (irregular mesh with coastal refinement) and the green curves the TPXO-8 values referring to the above-mentioned position research question is not as spatially selective as our coastal M2 amplitude. Lyard et al. (2006) have shown that the global simulation of tides generally benefits from topography-related mesh refinements. Other obvious cases of usefulness are transports through straits, processes along shelf breaks, ridges or continental slopes and many questions of regional oceanography.