# A non-hydrostatic algorithm for free-surface ocean modelling

Francis Auclair, Claude Estournel, Jochem Floor, Marine Herrmann, Cyril

Nguyen, Patrick Marsaleix

► To cite this version:

Francis Auclair, Claude Estournel, Jochem Floor, Marine Herrmann, Cyril Nguyen, et al. A non-hydrostatic algorithm for free-surface ocean modelling. Ocean Modelling, 2011, 36 (1-2), pp.49-70. 10.1016/j.ocemod.2010.09.006. hal-02110149

HAL Id: hal-02110149

https://hal.science/hal-02110149

Submitted on 28 Sep 2021

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
# A non-hydrostatic algorithm for free-surface ocean modelling

Francis Auclair*, Claude Estournel, Jochem W. Floor, Marine Herrmann, Cyril Nguyen, Patrick Marsaleix

Université de Toulouse, UPS, LA (Laboratoire d'Aérologie), 14 avenue Edouard Belin, F-31400 Toulouse, France
CNRS, LA (Laboratoire d'Aérologie), F-31400 Toulouse, France

An original implementation of a non-hydrostatic, free-surface algorithm based on a pressure correction method is proposed for ocean modelling. The free surface is implemented through an explicit scheme combined with a mode-spitting method but the depth-averaged velocity and the position of the free surface are updated at each non-hydrostatic iteration. The vertical momentum equation is also integrated up to the surface enabling a natural and accurate treatment of the surface layer. The consistent specification of the numerical schemes provides balanced transfers of potential and kinetic energy. This algorithm is well-suited for implementation as a non-hydrostatic kernel on originally hydrostatic free-surface ocean models such as Symphonie (http://poc.obs-mip.fr/pages/researchTopics/modelling/symphonie/symphonie.htm) for which it has originally been developed.

Energy balances associated with the propagation of short surface waves and solitary waves are presented for two dedicated well-documented configurations over closed domains. The buoyancy flux, the work rate of the pressure force together with the power of the advective terms are evaluated and discussed for the generation and the propagation of these two types of waves. The dissipation rate is in particular shown to be several orders of magnitude smaller than the work rates of the hydrostatic and non-hydrostatic pressure forces confirming the necessity for the exchanges of energy to be numerically balanced. The algorithm is subsequently applied to the complex generation of non-linear solitary internal waves by surface tides over Georges Bank, in the Gulf of Maine. The generation and the propagation of the observed non-linear and non-hydrostatic features in this region are correctly reproduced.

## 1. Introduction

The present algorithm tackles the problem of the removal of the hydrostatic assumption in ocean free-surface models. The challenges raised by the relaxation of this assumption under the Boussinesq approximation are of several types. Firstly, no prognostic equation is available to determine the non-hydrostatic pressure and the 3D Poisson system which must consequently be solved can be very demanding in terms of computing resources. Secondly, the presence of a free surface raises additional difficulties associated with the conservation of volume and heat or salt content or with the treatment of the non-hydrostatic dispersion of the fastest surface waves at affordable computing costs. Thirdly, the vertical velocity that is computed diagnostically under the hydrostatic assumption transforms into a prognostic variable with its own momentum equation and it must at the same time satisfies kinematic conditions at the surface and at the bottom. The present

algorithm gives adapted (and most often original) answers to these different problems.

The algorithm has been implemented and developed for Symphonie (Marsaleix et al., 2008, hereafter M08) but it can easily be adapted to any hydrostatic free-surface, σ-coordinate models when and where the hydrostatic assumption has to be relaxed. Indeed, even if hydrostatic modelling has encountered incontestable successes in the representation of the dynamics of the ocean, this assumption can also suffer from severe limitations. A sufficient condition for the hydrostatic assumption to be satisfied is that $U^2/L^2N^2 \ll 1$ where U, L and N stand for typical velocity, length scale and Brunt-Väisälä frequency of the flow motion (Marshall et al., 1997a). This can equivalently be restated as $\omega^2/N^2 \ll 1$ or $(h/L)^2/Ri \ll 1$ where $\omega$, $(h/L)$ and $Ri(=N^2h^2/U^2)$ are the pulsation, the aspect ratio and the Richardson number of the motion. These conditions are in particular not satisfied for instance in convective chimneys (Marshall and Schott, 1999) or in various types of small scale instabilities leading in particular to wave breaking and the wave induced mixing (Staquet and Sommeria, 2002), or the generation and propagation of non-linear (solitary) waves (New and Da Silva, 2002) cannot be accurately represented under a hydrostatic assumption. The hydrostatic assumption must also be

<footnote>* Corresponding author. Address: Laboratoire d'Aérologie, 14 Avenue Edouard Belin, 31400 Toulouse, France. Tel.: +33 561 332775; fax: +33 561 332790.

E-mail address: francis.auclair@aero.obs-mip.fr (F. Auclair).</footnote>
associated with the "traditional" approximation, i.e. the neglect of the horizontal component of the Coriolis parameter. Indeed to maintain the mechanical energy balance, the removal of this component in the vertical momentum equation must be associated with its cancellation in the horizontal momentum equation. This leads to a truncation of the angular momentum conservation principle and, as a consequence, to an approximate balance of enstrophy (Marshall et al., 1997a).

The treatment of the pressure components (whether hydrostatic, non-hydrostatic or associated with the position of the ocean free surface) is a central concern for any non-hydrostatic algorithms. Under the Boussinesq assumption, the equation of state only includes a bulk hydrostatic pressure associated to a reference density. It consequently does not depend on the non-hydrostatic component of the pressure and a splitting of the pressure between the hydrostatic and non-hydrostatic components can provide the necessary mathematical simplifications and substantially speeds up the computation. Taking advantage of this decoupling, the three components of the momentum equations can be stepped forward in time (ignoring the non-hydrostatic pressure force) before the non-hydrostatic component is adjusted to cancel the divergent part of the velocity field, leading to a "pure projection" algorithm (Chorin, 1967). An alternative is the "pressure correction" method: only an increment of non-hydrostatic pressure is updated at each time step and the most accurate available non-hydrostatic pressure is used to step forward the velocity. Such methods are sometimes formulated in the context of the "fractional step" algorithms in which the evolution operator is advantageously split (Dukowicz and Dvinsky, 1992). The projection operator can then be viewed as a "stiff" operator, the problem being that it might commute neither with the no-slip boundary conditions nor even with the implicit viscous diffusion operator, leading to "splitting errors". In this respect, Mahadevan et al. (1996) have demonstrated that this pressure correction method is preferable to pressure projection when implicit no-slip boundaries or free surface boundaries are implemented. Armfield and Street (2002) reach the same conclusion by reviewing the most popular non-hydrostatic projection algorithms in terms of accuracy and computational cost. The present algorithm is based on this "pressure-correction" method.

The presence of a free surface induces several additional difficulties. Its modelling requires in particular very short time steps and can consequently become intractable. To reduce computing costs while retaining explicit, fully non-linear evolution of the free surface, Blumberg and Mellor (1987) proposed a "mode splitting" algorithm based on a separated computation of the external (barotropic) and internal (baroclinic) modes. However, the extension of a mode-splitting algorithm to non-hydrostatic modelling is not straightforward since the update of the free surface does not commute with the enforcement of the non-divergence of the velocity field. The downward or upward acceleration at the surface tends in particular to modify (respectively to decrease or increase) the pressure gradient associated with a given surface elevation in comparison to that obtained for instance under a hydrostatic assumption. As a consequence, the resulting coupling of the surface elevation anomaly with the non-hydrostatic component of the pressure is sometimes simplified. Marshall et al. (1997a,b) make for instance a "rigid lid approximation" by considering that the propagation of surface waves is instantaneous which basically decouples free surface and non-hydrostatic pressure. To model the propagation of surface waves, Kanarska et al. (2007) conserve a mode splitting and construct their algorithm so that non-hydrostatic corrections vanish for the barotropic velocity and for the surface elevation, but they retain the non-hydrostatic pressure force in their "barotropic" mode. Using unconditionally stable semi-implicit approach, Casulli (1999) and Stelling and Zijlema (2003) propose to compute non-hydrostatic corrections for the

barotropic velocities and for the surface elevation. The non-hydrostatic adjustment of these variables is also retained in the present algorithm but it is associated with a mode-splitting approach (in the sense of Blumberg and Mellor (1987) and a Leap-Frog time stepping.

In Casulli (1999), the surface vertical velocity is obtained from the continuity equation in the upper layer which requires the enforcement of an additional equation (his Eq. (26)) to determine the non-hydrostatic pressure component in this layer. Stelling and Zijlema (2003) or Zijlema and Stelling (2005) chose a different option since they deduce the vertical velocity at the surface from its own momentum equation. This method has been applied to a mode splitting approach in the present algorithm. The representation of the free surface is rather natural when using the so-called sigma-coordinate since the upper layer can be treated as any other layer and the surface kinematic condition can be satisfied as a natural consequence of the depth integral of the continuity equation. However the apparent simplicity of treatment of the surface and bottom boundary conditions in $\sigma$-coordinates should also not hide some recurrent difficulties associated with their use in regions cumulating steep topography gradients and steep density stratification. These configurations can indeed lead to rather large truncation errors in the horizontal pressure gradient (Mellor et al., 1994) and Auclair et al. (2000) and Marsaleix et al. (2009) recently proposed adapted algorithms and numerical schemes that result in a drastic reduction of these errors. The formulation of the non-hydrostatic "Poisson" system is finally simpler (and usually symmetric positive definite) in Cartesian coordinates than in $\sigma$-coordinates. The resulting sparse system exhibits in particular only seven non-vanishing elements per lines versus 15 in $\sigma$-coordinates.

Marsaleix et al. (2009) additionally insist on the necessity for the various numerical schemes to be consistent with each other. They show in particular that a low order "rectangular" formulation of the hydrostatic pressure gradient could perform at least as well as a higher order Standard-Jacobian (in terms of truncation errors or vorticity balance) when it is associated with a "volume conserving" treatment of the density field. M08 also shows that the hydrostatic pressure gradient scheme must agree with the advection scheme for temperature and salinity to conserve mechanical energy. Energy conservation alone does not guaranty the stability of the numerical algorithm but Ferziger (1998) shows that a non-conservative model can undergo spurious variations of the kinetic energy and numerical oscillations. As importantly, when processes are to be modeled at widely spread scale range, energy conservation can become a necessity. Indeed Floor (2009) and Floor et al. (under review) show that the transfers of energy associated with diapycnal mixing induced by internal waves in the generation and propagation regions can be several orders of magnitude smaller than the buoyancy flux for these waves (i.e. the reversible exchange between the kinetic and potential energy compartments). As a consequence, an accurate energy balance needs to be enforced before mixing can be estimated. In the proposed algorithm, the numerical schemes for the continuity equation and the non-hydrostatic pressure force are consistent with each other and the "extrapolated" horizontal velocity at the surface and at the bottom (under a free-slip boundary conditions) must be equal to the closest available horizontal velocity. Although the problem of the conservation of energy that is dealt with in the present study is a general problem whatever the coordinate system, the various results are $\sigma$-coordinate specific.

To summarize, the present non-hydrostatic is based on a "pressure-correction", "free-surface", "generalized $\sigma$-coordinate" method. It shares with Casulli or Zijlema and Stelling's approaches (Casulli, 1999; Zijlema and Stelling, 2005) a similar treatment of the coupling of the surface elevation with the non-hydrostatic pressure. A non-hydrostatic correction of the free surface elevation
is indeed deduced from the adjustment of the non-hydrostatic pressure component but a mode-spitting approach enables the treatment of the free surface dynamics at affordable computational costs in the present algorithm. The set of analytical equations upon which is based the proposed algorithm is detailed in Section 2 and its numerical implementation is presented in Section 3. The algorithm is eventually used to model the generation of non-linear internal waves by surface tides in the region of Georges Bank in Section 4 and the transfers of energy more particularly associated with both surface and internal solitary waves are more specifically studied on dedicated closed configurations.

## 2. Analytical equations

### 2.1. Velocity and pressure fields in σ-coordinates

σ-coordinates can be related to Cartesian coordinates by:

$$ 
\begin{aligned}
z(x, y, \sigma; t) &= \sigma[H(x, y) + \zeta(x, y; t)] - H(x, y) \\
&= \sigma D(x, y; t) - H(x, y)
\end{aligned}
\quad (2.1) $$

where $H$ and $\zeta$ are respectively, the depth and the surface elevation anomaly and $D = H + \zeta$ is the total depth of the water column. At the surface ($z = \zeta$), $\sigma$ equals one while at the bottom ($z = -H$) $\sigma$ vanishes. Some useful relations are additionally given in Appendix A.

In $\sigma$-coordinates, the “true” vertical velocity $w$ can be written:

$$ 
\begin{aligned}
w(x, y, \sigma; t) &\equiv \frac{dz}{dt} = v_{\sigma} + \left. \frac{dz}{dt} \right|_{\sigma} = v_{\sigma} + \sigma \left. \frac{\partial \zeta}{\partial t} + v_{\zeta} \left. \frac{\partial z}{\partial x_{\zeta}} \right|_{\sigma} \right. \\
&= v_{\sigma} + \sigma \left. \frac{d\zeta}{dt} + (\sigma - 1) \frac{dH}{dt} \right.
\end{aligned}
\quad (2.2) $$

where the subscript “$\sigma$” refers to the derivation over a $\sigma$-layer $\alpha = 1$ (for $x$) or 2 (for $y$). In the remainder of this article, the indexes $z$ and $\sigma$ are more generally used when needed in order to specify that a derivative is respectively computed keeping $z$ or $\sigma$ constant. Expression (2.2) shows that the vertical velocity can be decomposed into three components:

* A vertical velocity across the $\sigma$-surface defined by $v_{\sigma} \equiv D \frac{d\sigma}{dt}$.

* A second component associated with the (space and time) variations of the free surface. This component is linearly decreasing from the surface to the bottom where it vanishes.

* A third component associated with the (space) variations of the bathymetry (with $\frac{dH}{dt} = v_{\zeta} \frac{\partial H}{\partial x_{\zeta}}$) which is also linearly decreasing this time from the bottom to the surface where it vanishes.

The density field can be decomposed into a background density and its anomaly: $\rho_0 + \rho$. The anomaly is a function of temperature ($T$) and salinity ($S$), whereas the pressure dependency of the density is neglected in the present study. The density is given by the equation of state which can simply be written $\rho = \rho(T,S)$. Several formulations of this equation are available in SYMPHONIE but only the linearized version is used for the present study.

As usual in a pressure correction method, the pressure is written: $p_0 + p + q$ where $p_0$ and $p$ are the hydrostatic components corresponding respectively to the background components of the density and to its anomaly: $\frac{\partial(p_0+p)}{\partial z} = \frac{\partial(p_0+p)}{\partial \sigma} = -g(\rho_0 + \rho)$, “q” stands for the remaining non-hydrostatic component of the pressure. Since the background density field is homogeneous, the horizontal derivative of the background pressure can be rewritten: $\frac{\partial p_0}{\partial x_{\zeta}} = \rho_0 g \frac{\partial \zeta}{\partial x_{\zeta}}$. This explicitly shows the coupling of the hydrostatic pressure force associated with the background density with the gradient of the free surface. It obviously does not imply that the two remaining components of the pressure force associated with density anomaly and non-hydrostatic pressure are decoupled from the surface elevation.

### 2.2. Dynamical and continuity equations

Under the Boussinesq approximation, the momentum equations for the horizontal velocity $v_{\zeta}$ and the vertical velocity ($w$) can be written:

$$ \left. \frac{\partial v_{\zeta}}{\partial t} \right|_{z} = -g \left. \frac{\partial \zeta}{\partial x_{\zeta}} \right|_{z} - \frac{1}{\rho_{0}} \left. \frac{\partial p}{\partial x_{\zeta}} \right|_{z} - \frac{1}{\rho_{0}} \left. \frac{\partial q}{\partial x_{\zeta}} \right|_{z} + A_{\zeta}^{z} \quad (2.3) $$

The operator “A” includes advection, viscosity and both components of the Coriolis terms (see M08 for details). After substitution of the pressure and flux relations (Appendix A), this can be rewritten in $\sigma$-coordinates:

$$ 
\begin{aligned}
\left. \frac{\partial D v_{\zeta}}{\partial t} \right|_{\sigma} = & -gD \left. \frac{\partial \zeta}{\partial x_{\zeta}} \right|_{\sigma} - \frac{1}{\rho_0} \left[ D \left. \frac{\partial p}{\partial x_{\zeta}} \right|_{\sigma} - \left. \frac{\partial z}{\partial x_{\zeta}} \right|_{\sigma} \frac{\partial p}{\partial \sigma} \right] \\
& - \frac{1}{\rho_0} \left[ D \left. \frac{\partial q}{\partial x_{\zeta}} \right|_{\sigma} - \left. \frac{\partial z}{\partial x_{\zeta}} \right|_{\sigma} \frac{\partial q}{\partial \sigma} \right] + D A_{\zeta}^{\sigma}
\end{aligned}
\quad (2.4) $$

Once the hydrostatic component of the pressure force has been simplified with the buoyancy force, the remaining equation for the vertical momentum is given by:

$$ \left. \frac{\partial w}{\partial t} \right|_{z} = - \frac{1}{\rho_0} \frac{\partial q}{\partial z} + A_3^z \quad (2.5) $$

In $\sigma$-coordinates, this leads to:

$$ \left. \frac{\partial D w}{\partial t} \right|_{\sigma} = - \frac{1}{\rho_0} \frac{\partial q}{\partial \sigma} + D A_3^{\sigma} \quad (2.6) $$

The continuity equation states that the divergence of the velocity field must vanish:

$$ \left. \frac{\partial v_{\zeta}}{\partial x_{\zeta}} \right|_{z} + \frac{\partial w}{\partial z} = 0 \quad (2.7) $$

which after substitution of the flux relation (Appendix A) can be rewritten in $\sigma$-coordinates

$$ \left. \frac{\partial D v_{\zeta}}{\partial x_{\zeta}} \right|_{\sigma} - \frac{\partial}{\partial \sigma} \left( \left. \frac{\partial z}{\partial x_{\zeta}} \right|_{\sigma} v_{\zeta} \right) + \frac{\partial w}{\partial \sigma} = 0 \quad (2.8) $$

here and throughout the entire article summation over repeated indexes over two horizontal direction is assumed. This equation is equivalent to Eq. (22) of M08 but it is expressed here for the “true” vertical velocity $w$ and not for the $\sigma$-vertical velocity $v_{\sigma}$ (see Appendix A). The evolution of the density field is prognosticated through the state equation and the heat and salt conservation equations:

$$ \left. \frac{\partial T}{\partial t} \right|_{z} = A_{T}^{z} \quad (2.9) $$

$$ \left. \frac{\partial S}{\partial t} \right|_{z} = A_{S}^{z} \quad (2.10) $$

### 2.3. Surface and bottom boundary conditions

The bottom boundary condition is classically specified for the velocity and/or its derivatives. In a free surface configuration, an additional boundary condition has to be specified at the interface between the water and the overlying atmosphere.

At the bottom, this boundary condition depends on the scales under study. If the bottom boundary layer is explicitly represented, then both the vertical and the horizontal velocities have to vanish for $\sigma = 0$, the bottom boundary condition being a no-slip condition:

$$ w(\sigma = 0) = v_{\zeta}(\sigma = 0) = 0 \quad (2.11) $$

If the bottom boundary layer is not explicitly represented then a kinematic condition has to be enforced claiming that the fluid cannot penetrate the motionless bottom:
$$w(\sigma = 0) = -v_{\alpha}(\sigma = 0) \frac{\partial H}{\partial x_{\alpha}} \qquad (2.12)$$

Both the free-slip and no-slip conditions satisfy (2.12) and in the following sections this relation should be considered as the generic bottom boundary condition without specifying whether it is a no-slip or free-slip boundary condition.

At the surface, an additional equation can also be obtained by claiming that a fluid particulate initially located at the surface remains at this location as time goes on. Along this surface, the continuity of the normal and tangential stress tensor and a surface kinematic condition must additionally be enforced. The resulting kinematic condition is:

$$w(\sigma = 1) = \frac{\partial \zeta}{\partial t} + v_{\alpha}(\sigma = 1) \frac{\partial \zeta}{\partial x_{\alpha}} \qquad (2.13)$$

This is clearly a particular case of Expression (2.2) for $\sigma = 1$. After substitution of the bottom kinematic relation (2.12), the depth integral of the continuity equation (2.8) can be written:

$$\left. \frac{\partial D \bar{v}_{\alpha}}{\partial x_{\alpha}} \right|_{\sigma} - v_{\alpha}(\sigma = 1) \frac{\partial \zeta}{\partial x_{\alpha}} + w(\sigma = 1) = 0 \qquad (2.14)$$

where the overlined velocity ($\bar{v}_{\alpha}$) is the depth-averaged velocity ($\bar{v}_{\alpha} = \frac{1}{D} \int_{-H}^{+H} v_{\alpha} dz$). Eqs. (2.13) and (2.14) can finally be combined to obtain:

$$\left. \frac{\partial D \bar{v}_{\alpha}}{\partial x_{\alpha}} \right|_{\sigma} + \frac{\partial \zeta}{\partial t} = 0 \qquad (2.15)$$

This equation relates the horizontal divergence of depth-integrated flow to the time variation of the surface elevation anomaly ($\partial \zeta / \partial t$). It is classically used as the prognostic equation for the surface elevation anomaly in free surface models (see for instance Blumberg and Mellor, 1987) and more specifically in the present algorithm.

Finally, many studies have been dedicated to the treatment of the boundary conditions for pressure under a non-hydrostatic assumption and Gresho and Sani (1987) have demonstrated that no boundary conditions for pressure should be implemented.

In the present study, no wind stress is specified at the free surface of the ocean, and the continuity of the stress tensor can be restricted to the continuity of the pressure field across the free surface which finally leads to:

$$q(\sigma = 1) = 0 \qquad (2.16)$$

## 2.4. Energy conservation

As stated in the introduction, the discrete formulation of the evolution of momentum and volume together with the discrete formulation of the continuity equation must lead to a consistent formulation of the evolution of a discrete equivalent of the kinetic and potential energy. This problem was studied in some detail by M08 under the assumption of hydrostaticity, so that kinetic energy was in this case truncated to its horizontal component. The authors concluded that the derivation of a balance of energy is first of all based on the coherence of the numerical scheme for the hydrostatic pressure gradient with the centered part of the advective scheme for density (i.e. temperature and salinity for a linear equation of state) and on a consistent coupling of the external and internal modes. The upstream components (if any) of the momentum, temperature and salinity fields are computed separately and are added to the “diffusive” terms.

A close inspection of the non-hydrostatic equation shows that the work rate of the non-hydrostatic pressure force must vanish for a closed domain. This follows from the identity:

$$v_{\alpha} \left. \frac{\partial q}{\partial x_{\alpha}} \right|_{z} + w \left. \frac{\partial q}{\partial z} \right| = \left( \left. \frac{\partial q}{\partial x_{\alpha}} \right|_{z} + \left. \frac{\partial q w}{\partial z} \right) - \left( \left. \frac{\partial v_{\alpha}}{\partial x_{\alpha}} \right|_{z} + \left. \frac{\partial w}{\partial z} \right) \right) q = 0 \qquad (2.17)$$

where the flux divergence terms in the first parenthesis of the right hand side vanish when integrated over a closed domain. The terms in the second parenthesis vanish due to the continuity equation (2.7). This property must obviously hold numerically which means that the numerical schemes for the non-hydrostatic pressure gradient and the continuity equation cannot be chosen independently.

## 3. Finite difference algorithm

Many ocean models and, in particular, most coastal ocean models (Blumberg and Mellor, 1987; Shchepetkin and McWilliams, 2005; M08) are based on an Arakawa C-grid and use $\sigma$-coordinates. The discretization stencils of the main non-hydrostatic operators are given in Fig. 1. Although this type of grid offers interesting properties in terms of energy conservation, it leads to difficulties when the kinematic relations must for instance be discretized over the upper and lower levels in $\sigma$-coordinates. Whereas the vertical velocity is available over these levels, its horizontal component is not. As a consequence, this component must be extrapolated which necessarily introduces numerical inconsistencies whatever the chosen scheme is. The amplitude of the related error should somehow go to zero as the vertical grid spacing close to the bottom and to the surface tends to zero.

To simplify the notations and without introducing any new restriction, only the terms associated with the motion along the x-axis are retained in the present section ($\alpha = 1$). The related velocity is written $u = v_1$. The numerical domain is thus restricted to $i \in [1, t_{\max}]$ for the remaining horizontal index and to $k \in [1, k_{\max}]$ for the vertical index. At density points $(i, k)$ and time $t$, the C-grid layers are defined by the total depth $D_i^t = H_i + \zeta_i^t$, the vertical position of the layers $z_{i,k}^t = (z_{i,k-1}^t + z_{i,k+1}^t)/2$ with $z_{i,k+1}^t = \sigma_{i,k+1/2} D_i^t - H_i$ (2.1) (M08). In this relation, $\sigma_{i,k+1/2}$ satisfies $0 = \sigma_{i,1/2} \le \sigma_{i,k+1/2} \le \sigma_{i,kmax+1/2} = 1$ and the thickness of the layers $\Delta z_{i,k}^t = z_{i,k+1/2}^t - z_{i,k-1/2}^t = D_i^t \Delta \sigma_{i,k} (k \in [1, k_{\max}])$ is thus directly proportional to $\Delta \sigma_{i,k}$. The position of the surface and bottom vertical-velocity layers satisfies $z_{i,kmax+1/2}^t = \zeta_i^t$ and $z_{i,k-1/2}^t = -H_i$. The resulting layer thickness at vertical velocity points $(i, k+1/2)$ is then given by $\Delta z_{i,k+1/2}^t = z_{i,k+1}^t - z_{i,k}^t$ for $k \in [1, k_{\max}]$, $\Delta z_{i,kmax+1/2}^t = z_{i,kmax+1/2}^t - z_{i,kmax-1/2}^t$ at the surface. The grid characteristics at horizontal velocity points $(i+1/2, k)$ are classically obtained by averaging the corresponding characteristics at neighboring density points.

In a pressure correction algorithm, the velocity and surface elevation are basically obtained in two steps. At time step $t + \Delta t$, “provisional” fields (written $\bar{u}_{i+1/2,k}^{t+\Delta t}$, $\bar{w}_{i,k+1/2}^{t+\Delta t}$ and $\bar{\zeta}_{i,k}^{t+\Delta t}$) are first obtained by solving the discrete analogue of the momentum equations (2.4) and (2.6) and of the volume conservation equation (2.15). Fields “increments” ($\delta u_{i+1/2,k}$, $\delta w_{i,k+1/2}$ and $\delta \zeta_{i,k}$) must then be computed so that the “final” fields ($u_{i+1/2,k}^{t+\Delta t} = \bar{u}_{i+1/2,k}^{t+\Delta t} + \delta u_{i+1/2,k}$, $w_{i,k+1/2}^{t+\Delta t} = \bar{w}_{i,k+1/2}^{t+\Delta t} + \delta w_{i,k+1/2}$ and $\zeta_{i,k}^{t+\Delta t} = \bar{\zeta}_{i,k}^{t+\Delta t} + \delta \zeta_{i,k}$) satisfy both the momentum equations, the continuity equation (2.8) and the kinematic conditions (2.12) and (2.13). In the present section, we consider that the algorithm has been integrated until time step $t$ and model fields are computed at $t + \Delta t$ and we present the different numerical schemes insisting more particularly on the coherent choices made for the discrete analogues of the continuity equation (2.8), the kinematic conditions (2.12) and (2.13) and the non-hydrostatic component of the pressure force.

### 3.1. Discrete formulation of the continuity equation

In $\sigma$-coordinates and at the density point $(i, k)$, the discrete form of Eq. (2.8) at $t + \Delta t$ can be written:
<figure><img src="image_2.png" /><figcaption>Fig. 1. Vertical section (at j) of the Arakawa C-grid stencil in σ-coordinates used for the discretization of the continuity equation (C) (at points (i,k)) and for the horizontal and vertical components of the non-hydrostatic pressure gradients (G) (respectively at points (i + 1/2,k) and (i,k + 1/2)). "Integer" levels (i and k) are represented by solid lines and "half-levels" (i + 1/2 and k + 1/2) by dashed lines. Three implementations of each operator are respectively shown at the surface, in the center of the water column and at the bottom. Red, grid and yellow lines indicate the dependency to a particular variable and α and β are the multiplying factor given in Section 3.1. Red stripped lines indicate dependencies of the continuity operator without computations of velocity increments. Pink and grey dashed rectangles indicate the variables implemented in the surface and bottom kinematic conditions. For the continuity and horizontal non-hydrostatic pressure gradient, the σ-corrective terms are represented independently on the right part of their respective panel while the vertical and horizontal derivatives (along σ-surfaces for the latter) are represented in the left part. Horizontal velocity and non-hydrostatic pressure in grey indicate variables at a non-classical position for a C-grid. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$$ \frac{\Delta z^{t+\Delta t} u^{t+\Delta t}_{i+\frac{1}{2},k} - \Delta z^{t+\Delta t} u^{t+\Delta t}_{i-\frac{1}{2},k}}{\Delta x} + w^{t+\Delta t}_{i,k+\frac{1}{2}} - w^{t+\Delta t}_{i,k-\frac{1}{2}} \\ - \left( \left[ \frac{\partial z}{\partial x} \right]_{\sigma} u \right]_{i,k+\frac{1}{2}}^{t+\Delta t} - \left[ \frac{\partial Z}{\partial x} \right]_{\sigma} u \left]_{i,k-\frac{1}{2}}^{t+\Delta t} \right) = 0 \quad (3.1) $$

where Δx is the horizontal grid resolution and Δt is the time step for the internal mode. However, unlike in the hydrostatic version of the model (M08), the “true” vertical velocity ($w^{t+\Delta t}_{i,k+\frac{1}{2}}$) rather than the σ-coordinate vertical velocity ($v_σ$) is used in this particular formulation of the continuity equation. $w^{t+\Delta t}_{i,k+\frac{1}{2}}$ is indeed the true prognostic variable for the vertical velocity in the non-hydrostatic formulation of the algorithm and it is calculated with the vertical momentum equation. At the opposite, $v_σ$ is a diagnostic variable that can be deduced from the “true” vertical velocity (see Appendix A) and it is exclusively used in the vertical advection of the momentum and tracers.

The discrete analogue of the analytical equation (2.2) for the true vertical velocity can finally be written:

$$ w_{i,k+\frac{1}{2}}^{t} = v_{\sigma,i,k+\frac{1}{2}}^{t} + \sigma_{i,k+\frac{1}{2}} \frac{\zeta_{i}^{t+\Delta t} - \zeta_{i}^{t-\Delta t}}{2\Delta t} + \left[ \frac{\partial z}{\partial x} \right]_{\sigma} u \left]_{i,k+\frac{1}{2}}^{t} \right) = 0 \quad (3.2) $$

