# Energy conservation issues in sigma-coordinate free-surface ocean models

Patrick Marsaleix, Francis Auclair, Jochem Willem Floor, Marine Julie Herrmann, Claude Estournel, Ivane Pairaud, Caroline Ulses

► To cite this version:

Patrick Marsaleix, Francis Auclair, Jochem Willem Floor, Marine Julie Herrmann, Claude Estournel, et al.. Energy conservation issues in sigma-coordinate free-surface ocean models. Ocean Modelling, 2008, 20 (1), pp.61-89. 10.1016/j.ocemod.2007.07.005. hal-02110201

HAL Id: hal-02110201

https://hal.science/hal-02110201

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
Energy conservation issues in sigma-coordinate
free-surface ocean models

Patrick Marsaleix *, Francis Auclair, Jochem Willem Floor, Marine Julie Herrmann, Claude Estournel, Ivane Pairaud, Caroline Ulses

Laboratoire d'Aérologie – CNRS et Université Paul Sabatier, 14, Avenue Edouard Belin, 31400 Toulouse, France

This paper focuses on the energy conservation properties of a hydrostatic, Boussinesq, coastal ocean model using a classic finite difference method. It is shown that the leapfrog time-stepping scheme, combined with the sigma-coordinate formalism and the motions of the free surface, prevents the momentum advection from exactly conserving energy. Because of the leapfrog scheme, the discrete form of the kinetic energy depends on the product of velocities at odd and even time steps and thus appears to be possibly negative when high-frequency modes develop. Besides, the study of the energy balance clarifies the numerical choices made for the computation of mixing processes. The time-splitting technique used to reduce the computation costs associated to the resolution of surface waves leads to the well-known external and internal mode equations. We show that these equations do not conserve energy if the coupling of these two modes is forward in time. Even if non-linear terms are negligible, this shortcoming can be significant regarding the pressure gradient term 'fro-zen' over a baroclinic time step. An alternative energy-conserving time-splitting technique is proposed in this paper. Discussion and conclusions are conducted in the light of a set of numerical experiments dedicated to surface and internal gravity waves.

Keywords: Energy conserving; Free surface; Sigma coordinate; Gravity waves

# 1. Introduction

The purpose of this study is to evaluate and enforce the energy conservation properties of the sigma-coor-
dinate free-surface ocean model Symphonie, described in Auclair et al. (2000). The model is inspired by the
POM model (Blumberg and Mellor, 1987; hereafter: BM87), widely used by the coastal ocean modelling

<footnote>* Corresponding author.

*E-mail addresses:* marp@aeroobs-mip.fr (P. Marsaleix), aucf@aeroobs-mip.fr (F. Auclair), floj@aeroobs-mip.fr (J.W. Floor), herm@aeroobs-mip.fr (M.J. Herrmann), estc@aeroobs-mip.fr (C. Estournel), ivane.pairaud@servmail.hmg.inpg.fr (I. Pairaud), c.ulses@nioo.knaw.nl (C. Ulses).</footnote>
community. As we will see in the following sections, some of the algorithmic choices (notably the pressure gradient force) are however significantly distinct from POM.

The incompressible Navier–Stokes equations are solved on a staggered C-grid using a classic finite difference method. In this context, the global conservation of primary quantities such as temperature, salinity or momentum is easily obtained as long as the advection problem is put in the form of a flux-divergence, by combining the continuity equation with tracer and momentum equations. The sigma-coordinate formalism is well adapted to the conservation issue posed by the free surface, since it leads to a set of equations in which the volume variation of grid cells is naturally included. The problem is slightly complicated by time-splitting techniques, used to limit calculus costs regarding surface gravity waves. The coupling of the external and internal modes (BM87) is partly conditioned by conservation issues, since the barotropic current involved in 3D advection fluxes is deduced from a time average of the external mode solution, consistent with volume variation of grid cells over one internal time step. As the integral of the flux-divergence over the domain volume is equivalent to the surface integral of the fluxes over the boundaries, global conservation of primary variables directly concerns the boundary conditions of the model. The global tracer balance related to advection and diffusion should remain intact in case of solid boundaries, since fluxes through the coastal or the bottom boundaries should equal zero. As far as velocities are concerned, bottom friction (and eventually coastline friction in case of no-slip condition) modifies the global momentum balance. Tracer and momentum balances strongly depend on open boundary conditions. Most open boundary schemes are based on local dynamics considerations only, such as the radiation of outgoing waves or the definition of incoming variables. Nonetheless, boundary conditions can also satisfy global constraints (Shulman et al., 1998; Marsaleix et al., 2006).

In a primitive equation model, the conservation of derived integrals of motions (e.g., net heat content, tracer quadratic variance, kinetic energy, etc.) is generally not straightforward. Since these quantities cannot be independently enforced, they actually are the consequence of the discretised momentum and tracer equations. The corresponding conservation properties are thus far from obvious and may not all be satisfied in some numerical methods. As a consequence, modellers possibly have to choose between for instance energy and enstrophy conservation (Sadourny, 1975; Madec et al., 1991). Nevertheless, powerful theories on oceanic circulations are based on vorticity or energy conservation properties. This largely explains the interest of model developers in energy-, vorticity-, and tracer-conserving methods (Arakawa and Lamb, 1977).

Although several studies have reported that energetically non-consistent models may perform better than energetically consistent models (Shchepetkin and McWilliams, 2003; Beckmann and Haidvogel, 1993), conservation of kinetic energy (hereafter KE) can also be viewed as a guarantee of numerical stability. As pointed out by Ferziger and Peric (2002, p. 161): “*if a numerical method is energy conservative and the net energy flux through the surface is zero the total kinetic energy in the domain does not grow with time*”, so that “*the velocity at every grid point in the domain must remain bounded*”. An energy conservative model must fulfil four main requirements. First, for closed domains, the momentum advection must not change the total KE. Second, the Coriolis term must not change the KE. Third, momentum-mixing terms must constitute an energy sink. Fourth, the variation of KE induced by the pressure gradient must be consistent with the variation of potential energy (hereafter PE) induced by the vertical advection of density. The momentum advection scheme proposed by Lilly (1965) and BM87 for the C-grid can be regarded as a particular case of the Arakawa and Lamb (1977, Eqs. (90)–(93)) energy-conserving scheme (note that the more general form of the Arakawa and Lamb’s scheme also enables enstrophy conservation for two-dimensional non-divergent flow). Concerning the Coriolis term, the discrete form used in BM87 is the same as in Arakawa and Lamb (1977). It is relatively easy to guarantee KE conservation for this simple term (see Eqs. (101)–(104) in Arakawa and Lamb, 1977), so that it is somehow of lesser interest than the other terms of the momentum equation. In the following, we neglect this process and consider a set of equations in the *Oxz* vertical plane. Concerning mixing terms, the expected consistency with the turbulent closure scheme is, to our knowledge, rather rarely mentioned in the literature (see, however, Burchard, 2002). This particular point is addressed in this paper. Finally, the pressure gradient force (PGF) is probably the most important term because of the large exchange between kinetic and potential energies through vertical advection of density. The importance of the consistency between the discrete forms of the PGF and the density advection is long-known in atmospheric modelling (Janjic, 1977; Arakawa and Suarez, 1983). In oceanic models using terrain following coordinates, PGF schemes are sensitive to truncation errors because of potentially steep bathymetric slopes, particularly in coastal areas (Haney, 1991). This issue
motivated several studies concerning ways to improve the accuracy of numerical schemes. Gerdes (1993), Song (1998), Song and Wright (1998), suggested replacing the straightforward PGF schemes (assuming density as a piecewise constant distribution within each grid box) by second-order density Jacobian type schemes, actually more accurate. Shchepetkin and McWilliams (2003, Eqs. (A14) and (A19)) however point out a minor lack of energetic consistency of the barotropic part of the PGF with the definition of discrete PE due to free-surface motion. To our knowledge, the effect of time differencing on conservation of quadratic quantities has been rarely discussed in the literature (see, however, the discussion about Crank–Nicholson scheme in Ferziger and Peric, 2002, p. 163). Much of the explanation may lie in the fact that time steps are usually much smaller than timescales associated with processes accessible to hydrostatic Boussinesq models (energy conservation of the time-stepping scheme could be reasonably expected because of its small truncation errors), whereas in regional or coastal models, the grid mesh is usually of the same order as the first baroclinic deformation radii. Similarly, few studies deal with energy conservation properties of time-splitting methods enabling to compute separately the baroclinic and barotropic parts of the flow (Higdon, 2005). In consequence, the present paper focuses on aspects rarely studied in the past: lack of energy conservation due to leapfrog temporal schemes, including free-surface motion, Asselin filter and time splitting of internal and external modes.

The plan of this paper is as follows: first, we describe the model. Secondly, we consider the details of the energy balance, starting with the time variation and velocity advection terms. Next, we discuss the mixing terms, including the Asselin filter (Asselin, 1972), considered here as a temporal mixing term and then we examine the pressure gradient terms. We subsequently address the PE balance, looking at the respective roles of advection and mixing tracer terms. Finally, we present a numerical experiment of internal wave generation over a submarine ridge, after which we present our conclusions.

## 2. Model description

The model grid is adjusted to the bathymetry thanks to a sigma-coordinate system leading to the following model equations:

$$ \frac{\partial \tilde{u}}{\partial t^*} + \frac{\partial \tilde{u}u}{\partial x^*} + \frac{\partial \tilde{u}v}{\partial y^*} + \frac{\partial \omega u}{\partial \sigma} - f \tilde{v} = -\frac{H}{\rho_0} \frac{\partial p}{\partial x} + \frac{\partial}{\partial x^*} \tilde{K}^x \frac{\partial u}{\partial x^*} + \frac{\partial}{\partial y^*} \tilde{K}^y \frac{\partial u}{\partial y^*} + \frac{\partial}{\partial \sigma} \frac{K^z}{H} \frac{\partial u}{\partial \sigma} \quad (1) $$

$$ \frac{\partial \tilde{v}}{\partial t^*} + \frac{\partial \tilde{v}u}{\partial x^*} + \frac{\partial \tilde{v}v}{\partial y^*} + \frac{\partial \omega v}{\partial \sigma} + f \tilde{u} = -\frac{H}{\rho_0} \frac{\partial p}{\partial y} + \frac{\partial}{\partial x^*} \tilde{K}^x \frac{\partial v}{\partial x^*} + \frac{\partial}{\partial y^*} \tilde{K}^y \frac{\partial v}{\partial y^*} + \frac{\partial}{\partial \sigma} \frac{K^z}{H} \frac{\partial v}{\partial \sigma} \quad (2) $$

where asterisks refer to the sigma-coordinate system and tildes indicate a multiplication by the scale factor $H = h + \eta$, which is the sum of the water column thickness at rest and the sea level anomaly, hereafter SLA (e.g., $\tilde{u} = Hu$). The appearance of the asterisk for time in partial differentiation with respect to time is explained by the fact that sigma levels actually move (due to SLA variations) with respect to the absolute system of coordinates. Details on the coordinate transformation leading to the formulation of derivatives in the sigma-coordinate system can be found in Johns et al. (1983). The division of the pressure gradient by a reference density comes from the Boussinesq approximation. Hydrostatic equilibrium and the expression of the horizontal derivative in the sigma-coordinate system lead to

