# Compact symmetric Poisson equation discretization for non-hydrostatic sigma coordinates ocean model

Roullet Guillaume, Molemaker M. Jeroen, Ducousso Nicolas, Dubos Thomas

► To cite this version:

Roullet Guillaume, Molemaker M. Jeroen, Ducousso Nicolas, Dubos Thomas. Compact symmetric Poisson equation discretization for non-hydrostatic sigma coordinates ocean model. Ocean Modelling, 2017, 118, pp.107-117. 10.1016/j.ocemod.2017.09.001. insu-03682732

HAL Id: insu-03682732

https://insu.hal.science/insu-03682732v1

Submitted on 4 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Compact symmetric Poisson equation discretization for non-hydrostatic sigma coordinates ocean model

Roullet Guillaume <fnref n="1" />,*, Molemaker Maarten Jeroen <fnref n="1" />,<fnref n="2" />, Ducousso Nicolas <fnref n="1" />, Dubos Thomas <fnref n="3" />

## Abstract :

In anticipation of relaxing the hydrostatic assumption in a sigma coordinates primitive equations ocean model, we show how a projection method can be designed with the use of a compact symmetric 15-point stencil for the Poisson equation. This is achieved by recognizing that, owing the non-orthogonality of the grid, the velocity has a contravariant and covariant set of components. The two sets play a different role in the primitive equations: the contravariant components enter the definition of the model fluxes, whereas the covariant components experience the forces and in particular the pressure gradient. By treating these two sets separately, the discretized gradient and divergence operators are simple finite differences. The two sets of components are related via a linear transformation, the metric tensor, which is entirely determined by the kinetic energy. We show how the spatial discretization of the kinetic energy fully controls the Poisson equation discretization, including its boundary conditions. The discretization of the Poisson equation is shown to converge at second order and to behave as well or better than alternative methods. This approach is a prerequisite to implement later an efficient Poisson solver, such as a multigrid algorithm.

## Highlights

* A new formulation of the non-hydrostatic primitive equations in sigma coordinates is proposed. * Momentum and flux are two different objects related via a linear transformation. * The Poisson equation for pressure is discretized with a compact symmetric stencil. * The approach is a prerequisite to implement a fast Poisson solver for a forthcoming non-hydrostatic version of ROMS.

**Keywords** : Non-hydrostatic, Sigma coordinates, Poisson equation

<footnote n="1">¹ Univ. Brest, CNRS, IRD, Ifremer, Laboratoire d'Océanographie Physique et Spatiale, IUEM, Brest, France.</footnote>

<footnote n="2">² UCLA, Los Angeles, USA</footnote>

<footnote n="3">³ IPSL/Lab. de Météorologie Dynamique, École Polytechnique, Palaiseau, France</footnote>

<footnote>* Corresponding author : Guillaume Roullet, email addresses : roullet@univ-brest.fr ; guillaume.roulet@univ-brest.fr</footnote>
1. Introduction

With the growing interest on processes at submesoscale (Nikurashin et al.,

2013; Molemaker et al., 2015; Callies et al., 2015; Gula et al., 2015) the ocean
3
4 modelling community needs non-hydrostatic (NH) ocean circulation models.

This will allow to investigate in realistic contexts the full range of internal

waves, symmetric instability, convection processes, etc. and to bridge the gap

with LES (Sullivan et al., 1994). A few circulation models do already handle

the non-hydrostatic physics properly: MITgcm (Marshall et al., 1997), POM

(Kanarska and Maderich, 2003), ROMS (Kanarska et al., 2007), Symphonie

(Auclair et al., 2011), GETM (Klingbeil and Burchard, 2013). However, NH
simulations on large grids, in realistic regional configurations are still awaited.

It is timely to resolve this. The reason for the NH circulation models rarity

is that relaxing the hydrostatic assumption, while keeping the Boussinesq

assumption, is not straightforward. At least two levels of difficulties must

be faced: i) the requirement to solve a 3D Poisson equation for the pressure

at each time step, made more complex in the case of sigma coordinates by

ii) the non-orthogonality of the grid. This second difficulty is absent in z-

coordinates models. Incidentally, for models with a time-splitting on the

free surface, a third difficulty is to properly handle the coupling between the
non-hydrostatic pressure and surface waves (Auclair et al., 2011).

The first problem amounts to solving a linear system of equations with
as many equations as number of grid points. For the intended purpose
of high-resolution turbulent simulations we have in mind grids as large as
2000 × 2000 × 1000 grid points, corresponding to $N \sim 10^9$ coupled equations.
It is a classical well known problem of High Performance Computing (HPC).
A possible way to circumvent the problem consists in relaxing the incom-
pressibility constraint, which forces to cope with sound waves. The pressure
can then be computed by integrating in time the compressible physics using a
time-splitting technique. In the atmospheric community, this is the approach
implemented in WRF (Skamarock and Klemp, 2008). In the ocean the ratio
between advective speed and the phase speed of sound waves is much more
unfavorable, leading to a required time split that is two orders of magnitude
larger. An approach that uses artificial impressibility to slow down the phase
speed of sound waves and hence reduce the required ratio in split time steps
is currently under development as part of the French regional ocean modeling
project CROCO, but not available in a published manuscript (Francis Au-
clair, private communication). Klingbeil and Burchard (2013) follow a similar
approach where a time-splitting on the baroclinic time step is introduced to
achieve the incompressibility constraint. These methods trade the price of
solving the Poisson problem with accepting a non-strict divergence-freeness
of the flow.

Otherwise, the numerical implementation necessarily requires the use of
an iterative solver to solve the Poisson equation. Many different numerical
methods exist but only a few of them succeed at maintaining a uniform
convergence rate for large grids. This is the case of the multigrid technique whose computational cost scales as $N \log N$, which is considered to be the optimal scaling. It has been shown recently to perform very well with large grids on massively parallel cluster, up to $10^{10}$ degrees of freedom and 65536 cores, in the context of atmospheric modelling (Müller and Scheichl, 2014).

To get the best performances, the stencil for the Poisson equation needs to be symmetric and compact. These two properties reduce the number of floating point operations, the amount of data storage and the data movement between the cpu and the memory. In modern computer architecture, data movement is the bottleneck for high performances (Williams et al., 2009).

Any numerical method that reduce data transfer is worth it. In the present study, the compactness means a 15-point stencil, compared to the usual 25-point one (Auclair et al., 2011). The symmetry allows to almost halve the number of coefficients for the matrix since only 8 coefficients are needed for a 15-point stencil (one main and 7 lower diagonals). Overall the present matrix offers a reduction by a factor 3 on the number of matrix coefficients compared to existing methods. In anticipation of implementing such solution into a sigma-coordinates model we present here a way to achieve compactness and symmetry for the Poisson equation.

The second problem is due to the non-orthogonality of the grid. One consequence is that the horizontal pressure gradient can not be computed as a simple finite difference between two adjacent horizontal grid cells. Chain rule and vertical interpolation should be used (Shchepetkin and McWilliams, 2003) resulting in wide stencils for the gradient operator. The divergence of vertical flux behaves similarly because the vertical flux involves horizontal
velocities. Since the Poisson equation arises as the successive action of the gradient and the divergence, the resulting stencil for the Poisson equation is tall in the vertical, involving overall 25 points (Auclair et al., 2011). The Poisson equation couples points over 5 levels in the vertical and 3 in both horizontal directions. The stencil can be made symmetric by defining the gradient as minus the adjoint of the divergence. This is the so-called compatible discretization (Taylor and Fournier, 2010). It is also possible to discretize the Poisson matrix directly and independently of how the gradient and the divergence are discretized (Kanarska et al., 2007). In that case, the stencil is compact with 15 points. The major drawback of such a method is that it fails at maintaining the compatibility between operators which impacts the energy conservation.

In this paper we show how the two problems can be solved jointly. The first step is to recognize that because of the non-orthogonality of the coordinates system, the velocity has two sets of components: a contravariant one and a covariant one. The two sets are related via a linear transformation, the metric tensor, which is invertible. The metric tensor is completely determined by the kinetic energy. The two sets of components play a different role and must be treated distinguishably. The contravariant components appear in the definition of fluxes (volume, tracer, momentum), whereas the covariant components appear in the force budget. Using the framework of the discrete differential geometry (Desbrun et al., 2008) we show that the contravariant components must be discretized at the cell faces of the primal grid whereas the covariant components should be discretized at the edges of the dual grid (the lines joining the cell centers). Such discretization is in line with the
C-grid staggering, and the finite volume discretization for tracers. It also completely adheres to Thuburn and Cotter (2012) way to discretize shallow water equations on non-orthogonal grids. The divergence and the gradient operators have then trivial discretization: they are simple finite differences.

In the spirit of Molemaker et al. (2005) and Dubos et al. (2015) we then deduce the discrete Poisson equation directly from the discretized kinetic energy and the discretized divergence operator. The price of this approach is to prognose either the covariant components or the contravariant ones. For sake of completeness we present the momentum equations written for the covariant components in both vector-invariant form and in flux-form. The vector-invariant form is similar to a standard Cartesian formulation whereas the flux-form involves additional pseudo-force terms. Up to now, sigma-coordinates models prognose Cartesian components of the velocity, namely horizontal and vertical, which are a blend of contravariant and covariant components. This causes unnecessary wide stencil for the divergence and the gradient and finally a too wide stencil for the Poisson equation. By explicitly and clarifying the subtleties induced by non-orthogonal coordinates, this paper offers a roadmap for an efficient implementation of a projection method.

The paper is organized as follows. We define in section 2 the contra- and covariant sets of components for the velocity and present the model continuous equations. In section 3, we present the spatial discretization of the variables, the divergence and the discretized kinetic energy and derive all the other discretizations from them. In section 4 we discuss the properties of the discretization on a test-case. A conclusion is given in section 5.
## 2. Formulation