The last two terms on the right-hand side of (3.1) are the σ-correction terms to the horizontal derivative. They have already been written here in a “flux”-like form that classically leads to the elimination of the terms associated with the “inner” layers ($k \in [1, k_{max}]$) when the continuity equation is integrated from the bottom to the surface of the ocean (Section 3.2). The discrete analogue of $[\frac{\partial z}{\partial x}]_{\sigma} u]_{i,k+\frac{1}{2}}^{t+\Delta t}$ is not unique but this choice determines the numerical schemes of the non-hydrostatic component of the pressure force (Appendix D) and of the bottom and surface kinematic conditions (Section 3.2). A low-order four-point scheme is chosen for $k \in [1, k_{max}]$:

$$ \left[ \frac{\partial z}{\partial x} \right]_{\sigma} u \left]_{i,k+\frac{1}{2}}^{t+\Delta t} \right) = \alpha_{i+\frac{1}{2},k+1}^{t+\Delta t} u_{i+\frac{1}{2},k+1}^{t+\Delta t} + \beta_{i+\frac{1}{2},k}^{t+\Delta t} u_{i+\frac{1}{2},k}^{t+\Delta t} + \alpha_{i-\frac{1}{2},k+1}^{t+\Delta t} u_{i-\frac{1}{2},k+1}^{t+\Delta t} \\ + \beta_{i-\frac{1}{2},k}^{t+\Delta t} u_{i-\frac{1}{2},k}^{t+\Delta t} \quad (3.3) $$

$$ \bullet \text{ with for } k \in [1, k_{\max}] : \alpha_{i+\frac{1}{2},k}^{t+\Delta t} = \frac{\Delta z_{i+\frac{1}{2},k}^{t+\Delta t}}{4\Delta z_{i+\frac{1}{2},k-\frac{1}{2}}^{t+\Delta t}} \frac{z_{i+1,k}^{t+\Delta t} - z_{i,k}^{t+\Delta t}}{\Delta x} \quad (3.4) $$

$$ \bullet \text{ and for } k \in [1, k_{\max}]: \beta_{i+\frac{1}{2},k}^{t+\Delta t} = \frac{\Delta z_{i+\frac{1}{2},k}^{t+\Delta t}}{4\Delta z_{i+\frac{1}{2},k+\frac{1}{2}}^{t+\Delta t}} \frac{z_{i+1,k}^{t+\Delta t} - z_{i,k}^{t+\Delta t}}{\Delta x} \quad (3.5) $$
At the lower and upper most levels (respectively), these terms take a particular form:

* at the bottom (in free-slip condition) :

$$ \left[ \frac{\partial z}{\partial x} \right]_{\sigma} \Big|_{i, \frac{1}{2}}^{t+\Delta t} = \alpha_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} u_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} + \alpha_{i-\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} u_{i-\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} \quad (3.6) $$

* at the surface :

$$ \left[ \frac{\partial z}{\partial x} \right]_{\sigma} \Big|_{i, k_{\max} + \frac{1}{2}}^{t+\Delta t} = \beta_{i+\frac{1}{2}, k_{\max} + \frac{1}{2}}^{t+\Delta t} u_{i+\frac{1}{2}, k_{\max} + \frac{1}{2}}^{t+\Delta t} + \beta_{i-\frac{1}{2}, k_{\max} + \frac{1}{2}}^{t+\Delta t} u_{i-\frac{1}{2}, k_{\max} + \frac{1}{2}}^{t+\Delta t} \quad (3.7) $$

where $\beta_{i+\frac{1}{2}, k_{\max} + \frac{1}{2}}^{t+\Delta t} = \frac{\zeta_{i+1}^{t+\Delta t} - \zeta_{i-1}^{t+\Delta t}}{2\Delta x}$ and $\alpha_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} = -\frac{H_{i+1} - H_{i-1}}{2\Delta x}$. Under a no-slip boundary condition, the right hand side of (3.6) and the velocity vanish trivially ($u_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} = w_{i, \frac{1}{2}}^{t+\Delta t} = 0$). The horizontal velocity finally needs to be “extrapolated” at the surface ($u_{i+\frac{1}{2}, k_{\max} + \frac{1}{2}}^{t+\Delta t}$) and at the bottom ($u_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t}$) under a free-slip assumption (see Section 3.5). No such relation is needed for the horizontal velocity under the hydrostatic assumption. Indeed, in this particular case, the continuity equation can be written in terms of the vertical velocity $v_\sigma$ (see Eq. (A.4) in Appendix A) leading to trivial kinematic conditions at the surface and at the bottom.

