THE DISCRETIZATION GRID AND THE DISCRETE OPERATORS

The discretization method is defined for concretenes on a special case of Delaunay triangulation on the sphere, i.e. the icosahedral geodesic grid described e.g. in [4]. The main reasons for the choice of this type of grid is its quasi-uniform coverage of the sphere, which solves automatically the pole problem and avoids very high Courant numbers at the poles. Furthermore its hierarchical structure provides a very natural setting for local grid refinement on nested grid hierarchies. Finite element approaches based on such geodesic grids have been introduced in [10], [11],[15]. Finite volume approaches were instead presented in [13], [25], [26]. The icosahedral construction yields a Delaunay triangulation of the sphere to which a Voronoi tessellation is naturally associated (see e.g. [24]), which consists of convex spherical polygons (either pentagons or hexagons, see figure 1). The Delaunay cells of the icosahedral grid are all triangles. For each side of a Voronoi cell there is a unique orthogonal side of the Delaunay cell associated to it.

The mass and vorticity preservation properties in ICON are achieved by use of triangular Delaunay cells on the sphere as control volumes and of the dual Voronoi cells (pentagons or hexagons) as control volumes for vorticity. The orthogonality of the primal and dual grid edges allows to use simple approximations of the gradient and rotation operators, in the framework of a C-type staggering of the discrete variables. This represents a major changes with respect to the discretization employed e.g. in GME (global model of Deutscher Wetterdienst, [20], GME shallow water version), where an A approach was used and discrete variables were defined at the vertices of the Delaunay grid and the orthogonality of primal and dual grids was not exploited.

Figure 1: Delaunay (triangles) and Voronoi (hexagons and pentagons) cells on the sphere obtained by diadic refinement of the regular icosahedron.

Some notation to describe the grid topology and geometry will now be introduced (see figure 2). Let then $i$ denote the generic cell of the Delaunay grid. Let ${\cal E }(i)$ denote the set of all edges of cell $i$ and ${\cal C }(i)$ the set of all cells which have edges in common with cell $i.$ The gridpoint associated to cell $i$ will also be referred to as the cell center. The generic vertex of a cell, which is also the center of a cell in the dual grid, is denoted by $v.$ ${\cal C }(v)$ denotes the set of all cells of which $v$ is a vertex and ${\cal E }(v)$ denotes the set of all edges of the dual cell whose center is vertex $v.$ The area of cell $i$ is denoted by $A_i,$ while the area of the dual cell is denoted by $A_v.$ Let then $l$ denote the generic edge of a cell. It is to be remarked that this index can be assigned at the same time to an edge of the primal grid and the edge of the dual grid, which by construction intersects the primal grid edge at its midpoint. The number of edges is actually equal for both grids.

Figure 2: Dual grid

The length of the edge $l$ of a cell is denoted by $\lambda_l$ and the distance between the centers of the cells adjacent to edge $l$ (i.e., the length of a edge of the dual cell) is denoted by $\delta _l.$ At each edge, a unit vector ${\bf N}_{l}$ normal to the edge $l$ is assigned. ${\bf T}_{l}$ denotes the unit vector tangential to the edge $l,$ chosen in such a way that ${\bf N}_{l} \times {\bf k}_l = {\bf T}_{l}$ holds, where ${\bf k}_l$ denotes the radial outgoing unit vector perpendicular to the tangent plane at the intersection of primal and dual edge $l.$ Furthermore, for each cell edge, the unit vector pointing in the outer normal direction with respect to cell $i$ is denoted by ${\bf n} _{i,l}.$ Unit vectors ${\bf n} _{v,l}$ are also introduced, as pointing in the outer normal direction with respect to the dual cell $v.$ The corresponding tangential vectors ${\bf t} _{v,l}$ are defined so that ${\bf n}_{v,l} \times {\bf t}_{v,l} = {\bf k}_l.$ It can be seen that, by simple geometric arguments, one has ${\bf N}_{l}\cdot {\bf t}_{v,l} = {\bf T}_{l}\cdot {\bf n}_{v,l} .$