As stated in the introduction, the objective of this study is to provide a discrete Poisson operator that can be solved efficiently by iterative methods on very large grids. In this section, we will introduce a dual representation of the velocities on the grid, which will prove to be advantageous in arriving at a compact, symmetric, discrete operator. Subsequently, we will present the governing equations of motion for these prognostic variables.

### 2.1. Flux and momentum

In orthogonal curvilinear coordinates $(\xi, \eta, z)$, the velocity of a fluid par- cell is naturally expressed as $\mathbf{u} = h_{\xi}\dot{\xi}\mathbf{i} + h_{\eta}\dot{\eta}\mathbf{j} + \dot{z}\mathbf{k}$, where the dot denotes the Lagrangian time derivative, $\mathbf{i}, \mathbf{j}, \mathbf{k}$ are the unit vectors in horizontal and vertical directions, and $(h_{\xi}, h_{\eta})$ are the Lamé coefficients in the horizontal di- rections. We define $(u_c, v_c, w_c) = (h_{\xi}\dot{\xi}, h_{\eta}\dot{\eta}, \dot{z})$ and denote them, with a slight abuse of usage, the Cartesian velocity components. The terrain following sigma coordinates are $(\xi, \eta, \sigma)$, for which the vertical position is $z(\xi, \eta, \sigma, t)$. The time dependency in $z$ is driven by the sea surface height variations and yields a breathing grid. Using the chain rule, this causes

$$ \mathbf{u} = u_c \mathbf{i} + v_c \mathbf{j} + (u_c s_{\xi} + v_c s_{\eta} + h_{\sigma} \dot{\sigma} + \partial_t z) \mathbf{k}, \quad (1) $$

where

$$ s_{\xi} = \frac{1}{h_{\xi}} \frac{\partial z}{\partial \xi}, \qquad s_{\eta} = \frac{1}{h_{\eta}} \frac{\partial z}{\partial \eta} \quad (2) $$

are the slopes of $\sigma$-surfaces, $h_{\sigma}$ is the Lamé coefficient in the vertical direction and $\partial_t z = \partial z / \partial t$ accounts for the breathing of grid. We denote $(U, V, W) =$
140 (h_ξ ˙ξ, h_η ˙η, h_σ ˙σ) the sigma coordinate velocity components. They are related to the Cartesian velocity components by

141

$$
\begin{align}
u_c &= U \nonumber \\
v_c &= V \tag{3}
\end{align}
$$

$$
w_c = U s_{\xi} + V s_{\eta} + W + \partial_t z .
$$

142 Hereafter we distinguish

$$
W^a = W + \partial_t z = w_c - U s_{\xi} - V s_{\eta}, \quad (4)
$$

the vertical velocity relative to a fixed grid, from W, the vertical velocity
through the breathing grid. In the sequel, we refer to (U, V, W<sup>a</sup>) as the con-
travariant components of velocity. Notice that it is a slight abuse of language,
since strictly speaking the contravariant components of velocity are just the
material derivatives of the coordinates, i.e. (ξ, ˙η, ˙σ), while (U, V, W<sup>a</sup>) incor-
porate the Lamé coefficients (h<sub>ξ</sub>, h<sub>η</sub>, h<sub>σ</sub>). With this definition, (U, V, W<sup>a</sup>) are
close to the Cartesian components typically manipulated by existing hydro-
static models, which simplifies the implementation of our method in such
models. A similar caveat applies to the covariant velocity components de-
fined below.

152

143

144

145

146

147

148

149

150

151

152

The velocity can be written in vectorial form as

153

$$
\mathbf{u} = U (\mathbf{i} + s_{\xi} \mathbf{k}) + V (\mathbf{j} + s_{\eta} \mathbf{k}) + W^{a} \mathbf{k}, \quad (5)
$$

for which it is clear that (i + s_ξ k, j + s_η k, k) are the basis vectors aligned with
the coordinates (ξ, η, σ). Interestingly there are two possible interpretations
for the contravariant components: i) according to (3) and (4), U and V are
154
155
156

154

155

156
horizontal components and $W^a$ is the vertical component normal to sigma iso-surfaces; ii) according to (5), $U$ and $V$ are along coordinates lines and $W^a$ is the vertical component. This dual interpretation is a source of confusion. The first interpretation is motivated by discrete differential geometry (Desbrun et al., 2008), the second interpretation is more in line with the physicist way of seeing the contravariant components. In this paper we will adopt the first interpretation and motivate it more thoroughly in section 3.

Now that we have the contravariant components of the velocity we can define the covariant components<fnref n="1" />. For that we use the kinetic energy $K$ as a natural inner product, allowing to define the metric tensor. In terms of $(U, V, W^a)$ $K$ reads

$$K = \frac{1}{2} [U^2 + V^2 + (Us_{\xi} + Vs_{\eta} + W^a)^2], \qquad (6)$$

where we have used $W^a$, not $W$, because the grid motion should not be accounted for in the kinetic energy. The covariant set of components $(u, v, w)$ is defined such that

$$K = \frac{1}{2}(uU + vV + wW^a). \qquad (7)$$

<footnote n="1">¹As for the contravariant velocity components, our definition of covariant components is slightly abusive. The standard covariant components are obtained by multiplying ours by the corresponding Lamé coefficients. Again, our definition results in $(u,v,w)$ being closer to the Cartesian components, simplifying the implementation in existing hydrostatic models while ensuring key relationships (7,8).</footnote>
171 We can find these components by differentiating K with respect to U, V, W<sup>a</sup>

$$
\begin{align}
u &= \partial K / \partial U = U(1 + s_{\xi}^{2}) + V s_{\xi} s_{\eta} + W^{a} s_{\xi} \nonumber \\
v &= \partial K / \partial V = U s_{\xi} s_{\eta} + V(1 + s_{\eta}^{2}) + W^{a} s_{\eta} \tag{8} \\
w &= \partial K / \partial W^{a} = U s_{\xi} + V s_{\eta} + W^{a}. \nonumber
\end{align}
$$

172 The covariant components can be related to the Cartesian ones by combining
173 equations (3) and (8)

$$
\begin{align}
u &= u_c + w_c s_\xi \nonumber \\
v &= v_c + w_c s_\eta \tag{9} \\
w &= w_c. \nonumber
\end{align}
$$

174 We may use these components to rewrite the velocity in vectorial form

$$
\mathbf{u} = u\mathbf{i} + v\mathbf{j} + w(\mathbf{k} - s_{\xi}\mathbf{i} - s_{\eta}\mathbf{j}). \quad (10)
$$

175 Again we have a double interpretation for the covariant components: i) ac-
176 cording to (9) (u, v, w) are aligned with the coordinates lines and in particular
177 w is the vertical component ; ii) according to (10), (u, v, w) are orthogonal
178 to coordinates lines. To be consistent with our previous choice, we need to
179 chose the first interpretation.

180

180 The system (8) defines the metric tensor relating contravariant to co-
181 variant components. It is a 3 × 3 full matrix. Interestingly the inverse
182 transformation has less entries

$$
\begin{align}
U &= u - w s_{\xi} \notag \\
V &= v - w s_{\eta} \tag{11} \\
W^{a} &= -u s_{\xi} - v s_{\eta} + w (1 + s_{\xi}^{2} + s_{\eta}^{2}) . \notag
\end{align}
$$
We define

$$ \mathbf{T} = \begin{pmatrix} 1 & 0 & -s_{\xi} \\ 0 & 1 & -s_{\eta} \\ -s_{\xi} & -s_{\eta} & 1 + s_{\xi}^{2} + s_{\eta}^{2} \end{pmatrix} \qquad (12) $$

the matrix transforming $(u, v, w)$ into $(U, V, W^a)$. It can be used to express the kinetic energy

$$ K = \frac{1}{2} (u, v, w) \mathbf{T} (u, v, w)^T, \qquad (13) $$

seen as an inner product between covariant components. The $\mathbf{T}$ matrix is the inverse of the metric tensor, it plays a central role in the discretization we propose.

At this stage we have set the framework. We have three different sets of velocity components available: Cartesian, contravariant and covariant. The Cartesian one is appealing but it does not correspond to either the coordinates lines or to the lines orthogonal to the coordinates hyperplanes. By contrast the contravariant and the covariant sets are natural. None of them is superior and they both have their interest. The covariant components are the natural prognostic variables of the momentum equations, while the contravariant components appear in all the fluxes and in particular into the incompressibility equation. For these reasons, $(u, v, w)$ will be called momentum and $(U, V, W^a)$ flux. In this approach flux and momentum are clearly two different objects. They handle their “horizontal” and “vertical” components very differently.

Finally because the two objects are related by the kinetic energy, their relation depends on whether the hydrostatic assumption is made or not.
203 Within the hydrostatic assumption, K simplifies into

$$K_{hydro} = \frac{1}{2} [U^2 + V^2]. \qquad (14)$$

This implies $w = 0$, which expresses the fact that the vertical momentum equation is a diagnostic equation. It also implies $(u, v) = (U, V) = (u_c, v_c)$, which means that the distinction between contravariant and covariant vanishes. Therefore everything can be cast in terms of Cartesian components and these subtleties can be ignored. The case of $W^a$ is special because it keeps its definition (4), albeit with $w = 0$. $W$, the flux through the moving grid, is called the omega velocity. The distinction between Cartesian, contravariant and covariant components is thus essentially a non-hydrostatic issue.

## 2.2. Equations in momentum variables

We now derive the momentum equations for the newly defined momentum components $(u, v, w)$, used as prognostic variables. As starting point, we will repeat equations (3.8 – 3.10) together with equations (3.13) and (3.16) from Kanarska et al. (2007) that define the momentum equations in terms of substantial derivatives of Cartesian velocity components, i.e. in flux form,