The red panel in Fig. 1 gives an overall description of the discretization stencil for Eq. (3.1). The horizontal along-$\sigma$ and vertical terms (first four terms on the left hand side of Eq. (3.1) are shown on the left side of the panel whereas the $\sigma$-correction (last two terms in (3.1)) are shown on the right part.

## 3.2. Provisional free surface anomaly and kinematic conditions

The provisional evolution of the free surface anomaly is computed by means of an explicit scheme and the chosen implementation is fully non-linear and conservative (see M08 for a detailed description of the hydrostatic implementation). The CFL criterion for surface motions such as surface gravity waves is however rather restrictive and a mode-spitting method is consequently implemented, the external-mode time ($\Delta t_e$) satisfying $\Delta t_e \le \Delta t$. A discrete analogue of Eq. (2.15) for the external mode is thus:

$$ \tilde{\zeta}_{i}^{t+\Delta t_e} = \tilde{\zeta}_{i}^{t-\Delta t_e} + \frac{2\Delta t_e}{\Delta x} \left( \left[ \tilde{\mathbf{D}}\tilde{\mathbf{u}} \right]_{i+\frac{1}{2}}^{t} - \left[ \mathbf{D}\tilde{\mathbf{u}} \right]_{i-\frac{1}{2}}^{t} \right) \quad (3.8) $$

where bolded variables are associated with the external mode. The internal-mode and external-mode surface elevation anomaly are related by:

$$ \tilde{\zeta}_{i}^{t+2\Delta t} = \frac{1}{2} (\tilde{\zeta}_{i}^{t+2\Delta t} + \tilde{\zeta}_{i}^{t+2\Delta t+\Delta t_e}) \quad (3.9) $$

The internal-mode and external-mode depth averaged velocities are related by:

$$ [\tilde{D}\tilde{u}]_{i}^{t+\Delta t} = \langle [\tilde{D}\tilde{u}]_{i}^{t\epsilon} \rangle_{t\epsilon \in [t, t+\Delta t]} + \langle [\tilde{D}\tilde{u}]_{i}^{t\epsilon} \rangle_{t\epsilon \in [t+\Delta t, t+2\Delta t]} \quad (3.10) $$

where $\langle \cdot \rangle_{t\epsilon \in [t, t+\Delta t]}$ refers to the discrete time average between $t$ and $t+\Delta t$. The discrete analogue of Eq. (2.15) for the internal-mode can be obtained by summing (3.8) between $t$ and $t+2\Delta t$ and the internal-mode provisional free surface anomaly ($\zeta_{i}^{t+2\Delta t}$) consequently satisfies:

$$ \tilde{\zeta}_{i}^{t+2\Delta t} = \zeta_{i}^{t} - \frac{2\Delta t}{\Delta x} \left( [\tilde{D}\tilde{u}]_{i+\frac{1}{2}}^{t+\Delta t} - [\tilde{D}\tilde{u}]_{i-\frac{1}{2}}^{t+\Delta t} \right) \quad (3.11) $$

Unlike in Kanarska et al. (2007), a non-hydrostatic increment is computed for the internal-mode surface elevation at each time step

(see Section 3.4). A non-hydrostatic increment is also added to the external-mode surface elevation so that relation (3.9) remains valid for true surface elevations.

The discretization schemes for the bottom and surface kinematic conditions must agree with those chosen for the continuity equation and cannot be chosen independently. By integrating vertically Eq. (3.1), it becomes indeed clear that the discretization scheme of the bottom kinematic condition must more specifically agree with the continuity equation if one wishes that the surface kinematic condition be simultaneously satisfied. This can be shown by summing the discrete continuity equation (3.1) over the depth for final the velocity field:

$$ \frac{1}{\Delta x} \left( [\tilde{D}\tilde{u}]_{i+\frac{1}{2}}^{t+\Delta t} - [\tilde{D}\tilde{u}]_{i-\frac{1}{2}}^{t+\Delta t} \right) + w_{i, k_{\max}+\frac{1}{2}}^{t+\Delta t} - \left( \beta_{i+\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t} u_{i+\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t} + \beta_{i-\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t} u_{i-\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t} \right) \\ - w_{i, \frac{1}{2}}^{t+\Delta t} + \left( \alpha_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} u_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} + \alpha_{i-\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} u_{i-\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} \right) = 0 \quad (3.12) $$

where the provisional horizontal velocities at the bottom ($\tilde{u}_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t}$) and at the surface ($\tilde{u}_{i-\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t}$) are once more not considered at their usual location on a C-grid and must thus be deduced from the provisional velocity in the neighbourhood (see Section 3.5). The chosen bilinear scheme (in the $x-z$ section) for $[\partial z/\partial x]_{i,k+\frac{1}{2}}^{t+\Delta t}$ leads to a natural formulation of the surface and bottom kinematic conditions. Indeed, the last three terms in (3.12) provide the discrete formulation of the kinematic boundary condition (2.12) that is first implemented to compute the bottom provisional vertical velocity:

$$ w_{i, \frac{1}{2}}^{t+\Delta t} = -\frac{1}{2} \left( \frac{(H_{i+1} - H_i)}{\Delta x} u_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} + \frac{(H_i - H_{i-1})}{\Delta x} u_{i-\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} \right) \quad (3.13) $$

The implementation of Eq. (3.13) is represented by a grey dashed rectangle in Fig. 1. Under the no-slip boundary condition, this equation trivially simplifies to: $\tilde{w}_{i, \frac{1}{2}}^{t+\Delta t} = \tilde{u}_{i-\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} = 0$.

The discrete formulation of the surface kinematic boundary condition can finally be recovered by substituting Eqs. (3.11) and (3.13) into (3.12), leading to the following formulation of the kinematic condition at the surface:

$$ -\frac{\zeta_{i}^{t+2\Delta t} - \zeta_{i}^{t}}{2\Delta t} + w_{i, k_{\max}+\frac{1}{2}}^{t+\Delta t} - \frac{1}{2} \left( \frac{(\zeta_{i+1}^{t+\Delta t} - \zeta_{i}^{t+\Delta t})}{\Delta x} u_{i+\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t} + \frac{(\zeta_{i}^{t+\Delta t} - \zeta_{i-1}^{t+\Delta t})}{\Delta x} u_{i-\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t} \right) = 0 \quad (3.14) $$

The implementation of (3.14) is given by a red dashed rectangle in Fig. 1. The surface and bottom boundary conditions are consequently treated on very different dynamical grounds and the vertical momentum equation is integrated all the way to the surface layer whereas the bottom vertical velocity is deduced from the bottom kinematic condition. This is a consequence of the very different dynamical behaviors of these layers due to the possible propagation of surface waves which has no equivalent at the bottom. No non-hydrostatic correction is further added to the vertical velocity along this bottom layer and its non-hydrostatic increment vanishes leading to $w_{i, \frac{1}{2}}^{t+\Delta t} = \tilde{w}_{i, \frac{1}{2}}^{t+\Delta t}$ whereas at the surface $w_{i, k_{\max}+\frac{1}{2}}^{t+\Delta t} \neq \tilde{w}_{i, k_{\max}+\frac{1}{2}}^{t+\Delta t}$. For the same reason we shall impose that $u_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t} = \tilde{u}_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t}$ so that the final velocity at time step $t+\Delta t$ ($u_{i+\frac{1}{2}, \frac{1}{2}}^{t+\Delta t}, w_{i, \frac{1}{2}}^{t+\Delta t}$) satisfies the kinematic condition at the bottom (3.13).

## 3.3. Provisional horizontal and vertical velocity

At time step $t+\Delta t$, the provisional velocity ($\tilde{u}_{i+\frac{1}{2}, k}^{t+\Delta t}, \tilde{w}_{i, k+\frac{1}{2}}^{t+\Delta t}$) can be obtained by integrating the momentum equations (2.4) and (2.6). This provisional field has no reason to also satisfy the continuity
equation (3.1) and a non-hydrostatic correction (the incremental velocity: $\delta u_{i+\frac{1}{2},k}$, $\delta w_{i+\frac{1}{2}}$) must then be computed so that the final velocity field $(u_{i+\frac{1}{2},k}^{+\Delta t}, w_{i,k+\frac{1}{2}}^{+\Delta t})$ satisfies both the momentum and continuity equations. The depth-averaged velocity is written $(\tilde{u}_{i+\frac{1}{2},k}^{+\Delta t}, \tilde{w}_{i,k+\frac{1}{2}}^{+\Delta t})$. The momentum equations for a Leap-Frog time stepping consequently reads:

$$ \bullet \text{ for } k \in [1, k_{\text{max}}]: [\Delta z \tilde{u}]_{i+\frac{1}{2},k}^{+\Delta t} = [\Delta z u]_{i+\frac{1}{2},k}^{+\Delta t} + 2\Delta t \left( -\frac{\Delta z_{i+\frac{1}{2},k}^{t} \partial q}{\rho_0 \partial x} \Bigg|_{z,i+\frac{1}{2},k}^{t} + f_{u,i+\frac{1}{2},k}^{t} \right) \quad (3.15) $$

$$ \bullet \text{ for } k \in [1, k_{\text{max}}]: [\Delta z \tilde{w}]_{i,k+\frac{1}{2}}^{+\Delta t} = [\Delta z w]_{i,k+\frac{1}{2}}^{+\Delta t} + 2\Delta t \left( -\frac{\Delta z_{i,k+\frac{1}{2}}^{t} \partial q}{\rho_0 \partial z} \Bigg|_{z,i,k+\frac{1}{2}}^{t} + f_{w,i,k+\frac{1}{2}}^{t} \right) \quad (3.16) $$

where $f_u$ and $f_w$ are non-linear operators including momentum advection and viscosity, the Coriolis pseudo-force and the hydrostatic component of the pressure force. A complete description of the hydrostatic kernel is given in M08 and the discrete analogue of the vertical momentum equation and Coriolis pseudo-force are detailed in Appendix B.

Since the algorithm is based on a "pressure correction" method, the non-hydrostatic pressure force is implemented both in the provisional horizontal and vertical momentum equations. No additional boundary condition for the vertical velocity is necessary at the surface since the vertical momentum equation (3.16) can be integrated all the way to the free surface. The provisional (and final) vertical velocity at the lower layer $(\tilde{w}_{i,\frac{1}{2}}^{+\Delta t} = w_{i,\frac{1}{2}}^{+\Delta t})$ is computed through (3.13) under the free-slip boundary condition.

To prognosticate the external-mode depth-averaged horizontal velocity $(\tilde{u})$, the horizontal momentum equation (3.15) is depth integrated from the lower ($k=1$) to the upper ($k=k_{\text{max}}$) layer and solved with a smaller time step $\Delta t_e$:

$$ [\tilde{\mathbf{D}}\tilde{\mathbf{u}}]_{i+\frac{1}{2}}^{t+\Delta t_e} = [\tilde{\mathbf{D}}\tilde{\mathbf{u}}]_{i+\frac{1}{2}}^{t-\Delta t_e} + 2\Delta t_e - \frac{1}{\rho_0} \sum_{k=1}^{k_{\text{max}}} \Delta \tilde{\mathbf{z}}_{i+\frac{1}{2},k}^{t} \frac{\partial q}{\partial x} \Bigg|_{z,i+\frac{1}{2},k}^{t-\frac{1}{2}\Delta t} + \mathbf{F}_{u,i+\frac{1}{2},k}^{t} \quad (3.17) $$

The baroclinic and non-linear terms (non-linear fluxes, baroclinic and non-hydrostatic pressure forces) are frozen.

We have shown in Section 2.4 that the work rate of the non-hydrostatic component of the pressure force must vanish over a closed domain (2.17). A consequence is that the numerical scheme for the non-hydrostatic pressure gradient cannot be chosen independently from the one used for the continuity equation (see Appendix D for details). For $k \in [1, k_{\text{max}}]$ and at time step $t$, the non-hydrostatic pressure force corresponding to the chosen numerical scheme for the continuity equations (3.1), (3.3)-(3.6) and (3.7) must consequently be written as a linear time-dependent function of the non-hydrostatic pressure:

$$ \frac{\partial q}{\partial x} \Bigg|_{z,i+\frac{1}{2},k}^{t} = \frac{q_{i+1,k}^{t} - q_{i,k}^{t}}{\Delta x} - \frac{z_{i+1,k}^{t} - z_{i,k}^{t}}{\Delta x} \frac{\partial q}{\partial z} \Bigg|_{i+\frac{1}{2},k}^{t} \quad (3.18) $$

$$ \frac{\partial q}{\partial z} \Bigg|_{i,k+1/2}^{t} = \frac{q_{i,k+1}^{t} - q_{i,k}^{t}}{\Delta z_{i,k+1/2}^{t}} \quad (3.19) $$

where the last term of (3.19) is a $\sigma$-coordinate "correction". The vertical gradient at point $(i+\frac{1}{2}, k)$ must then satisfy:

$$ \frac{\partial q}{\partial z} \Bigg|_{i+1/2, k}^{t} = \frac{1}{4} \frac{q_{i+1, k+1}^{t} + q_{i,k+1}^{t} - q_{i+1, k}^{t} - q_{i,k}^{t} + q_{i+1, k}^{t} + q_{i,k}^{t} - q_{i+1, k-1}^{t} - q_{i,k-1}^{t}}{\Delta z_{i+1/2, k+1/2}^{t} \Delta z_{i+1/2, k-1/2}^{t}} \quad (3.20) $$

for $k \in [1, k_{\text{max}}]$. Taking further advantage of the fact that the non-hydrostatic pressure must vanish at the surface (2.16), a backward formulation of the pressure gradient matrix is imposed over the

upper layer: $(\frac{\partial q}{\partial z})_{i,k_{\text{max}}+\frac{1}{2}}^{t} = \frac{0-q_{i,k_{\text{max}}}^{t}}{\Delta z_{i,k_{\text{max}}+\frac{1}{2}}^{t}}$ and $(\frac{\partial q}{\partial z})_{i+\frac{1}{2},k_{\text{max}}+\frac{1}{2}}^{t} = \frac{0-q_{i+1,k_{\text{max}}-\frac{1}{2}}^{t}-q_{i,k_{\text{max}}}^{t}}{\Delta z_{i+\frac{1}{2},k_{\text{max}}+\frac{1}{2}}^{t}}$.

The green and yellow panels of Fig. 1 show the implementation of the non-hydrostatic pressure gradient in a vertical section. The along-$\sigma$ and $\sigma$-corrective components of the non-hydrostatic pressure gradient (first and second terms on the right hand side of Eq. (3.18)) are respectively, shown in the left and right part of green panel.

### 3.4. Poisson equation and non-hydrostatic increments

As the provisional velocity does not satisfy the continuity equation, it must be corrected or "projected" onto the divergence-free solution submanifold. An increment of non-hydrostatic pressure ($\delta q$) is thus calculated and it finally provides the needed velocity increment following:

$$ \bullet \text{ for } k \in [1, k_{\text{max}}]: \delta u_{i+\frac{1}{2},k} = -\frac{2\Delta t}{\rho_0} \frac{\partial \delta q}{\partial x} \Bigg|_{z,i+\frac{1}{2},k}^{t+\Delta t} \quad (3.21) $$

$$ \bullet \text{ for } k \in [1, k_{\text{max}}]: \delta w_{i,k+\frac{1}{2}} = -\frac{2\Delta t}{\rho_0} \frac{\partial \delta q}{\partial z} \Bigg|_{z,i+\frac{1}{2},k}^{t+\Delta t} \quad \text{and} \quad \delta w_{i,\frac{1}{2}} = 0 \quad (3.22) $$

The time positioning of the non-hydrostatic increments is rather arbitrary and is not specified. To cancel the divergence of the velocity, the non-hydrostatic increment ($\delta q$) must then satisfy:

$$ \nabla_z^2(\delta q) = \rho_0 \frac{(\vec{\nabla} \cdot \vec{\nabla}^{t+\Delta t})}{2\Delta t} \quad (3.23) $$

where $\nabla_z^2$ is the Laplacian of the non-hydrostatic pressure field in z-coordinates and $\vec{v} = (u, w)$. This equation is obtained by substituting (3.21) and (3.22) into (3.1). The discretization over the C-grid of the resulting Poisson equation (3.23) must consequently be deduced from the discrete continuity equation (3.1) and from the non-hydrostatic pressure gradient force (based on (3.18)-(3.20)). This is further detailed in Appendix C.3.

Once the non-hydrostatic pressure increment has been calculated, the final non-hydrostatic pressure can be computed $(q_{i,k}^{t+\Delta t} = \tilde{q}_{i,k}^{t+\Delta t} + \delta q_{i,k})$. The incremental velocity can be deduced from Eqs. (3.21) and (3.22) to obtain the final internal-mode velocity (Fig. 2): $u_{i+\frac{1}{2},k}^{t+\Delta t} = \tilde{u}_{i+\frac{1}{2},k}^{t+\Delta t} + \delta u_{i+\frac{1}{2},k}$, $w_{i,k+\frac{1}{2}}^{t+\Delta t} = \tilde{w}_{i,k+\frac{1}{2}}^{t+\Delta t} + \delta w_{i,k+\frac{1}{2}}$. The non-hydrostatic pressure component and the hydrostatic surface pressure component are coupled, which implies that the horizontal divergence of the depth-integrated velocity increment $[D^{t+\Delta t} \delta \bar{u}]_{i+\frac{1}{2}}$ does not vanish and an incremental surface elevation anomaly must be computed:

$$ \delta \zeta_i = -\frac{2\Delta t}{\Delta x} \left( [D^{t+\Delta t} \delta \bar{u}]_{i+\frac{1}{2}} - [D^{t+\Delta t} \delta \bar{u}]_{i-\frac{1}{2}} \right) \quad (3.24) $$

We have shown in Section 3.1 that the surface elevation and velocity fields automatically satisfy the surface kinematic condition (3.14) and, unlike at the bottom boundary, their increments do not need to vanish. The final depth of the water column ($D_i^{t+2\Delta t}$) together with the grid position ($z_{i,k}^{t+2\Delta t}$) and thickness ($\Delta z_{i,k}^{t+2\Delta t}$) can in turn be updated. As both the surface elevation anomaly and the depth-averaged velocity have non vanishing non-hydrostatic increments, a correction must also be added to the corresponding
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Schematic representation of the mode-spitting algorithm used in Symphonie-NH. (Yellow): internal mode, (blue): external mode, (red): non-hydrostatic pressure correction. Blue and yellow horizontal arrows stand each for one integration of respectively the external and internal mode kernels, the two lines of arrows for each mode represent the overlapping Leap-Frog integrations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<caption>**Table 1**</caption>

CPU time in seconds for several preconditioning methods and solvers for 1000 iterations.

<table><thead><tr><th></th><th>TEST 2D</th><th>TEST 2D-large</th><th>TEST 3D</th></tr></thead><tbody><tr><td>MUMPS</td><td>9.11</td><td>5.82 × 10<sup>4</sup></td><td>4.12 × 10<sup>3</sup></td></tr><tr><td>ILUT and FGMRES</td><td>3.25</td><td>1.24 × 10<sup>3</sup></td><td>9.35 × 10<sup>2</sup></td></tr><tr><td>ILUK and FGMRES</td><td>4.10</td><td>1.65 × 10<sup>3</sup></td><td>No convergence</td></tr><tr><td>ILUT and DQGMRES</td><td>3.28</td><td>1.22 × 10<sup>3</sup></td><td>1.04 × 10<sup>2</sup></td></tr></tbody></table>

variables in the external mode. This correction is applied uniformly between $t$ and $t + 2\Delta t$ so that the relation coupling the internal and external modes (3.9) and (3.10) remains true for the final fields.

Several algorithms from the HIPS,<fnref n="1" /> HYPRE<fnref n="2" /> and SPARSKIT<fnref n="3" /> software libraries have been tested in order to solve the linear system (3.23). To this purpose, three benchmarks associated with the modelling of the free oscillations of the 10 m by 10 m box presented in the Section 4.1 have been implemented. For the first TEST (TEST 2D), the computational grid is two-dimensional (64 × 1 × 64), for the second TEST (TEST 2D-Large) a larger grid is implemented (1000 × 1 × 1000) while the remaining (TEST 3D) is three-dimensional (64 × 64 × 64). The results are given in Table 1 in terms of CPU time. To start with, a direct solver (MUMPS<fnref n="4" />) has been implemented and is used as a basis of comparison. Several types of preconditioning have first been evaluated with the FGMRES solver, then the best conditioner has been used with several solvers (Saad, 2003). For each implementation, 1000 model iterations have been completed. These tests have shown that the best configuration is the combination of a ILUT preconditioning with the FGMRES solver which appeared to be more robust than the DQGMRES. The solver converges on average in less than two iterations in the 2D case indicating that the numerical system is rather well-conditioned. In these simple configurations, computational costs have also been drastically reduced by recomputing the preconditioning only when the number of iterations of the solver reaches a chosen threshold and by optimally reorganizing the columns and the lines of the matrix system. In the different test cases shown here, only one preconditioning at the first time step was necessary. The 3D case does not show radically different behavior compared to their corresponding