$$ -\frac{1}{\rho_0} \frac{\partial p}{\partial x} = -g \frac{\partial \eta}{\partial x} - \frac{1}{\rho_0} \left( \frac{\partial p'}{\partial x^*} + \frac{\partial z}{\partial x^*} g \rho' \right) \quad (3) $$

where $p' = g \int_{z}^{\eta} \rho' dz'$ is the hydrostatic pressure anomaly associated with the density anomaly, $\rho' = \rho - \rho_0$. This formulation of the pressure gradient is characterised by the appearance of a term depending on the SLA only (first term on the right hand side of expression (3)). This is motivated by considerations on the free-surface computation that will be detailed later on. For sake of clarity of the arguments developed in this paper, density is related to temperature and salinity through a linear equation of state. Thus, the density variations are governed by a simple advection diffusion tracer equation, linearly deduced from those for temperature and salinity, namely

$$ \frac{\partial \tilde{p}}{\partial t^*} + \frac{\partial \tilde{u} \rho}{\partial x^*} + \frac{\partial \tilde{v} \rho}{\partial y^*} + \frac{\partial \omega \rho}{\partial \sigma} = \frac{\partial}{\partial x^*} \tilde{K}_T^x \frac{\partial \rho}{\partial x^*} + \frac{\partial}{\partial y^*} \tilde{K}_T^y \frac{\partial \rho}{\partial y^*} + \frac{\partial}{\partial \sigma} \frac{K_T^z}{H} \frac{\partial \rho}{\partial \sigma} \quad (4) $$
Using a linear equation of state is reasonable in a coastal context, since pressure effects on water density remain limited as long as we consider the shallow water approximation. We also note that a non-linear equation possibly destroys the energetic consistency of the pressure gradient force (Gerdes, 1993), but the limitation of this assumption should be kept in mind. The parameterisation of sub-grid processes is given by diffusion terms on the right hand side (hereafter: RHS) of (4), involving mixing coefficients that can be calculated in different ways. Although $(K^x, K^z)$ are generally calculated by a turbulent closure scheme representing the energy cascade toward small scales (Mellor and Yamada, 1982), the other mixing coefficients $(\tilde{K}^x, \tilde{K}^y, \tilde{K}_T^x, \tilde{K}_T^y)$ are usually given by more basic techniques intended to reduce numerical noise (Smagorinsky, 1993), if not simply reduced to a fixed value. Momentum and tracer advection have the form of a flux-divergence, resulting from the combination of the primitive equations with the continuity equation. The latter is given by

$$ \frac{\partial \eta}{\partial t^*} + \frac{\partial \tilde{u}}{\partial x^*} + \frac{\partial \tilde{v}}{\partial y^*} + \frac{\partial \omega}{\partial \sigma} = 0 \qquad (5) $$

This equation is used to compute $\omega$, the vertical velocity in sigma coordinates (BM87). We will just retain here that $\omega$ vanishes at the bottom and at the sea surface, so that the depth integral of Eq. (5) leads to a simple relationship linking the SLA variations to the depth-averaged component of the current, $(\hat{u}, \hat{v}) = H^{-1} \int_{-h}^{h} (u, v) dz$. This relationship reads

$$ \frac{\partial \eta}{\partial t^*} + \frac{\partial \tilde{u}}{\partial x^*} + \frac{\partial \tilde{v}}{\partial y^*} = 0 \qquad (6) $$

with $(\tilde{u}, \tilde{v}) = H(\hat{u}, \hat{v})$. In the case of a finite difference method using a leapfrog scheme for the time-stepping procedure, Eq. (6) is computed first, i.e. the SLA variation term is computed from the divergence of the depth-averaged current known from the previous iteration. Then Eq. (5) gives $\omega$, using the horizontal current components from the previous iteration and the SLA variation term provided by Eq. (6). Momentum Eqs. (1) and (2) and density equation (4) are then computed using mixing coefficient values from the previous iteration. In the case of an explicit method, the barotropic part of the pressure gradient (first term at the RHS of Eq. (3)) poses a problem of numerical stability, generally leading to a separate treatment of the barotropic variable, using a time step specially suited for surface wave computation. This time step is normally much smaller than what we would obtain if baroclinic physics was only considered. The so-called time-splitting technique, presented in BM87, is based on Eq. (6), completed with two equations for the depth-averaged current, deduced from the vertical integral of momentum equations (1) and (2), i.e.,

$$ \frac{\partial \tilde{u}}{\partial t^*} + \frac{\partial \tilde{u} \tilde{u}}{\partial x^*} + \frac{\partial \tilde{v} \tilde{u}}{\partial y^*} - f \tilde{v} = -gH \frac{\partial \eta}{\partial x} + \frac{\partial}{\partial x^*} \tilde{K}^x \frac{\partial \tilde{u}}{\partial x^*} + \frac{\partial}{\partial y^*} \tilde{K}^y \frac{\partial \tilde{u}}{\partial y^*} + \left[ \frac{K^z}{H} \frac{\partial u}{\partial \sigma} \right]_{\sigma_0}^{\sigma_1} \\ - H \int_{\sigma_0}^{\sigma_1} \left( \frac{1}{\rho_0} \frac{\partial p'}{\partial x} + \frac{\partial \tilde{u}' u'}{\partial x^*} + \frac{\partial \tilde{v}' u'}{\partial y^*} \right) d\sigma \qquad (7) $$

$$ \frac{\partial \tilde{v}}{\partial t^*} + \frac{\partial \tilde{u} \tilde{v}}{\partial x^*} + \frac{\partial \tilde{v} \tilde{v}}{\partial y^*} + f \tilde{u} = -gH \frac{\partial \eta}{\partial y} + \frac{\partial}{\partial x^*} \tilde{K}^x \frac{\partial \tilde{v}}{\partial x^*} + \frac{\partial}{\partial y^*} \tilde{K}^y \frac{\partial \tilde{v}}{\partial y^*} + \left[ \frac{K^z}{H} \frac{\partial v}{\partial \sigma} \right]_{\sigma_0}^{\sigma_1} \\ - H \int_{\sigma_0}^{\sigma_1} \left( \frac{1}{\rho_0} \frac{\partial p'}{\partial y} + \frac{\partial \tilde{u}' v'}{\partial x^*} + \frac{\partial \tilde{v}' v'}{\partial y^*} \right) d\sigma \qquad (8) $$

where $(\sigma_0, \sigma_1)$ are the bottom and surface values of $\sigma$. Our model uses the sigma-coordinate transformation of Johns and Oguz (1987) with $(\sigma_0, \sigma_1) = (0, 1)$. Let us note however that this is just a special case of the more general terrain following coordinate transformation. Besides the numerical aspects developed in this paper can virtually be transposed to a wide range of coordinate transformations, like for instance the hybrid terrain pressure coordinate system often used in atmospheric models (see for instance Arakawa and Suarez, 1983). Note that $(\tilde{K}^x, \tilde{K}^y)$ correspond to the depth-averaged values of $(K^x, K^y)$. Eqs. (6)–(8) (hereafter external mode equations) are used to compute $(\eta, \hat{u}, \hat{v})$ with a small time step. The terms under the integral at the RHS of Eqs. (7) and (8) are frozen over a time period equivalent to one time step of the so-called internal mode equations. The latter consist of momentum equations (1) and (2), and tracer equation (4). The time-splitting technique works in the following manner: at the beginning of each iteration, a baroclinic current, $(u', v')$, is
obtained by removing the depth-averaged value from the absolute current provided by the internal mode solution. This baroclinic current and the three-dimensional pressure anomaly, $p'$, provide the last term at the RHS of Eqs. (7) and (8). The latter and Eq. (6) are then integrated over one external mode sequence with a small time step. At the end of the external mode sequence, the depth-averaged component of the absolute current is updated using a time-averaged value of the external mode solution. Three-dimensional equations, e.g. continuity equation (5), momentum equations (1) and (2) and tracer equation are then computed. We note that Eqs. (1) and (2) can actually be viewed as a pair of equations for the baroclinic components of the current ($u'$, $v'$) = ($u - \hat{u}$, $v - \hat{v}$).

### 3. Kinetic energy balance

#### 3.1. Time variation and advection terms

##### 3.1.1. Analytical formulation

An equation for the KE rate-of-change is obtained by multiplying the momentum equations by their respective velocity components. A global balance is obtained by integrating the latter over the numerical domain and over the duration of the simulation. The time splitting of the model into two modes, the so-called external and internal modes, naturally leads to summing the barotropic and baroclinic balances (since the orthogonality of the modes leads to $\int_{\sigma_0}^{\sigma_1} u \frac{\partial u}{\partial t} d\sigma = \hat{u} \frac{\partial \hat{u}}{\partial t} + \int_{\sigma_0}^{\sigma_1} u' \frac{\partial u'}{\partial t'} d\sigma$), separately obtained on one hand by multiplying the external mode Eqs. (7) and (8) by the depth-averaged components of current, ($\hat{u}$, $\hat{v}$), and on the other hand by multiplying the internal mode equations (1) and (2) by the baroclinic components of current, ($u'$, $v'$). Let us recall that, although it appears to be a set of equations for the total current, the internal mode is only used to update the baroclinic component of velocities (BM87). We also note that, as we have $\int_{\sigma_0}^{\sigma_1} u' \frac{\partial \hat{u}}{\partial t'} d\sigma = \int_{\sigma_0}^{\sigma_1} u' \frac{\partial \hat{u}'}{\partial t'} d\sigma$, it is in fact no different to compute the baroclinic KE balance from Eqs. (1) and (2) or from a pair of equations that would resolve $u'$ and $v'$ specifically, a property largely used in the following. Moreover, in order to simplify the presentation of the energy balance without reducing the scope of our conclusions, we now consider the previous equations in a vertical $Oxz$-plane. In other words, we neglect current components and derivatives in the $Oy$ direction. After some simple algebraic manipulations, the contribution of the time variation terms to the global energy balance is written:

$$ \int_{t0}^{t1} \int_{x0}^{x1} \left( \hat{u} \frac{\partial \hat{u}}{\partial t^*} + \int_{\sigma_0}^{\sigma_1} u' \frac{\partial \hat{u}'}{\partial t^*} d\sigma \right) dx^* dt^* = \int_{x0}^{x1} \left( \Delta \tilde{E}_C + \int_{\sigma_0}^{\sigma_1} \Delta \tilde{E}'_C d\sigma \right) dx^* \\ + \int_{t0}^{t1} \int_{x0}^{x1} \left( \hat{E}_C + \int_{\sigma_0}^{\sigma_1} E'_C d\sigma \right) \frac{\partial \eta}{\partial t^*} dx^* dt^* \quad (9) $$

where $\hat{E}_C = \hat{u}^2/2$ and $E'_C = u'^2/2$ respectively represent the KE of the external and internal mode and $\Delta$ stands for a difference between the final and initial times, $t1$ and $t0$. The last term at the RHS of Eq. (9), depending on SLA variations, should vanish when the contribution of advection terms is taken into account. The latter is given by

$$ \int_{t0}^{t1} \int_{x0}^{x1} \left[ \hat{u} \left( \frac{\partial \tilde{u}\hat{u}}{\partial x^*} + \int_{\sigma_0}^{\sigma_1} \frac{\partial \tilde{u}'u'}{\partial x^*} d\sigma \right) + \int_{\sigma_0}^{\sigma_1} u' \left( \frac{\partial (\tilde{u}u' + \tilde{u}'\hat{u} + \hat{u}'u')}{\partial x^*} + \frac{\partial \omega(\hat{u} + u')}{\partial \sigma} \right) d\sigma \right] dx^* dt^* \quad (10) $$

##### 3.1.2. Discretised formulation

The discrete formulation makes use of an Arakawa C-grid (Arakawa and Lamb, 1977). Figs. 1 and 2 show the location of the different variables on the model grid. Momentum equations are computed at horizontal grid half-integer indexes, $i$, ranging from $1 + 1/2$ to $M - 1/2$ and vertical integer indexes, $k$, increasing upward from 1 to $R - 1$. Lateral boundary conditions are applied on $i = 1/2$ and $M + 1/2$. Tracer equations are computed at horizontal grid integer indexes, $i$, ranging from 1 to $M$ and vertical integer indexes, $k$. We use the following sigma-coordinate transformation:

$$ \sigma = (z + h)/H \quad (11) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of the C grid in the Ox direction. Triangles indicate location of horizontal velocities. Circles indicate location of tracers, surface pressure, bathymetry, vertical velocities and TKE variables.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Sketch of the C grid in the Oz direction. Rhombuses indicate the location of horizontal velocities and tracers; squares indicate location of vertical velocities and TKE variables.</figcaption></figure>

Half-integer index sigmas, ranging from 0 at the bottom to 1 at the free surface, are defined first. Note that a non-linear distribution can be employed at this stage. Integer index sigmas are then deduced from the latter according to

$$ \sigma_k = (\sigma_{k+1/2} + \sigma_{k-1/2})/2 \quad (12) $$

from which it follows that tracer nodal value at integer indexes $(i, k)$ represents the mean over the cell surface, $\Delta x \Delta z_{i,k}$, where $\Delta x$ is the distance between horizontal grid-points and the thickness of the cell is defined as

$$ \Delta z_{i,k} = H_i \Delta \sigma_k = H_i (\sigma_{k+1/2} - \sigma_{k-1/2}) \quad (13) $$

It follows from (11)-(13) that

$$ \sigma_k = \sum_{q=1,k-1} \Delta\sigma_q + 0.5\Delta\sigma_k \quad (14) $$

and

$$ z_{i,k} - z_{i,k-1} = H_i (\Delta\sigma_k + \Delta\sigma_{k-1})/2 \quad (15) $$

We recall that $(\tilde{u}, \tilde{u}') = H(\hat{u}, u')$ and that $H^{t+\Delta t} - H^t = \eta^{t+\Delta t} - \eta^t$. To derive discrete schemes we will use the following averaging operators:

$$ \begin{aligned} \tilde{a}_{i,k}^x &= \frac{1}{2}(a_{i+1/2,k} + a_{i-1/2,k}) \\ \tilde{a}_{i+1/2,k}^x &= \frac{1}{2}(a_{i,k} + a_{i+1,k}) \\ \tilde{a}_{i,k}^z &= \frac{1}{2}(a_{i,k+1/2} + a_{i,k-1/2}) \\ \tilde{a}_{i,k+1/2}^z &= \frac{1}{2}(a_{i,k} + a_{i,k+1}) \end{aligned} \quad (16) $$
Eq. (9) represents the sum of barotropic and baroclinic relations, namely $\int_{t0}^{t1} \int_{x0}^{x1} \hat{u} \frac{\hat{\bar{u}}}{\hat{c}^t} dx^* dt^* = \int_{x0}^{x1} \Delta \tilde{E}_C dx^* + \int_{t0}^{t1} \int_{x0}^{x1} \hat{E}_C \frac{\hat{\bar{\eta}}}{\hat{c}^t} dx^* dt^*$ and $\int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} u' \frac{\hat{\bar{u}}'}{\hat{c}^t} d\sigma dx^* dt^* = \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} \Delta \tilde{E}'_C d\sigma dx^* + \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} E'_C d\sigma \frac{\hat{\bar{\eta}}}{\hat{c}^t} dx^* dt^*$ respectively. Their finite difference counterpart, based on our common leapfrog scheme, is then given by

$$ \sum_{\substack{i=1,M-1 \\ t=t0,t1}} \left( \hat{u}' \frac{\tilde{u}^{t+\Delta t_e} - \tilde{u}^{t-\Delta t_e}}{2\Delta t_e} \Delta x \Delta t_e \right)_{i+1/2} = \sum_{i=1,M-1} \left( \frac{\Delta x}{2} \left( \hat{u}^{t1} \tilde{u}^{t1+\Delta t_e} - \hat{u}^{t0} \tilde{u}^{t0-\Delta t_e} \right) \right)_{i+1/2} \\ + \sum_{\substack{i=1,M-1 \\ t=t0,t1-\Delta t_e}} \left( \frac{\hat{u}' \hat{u}^{t+\Delta t_e}}{2} \frac{\bar{\eta}^{x+\Delta t_e} - \bar{\eta}^x}{\Delta t_e} \Delta x \Delta t_e \right)_{i+1/2} \quad (17) $$

and

$$ \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t0,t1}} \left( u'' \frac{\tilde{u}''^{t+\Delta t} - \tilde{u}''^{t-\Delta t}}{2\Delta t} \Delta x \Delta \sigma \Delta t \right)_{i+1/2,k} = \sum_{\substack{i=1,M-1 \\ k=1,R-1}} \left( \frac{\Delta x \Delta \sigma}{2} \left( u''^1 \tilde{u}''^{1+\Delta t} - u''^0 \tilde{u}''^{0-\Delta t} \right) \right)_{i+1/2,k} \\ + \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t0,t1-\Delta t}} \left( \frac{u'' u''^{t+\Delta t}}{2} \frac{\bar{\eta}^{x+\Delta t} - \bar{\eta}^x}{\Delta t} \Delta x \Delta \sigma \Delta t \right)_{i+1/2,k} \quad (18) $$