$$ \begin{aligned} \frac{D u_c}{Dt} &= -A_{\xi} \left. \frac{\partial(p^h + q)}{\partial \xi} \right|_z + \mathcal{V}(C_{u_c} + G_{u_c} + F_{u_c}) \\ \frac{D v_c}{Dt} &= -A_{\eta} \left. \frac{\partial(p^h + q)}{\partial \eta} \right|_z + \mathcal{V}(C_{v_c} + G_{v_c} + F_{v_c}) \\ \frac{D w_c}{Dt} &= -A_{\sigma} \frac{\partial q}{\partial \sigma} + \mathcal{V}(C_{w_c} + F_{w_c}) \end{aligned} \qquad (15) $$

where $(A_{\xi}, A_{\eta}, A_{\sigma}) = (h_{\eta}h_{\sigma}, h_{\xi}h_{\sigma}, h_{\xi}h_{\eta})$ are elementary surfaces and $\mathcal{V} = h_{\xi}h_{\eta}h_{\sigma}$ is the elementary volume. The gravitational term and the vertical
hydrostatic pressure gradient are absent in (15) because they balance. The two contributions $p^h$ and $q$ to the total pressure thus play a completely different role: $p^h$ accounts for the gravitational effects whereas $q$ enforces the incompressibility. The additional subscripts $c$ emphasize that the momenta are the Cartesian velocity components. $C_{u_c,v_c,w_c}$ represent the Coriolis force, $G_{u_c,v_c}$ represent the curvilinear terms arising from the curvature of the grid in the horizontal plane, and $F_{u_c,v_c,w_c}$ represent viscous dissipation, mixing and forcing terms. Note that there is no curvature term in the vertical direction and that $C_{w_c} = 0$ if the traditional assumption is made. The substantial derivative

$$ \frac{\mathcal{D}*}{\mathcal{D}t} = \frac{\partial}{\partial t}(\mathcal{V}*) + \frac{\partial}{\partial \xi}(\mathcal{A}_{\xi}U*) + \frac{\partial}{\partial \eta}(\mathcal{A}_{\eta}V*) + \frac{\partial}{\partial \sigma}(\mathcal{A}_{\sigma}W*) \quad (16) $$

uses the flux through the breathing mesh. Note that in the present formalism (16) is not impacted by the new definition of momentum because the flux defined in (3) is identical to the one in Kanarska et al. (2007), up to the area metric factors. Consequently the continuity equation remains identical

$$ \frac{\partial}{\partial t}\mathcal{V} + \frac{\partial}{\partial \xi}(\mathcal{A}_{\xi}U) + \frac{\partial}{\partial \eta}(\mathcal{A}_{\eta}V) + \frac{\partial}{\partial \sigma}(\mathcal{A}_{\sigma}W) = 0. \quad (17) $$

We can now recast the equations in terms of momentum. First we apply the total derivative (16) to (9)

$$ \begin{aligned} \frac{\mathcal{D}u}{\mathcal{D}t} &= \frac{\mathcal{D}u_c}{\mathcal{D}t} + s_{\xi} \frac{\mathcal{D}w_c}{\mathcal{D}t} + w_c \frac{\mathcal{D}s_{\xi}}{\mathcal{D}t} \\ \frac{\mathcal{D}v}{\mathcal{D}t} &= \frac{\mathcal{D}v_c}{\mathcal{D}t} + s_{\eta} \frac{\mathcal{D}w_c}{\mathcal{D}t} + w_c \frac{\mathcal{D}s_{\eta}}{\mathcal{D}t} \\ \frac{\mathcal{D}w}{\mathcal{D}t} &= \frac{\mathcal{D}w_c}{\mathcal{D}t}, \end{aligned} \quad (18) $$
237 then substitute the material derivatives of $(u_c, v_c, w_c)$ with (15) to get

$$
\begin{align}
\frac{\mathcal{D}u}{\mathcal{D}t} &= -\mathcal{A}_{\xi} \left. \frac{\partial(p^h + q)}{\partial\xi} \right|_z - \mathcal{A}_{\sigma} s_{\xi} \frac{\partial q}{\partial\sigma} + w \frac{\mathcal{D}s_{\xi}}{\mathcal{D}t} + \mathcal{V}(C_u + G_u + F_u) \nonumber \\
\frac{\mathcal{D}v}{\mathcal{D}t} &= -\mathcal{A}_{\eta} \left. \frac{\partial(p^h + q)}{\partial\eta} \right|_z - \mathcal{A}_{\sigma} s_{\eta} \frac{\partial q}{\partial\sigma} + w \frac{\mathcal{D}s_{\eta}}{\mathcal{D}t} + \mathcal{V}(C_v + G_v + F_v) \tag{19} \\
\frac{\mathcal{D}w}{\mathcal{D}t} &= -\mathcal{A}_{\sigma} \frac{\partial q}{\partial\sigma} + \mathcal{V}(C_w + F_w), \nonumber
\end{align}
$$

238 where $(F_u, F_v, F_w) = (F_{u_c} - s_{\xi}F_{w_c}, F_{v_c} - s_{\eta}F_{w_c}, F_{w_c})$ and $(G_u, G_v) = (G_{u_c}, G_{v_c})$.

239 The Coriolis terms $C_{u,v,w}$ are detailed in Appendix A. Using the chain rule

$$
\mathcal{A}_{\xi} \left. \frac{\partial *}{\partial \xi} \right|_{z} = \mathcal{A}_{\xi} \frac{\partial *}{\partial \xi} - \mathcal{A}_{\sigma} s_{\xi} \frac{\partial *}{\partial \sigma} \quad (20)
$$

240 simplifies the expression for the non-hydrostatic pressure gradient

$$
\begin{align}
\frac{\mathcal{D}u}{\mathcal{D}t} &= -\mathcal{A}_{\xi} \left. \frac{\partial p^h}{\partial \xi} \right|_z - \mathcal{A}_{\xi} \frac{\partial q}{\partial \xi} + w \frac{\mathcal{D}s_{\xi}}{\mathcal{D}t} + \mathcal{V}(C_u + G_u + F_u) \nonumber \\
\frac{\mathcal{D}v}{\mathcal{D}t} &= -\mathcal{A}_{\eta} \left. \frac{\partial p^h}{\partial \eta} \right|_z - \mathcal{A}_{\eta} \frac{\partial q}{\partial \eta} + w \frac{\mathcal{D}s_{\eta}}{\mathcal{D}t} + \mathcal{V}(C_v + G_v + F_v) \tag{21} \\
\frac{\mathcal{D}w}{\mathcal{D}t} &= -\mathcal{A}_{\sigma} \frac{\partial q}{\partial \sigma} + \mathcal{V}(C_w + F_w). \nonumber
\end{align}
$$

241 (21) is the new system of momentum equations we propose to use to integrate
242 the primitive equations in sigma-coordinates model when the hydrostatic
243 assumption is relaxed. The terms in (21) involving the substantial derivatives
244 of $s_{\xi}$ and $s_{\eta}$ can also be written as

$$
\begin{align}
w \frac{\mathcal{D} s_{\xi}}{\mathcal{D} t} &= w \mathcal{V} \frac{\partial s_{\xi}}{\partial t} + w U \mathcal{A}_{\xi} \frac{\partial s_{\xi}}{\partial \xi} + w V \mathcal{A}_{\eta} \frac{\partial s_{\xi}}{\partial \eta} + w W \mathcal{A}_{\sigma} \frac{\partial s_{\xi}}{\partial \sigma} \tag{22} \\
w \frac{\mathcal{D} s_{\eta}}{\mathcal{D} t} &= w \mathcal{V} \frac{\partial s_{\eta}}{\partial t} + w U \mathcal{A}_{\xi} \frac{\partial s_{\eta}}{\partial \xi} + w V \mathcal{A}_{\eta} \frac{\partial s_{\eta}}{\partial \eta} + w W \mathcal{A}_{\sigma} \frac{\partial s_{\eta}}{\partial \sigma} . \nonumber
\end{align}
$$

245 This derivation makes (22) inherently new terms compared to Kanarska et al.

246 (2007). The system (21) can be derived independently by application of
247 Hamilton's principle (Appendix C). It yields the vector-invariant form.
248 After a bit of algebra the flux form can be recovered along with the expression
249 for the curvilinear terms. This alternative derivation blends the contributions
250 of (22) and $\nabla G_{u,v}$ into a same term. Finally note that the chain rule (20)
251 can also be applied to the hydrostatic pressure gradient term in (21). This
252 introduces a buoyancy term in the horizontal components of the momentum
253 equations, like in the variational approach (Appendix C). In practice the
254 hydrostatic pressure gradient and this buoyancy term are computed together
255 to limit the pressure gradient error (Shchepetkin and McWilliams, 2003).

256 Finally we want to emphasize that there is no necessity to prognose the
257 covariant components. A perfectly valid choice could be to prognose the
258 contravariant ones ($U, V, W^a$). In fact, it seems that the good choice of prog-
259 nostic variables depend on the formulation of the momentum equations: in
260 vector-invariant, the covariant components appear to be natural (see Eq C.4),
261 in flux-form, the contravariant components are more natural.

### 3. Discrete operators for the non-hydrostatic pressure

263 In this section we first show how the Poisson equation for pressure can be
264 derived from the definitions of the kinetic energy and the divergence operator.
265 We then discretize the kinetic energy and deduce the discretizations for the
266 $T$ matrix and the Poisson equation.

#### 3.1. *Pressure projection in matrix formulation*

268 Following Kanarska et al. (2007), we will use the language of pressure
269 projection methods (Chorin, 1968; Bell et al., 1989). These methods refer to
270 a class of fractional step, implicit algorithms for unsteady flows that satisfy a
non-divergence constraint. In the fractional step method a provisional velocity field is constructed first by advancing the momentum equations in time without enforcing the no-divergence constraint. Subsequently, the divergent part of velocity field is removed by a projection to the no-divergence subspace. Here, we will present a derivation following Molemaker et al. (2005) to arrive, first at an expression for the operators in matrix form and second, specifically in terms of discrete coefficients.