2D counterparts and the ILUT preconditioner remains the best choice while the FGMRES is still preferred.

### 3.5. Boundary conditions

The boundary conditions for the horizontal component of the provisional velocity ($\bar{u}^{t+\Delta t}_{i+\frac{1}{2},k}$ at $i = 1/2$ and $i = i_{max} + 1/2$) along the open boundaries are described in Marsaleix et al. (2006). The remaining boundary conditions are best described in terms of conservation of volume, i.e. considering the implementation of the continuity equation (3.1). Along the boundary, this equation is indeed satisfied directly by adjusting the vertical velocity ($\bar{w}^{t+\Delta t}_{i,k+\frac{1}{2}}$) at the outermost points ($i = 1$ and $i_{max}$). At points inside the domain, i.e. for $i \in [2, i_{max} - 1]$, the continuity equation is satisfied by computing non-hydrostatic increments thought (3.23). To do so, several options are available to specify the non-hydrostatic pressure increment along the boundary. If the outer field is supposed to be hydrostatic, the non-hydrostatic pressure must be equal to zero along the boundary, if not, a vanishing gradient across the boundary can for instance be specified. The latter scheme is used in the Georges Bank Experiment along the open boundaries and the non-hydrostatic increments of both the horizontal and vertical velocities are canceled at outermost grid points: $\delta u_{\frac{1}{2},k} = \delta u_{i_{max}+\frac{1}{2},k} = 0$ and $\delta w_{1,k+\frac{1}{2}} = \delta w_{i_{max},k+\frac{1}{2}} = 0$.

In the particular case of closed lateral boundary conditions (such as in Experiments 1 and 3) the provisional and incremental horizontal velocities vanish at $i = 1/2$ and $i = i_{max} + 1/2$ while the provisional and incremental vertical velocities vanish at $i = 1$ and $i = i_{max}$.

The non-hydrostatic pressure ($q^{t+\Delta t}_{i,k_{max}+\frac{1}{2}}$) at the surface must additionally vanish following (2.16). This obviously does not imply that the pressure is hydrostatic in the surface layer since the non-hydrostatic pressure also satisfies $q^{t+\Delta t}_{i,k_{max}} \neq 0$.

Relations (3.13) and (3.14) additionally include extrapolated horizontal velocities. At the first glance, their implementation seems arbitrary and is only a matter of choice of an accurate relation. Although the numerical implementation of the continuity equation does not lead to restrictive conditions, the conservation of energy does. As Eq. (2.17) must indeed be satisfied numerically, the following relations hold for the horizontal velocity over the surface and bottom layers:

<footnote n="1">¹ http://hips.gforge.inria.fr/.</footnote>

<footnote n="2">² https://computation.llnl.gov/casc/sc2001_fliers/hypre/hypre01.html.</footnote>

<footnote n="3">³ http://www-users.cs.umn.edu/~saad/software/SPARSKIT/sparskit.html.</footnote>

<footnote n="4">⁴ http://mumps.enseeiht.fr/index.php.</footnote>
$$ \tilde{u}_{i+\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t}=\tilde{u}_{i+\frac{1}{2}, k_{\max}}^{t+\Delta t} \quad \text{and} \quad u_{i+\frac{1}{2}, k_{\max}+\frac{1}{2}}^{t+\Delta t}=u_{i+\frac{1}{2}, k_{\max}}^{t+\Delta t} \quad (3.25) $$

$$ u_{i+\frac{1}{2}, 2}^{t+\Delta t}=\tilde{u}_{i+\frac{1}{2}, 2}^{t+\Delta t}=\tilde{u}_{i+\frac{1}{2}, 1}^{t+\Delta t} \quad (\text{free-slip}) \quad \text{or} $$

$$ u_{i+\frac{1}{2}, 2}^{t+\Delta t}=\tilde{u}_{i+\frac{1}{2}, 2}^{t+\Delta t}=0 \quad (\text{no-slip}) \quad (3.26) $$

The second relation at the bottom holds for the provisional velocity. Indeed, in any case, the incremental velocity at the bottom vanishes ($\delta u_{i+\frac{1}{2}, 2} = 0$). A complete derivation of the work rate of the non-hydrostatic pressure force is given in Appendix D. In $\sigma$-coordinates, the energy conservation requirement thus leads to a trivial extrapolation scheme for the horizontal velocity in the surface and bottom layers (in free-slip bottom boundary). It is for instance possible to show (the details of this fastidious calculus are not reproduced here) that higher order extrapolation schemes based on the two closest non-hydrostatic pressure force and horizontal velocity points do not conserve energy in the surface layer and leads to non-physical extrapolation coefficients in the bottom layer.

### 3.6. Three-step algorithm

A step-by-step summary of the proposed non-hydrostatic free-surface mode-spitting $\sigma$-coordinate algorithm is proposed in the present section. The final velocity ($u_{i+\frac{1}{2}, k}^t, w_{i,k+\frac{1}{2}}^t$), temperature and salinity ($T_{i,k}^t$ and $S_{i,k}^t$), hydrostatic and non-hydrostatic components of the pressure field ($p_{i,k}^t$ and $q_{i,k}^t$) and density anomaly ($\rho_{i,k}^t$) are supposed to be known until time $t$. For this to be true, the external mode must have been integrated until $t + \Delta t$ (Fig. 2), the external-mode depth-averaged velocity $[Du]_i^t$ and surface elevation anomaly ($\zeta_i^t$) are consequently supposed to be available until time step $t + \Delta t$. Indeed, the final internal-mode depth-integrated mean velocity $[Du]_i^t$ is defined as the time average of its external mode counterpart between $t - \Delta t$ and $t + \Delta t$ (3.10) and the final position of the internal-mode surface elevation ($\zeta_i^{t+\Delta t}$) can thus be deduced from the external-mode surface elevation at $t + \Delta t$ (3.9). The time-dependent computational grid, i.e. the “final” position of the grid layers ($z_{i,k}^{t+\Delta t}$) and their thickness ($\Delta z_{i,k}^{t+\Delta t}$) are consequently available until time step $t + \Delta t$.

**Step 1:** Internal-mode provisional velocity and mass field (at $t + \Delta t$).

The internal-mode momentum equations (3.15) and (3.16) are integrated from time step $t - \Delta t$ to $t + \Delta t$ to compute the provisional (baroclinic) velocity at $t + \Delta t$. There is no need for additional boundary conditions at the upper level and the vertical velocity at the bottom is computed through the bottom kinematic relation (3.13) in which the horizontal velocity ($u_{i+\frac{1}{2}, 2}^{t+\Delta t}$) satisfies (3.26). The heat and salt conservation equations are integrated from $t - \Delta t$ and $t + \Delta t$ to update the temperature and salinity at $t + \Delta t$ before the density diagnostic field at the same time step is obtained from the equation of state.

**Input:** These equations require the knowledge of the final velocity, temperature and salinity at time step $t - \Delta t$ and $t$ together with the final pressure field (hydrostatic, non-hydrostatic and surface components) at time step $t$ and the final computational grid at $t + \Delta t$.

**Output:** The provisional velocity ($\tilde{u}_{i+\frac{1}{2}, k}^{t+\Delta t}, \tilde{w}_{i,k+\frac{1}{2}}^{t+\Delta t}$) and final temperature, salinity and density fields ($T_{i,k}^{t+\Delta t}, S_{i,k}^{t+\Delta t}, \rho_{i,k}^{t+\Delta t}$) are available at the end of Step 1.

**Step 2:** External-mode provisional depth-averaged velocity and surface elevation (at $t + 2\Delta t$).

The external-mode Eqs. (3.8) and (3.17) are integrated from $t + \Delta t$ to $t + 2\Delta t$ to calculate the provisional surface elevation anomaly and the provisional depth averaged horizontal current. These equations are integrated with a smaller time step ($\Delta t_e$) to satisfy the CFL criterion associated in particular with the propagation of surface wave.

**Input:** The integration of the external-mode equation is based on the provisional internal-mode (velocity, depth-integrated advective fluxes, hydrostatic and non-hydrostatic components of the pressure force...) at $t + \Delta t$. The external-mode depth-averaged provisional velocity and surface elevation anomaly are also needed at $t + \Delta t$.

**Output:** The internal-mode depth-averaged provisional velocity at $t + \Delta t$ and the provisional surface elevation at $t + 2\Delta t$ can be deduced from Eqs. (3.10) and (3.9). At the end of Step 2, the barotropic and baroclinic components of the provisional velocity are consequently known at $t + \Delta t$. They satisfy the momentum equations and the bottom kinematic condition but the provisional velocity satisfies neither the continuity equation nor the surface kinematic solution.

**Step 3:** Incremental non-hydrostatic pressure and velocity at $t + \Delta t$ and surface elevation at $t + 2\Delta t$.

Non-hydrostatic increments are then computed for the non-hydrostatic pressure component, the velocity and the surface elevation anomaly so that the resulting final velocity satisfies both the momentum and the continuity equations. The non-hydrostatic pressure increment ($\delta q$) is first obtained by solving the Poisson (Eq. (3.23) see Appendix C for a complete discrete formulation). The Poisson operator, the incremental velocities and the non-hydrostatic pressure force are computed by matrix–matrix and matrix–vector multiplications.

**Input:** The provisional velocity at $t + \Delta t$ (calculated at Steps 1 and 2) is needed to compute the right-hand side of the Poisson system (3.23) while the surface elevation, the layer position and thickness at $t + \Delta t$ are used to compute the matrix operators for the continuity equation (Appendix C.1) and non-hydrostatic pressure force (Appendix C.2).

**Output:** The incremental velocity at $t + \Delta t$ can finally be obtained from relations (3.21) and (3.22). The $\sigma$-vertical velocity ($v_\sigma|_{i,k+\frac{1}{2}}^{t+\Delta t}$) can be deduced from the true vertical velocity ($w_{i,k+\frac{1}{2}}^{t+\Delta t}$) through Eq. (3.2). The vertical velocity $v_\sigma|_{i,k+\frac{1}{2}}^{t+\Delta t}$ is then used in the vertical advection of the momentum and tracers. The non-hydrostatic pressure force at time step $t + \Delta t$ ($-\frac{1}{\rho_0} \frac{\partial q}{\partial x}|_{z=i+\frac{1}{2}, k}$) can also be calculated by updating the corresponding matrix operator ($G^{t+\Delta t}$, see Appendix C). The final surface elevation anomaly and the final depth-averaged velocity for the internal mode at $t + 2\Delta t$ satisfy respectively $\zeta_i^{t+2\Delta t} = \frac{1}{2}(\tilde{\zeta}_i^{t+2\Delta t} + \tilde{\zeta}_i^{t+2\Delta t+\Delta t_e}) + \delta\zeta_i$ and $[D\tilde{u}]_i^{t+\Delta t} = \langle [\tilde{D}\tilde{u}]_i^{t_e} \rangle_{t_e \in [t, t+2\Delta t]} + \delta[D^{t+\Delta t}\tilde{u}]_i$. The final velocity at time step $t + \Delta t$ satisfies both the continuity equation (3.1) and the following formulation of the momentum equations:

$$ [\Delta zu]_{i+\frac{1}{2}, k}^{t+\Delta t} = [\Delta zu]_{i+\frac{1}{2}, k}^{t-\Delta t} + 2\Delta t \left[ -\frac{\Delta z_{i+\frac{1}{2}, k}^t \partial q}{\rho_0 \partial x} \Bigg|_{z=i+\frac{1}{2}, k}^t - \frac{\Delta z_{i+\frac{1}{2}, k}^{t+\Delta t} \partial \delta q}{\rho_0 \partial x} \Bigg|_{z=i+\frac{1}{2}, k}^{t+\Delta t} + f_{u=i+\frac{1}{2}, k}^t \right] \quad (3.27) $$
$$[\Delta z w]_{i,k+\frac{1}{2}}^{t+\Delta t} = [\Delta z w]_{i,k+\frac{1}{2}}^{t-\Delta t} + 2\Delta t \left[ -\frac{\Delta z_{i,k+\frac{1}{2}}^t}{\rho_0} \frac{\partial q}{\partial z} \Big|_{i,k+\frac{1}{2}}^t - \frac{\Delta z_{i,k+\frac{1}{2}}^{t+\Delta t}}{\rho_0} \frac{\partial \delta q}{\partial z} \Big|_{i,k+\frac{1}{2}}^{t+\Delta t} + f_{w,i,k+\frac{1}{2}}^t \right] \quad (3.28)$$

As such, this algorithm presents several well-known drawbacks. First of all, the implicit operator for viscosity computed at Step 1 does not commute with the projection operator applied at Step 3. Indeed this operator cannot recover its pure implicit formulation (for vertical viscosity) and the provisional velocity at time step $t + \Delta t$ remains in the non-linear operator "f" (3.27) and (3.28) in place of the total velocity. A possible solution would be to implement a Lagged viscous diffusion scheme (Schlesinger et al., 1983). It would indeed conserve the filtering properties of the implicit scheme with respect to the Leap-Frog scheme and commutes with the projection operator.

## 4. Modelling surface and non-linear internal gravity waves

The propagation of non-hydrostatic surface waves is first studied in the case of the natural barotropic linear oscillations of a 10 m by 10 m box (Chen, 2003). Then the generation of non-linear internal tides over George Bank (Gulf of Maine) is investigated. This complex region cumulates large amplitude M2 tidal currents, a well-marked pycnocline and steep bathymetry slopes. The idealized tidal forcing, stratification and bathymetry together with the Coriolis dispersion are similar to those used by Lamb (1994, hereafter L94). The dynamics of the internal waves in this region is investigated in detail, deferring to future studies the modelling of permanent along bank currents associated in particular with tidal rectification (Loder, 1980). The main objective is indeed to validate the present algorithm in a configuration cumulating large topography gradients and large amplitude surface and internal waves based on in situ observations and on L94's modelling. The generation and the propagation of non-linear internal waves are finally studied in some details in the case of the oscillations of a two-layer fluid in a box. In each case, a complete balance of kinetic energy is achieved globally for the whole domain and the non-vanishing rates of transfers are analyzed giving both some numerical and physical insights.

The Georges Bank experiment can thus be considered as LES (Large Eddy Simulation) while the others are DNS (Direct Numerical Simulation) and the balances and inner transfers of energy for the different waves encountered in this realistic implementation are more specifically detailed in the first and third experiments specifically implemented over closed domains, i.e. with no complex exchange through out the open boundaries.

### 4.1. Natural oscillations in a 10 m by 10 m box (Experiment 1)

The first experiment is dedicated to a detailed inspection of the representation of short (non-hydrostatic) surface waves: the evolution of the RMS errors, the characteristic relation and the energy transfers are more particularly studied. To achieve this, the natural free surface oscillations of a homogeneous fluid in a 10 m by 10 m deep box are classically simulated.

The conditions of this experiment are similar to the ones proposed in the literature and are given in Table 2. Following Chen (2003), the free surface anomaly at $t = 0$ satisfies:

$$\eta(x, t) = 10 + 0.15 \cos \left( \frac{\pi}{10} x \right) \quad (4.1)$$

No initial velocity is specified and the density is homogeneous. The horizontal and vertical grid resolutions are equal to 20 cm. The external and internal mode time steps are chosen equal (a test of

<table><caption>**Table 2**
Model parameters for Experiment 1.</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Domain size (m)</td><td>(10, 10)</td></tr><tr><td>Domain size in points: (i<sub>max</sub>, k<sub>max</sub>)</td><td>(50, 50)</td></tr><tr><td>Δx = Δy</td><td>20 cm</td></tr><tr><td>External mode time step: Δt<sub>e</sub></td><td>6.43 × 10<sup>-3</sup> s</td></tr><tr><td>Internal mode time step: Δt</td><td>6.43 × 10<sup>-3</sup> s</td></tr><tr><td>Viscosity</td><td>10<sup>-6</sup> m<sup>2</sup>s<sup>-1</sup></td></tr><tr><td>Diffusivity</td><td>10<sup>-7</sup> m<sup>2</sup>s<sup>-1</sup></td></tr><tr><td>Advection scheme</td><td>Centered</td></tr><tr><td>Asselin filter parameter</td><td>0.1</td></tr><tr><td>Bottom boundary condition</td><td>No slip</td></tr><tr><td>Coriolis parameter</td><td>0</td></tr><tr><td>Equation of state</td><td>Linear</td></tr></tbody></table>

sensitivity is issued concerning the consequences of a larger internal mode time step). The viscosity and diffusivity are molecular, no turbulent scheme is used and the advection schemes are all centered. A no-slip boundary condition is implemented.

The period of the natural oscillations of the box is equal to $T_{NH} = 3.58$ s. However, under the hydrostatic assumption, which is in this case equivalent to the long wave approximation, this period reduces to $T_H = 2.02$ s. This difference is due to the fact that long (hydrostatic) surface waves have a higher phase velocity than shorter (non-hydrostatic) waves.

Fig. 3 shows the evolution with time of the kinetic energy together with the main energy transfers during the first 20 s. Potential energy is evaluated as in (M08) and, in the non-hydrostatic context, a vertical component is added to the total kinetic energy per unit mass which becomes:

$$K_t^{t+\frac{1}{2}} = \sum_{i=2}^{i_{max}-1} \sum_{k=1}^{k_{max}} \left( \frac{1}{2} u_{i+\frac{1}{2},k}^t \bar{u}_{i+\frac{1}{2},k}^{t+\Delta t} + \frac{1}{2} w_{i,k+\frac{1}{2}}^t \bar{w}_{i,k+\frac{1}{2}}^{t+\Delta t} \right) \quad (4.2)$$

The evolution equation of the vertical component of the kinetic energy is obtained by multiplying the vertical momentum equation (3.28) written for $k \in [1/2, k_{max} + 1/2]$ (all the way to the surface layer but not at the bottom) by the vertical component of the current ($w_{i,k+1/2}^t$). The global power evolutions are obtained by integrating the local rates over the whole computational domain and are given per unit length in the y-direction. They are shown descending in order from the largest (Fig. 3a) to the smallest (Fig. 3e).