Eqs. (17) and (18) respectively represent the external and internal contributions. In Eq. (17), $t$ varies from $t0$ to $t1$ in increments of $\Delta t_e$, the external time step, and in increments of $\Delta t$, the internal time step, in Eq. (18). The first term on the RHS of Eqs. (17) and (18) respectively, corresponds to the difference of barotropic (respectively baroclinic) KE between the final and initial times ($t_1, t_0$). Let us note that, according to Eqs. (17) and (18), the discrete form of the kinetic barotropic and baroclinic energy, $(\hat{u}'\tilde{u}^{t+\Delta t_e}, u''\tilde{u}''^{t+\Delta t})/2$, can be negative, since velocities at iterations $t$ and $t + \Delta t_e$ or $t + \Delta t$ can have opposite signs.

### 3.1.3. Surface motion

The last term of the RHS of Eq. (17) should equilibrate with the energy balance related to advection terms, since, in case of closed boundaries, we actually have: $\int_{x0}^{x1} \hat{u} \frac{\hat{\bar{u}}}{\hat{c}^x} dx^* = - \int_{x0}^{x1} \frac{\hat{u}^2}{2} \frac{\hat{\bar{\eta}}}{\hat{c}^x} dx^*$. Using the numerical counterpart of the barotropic continuity equation (6),

$$ \left( \tilde{u}'_{i+1/2} - \tilde{u}'_{i-1/2} \right) / \Delta x = -(\eta_i^{t+\Delta t_e} - \eta_i^{t-\Delta t_e}) / 2\Delta t_e, \quad (19) $$

the numerical rule of integration by parts given by (A.1) and (A.2) (in Appendix) and assuming closed boundary ($\hat{u}_{1/2} = \hat{u}_{M+1/2} = 0$), one can show that the finite difference counterpart, based on our common centred scheme ($F_i = \bar{\tilde{u}}_i \bar{\tilde{u}}_i^x$), verifies

$$ \sum_{\substack{i=1,M-1 \\ t=t0,t1}} \hat{u}'_{i+1/2} (F_{i+1} - F_i) \Delta t_e = - \sum_{\substack{i=1,M-1 \\ t=t0,t1}} \frac{\hat{u}^2}{2} \frac{\bar{\eta}^{x+\Delta t_e} - \bar{\eta}^{x-\Delta t_e}}{2\Delta t_e} \Delta x \Delta t_e)_{i+1/2} \quad (20) $$

Note that the expected equilibrium between (20) and the second term of the RHS of (17) is not exactly achieved, since on the one hand $\hat{u}'\hat{u}^{t+\Delta t_e}/2$ is not exactly equivalent to $\hat{u}^2/2$, while on the other hand the forward derivative, $(\bar{\eta}^{x+\Delta t_e} - \bar{\eta}^x)/\Delta t_e$, is not strictly equivalent to the centred derivative, $(\bar{\eta}^{x+\Delta t_e} - \bar{\eta}^{x-\Delta t_e})/2\Delta t_e$. Nonetheless, we assume this shortcoming to be negligible if the time step is small enough. In a similar manner, we can show that the finite difference counterpart for the energy balance related to $u'\frac{\hat{\bar{u}}u'}{\hat{c}^x}$ (with momentum fluxes given by $F_i = \bar{\tilde{u}}_i \bar{\tilde{u}}_i^x$ and closed boundaries) verifies

$$ \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t0,t1}} u''_{i+1/2,k} (F_{i+1,k} - F_{i,k}) \Delta t \Delta \sigma_k = - \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t0,t1}} \frac{u'^2}{2} \frac{\bar{\eta}^{x+\Delta t} - \bar{\eta}^{x-\Delta t}}{2} \Delta x \Delta \sigma)_{i+1/2,k} \quad (21) $$
For the reasons similar to previously mentioned, (21) does not exactly equilibrate the second term at the RHS of Eq. (18). Some details concerning the use of barotropic currents in the internal mode equations may help in understanding Eq. (21). Let us note that, because of the time-splitting technique, the barotropic current involved in the momentum fluxes of (21) corresponds to a time average of the external mode solution over the period bounded by internal time steps $t - \Delta t$ and $t + \Delta t$ (Mellor, 2003). The transport divergence term, emerging due to the rule of integration by parts (A.1), is then balanced by a variation of the SLA between time steps $t - \Delta t$ and $t + \Delta t$, actually resulting in the RHS of (21).

### 3.1.4. Purely baroclinic terms of advection

Purely baroclinic terms of advection, $u'(\frac{\partial \tilde{u}'}{\partial x'} + \frac{\partial \omega u'}{\partial \sigma})$, do not pose the same problem, since numerous authors (Arakawa and Lamb, 1977; Haltiner and Williams, 1980) have long shown they conserve the global KE, as long as the momentum fluxes and the continuity equation are formulated consistently. Indeed, defining the momentum fluxes as in BM87 ($F_{i,k}^X = \tilde{u}_{i,k}^x \tilde{u}_{i,k}^x$, $F_{i+1/2,k+1/2}^Z = (\bar{\omega}^x \tilde{u}^z)_{i+1/2,k+1/2}$), using the finite difference counterpart of the baroclinic continuity equation (obtained by subtracting Eqs. (5) and (6)), i.e.,

$$ (\tilde{u}'_{i+1/2,k} - \tilde{u}'_{i-1/2,k})/\Delta x = -(\omega_{i,k+1/2} - \omega_{i,k-1/2})/\Delta \sigma_k, \quad (22) $$

applying the numerical rule (A.1) and (A.2), and assuming closed boundaries ($u'_{1/2,k} = u'_{M+1/2,k} = 0$), one can show that the numerical counterpart of the corresponding energy balance, $\sum_{\substack{k=1,R-1 \\ t=0,t1}} u_{i+1/2,k}^t (F_{i+1,k}^X - F_{i-1,k}^X)\Delta\sigma_k + (F_{i+1/2,k+1/2}^Z - F_{i+1/2,k-1/2}^Z)\Delta x\Delta t$, does not change the global KE.

### 3.1.5. External/internal modes energy transfer

Terms in expression (10) that we have not detailed yet are

$$ \int_{t0}^{t1} \int_{x0}^{x1} \left[ \hat{u} \int_{\sigma_0}^{\sigma_1} \frac{\partial \tilde{u}' u'}{\partial x^*} d\sigma + \int_{\sigma_0}^{\sigma_1} u' \left( \frac{\partial \tilde{u}' \hat{u}}{\partial x^*} + \frac{\partial \omega \hat{u}}{\partial \sigma} \right) d\sigma \right] dx^* dt^* \quad (23) $$

Let us recall that the term under the first vertical integral has been obtained by multiplying the non-linear terms at the RHS of (7) by $\hat{u}$ and thus is a contribution to the barotropic KE balance. Conversely, the term under the second integral was obtained by multiplying the non-linear terms (a part of them actually) at the left hand side (hereafter LHS) of (1) by $u'$ and thus is a contribution to the baroclinic energy balance. Using the continuity equations (5) and (6), we see that the latter satisfies

$$ u' \left( \frac{\partial \tilde{u}' \hat{u}}{\partial x^*} + \frac{\partial \omega \hat{u}}{\partial \sigma} \right) = \frac{\partial u' \tilde{u}' \hat{u}}{\partial x^*} - \hat{u} \frac{\partial \tilde{u}' u'}{\partial x^*} \quad (24) $$

If, as previously, we consider a closed domain, the first term at the RHS of (24) disappears when the global balance is calculated. In that case, the second vertical integral of (23) equilibrates the first one, indicating an energy transfer between the external and internal modes, through the advection terms. In order to write a numerical expression for (23), let us recall that the numerical counterpart for the momentum fluxes, $(\tilde{u}'u', \tilde{u}'\hat{u}, \omega\hat{u})$, is based on a common centred scheme, namely: $F_{i,k}^{XE} = \tilde{u}_{i,k}^x \tilde{u}_{i,k}^x$, $F_{i,k}^{XI} = \tilde{u}_{i,k}^x \tilde{u}_{i,k}^x$ and $F_{i+1/2,k+1/2}^{ZI} = (\bar{\omega}^x \hat{u})_{i+1/2,k+1/2}$. We will examine separately the numerical expression of the three terms involved in (23). Assuming closed boundaries, we note that the balance related to the horizontal fluxes satisfies

$$ \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=0,t1}} \hat{u}_{i+1/2}^t \frac{F_{i+1,k}^{XE} - F_{i,k}^{XE}}{\Delta x} \Delta \sigma_k \Delta t \Delta x = - \sum_{\substack{i=2,M-1 \\ k=1,R-1 \\ t=0,t1}} (\tilde{u}^{\prime \prime \prime} x \tilde{u}^{\prime \prime \prime})_{i,k} (\hat{u}_{i+1/2}^t - \hat{u}_{i-1/2}^t) \Delta \sigma_k \Delta t \quad (25) $$

$$ \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=0,t1}} u_{i+1/2,k}^t \frac{F_{i+1,k}^{XI} - F_{i,k}^{XI}}{\Delta x} \Delta \sigma_k \Delta t \Delta x = - \sum_{\substack{i=2,M-1 \\ k=1,R-1 \\ t=0,t1}} (\tilde{u}^{\prime \prime \prime \prime} x \tilde{u}^{\prime \prime \prime \prime})_{i,k} (u_{i+1/2,k}^t - u_{i-1/2,k}^t) \Delta \sigma_k \Delta t \quad (26) $$
A few comments are now helpful. The baroclinic momentum flux, $F_{i,k}^{XE}$, is involved in the coupling of the external and internal modes since it also appears at the RHS of the barotropic momentum equation (7). As we use the time-splitting technique described in Mellor (2003), this term is kept frozen over the external mode sequence bounded by two successive time steps, $t$ and $t + \Delta t$, of the internal mode.

This means that the barotropic current, $\hat{u}$, appearing in expression (25) is actually equivalent to a time average of the external mode solution between internal time steps $t$ and $t + \Delta t$, besides explaining the use of time step $\Delta t$ instead of $\Delta t_c$. In contrast, the barotropic current appearing in the momentum flux $F_{i,k}^{XI}$ is involved in the internal mode computation. Since we use a leapfrog scheme, baroclinic velocities are obtained by integrating momentum equations from time step $t - \Delta t$ to $t + \Delta t$. Some conservation considerations, related to the use of the time-splitting technique, normally affect the way we compute the barotropic component of the current involved in advection terms of the internal mode (Deleersnijder, 1993). As in Mellor (2003), the barotropic current is here provided by the external mode solution, time averaged over the two external mode sequences bounded by internal time steps $t - \Delta t$ and $t + \Delta t$. Therefore, barotropic currents appearing in (25) and (26) are not exactly equivalent: the former have been averaged over one external mode sequence and the latter over two. As we will see now, this prevents the energy transfer between external and internal mode expressed by (23) from being numerically balanced. To achieve this, we temporarily assume that barotropic currents of (25) and (26) were strictly the same. If so, and noting that we have $\overline{u}_{i,k}^{x}(\hat{u}_{i+1/2} - \hat{u}_{i-1/2}) + \overline{\hat{u}}_{i,k}^{x}(u'_{i+1/2,k} - u'_{i-1/2,k}) = u'_{i+1/2,k}\hat{u}_{i+1/2} - u'_{i-1/2,k}\hat{u}_{i-1/2}$, one can show that the sum of (25) and (26) gives

$$ \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t_0,t_1}} (u^{\prime\prime}\hat{u}^{\prime})_{i+1/2,k} \frac{(\overline{\hat{u}^{\prime\prime}}^{x})_{i+1,k} - (\overline{\hat{u}^{\prime\prime}}^{x})_{i,k}}{\Delta x} \Delta\sigma_k \Delta t \Delta x \quad (27) $$

Now, if we consider the vertical flux term of (23) we obtain

$$ \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t_0,t_1}} u''_{i+1/2,k} \frac{F_{i+1/2,k+1/2}^{ZI} - F_{i+1/2,k-1/2}^{ZI}}{\Delta\sigma_k} \Delta\sigma_k \Delta t \Delta x \\ = \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t_0,t_1}} (u''\hat{u}')_{i+1/2,k} \frac{(\overline{\omega}^{\prime\prime})_{i+1/2,k+1/2} - (\overline{\omega}^{\prime\prime})_{i+1/2,k-1/2}}{\Delta\sigma_k} \Delta\sigma_k \Delta t \Delta x \quad (28) $$

According to the baroclinic continuity scheme (22), we can substitute $(\overline{\omega}^{\prime\prime})_{i+1/2,k+1/2} - (\overline{\omega}^{\prime\prime})_{i+1/2,k-1/2}/\Delta\sigma_k = -(\overline{\hat{u}}_{i+1,k}^{x} - \overline{\hat{u}}_{i,k}^{x})/\Delta x$ in the RHS of expression (28), and see that the sum of (27) and (28) is exactly zero. In other words, the equilibrium expressed by (23) and (24) is actually numerically achieved. Unfortunately, this does not happen when barotropic currents appearing in (25) and (26) are different. Nonetheless, we suppose that the value of the barotropic value is not very sensitive to the way the time average is computed (over one or two external mode sequences), notably in low frequency cases. Therefore, the associated energy conservation defect may be regarded as negligible in most cases. We will make this point clear by means of the numerical tests presented in Section 6.

### 3.1.6. Partial conclusion

At this stage, we retain that the discrete form of the KE can be negative because of the leapfrog time-splitting procedure. Furthermore, the time variation term of KE contains a time variation term of SLA that is not exactly balanced by advection terms, evidencing an energy-conserving failure of the leapfrog scheme. Finally, we note that the separation of the model equations in two (external and internal) modes must be dealt with carefully, especially regarding the exchanges between the barotropic and baroclinic forms of the KE, notably those involving advection terms. These processes are not balanced if a classic time-splitting technique is used, because of the inconsistency between the leapfrog scheme used to update internal mode equations, and the forward-like nature of the frozen terms in the external mode equations. We will see that the twin external mode method is a possible alternative to that particular problem.
## 3.2. Dissipation of kinetic energy induced by diffusion terms

### 3.2.1. Analytical formulation

The effect of the diffusion terms on the energy balance is given by

$$ \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} u \left( \frac{\partial}{\partial x^*} \tilde{K}^x \frac{\partial u}{\partial x^*} + \frac{\partial}{\partial \sigma} \frac{K^z}{H} \frac{\partial u}{\partial \sigma} \right) dx^* d\sigma dt^* = \int_{t0}^{t1} \int_{\sigma_0}^{\sigma_1} \left[ u \tilde{K}^x \frac{\partial u}{\partial x^*} \right]_{x0}^{x1} d\sigma dt^* + \int_{t0}^{t1} \int_{x0}^{x1} \left[ u \frac{K^z}{H} \frac{\partial u}{\partial \sigma} \right]_{\sigma_0}^{\sigma_1} dx^* dt^* \\ - \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} \left( \tilde{K}^x \left( \frac{\partial u}{\partial x^*} \right)^2 + \frac{K^z}{H} \left( \frac{\partial u}{\partial \sigma} \right)^2 \right) dx^* d\sigma dt^* \quad (29) $$