We define the diagonal matrices $L$, $A$, $V$ corresponding to grid lengths, cell face areas and cell volumes working on velocity points. The elements of $L$ are the Lamé coefficients. The matrices are related by $V = AL$. We define $d$ the differentiation matrix acting on cell centers (rho-points) and returning gradient components at face centers (velocity points). $d$ carries no dimension. The divergence is $D = -d^T A$. Note the minus and the transpose. The divergence works on fluxes and returns finite volume scalars (scalar times the elementary volume).

Using (12), we can write the relation between fluxes and momenta compactly in matrix form as:

$$X^a = T x \qquad (23)$$

where $x = (u, v, w)^T$ is the solution vector in terms of momenta and $X^a = (U, V, W^a)^T$ the equivalent in terms of fluxes. We use the domain integrated kinetic energy $K = \frac{1}{2}x^T V X^a$ to define the inner product of two momenta $x$ and $y$

$$<x, y> = x^T V T y. \qquad (24)$$

To be an inner product $VT$ needs to be symmetric. This is achieved by
discretizing the kinetic energy in such way that

$$ \mathcal{K} = \frac{1}{2} (\mathbf{L}x)^T \tilde{\mathbf{T}} (\mathbf{L}x), \qquad (25) $$

where $\tilde{\mathbf{T}}$ is a symmetric matrix. The discretization of $\mathbf{T}$ is then deduced

$$ \mathbf{T} = \mathbf{A}^{-1} \tilde{\mathbf{T}} \mathbf{L}. \qquad (26) $$

The non-divergence condition is a linear constraint of the flow that, using (17), can now be written in terms of the momenta as:

$$ \mathbf{D}\mathbf{T} x_D = 0 \qquad (27) $$

where the subscript $D$ indicates a solution vector that satisfies the no-divergence constraint (27).

We are now sufficiently equipped to project from a provisional velocity field $x$, to a non-divergent field $x_D$. If $\epsilon$ is the difference vector between the provisional solution and its projection onto the no-divergence divergence subspace:

$$ \epsilon = x - x_D, \qquad (28) $$

we demand that this correction vector is orthogonal to the subspace of no-divergence solutions, or:

$$ \langle x_D, \epsilon \rangle = 0. \qquad (29) $$

The correction is of the form

$$ \epsilon = \mathbf{V}^{-1} \mathbf{D}^T p, \qquad (30) $$

with $p$ a scalar field, because

$$ \langle x_D, \mathbf{V}^{-1} \mathbf{D}^T p \rangle = x_D^T \mathbf{T}^T \mathbf{D}^T p = (\mathbf{D} \mathbf{T} x_D)^T p = 0. \qquad (31) $$
$p$ acts as a Lagrange multiplier. Its function is to help satisfy the linear constraint of non-divergence. Physically speaking, the pressure functions as such for an incompressible flow, and we see that the correction vector $\epsilon$ is, in fact, the pressure gradient. We now have the following system:

$$VTx - VTx_D = (DT)^T p \quad (32)$$

$$DT x_D = 0. \quad (33)$$

From which we can easily derive an equation for $p$

$$(DV^{-1}T^T D^T) p = DT x. \quad (34)$$

(34) is the Poisson equation for pressure and

$$P = DV^{-1}T^T D^T \quad (35)$$

is the Poisson matrix. The matrix can be simplified further into

$$P = d^T \tilde{T} d, \quad (36)$$

which is symmetric if $\tilde{T}$ is. We have done all this effort to arrive at the Poisson equation for the pressure (34). We see that once the kinetic energy is discretized in a symmetric form using the $\tilde{T}$ matrix, the divergence operator $D$, and the volume matrix $V$ are fixed, the definition of the other operators follow. The flux-momentum approach yields a compact 15-point stencil because $P$ can be cast into the simplified form (36) for which the differential operator $d$ has a minimal form.

### 3.2. Variable staggering

We describe here how the momentum and flux can be discretized in such way that the operators grad and div operating on them acquire the simplest
possible discretization, namely basic finite differences. For that we use the
framework of discrete differential geometry (DDG) (Desbrun et al., 2008;
Grady and Polimeni, 2010), which extends the common concepts of finite
volume and C-grid staggering to a higher degree of generality. The frame-
work also allows to rationalize why the contravariant components should be
interpreted as normal to coordinate hypersurfaces and not as aligned with
the coordinate lines. The interpretation we present is also completely aligned
with Thuburn and Cotter (2012) for non-orthogonal grids in the context of
shallow water equations.

The goal of this section is to introduce all the elements present in figure 1,
adapted from Desbrun et al. (2008), which summarizes the various aspects of
the spatial discretization. The (ξ, η, σ) coordinates are discretized in logically
rectangular cells, with indices (i, j, k) pointing to cell vertices. The cell edges
are along coordinates lines, cell faces are along coordinates hyperplanes and
finally each cell form a volume. This defines the primal grid. The dual grid
is defined by setting vertices at the primal cell centers, then dual grid edges
connect these dual grid vertices through primal faces. Dual faces connect
dual edges. Finally the dual volume cell is defined as the interior of dual
faces, it surrounds each primal vertex.

The deep breakthrough of DDG is to see the usual operators (gradient,
curl and div) as transformations between objects: the gradient transforms
vertices to edges, the curl edges to faces and the divergence faces to volume.
Each operator is then simply a finite difference between quantities. We detail
this below.

Scalar quantities are naturally discretized at either cell centers or cell
<figure><img src="image_2.png" /><figcaption>Figure 1: Staggering of variables on the primal and dual grids according to their nature: velocity, pressure or tracer. Contravariant components are defined on the primal, covariant components on the dual. They are related by the T transformation. By interpreting the contravariant as flux, living on faces, and the covariant as gradient terms, living on the edges, this makes the pressure gradient operation and the divergence of the flux very simple. This interpretation is consistent with a discretization of the pressure at the primal cell centers (dual vertices) and a finite volume interpretation of tracers (primal volume). Each object is connected to its right-neighbour (primal) and left-neighbour (dual) by a differentiation operator which, depending on the object, is the gradient, the curl or the divergence. The vectorial interpretation is illustrated by the cyan dashed box. The model's objects and operators (bold) are delimited by the black dashed box.</figcaption></figure>

volumes depending on whether the quantity is a true scalar, like pressure, or an integrated quantity, like tracer content. Vectors are naturally discretized at either cell edges or cell faces. Contravariant components are defined on the primal whereas the covariant components are on the dual. The two interpretations for the contravariant components can now be enriched by this geometric framework: the contravariant components are either i) on the faces of the primal or ii) on the edges of the primal, which is the vectorial inter-

345

346

347

348

349

350

351
pretation. Interpretation i) is the one we made, it makes the flux naturally
transformed into a volume via the divergence operator. With interpretation
ii), the divergence of the contravariant components is not a trivial operation.
Similarly, the momentum should behave like a gradient, because it is covari-
ant. It can be discretized either i) on the edges of the dual or on ii) the
faces. To be consistent, it has to be one edges. It makes the computation of
the pressure gradient a simple finite difference. By using discretization ii),
the expression of the pressure gradient acting on the momentum would be
complicated.

The consistency of i) with the C-grid discretization practice is even deeper:
the finite volume discretization of tracers implies that the flux sits at the pri-
mal faces, and the discretization of pressure at cell centers implies that the
covariant velocity sits at the dual edges. Finally, having the contravariant
components on the primal faces makes the operation div **U** covariant, i.e.
coordinates independent.

The **T** matrix transforms objects from the dual edges to the primal faces.
Lastly, from figure 1 it is clear that the vorticity, defined as the curl of
the momentum lives at the dual faces. It inherits covariant components.
The vorticity appears explicitly in the vector-invariant form. It is only a
diagnostic quantity in flux-form.

From now on, each variable is indexed with a triplet of indices (i, j, k).
Faces and edges with index (i, j, k) are located to the west, south and bottom
of the cell (i, j, k). To lighten the notation the indices will be indicated only
when they differ from i, j or k, e.g. φᵢ₋₁ will stand for φᵢ₋₁,ⱼ,ₖ. There is
nᵣ cells in the vertical and nᵣ + 1 w-levels. We define the right difference
operator $\delta_{i,r}$ acting on the field $\phi_{i,j,k}$ in the direction $\xi$ as $\delta_{i,r}[\phi] = \phi_{i+1} - \phi$ and the left difference operator $\delta_{i,l}[\phi] = \phi - \phi_{i-1}$. We define the right averaging operator $\bar{\phi}^{i,r} = (\phi_{i+1} + \phi)/2$ and the left averaging one $\bar{\phi}^{i,l} = (\phi + \phi_{i-1})/2$. Likewise we define the same operators acting in direction $\eta$ with index $j$, and in direction $s$ with index $k$.

Grid spacing $(h_\xi, h_\eta, h_\sigma)$, and local slopes $(s_\xi, s_\eta)$ are also defined at cell centers with

$$s_{\xi} = \left[ \frac{\delta_{i,l}[z]}{h_{\xi}^{u}} \right]^{i,r} \quad \text{and} \quad s_{\eta} = \left[ \frac{\delta_{j,l}[z]}{h_{\eta}^{v}} \right]^{j,r}. \qquad (37)$$

A superscript of $u, v, w$ denotes a quantity that has been averaged to the position of its respective velocity, e.g. $h_\xi^u = \delta_{i,r}[h_\xi]$. In particular $h_\sigma^w$ is only half of its value at the bottom ($k=1$) and top level ($k=n_z+1$). Faces areas are defined at velocity points using product of two lengths, e.g. $A_\xi^u = h_\eta^u h_\sigma^u$. Volumes are defined at velocity points using the product of the lengths in each direction, e.g. $V^u = h_\xi^u h_\eta^u h_\sigma^u$.