The evolution of the horizontal and vertical kinetic energy is similar (Fig. 3a) and classically exhibits a period equal to half the period of the oscillations due to its non-linear nature. The evolution during the first 20 s does not exhibit any decreasing nor increasing tendency of the amplitude of these oscillations. The work rate of the surface elevation pressure force matches the total variation rate of the kinetic energy within a few percents (about 2%). This is a consequence of the perfect canceling (i.e. to computer precision) of the work rate of the horizontal and vertical components of the non-hydrostatic pressure force which are shown on Fig. 3b and to the small dissipation rates shown on (Fig. 3d and e). The canceling of the non-hydrostatic pressure gradient for a closed domain was discussed in Section 2.4 (Eq. (2.17)) and is thus recovered numerically. Fig. 3c shows the evolution of the non-hydrostatic increments for respectively the horizontal and vertical components of the velocity. The horizontal increment is lagging the variation of the kinetic energy (Fig. 3a) by a quarter of a period. Interestingly enough the horizontal and vertical increments have opposite signs but rigorously the same amplitude as a consequence of the symmetric nature of the flow.

The evolution of the power of the advection terms in the horizontal and vertical momentum equations together with the work rates of the horizontal and vertical components of the Asselin filter
<figure><img src="image_4.png" /><figcaption>Fig. 3. main energy transfers in Experiment 1 per unit length in the y-direction. (a): kinetic energy. Dashed line: variation rate of the kinetic energy (Kt), Solid line: variation rate of the horizontal kinetic energy (DdKl/dt, x), (+): variation rate of the vertical kinetic energy (DdKl/dt, z), (×): work rate of surface pressure force (Surf, x). (b) Solid line: work rate of the non-hydrostatic pressure force in the horizontal direction (Nh pgrad, x), (+): minus the work rate of the non-hydrostatic pressure force in the vertical direction (−Nh pgrad, z). (c) non-hydrostatic correction: Solid line: work rate of the non-hydrostatic correction of the horizontal momentum (Corr NH, x), (+): Minus the work rate of the non-hydrostatic correction of the vertical momentum (−Corr. NH, z). (d) Dashed line: power of the horizontal advection (Adv, x), Dash-dotted line: power of the vertical advection (Adv, z), Solid line: Power of the Asselin filter for the horizontal velocity, (+): power of the Asselin filter for the vertical velocity. (e) viscous diffusion: Solid line and (+) signs: power of the respectively the horizontal and vertical viscous diffusions of horizontal momentum (hor dif, x and vert dif, x), Dashed and dash-dotted lines: power of the respectively the horizontal and vertical viscous diffusions of vertical momentum (hor dif, z and vert dif, z).</figcaption></figure>

are presented on Fig. 3d. The advective terms are about one order of magnitude larger in the vertical momentum equation than in the horizontal momentum equation and both advection terms exhibit an irregular evolution due to their non-linear nature and to the rather low spatial resolution (20 cm). A period of transition of about 1 s can be observed for the work rates of the Asselin filter. Indeed during the very first second, the amplitude of these rates are one order of magnitude larger than during the remainder of the simulation (these values have been truncated in Fig. 3d). After this period, the energy transfer associated with the Asselin filter decreases, becomes periodic and its amplitude is more than three orders of magnitude smaller than the work rate of the surface pressure force. Interestingly, during the transition periods this transfer of energy can be positive while after about four periods it remains negative which could have been expected for a damping filter. These positive values are associated with an increase of the kinetic energy by the Asselin filter. As shown by M08, this is one of the interesting properties of the Leap Frog time stepping scheme in the case of high frequency variations of the velocity as energy is paradoxically injected in order to stabilize the numerical simulation.

The power of the horizontal and vertical viscous terms in the horizontal and vertical momentum equations is finally given in Fig. 3e. It is about one million times smaller than the work rate of the surface pressure force and is out of phase by about a quarter of a period. The horizontal and vertical viscous diffusions of the vertical momentum are one order of magnitude higher than the

diffusion of the horizontal component. Whereas the horizontal and vertical components of the diffusion rate in the horizontal momentum equation are equal and negative, the vertical component of the diffusion rate in the vertical momentum equation is positive meaning that it is globally responsible for an increase of the vertical kinetic energy. This can be explained by considering the following decomposition of the vertical viscous diffusion rate for a constant and homogeneous vertical diffusivity $K_z$ under the no-slip bottom boundary condition:

$$ \int_{-H}^{\eta} K_z w \frac{\partial^2 w}{\partial z^2} dz = K_z w(\eta) \left. \frac{\partial w}{\partial z} \right|_{z=\eta} - K_z \int_{-H}^{\eta} \left( \frac{\partial w}{\partial z} \right)^2 dz \quad (4.3) $$

The second term of the right hand side is the local dissipation rate and is negative. However, the first term does not vanish and can be positive at the free surface. The right hand side can consequently be either positive or negative depending in particular on the amplitude of the surface wave.

This experiment has been repeated for several ratios of internal to external mode time steps. Each experiment lasts 20 s and the time step for the external mode is based on CFL criteria (it remains equal to 6.42 ms) while the time step for the internal mode is varied to obtain successively time-step ratios ($\Delta t/\Delta t_e$) of 2 and 10. To evaluate the evolution of the amplitude of the oscillations, the value of the last maximum along the left wall is computed while the period is calculated by considering the oscillations of the same
point. The amplitude of the initial tilt is reduced to 1 mm to avoid the growing of any non-linear instability.

When the external and internal modes are integrated simultaneously with a time step $\Delta t$ of 6.42 ms (this case corresponds to a ratio of 1), the amplitude of the oscillations deviates after 20 s by only 0.17% from the amplitude of the initial tilt and the period approaches the theoretical non-hydrostatic period $T_{NH}$ by $4.1 \times 10^{-3}$%.

When a ratio of 2 is imposed leading to an internal mode time step of 12.84 ms, the amplitude deviates by 0.11% from its initial value and the period approaches this time $T_{NH}$ by 0.059%. For a time-step ratio of 10, the amplitude changes finally by 2.1% and the period by 0.78%. As could have been expected, the precision is reduced when the internal mode time step is enhanced but the modification of the period and amplitude of the oscillations is rather small but cannot be neglected.

This balance of energy gives some physical insights to understand the mechanisms controlling the free oscillations. To start with, the variation rate of the kinetic energy is clearly associated with the work rate of the surface pressure force which confirms that the free oscillations are due to a transfer of energy between the kinetic and the potential compartments, the potential energy being associated in this case with the variation of the free surface. The total dissipation rate can be precisely evaluated and has been shown to be one million times smaller than the total variation rate. This clearly shows that the evaluation of the mixing rate in the propagation region of such waves must be associated with an accurate, energetically controlled modelling approach. Indeed any spurious numerical energy lost larger than 1 per million times the evolution rate of the kinetic energy prevents any evaluation of the viscous dissipation rate. This is true whether the free oscillations are hydrostatic (not shown) or non-hydrostatic. Indeed, as the work rate of the non-hydrostatic pressure force vanishes, it has no impact on the global kinetic energy balance. It is rather

associated with a redistribution of the energy between the vertical and horizontal components of the kinetic energy.

The evolution of the surface elevation anomaly exhibits close agreement with theory and the curves cannot be distinguished (not shown). Fig. 4 shows the evolution of the RMS errors of the main wave parameters with the grid resolution and the time step. The RMS errors are presented for the surface elevation (a), the horizontal (b) and vertical (c) components of the velocity and the non-hydrostatic pressure (d) against the grid scale. The numerical solutions are compared to the analytical solutions for the same configuration given by Chen (2003, Eqs. (43) and (44)). The grid scale ($\Delta x$) is varied from $(10/2^6)$ to $(10/2^3)$ m. For each grid scale, the RMS errors are shown for several time step spreading from $(\Delta x/2^7)$ to $(\Delta x/2^2)$ s. The RMS errors for the four parameters increase linearly with the grid resolution except for the largest time step $(\Delta x/2^2)$ s. In this case, the RMS errors increase much more quickly and are not shown for $\Delta x$ higher than $(10/2^4)$ m. The order of magnitude of the RMS errors remains rather low but is about 1000 times higher for the non-hydrostatic pressure than for the other three parameters. The evolution of the errors is first-order and indicates accurate behavior of the algorithm together with a rather low level of error.

## 4.2. Formation of nonlinear waves over George Bank (Experiment 2)

The formation of non-linear solitary waves over a steep continental slope is now investigated. The chosen configuration corresponds to Lamb's idealized representation of George Bank (L94) and the bottom topography is consequently given by:

$$h(x) = 65 + 97.5[1 - \tanh(2.2510^{-4}(x - x_0))] \quad (4.4)$$

while the density stratification corresponds to the first profile proposed by Lamb (L94, Eqs. 7a-b). A monochromatic barotropic tidal wave is forced at the eastern boundary with a period $T = 12.4$ h (M2)

<figure><img src="image_5.png" /><figcaption>Fig. 4. RMS Errors for the surface elevation (a), the horizontal (b) and vertical (c) components of the velocity and the non-hydrostatic pressure (d) against the grid scale. The reference is the analytical solution by Chen (2003). Time step is equal to $\Delta x/2^2$ (circles), $\Delta x/2^3$ (crosses), $\Delta x/2^4$ (points), $\Delta x/2^5$ (stars)) and $\Delta x/2^6$ (squares). Dotted line has a slope equal to 1.</figcaption></figure>
and an amplitude corresponding to a volume flux of 59.8 m² s⁻¹. This induces a maximum tidal current of about 92 cm/s over the shelf. At t = 0, no initial field is specified so that an off-bank tidal flow can freely build up in so far as the free surface assumption permits the propagation of surface waves. This however introduces some differences with L94's initialization procedure. First of all, it takes about half an hour for the surface tide to propagate from the eastern boundary to the bank edge. During this period the isopycnal surfaces above the bank remain quasi-still whereas they start to move at the initial time step in L94. No time offset was yet introduced but the stability of the algorithm permits the investigation of at least the first three periods. To speed up the integration, the internal mode time step is chosen to be 6-times as large as the external mode time step. Larger time steps lead to small scale instabilities after a few tidal periods. The amplitude of the Asselin filter is twice the amplitude used for the remaining two Experiments.

Cyclic boundary conditions are implemented in the along-slope direction in order to simulate the along slope currents due to the non-vanishing Coriolis pseudo-force. The main characteristics of this realistic model implementation are summarized in Table 3. The horizontal resolution is 15 m with 100 σ-levels and the cosine component of the Coriolis force is neglected. Unlike in L94, the flow is viscous: a turbulent closure scheme is implemented in the vertical direction and the advection scheme for momentum is centered vertically and upstream horizontally. A rather low background isotropic eddy diffusivity of 10⁻⁵ m²s⁻¹ is used for the temperature and salinity equations. A free-slip boundary condition is implemented at the bottom.

As observed by L94, the first tidal period is a numerical spin-up and we shall show that the dynamics of the internal waves takes a few periods to build up. The dynamically more realistic second and third tidal periods are thus more particularly investigated. Fig. 5a shows the position of the isopycnals at t = 1.625T and can be compared to L94's Fig. 6m while Fig. 5b shows the same isopycnals one period later (i.e. at t = 2.625T). At t = 1.625T first, the large depression formed by the off-bank tidal flow between T and T + T/4 has separated into two depressions labeled A₂ and C₂ propagating respectively, on and off-bank. As far as the on-bank dynamics is concerned, A₂ is about to merge with another on-bank propagating depression (B₁) generated during the second part of first tidal period in agreement with the observations reported by Brickman and Loder (1993). Depression B₁ has been advected back to the bank area by the off-bank flow. Between 10 and 20-m deep, depression A₂ has an amplitude of about 10 m at t = 1.625T and 15 m at t = 2.625T versus about 20 m in L94. At t = 1.625T, depression B₁ and to a lesser extent A₂ have also penetrated 2.5 km further onto

<table><caption>**Table 3**
Model parameters for Georges Bank experiment (Experiment 2).</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Domain size (m)</td><td>(67,500, 45, 260)</td></tr><tr><td>Domain size in points: (i<sub>max</sub>, j<sub>max</sub>, k<sub>max</sub>)</td><td>(4500, 3, 100) Cyclic boundary conditions in the y-direction</td></tr><tr><td>Δx = Δy</td><td>15 m</td></tr><tr><td>External mode time step: Δt<sub>e</sub></td><td>0.0945 s</td></tr><tr><td>Internal mode time step: Δt</td><td>0.567 s</td></tr><tr><td>Viscosity</td><td>(Horizontally: diffusive part of the upstream advective scheme, vertically: turbulent scheme)</td></tr><tr><td>Diffusivity</td><td>10<sup>-5</sup> m<sup>2</sup>s<sup>-1</sup></td></tr><tr><td>Advection scheme</td><td>Upstream</td></tr><tr><td>Asselin filter parameter</td><td>0.2</td></tr><tr><td>Bottom boundary condition</td><td>Free slip</td></tr><tr><td>Coriolis parameter</td><td>f = 10<sup>-4</sup> s<sup>-1</sup></td></tr><tr><td>Equation of state</td><td>Linear</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Cross-isobath vertical sections of the density fields in kg/m² at t = 1.625 T (a) and t = 2.625 T (b). Letters A–C indicate the location of the depressions of the pycnocline and E indicates the position of a positive elevation of the pycnocline, the subscripts indicate at which period these depressions or elevations are generated.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Cross-isobath vertical sections of the horizontal (a) and vertical (b) components of the velocity at *t* = *T*. Contours are drawn every 0.04 m/s for the horizontal velocity and every 0.0003 m/s for the vertical velocity. Letter F indicates the position of a positive elevation of the pycnocline.</figcaption></figure>

the bank than in L94 (Fig. 5a). At $t = 2.625T$, $B_2$ and $A_3$ are however located in the same area (Fig. 5b). These two depressions are separated by a distance close to the 5 km observed by Loder et al. (1992) (after correcting for ship steaming time (L94)): 4.9 km at $t = 1.625T$ and about 4 km at $t = 2.625T$. This also corresponds to the distance found by L94. Deeper in the water column, in the area of the bank edge, the amplitude of the depression reaches 32 m at $t = 1.625T$ and 40 m at $t = 2.625T$ versus 40–45 m in L94. An older, now on-bank propagating depression labeled $A_1$ is located around $x = -28$ km close to the position found by L94 at $t = 1.625T$. At $t = 2.625T$, the corresponding depression $A_2$ is much closer to the bank (between $x = -25$ m and $x = -22.5$ m) and the position of the front is then depth dependent. This “older” depression has been advected back and forth to the bank by the off-bank and on-bank tidal currents and its position consequently appears to be strongly dependent on the strength of the tidal current.

During the second period and over the plain, the depression $C_2$ has evolved into a solibore and is located at $x = 11.8$ km, i.e. 2 km closer to the bank than in L94. It also shows a steeper front than in L94 and it is followed by a smaller number of (smaller amplitude) secondary oscillations. These oscillations still appear to be quickly growing and they only reach an amplitude close to the one of the front between $t = 1.625T$ and $t = 1.750T$ versus $t = 1.375T$ and $t = 1.5T$ in L94. The evolution of this depression during the third period is however similar to the one found in L94 (2.b), and the amplitude of the secondary oscillations is already higher at $t = 2.625T$ than at $t = 1.625T$. This is consistent with the space shuttle observations reported by La Violette et al. (1990) who reported the presence of “several packets of internal waves propagating away from the northern edge of Georges Bank into the Gulf of Maine”. The number and the amplitude of the secondary maxima are more specifically studied in Section 4.3 and they are shown to be very dependent on the internal mode time step.

At $t = 1.625T$, an elevation ($E_2$) is also being generated above the bank edge and is followed by another depression ($B_2$). The slopes of the isopycnals between $A_3$ and $E_3$ are steeper at $t = 2.625T$ and their structure in this area looks like that found by L94. Due to the viscous nature of the flow, the instabilities following depression $B_2$ have been dissipated. A closer inspection of the turbulent dissipation (not shown) shows higher amplitudes in this region confirming the presence of small scale instabilities. The free surface hypothesis also enables the evaluation of the surface signature of the depression (about 2 cm) between the plain and the shelf.

The off-bank horizontal and vertical velocities are given at $t = T$ on Fig. 6. These plots can be compared to L94’s Fig. 7a and b. The structure of the velocity fields is in close agreement with L94. Two vortices can indeed be identified on each side of the bank slope. The vortex located at the off-bank side is anticlockwise while the one at the on-bank side is of opposite sign. These vortices were obviously not present at the beginning of the first tidal period which explains at least part of the differences observed between the first and second tidal period. The strength of the vortexes is somehow larger in L94 with horizontal velocities reaching as much as 40 cm/s close to the surface over the bank while they remain lower than 13 cm/s at $t = T$ (but are still increasing) and only approach 30 cm/s one period latter. Fig. 6 also confirms the Mode-2 structure of depression $B_1$ while depression $A_1$ clearly looks like a Mode-1 wave.

As a conclusion, the proposed algorithm provides in this complex test case very similar dynamics as the non-hydrostatics experiment by L94 confirming, if needed, the dynamics of the non-linear internal waves found in this study, corroborating also the available observations. This result is particularly satisfactory for several reasons associated with the basic assumptions on which the two algorithms are based. The present algorithm is “free surface” while that used by L94 is “rigid lid” which basically means that the dynamics
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** position of the pycnocline and energy transfers in Experiment 3 per unit length in the y-direction. (a) Position of the pycnocline in meters at three locations. Dashed line: at the left wall, Dash-dotted line: at the right wall, Solid line: at the center of the basin. (b) Solid line: variation rate of the kinetic energy (DdKI/dt), (o): work rate of the hydrostatic pressure force, (+): buoyancy flux (Phi z). (c): Solid line: work rate of the horizontal component of the non-hydrostatic pressure force (Nh pgrad, x). (+): minus the work rate of the vertical component of the non-hydrostatic pressure force (Nh pgrad, z). Dashed line: power of the sum of the mixing terms (horizontal and vertical viscous terms and Asselin filter). (d) Dash-dotted line: power of the advection of the horizontal momentum (Adv, x), Solid line: power of the advection of the vertical momentum (Adv, z). Dashed line: power of the Asselin filter.</figcaption></figure>

of the forcing surface tidal wave is not represented in L94 while it is in the present study. No initial field had to be specified and the barotropic tide propagates from the boundaries. This introduces a small offset in the establishment of the dynamics with regard to the configuration by L94. The fields show similar structures but the amplitude of the velocity anomalies are for instance higher during the third period in the present study but do not keep on increasing during the following periods, unlike in L94 (Section 4.2) where "the ensuing large velocities and small-scale features necessitated the termination of the model run".