In order to develop an analog of the rectangular C-type staggering (see e.g. [2]) on the Delaunay grids, the mass points are defined as the centers of the grid cells, while the velocity points are defined for each cell edge as the intersection between the edges of the Voronoi and Delaunay cells (see figure 2). By construction, each of these points is equidistant from the centers of the Voronoi cells at the ends of that edge. A velocity point is also the intersection of the edge of the cell with the arc connecting the centers of the cells adjacent to that edge. These points are the locations of the discrete normal velocity components with respect to the cell edge. Given the edge $l$ of a cell, the adjacent cells are denoted by the indexes $i(l,1)$ and $i(l,2),$ respectively. The indexes are chosen so that the direction from $i(l,1)$ to $i(l,2)$ is the positive direction of the normal vector ${\bf N}_{l}.$ Vertex indexes $v(l,1)$ and $v(l,2)$ can also be defined analogously, so that the direction from $v(l,1)$ to $v(l,2)$ is the positive direction of the vector ${\bf T}_{l}.$ Given a generic discrete vector field ${\bf g}$ on the sphere, its value at a velocity point can be represented as ${\bf g}_l= g_l{\bf N}_l +\hat g_l {\bf T}_l ,$ where $g_l, \hat g_l$ denote the normal and the tangential components, respectively. In a C grid discretization approach, the discrete prognostic variables considered are the value of the height field $h_i$ at the mass points (interpreted as a cell averaged value) and the normal velocity components $u_l.$ The tangential velocity components, which are needed e.g. for the computation of the Coriolis force term, must be reconstructed. It is to be remarked that this type of grid arrangement yields a spatial discretization which is very similar to the Raviart-Thomas finite element of order 0 (see e.g. [23]).

Given these definitions, discrete operators can be introduced, which will be employed then to define the proposed numerical algorithm. Thanks to the orthogonality between the cell edges and the arcs connecting the cell centers, the directional derivatives in the normal and tangential directions are easily approximated as

$$( \delta_{\nu}h)_l=( h_{i(l,2)} - h_{i(l,1)})/ \delta_{l} , ( \delta_{\tau}h)_l=( h_{v(l,2)} - h_{v(l,1)})/ \lambda_{l}.$$ (5)

Again because of the same orthogonality, discrete divergence and curl operators are now introduced in the context of the C grid staggering outlined above. These operators are defined as acting on a set of values $g_l$ assigned at the edges of the Voronoi-Delaunay grid. These are to be interpreted as the components of a vector field ${\bf g}_l$ normal to the cell edges, i.e. $g_l={\bf g}_l\cdot {\bf N}_l.$ The discrete divergence and curl operator can be naturally defined as follows:

$${\mathrm div}(g)_i = \frac {1}{A_i}\sum ^{} _{l\in {\cal E}(i) } ... ...{A_v}\sum ^{} _{l\in {\cal E}(v) }{g_l} {\bf N}_l \cdot {\bf t}_{v,l}\delta_l .$$ (6)

The discrete analogous of the Helmholtz decomposition theorem was proven in [22] for this type of grid arrangement and discrete operators. It should be observed that the velocity points are not exactly equidistant from the adjacent Delaunay grid cell centers. As a result, the difference operators described above are only first order accurate. However, grid optimization procedures such as those introduced in [13,14] can partly cure this problem, by reducing the off-centering to rather small values.

In table 1 some characteristic quantities of the triangular icosahedral grid at various resolutions after Heikes-Randall optimization [13,14] can be found.

Table 1: The triangular icosahedral grid at various resolutions after Heikes-Randall optimization

grid level	N. of triangular cells	N. of edges	Av. edge(deg)	Av $\lambda_{l}$(km)	Av $\delta_{l}$(km)
1	20	30	63.435	7053.888	4649.255
2	80	140	33.859	3765.050	2251.109
3	320	480	17.232	1916.228	1116.204
4	1280	1920	8.654	962.284	556.848
5	5120	7680	4.331	481.643	278.255
6	20480	30720	2.166	240.882	139.105
7	81920	122880	1.083	120.449	69.550
8	327680	491520	0.542	60.225	34.755