### 3.3. Discretization of the gradient and divergence

As announced, the discretization of the gradient and divergence boil down to simple finite differences. If we denote $X = (U, V, W)^T$ the flux vector and $\delta = DX$ the result of the action of $D$, then the discretization of the divergence reads

$$\delta_{i,r}[A_{\xi}^{u}U] + \delta_{j,r}[A_{\eta}^{v}V] + \delta_{k,r}[A_{\sigma}^{w}W] = \delta. \qquad (38)$$

The definition is valid at all grid points, including at the boundaries. The bottom boundary condition is handled outside of the divergence by simply setting $W=0$. The gradient is defined as the skew-adjoint of $D$. For interior
398 points it is easy to show that

$$h_i^u u_i = \delta_{i,l}[q], \quad h_i^v v_i = \delta_{j,l}[q], \quad h_i^w w_i = \delta_{k,l}[q], \qquad (39)$$

where $q$ is a scalar defined on the dual vertices and $(u, v, w)$ are the (covariant) gradient components. The expression for the curl is given in Appendix C. It is basically computed by estimating the circulation around each face.

### 3.4. Discretization of the kinetic energy

As said above, the discretization of the $\mathbf{T}$ matrix derives from $\mathcal{K}$ written in the form (25). To discretize $\mathcal{K}$ we start from the continuous expression

$$\mathcal{K} = \iiint \left[ \frac{1}{2}(u^2 + v^2 + \gamma w^2) - (s_{\xi}uw + s_{\eta}vw) \right] \mathcal{V} d\xi d\eta d\sigma, \quad (40)$$

with $\gamma = 1 + s_{\xi}^2 + s_{\eta}^2$, and split it into several contributions

$$\mathcal{K} = \mathcal{K}_{\text{int}} + \mathcal{K}_{\text{cross}} + \mathcal{K}_{\text{bottom}}, \qquad (41)$$

where $\mathcal{K}_{\text{bottom}}$ accounts for all terms involving the bottom vertical momentum $w_b$ (w at $k=1$), $\mathcal{K}_{\text{cross}}$ accounts for the cross terms $uw$ and $vw$, and $\mathcal{K}_{\text{int}}$ is the remainder. The interior kinetic energy

$$\mathcal{K}_{\text{int}} = \frac{1}{2} \sum_{i,j} \left\{ \sum_{k=1}^{n_z} u^2 \mathcal{V}^u + \sum_{k=1}^{n_z} v^2 \mathcal{V}^v + \sum_{k=2}^{n_z+1} \gamma^w w^2 \mathcal{V}^w \right\} \qquad (42)$$

is the simplest term. Each component is weighted with its volume defined at the component location. This part of the kinetic energy encodes the diagonal terms of the $\tilde{\mathbf{T}}$ matrix, which are $(A^u/h_\xi^u, A^v/h_\eta^v, c^w A^w/h_\sigma^w)$. The cross terms
contribution reads

$$
\begin{equation}
\begin{split}
K_{cross} = -\frac{1}{2} \sum_{i,j} & \left[ s_{\xi} h_{\eta} \overline{[h_{\xi}^{u} u]}^{i,r} + s_{\eta} h_{\xi} \overline{[h_{\eta}^{v} v]}^{j,r} \right]_{k=1} [h_{\sigma}^{w} w]_{k=2} \\
& - \sum_{i,j} \sum_{k=2}^{n_z-1} \left[ s_{\xi} h_{\eta} \overline{[h_{\xi}^{u} u]}^{i,r} + s_{\eta} h_{\xi} \overline{[h_{\eta}^{v} v]}^{j,r} \right] \overline{[h_{\sigma}^{w} w]}^{w,r} \\
& - \frac{1}{2} \sum_{i,j} \left[ s_{\xi} h_{\eta} \overline{[h_{\xi}^{u} u]}^{i,r} + s_{\eta} h_{\xi} \overline{[h_{\eta}^{v} v]}^{j,r} \right]_{k=nz} [h_{\sigma}^{w} w]_{k=nz} \\
& - \sum_{i,j} \left[ s_{\xi} h_{\eta} \overline{[h_{\xi}^{u} u]}^{i,r} + s_{\eta} h_{\xi} \overline{[h_{\eta}^{v} v]}^{j,r} \right]_{k=nz} \overline{[h_{\sigma}^{w} w]}_{k=nz+1}.
\end{split}
\tag{43}
\end{equation}
$$

All terms are defined at cell centers (rho points). The volume elements have been carefully split into three length elements to form elementary circulations to ensure (25). The choice to average the terms at cell centers allows for the maximum of terms cancellation in the resulting pressure matrix. The sum has been split in horizontal slabs: in the bottom level ($k = 1$) we have removed the part involving $w_b$, in the top level ($k = n_z$) we have introduced an extra factor 2 to compensate for $h_\sigma^w$ being halved.

409

410

411

412

413

414

415

The bottom contribution accounts for the kinetic energy density $K = \frac{1}{2}\gamma w^2 - (s_\xi u w + s_\eta v w)$ integrated over the bottom cells lower half volume. The $u^2$ and $v^2$ are already accounted in $K_{int}$. Because $w$ at the bottom is not a prognostic variable, it should be substituted with its diagnostic relation $\gamma w = s_\xi u + s_\eta v$, obtained from (11), the definition of $W^a$, and the fact that the breathing vanishes at the bottom boundary ($\partial_t z = 0$). After substitution, the kinetic energy density reads $K = -\frac{1}{2\gamma}[(s_\xi u)^2 + (s_\eta v)^2 + 2s_\xi s_\eta uv]$, which is the continuous density we use for the discretization at the bottom. The
discretized volume-integrated bottom kinetic energy then reads

$$ K_{bottom} = -\frac{1}{4} \sum_{i,j,k=1} \left\{ \left[ \frac{s_{\xi}^2}{\gamma} \right]^{i,l} u^2 \mathcal{V}^u + \left[ \frac{s_{\eta}^2}{\gamma} \right]^{j,l} v^2 \mathcal{V}^v + 2 \frac{s_{\xi}s_{\eta}h_{\sigma}}{\gamma} [\overline{h_{\xi}^u u}]^{i,r} [\overline{h_{\eta}^v v}]^{j,r} \right\}, \quad (44) $$

where the $u^2$ and $v^2$ have been discretized at their velocity point and the $uv$ cross term has been averaged at cell center. The whole term is indeed weighted with a half volume cell because of the extra factor 1/2 that compensates for $\mathcal{V}^u$, $\mathcal{V}^v$ and $h_{\sigma}$ being full cell metrics.

### 3.5. Discretization of the fluxes

The flux discretization, that is the action of $\mathcal{T}$ on the momentum, can now be deduced either upon identification of the $\tilde{\mathcal{T}}$ with (25) and then application of (26) or by differentiating $\mathcal{K}$ with respect to the momentum, e.g.

$$ U_{i,j,k} = \frac{1}{\mathcal{A}_{\xi}^{u}} \frac{\partial \mathcal{K}}{\partial(h_{\xi}^{u} u_{i,j,k})}. \quad (45) $$

We have for the horizontal flux $U$

* In the interior $(1 < k < n_z)$

$$ U = u - \frac{1}{\mathcal{A}_{\xi}^{u}} \overline{s_{\xi} h_{\eta} [\overline{h_{\sigma}^{w} w}]^{w, r, i, l}} \quad (46) $$

* At the bottom $k=1$

$$ U = \left( 1 - \left[ \frac{s_{\xi}^2}{2\gamma} \right]^{i,l} \right) u - \frac{1}{2\mathcal{A}_{\xi}^u} \left( \frac{s_{\xi}s_{\eta}h_{\sigma}}{\gamma} \overline{[h_{\eta}^v v]^{j,r}} + \frac{s_{\xi}h_{\eta}}{\sigma} [h_{\sigma}^w w]_{k=2}^{i,l} \right) \quad (47) $$
427 • At the top $k = n_z$

$$U = u - \frac{1}{2\mathcal{A}_{\xi}^{u}} \left( \overline{s_{\xi}h_{\eta}[h_{\sigma}^{w}w]_{k=n_z}} + 2s_{\xi}h_{\eta}[h_{\sigma}^{w}w]_{k=n_z+1} \right)^{i,l} \quad (48)$$

428 The discretization for $V$ is similar. It is obtained by applying the permutation $(\xi, u, i) \leftrightarrow (\eta, v, j)$ on the above expressions.

429 The discretization of the vertical flux $W^a$ reads

431 • at the bottom $k = 1$

$$W^a = 0, \quad (49)$$

432 • in the interior $1 < k < n_z + 1$

$$W^a = \gamma^w w - \frac{1}{\mathcal{A}_{\sigma}^w} \left( \overline{s_{\xi}h_{\eta}[h_{\xi}^u u]^{i,r}} + \overline{s_{\eta}h_{\xi}[h_{\eta}^v v]^{j,r}} \right), \quad (50)$$

433 • and at the surface $k = n_z + 1$

$$W^a = \gamma^w w - \frac{1}{\mathcal{A}_{\sigma}^w} \left[ \overline{s_{\xi}h_{\eta}[h_{\xi}^u u]^{i,r}} + \overline{s_{\eta}h_{\xi}[h_{\eta}^v v]^{j,r}} \right]_{k=nz} . \quad (51)$$

434 These expressions complete the definition for the discrete $T$ matrix, that transforms momentum into flux.

### 3.6. Discretization of the Poisson operator