This is a crucial point as the present model is supposed to be implemented in realistic configurations and as a consequence, includes dissipative processes. However, in spite of these fundamental differences, the models show rather close agreements, the most striking discrepancies being associated with the strength of the various dynamical structures. In the present simulation, the smoother character of the flow (even during the third period) is indeed to be associated with its viscous nature. In the present algorithm, the dissipation takes four different forms. In the momentum equations, vertical viscous diffusion is associated with the turbulent closure scheme (Gaspar et al., 1990) while horizontal diffusion is in the present case restricted to the diffusive part of the upstream advective scheme. The evolution of the mass field (through the density) has a small background diffusivity equal to $10^{-5} \text{ m}^2/\text{s}$. To finish, the Asselin (time) filter has a classic value of 0.2. Several sensitivity tests have been carried out. They show that a reduction of the dissipation results in larger gradients and steeper fronts. To go further the energy transfers associated with these dissipation schemes have to be quantified. Lamb (2007) showed for a rather

similar configuration of non-linear internal waves in the region of Georges Bank that the fluxes of kinetic and available potential energy were of the same order of magnitude as the traditional energy flux radiated by internal waves. This conclusion was reached for a stronger thermocline and consequently larger amplitude internal solitary waves. However, as the available observations do not permit to investigate the details of the various fluxes of energy, the latter are to the present time more advantageously evaluated on the simpler “academic” test cases presented in the following sections. As far as the present configuration is concerned, the choice was thus made to implement a configuration as close as possible to an operational configuration with well-known and well-identified numerical schemes. This means in particular that, unlike in L94, the present model can be integrated for several tidal periods, with the only restriction that the density profile is eroded with time as no restratification mechanism (other than advection through the open boundaries) has been included. Steeper fronts and gradients can then be obtained with the present algorithm by increasing the horizontal and vertical resolutions.

To further investigate the energy transfers and the sensitivity of the properties of non-linear internal waves to the modelling process, we now more specifically focus on the non-linear and non-hydrostatic oscillations of a pycnocline in a closed basin.

### 4.3. Non-linear and non-hydrostatic oscillations of a pycnocline (Experiment 3)

In this third experiment, the formation of interfacial solitons encountered in the previous experiment 2 is revisited in a
two-layer fluid (Horn et al., 2001). The pycnocline is initially submitted to a small perturbation before it freely oscillates during 400 s. The main numerical parameters of the experiment are given in Table 4. The upper layer depth is notably equal to 30% of the total depth while the amplitude of the initial (linear) perturbation is equal to 45% of the upper layer depth. The resulting regime corresponds to the formation of solitons (this corresponds to Regime 2 of Horn et al.). The external and internal mode time steps are equal and the grid resolution is 10 cm with 74 σ-levels. As in the previous experiment, no initial velocity is specified and a no-slip boundary condition is implemented. The viscosity and the diffusivity are also molecular: no turbulent closure scheme is used and the advection schemes are centered.

Fig. 7a shows the evolution with time of the pycnocline depth at three locations: one along each wall and one at the center of the domain. A forming soliton reaches the center of the domain after about 230 s and the amplitude of the soliton is more than twice as large at the points located along the walls due to the reflection. Its formation is thus a little slower than in the observations made by Horn et al. (2001). This difference can be attributed to the initialization procedure. Horn et al. generate indeed the initial oscillation by moving their basin while a linear slope is imposed in the present implementation. Several initial fields have been tested, the specification of a sinusoidal slope at the initial time resulting for instance in a development of the soliton structure in only 120 s, i.e. quicker than in the observations made by Horn et al. The same type of conclusion is reached by Hodges et al. (2006).

Fig. 7b-d show the evolution of the non-vanishing terms of the kinetic energy equation. Three groups of terms are of particular interest for the present study: the work rate of the pressure force (both hydrostatic and non-hydrostatic), the power of the horizontal and vertical advections and the power of dissipation.

During the first minute of simulation, the variation rate of the kinetic energy (Fig. 7b) is approximately equal to the power of the hydrostatic pressure force and the remainder is to be attributed to the work rate of the surface pressure force and for a minor part to the dissipation rate. The curve representing the hydrostatic pressure force cannot be distinguished from the one representing the buoyancy flux. The dates when the front of the oscillation reaches one of the walls can easily be observed and correspond to large oscillations of the rate of change of the kinetic energy upon which are superimposed the high-frequency oscillations of the solitons after the fourth reflection. The buoyancy flux is responsible for the oscillatory exchange of kinetic and potential energy and, in the case of a closed basin it approximately cancels the power of the advection of potential energy.

The horizontal and vertical components of the non-hydrostatic pressure force (Fig. 7c) exactly cancel each other. After about

260 s, the non-hydrostatic contribution to the pressure force undergoes large amplitude oscillations right after the reflection of the wave front at one of the walls confirming the importance of the non-hydrostatic dispersion, the amplitude being particularly large when the amplitude of the solitons is large, i.e. at the walls. This clearly confirms the dispersive nature of the non-hydrostatic pressure contribution. At the same time, the work rate of advection (Fig. 7d) is also larger each time the front of the interfacial oscillation reaches a wall, but it remains about 1000 times smaller than the rate change of the kinetic energy. This periodic increase corresponds to the non-linear steepening of the front. Fig. 7c and d consequently provide a clear proof of the presence of solitary waves associated with a non-linear steepening of the wave front and to a non-hydrostatic dispersion of this front.

The energy rate of change due to the mixing terms, i.e. the sum of the horizontal and vertical viscous diffusions and of the Asselin filter, is at most one order of magnitude smaller than the rate of change of the kinetic energy and shows larger values when the front of the oscillations hits the walls. It is to be noted that the amplitude of the dissipation fluxes is reduced when the soliton are formed. The dissipation associated with the Asselin filter is also shown to contribute for less than one percent to the total mixing rate (Fig. 7d). After a few seconds, these terms lower to about 0.01% of the total rate of change of the kinetic energy.

The same experiment is finally presented for varying ratios of the external to the internal mode time-steps and Fig. 8 shows the resulting evolution of the position of the pycnocline at a point located at the center of the basin for $t \in [280 \text{ s}, 350 \text{ s}]$. In all these sensitivity tests, the external mode time-step remains equal to 0.0189 s and, to facilitate the comparison, the viscous diffusion schemes are similar to those used in Georges Bank experiment (vertically: turbulent scheme, horizontally: background diffusivity of $10^{-5} \text{ m}^2\text{s}^{-1}$ and diffusive part of the upstream advective schemes). All the curves show the soliton-like evolution of the pycnocline with exactly the same periodicity. Before the appearance of the first soliton, the curves for the different ratios cannot be distinguished (not shown) but the representation of the solitons differs afterwards: the number of secondary maxima and their amplitude are indeed drastically reduced when the time-step ratio is increased. With a ratio greater or equal than 8 only 3 maxima are present. If this ratio is decreased to 4, one additional secondary maximum can be modeled while a maximum ratio of 2 is in particular needed to accurately simulate the 5th maximum. Even if no direct comparison with the time step used by L94 is possible (this information is not available for a direct comparison), this explains the lower amplitude of the secondary maxima in Georges Bank experiment which is integrated with a time-step ratio of 6. This somehow restrains the use of the mode-splitting algorithm when small scale solitons are present in the domain and an embedded modeling strategy must then be implemented (Auclair et al., 2006).

A complete plot of the various components of the kinetic energy balance provides thus both a qualitative and quantitative description of the non-linear steepening and non-hydrostatic dispersion mechanisms associated with the formation of solitary internal waves. As a comparison, it is then interesting to examine the same simulation under the hydrostatic assumption (Fig. 9). Not surprisingly, the first 250 s are rather similar and the buoyancy flux cancels approximately the hydrostatic pressure force (Fig. 9b). However, some oscillations looking like solitary waves with a much higher frequency appear between 250 s and 300 s. These oscillations are obviously not sustained by a dispersive non-hydrostatic mechanism and only mimic the oscillations due the propagation of a soliton. The order of magnitude of the mixing rate is rather similar to the non-hydrostatic case and the contribution of the Asselin filter remains low.

<table><caption>**Table 4**
Model parameters for Experiment 3.</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Domain size (m)</td><td>(6, 0.29)</td></tr><tr><td>Domain size in points: (i<sub>max</sub>, k<sub>max</sub>)</td><td>(60, 74)</td></tr><tr><td>Δx = Δy</td><td>10 cm</td></tr><tr><td>External mode time step: Δt<sub>e</sub></td><td>0.018 s</td></tr><tr><td>Internal mode time step: Δt</td><td>0.018 s</td></tr><tr><td>Viscosity</td><td>10<sup>-6</sup> m<sup>2</sup>s<sup>-1</sup></td></tr><tr><td>Diffusivity</td><td>10<sup>-7</sup> m<sup>2</sup>s<sup>-1</sup></td></tr><tr><td>Advection scheme</td><td>Centered</td></tr><tr><td>Asselin filter parameter</td><td>0.1</td></tr><tr><td>Bottom boundary condition</td><td>No slip</td></tr><tr><td>Horn's Experiment Δ parameter</td><td>0.45</td></tr><tr><td>Horn's Experiment γ parameter</td><td>0.3</td></tr><tr><td>Coriolis parameter</td><td>0</td></tr><tr><td>Equation of state</td><td>Linear</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** position of the pycnocline at the center of basin for $t \in [280s, 350s]$. The evolution with time of the position of the pycnocline is given for several external to internal mode time-step ratios ($r = \Delta t_e/\Delta t$). Ratios of 1, 2, 4, 6, 8 and 10 are plotted.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** position of the pycnocline and energy transfers in the hydrostatic implementation of Experiment 3 per unit length in the y-direction. (a) Position of the pycnocline in meters at three locations. Dashed line: at the left wall, Dash-dotted line: at the right wall, Solid line: at the center of the basin. (b) Solid line: variation rate of the kinetic energy (DdKI/dt), (o): work rate of the hydrostatic pressure force, (+): buoyancy flux (Phi z). (c): Solid line: Dashed line: power of the sum of the mixing terms (horizontal and vertical viscous terms and Asselin filter). (d) Dash-dotted line: power of the advection of the horizontal momentum (Adv, x), Solid line: power of the advection of the vertical momentum (Adv, z). Dashed line: power of the Asselin filter.</figcaption></figure>

Interestingly enough, it can be noted that since the work rates of the horizontal and vertical non-hydrostatic pressure forces cancel, the hydrostatic and non-hydrostatic balances of the kinetic energy look misleadingly rather similar and the contribution of the advection terms remains low. The non-linear exchange of energy

they are responsible for indeed only equal to a few percents of the exchange of energy. This confirms, if needed, the importance of the energy balances.

These results can be compared to the modelling of the same experiments by Hodges et al. (2006). They indeed evaluate
indirectly the dissipation of kinetic and potential energy in order to assess numerical errors. In the present work the various contributions to the energy transfers are separated. The main advantage of such quantitative evaluations is to investigate in detail the dynamics of the studied mechanisms and, if needed, to further control them.

## 5. Discussion and conclusion

A Boussinesq, free-surface, non-hydrostatic, mode-splitting algorithm in $\sigma$-coordinate has been proposed for ocean modelling. This algorithm is designed to be easily implemented as a non-hydrostatic kernel for originally hydrostatic ocean models. The "pressure correction" kernel exhibits several original aspects such as the implementation of the free surface in $\sigma$-coordinates and the treatment of the vertical velocity in the bottom and surface layers. The dynamical treatment of the free surface and its coupling with the non-hydrostatic pressure is somehow similar to the algorithms proposed by Casulli (1999) and Zijlema and Stelling (2005). A noteworthy difference with these studies is the treatment of the free surface which is integrated using a time-splitting explicit scheme. This enables the propagation of non-linear waves without requiring prohibitively small time steps. The vertical velocity is integrated using the vertical momentum equation all to way to the surface and, at the lower level, the kinematic condition is used as a boundary condition. No incremental vertical or horizontal velocity is computed along this bottom layer and the surface kinematic condition can be implemented as a consequence of the depth-integrated continuity equation and of this bottom kinematic condition. The projection operator (associated with the non-hydrostatic correction) consequently commutes with the bottom boundary conditions. Another original aspect is the strict conservation of the energy fluxes. In this respect, the proposed approach is a non-hydrostatic extension of the conservative algorithm detailed by M08. We have additionally shown in the present study that the discretization of the work rate of the non-hydrostatic pressure gradient should vanish in a closed domain which imposes a simple extrapolation scheme for the horizontal velocity at the top and at the bottom of the water column. Simple low-order numerical schemes have been favored. The formulation of the $\sigma$-correction term of the continuity equation (3.1) leads to natural expressions of the bottom and surface kinematic conditions (respectively, Eqs. (3.13) and (3.14)). Once this first scheme is chosen, energy conservation imposes the discretization scheme for the non-hydrostatic pressure gradient.

To illustrate the properties of the algorithm, the propagation of short surface waves and non-linear internal waves in closed boxes have been studied in Experiments (1) and (3). For both experiments, a complete energy balance has been carried out and the evolution of the main energy transfers has been presented. In addition to the numerical confirmation of some well-known energy transfers, several original results have been obtained. The study of the energy transfers associated with the Asselin filter has further confirmed that such a filter may not only dissipate but also increase the kinetic energy (M08). In the third experiment it was also possible to quantify the power of the advection and of the non-hydrostatic dispersion in the propagation of non-linear solitary waves. It has finally been shown that a large increase of the external to internal time step ratio associated with an increase of the internal mode time step can lead to a reduction of the number of secondary maxima and to a lowering of their amplitude.

The algorithm has also been implemented to model the generation and the propagation of large amplitude, non-linear internal tides in the region of Georges Bank (Experiment 2). The dynamics of the non-linear waves were similar to the observations and

numerical modelling reported in L94. This model provided indeed an interesting comparison since it was based on rather different assumptions and technical choices: it is in particular rigid-lid and inviscid. The characteristics and the position of the depressions and elevations of the isopycnals are in close agreement. The generation and splitting of the depression together with the off-bank propagation of a solibore have also been confirmed. The major differences between the present modelling and that of L94 are associated with the inviscid nature of the flow in the latter configuration and to a lesser extent to the explicit modelling of the free-surface dynamics. Although the absence of any dissipation was indeed interesting to investigate the generation and the propagation of internal waves, it did not prevent the generation and the amplification of small-scale instabilities which limits the length of the run in L94. In contrast, the present configuration was run with a realistic viscous, diffusive and turbulent algorithm. The dissipation schemes clearly avoid the accumulation of energy at small scales. An associated consequence is the limited and controlled smoothing of the main structures, while steeper gradients are obtained only by increasing the horizontal and vertical resolution.

## Acknowledgments

We gratefully thank the computer team of the Laboratoire d'Aérologie for its support. This research was supported by the LEFE-IDAO "Ondes internes" and the ANR-PIWO (08-BLAN-0113) French National programs.

## Appendix A. $\sigma$-coordinate system

Two expressions are worth being acknowledged in $\sigma$-coordinates when dealing with horizontal derivatives. The "gradient" formulation is first given by:

$$ \frac{\partial q}{\partial x_z} = \frac{\partial q}{\partial x_z}_\sigma - \frac{1}{D} \frac{\partial z}{\partial x_z}_\sigma \frac{\partial q}{\partial \sigma} \quad (A.1) $$

and the "flux" formulation can then be written:

$$ \frac{\partial v_z}{\partial x_z}_z = \frac{1}{D} \frac{\partial D v_z}{\partial x_z}_\sigma - \frac{1}{D} \frac{\partial}{\partial \sigma} \left( \frac{\partial z}{\partial x_z} \Big|_{\sigma} \right) \quad (A.2) $$

The time derivative in $\sigma$-coordinates is given by: $\frac{\partial q}{\partial t}|_z = \frac{\partial q}{\partial t}|_\sigma - \frac{\sigma}{D} \frac{\partial q}{\partial \sigma} \frac{\partial z}{\partial \sigma}$ and the vertical derivative of the vertical velocity satisfies the relations:

$$ \frac{\partial w}{\partial \sigma} - \frac{\partial}{\partial \sigma} \left[ v_z \frac{\partial z}{\partial x_z} \Big|_{\sigma} \right] = \frac{\partial v_z}{\partial \sigma} + \frac{\partial \zeta}{\partial t} \quad (A.3) $$

The formulation of the continuity equation used in the hydrostatic version of the model (M08) can be then be recovered from (2.8):

$$ \frac{\partial D v_z}{\partial x_z} \Big|_{\sigma} + \frac{\partial v_z}{\partial \sigma} + \frac{\partial \zeta}{\partial t} = 0 \quad (A.4) $$

The "true" vertical velocity (w) being the non-hydrostatic model prognostic variable, Eq. (2.8) rather than (A.4) is implemented. The $\sigma$-vertical velocity ($v_\sigma$) is a diagnostic variable and is only used to compute the vertical component of the advection of momentum and tracers. It can be derived either from (A.3) or (A.4).

## Appendix B. Momentum equations

The horizontal and vertical provisional momentum equations are respectively given by (3.15) and (3.16) (the bottom vertical velocity ($w_{i\frac{1}{2}}^{\pm\Delta t}$ = $w_{i\frac{1}{2}}^{\pm\Delta t}$) is computed through the kinematic condition (3.13)):
$$
\text{For } k \in [1, k_{\text{max}}]: \left[ \Delta z \tilde{u} \right]_{i+\frac{1}{2}, k}^{t+\Delta t} = \left[ \Delta z u \right]_{i+\frac{1}{2}, k}^{t-\Delta t} + 2\Delta t \left( -\frac{\Delta z_t}{\rho_0} \frac{\partial q_t}{\partial x} \right)_{z, \frac{1}{2}, k}^t + f_{u,i+\frac{1}{2},k}^t \right) \quad (3.15)
$$

$$
\text{For } k \in [1, k_{\text{max}}]: \left[ \Delta z \tilde{w} \right]_{i+\frac{1}{2}, k}^{t+M} = \left[ \Delta z w \right]_{i+\frac{1}{2}, k}^{t-\Delta t} + 2\Delta t \left( -\frac{\Delta z_t}{\rho_0} \frac{\partial q_t}{\partial z} \right)_{z, i+\frac{1}{2}}^t + f_{w,i+\frac{1}{2}, k}^t \quad (3.16)
$$

where $f_u$ and $f_w$ include the momentum advection and viscosity to-
gether with the horizontal hydrostatic pressure and Coriolis forces.
The discretization of the former operator in the horizontal momen-
tum equation is detailed in M08 except for the non-traditional Cori-
olis contribution given by:

$$
-\frac{1}{4} \left( f'_{i} \Delta z_{i,k}^{t} \left( w_{i,k+\frac{1}{2}}^{t} + w_{i,k-\frac{1}{2}}^{t} \right) + f_{i+1}^{t} \Delta z_{i+1,k}^{t} \left( w_{i+1,k+\frac{1}{2}}^{t} + w_{i+1,k-\frac{1}{2}}^{t} \right) \right) \quad (B.1)
$$

where $f' = 2\Omega\cos\phi$ and $\phi$ is the latitude.

In the vertical momentum equation, the operator $(f_t^w)_{i,k+\frac{1}{2}}$ is the sum of the following contributions:

* Advection of vertical momentum (a centered scheme is implemented for the present study): $-(\text{Adv}v_{i+\frac{1}{2},k}^{H,t} - \text{Adv}v_{i-\frac{1}{2},k}^{H,t} + \text{Adv}v_{i,k+1}^{V,t} - \text{Adv}v_{i,k-1}^{V,t})$ where the horizontal component of the advective flux is:

$$
\text{Adv}_{i+\frac{1}{2},k+\frac{1}{2}}^{H,t} = \frac{1}{4\Delta x} \left( [\Delta z u]_{i+\frac{1}{2},k+1}^{t} + [\Delta z u]_{i+\frac{1}{2},k}^{t} \right) \left( w_{i+1,k+\frac{1}{2}}^{t} + w_{i,k+\frac{1}{2}}^{t} \right) \quad \text{for } k \in [1, k_{\text{max}}] \quad (B.2)
$$

$$
\text{Adv}_{i+\frac{1}{2},k+\frac{1}{2}}^{H,t} = \frac{1}{2\Delta x} [\Delta z u]_{i+\frac{1}{2},k_{\text{max}}}^{t} (w_{i+1,k+\frac{1}{2}}^{t} + w_{i,k+\frac{1}{2}}^{t}) \quad \text{for } k=k_{\text{max}} \qquad (B.3)
$$

and the vertical component is:

$$
\text{Adv}_{i,k}^{V,t} = \frac{1}{4} (\omega_{i,k+\frac{1}{2}}^t + \omega_{i,k-\frac{1}{2}}^t) (w_{i,k+\frac{1}{2}}^t + w_{i,k-\frac{1}{2}}^t) \quad \text{for } k \in [1, k_{\text{max}}] \qquad (B.4)
$$

$$
\mathrm{Adv}_{i,k_{\mathrm{max}}+1}^{V,t} = 0 \quad \text{for } k = k_{\mathrm{max}} + 1 \qquad (\mathrm{B}.5)
$$

* Viscous dissipation of vertical momentum: $Visc_{i+\frac{1}{2},k+\frac{1}{2}}^{H,t} - Visc_{i-\frac{1}{2},k+\frac{1}{2}}^{H,t} + Visc_{i,k+1}^{V,t} - Visc_{i,k}^{V,t}$ where the horizontal component of the viscous flux is:

$$
\text{Visc}_{i+\frac{1}{2},k+\frac{1}{2}}^{H,t} = K_{i+\frac{1}{2},k+\frac{1}{2}}^{H} \frac{\Delta z_t}{\Delta x^2} \left( w_{i+1,k+\frac{1}{2}}^t - w_{i,k+\frac{1}{2}}^t \right) \quad \text{for } k \in [1, k_{\text{max}}] \qquad (\text{B.6})
$$

and the vertical component is:

$$
\text{Visc}_{i,k}^{V,t} = \frac{K_{i,k}^{V}}{\Delta z_{i,k}^{t}} (\tilde{w}_{i,k+\frac{1}{2}}^{t+\Delta t} - \tilde{w}_{i,k-\frac{1}{2}}^{t+\Delta t}) \quad \text{for } k \in [1, k_{\text{max}}] \quad (B.7)
$$

As in the horizontal direction, only the vertical component of the
viscous dissipation scheme is implicit.

* Non-traditional Coriolis pseudo-force:

$$
\frac{1}{4} \left( f'_{i} \Delta z_{i,k}^{t} \left( u_{i+\frac{1}{2},k}^{t} + u_{i-\frac{1}{2},k}^{t} \right) + f'_{i} \Delta z_{i,k+1}^{t} \left( u_{i+\frac{1}{2},k+1}^{t} + u_{i-\frac{1}{2},k+1}^{t} \right) \right) \quad \text{for } k \in [1, k_{\text{max}}] \qquad (B.8)
$$

and

$$
\frac{1}{4} \left( f'_{i} \Delta z_{i,k_{\text{max}}}^{t} \left( u_{i+\frac{1}{2},k_{\text{max}}}^{t} + u_{i-\frac{1}{2},k_{\text{max}}}^{t} \right) \right) \quad \text{for } k = k_{\text{max}} \qquad (B.9)
$$

**Appendix C. Poisson operator**

A sparse-matrix formulation of the non-hydrostatic algorithm
has been implemented for several reasons. The “continuity” and
“non-hydrostatic pressure force” matrix operators are used not only
to construct the Laplacian operator appearing in the non-hydro-

static Poisson equation (3.23) but also to compute respectively the
right-hand side of this equation and the non-hydrostatic velocity
increment. The non-hydrostatic pressure gradient operator is also
used to compute the non-hydrostatic pressure force in the barotrop-
ic (3.17) and baroclinic momentum equations (3.15) and (3.16). On
top of the possible reduction of computing time, the matrix imple-
mentation leads de facto to a perfect numerical compatibility be-
tween the continuity equations and Laplacian operator.

At time step *t*, the non-hydrostatic pressure and velocity vectors for the region of interest (*R*) can be defined respectively by
$$Q^t = \{q_{i,k}^t\}_{(i,k) \in \mathcal{R}} \quad \text{and} \quad V^t = \{u_{i+\frac{1}{2},k}^t, w_{i,k+\frac{1}{2}}^t\}_{(i,k) \in \mathcal{R}}$$

C.1. Continuity equation

At time step $t + \Delta t$, the continuity equation (3.1) can be
rewritten:

$$
C^{t+\Delta t} V^{t+\Delta t} = 0 \qquad (C.1)
$$

In the following, the time indexes are all related to time step $t + \Delta t$ and are dropped to simplify the notations. At point $(i, k)$ and for $k \neq 1$ (strictly above the bottom layer) this leads to:

$$
c_1 u_{l,k} + c_2 u_{l-1,k} + c_3 u_{l,k+1} + c_4 u_{l-1,k+1} + c_5 u_{l,k-1} + c_6 u_{l-1,k-1} \\
+ c_7 w_{l,k+1/2} + c_8 w_{l,k-1/2} = 0 \qquad (C.2a)
$$

where $I = i + 1/2$. With $\alpha_{l,k}$ and $\beta_{l,k}$ defined by Eqs. (3.4) and (3.5), the coefficients are given by:

* for $k \in [1, k_{\text{max}}]$:

$$
c_1 = \Delta z_{l,k} / \Delta x + \alpha_{l,k} - \beta_{l,k}
$$

$$
c_2 = -\Delta z_{l-1,k} / \Delta x + \alpha_{l-1,k} - \beta_{l-1,k}
$$

$$
c_3 = -\alpha_{l,k+1}
$$

$$
c_4 = -\alpha_{l-1,k+1} c_5 = \beta_{l,k-1}
$$

$$
c_6 = \beta_{l-1,k-1} \quad c_7 = 1 \quad c_8 = -1
$$

* for $k = k_{\text{max}}$:

$$
c_1 = \Delta z_{l,k_{\text{max}}}/\Delta x + \alpha_{l,k_{\text{max}}}- \beta_{l,k_{\text{max}}+1/2}
$$

$$
c_2 = -\Delta z_{l-1,k_{\text{max}}}/\Delta x + \alpha_{l-1,k_{\text{max}}}- \beta_{l-1,k_{\text{max}}+1/2}
$$

$$
c_5 = \beta_{l,k_{\text{max}}-1}
$$

$$
c_6 = \beta_{l,k_{\text{max}}-1}
$$

$c_7 = 1$

$c_8 = -1$

$c_3 = c_4 = 0$

* For $k=1$: the continuity equation is given by:

$$
c_1 u_{l,1} + c_{1b} u_{l,1/2} + c_2 u_{l-1,1} + c_{2b} u_{l-1,1} + c_3 u_{l,2} + c_4 u_{l-1,2} \\
+ c_7 w_{l,3/2} + c_8 w_{l,1/2} = 0 \qquad (C.2b)
$$

where

$c_1 = \Delta z_{l,1} / \Delta x - \beta_{l,1}$ and $c_{1b} = \alpha_{l,1/2}$

$c_2 = -\Delta z_{l-1,1} / \Delta x - \beta_{l-1,1}$ and $c_{2b} = \alpha_{l-1,1/2}$

$c_3 = -\alpha_{l,2}$

$c_4 = -\alpha_{l-1,2}$

$c_5 = c_6 = 0$

$c_7 = 1 c_8 = -1$

C.2. Non-hydrostatic pressure force

The non-hydrostatic velocity increments (3.21) and (3.22) can
also be rewritten in matrix notations as:
$$ \delta V = - \frac{2 \Delta t}{\rho_0} G^{t+\Delta t} \delta Q \qquad (C.3) $$

For $k \in [1, k_{max}]$, this leads to:

$$ \delta u_{l,k} = - \frac{2 \Delta t}{\rho_0} (g_1 \delta q_{i,k} + g_2 \delta q_{i+1,k} + g_3 \delta q_{i,k+1} + g_4 \delta q_{i+1,k+1} + g_5 \delta q_{i,k-1} + g_6 \delta q_{i+1,k-1}) \qquad (C.4) $$

$$ \delta w_{i,k+1/2} = - \frac{2 \Delta t}{\rho_0} (g_7 \delta q_{i,k+1} + g_8 \delta q_{i,k}) \qquad (C.5) $$

The coefficients are given by:

* for $k \in [1, k_{max}]$:

$$ g_1 = -1/\Delta x + (\beta_{1,k} - \alpha_{1,k})/\Delta z_{1,k} \quad g_2 = 1/\Delta x + (\beta_{1,k} - \alpha_{1,k})/\Delta z_{1,k} $$

$$ g_3 = -\beta_{1,k}/\Delta z_{1,k} \quad g_4 = -\beta_{1,k}/\Delta z_{1,k} \quad g_5 = \alpha_{1,k}/\Delta z_{1,k} \quad g_6 = \alpha_{1,k}/\Delta z_{1,k} $$

$$ g_7 = 1/\Delta z_{i,k+1/2} \quad g_8 = -1/\Delta z_{i,k+1/2} $$

* for $k = k_{max}$:

$$ g_1 = -1/\Delta x + (\beta_{1,k_{max}+1/2} - \alpha_{1,k_{max}})/\Delta z_{1,k_{max}} $$

$$ g_2 = 1/\Delta x + (\beta_{1,k_{max}+1/2} - \alpha_{1,k_{max}})/\Delta z_{1,k_{max}} $$

$$ g_5 = \alpha_{1,k_{max}}/\Delta z_{1,k_{max}} \quad g_6 = \alpha_{1,k_{max}}/\Delta z_{1,k_{max}} \quad g_7 = 1/\Delta z_{i,k_{max}+1/2} $$

$$ g_8 = -1/\Delta z_{i,k_{max}+1/2} \quad g_9 = g_4 = 0 $$

* For $k = 1$:

$$ g_1 = -1/\Delta x + \beta_{1,1}/\Delta z_{1,1} \quad g_2 = 1/\Delta x + \beta_{1,1}/\Delta z_{1,1} $$

$$ g_3 = -\beta_{1,1}/\Delta z_{1,1} \quad g_4 = -\beta_{1,1}/\Delta z_{1,1} \quad g_7 = 1/\Delta z_{i,3/2} \quad g_8 = -1/\Delta z_{i,3/2} $$

$$ g_5 = g_6 = 0 $$

### C.3. Poisson equation

In matrix notation, the Poisson system is obtained by multiplying the continuity and pressure gradient matrix operators:

$$ C^{t+\Delta t} \delta Q = \frac{\rho_0}{2\Delta t} C^{t+\Delta t} \tilde{V}^{t+\Delta t} \qquad (C.6) $$

At point $(i, k)$ and for $k \neq 1$ (strictly above the bottom layer) this leads to:

$$ \begin{aligned} a_1 \delta q_{i,k} + a_2 \delta q_{i+1,k} + a_3 \delta q_{i-1,k} + a_4 \delta q_{i,k+1} + a_5 \delta q_{i+1,k+1} & \\ + a_6 \delta q_{i-1,k+1} + a_7 \delta q_{i,k-1} + a_8 \delta q_{i-1,k-1} + a_9 \delta q_{i-1,k-1} & \\ + a_{10} \delta q_{i,k+2} + a_{11} \delta q_{i+1,k+2} + a_{12} \delta q_{i-1,k+2} + a_{13} \delta q_{i,k-2} & \\ + a_{14} \delta q_{i+1,k-2} + a_{15} \delta q_{i-1,k-2} & \\ = \frac{\rho_0}{2\Delta t} \left( c_1 \tilde{u}_{i,k}^{t+\Delta t} + c_2 \tilde{u}_{i-1,k}^{t+\Delta t} + c_3 \tilde{u}_{i,k+1}^{t+\Delta t} + c_4 \tilde{u}_{i-1,k+1}^{t+\Delta t} + c_5 \tilde{u}_{i,k-1}^{t+\Delta t} \right. & \\ \left. + c_6 \tilde{u}_{i-1,k-1}^{t+\Delta t} + c_7 \tilde{w}_{i,k+1/2}^{t+\Delta t} + c_8 \tilde{w}_{i,k-1/2}^{t+\Delta t} \right) \end{aligned} \qquad (C.7) $$

where

$$ \begin{aligned} a_1 = & -\frac{\Delta z_{1,k}}{\Delta x^2} - \frac{2\alpha_{1,k}}{\Delta x} + \frac{2\beta_{1,k}}{\Delta x} + \frac{2\alpha_{1,k}\beta_{1,k}}{\Delta z_{1,k}} - \frac{\alpha_{1,k}^2}{\Delta z_{1,k}} - \frac{\beta_{1,k}^2}{\Delta z_{1,k}} - \frac{\delta z_{1,k-1}}{\Delta x^2} \\ & + \frac{2\alpha_{1-1,k}}{\Delta x} - \frac{2\beta_{1-1,k}}{\Delta x} + \frac{2\alpha_{1-1,k}\beta_{1-1,k}}{\Delta z_{1-1,k}} - \frac{\alpha_{1-1,k}^2}{\Delta z_{1-1,k}} - \frac{\beta_{1-1,k}^2}{\Delta z_{1-1,k}} - \frac{\alpha_{1-1,k}^2}{\Delta z_{1,k+1}} \\ & - \frac{\alpha_{1-1,k+1}^2}{\Delta z_{1-1,k+1}} - \frac{\beta_{1-1,k+1}^2}{\Delta z_{1-1,k+1}} - \frac{\alpha_{1-1,k}^2}{\Delta z_{1,k+1}} - \frac{\beta_{1-1,k+1}^2}{\Delta z_{1,k+1}} - \frac{1}{\Delta z_{i,k+1/2}} - \frac{1}{\Delta z_{i,k-1/2}} \end{aligned} $$

$$ a_2 = \frac{\Delta z_{1,k}}{\Delta x^2} + \frac{2\alpha_{1,k}\beta_{1,k}}{\Delta z_{1,k}} - \frac{\alpha_{1,k}^2}{\Delta z_{1,k}} - \frac{\beta_{1,k}^2}{\Delta z_{1,k}} - \frac{\alpha_{1,k}^2}{\Delta z_{1,k+1}} - \frac{\beta_{1,k}^2}{\Delta z_{1,k-1}} - \frac{1}{\Delta z_{i,k+1/2}} - \frac{1}{\Delta z_{i,k-1/2}} $$

$$ a_3 = \frac{\delta z_{1,k}}{\Delta x^2} + \frac{2\alpha_{1,k}\beta_{1,k}}{\Delta z_{1,k}} - \frac{\alpha_{1,k}^2}{\Delta z_{1,k}} - \frac{\beta_{1,k}^2}{\Delta z_{1,k}} - \frac{\alpha_{1,k}^2}{\Delta z_{1,k+1}} - \frac{\beta_{1,k}^2}{\Delta z_{1,k-1}} - \frac{1}{\Delta z_{i,k+1/2}} - \frac{1}{\Delta z_{i,k-1/2}} $$

$$ \begin{aligned} a_4 = & -\frac{\beta_{1,k}}{\Delta x} - \frac{\alpha_{1,k}\beta_{1,k}}{\Delta z_{1,k}} + \frac{\beta_{1,k}^2}{\Delta z_{1,k}} + \frac{\beta_{1-1,k}}{\Delta x} - \frac{\alpha_{1-1,k}\beta_{1-1,k}}{\Delta z_{1-1,k}} + \frac{\beta_{1-1,k}^2}{\Delta z_{1-1,k}} + \frac{\alpha_{1,k+1}}{\Delta x} \\ & - \frac{\alpha_{1,k+1}\beta_{1,k+1}}{\Delta z_{1,k+1}} + \frac{\alpha_{1,k+1}^2}{\Delta z_{1,k+1}} - \frac{\alpha_{1-1,k+1}\beta_{1-1,k+1}}{\Delta z_{1-1,k+1}} + \frac{\alpha_{1-1,k+1}^2}{\Delta z_{1-1,k+1}} \\ & + \frac{1}{\Delta z_{i,k+1/2}} \end{aligned} $$

$$ a_5 = -\frac{\beta_{1,k}}{\Delta x} - \frac{\alpha_{1,k}\beta_{1,k}}{\Delta z_{1,k}} + \frac{\beta_{1,k}^2}{\Delta z_{1,k}} - \frac{\alpha_{1,k+1}}{\Delta x} - \frac{\alpha_{1,k+1}\beta_{1,k+1}}{\Delta z_{1,k+1}} + \frac{\alpha_{1,k+1}^2}{\Delta z_{1,k+1}} $$

$$ a_6 = \frac{\beta_{1-1,k}}{\Delta x} - \frac{\alpha_{1-1,k}\beta_{1-1,k}}{\Delta z_{1-1,k}} + \frac{\beta_{1-1,k}^2}{\Delta z_{1-1,k}} + \frac{\alpha_{1-1,k+1}}{\Delta x} - \frac{\alpha_{1-1,k+1}\beta_{1-1,k+1}}{\Delta z_{1-1,k+1}} \\ + \frac{\alpha_{1-1,k+1}^2}{\Delta z_{1-1,k+1}} $$

$$ \begin{aligned} a_7 = & \frac{\alpha_{1,k}}{\Delta x} - \frac{\alpha_{1,k}\beta_{1,k}}{\Delta z_{1,k}} + \frac{\alpha_{1,k}^2}{\Delta z_{1,k}} - \frac{\alpha_{1-1,k}}{\Delta x} - \frac{\alpha_{1-1,k}\beta_{1-1,k}}{\Delta z_{1-1,k}} + \frac{\alpha_{1-1,k}^2}{\Delta z_{1-1,k}} - \frac{\beta_{1,k-1}}{\Delta x} \\ & - \frac{\alpha_{1,k-1}\beta_{1,k-1}}{\Delta z_{1,k-1}} + \frac{\beta_{1,k-1}^2}{\Delta z_{1,k-1}} - \frac{\beta_{1-1,k-1}}{\Delta x} - \frac{\alpha_{1-1,k-1}\beta_{1-1,k-1}}{\Delta z_{1-1,k+1}} + \frac{\beta_{1-1,k-1}^2}{\Delta z_{1-1,k-1}} \\ & + \frac{1}{\Delta z_{i,k-1/2}} \end{aligned} $$

$$ a_8 = \frac{\alpha_{1,k}}{\Delta x} - \frac{\alpha_{1,k}\beta_{1,k}}{\Delta z_{1,k}} + \frac{\alpha_{1,k}^2}{\Delta z_{1,k}} + \frac{\beta_{1,k-1}}{\Delta x} - \frac{\alpha_{1,k-1}\beta_{1,k-1}}{\Delta z_{1,k-1}} + \frac{\beta_{1,k-1}^2}{\Delta z_{1,k-1}} $$

$$ a_9 = -\frac{\alpha_{1-1,k}}{\Delta x} - \frac{\alpha_{1-1,k}\beta_{1-1,k}}{\Delta z_{1-1,k}} + \frac{\alpha_{1-1,k}^2}{\Delta z_{1-1,k}} - \frac{\beta_{1-1,k-1}}{\Delta x} - \frac{\alpha_{1-1,k-1}\beta_{1-1,k-1}}{\Delta z_{1-1,k-1}} \\ + \frac{\beta_{1-1,k-1}^2}{\Delta z_{1-1,k-1}} $$

$$ a_{10} = \frac{\alpha_{1,k+1}\beta_{1,k+1}}{\Delta z_{1,k+1}} + \frac{\alpha_{1-1,k+1}\beta_{1-1,k+1}}{\Delta z_{1-1,k+1}} \quad a_{11} = \frac{\alpha_{1,k+1}\beta_{1,k+1}}{\Delta z_{1,k+1}} $$

$$ a_{12} = \frac{\alpha_{1-1,k+1}\beta_{1-1,k+1}}{\Delta z_{1-1,k+1}} \quad a_{13} = \frac{\alpha_{1,k-1}\beta_{1,k-1}}{\Delta z_{1,k-1}} + \frac{\alpha_{1-1,k-1}\beta_{1-1,k-1}}{\Delta z_{1-1,k-1}} $$

$$ a_{14} = \frac{\alpha_{1,k-1}\beta_{1,k-1}}{\Delta z_{1,k-1}} \quad a_{15} = \frac{\alpha_{1-1,k-1}\beta_{1-1,k-1}}{\Delta z_{1-1,k-1}} $$

The coefficients have a special form at the surface (for $k = k_{max}$) (Sections C.1 and C.2).

At the bottom ($k = 1$):

$$ \begin{aligned} a_1 \delta q_{i,1} + a_2 \delta q_{i+1,1} + a_3 \delta q_{i-1,1} + a_4 \delta q_{i,2} + a_5 \delta q_{i+1,2} + a_6 \delta q_{i-1,2} & \\ + a_{10} \delta q_{i,3} + a_{11} \delta q_{i+1,3} + a_{12} \delta q_{i-1,3} & \\ = \frac{\rho_0}{2\Delta t} \left( c_1 \tilde{u}_{i,1}^{t+\Delta t} + c_1 \tilde{u}_{i,1}^{t+\Delta t} + c_2 \tilde{u}_{i,1/2}^{t+\Delta t} + c_2 \tilde{u}_{i-1,1}^{t+\Delta t} + c_2 \tilde{u}_{i-1,1/2}^{t+\Delta t} \right. & \\ \left. + c_3 \tilde{u}_{i,2}^{t+\Delta t} + c_4 \tilde{u}_{i-1,2}^{t+\Delta t} + c_7 \tilde{w}_{i,k+1/2}^{t+\Delta t} \right) \end{aligned} \qquad (C.7) $$
where

$$a_1 = -\frac{\delta z_{1,1}}{\Delta x^2} + \frac{2\beta_{1,1}}{\Delta x} - \frac{\beta_{1,1}^2}{\Delta z_{1,1}} - \frac{\delta z_{1,1}}{\Delta x^2} - \frac{2\beta_{1,1}}{\Delta x} - \frac{\beta_{1,1}^2}{\Delta z_{1,1}} - \frac{1}{\Delta z_{1,3/2}}$$

$$a_2 = \frac{\delta z_{1,1}}{\Delta x^2} - \frac{\beta_{1,1}^2}{\Delta z_{1,1}} \qquad a_3 = \frac{\delta z_{1,1}}{\Delta x^2} - \frac{\beta_{1,1}^2}{\Delta z_{1,1}}$$

$$a_4 = -\frac{\beta_{1,1}}{\Delta x} + \frac{\beta_{1,1}^2}{\Delta z_{1,1}} + \frac{\beta_{1,1}}{\Delta x} + \frac{\beta_{1,1}^2}{\Delta z_{1,1}} + \frac{\alpha_{1,2}}{\Delta x} - \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}} + \frac{\alpha_{1,2}^2}{\Delta z_{1,2}} - \frac{\alpha_{1,2}}{\Delta x} - \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}} + \frac{1}{\Delta z_{1,3/2}}$$