The first term on the RHS of (29) represents energy fluxes through lateral boundaries, which vanish in the case of closed boundaries, as we suppose now. The second term represents energy fluxes related to surface and bottom stresses. The last term is always negative, indicating diffusion is an energy sink. Models of the POM-type generally use a turbulence closure scheme based on a prognostic equation for the turbulence kinetic energy (hereafter TKE). The latter normally contains a production term, $\tilde{K}^z(\frac{\partial u}{\partial \sigma})^2$, (BM87; Gaspar et al., 1990), so that the loss of energy induced by vertical mixing in (29) in fact involves an energy transfer process from scales fully resolved by the model grid to the scale of sub-grid turbulence, parameterised by the turbulence closure scheme. Practically, the energy lost by vertical diffusion is counterbalanced by an equivalent gain of TKE. Horizontal turbulence is generally not accurately represented in such models, so that dissipation induced by horizontal mixing, $-\tilde{K}^x(\frac{\partial u}{\partial x^*})^2$, is directly lost by the system.

### 3.2.2. Asselin filter

Although it does not appear in the mathematical formulation of the model equations, we also consider a time diffusion effect, induced by the use of an Asselin filter (Asselin, 1972). The Asselin filter is implemented as in BM87, that is

$$ \phi^{*t} = \phi^t + \frac{\alpha}{2} (\phi^{t+\Delta t} - 2\phi^t + \phi^{*t-\Delta t}) \quad (30) $$

where $\phi^t$ is any of the model variables and $\phi^{*t}$ is the corresponding filtered variable. Combined with the leap-frog scheme, the time-stepping procedure now is

$$ \frac{\phi^{t+\Delta t} - \phi^{*t-\Delta t}}{2\Delta t} = \text{RHS}^t \quad (31) $$

where RHS$^t$ is the numerical counterpart of the various model processes, except the Asselin filter term.

An inconvenience of expression (31) is that we lose, apparently, a basic property of local conservation, since the numerical integral over the simulation duration of the LHS term of (31) is no longer equal to the difference between the final and initial states. In fact, expression (31) can be rewritten using (30) at time $t - \Delta t$, namely

$$ \phi^{t+\Delta t} - \phi^{t-\Delta t} = 2\Delta t \text{RHS}^t + \frac{\alpha}{2} (\phi^t - 2\phi^{t-\Delta t} + \phi^{*t-2\Delta t}) \quad (32) $$

where the filtered variable, $\phi^{*t-2\Delta t}$, can in turn be replaced by $\phi^{t-2\Delta t}$ owing to (30), thus introducing the additional term $(\frac{\alpha}{2})^2(\phi^{t-\Delta t} - 2\phi^{t-2\Delta t} + \phi^{*t-3\Delta t})$ on the RHS of (32). This operation can be repeated for $\phi^{*t-3\Delta t}$ and so on, so that we finally obtain that the RHS of (31) is equivalent to $2\Delta t \text{RHS}^t + \sum_{q=1,q_{\text{max}}}(\frac{\alpha}{2})^q(\phi^{t-(q-1)\Delta t} - 2\phi^{t-q\Delta t} + \phi^{t-(q+1)\Delta t})$ where $q_{\text{max}}$ is the total number of iterations since the beginning of the simulation. As the Asselin coefficient, $\alpha$, generally is a positive number much smaller than 1 (BM87), the terms of the polynomial series with $q \ge 2$ can be neglected, so that (32) is virtually equivalent to

$$ \phi^{t+\Delta t} - \phi^{t-\Delta t} = 2\Delta t \text{RHS}^t + \frac{\alpha}{2} (\phi^t - 2\phi^{t-\Delta t} + \phi^{t-2\Delta t}) \quad (33) $$

We will finally retain (33), because its conservation properties are more interesting than the ones of the classic Asselin scheme. Indeed, the numerical equivalent of the time integral of (33) gives
$$ \frac{1}{2} (\phi^{t+\Delta t} + \phi^{t1} - \frac{\alpha}{2}(\phi^{t1} - \phi^{t1-\Delta t})) = \sum_{t=t0,t1} \text{RHS}^t \Delta t + \frac{1}{2} (\phi^{t0-\Delta t} + \phi^{t0} - \frac{\alpha}{2}(\phi^{t0-\Delta t} - \phi^{t0-2\Delta t})) \quad (34) $$

showing that the final state of the variables now only depends on the initial state and the time integral of the forcing terms, as may be expected. Also note that expression (33) cannot be implemented as such, due to the sigma formulation of our model, so that we actually use the following scheme:

$$ \tilde{\phi}^{t+\Delta t} - \tilde{\phi}^{t-\Delta t} = 2 \mathrm{d}t \mathrm{R} \tilde{\mathrm{H}} \mathrm{S}^t + \tilde{\mathrm{F}}_{\mathrm{ass}}^{t-0.5\Delta t} - \tilde{\mathrm{F}}_{\mathrm{ass}}^{t-1.5\Delta t} \quad (35) $$

with

$$ (\tilde{F}_{\text{ass}}^{t-0.5\Delta t}, \tilde{F}_{\text{ass}}^{t-1.5\Delta t}) = \frac{\alpha}{2} \left( \frac{H^t + H^{t-\Delta t}}{2} (\phi^t - \phi^{t-\Delta t}), \frac{H^{t-\Delta t} + H^{t-2\Delta t}}{2} (\phi^{t-\Delta t} - \phi^{t-2\Delta t}) \right) \quad (36) $$

The Asselin filter in its standard form is known to destroy conservation properties for tracer variables (e.g. Griffies et al., 2001). The advantage of (35) and (36) is that the filter leaves homogeneous fields unchanged despite possible motions of the SLA. In fact, the filter can be seen as a finite difference counterpart for a time diffusion process, $\frac{\partial}{\partial t} \tilde{K}^A \frac{\partial \phi}{\partial t}$ (with $\tilde{K}^A/\Delta t^2 = \tilde{\alpha}/2$), comparable to horizontal and vertical mixing processes. Regarding velocity, for example, we expect that the KE dissipated by the filter is given by the integral over time of $-\tilde{K}^A(\frac{\partial u}{\partial t})^2$. Now, let us note that the computation of the filter is centred on the previous time step, $t - \Delta t$. The energy balance approach emphasises the interest of this particular feature. Indeed, using A1, the KE balance related to the filter satisfies

$$ \sum_{t=t0,t1} u^t (\tilde{F}_{\text{ass}}^{t-0.5\Delta t} - \tilde{F}_{\text{ass}}^{t-1.5\Delta t}) \mathrm{d}t = u^{t1} \tilde{F}_{\text{ass}}^{t1-0.5\Delta t} - u^{t0} \tilde{F}_{\text{ass}}^{t0-1.5\Delta t} - \sum_{t=t0,t1-\Delta t} \tilde{K}^A \frac{(u^t - u^{t-\Delta t})(u^{t+\Delta t} - u^t)}{\Delta t} \quad (37) $$

We clearly see that the product $(u^t - u^{t-\Delta t})(u^{t+\Delta t} - u^t)$, at the RHS of (37), should be positive most of the time (so that the global effect of the filter is to dissipate energy, as expected). However, it is possibly negative, for instance if the simulation develops some high-frequency instabilities reversing the signs of the current at successive time steps. Although we might conclude that the KE produced by the filter amplifies the numerical noise, in fact it does not. Indeed, in accordance with (17) and (18), we have seen that KE of high-frequency instabilities is negative, so that the apparent production of energy by the filter reduces the level of numerical noise. We would not have this interesting property if the computation of the filter was centred on the present time step, $t$, as an intuitive approach would suggest.

### 3.2.3. Horizontal and vertical mixing

These considerations also apply to the other diffusion processes. Horizontal mixing is therefore computed with velocities obtained at the previous time step $t - \Delta t$. The related KE balance, supposing a closed domain and using (A.1), satisfies

$$ \begin{aligned} & \sum_{\substack{t=t0,t1 \\ k=1,R-1 \\ i=1,M-1}} \frac{u^t}{\Delta x} (\tilde{K}^x(u^{t-\Delta t}_{i+1/2,k} - u^{t-\Delta t}_{i+1/2,k}) - \tilde{K}^x(u^{t-\Delta t}_{i+1/2,k} - u^{t-\Delta t}_{i-1/2,k}))\Delta\sigma_k \Delta t \\ & = -\sum_{\substack{t=t0,t1 \\ k=1,R-1 \\ i=2,M-1}} \frac{\tilde{K}^x}{\Delta x} (u^{t-\Delta t}_{i+1/2,k} - u^{t-\Delta t}_{i-1/2,k})(u^t_{i+1/2,k} - u^t_{i-1/2,k})\Delta\sigma_k \Delta t \end{aligned} \quad (38) $$

In most cases, the RHS of (38) is negative, evidencing the expected dissipation effect of the horizontal mixing. In the particular case of a numerical mode having a $2\Delta t$ period the RHS of (38) becomes positive, counteracting the development of the former, since KE of numerical noise is negative. Similar comments apply for vertical mixing, except that its computation is now centred on the following time step, $t + \Delta t$ (BM87). As an inconvenient consequence, at a given grid node $(i+1/2,k)$, any model equation contains several unknown variables inducing a vertical dependence on the other equations, requiring a linear system with the equations
of a same grid column be solved. On the other hand, implicit schemes are known to be more stable. Using A1, the energy balance associated with vertical mixing verifies

$$ 
\begin{aligned}
& \sum_{\substack{t=t_0,t1 \\ k=1,R-1 \\ i=1,M-1}} u_{i+1/2,k}^t \left( \frac{\tau_{i+1/2,k+1/2} - \tau_{i+1/2,k-1/2}}{\Delta\sigma_k} \right) \Delta x \Delta t \Delta \sigma_k \\
&= \sum_{\substack{t=t_0,t1 \\ i=1,M-1}} (u_{i+1/2,R-1}^t \tau_{i+1/2}^{\text{Surf}} - u_{i+1/2,1}^t \tau_{i+1/2}^{\text{Bot}}) \Delta x \Delta t - \sum_{\substack{t=t_0,t1 \\ k=1,R-2 \\ i=1,M-1}} K^z \frac{u_{i+1/2,k+1}^{t+\Delta t} - u_{i+1/2,k}^{t+\Delta t}}{z_{i+1/2,k+1} - z_{i+1/2,k}} (u_{i+1/2,k+1}^t - u_{i+1/2,k}^t) \Delta x \Delta t
\end{aligned} 
\quad (39) $$

where $\tau_{i+1/2,k+1/2} = K^z(u_{i+1/2,k+1}^{t+\Delta t} - u_{i+1/2,k}^{t+\Delta t})/(z_{i+1/2,k+1} - z_{i+1/2,k})$ and $(\tau_{i+1/2}^{\text{Surf}}, \tau_{i+1/2}^{\text{Bot}})$ are the surface and bottom stresses, acting as source or sink of energy. The last term at the RHS of (39) is a term of energy dissipation.

For the sake of simplicity we have not detailed the separated contribution of internal and external mode diffusion terms to the energy balance. This can easily be deduced from Eqs. (29) to (39), since we use a constant mixing coefficient, $K^x$, so that the depth-averaged horizontal diffusion process of Eq. (1) is simply equivalent to the horizontal diffusion of the depth-averaged current. In other words, we verify that

$$ \int_{\sigma_0}^{\sigma_1} \left( \frac{\partial}{\partial x^*} \tilde{K}^x \frac{\partial u}{\partial x^*} \right) d\sigma = \frac{\partial}{\partial x^*} \tilde{K}^x \frac{\partial \hat{u}}{\partial x^*} \quad (40) $$

In case of a non-constant mixing coefficient, Eq. (7) should take into account an additional term, $\int_{\sigma_0}^{\sigma_1} (\frac{\partial}{\partial x^*} \tilde{K}' \frac{\partial \hat{u}}{\partial x^*}) d\sigma$, with $K' = K^x - \int_{\sigma_0}^{\sigma_1} K^x d\sigma$.

### 3.2.4. Partial conclusion

At this stage, we retain that the Asselin filter can be approximated as a time diffusion process, and thus considered as a sink of energy. We also retain that a backward centring helps the scheme removing high-frequency instabilities. For the same reasons, horizontal and vertical mixing computations are respectively centred on the previous and following time steps.

## 3.3. Pressure gradient

### 3.3.1. SLA slope

The pressure gradient is given by expression (3). Note that the latter contains a purely barotropic term associated to the slope of the SLA. The contribution to the barotropic KE is obtained by multiplying the SLA slope term, $gH \frac{\partial \eta}{\partial x^*}$, of Eq. (7) by the barotropic current, $\hat{u}$. Using the barotropic continuity equation (6), and the divergence theorem, we obtain the related energy balance:

$$ - \int_{t_0}^{t1} \int_{x0}^{x1} g \hat{u} H \frac{\partial \eta}{\partial x} dx^* dt^* = - \int_{x0}^{x1} \left[ \frac{g \eta^2}{2} \right]_{t0}^{t1} dx^* - \int_{t0}^{t1} g [H \hat{u} \eta]_{x0}^{x1} dt^* \quad (41) $$

where the first term at the RHS of (41) represents a variation of barotropic PE between initial and final times, $t0$ and $t1$. The second term at the RHS of (41) vanishes in case of closed boundaries. The discrete form of the surface pressure gradient is simply $-g\tilde{H}_{i+1/2}^x(\eta_{i+1} - \eta_i)/\Delta x$. Energy conservation is easily demonstrated (e.g. Arakawa and Lamb, 1977) using the discrete barotropic continuity equation (19) and the numerical rule of integration by parts (A.1). In the case of closed boundaries the finite difference counterpart of (41) is

$$ - \sum_{\substack{t=t_0,t1 \\ i=1,M-1}} g \tilde{u}_{i+1/2}^t \frac{\eta_{i+1}^t - \eta_i^t}{\Delta x} \Delta t_e = - \sum_{i=2,M-1} \frac{g}{2} (\eta_i^{t1} \eta_i^{t1+\Delta t_e} - \eta_i^{t0} \eta_i^{t0-\Delta t_e}) \Delta x \quad (42) $$

where the numerical counterpart of the barotropic PE anomaly is given by $g\eta^t \eta^{t+\Delta t_e}/2$. While the continuous form of the latter is always positive, we see that it may become negative in certain cases (in particular if the model becomes unstable and develops numerical noise), due to the leapfrog time-stepping scheme. This issue
may be taken into account in order to obtain a consistent construction of the other numerical schemes in the external mode, such as the boundary conditions. In case of open boundaries, the time integral of the PE flux through the boundaries should be added to the RHS of (42), $g(\overline{H}^x \hat{u}')_{1+1/2} \eta'_1 - g(\overline{H}^x \hat{u}'')_{M-1/2} \eta'_M$. Using a radiative boundary condition that is decentred in time and space as described in Marsaleix et al. (2006), these energy fluxes become $-g c \eta'_1^{t+\Delta t_e} \eta'_1 - g c \eta'_M^{t+\Delta t_e} \eta'_M$, where c is the phase velocity of surface waves. Clearly, in most cases, these fluxes are negative, since $\eta'_t$ and $\eta'^{t+\Delta t_e}$ share the same sign. In other words, the radiative boundary conditions in question subtract energy from the system, a property that may be considered to enhance the model stability. When case the model develops numerical noise, $\eta'_t$ and $\eta'^{t+\Delta t_e}$ are of opposite sign, so that the discrete form of the PE anomaly is negative. The fluxes through the open boundaries remain consistent, since they become positive, and help reduce the noise level.

### 3.3.2. Baroclinic contribution of the density anomaly gradient

We consider the pressure gradient terms depending on density gradients, which are the two last terms at the RHS of Eq. (3). These terms both contribute to the baroclinic and barotropic energy balance. Let us first examine the baroclinic contribution which is obtained by multiplying, in Eq. (1), the part of the pressure that depends on density gradients by the baroclinic current, $u'$. Using the divergence theorem, the hydrostatic equilibrium and the baroclinic continuity equations, we can show that we have

$$ -\int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} u' \left( \frac{\partial p'}{\partial x^*} + g \frac{\partial z}{\partial x^*} \rho' \right) H \, d\sigma \, dx^* \, dt^* = -g \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} \left( \omega \rho' H + \frac{\partial z}{\partial x^*} \tilde{u}' \rho' \right) d\sigma \, dx^* \, dt^* \quad (43) $$