435 The Poisson matrix is obtained from (35) with the use of the discretization of the fluxes and the divergence. Its discretization is given in appendix Appendix B. Incidentally the compactness of the stencil does not change the matrix condition number, compared to any other discretization. Indeed, the condition number is controlled by the grid characteristics: grid steps and number of grid points.
The combination of the discrete divergence and the discrete kinetic energy fully determine the boundary conditions on the pressure equation. At the surface the boundary condition is of Dirichlet type ($p = 0$) whereas at the bottom and laterally it is of Neumann type. To exhibit them we use the property that the flux is divergent-free $\langle p, \mathbf{D}\mathbf{T}x_D \rangle = 0$, then make use of the adjoint to find $\langle \mathbf{D}^T p, \mathbf{T}x_D \rangle = 0$ which expresses that the pressure force does not work. By expanding this inner product as a sum on velocity points we recover the discrete kinetic energy for which one of the momentum component is replaced with the pressure gradient. In particular the $w^2$ term at $(i, j, n_z + 1)$ in $K_{int}$ is now $-\gamma^w \mathcal{V}^w w p_{i,j,n}/h_{i,j,n_z+1}^w$, from which we deduce the vertical pressure gradient at the surface is $\partial_z p = -p_{i,j,n}/h_{i,j,n_z+1}^w$. This discretization is indeed a vertical gradient if $p_{i,j,n_z+1/2} = 0$, thus proving that the surface boundary condition is of Dirichlet type. The surface boundary condition is realized in the Poisson stencil by the $A_\sigma^w/h_\sigma^w$ term at $k = n_z+1$ in the diagonal coefficient $P_1$ (see (B.2) in Appendix B). The bottom boundary condition is handled by $K_{bot}$ and behaves such that $W = 0$ at the bottom. The pressure gradient at the bottom is not defined but its associated vertical flux ($\mathbf{TD}^T p$) is zero there, which corresponds to a Neumann boundary condition. The lateral boundary conditions is treated by assuming the covariant velocity, either $u$ or $v$, vanishes. This removes in the Poisson matrix all connections between interior and outer pressure points.

## 4. Discussion

We now assess the properties of the discretization on a test-case. The goal is to test the discretized projection in a configuration simple enough, yet
ensuring that none of the 15 coefficients of the Poisson matrix vanishes. The
comparison with an analytical solution would have been ideal but we have
found none with a genuinely three-dimensional topography.

The test-case, expressed in Cartesian coordinates, consists in projecting a
470
divergent velocity field (u, v, w) = (0, 0, cos(π x / L)) onto the divergent-free
471
subspace in the closed domain (x, y, z) ∈ [0, L] × [0, L] × [-H + h(x, y), 0],
472
where h(x, y) = h₀ exp(-α (x-L/2)² - α (y-L/2)²) is a sea-mount topogra-
473
phy. The numerical values are L = 10 km, H = 4 km, h₀ = 2 km, α = 25/L².
474
The test-case is particularly severe for the bottom boundary condition be-
475
cause the initial flow field is already divergence-free but it is far from being
476
the solution because it does not satisfy the no-flow bottom boundary condi-
477
tion. The right hand side for the Poisson equation is thus zero everywhere
478
except along the topography.
479

Even though the test-case is relatively simple, it has no analytical solu-
tion. It is therefore difficult to prove that the numerical solution obtained
with the 15-point stencil is correct. To circumvent this issue we have com-
puted the solution using two other methods: the tall 25-point stencil obtained
with Cartesian components and the 7 points stencil derived using a z-level dis-
cretization. The 25-point stencil matrix can be written as P<sub>w</sub> = D<sub>w</sub>V<sup>-1</sup>D<sub>w</sub><sup>T</sup>,
where the divergence operator D<sub>w</sub> has a wider stencil than D because of
the extra derivatives coming from the chain rule. The approach closely fol-
lows Auclair et al. (2011) but it differs by explicitly obtaining a symmetric
matrix. In the z-levels discretization, the topography is handled via a mask
array that flags cells outside the domain. The resulting matrix is similar
to what is implemented in MITgcm (Marshall et al., 1997). In both cases,
480
481
482
483
484
485
486
487
488
489
490
491
the bottom boundary condition is trivially satisfied by resorting a divergence matrix that does not uses the *w* of the bottom level to compute the divergence in the first level. In addition to these three matrices we also test the impact of the extra-care needed to handle to bottom boundary condition in the 15-point stencil case by implementing a degraded version of it in which we simply remove in the **T** matrix all the contributions from the $K_{\text{bot}}$ term. This matrix is obtained by setting the flag $\mu = 0$ in the Poisson stencil given in appendix.

We solve the problem on a series of grids $(n_x, n_y, n_z) = 2^n \times (25, 25, 10)$ with doubling resolution. The matrices are implemented in Matlab. Because the grids need to fit in RAM of a single chipset, we were able to test the range $n = 0, 1, 2, 3$ for the four matrices. The solution is a non trivial flow. Figure 2 shows the Cartesian $u_c$ and $w_c$ components of the flow in the mid-plane $y = L/2$ for $n = 3$. The flow has on one side a down-welling branch emanating from the surface and an upwelling branch on the other. In between, the flow is passing around (not shown) and over the sea-mount. The four implementations give an identical flow. At this resolution the differences are several orders of magnitude smaller than the solution. To quantitatively measure the convergence we compute the kinetic energy $K$ of each solution. Results are shown in Figure 3. They clearly converge to the same value $K_\infty$, indicated by the tick black line, albeit with a different convergence rate. The 15- and the 25-point stencils converge more rapidly to $K_\infty$. Specifically their increment of accuracy $\delta K = K(n) - K(n-1)$ decreases by a factor of 4 when the resolution is doubled, implying that $K$ is of the form $K = K_\infty + c 4^{-n}$, where $c$ is a constant. We can then use the Richardson extrapolation on the
517 15-point results to estimate this limit value as $K_\infty \approx K(n) + \delta K(n)/3$.

518 This allows to monitor the absolute error $|K - K_\infty|$ (inset of Fig. 3). This
519 figure confirms that the 15- and 25-point stencil converge to the solution at
520 second order, as indicated by the -2 slope. Interestingly, the 15-point sten-
521 cil is slightly more accurate than the 25-point stencil. This can be traced
522 back to slightly smoother vertical shears near the bottom. The z-levels only
523 converges at first order, which is consistent with the discretization of the to-
524 pography that is only first order. For that case the solution has the footprint
525 of steps near the topography (not shown). The degraded 15-point stencil
526 also has a global order 1 convergence. Comparisons with the full 15-point
527 stencil reveal that the horizontal flux vertical profile is discontinuous at the
528 bottom level, generating a maximum of error there. Because the problem
529 is essentially global this impacts the solution everywhere. The contribution
530 of the horizontal cross terms in the bottom level ensures a uniform second
531 order estimate of the velocities everywhere.

532 **5. Conclusion**

533 We have presented an alternative system of equations for the primitive
534 equations written in sigma coordinates with a breathing mesh. It is based on
535 the recognition of two sets of components for the velocity: a contravariant
536 set, involved in the flux definitions, and a covariant one, involved in the
537 momentum definition. The two sets are related by a linear transformation,
538 the metric tensor, fully determined by the kinetic energy. We provide the
539 momentum equations written for the covariant components in both vector-
540 invariant and flux-form. We show that within the hydrostatic assumption,
<figure><img src="image_3.png" /><figcaption>Figure 2: Vertical section of the velocity Cartesian components $w_c$ (a) and $u_c$ (b) at $y = L/2$, for the full solution obtained with the 15-point stencil matrix.</figcaption></figure>

the distinction between contra- and covariant vanishes, which explains why the distinction has been largely ignored so far and why it is genuinely a non-hydrostatic issue. In flux-form the new system of momentum equation possesses extra pseudo-force-like terms in the horizontal. They involve the partial derivatives of the sigma slopes. The chief gain of this system, and the reason we present it, is to remove chain rule operations in the gradient and divergence operators. In the NH case this allows, in turn, to express the Poisson equation for the pressure very compactly. We show how this property can be implemented numerically yielding to a compact symmetric stencil for the Poisson matrix with only 15 points. The discretization is done at the level
<figure><img src="image_4.png" /><figcaption>Figure 3: Kinetic energy as a function of resolution *n* for four different Poisson matrices: the 15-point stencil **P** matrix (blue), the 25 points Cartesian matrix **P**<sub>w</sub> (red), the 15-point matrix with the bottom boundary condition removed (yellow), the 7 points z-levels matrix (purple). The thick line is the limit value *K*<sub>∞</sub> obtained by Richardson extrapolation. The inset shows the absolute error |*K* − *K*<sub>∞</sub>| as a function of resolution.</figcaption></figure>

of the kinetic energy, which ensures that the pressure work vanishes globally.
It also ensures the symmetry of the stencil. We show how the no-flow bottom
boundary condition can be discretized with a second order accuracy. This
discretization is also done directly at the level of the kinetic energy.

This approach opens a promising way of solving NH physics in a sigma
coordinates model. The full implementation of this system of equations in
an ocean model will be reported in a forthcoming paper.
## 558 Appendix A. Coriolis terms

In the non-traditional approximation the Coriolis terms involve the three components of the Earth angular speed $\Omega$, namely

$$ 
\begin{aligned}
f_{\xi} &= 2\Omega \cos \phi \sin \alpha \\
f_{\eta} &= 2\Omega \cos \phi \cos \alpha \\
f_{z} &= 2\Omega \sin \phi,
\end{aligned} 
\quad (A.1) $$

where $\alpha$ is the angle between the $\xi$ and the zonal direction and $\phi$ is the latitude. From Appendix B of Kanarska et al. (2007) the Coriolis terms for Cartesian velocity components are

$$ 
\begin{aligned}
C_{u_c} &= -v_c f_z + w_c f_{\eta} \\
C_{v_c} &= +u_c f_z - w_c f_{\xi} \\
C_{w_c} &= -u_c f_{\eta} + v_c f_{\xi}.
\end{aligned} 
\quad (A.2) $$

Using the relations (9) yields the momentum based expression