$$a_5 = -\frac{\beta_{1,1}}{\Delta x} + \frac{\beta_{1,1}^2}{\Delta z_{1,1}} - \frac{\alpha_{1,2}}{\Delta x} - \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}} \qquad a_6 = \frac{\beta_{1,1}}{\Delta x} + \frac{\beta_{1,1}^2}{\Delta z_{1,1}} + \frac{\alpha_{1,2}}{\Delta x} - \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}}$$

$$a_7 = a_8 = a_9 = 0 \qquad a_{10} = \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}} + \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}}$$

$$a_{11} = \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}} \qquad a_{12} = \frac{\alpha_{1,2}\beta_{1,2}}{\Delta z_{1,2}}$$

## Appendix D. Non-hydrostatic pressure work

With the notations defined in Appendix C, the work rate of the non-hydrostatic pressure force in a closed domain satisfying $(u_{\frac{1}{2},k}^{t+\Delta t} = u_{\frac{1}{2},k}^{t+\Delta t} = q_{1,k}^{t+\Delta t} = q_{\frac{1}{2},k}^{t+\Delta t} = w_{1,k+\frac{1}{2}}^{t+\Delta t} = w_{\frac{1}{2},k+\frac{1}{2}}^{t+\Delta t} = 0)$ can now be written (at $t + \Delta t$):

At the left hand side, Lines (I)-(IV) correspond to the work rate of the horizontal component of the non-hydrostatic pressure gradient whereas Lines (V)-(VI) correspond to the work rate of its vertical component. The expression of the bottom (Line IV) and surface layer (Lines III and VI) have been written separately. At the right hand side, the non-hydrostatic component of the pressure at $(i, k)$ for $k \in [1, k_{max}]$ times the continuity equation at the same location is recovered at Lines (i), (ii), (iii), (iv) and (vi) whereas the non-hydrostatic pressure component at point $(i, 1)$ times the bottom kinematic condition at the same location is recovered at Lines (v) and (vii) (these terms are underlined). The dissymmetry between the bottom and surface layer is a consequence of the enforcement of the vertical momentum equation in the surface layer ($k = k_{max} + 1/2$) while the vertical velocity at the bottom is computed through the bottom kinematic condition.

## References

Armfield, S., Street, R., 2002. An analysis and comparison of the time accuracy of the fractional-step methods for the Navier-Stokes equations on staggered grids. Int. J. Numer Methods Fluids 38, 255-282.

Auclair, F., Marsaleix, P., Estournel, C., 2000. Sigma coordinate pressure gradient errors: evaluation and reduction by an inverse method. J. Atmos. Ocean Technol. 17, 1348-1367.

Auclair, F., Estournel, C., Marsaleix, P., Pairaud, I., 2006. On coastal embedded modelling. Geophys. Res. Lett. 33 (14), L1 4602. doi:10.1029/2006GL026099.

$$\sum_{i=2}^{i_{max}-1} \sum_{k=1}^{k_{max}} [\Delta z u]_{i,k} \left( \frac{q_{i+1,k} - q_{i,k}}{\Delta x} \right) \qquad (I)$$

$$- \sum_{i=2}^{i_{max}-1} \sum_{k=2}^{k_{max}-1} [\Delta z u]_{i,k} \left[ \alpha_{i,k} \frac{q_{i,k} + q_{i+1,k} - q_{i,k-1} - q_{i+1,k-1}}{\Delta z_{i,k-1}} + \beta_{i,k} \frac{q_{i,k+1} + q_{i+1,k+1} - q_{i,k} - q_{i+1,k}}{\Delta z_{i,k}} \right] \qquad (II)$$

$$- \sum_{i=2}^{i_{max}-1} [\Delta z u]_{i,k_{max}} \left[ \alpha_{i,k_{max}} \frac{q_{i,k_{max}} + q_{i+1,k_{max}} - q_{i,k_{max}-1} - q_{i+1,k_{max}-1}}{\Delta z_{i,k_{max}-1}} + \beta_{i,k_{max}} \frac{-q_{i,k_{max}} - q_{i+1,k_{max}}}{\Delta z_{i,k_{max}}} \right] \qquad (III)$$

$$- \sum_{i=2}^{i_{max}-1} [\Delta z u]_{i,1} \left[ \beta_{1,1} \frac{q_{i,2} + q_{i+1,2} - q_{i,1} - q_{i+1,1}}{\Delta z_{i,1}^{t+\Delta t}} \right] \qquad (IV)$$

$$+ \sum_{i=2}^{i_{max}-1} \sum_{k=1}^{k_{max}-1} [\Delta z w]_{i,k+\frac{1}{2}} \frac{q_{i,k+1} - q_{i,k}}{\Delta z_{i,k+\frac{1}{2}}} \qquad (V)$$

$$+ [\Delta z w]_{i,k_{max}+\frac{1}{2}} \frac{0 - q_{i,k_{max}}}{\Delta z_{i,k_{max}+\frac{1}{2}}} \qquad (VI)$$

$$= - \sum_{i=2}^{i_{max}-1} \sum_{k=1}^{k_{max}} q_{i,k} \left( \frac{[\Delta z u]_{i,k} - [\Delta z u]_{i-1,k}}{\Delta x} \right) \qquad (i)$$

$$- \sum_{i=2}^{i_{max}-1} \sum_{k=2}^{k_{max}-1} q_{i,k} [(\alpha_{i,k} - \beta_{i,k})u_{i,k} + (\alpha_{i-1,k} - \beta_{i-1,k})u_{i-1,k} - \alpha_{i,k+1}u_{i,k+1} - \alpha_{i-1,k+1}u_{i-1,k+1} + \beta_{i,k-1}u_{i,k-1} + \beta_{i-1,k-1}u_{i-1,k-1}] \qquad (ii)$$

$$- \sum_{i=2}^{i_{max}-1} q_{i,k_{max}} [(\alpha_{i,k_{max}} - \beta_{i,k_{max}})u_{i,k_{max}} + (\alpha_{i-1,k_{max}} - \beta_{i-1,k_{max}})u_{i-1,k_{max}} + \beta_{i,k_{max}-1}u_{i,k_{max}-1} + \beta_{i-1,k_{max}-1}u_{i-1,k_{max}-1}] \qquad (iii)$$

$$- \sum_{i=2}^{i_{max}-1} q_{i,1} [\alpha_{i,\frac{1}{2}} u_{i,\frac{1}{2}} + \alpha_{i-1,\frac{1}{2}} u_{i-1,\frac{1}{2}} - \beta_{i,1} u_{i,1} - \beta_{i-1,1} u_{i-1,1} - \alpha_{i,2} u_{i,2} - \alpha_{i-1,2} u_{i-1,2}] \qquad (iv)$$

$$+ q_{i,1} (\alpha_{i,\frac{1}{2}} u_{i,\frac{1}{2}} + \alpha_{i-1,\frac{1}{2}} u_{i-1,\frac{1}{2}}) \qquad (v)$$

$$- \sum_{i=2}^{i_{max}-1} \sum_{k=1}^{k_{max}} q_{i,k} (w_{i,k+\frac{1}{2}} - w_{i,k-\frac{1}{2}}) \qquad (vi)$$

$$- q_{i,1} w_{i,\frac{1}{2}} (\text{vii}) \qquad (D.1)$$
Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal ocean circulation model. In: Heaps, Norman S. (Ed.), Three-dimensional Coastal Ocean Modell. Coastal and Estuarine Sciences, vol. 4. AGU, Washington D.C., pp. 1–16.

Brickman, D., Loder, J.W., 1993. The energetics of the internal tide on northern Georges Bank. J. Phys. Oceanogr. 23, 409–424.

Casulli, V., 1999. A semi-implicit finite difference method for non-hydrostatic, free-surface flows. Int. J. Numer. Methods Fluids 30, 425–440.

Chen, X.J., 2003. A fully hydrodynamic model for three-dimensional, free-surface flows. Int. J. Numer. Methods Fluids 42, 929–952.

Chorin, A., 1967. A numerical method for solving incompressible viscous flow problems. J. Comput. Phys. 2, 12–26.

Dukowicz, J.K., Dvinsky, A.S., 1992. Approximate factorization as a high order splitting for the implicit incompressible flow equations. J. Comput. Phys. 102 (2), 336–347.

Ferziger, J.H., 1998. Numerical Methods for Engineering Applications. John Wiley & Sons. pp. 400.

Floor, 2009. Energetics of Internal Tide Generation, Propagation and Dissipation. Ph.D. Thesis, Univ. of Toulouse Paul Sabatier. <http://pocobs-mip.fr/DOCUMENTS/FLOOR_jochem_w_200910_phd.pdf>.

Floor J.W., Auclair, F., Marsaleix, P., under review. Energy transfers in internal tide generation, propagation and dissipation in the deep ocean.

Gaspar, P., Gregoris, Y., Lefevre, J.M., 1990. A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: tests at station Papa and long-term upper ocean study site. J. Geophys. Res. 95, 16,179–16,193.

Gresho, P.M., Sani, R.L., 1987. On pressure boundary conditions for the incompressible NavierStokes equations. Int. J. Numer. Methods Fluids 7, 1111.

Hodges, B.R., Laval, B., Wadzuk, B.M., 2006. Numerical error assessment and a temporal horizon for internal waves in a hydrostatic model. Ocean Modell. 13 (1), 44–64. doi:10.1016/j.ocemod.2005.09.005.

Horn, D.A., Imberger, J., Ivey, G.N., 2001. The degeneration of large-scale interfacial gravity waves in lakes. J. Fluid Mech. 434, 181–207.

Kanarska, Y., Shchepetkin, A., McWilliams, J.C., 2007. Algorithm for non-hydrostatic dynamics in the Regional Oceanic Modelling System. Ocean Modell. 18, 143–174.

Lamb, K.G., 1994. (L94): Numerical experiments of internal wave generation by strong tidal flow across a finite amplitude bank edge. J. Geophys. Res. 99 (C1), 843–864.

Lamb, K.G., 2007. Energy and pseudoenergy flux in the internal wave field generated by tidal flow over topography. Cont. Shelf Res. 27, 1208–1232.

La Violette, P.E., Johnson, D.R., Brooks, D.A., 1990. Sun-glitter photographs of Georges Bank and the Gulf of Maine from the space shuttle. Oceanography 3, 43–49.

Loder, J.W., 1980. Topographic rectification of tidal currents on the sides of Georges Bank. J. Phys. Oceanogr. 10, 1399–1416.

Loder, J.W., Brickman, D., Horne, E.P.W., 1992. Detailed structure of currents and hydrography on the northern side of Georges Banks. J. Geophys. Res. 97, 14,331–14,351.

Mahadevan, A., Oliger, J., Street, R., 1996. A non-hydrostatic mesoscale ocean model. Part I: Well-posedness and scaling. J. Phys. Oceanogr. 26, 1868–1880.

Marsaleix, P., Auclair, F., Estournel, C., 2006. Considerations on open boundary conditions for regional and coastal ocean models. J. Atmos. Oceanic. Technol. 23, 1604–1612.

Marsaleix, P., Auclair, F., Floor, J.W., Herrmann, M.J., Estournel, C., Pairaud, I., Ulses, C., 2008. (M08): Energy conservation issues in sigma-coordinate free-surface ocean models. Ocean Modell. 20 (1), 61–89.

Marsaleix, P., Auclair, F., Floor, J.W., Herrmann, M.J., Estournel, C., 2009. Low-order pressure gradient schemes in sigma coordinate models: the seamount test revisited. Ocean Modell. 30 (2–3), 169–177. doi:10.1016/j.ocemod.2009.06.011.

Marshall, J., Schott, F., 1999. Open-ocean convection: observations, theory, and models. Rev. Geophys. 37, 1–64.

Marshall, J., Hill, C., Perelman, L., Adcroft, A., 1997a. Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. J. Geophys. Res. 102 (c3), 5733–5753.

Marshall, J., Adcroft, A., Hill, C., Perelman, L., Heisey, C., 1997b. A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers. J. Geophys. Res. 102 (c3), 5753–5766.

Mellor, G.L., Ezer, T., Oey, L.Y., 1994. The pressure gradient conundrum of sigma coordinate ocean models. J. Atmos. Oceanic. Technol. 11 (4), 1126–1134. Part 2.

New, A.L., Da Silva, J.C.B., 2002. Remote-sensing evidence for the local generation of internal soliton packets in the central Bay of Biscay. Deep-Sea Res. I 49 (2002), 915–934.

Saad, Y., 2003. Iterative Methods for Sparse Linear Systems. SIAMS Editions. ISBN 0-89871-534-2.

Schlesinger, R.E., Uccellini, L.W., Johnson, D.R., 1983. The effects of the Asselin filter on numerical solutions to the linearized shallow-water wave equations. Mon. Weather Rev. 111, 455–467.

Shchepetkin, A.F., McWilliams, J.C., 2005. The Regional Ocean Modeling System (ROMS): a split-explicit, free-surface, topography-following coordinates ocean model. Ocean Modell. 9, 347–404.

Staquet, C., Sommeria, C.J., 2002. Mixing by breaking internal gravity waves: from instabilities to turbulence. Ann. Rev. Fluid Mech. 34, 559–593.

Stelling, G., Zijlema, M., 2003. An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation. Int. J. Numer. Methods Fluids 43, 1–23. doi:10.1002/d.595.

Zijlema, M., Stelling, G.S., 2005. Further experiences with computing non-hydrostatic free-surface flows involving water waves. Int. J. Numer. Methods Fluids 48, 169–197.