We will see in the following section that PE balance related to the advection of density contains the same terms as those at the RHS of (43), thus evidencing a mechanism of conversion between the kinetic and potential forms of energy. Here, the pressure gradient force is given by a straightforward discretisation, namely

$$ \frac{\partial p'}{\partial x^*} + g \frac{\partial z}{\partial x^*} \rho' = \frac{p'_{i+1,k} - p'_{i,k}}{\Delta x} + g \frac{z_{i+1,k} - z_{i,k}}{\Delta x} \rho'^{t,x}_{i+1/2,k} \quad (44) $$

with the hydrostatic pressure computed according to

$$ p'_{i,k} = gH_i \left[ \sum_{q=k+1,R-1} \rho'_{i,q} \Delta\sigma_q + 0.5 \rho'_{i,k} \Delta\sigma_k \right] \quad (45) $$

Let us recall that numerous authors have considered the discretisation of the PGF. These studies were often motivated by the problem of truncation errors induced by the sigma-coordinate scheme. Using a straightforward scheme similar to the one used here, Beckmann and Haidvogel (1993) and Haney (1991) showed that the PGF accuracy can be improved by eliminating a mean density profile before calculating the pressure gradient. However, this procedure is not retained here since it is known that this destroys the energetic consistency of the PGF with vertical density advection (Gerdes, 1993). An alternative to the schemes (44) and (45) is the second-order POM density Jacobian. Shchepetkin and McWilliams (2003) showed that this scheme is algebraically equivalent to (44) provided that the rectangular integration in (45) is replaced by a trapezoidal integration. As a consequence the density Jacobian PGF vanishes identically, regardless the vertical structure of the grid, if density is simply a linear function of z, a property that is not exactly satisfied by (44) and (45). On the other hand, the rectangular integration (45) does not require any particular boundary condition at the top of the grid while the density Jacobian scheme does (Shchepetkin and McWilliams, 2003). The energetic consistency of the density Jacobian PGF with density advection terms is discussed in Song and Wright (1998). Their study, however, neglected the effects of sea surface elevation variations. This was done by Shchepetkin and McWilliams (2003) who showed that the energetic consistency of the barotropic part of the PGF with the time differentiation term in the discrete density equation was possibly compromised by free-surface variations. Opposingly we will see in Section 4.1 that the rectangular integration in (45), that assumes a piece-wise constant distribution within each vertical cell, is consistent with the definition of the discrete PE. The scheme (44) and (45) was actually chosen for this particular reason but we keep in mind that the density Jacobian scheme is virtually more accurate.
Let us now turn to the numerical counterpart of expression (43). It follows from (45) that

$$p'_{i,k} - p'_{i,k-1} = -\frac{gH_i}{2}(\Delta\sigma_k \rho'_{i,k} + \Delta\sigma_{k-1} \rho'_{i,k-1}) \quad (46)$$

Using (44) and (46), the continuity equation (22) and the closed boundary condition, $u_{1/2,k} = u_{M+1/2,k} = 0$, we can show that the numerical counterpart expression (43) verifies

$$\begin{aligned} & -\sum_{\substack{t=0,t1 \\ i=1,M-1 \\ k=1,R-1}} u'_{i+1/2,k} \left( \frac{p'_{i+1,k} - p'_{i,k}}{\Delta x} + g \frac{z_{i+1,k} - z_{i,k}}{\Delta x} \overline{\rho'}^x_{i+1/2,k} \right) \overline{H}^x_{i+1/2} \Delta \sigma_k \Delta x \Delta t \\ &= -\frac{g}{2} \sum_{\substack{t=0,t1 \\ i=1,M \\ k=2,R-1}} \omega_{i,k} (\Delta \sigma_k \rho'_{i,k} + \Delta \sigma_{k-1} \rho'_{i,k-1}) H_i \Delta x \Delta t - g \sum_{\substack{t=t0,t1 \\ i=1,M-1 \\ k=1,R-1}} \tilde{u}'_{i+1/2,k} (z_{i+1,k} - z_{i,k}) \overline{\rho'}^x_{i+1/2,k} \Delta \sigma_k \Delta t \end{aligned} \quad (47)$$

The RHS of (47) will be compared to the corresponding term of the PE balance in the following section.

### 3.3.3. Barotropic contribution of the density anomaly gradient

Here, we examine contribution of pressure gradient force associated with density anomaly to the barotropic kinetic energy balance. The latter is obtained by multiplying, in the barotropic momentum equation (7), the part of the pressure gradient depending on density gradients by the barotropic current $\hat{u}$. Using the property of the integral of a divergence over a closed domain and the barotropic continuity equation (6), we can show that we have

$$-\int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} \hat{u} \left( \frac{\partial p'}{\partial x^*} + g \frac{\partial z}{\partial x^*} \rho' \right) H \, d\sigma \, dx^* \, dt^* = -\int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma_0}^{\sigma_1} \left( p' \frac{\partial \eta}{\partial t} + g \frac{\partial z}{\partial x^*} \tilde{u} \rho' \right) d\sigma \, dx^* \, dt^* \quad (48)$$

We will see in the following section, focussing on the PE balance, that the RHS of (48) evidences a mechanism of conversion between the kinetic and potential forms of energy. Now, we examine the numerical counterpart of the LHS of (48). Using the closed boundary condition, $\hat{u}_{1/2} = \hat{u}_{M+1/2} = 0$, the latter verifies