$$ 
\begin{aligned}
C_u &= C_{u_c} + s_{\xi} C_{w_c} \\
C_v &= C_{v_c} + s_{\eta} C_{w_c} \\
C_w &= C_{w_c}.
\end{aligned} 
\quad (A.3) $$

Rearranging using (3) and (9) again, then defining

$$ f_{\sigma} = f_z - s_{\xi} f_{\xi} - s_{\eta} f_{\eta}, \quad (A.4) $$

yields

$$ 
\begin{aligned}
C_u &= -V f_{\sigma} + W^a f_{\eta} \\
C_v &= +U f_{\sigma} - W^a f_{\xi} \\
C_w &= -U f_{\eta} + V f_{\xi}.
\end{aligned} 
\quad (A.5) $$
Note that the Coriolis force involves the vertical flux $W^a$ relative to a fix grid, not $W$. These terms are obviously energy neutral because the energy equation is obtained by forming the inner product of the momentum equations with the flux $(U, V, W^a)$. The traditional approximation case is obtained by setting $f_\xi = f_\eta = 0$.

## 572 Appendix B. Poisson matrix

The discretization of the Poisson matrix follows from (35). Its action on a discrete scalar field $q_{i,j,k}$ can be written as

$$ 
\begin{aligned}
& P_1 q + P_2 q_{k-1} + P_{2,k+1} q_{k+1} + P_3 q_{j-1,k+1} + P_{3,j+1,k-1} q_{j+1,k-1} \\
& \qquad + P_4 q_{j-1} + P_{4,j+1} q_{j+1} + P_5 q_{j-1,k-1} + P_{5,j+1,k+1} q_{j+1,k+1} \\
& \qquad + P_6 q_{i-1,k+1} + P_{6,i+1,k-1} q_{i+1,k-1} + P_7 q_{i-1} + P_{7,i+1} q_{i+1} \\
& \qquad + P_8 q_{i-1,k-1} + P_{8,i+1,k+1} q_{i+1,k+1} + P_9 q_{i-1,j-1} + P_{9,i+1,j+1} q_{i+1,j+1} \\
& \qquad + P_{10} q_{i-1,j+1} + P_{10,i+1,j-1} q_{i+1,j-1} = r_{i,j,k} \quad (B.1)
\end{aligned}
$$

where the coefficients $P_l, l = 1, 2, \dots, 10$, define the stencil of the matrix and $r_{i,j,k}$ is the result. Coefficients pointing to indices outside of the index range are set to 0. To lighten the notation the indices are indicated only when they differ from $i, j$ or $k$, e.g. $P_2$ stands for $P_{2,i,j,k}$. Each coefficient has a general expression for interior levels $1 < k < n_z$ but differ at the top and bottom levels. We define the indicatrix function $\Phi_{l,i,j,k}$ that is 0 for $k \neq l$ and 1 for $k = l$. We also define a flag $\mu = 1$ for the full 15-point Poisson matrix and $\mu = 0$ for the degraded 15-point matrix in which the bottom boundary
581 condition is not explicitly handled. The diagonal coefficient reads

$$P_1 = \frac{A_{\sigma}^{w}}{h_{\sigma}^{w}}(1 - \Phi_1) + \left. \frac{A_{\sigma}^{w}}{h_{\sigma}^{w}} \right|_{k+1} + \left. \frac{A_{\eta}^{v}}{h_{\eta}^{v}} + \frac{A_{\eta}^{v}}{h_{\eta}^{v}} \right|_{j+1} + \left. \frac{A_{\xi}^{u}}{h_{\xi}^{u}} + \frac{A_{\xi}^{u}}{h_{\xi}^{u}} \right|_{i+1} . \quad (B.2)$$

Recall that $A_{\sigma}^{w}/h_{\sigma}^{w}$ is well defined at $k = n_z + 1$. The other coefficients read

$$P_2 = \frac{A_{\sigma}^{w}}{h_{\sigma}^{w}}\gamma^{w}(1 - \Phi_1) \qquad (B.3)$$

$$P_3 = \frac{1}{4}(s_{\eta}h_{\xi} + s_{\eta,j-1}h_{\xi,j-1})(1 - \Phi_{nz}) \qquad (B.4)$$

$$P_4 = \frac{A_{\eta}^{v}}{h_{\eta}^{v}}\left(1 - \frac{1}{4}\left[\frac{s_{\eta}^{2}}{\gamma} + \frac{s_{\eta,j-1}^{2}}{\gamma_{j-1}}\right]\Phi_{1}\mu\right) + \frac{1}{4}(s_{\eta}h_{\xi} - s_{\eta,j-1}h_{\xi,j-1})(\Phi_{nz} + \Phi_{1}) \qquad (B.5)$$

$$P_5 = -\frac{1}{4}(s_{\eta}h_{\xi} + s_{\eta,j-1}h_{\xi,j-1})(1 - \Phi_1) \qquad (B.6)$$

$$P_6 = \frac{1}{4}(s_{\xi}h_{\eta} + s_{\xi,i-1}h_{\eta,i-1})(1 - \Phi_{nz}) \qquad (B.7)$$

$$P_7 = \frac{A_{\xi}^{u}}{h_{\xi}^{u}}\left(1 - \frac{1}{4}\left[\frac{s_{\xi}^{2}}{\gamma} + \frac{s_{\xi,i-1}^{2}}{\gamma_{i-1}}\right]\Phi_{1}\mu\right) + \frac{1}{4}(s_{\xi}h_{\eta} - s_{\xi,i-1}h_{\eta,i-1})(\Phi_{nz} + \Phi_{1}) \qquad (B.8)$$

$$P_8 = -\frac{1}{4}(s_{\xi}h_{\eta} + s_{\xi,i-1}h_{\eta,i-1})(1 - \Phi_1) \qquad (B.9)$$

$$P_9 = -\frac{1}{8}\left(\frac{s_{\xi,i-1}s_{\eta,i-1}h_{\sigma,i-1}}{\gamma_{i-1}} + \frac{s_{\xi,j-1}s_{\eta,j-1}h_{\sigma,j-1}}{\gamma_{j-1}}\right)\Phi_1\mu \qquad (B.10)$$

$$P_{10} = \frac{1}{8}\left(\frac{s_{\xi,i-1}s_{\eta,i-1}h_{\sigma,i-1}}{\gamma_{i-1}} + \frac{s_{\xi,j+1}s_{\eta,j+1}h_{\sigma,j+1}}{\gamma_{j+1}}\right)\Phi_1\mu \qquad (B.11)$$

582 As it is seen, only 8 coefficients are needed away from the top and bottom boundaries. This is a manifestation of the both the compactness and the symmetry of the stencil.
## 585 Appendix C. Variational derivation

In this section we recover the momentum equations (21) with a variational derivation using Hamilton’s principle. However we were able to do it only for the case of a non-breathing grid, so that all terms involving $\partial_t z$, $\partial_t s_\xi$ and $\partial_t s_\eta$ vanish. We also remove the full Coriolis force because the variational approach requires to have horizontal Lamé coefficients ($h_\xi, h_\eta$) that vary in the vertical, which corresponds to the deep atmosphere assumption in atmospheric models (Tort and Dubos, 2014).

To derive the momentum equations we follow Holm et al. (1998) variational derivation and adapt it to the case of non-orthogonal curvilinear coordinates. We start from the Lagrangian

$$ \mathcal{L}(U_i, \rho, b, p) = \iiint \left[ \rho \left( \frac{1}{2} u_i U_i + bz \right) - p \left( \frac{\rho}{\rho_0} - 1 \right) \right] V d\xi d\eta d\sigma, \quad (C.1) $$

where the index notation $i = 1, 2, 3$ stands for $(\xi, \eta, \sigma)$, $U_i$'s are the components of the flux $\mathbf{U}$, $u_i$'s are the component of the momentum $\mathbf{u}$ and $b$ is the buoyancy relative to the Boussinesq reference density $\rho_0$. The variations of the Lagrangian are

$$ \begin{aligned} \frac{1}{\rho} \frac{\delta \mathcal{L}}{\delta U_i} &= u_i \\ \frac{1}{\rho} \frac{\delta \mathcal{L}}{\delta b} &= z \\ \frac{\delta \mathcal{L}}{\delta \rho} &= \frac{1}{2} u_i U_i + bz - \frac{p}{\rho_0} \\ \frac{\delta \mathcal{L}}{\delta p} &= \frac{\rho}{\rho_0} - 1. \end{aligned} \quad (C.2) $$

The total pressure $p = p^h + q$ arises as a Lagrange multiplier enforcing the incompressibility constraint $\rho = \rho_0$. The Euler Poincaré equation, or
602 momentum equation, reads (Holm et al., 1998)

$$ \partial_t \mathbf{u} + \mathrm{curl} \mathbf{u} \times \mathbf{U} + \nabla (\mathbf{u} \cdot \mathbf{U}) - \nabla \frac{\delta \mathcal{L}}{\delta \rho} + \frac{1}{\rho} \frac{\delta \mathcal{L}}{\delta b} \nabla b = 0, \quad (C.3) $$

603 namely

$$ \partial_t \mathbf{u} + \mathrm{curl} \mathbf{u} \times \mathbf{U} = - \nabla \left( \frac{1}{2} \mathbf{U} \cdot \mathbf{u} + \frac{p}{\rho_0} \right) + b \nabla z. \quad (C.4) $$

604 (C.4) is the momentum equation in vector-invariant form. It is identical to 605 Tort and Dubos (2014). The formalism comes with the tracer conservation 606 equation $\partial_t \phi + \mathrm{div}(\mathbf{U}\phi) = 0$, that holds for any tracer. In particular it 607 translates into the continuity equation (17) for $\phi = \rho$. We now show how 608 (C.4) can be transformed into the flux form (21) plus a pseudo force $\tilde{\mathbf{G}}$. The 609 idea is to prove that

$$ \left[ \mathrm{curl} \mathbf{u} \times \mathbf{U} + \nabla \left( \frac{1}{2} \mathbf{U} \cdot \mathbf{u} \right) \right] \mathcal{V} = \partial_j (\mathcal{A}_j U_j \mathbf{u}) - \tilde{\mathbf{G}} \mathcal{V}. \quad (C.5) $$

610 In components the curl term reads

$$ \begin{aligned} \mathcal{V} (\mathrm{curl} \mathbf{u} \times \mathbf{U})|_i &= \epsilon_{ijk} \epsilon_{jlm} \partial_l (u_m h_m) U_k \frac{\mathcal{A}_k}{h_i} \\ &= (\delta_{kl} \delta_{im} - \delta_{km} \delta_{il}) \partial_l (u_m h_m) U_k \frac{\mathcal{A}_k}{h_i} \\ &= \partial_k (U_k \mathcal{A}_k u_i) - U_k \mathcal{A}_i \partial_i u_k + U_k \frac{\mathcal{A}_k}{h_i} (u_i \partial_k h_i - u_k \partial_i h_k), \end{aligned} \quad (C.6) $$

611 where we have used $\partial_k (U_k \mathcal{A}_k) = 0$. The next step is to rearrange

$$ \begin{aligned} \mathcal{V} \nabla \left( \frac{1}{2} \mathbf{U} \cdot \mathbf{u} \right) \Big|_i &= \frac{1}{2} \mathcal{A}_i \partial_i (U_k u_k) \\ &= \mathcal{A}_i U_k \partial_i u_k + \frac{1}{2} \mathcal{A}_i (u_k \partial_i U_k - U_k \partial_i u_k) \quad (C.7) \\ &= \mathcal{A}_i U_k \partial_i u_k + \frac{1}{2} \mathcal{A}_i u_k u_j \partial_i T_{jk}. \end{aligned} $$
Summing (C.6) and (C.7) gives (C.5) with the general expression for the pseudo-force

$$ \tilde{\mathbf{G}}_i = \frac{U_k}{h_i h_k} (u_k \partial_i h_k - u_i \partial_k h_i) - \frac{1}{2h_i} u_k u_j \partial_i T_{jk}. \quad (C.8) $$

Using the expression for $T_{ij}$ allows to explicit each component

$$ \begin{aligned} \tilde{G}_u &= + (v \partial_{\xi} h_{\eta} - u \partial_{\eta} h_{\xi}) \frac{V}{h_{\xi} h_{\eta}} + \frac{w W}{h_{\xi} h_{\sigma}} \partial_{\xi} h_{\sigma} + \frac{w U}{h_{\xi}} \partial_{\xi} s_{\xi} + \frac{w V}{h_{\xi}} \partial_{\xi} s_{\eta} \\ \tilde{G}_v &= - (v \partial_{\xi} h_{\eta} - u \partial_{\eta} h_{\xi}) \frac{U}{h_{\xi} h_{\eta}} + \frac{w W}{h_{\eta} h_{\sigma}} \partial_{\eta} h_{\sigma} + \frac{w U}{h_{\eta}} \partial_{\eta} s_{\xi} + \frac{w V}{h_{\eta}} \partial_{\eta} s_{\eta} \end{aligned} \quad (C.9) \\ \tilde{G}_w = 0, \end{aligned} $$

where we have used the shallow fluid approximation $\partial_{\sigma}h_{\xi} = \partial_{\sigma}h_{\eta} = 0$ and $\partial_{\xi}h_{\sigma} = h_{\xi}\partial_{\sigma}s_{\xi}$, $\partial_{\eta}h_{\sigma} = h_{\eta}\partial_{\sigma}s_{\eta}$. This pseudo-force can be transformed further into the form of (21–22), namely

$$ \begin{aligned} \tilde{G}_u &= G_u + \frac{wU}{h_\xi} \partial_\xi s_\xi + \frac{wV}{h_\eta} \partial_\eta s_\xi + \frac{wW}{h_\sigma} \partial_\sigma s_\xi \\ \tilde{G}_v &= G_v + \frac{wU}{h_\xi} \partial_\xi s_\eta + \frac{wV}{h_\eta} \partial_\eta s_\eta + \frac{wW}{h_\sigma} \partial_\sigma s_\eta, \end{aligned} \quad (C.10) $$

with

$$ \begin{aligned} G_u &= + (V \partial_{\xi} h_{\eta} - U \partial_{\eta} h_{\xi}) \frac{V}{h_{\xi} h_{\eta}} \\ G_v &= - (V \partial_{\xi} h_{\eta} - U \partial_{\eta} h_{\xi}) \frac{U}{h_{\eta} h_{\xi}}, \end{aligned} \quad (C.11) $$

where we have used the property

$$ \frac{1}{h_{\xi}} \partial_{\xi} s_{\eta} - \frac{1}{h_{\eta}} \partial_{\eta} s_{\xi} = \frac{1}{h_{\xi} h_{\eta}} (s_{\xi} \partial_{\eta} h_{\xi} - s_{\eta} \partial_{\xi} h_{\eta}). \quad (C.12) $$

The pseudo-force differs from the hydrostatic case by the presence of terms proportional to $w$, the vertical momentum. They vanish in the hydrostatic case because the vertical momentum is identically zero.
Acknowledgments. We thank Labex Mer for supporting this study and the two reviewers, Francis Auclair and an anonymous one, for helping improving the clarity of the paper.

626 References

627 M. Nikurashin, G. K. Vallis, A. Adcroft, Routes to energy dissipation for geostrophic flows in the Southern Ocean, Nature Geoscience 6 (1) (2013) 48–51.

630 M. J. Molemaker, J. C. McWilliams, W. K. Dewar, Submesoscale instability and generation of mesoscale anticyclones near a separation of the California Undercurrent, Journal of Physical Oceanography 45 (3) (2015) 613–629.

633 J. Callies, R. Ferrari, J. M. Klymak, J. Gula, Seasonality in submesoscale turbulence, Nature communications 6.

635 J. Gula, M. J. Molemaker, J. C. McWilliams, Gulf Stream dynamics along the southeastern US seaboard, Journal of Physical Oceanography 45 (3) (2015) 690–715.

638 P. P. Sullivan, J. C. McWilliams, C.-H. Moeng, A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows, Boundary-Layer Meteorology 71 (3) (1994) 247–276.

641 J. Marshall, C. Hill, L. Perelman, A. Adcroft, Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling, J. Geophys. Res. 102 (1997) 5733–5752.
644 Y. Kanarska, V. Maderich, A non-hydrostatic numerical model for calculating free-surface stratified flows, Ocean Dynamics 53 (3) (2003) 176–185.

645 Y. Kanarska, A. Shchepetkin, J. McWilliams, Algorithm for non-hydrostatic dynamics in the regional oceanic modeling system, Ocean Modelling 18 (3) (2007) 143–174.

649 F. Auclair, C. Estournel, J. W. Floor, M. Herrmann, C. Nguyen, P. Marsaleix, A non-hydrostatic algorithm for free-surface ocean modelling, Ocean modelling 36 (1) (2011) 49–70.

652 K. Klingbeil, H. Burchard, Implementation of a direct nonhydrostatic pressure gradient discretisation into a layered ocean model, Ocean Modelling 65 (2013) 64–77.

655 W. C. Skamarock, J. B. Klemp, A time-split nonhydrostatic atmospheric model for weather research and forecasting applications, Journal of Computational Physics 227 (7) (2008) 3465–3485.

658 E. H. Müller, R. Scheichl, Massively parallel solvers for elliptic partial differential equations in numerical weather and climate prediction, Quarterly Journal of the Royal Meteorological Society 140 (685) (2014) 2608–2624.

661 S. Williams, A. Waterman, D. Patterson, Roofline: an insightful visual performance model for multicore architectures, Communications of the ACM 52 (4) (2009) 65–76.

664 A. F. Shchepetkin, J. C. McWilliams, A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate, Journal of Geophysical Research: Oceans 108 (C3).
667 M. A. Taylor, A. Fournier, A compatible and conservative spectral element method on unstructured grids, Journal of Computational Physics 229 (17) (2010) 5879–5895.
668
669

670 M. Desbrun, E. Kanso, Y. Tong, Discrete differential forms for computational modeling, in: Discrete differential geometry, Springer, 287–324, 2008.
671

672 J. Thuburn, C. J. Cotter, A framework for mimetic discretization of the rotating shallow-water equations on arbitrary polygonal grids, SIAM Journal on Scientific Computing 34 (3) (2012) B203–B225.
673
674

675 M. J. Molemaker, J. C. McWilliams, I. Yavneh, Baroclinic instability and loss of balance, Journal of Physical Oceanography 35 (9) (2005) 1505–1517.
676

677 T. Dubos, S. Dubey, M. Tort, R. Mittal, Y. Meurdesoif, F. Hourdin, DY-NAMICO, an icosahedral hydrostatic dynamical core designed for consistency and versatility, Geoscientific Model Development Discussions 8 (2) (2015) 1749–1800.
678
679
680

681 A. J. Chorin, Numerical solution of the Navier-Stokes equations, Mathematics of computation 22 (104) (1968) 745–762.
682

683 J. B. Bell, P. Colella, H. M. Glaz, A second-order projection method for the incompressible Navier-Stokes equations, Journal of Computational Physics 85 (2) (1989) 257–283.
684
685

686 L. J. Grady, J. Polimeni, Discrete calculus: Applied analysis on graphs for computational science, Springer Science & Business Media, 2010.
687
688 M. Tort, T. Dubos, Dynamically consistent shallow-atmosphere equations with a complete Coriolis force, Quarterly Journal of the Royal Meteorological Society 140 (684) (2014) 2388–2392.

691 D. D. Holm, J. E. Marsden, T. S. Ratiu, The Euler-Poincaré equations and semidirect products with applications to continuum theories, Adv. in Math. 137 (1998) 1–81.