$$\begin{aligned} & -\sum_{\substack{t=t0,t1 \\ i=1,M-1 \\ k=1,R-1}} \hat{u}_{i+1/2} \left( \frac{p'_{i+1,k} - p'_{i,k}}{\Delta x} + g \frac{z_{i+1,k} - z_{i,k}}{\Delta x} \overline{\rho'}^x_{i+1/2,k} \right) \overline{H}^x_{i+1/2} \Delta \sigma_k \Delta x \Delta t \\ &= \sum_{\substack{t=t0,t1 \\ i=1,M \\ k=1,R-1}} p'_{i,k} \frac{\tilde{u}_{i+1/2} - \tilde{u}_{i-1/2}}{\Delta x} \Delta \sigma_k \Delta x \Delta t - g \sum_{\substack{t=t0,t1 \\ i=1,M-1 \\ k=1,R-1}} \tilde{u}_{i+1/2} (z_{i+1,k} - z_{i,k}) \overline{\rho'}^x_{i+1/2,k} \Delta \sigma_k \Delta t \end{aligned} \quad (49)$$

We can see that the RHS of (49) represents a numerical counterpart of the RHS of (48), especially if we replace the velocity gradient by the time variation of SLA given by the barotropic continuity equation (19). However, this point needs some comments concerning the coupling of the internal and external modes. The barotropic KE balance is indeed based on the external mode equations, the latter being computed with a time step smaller than the time step of the internal mode. The barotropic current appearing at the RHS of (49) is combined with some terms provided by the internal mode, $(p', \rho', z)$, which are kept frozen over a sequence of computation of the external mode, bounded by the internal time steps $t$ and $t + \Delta t$. In other words, the numerical time integral of (49), that should normally use the time-stepping scheme of the external mode, can simply be computed with the time-stepping scheme of the internal mode (as for the numerical integral of (47)), insofar as the barotropic current value has been averaged over the duration of the external mode sequence bounded by internal time steps $t$ and $t + \Delta t$. If so, the continuity equation (19) requires that the barotropic current divergence equilibrates the time variation of the SLA between time step $t$ and $t + \Delta t$. In other words, the first term at the RHS of (49) can be rewritten as
$$ \sum_{\substack{t=t_0,t_1 \\ i=1,M \\ k=1,R-1}} p'_{i,k} \frac{\tilde{u}_{i+1/2} - \tilde{u}_{i-1/2}}{\Delta x} \Delta \sigma_k \Delta x \Delta t = - \sum_{\substack{t=t_0,t_1 \\ i=1,M \\ k=1,R-1}} p''_{i,k} \frac{\eta_i^{t+\Delta t} - \eta_i^t}{\Delta t} \Delta \sigma_k \Delta x \Delta t \quad (50) $$

### 3.3.4. Partial conclusion

At this stage we retain, according to the relationship (41), that the variations of barotropic KE induced by the slope of the SLA must be equilibrated by an equivalent variation of the barotropic PE, of which the discrete form may become negative due to the leapfrog scheme. We have also seen that this process is exactly reproduced by the model, numerically speaking, provided that the barotropic PE is correctly defined, as stated by expression (42).

## 4. Potential energy balance

As far as the effects of density are concerned, a PE balance is obtained by multiplying the density equation (4) by gz (Winters et al., 1995). Considering a closed domain and no turbulent fluxes through the bottom and the surface, this gives

$$ \int_{t_0}^{t_1} \int_{x_0}^{x_1} \int_{\sigma_0}^{\sigma_1} \frac{\partial gz \tilde{\rho}'}{\partial t^*} d\sigma dx^* dt^* = \int_{t_0}^{t_1} \int_{x_0}^{x_1} \int_{\sigma_0}^{\sigma_1} \left( g \tilde{\rho}' \frac{\partial z}{\partial t^*} + g \rho' \tilde{u} \frac{\partial z}{\partial x^*} + g \rho' \omega H \right) d\sigma dx^* dt^* \\ - \int_{t_0}^{t_1} \int_{x_0}^{x_1} \int_{\sigma_0}^{\sigma_1} \left( g \tilde{K}_T^x \frac{\partial \rho'}{\partial x^*} \frac{\partial z}{\partial x^*} + g K_T^z \frac{\partial \rho'}{\partial \sigma} \right) d\sigma dx^* dt^* \quad (51) $$

The LHS of (51) is a time variation term of PE. The first integral at the RHS of (51) represents the effect of density advection on the PE balance. Considering that in the sigma-coordinate system the vertical velocity is expressed as $w = \frac{\partial z}{\partial t^*} + u \frac{\partial z}{\partial x^*} + \omega$, we recover the classic term, $-g\rho'w$ (Gill, 1982), expressing the effect of vertical motions on the PE level, as expected. The second integral at the RHS of (51) represents the effects of mixing.

### 4.1. Time variation term

The time variation term of PE, $\frac{\partial gz \tilde{\rho}'}{\partial t^*}$, and the term related to the variation of the vertical grid levels, $gz \frac{\partial \tilde{\rho}'}{\partial t^*}$, both come from the time variation term of the density equation (4), since we have $gz \frac{\partial \tilde{\rho}'}{\partial t^*} = \frac{\partial gz \tilde{\rho}'}{\partial t^*} - g \tilde{\rho}' \frac{\partial z}{\partial t^*}$. As far as the PE balance is concerned, the corresponding numerical counterpart verifies

$$ \sum_{\substack{t=t_0,t_1 \\ i=1,M \\ k=1,R-1}} gz'_{i,k} \frac{\tilde{\rho}''_{i,k}^{t+\Delta t} - \tilde{\rho}''_{i,k}^{t-\Delta t}}{2\Delta t} \Delta t \Delta \sigma_k \Delta x = \sum_{\substack{i=1,M \\ k=1,R-1}} \frac{g}{2} \left( z_{i,k}^{t1} \tilde{\rho}''_{i,k}^{t1+\Delta t} + z_{i,k}^{t1-\Delta t} \tilde{\rho}''_{i,k}^{t1-\Delta t} - z_{i,k}^{t0} \tilde{\rho}''_{i,k}^{t0-\Delta t} - z_{i,k}^{t0+\Delta t} \tilde{\rho}''_{i,k}^{t0+\Delta t} \right) \Delta \sigma_k \Delta x \\ - \sum_{\substack{i=t_0+\Delta t,t_1-\Delta t \\ i=1,M \\ k=1,R-1}} g H'_{i,k} \rho''_{i,k} \sigma_k \frac{\eta_i^{t+\Delta t} - \eta_i^{t-\Delta t}}{2\Delta t} \Delta t \Delta \sigma_k \Delta x \quad (52) $$

The first term at the RHS of (52) represents the difference of PE between the beginning, $\frac{g}{2}(z_{i,k}^{t0} \tilde{\rho}''_{i,k}^{t0-\Delta t} + z_{i,k}^{t0+\Delta t} \tilde{\rho}''_{i,k}^{t0+\Delta t})$, and the end, $\frac{g}{2}(z_{i,k}^{t1} \tilde{\rho}''_{i,k}^{t1-\Delta t} + z_{i,k}^{t1+\Delta t} \tilde{\rho}''_{i,k}^{t1+\Delta t})$, of the simulation. In the second term of the RHS of (52), we have used the sigma-coordinate transformation (11) to express the variations of grid levels as a function of the SLA, namely

$$ \frac{\partial z}{\partial t^*} = \sigma \frac{\partial \eta}{\partial t^*} \quad (53) $$

Also recall that we have $\tilde{\rho}' = H\rho'$. Expressed thus, the second term at the RHS of (52) allows to better evidence a mechanism of energy conversion involving the pressure gradient since, using the hydrostatic equilibrium and (53), one can show that the term related to grid level variations actually satisfies
$$ \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma0}^{\sigma1} g \tilde{\rho}' \frac{\partial z}{\partial t^*} d\sigma dx^* dt^* = \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma0}^{\sigma1} p' \frac{\partial \eta}{\partial t^*} d\sigma dx^* dt^* \quad (54) $$

We see that the RHS of (54) also appears in expression (48), but with an opposite sign. Additionally, (54) shows that the first term at the RHS of (51) should be equilibrated by the contribution of the pressure anomaly gradient to the baroclinic and barotropic kinetic energies, namely (43) plus (48). Let us now examine the numerical counterpart of (54), in order to compare it to the corresponding numerical counterpart of Eq. (48), which is (50).

Combining (45) and (14) we obtain

$$ \sum_{k=1,R-1} p_k'' \Delta \sigma_k = g H^t \sum_{k=1,R-1} \rho_k'' \sigma_k \Delta \sigma_k \quad (55) $$

It follows from (55) that the RHS of (50) actually verifies

$$ - \sum_{\substack{t=t0,t1 \\ i=1,M \\ k=1,R-1}} p_{i,k}'' \frac{\eta_i^{t+\Delta t} - \eta_i^t}{\Delta t} \Delta \sigma_k \Delta x \Delta t = -g \sum_{\substack{t=t0,t1 \\ i=1,M \\ k=1,R-1}} H_i^t \rho_{i,k}'' \sigma_k \frac{\eta_i^{t+\Delta t} - \eta_i^t}{\Delta t} \Delta t \Delta \sigma_k \Delta x \quad (56) $$

Clearly, the RHS of (56) is not strictly equivalent to the second term at the RHS of (52), since the time derivative of the SLA is computed with values obtained at time steps $t$ and $t + \Delta t$ for the former, and time steps $t - \Delta t$ and $t + \Delta t$ for the latter. Therefore, the equilibrium expressed by (54) cannot be faithfully reproduced by the model because of the inconsistency between, on one the hand, the leapfrog time-stepping procedure used for the density equation, and the forward-like nature of the frozen terms coupling the internal and external modes on the other hand. We note that except the time differentiation term of SLA, the RHS of (56) and (52) are identical, as a result of the energetic consistency of the rectangular scheme used to compute the hydrostatic pressure (45) with the definition of the discrete PE.

## 4.2. Density advection terms

This part of the balance concerns the terms $g p' \tilde{u} \frac{\partial z}{\partial x^*} + g p' \omega H$ appearing at the RHS of (51), which stem from density advection, $\frac{\partial \tilde{u} \rho'}{\partial x^*} + \frac{\partial \omega \rho'}{\partial \sigma}$. In order to build the numerical counterpart of the corresponding balance, we use the classic centred advection scheme of the model. Meanwhile, the current is expressed as the sum of its external and internal parts, ($\tilde{u} = \hat{u} + \tilde{u}'$), in order to separate the barotropic and baroclinic contributions. Using (15), we obtain the corresponding PE balance:

$$ \begin{aligned} & \sum_{\substack{t=t0,t1 \\ i=1,M \\ k=1,R-1}} g z_{i,k} \left( \frac{F_{i+1/2,k}^x - F_{i-1/2,k}^x}{\Delta x} + \frac{F_{i,k+1/2}^z - F_{i,k-1/2}^z}{\Delta \sigma_k} \right) \Delta \sigma_k \Delta x \Delta t \\ & = - \sum_{\substack{t=t0,t1 \\ i=1,M-1 \\ k=1,R-1}} g \tilde{u}_{i+1/2} \overline{\rho'}_{i+1/2,k}^x (z_{i+1,k} - z_{i,k}) \Delta \sigma_k \Delta t - \sum_{\substack{t=t0,t1 \\ i=1,M-1 \\ k=1,R-1}} g \tilde{u}'_{i+1/2,k} \overline{\rho'}_{i+1/2,k}^x (z_{i+1,k} - z_{i,k}) \Delta \sigma_k \Delta t \\ & \quad - \sum_{\substack{t=t0,t1 \\ i=1,M \\ k=2,R-1}} g \omega_{i,k} \frac{\Delta \sigma_k \rho'_{i,k} + \Delta \sigma_{k-1} \rho'_{i,k-1}}{2} H_i \Delta x \Delta t \end{aligned} \quad (57) $$

with $F_{i+1/2,k}^x = (\tilde{u} \overline{\rho'}^x)_{i+1/2,k}$ and $F_{i,k+1/2}^z = \omega_{i,k+1/2} \frac{\Delta \sigma_k \rho'_{i,k} + \Delta \sigma_{k+1} \rho'_{i,k+1}}{\Delta \sigma_k + \Delta \sigma_{k+1}}$ the advection fluxes. In the vertical advection flux, the interface density value is given by a particular weighted average in order to maintain a consistent energy balance with the hydrostatic pressure scheme (45). Let us note, that, from the point of view of the
advection problem, this is a priori not the best estimate of the mid-point value, $\rho_{i,k+1/2}$, since a second-order scheme would rather reverse the $\Delta\sigma$'s order, e.g. $\rho_{i,k+1/2} = (\Delta\sigma_{k+1}\rho'_{i,k} + \Delta\sigma_k\rho'_{i,k+1})/(\Delta\sigma_k + \Delta\sigma_{k+1})$.

Thus, we see that, as expected, the two last terms of (57) exactly equilibrate the numerical counterpart of (43), i.e., (47).

The first term at the RHS of (57) is similar to the second term at the RHS of (49), but the expected equilibrium is not exactly achieved. On one hand, the barotropic current appearing in (49) corresponds to a time averaged value of the external sequence bounded by time steps $t$ and $t + \Delta t$. On the other hand, the barotropic current in (57) corresponds to a time averaged value of the two external sequences bounded by time steps $t - \Delta t$ and $t + \Delta t$. As discussed in previous sections, this is the consequence, first of the forward-like nature of the frozen terms in the external mode equations, and second of the leapfrog time-stepping procedure used to compute the density equation.

### 4.3. Partial conclusion

At this stage, let us recapitulate the results just obtained; thanks to (43), (48) and (54) we know that the KE variations related to the pressure anomaly gradient must equilibrate the PE variations related to density advection represented by the first term at the RHS of (51), namely

$$ - \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma0}^{\sigma1} u \left( \frac{\partial p'}{\partial x^*} + g \frac{\partial z}{\partial x^*} \rho' \right) H \, d\sigma \, dx^* \, dt^* = \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma0}^{\sigma1} \left( g \tilde{\rho}' \frac{\partial z}{\partial t^*} + g \rho' \tilde{u} \frac{\partial z}{\partial x^*} + g \rho' \omega H \right) d\sigma \, dx^* \, dt^* \quad (58) $$

In order to check whether this property is numerically reproduced by the model, we first calculated the numerical counterpart of the LHS of (58), considering the barotropic and baroclinic contribution separately, leading to numerical expressions (47) and (49). Then, we calculated the numerical counterpart of the RHS of (58), considering the terms $g\tilde{\rho}'\frac{\partial z}{\partial t^*}$ and $g\rho'\tilde{u}\frac{\partial z}{\partial x^*} + g\rho'\omega H$ separately, arriving at numerical expressions (52) and (57). At first sight, expressions (47), (49), (52) and (57) seemed numerically consistent, but extensive examination revealed that the balance expressed by (58) is in fact not exactly recovered, because of the method used to couple the external and internal modes. Indeed, the barotropic current of (49) has been averaged over one external sequence only, when a perfect equilibrium would actually require a time average over two sequences bounded by time steps $t - \Delta t$ and $t + \Delta t$.

### 4.4. Density diffusion terms

The contribution of density diffusion terms to the PE balance is obtained by multiplying the mixing terms by $gz$. Taking into account that we suppose lateral, surface and bottom fluxes to be zero, the corresponding balance satisfies

$$ \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma0}^{\sigma1} gz \left( \frac{\partial}{\partial x^*} \tilde{K}_T^x \frac{\partial \rho'}{\partial x^*} + \frac{\partial}{\partial \sigma} \frac{K_T^z}{H} \frac{\partial \rho'}{\partial \sigma} \right) d\sigma \, dx^* \, dt^* = \int_{t0}^{t1} \int_{x0}^{x1} \int_{\sigma0}^{\sigma1} \left( -g \tilde{K}_T^x \frac{\partial \rho'}{\partial x^*} \frac{\partial z}{\partial x^*} - g K_T^z \frac{\partial \rho'}{\partial \sigma} \right) d\sigma \, dx^* \, dt^* \quad (59) $$

The last term in brackets on the RHS of (59) is generally positive, since the vertical density gradient, $\frac{\partial\rho'}{\partial\sigma}$, is generally negative. This means that vertical mixing tends to increase the PE level, in other words: energy is consumed by vertical mixing. In general, this term is retrieved, with an opposite sign, in turbulence closure schemes based on a prognostic equation of the TKE (Mellor and Yamada, 1982; Gaspar et al., 1990). This highlights the fact that the energy amount required to mix water should be taken from the TKE tank. As far as our turbulence scheme is concerned (Gaspar et al., 1990), the vertical mixing coefficient is proportional to the square root of the TKE and thus mixing does not occur if the TKE level is low, and inversely. The first term between brackets at the RHS of (59) corresponds to the effect of horizontal mixing. Strictly speaking, the word “horizontal” is not appropriate since the slope of iso-sigma levels, $\frac{\partial z}{\partial x^*}$, causes this mixing term to contain a part of vertical mixing. As a consequence, this “horizontal” mixing induces a variation of the PE level. As far
as we know, the latter is not taken into account by the turbulence closure schemes mentioned previously. This point will be discussed in the following. These considerations on mixing processes may raise questions regarding the advection scheme. Indeed, the centred form of the density advection scheme presented in (57) is not always suitable, for instance because of numerical instabilities that can possibly develop in the wake of frontal structures (Reffray et al., 2004). As suggested by some authors, it is possible to counter this drawback by adding a diffusion term to the advection scheme, which prevents the development of numerical noise, using for instance a Lax-Wendroff method or a TVD scheme (James, 1996). This diffusion operator has a form similar to the scheme presented in Section 3.2, except that the mixing coefficient computation depends on numerical concerns, such as the limitation of overshoot effects, rather than on physical considerations regarding turbulence. This additional diffusion term is likely to induce significant vertical mixing, which may eventually raise the question of a corresponding balance in the TKE scheme. In other words, whether any vertical diffusion process should be equilibrated in the TKE equation by an equivalent buoyancy term, which will be discussed in the following.

As far as the PE balance is concerned, the numerical counterpart of (59) is

$$ 
\begin{aligned}
& \sum_{\substack{t=t_0,t1 \\ i=1,M \\ k=1,R-1}} g z_{i,k} \left( \frac{F_{i+1/2,k}^x - F_{i-1/2,k}^x}{\Delta x} + \frac{F_{i,k+1/2}^z - F_{i,k-1/2}^z}{\Delta \sigma_k} \right) \Delta \sigma_k \Delta x \Delta t \\
& = - \sum_{\substack{t=t_0,t1 \\ i=1,M-1 \\ k=1,R-1}} g (\tilde{K}_{Ti+1/2,k}^x + \tilde{A}_{i+1/2,k}^x) (\rho_{i+1,k}^{t-\Delta t} - \rho_{i,k}^{t-\Delta t}) \frac{z_{i+1,k} - z_{i,k}}{\Delta x} \Delta \sigma_k \Delta t \\
& \quad - \sum_{\substack{t=t_0,t1 \\ i=1,M \\ k=1,R-2}} g (K_{Ti,k+1/2}^z + A_{i,k+1/2}^z) (\rho_{i,k+1}^{t+\Delta t} - \rho_{i,k}^{t+\Delta t}) \Delta x \Delta t
\end{aligned} 
\quad (60) $$

with $F_{i+1/2,k}^x = (\tilde{K}_{Ti+1/2,k}^x + \tilde{A}_{i+1/2,k}^x) \frac{\rho_{i+1,k}^{t-\Delta t} - \rho_{i,k}^{t-\Delta t}}{\Delta x}$ and $F_{i,k+1/2}^z = (K_{Ti,k+1/2}^z + A_{i,k+1/2}^z) \frac{\rho_{i,k+1}^{t+\Delta t} - \rho_{i,k}^{t+\Delta t}}{\Delta x}$ where $(\tilde{K}_T^x, K_T^z)$ are mixing coefficients given by physical arguments on turbulence, and $(\tilde{A}^x, A^z)$ those deduced from considerations on numerical stability properties of the advection scheme. As for the momentum equations, we note that the horizontal density-mixing scheme is centred on the previous time step, $t - \Delta t$, while the vertical density-mixing is centred on the following time step, i.e. $t + \Delta t$, for numerical stability reasons. We also note that the first term at the RHS of (60) vanishes when the slope of iso-sigma levels is zero. Finally, we retain that the RHS of (60) can be possibly used to build a numerical scheme for the buoyancy term of the TKE equation, in order to retrieve an expected energy-conserving property. Namely, that through the buoyancy term the TKE tank loses (gains) an amount of energy exactly equivalent to what was gained (lost) by the global PE through density mixing processes. This includes vertical mixing occurring indirectly because of the sigma level slope, and possibly additional diffusion required to stabilize the advection scheme. In the same range of ideas, we note that expression (39) can suggest an energy-conserving scheme for the production term of a TKE equation. The interest of using a conservative form of the TKE equation is presented in Burchard (2002).

## 5. Time splitting: an energy-conserving scheme

It is clear that energy exchanges involving frozen terms of advection and pressure anomaly gradient in external mode equations are not exactly equilibrated, because of the forward-like nature of the coupling of external and internal modes. Here, an alternative scheme based on centred-like way of coupling the two modes is proposed. We aim to use the frozen terms in a centred manner; in practice, the frozen terms computed by the internal mode at time step $t$ should be involved in the two external mode sequences bounded by time step $t - \Delta t$ and $t + \Delta t$. Because of this constraint two external modes are used (hereafter: twin external modes). The coupling mechanism of the twin external modes with the internal mode is shown by Fig. 3. This figure can be explained as follows:
<figure><img src="image_4.png" /><figcaption>Fig. 3. Sketch of the twin external modes method.</figcaption></figure>

Step 1: At time step $t$, the internal mode solution allows computation of the frozen terms for external mode 1:

$$ \int_{\sigma_0}^{\sigma_1} \left( -\frac{1}{\rho_0} \frac{\partial p'}{\partial x^*} + \frac{g \rho'}{\rho_0} \frac{\partial z}{\partial x^*} - \frac{\partial \tilde{u}' u'}{\partial x^*} \right) d\sigma \quad (61) $$

Step 2: External mode 1 allows integration of SLA and transport equations from internal time step $t - \Delta t$ to $t + \Delta t$.

Step 3: The transport computed by external mode 1 is time averaged over the period between $t - \Delta t$ and $t + \Delta t$. This enables to update the depth-averaged component of the current in the internal mode equations, similar to the time-splitting technique of BM87.

Step 4: The internal mode is integrated from time step $t - \Delta t$ to $t + \Delta t$.

Step 5: At time step $t + \Delta t$, the internal mode solution allows for computation of the frozen terms for external mode 2.

Step 6: External mode 2 allows for integration of SLA and transport equations from internal time step $t$ to $t + 2\Delta t$.

Step 7: The transport computed by external mode 2 is time averaged over the period between $t$ and $t + 2\Delta t$. This enables to update the depth-averaged component of the current in the internal mode equations.

Step 8: The internal mode is integrated from time step $t$ to $t + 2\Delta t$.

Steps 9–16 are the same as steps 1–8, but for the two following internal time steps.

The barotropic kinetic and PE balances are now computed using the average of solutions of external modes 1 and 2. On the other hand, the energy balance related to the internal mode equations is computed in the same fashion as before. Thus, barotropic KE related to the pressure anomaly gradient is numerically equivalent to expression (49), except that now the barotropic current is averaged from time step $t - \Delta t$ to $t + \Delta t$. As a consequence, the RHS of expression (50) can be rewritten as

$$ -\sum_{\substack{t=t_0, t_1 \\ i=1, M \\ k=1, R-1}} p_{i,k}^{t} \frac{\eta_{i}^{t+\Delta t} - \eta_{i}^{t-\Delta t}}{2} \Delta\sigma_{k} \Delta x \quad (62) $$

which, taking (55) into account, is now numerically equivalent to the second term at the RHS of (52). However, the barotropic component of the current appearing in (49) is exactly the same as the one involved in the density advection, so that the first term at the RHS of (57) now perfectly equilibrates the second term at the RHS of (49). In the same vein, concerning KE exchanges between external and internal modes through advection terms, the twin external modes method causes the barotropic current appearing in expressions (25) and (26) to be numerically the same, so ensuring a perfect equilibrium.

Among the possible shortcomings of the twin external mode method, one might fear the separation of the two external mode solutions, through a process similar to the separation of the solutions at even and odd
iterations that might be induced by leapfrog time-stepping procedures. Similarly to the Asselin filter, intended
to bring even and odd trajectories closer, one could imagine bridging an eventual gap between the two external
modes by a filtering procedure. However, this appears unnecessary, at least in the frame of the numerical
experiments presented in the following section. Indeed, the twin external modes seem linked closely enough,
driven along a common trajectory by the internal mode, so that finally no filtering procedure was used. Also
note that the method of the twin external modes technique still ensures global conservation properties for pri-
mary variables such as momentum, heat or salt, as confirmed by the numerical experiments presented
promptly.

**6. Numerical experiments and discussion**

The results of the previous section are now discussed by means of a set of numerical experiments. Our
reference experiment comprises the generation of internal gravity waves by surface gravity (external) waves
over a submarine ridge. The aim is to construct a complete energy balance, and to quantify the shortcomings
identified in the previous sections. First of all, the energy conservation errors related to the leapfrog time-
stepping procedure in the momentum equations are considered, since we have shown that the SLA time
variation term emerging from the momentum time variation term (a particularity of the sigma-coordinate
formalism) is not *exactly* equilibrated by the term issued from advection terms as it should however be.
Let us recall that this shortcoming is evidenced by the inconsistency between, on one hand, expressions
(17) and (20), and expressions (18) and (21) on the other hand. Secondly, we examine the other shortcoming
of the momentum advection scheme, appearing when the internal to external mode coupling method leads to
inconsistent values of the barotropic current in expressions (25) and (26), thereby preventing equilibrium KE
exchanges between the barotropic and baroclinic modes. Subsequently, we consider energy transfers that play
an important role in internal wave processes. That is, the conversion of mechanical barotropic energy into
baroclinic PE is considered, which involves on one hand density advection, and the frozen term of pressure
anomaly gradient in the external mode on the other hand. This latter process is of particular interest, since,
firstly, we may expect it to be quantitatively significant, and secondly, we know that this conversion mech-
anism is not numerically exact under the classic forward-like way of coupling the external and internal
modes. Two model configurations are then tested. The reference version uses the classic time-splitting tech-
nique, based on the forward-like coupling of the two modes suggested by BM87. In other words, the frozen
terms, provided by the internal mode solution at time step *t*, force the external mode equations between time
steps *t* and *t* + *Δt*. We compare this version of the model to that based on the twin external modes method
presented in the previous section, illustrating that it ensures an exact balance of the energy exchanges with the
internal mode. Furthermore, the contribution of the various dissipation terms to the energy balance is
considered, using standard values for the Asselin and horizontal mixing coefficients, while the vertical mixing
coefficient is computed using the turbulence closure scheme proposed by Gaspar et al. (1990). Finally, the
density mixing processes are addressed.

The reference simulation is described as follows. As in the previous sections, a two-dimensional $Oxz$ numerical domain with closed boundaries is used. The bathymetry is constant ($h = 5000$ m) except in the central part of the domain where we have introduced a submarine ridge, 1000 m high and about 60 km wide. The bathymetry is defined as $h(x) = 5000 - 1000e^{\left(-\frac{x-100}{15}\right)^2}$ where $x-100$ corresponds to the distance, in km, from the middle of the domain. The total length of the domain is 4800 km. At the initial state, a constant surface pressure gradient is prescribed, such that the SLA gap between the two closed boundaries is equal to one meter. The initial density field is horizontally homogeneous, with a constant vertical gradient, $2 \times 10^{-4} \text{ kg m}^{-4}$, while the reference density equals $\rho_0 = 1000 \text{ kg m}^{-3}$. Initial velocities are zero. The geometric features of the domain will lead to the generation of a barotropic seiche at a 12-h period. The seiche will lose a part of its mechanical energy in the production of internal waves over the ridge. Conclusions of these experiments will be helpful when, in future studies, we will apply the model to quantify the energy lost by M2 tidal waves in the generation of internal tides.

The external time step, deduced from a numerical stability criterion considering the grid resolution and the
theoretical surface wave celerity (BM87), is taken to be 4.31 s. In the reference simulation, the internal time
step is 50 times larger than the external one, about 215 s. The left column of Table 1 summarises the various
<table><caption>Table 1
Model parameters of the reference simulation and sensitivity tests</caption><thead><tr><td>Model parameters</td><td>Reference simulation</td><td>Sensitivity tests</td></tr></thead><tbody><tr><td>Time-splitting</td><td>Single external mode</td><td>Twin external modes</td></tr><tr><td>Internal/external time step</td><td>50</td><td>10, 100</td></tr><tr><td>Asselin coefficient (momentum, tracer)</td><td>(0.1, 0)</td><td></td></tr><tr><td>Horizontal mix. coef. (momentum, tracer)</td><td>(15, 0) m<sup>2</sup> s<sup>-1</sup></td><td></td></tr><tr><td>Turbulence closure</td><td>Gaspar et al. (1990)</td><td></td></tr><tr><td>Period of the barotropic seiche</td><td>12 h</td><td>6 h, 24 h</td></tr></tbody></table>

parameters of the reference simulation. The right column shows the parameter modifications considered in order to perform a set of sensitivity tests.

In Table 1, for the sake of clarity, the classical expression “Asselin filter” appears, while in fact the modified version of the Asselin filter, i.e., the time diffusion scheme (35), proposed in Section 3.2 was used. Also note that the reference simulation does not include any horizontal mixing processes in the density equation.

The simulation quickly exhibits a regular barotropic oscillation at a 12-h period. Internal waves are excited above the ridge at the same period. We now briefly examine their characteristics, using a vector projection of the model solution onto the wave propagation modes analytically predicted for the case of a constant Brunt-Väisälä frequency. As far as the vertical velocity is concerned, these propagating modes have a sinusoidal shape (Gill, 1982, [p. 155]). For coherence with the remaining of the paper, here the horizontal (rather than the vertical) component of the velocity is analysed, considering the KE related to each propagation mode separately. Thus, as far as the horizontal velocity is concerned, the z dependence of the modes is of the type $\cos(m\pi z/h)$ where m is an integer number. Therefore, the horizontal velocity is converted into the cosinusoidal series $A_m(x, t) \cos(m\pi z/h)$. The amplitudes, $A_m(x, t)$, of the seven first modes are shown in Fig. 4. The wavelength of the first ($m=1$) baroclinic propagation mode of the computed internal waves is about 95 km. Fig. 4a shows that the corresponding phase speed is about 2.2 m/s, in good agreement with the analytical value (Gill, 1982), obtained in the idealistic case of waves propagating in a linearly stratified medium with a Brunt-Väisälä frequency, $(g\rho_0^{-1}d\rho/dz)^{0.5}$, equal to $1.4 \times 10^{-3} s^{-1}$. Fig. 4b–g also shows that higher ($m>1$) modes verify the expected relationship $C_m = C_1/m$. Nevertheless, the phase speed seems to be underestimated when m is greater than 8. For instance, at the beginning of the simulation, the 10th mode (not shown) seems to propagate with a speed of 0.18 m/s instead of the expected analytical value of 0.22 m/s. It should be noted, considering the 12-h period, that the associated wavelength should be close to 9 km, while such small patterns cannot be correctly represented by a grid with a 3 km horizontal resolution. Nonetheless, Mode 7, with a related wavelength roughly four times the horizontal grid mesh, is correctly modelled (additionally, note that a grid with 30 vertical levels can accurately represent the vertical structure of Mode 7).

Fig. 5 shows the evolution of the global barotropic mechanical energy over the first 12 days. This period corresponds to the time required for the fastest baroclinic mode to reach the closed boundaries, roughly 2400 km from the ridge where internal waves are excited. Since absolute magnitudes have little interest in an academic two-dimensional model, curves of energy balance presented in the following figures have been normalized by the global initial mechanical energy, $E_0$. The energy balance curves each represent a fraction of the barotropic PE available at the initial state. Fig. 5 shows that $0.12 \times E_0$ (or 12% of the global initial barotropic energy) has been dissipated and/or converted into baroclinic energy after 12 days, and that in the same time the baroclinic energy has increased by $0.045 \times E_0$. At this stage, no baroclinic PE curves are presented, because we have not yet introduced the distinction between the available PE and irreversible transformations of PE related to mixing processes. This point is discussed later on.

As mentioned previously, we will first examine the inconsistency between, on one hand, the SLA time variation term emerging from the time variation term of momentum, and its counterpart emerging from advection terms on the other hand, as evidenced by expressions (17), (18), (20) and (21). Using the latter, the time evolution of the corresponding errors for the barotropic and baroclinic balances was computed, which are respectively represented by the quantities $\epsilon_e$ and $\epsilon_t$, defined as follows:
<figure><img src="image_5.png" /><figcaption>Fig. 4. Amplitude (in m/s) of the baroclinic propagation modes as a function of time (in days) and distance (in km) from the ridge. (a) mode 1, (b) mode 2, (c) mode 3, (d) mode 4, (e) mode 5, (f) mode 6, (g) mode 7.</figcaption></figure>

$$ \epsilon_e = \sum_{\substack{i=1,M-1 \\ t=t_0,t_1}} \left( \frac{\hat{u}^t \hat{u}^{t+\Delta t_e}}{2} \frac{\overline{\eta^{t+\Delta t_e}}^x - \overline{\eta^t}^x}{\Delta t_e} - \frac{\hat{u}^{t2}}{2} \frac{\overline{\eta^{t+\Delta t_e}}^x - \overline{\eta^{t-\Delta t_e}}^x}{2\Delta t_e} \right)_{i+1/2} \Delta x \Delta t_e \quad (63) $$

$$ \epsilon_i = \sum_{\substack{i=1,M-1 \\ k=1,R-1 \\ t=t_0,t_1}} \left( \frac{u^t u^{t+\Delta t}}{2} \frac{\overline{\eta^{t+\Delta t}}^x - \overline{\eta^t}^x}{\Delta t} - \frac{u^{t2}}{2} \frac{\overline{\eta^{t+\Delta t}}^x - \overline{\eta^{t-\Delta t}}^x}{2\Delta t} \right)_{i+1/2,k} \Delta x \Delta \sigma_k \Delta t \quad (64) $$

The errors are observed to be rather small (figure not shown), since $\epsilon_e$ and $\epsilon_i$ are respectively of the order of $10^{-12} \times E_0$ and $10^{-10} \times E_0$ after 12 days of simulation, which is several orders of magnitude smaller than the energy variations at play in the current simulation. Thus, these shortcomings have no influence on the accuracy of the global energy balance, at least in the context of the current numerical experiment.

Now, let us examine the equilibrium of exchanges of KE between external and internal modes through momentum advection terms. In a closed domain, due to this equilibrium expression (23) should vanish. The numerical counterpart of expression (23) was shown to be obtained by adding expressions (25), (26) and (28). In addition, it was shown that, because of the forward-like nature of the coupling of the external and internal modes, the inconsistency between values of the barotropic current appearing in expression (25)
<figure><img src="image_6.png" /><figcaption>Fig. 4 (continued)</figcaption></figure>

and (26) is likely to prevent $\epsilon = (25) + (26) + (28)$ from vanishing. Nevertheless, this error remains negligible (figure not shown), since after 12 days simulation we have $\epsilon \approx 10^{-8} \times E_0$. Therefore, it does not change the global energy balance, at least, once more, as far as this particular experiment is concerned.

The third shortcoming identified in the previous sections concerns the exchanges between the barotropic kinetic and potential forms of energy through, on one hand, the frozen term of pressure anomaly gradient of the external mode, and, on the other hand, the density advection by the barotropic component of the velocity. This drawback arises from the forward-like nature of the coupling of the external and internal modes, which causes the numerical expressions (47), (49), (52) and (57) to be numerically inconsistent. This inconsistency prevents the energy balance expressed by (58) to be numerically exact. Indeed, Fig. 6 shows that the two curves respectively representing the time evolution of the RHS and LHS numerical counterpart of (58) are not superimposed. The relative error is defined as

$$ \epsilon = \frac{|A - B|}{0.5 \times (|A| + |B|)} \qquad (65) $$

where A and B are respectively the numerical counterparts of the RHS and LHS of expression (58). After 12 days, the LHS of (58) is equal to $0.0447 \times E_0$. This value is quite significant, since it is comparable to the increase of the baroclinic KE over the same period. The corresponding value at the RHS of (58) is about $0.0456 \times E_0$ (Fig. 6). According to (65), this represents a relative error of 2%. Recall that this error is caused
<figure><img src="image_7.png" /><figcaption>Fig. 5. Global barotropic mechanical energy (solid line) and global kinetic baroclinic energy (dashed line).</figcaption></figure>

by the inconsistency between barotropic current values involved in the numerical counterpart of (58) because of
the forward-like nature of the frozen term of the external mode: the barotropic current in (49) corresponding to
a time average over the period [t, t + Δt] when an exact equilibrium with (52) and (57) would actually require a
time average over two external mode sequences, i.e., [t - Δt, t + Δt]. The error margin of the barotropic current
value, and thus its impact on the imbalance of (65), is likely to be smaller for a signal with a longer period. Alter-
natively, for a given periodicity, one expects the error to be smaller when using a shorter time step. To check
these two properties, sensitivity tests regarding the periodicity of the barotropic signal and the model time step
were performed. First of all, the reference time step (Table 1) is maintained, and periods of 6 h and 24 h are
considered. Note that these new periods are simply obtained by changing the size of the numerical domain.
The relative error, given by (65), is assessed at the end of each simulation. To maintain some coherence with
the previous 12-h period simulation, the end of these two simulations corresponds to the time required for

<figure><img src="image_8.png" /><figcaption>Fig. 6. Numerical counterparts of the LHS (solid line) and RHS (dashed line) terms of expression (62).</figcaption></figure>
the fastest baroclinic waves to reach the closed boundaries of the numerical domain, i.e., respectively 6 and 24 days. In the case of a 6-h periodic signal, the relative error is about 12%, which actually exceeds the imbalance of the reference simulation. Conversely, in the case of the 24-h periodic signal, the relative error is about 1%, which is smaller than for the imbalance of the reference simulation, as expected. The second test focuses on the internal time step, considering a given periodicity of the barotropic signal. The internal time step of the reference simulation corresponds to 50 external time steps. As far as the 12-h periodic signal is concerned, the relative error is smaller (0.4%) with a shorter time step (equivalent to 10 external time steps) and greater (5%) with a longer internal time step (equivalent to 100 external time steps), as expected. Table 2 summarises the results obtained for all tests performed regarding the periodicity of the barotropic signal and the internal time step. Clearly, in the least favourable case (longer time step and shorter period) the relative error can reach 30%. Conversely, the equilibrium shortcoming is negligible in the most favourable case.

The imbalance of the energy transfer between terms involving density advection and the pressure anomaly gradient has undesirable numerical repercussions on the amount of energy removed from the barotropic energy tank through the internal wave generation process. Among them, we note a strong dependency on the value of the internal time step. Fig. 7 shows, for the 6-h periodic signal case, the decrease of the barotropic energy with time. After 6 days, which is the time required for waves to reach the boundaries of the domain, the decrease of the global mechanical barotropic energy is about $0.068 \times E_0$ (respectively $0.074 \times E_0$) for an internal time step equal to 100 (respectively 10) external time steps. Alternatively, the same simulations performed with the conserving method of the twin external modes appear to be much less sensitive to the value of the internal time step. Indeed, after 6 days, the decrease of barotropic mechanical energy is very close to $0.074 \times E_0$ regardless of the value of internal time step (Fig. 8). Therefore, the decrease of energy is close to that obtained with both the classic time-splitting method and a very small internal time step, suggesting

<table><caption>Table 2
Balance error (%) given by Eq. (61), at the end of simulation, depending on wave-periodicity and internal over external time step ratio</caption><thead><tr><th></th><th>6 h</th><th>12 h</th><th>24 h</th></tr></thead><tbody><tr><td>Int/ext time step: 100 (%)</td><td>30</td><td>5</td><td>2</td></tr><tr><td>Int/ext time step: 50 (%)</td><td>12</td><td>2</td><td>1</td></tr><tr><td>Int/ext time step: 10 (%)</td><td>2</td><td>0.4</td><td>0.2</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Global barotropic mechanical energy using an internal/external time step ratio equal to 10 (solid line), 50 (dashed line), 100 (dotted line).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Same as Fig. 7 but with the twin external mode method.</figcaption></figure>

that the shortcoming of the classic technique in terms of energy conservation is to underestimate the loss of
barotropic energy for large internal time steps.

7. Conclusions

The numerical schemes of an ocean circulation model are considered in details in this work, from the point of view of the global conservation of energy. An overview of the kinetic and PE transfers is shown in Fig. 9: solid and dashed arrows recapitulate respectively the exact and approximate numerical balances associated to these transfers together with the section numbers where they are discussed.

First, because of the free-surface method, the balance between the advection and tendency terms in the momentum equation is not achieved anymore. However, numerical experiments indicate that this shortcoming is negligible. Second, it is shown that some energy conservation properties are lost due to the time-splitting method used to limit computational costs associated to the resolution of surface waves. This arises when the coupling of the external and internal modes is forward in time. In other words, when the coupling terms, provided at time step *t* by the internal mode solution, force the external mode equations over a period of time bounded by time steps *t* and *t* + *Δ*t. Numerical tests show that the lack of energy conservation regarding the gradient of pressure anomaly can be important, as non-linear terms appear to be negligible in the framework of the considered test case. This drawback is shown to be more significant for large internal time steps, or alternatively, when the internal wave periods are short. The principal effect is that the amount of mechanical energy removed from the external mode by the generation of internal waves can be biased or at least dependent on the internal mode time step. An energy-conserving alternative consists of using a leapfrog type coupling of the external and internal mode, where the frozen terms, which force the external mode equations over the period of time bounded by internal time steps *t* − *Δ*t and *t* + *Δ*t, are provided at time step *t* by the internal mode equations. The main consequence is that the dependency of the barotropic mechanical energy level on the model internal time step is reduced. This conservative time-splitting technique is based on two sets of twin external mode equations. Numerical tests show that the two barotropic solutions remain close to each other, since the single internal mode that the two external modes have in common binds the two barotropic solutions to the same trajectory. This, however, may cease to be the case in realistic three-dimensional simulations; future research might involve defining a way to ensure this property. Also note, that the twin external mode technique has the inconvenience of doubling the computational cost of the barotropic mode.

Moreover, the study of the energy balance sheds some light on numerical options concerning the mixing
processes of the model. Considering the Asselin filter, which can be approximated as a classic mixing
<figure><img src="image_11.png" /><figcaption>Fig. 9. Overview of the kinetic and potential energy transfers in a free-surface s coordinate model. Red arrows: satisfied numerical balances, red dashed arrows: approximated balances due to leapfrog scheme, blue dashed arrows: approximated balances due to time splitting, black arrows: transfers toward TKE, external and internal energy tanks. Note that “external” and “internal” energies refer here to the fundamental quantities defined in Winters et al. (1995). Numbers in parentheses indicate the section number where the related balance is discussed. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

second-order operator, or considering the horizontal diffusion of velocities, the energy balance is a way to explain why computing mixing processes using the variables provided at the previous time step improve the damping of high-frequency numerical modes.

Not all details of the role of tracer mixing in the energy balance were considered, because, although a full-fledged and non-trivial problem, the question of background PE is beyond the scope of this paper. Nevertheless, the energy downscaling problem raised by the turbulence closure is addressed. Insomuch as the increase of PE induced by density mixing should be associated to an equivalent decrease term in the TKE equation, the numerical analysis of the energy balance calls for an appropriate (i.e., conservative) form of the buoyancy term of the turbulence closure scheme, as recommended by Burchard (2002). In a sigma-coordinate model, horizontal tracer mixing operators, explicitly formulated or possibly hidden in an upstream type advection scheme, actually contain a part of vertical mixing which is, a priori, not taken into account in the TKE equation. However, according to preliminary tests (not shown), it seems illusory to think that the TKE equation can be adapted, because the amount of energy involved in these unexpected vertical mixing processes likely exceeds the amount of available TKE.

The conclusions of this work form the foundation of a model well adapted to the study of the energy balance of internal waves. The estimate of the energy transferred from barotropic tides to internal waves in the Bay of Biscay (North East Atlantic) will be the subject of future studies.
Acknowledgements

This study was funded by the European MFSTEP Project (EU Contract EVK3-CT-2002-00075), the European INSEA Project (SST4-CT-2005-012336), the “LEFE-IDAO-2006” and “OPA Symphonique” french projects. The authors thank Cyril Nguyen and the Laboratoire d’Aérologie computer team, Serge Prieur, Laurent Cabanas, Jérémy Leclercq, Didier Gazen, and Juan Escobar for their support. We thank the anonymous reviewers for their thoughtful and constructive comments.

Appendix A

The numerical rule of integration by parts is given by

$$ \sum_{n=N_0, N_1} A_{n+1/2} (B_{n+1} - B_n) = - \sum_{n=N_0+1, N_1} B_n (A_{n+1/2} - A_{n-1/2}) + B_{N_1+1} A_{N_1+1/2} - B_{N_0} A_{N_0+1/2}. \quad (A.1) $$

When $B_n$ is a momentum advective flux of the form $B_n = C_n(A_{n+1/2} + A_{n-1/2})$, where $A_{n+1/2}$ is for instance a velocity component with closed boundary conditions ($A_{N_0+1/2} = A_{N_1+1/2} = 0$), (A.1) becomes

$$ \sum_{n=N_0, N_1} A_{n+1/2} (B_{n+1} - B_n) = \sum_{n=N_0, N_1} A_{n+1/2}^2 (C_{n+1} - C_n) \quad (A.2) $$

References

Arakawa, A., Lamb, V.R., 1977. Computational design of the basic dynamical processes of the UCLA general circulation model. Methods in Computational Physics 17, 174–267.

Arakawa, A., Suarez, M.J., 1983. Vertical differencing of the primitive equations in sigma coordinates. Monthly Weather Review 111, 34–45.

Asselin, R., 1972. Frequency filter for time integrations. Monthly Weather Review 100, 487–490.

Auclair, F., Casitas, S., Marsaleix, P., 2000. Application of an inverse method to coastal modelling. Journal of Atmospheric and Oceanic Technology 17, 1368–1391.

Beckmann, A., Haidvogel, D.B., 1993. Numerical simulation of flow around a tall isolated seamount. Part I: problem formulation and model accuracy. Journal of Physical Oceanography 23, 1736–1753.

Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal circulation model. In: Heaps, N. (Ed.), Three-Dimensional Coastal Ocean Models, In: Coastal Estuarine Science, vol. 4. American Geophysical Union, pp. 1–16.

Burchard, H., 2002. Energy-conserving discretisation of turbulent shear and buoyancy production. Ocean Modelling 4, 347–361.

Deleersnijder, E., 1993. Numerical mass conservation in a free-surface coordinate marine model with mode splitting. Journal of Marine Systems 4, 365–370.

Ferziger, J.H., Peric, M., 2002. Computational methods for fluid dynamics, third ed. Springer, 423 pp.

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. Journal of Geophysical Research 95, 16,179–16,193.

Gerdes, R., 1993. A primitive equation ocean circulation model using a general vertical coordinate transformation. 1. Description and testing the model. Journal of Geophysical Research 98, 14,683–14,701.

Gill, A.E., 1982. Atmosphere Ocean Dynamics. Academic Press, 662 pp.

Griffies, S.M., Pacanowski, R.C., Schmidt, M., Balaji, V., 2001. Tracer conservation with an explicit free surface method for z-coordinate ocean models. Monthly Weather Review 5, 1081–1098.

Haltiner, G.J., Williams, R.T., 1980. Numerical prediction and dynamic meteorology, second ed. John Wiley & Sons, 477 pp.

Haney, R.L., 1991. On the pressure gradient force over steep topography in sigma coordinate ocean models. Journal of Physical Oceanography 21, 610–619.

Higdon, R.L., 2005. A two-level time-stepping method for layered ocean circulation models: further development and testing. Journal of Computational Physics 206, 463–504.

James, I.D., 1996. Advection schemes for shelf sea models. Journal of Marine Systems 8, 237–254.

Janjic, Z.I., 1977. Pressure gradient force and advection scheme used for forecasting with steep and small scale topography. Contributions to atmospheric physics 50, 186–199.

Johns, B., Oguz, T., 1987. Turbulent energy closure schemes. In: Heaps, N. (Ed.), Three-Dimensional Coastal Ocean Models, In: Coastal Estuarine Science, vol. 4. American Geophysical Union, pp. 17–39.

Johns, B., Sinha, P.C., Dube, S.K., Mohanty, U.C., Rao, A.D., 1983. Simulation of storm surges using a three-dimensional numerical model: an application to the 1977 Andhra Cyclone. Quarterly Journal of the Royal Meteorological Society 109, 211–224.

Lilly, D.K., 1965. On the computational stability of numerical solutions of time-dependent non-linear geophysical fluid dynamics problems. Monthly Weather Review 93, 11–26.
Madec, G., Delecluse, P., Crepon, M., Chartier, M., 1991. A three-dimensional numerical study of deep-water formation in the North-western Mediterranean Sea. Journal of Physical Oceanography 21, 1349-1371.

Marsaleix, P., Auclair, F., Estournel, C., 2006. Considerations on open boundary conditions for regional and coastal ocean models. Journal of Atmospheric and Oceanic Technology 23, 1604-1613.

Mellor, G.L., 2003. Users guide for a three-dimensional, primitive equation, numerical ocean model. Program in Atmospheric and Oceanic Sciences. Princeton University, 53 pp.

Mellor, G.L., Yamada, T.Y., 1982. Development of a turbulent closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics 20, 851-875.

Reffray, G., Fraunié, P., Marsaleix, P., 2004. Secondary flows induced by wind forcing in the Rhone region of freshwater influence. Ocean Dynamics 54, 179-196.

Sadourny, R., 1975. The dynamics of finite-difference models of the shallow-water equations. Journal of Atmospheric Sciences 32, 680-689.

Shechepetkin, A.F., McWilliams, J.C., 2003. A method for computing horizontal pressure-gradient force in an oceanic model with non-aligned vertical coordinate. Journal of Geophysical Research 108 (3090), 35.1-35.34.

Shulman, I., Lewis, J.K., Blumberg, A.F., Nicholas Kim, B., 1998. Optimized boundary conditions and data assimilation with application to the M2 tide in the Yellow Sea. Journal of Atmospheric and Oceanic Technology 15, 1066-1071.

Smagorinsky, J., 1993. Some historical remarks on the use of nonlinear viscosities. In: Galperin, B., Orszag, S. (Eds.), Large Eddy Simulations of Complex Engineering and Geophysical Flows. Cambridge University Press, Cambridge, UK, pp. 4-36.

Song, Y.T., 1998. A general pressure gradient formulation for ocean models. Part I: scheme design and diagnostic analysis. Monthly Weather Review 126, 3213-3230.

Song, Y.T., Wright, D.G., 1998. A general pressure gradient formulation for ocean models. Part II: energy, momentum, and bottom torque consistency. Monthly Weather Review 126, 3231-3247.

Winters, K.B., Lombard, P.N., Riley, J.J., D'Asaro, E., 1995. Available potential energy and mixing in density-stratified fluids. Journal of Fluid Mechanics 289, 115-128.