ABSTRACT

We present a set of equations describing the nonlinear dynamics of flows constrained by environmental rotation and stratification (Rossby numbers $Ro\in \left[\mathrm{0.1,0.5}\right]$ and Burger numbers of order unity). The fluid is assumed incompressible, adiabatic, inviscid and in hydrostatic balance. This set of equations is derived from the Navier Stokes equations (with the above properties), using a Rossby number expansion with second order truncation. The resulting model has the following properties: 1) it can represent motions with moderate Rossby numbers and a Burger number of order unity; 2) it filters inertia-gravity waves by assuming that the divergence of horizontal velocity remains small; 3) it is written in terms of a single function of space and time (pressure, generalized streamfunction or Bernoulli function); 4) it conserves total (Ertel) vorticity in a Lagrangian form, and its quadratic norm (potential enstrophy) at the model order in Rossby number; 5) it also conserves total energy at the same order if the work of pressure forces vanishes when integrated over the fluid domain. The layerwise version of the model is finally presented, written in terms of pressure. Integral properties (energy, enstrophy) are conserved by these layerwise equations. The model equations agree with the generalized geostrophy equations in the appropriate parameter regime. Application to vortex dynamics are mentioned.

Keywords:

Incompressible and Stratified Fluids, Intermediate Model, Rossby Number, Burger Number, Conservation Properties, Energy, Potential Vorticity

1. Introduction

Flows of incompressible fluids, strongly constrained in the horizontal plane by stratification and by background rotation, are often found in geophysical environments and in laboratory experiments (see [1] and references therein). In particular, motions at mesoscale and large scale in the Earth oceans and atmosphere are mostly horizontal due to the combined action of the Coriolis and buoyancy forces. Under these conditions, intense flow structures, such as jets and vortices, can be long-lived and ensure a significant part of the mass, momentum and heat transfers within the fluid (see [2] [3] ; and their references). In the ocean, the intense western boundary currents (e.g. the Gulf Stream), but also the numerous and less intense eastern boundary currents (e.g. the Mediterranean Water Undercurrents along the Iberian continental slope), are unstable at various locations, and shed such long-lived vortices (the Gulf-Stream Cold-Core and Warm-Core rings, or the Meddies, Mediterranean Water Eddies). Such eddies have therefore small, but finite Rossby numbers ( $Ro\in \left[\mathrm{0.1,0.5}\right]$ ) and Burger numbers most often on the order of unity ( $Bu\in \left[\mathrm{0.5,5.0}\right]$ ).

Modeling these flows has therefore been the subject of many studies, using the “primitive equations”, i.e. the three-dimensional Navier-Stokes equations with hydrostatic balance, Boussinesq approximation and incompressibility. These primitive equations, describe both fast (wave-like) and slow (vortex-like) motions. These two subspaces of solutions of the primitive equations are independent only for linear dynamics. In nonlinear dynamics, fast and slow motions can interact : initially slow motions, free of fast waves, remain in the slow subspace, but nonlinear interaction of fast waves can contribute to slow motions. Therefore the slow motions subspace is called a quasi manifold.

The importance of quasi-2D jets and vortices in geophysical flows and laboratory experiments has motivated the search for asymptotic limits of the primitive equations on the slow quasimanifold (see for instance [4] ). This regime is governed by the time evolution of potential vorticity (also called Ertel vorticity, a combination of vertical vorticity with vertical density gradient). For incompressible and inviscid fluids, potential vorticity is conserved by each fluid element (see [5] ). To filter out fast (inertia-gravity) waves, zero or weak divergence of horizontal velocity is usually imposed. The most prominent example of slow regime is the horizontal balance between the Coriolis acceleration and the pressure gradient (for parallel flows). This equilibrium is known as geostrophic balance, for which horizontal velocity has zero divergence. The following model in complexity, with small horizontal velocity divergence, is the quasi-geostrophic model [6] . This model is derived from the primitive equations, by assuming dominant Coriolis acceleration compared with relative acceleration (i.e. small Rossby number Ro) and comparable influence of buoyancy and Coriolis effects (i.e. Burger number Bu of order unity). The quasi-geostrophic model conserves energy and potential vorticity up to ( $R{o}^2$ ).

At the following order in Rossby number, expansions of the primitive equations yield so-called “intermediate models” which fall into two general classes: the balance equations with horizontal velocity splitting into rotational and divergent parts ( [7] [8] [9] ) and the generalized geostrophic equations with expansion of horizontal velocity in horizontal derivatives of pressure ( [10] [11] ).

The first model has a wide range of applicability (in Rossby and Burger numbers) but relies numerically on a complex and time-consuming iterative method [12] . The second model is much easier to implement but presents an excessive bias between cyclones and anticyclones ( [13] [14] ).

The purpose of the present work is to propose a set of equations, derived along the same principle as the generalized geostrophic equations, but assuring the conservation of energy and of potential (Ertel) vorticity (and therefore of its L²-norm, total enstrophy) in the derivation. The derivation is also carried out directly from the primitive equations for continuously stratified flows, contrary to the generalized geostrophic model, derived from the shallow-water equations (layerwise primitive equations). Then, the projection of these equations onto a set of homogeneous layers is shown to be in agreement with equations directly derived from the shallow-water equations. Such a model is important for the study of mesoscale vortex dynamics in the ocean; for these vortices, Ro is not very small but often remains smaller than 1. When modeling the stability or the interaction of such vortices, it is essential that the model preserves their potential vorticity, and their energy when both are physically conserved ( [15] [16] [17] ).

Here, we present the primitive equations and the physical scaling (Section 2). In Section 3, the intermediate model is derived and its conservation properties at $O\left(R{o}^3\right)$ are presented: Lagrangian invariance of Ertel vorticity, energy conservation in any volume for which the integral flux of the Bernoulli function across the boundary vanishes (as for primitive equations). The model equations are written in terms of a generalized streamfunction or of Bernoulli function. In Section 4, the layerwise model is presented and compared with previously derived models (frontal geostrophic equations, generalized geostrophic equations). Finally, conclusions are drawn.

2. Mathematical Framework

2.1. Dimensional Primitive Equations

The fundamental equations for fluid motions strongly constrained by rotation and stratification are the primitive equations, derived from the three-dimensional Navier-Stokes equations with background rotation, under the following assumptions:

・ incompressibility (moderate particle velocities with respect to the speed of sound);

・ Boussinesq approximation (invariance of density in the horizontal momentum equations);

・ negligible vertical accelerations leading to hydrostatic balance;

・ moderate vertical velocities, implying that the Coriolis force mostly acts on the horizontal velocity;

・ background rotation a linear function of the y coordinate $f=f_0+\beta y$ (beta-plane approximation). On Earth, this corresponds to a projection of the motion on the tangential plane to Earth at a given latitude: $f=2\Omega \sin \left(\theta \right)=2\Omega \left[\sin \left({\theta }_0\right)+\left(\theta -{\theta }_0\right)\cos \left({\theta }_0\right)\right]$ .

Under these hypotheses, assuming adiabatic and inviscid evolutions, the primitive equations are written in dimensional form ( [5] ):

$Q{M}_x\mathrm{<mo>\; :\; </mo>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\left({\partial }_t+u^{\mathrm{<mo>\; *\; </mo>}}{\partial }_x+v^{\mathrm{<mo>\; *\; </mo>}}{\partial }_y+w^{\mathrm{<mo>\; *\; </mo>}}{\partial }_z\right)u^{\mathrm{<mo>\; *\; </mo>}}-f{v}^{\mathrm{<mo>\; *\; </mo>}}=-\frac{1}{{\rho }_s^{\mathrm{<mo>\; *\; </mo>}}}{\partial }_x{p}_T^{\mathrm{<mo>\; *\; </mo>}}$ (1)

$Q{M}_y\mathrm{<mo>\; :\; </mo>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\left({\partial }_t+u^{\mathrm{<mo>\; *\; </mo>}}{\partial }_x+v^{\mathrm{<mo>\; *\; </mo>}}{\partial }_y+w^{\mathrm{<mo>\; *\; </mo>}}{\partial }_z\right)v^{\mathrm{<mo>\; *\; </mo>}}+f{u}^{\mathrm{<mo>\; *\; </mo>}}=-\frac{1}{{\rho }_s^{\mathrm{<mo>\; *\; </mo>}}}{\partial }_y{p}_T^{\mathrm{<mo>\; *\; </mo>}}$ (2)

$Q{M}_z\mathrm{<mo>\; :\; </mo>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\partial }_z{p}_T^{\mathrm{<mo>\; *\; </mo>}}=-{\rho }_T^{\mathrm{<mo>\; *\; </mo>}}g$ (3)

$C\mathrm{<mo>\; :\; </mo>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\partial }_x{u}^{\mathrm{<mo>\; *\; </mo>}}+{\partial }_y{v}^{\mathrm{<mo>\; *\; </mo>}}+{\partial }_z{w}^{\mathrm{<mo>\; *\; </mo>}}=0$ (4)

$H\mathrm{<mo>\; :\; </mo>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\left({\partial }_t+u^{\mathrm{<mo>\; *\; </mo>}}{\partial }_x+v^{\mathrm{<mo>\; *\; </mo>}}{\partial }_y+w^{\mathrm{<mo>\; *\; </mo>}}{\partial }_z\right){\rho }_T^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>\; 0}$ (5)

where $u^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{w}^{\mathrm{<mo>\; *\; </mo>}}$ are the velocity components (functions of $x\mathrm{,}y\mathrm{,}z\mathrm{,}t$ ), g is the local gravity intensity; $p_T^{\mathrm{<mo>\; *\; </mo>}}$ and ${\rho }_T^{\mathrm{<mo>\; *\; </mo>}}$ are the total pressure and density; pressure is split into static and dynamic parts and density into mean and anomaly:

${\rho }_T^{*}={\rho }_s^{*}\left(z\right)\left[1+{\rho }^{*}\left(x,y,z,t\right)\right],\mathrm{<mi>\; </mi>}{p}_T^{*}=p_s^{*}\left(z\right)+{\rho }_s^{*}\left(z\right)p^{*}\left(x,y,z,t\right),$

assuming that the vertical variations dominate the horizontal ones. Viscous dissipation of momentum, forces, and diabatic terms can be added to the model if necessary.

Stars distinguish these dimensional variables (and their relations) from their dimensionless counterparts in Section 2.3 and after. The hydrostatic relation holds for both parts:

${\partial }_z{p}_s^{*}=-{\rho }_s^{*}g,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\partial }_z\left({\rho }_s^{*}{p}^{*}\right)=\mathrm{<mi>\; </mi>}-{\rho }_s^{*}{\rho }^{*}g.$

Note that Equation (H) can also be written

$\frac{\text{d}{\rho }^{*}}{\text{d}t}+w^{*}\left[\frac{1}{{\rho }_s^{*}}{\partial }_z{\rho }_s^{*}\left(1+{\rho }^{*}\right)+{\partial }_z{\rho }^{*}\right]=0$ (6)

with $\text{d}/\text{d}t={\partial }_t+u^{*}{\partial }_x+v^{*}{\partial }_y$ . The original five equations relate the five unknowns $u^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{w}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\rho }^{\mathrm{<mo>\; *\; </mo>}}$ . Note that ${\rho }_s^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{p}_s^{\mathrm{<mo>\; *\; </mo>}}$ are known when the state of the fluid at rest is given.

2.2. Scale Analysis

The horizontal velocity $u^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}$ is now scaled by U, the vertical velocity $w^{\mathrm{<mo>\; *\; </mo>}}$ by W, horizontal lengths $x\mathrm{,}y$ by L, heights by H, time by T, dynamical pressure ( $p^{\mathrm{<mo>\; *\; </mo>}}$ ) by P and mean density ( ${\rho }_s^{\mathrm{<mo>\; *\; </mo>}}$ ) by ${\rho }_0$ .

As density mainly varies along the vertical axis, the Brunt-Vaisala frequency is defined as $N^2=-g{\partial }_z{\rho }_s^{*}/{\rho }_s^{*}$ on average on the horizontal domain; it represents the oscillation frequency of fluid particles around a static equilibrium in a stable configuration. $N_0^2$ is the vertical average of $N^2\left(z\right)$ and $S\left(z\right)=N^2\left(z\right)/N_0^2$ .

From these scales, several dimensionless parameters are defined:

・ the spatial Rossby number $Ro=U/f_{0}L$ ,

・ the temporal Rossby number $R{o}_T=1/f_{0}T$ ,

・ the beta-effect intensity $\beta ={\beta }^{*}L/f_0$ ,

・ the Burger number $Bu={\left(N_{0}H/f_{0}L\right)}^2$ ,

・ the aspect ratio $\delta =H/L$ .

From these numbers and dynamical relations, we scale:

・ pressure via geostrophy, $P={\rho }_{0}fLU$ ;

・ mean density and anomaly from the the pressure scaling and via hydrostatics (QM_z):

${\rho }_T^{*}={\rho }_0{\rho }_s\left(z\right)\left[1+\left(Ro/Bu\right)\left(N_0^{2}H/g\right)\rho \left(x,y,z,t\right)\right];$

・ the vertical velocity from the previous scalings and with equation (H), $W=U\delta Ro/Bu$ .

To illustrate this scaling, consider vortices in the ocean (at mid latitudes), for which $L=50\text{km}$ , $H=1\text{km}$ , $U=1.0\text{m}/\text{s}$ , $T={10}^5\text{s}$ , $f_0={10}^{-4}{\text{s}}^{-1}$ , ${\beta }^{*}=2\times {10}^{-11}{\text{m}}^{-1}\cdot {\text{s}}^{-1}$ , $N_0^2={10}^{-5}{\text{s}}^{-2}$ , leading to $Bu=2.5$ , $Ro=0.2$ , $R{o}_T=0.1$ , $\beta =0.01$ . The following model is designed to study such intense flows (jets or vortices), with $Ro\ge 0.1$ , beyond the quasi-geostrophic approximation. The upper limit of validity of the following model is determined a posteriori from estimates of the first truncated terms in the equations, as well as with numerical tests: the model loses accuracy near $Ro=0.5$ , for $Bu~O\left(1\right)$ . For $Ro~1$ , the primitive equations should be used. Moreover, as $\beta$ is much smaller than $Ro$ , the beta-plane approximation is justified.

2.3. Dimensionless Primitive Equations

With this scaling and dropping stars for dimensional variables, equation (QM_x) is:

$Ro\frac{Du}{Dt}-\left(1+\beta y\right)v=\mathrm{<mi>\; </mi>}-{\partial }_{x}p\mathrm{,}$ (7)

Equation QM_y is:

$Ro\frac{Dv}{Dt}+\left(1+\beta y\right)u=-{\partial }_{y}p\mathrm{,}$ (8)

Equation QM_z is:

${\partial }_{z}p=-\rho ,$ (9)

Equation C is:

${\partial }_{x}u+{\partial }_{y}v+\frac{Ro}{Bu}{\partial }_{z}w=0,$ (10)

and Equation H is

$\frac{D\rho }{Dt}=wS\left(z\right)\mathrm{,}$ (11)

where the total (three-dimensional) derivative is written in scaled form:

$\frac{D}{Dt}=\frac{R{o}_T}{Ro}{\partial }_t+u{\partial }_x+v{\partial }_y+\frac{Ro}{Bu}w{\partial }_z,$ (12)

for which no truncation of $u\mathrm{,}v\mathrm{,}w$ in $Ro$ has yet been made.

Following [5] , it is straightforward to show that these equations conserve total (potential) vorticity

$\frac{D\Pi }{Dt}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>\; 0},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\Pi \mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left(\left(1+\beta y\right)n_z+Ro\nabla \times u\right)\cdot \left(-S\left(z\right)n_z+RoBu\nabla \rho \right)$

where $n_z$ is the unit vertical vector.

3. Derivation of the Model Equations

3.1. Generalized Vorticity Equation

We first rewrite equation H to extract w as an exact derivative:

$\frac{Ro}{Bu}w\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{Ro}{Bu}\frac{D}{Dt}\left[\frac{\rho }{S}-\frac{Ro}{Bu}\frac{{\rho }^2}{2S}{\partial }_z\frac{1}{S}\right]+O\left(R{o}^3\right)=\frac{Ro}{Bu}\frac{DX}{Dt}+O\left(R{o}^3\right).$ (13)

This defines X, a quantity physically related to the static and dynamic densities. This expression of the vertical velocity w is central to the subsequent algebra. We note that this expression is similar in form to the quasi-geostrophic expression of w, for which the Lagrangian derivative is simply two-dimensional:

$w=\text{d}\left(\rho /S\right)/\text{d}t$ .

We also introduce the absolute vorticity Y, the sum of the relative vorticity and of the Coriolis parameter:

$RoY=\beta y+Ro\left({\partial }_{x}v-{\partial }_{y}u\right)\mathrm{.}$ (14)

The generalized vorticity equation (GV), that we derive hereafter, can be obtained by successive steps:

・ first by taking the curl of (QM_x, QM_y),

・ then by substituting $u\mathrm{,}v$ via (QM_x, QM_y), ${\partial }_{z}p$ via (QM_z) and w via (H),

・ and by canceling the three-dimensional velocity divergence via equation (C) (without any approximation).

・ Consistent $O\left(R{o}^3\right)$ eliminations finally lead to the generalized vorticity equation (Equation (15) hereafter).

This rather tedious algebra can be more elegantly replaced by the following symbolic combination leading to (GV):

$\begin{array}{c}\left(GV\right)\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left(\nabla \times \left(QM\right)\right)\cdot n_z+\frac{Ro}{Bu}\left(1+\frac{Ro}{Bu}{\partial }_{z}X\right)\left[\nabla X\times \left({\partial }_z\left(QM\right)-\nabla \left(Q{M}_z\right)\right)\right]\cdot n_z\\ -\left[\left(1+RoY\right)+\frac{R{o}^2}{Bu}\left[\nabla X\times {\partial }_{z}u\right]\cdot n_z\right]\left(C\right)\\ +\frac{Ro}{Bu}\left(1+\frac{Ro}{Bu}{\partial }_{z}X\right)\left[\left[\left(1+RoY\right)+\left(Ro{\partial }_{z}u\right)\times \nabla \left(H\right)\right]\cdot n_z\right],\end{array}$

where $QM=\left(Q{M}_x\mathrm{,}Q{M}_y\mathrm{,}Q{M}_z\right)$ is a symbolic writing of Equations (7)-(9).

Physically, in this expression:

・ the first term on the first line is the absolute vorticity;

・ the rest of the first line represents the coupling of the baroclinicity vector with the vortex stretching;

・ the first term on the second line is the coupling of the velocity divergence with the absolute vorticity;

・ the second part of the second line is the coupling of that velocity divergence with baroclinicity;

・ finally, the third line contains the vortex stretching, its coupling with the absolute vorticity and with the baroclinicity vector.

At first order in Ro, the total vorticity is that of the quasi-geostrophic (QG) model, and at the following order in Ro all terms are nonlinear and represent couplings between the QG terms.

Following the formal calculation, the (GV) equation is obtained:

$\left(1+\frac{Ro}{Bu}{\partial }_{z}X\right)\frac{D}{Dt}\left[\left(1+RoY\right)\left(1-\frac{Ro}{Bu}{\partial }_{z}X\right)+\frac{R{o}^2}{Bu}\left(\nabla X\times {\partial }_{z}u\right)\cdot n_z\right]=\mathrm{<mi>\; </mi>}O\left(R{o}^3\right)\mathrm{.}$

Equation (GV) simply expresses the three-dimensional Lagrangian conservation of the generalized vorticity Q up to order $R{o}^3$ :

$\frac{DQ}{Dt}=\mathrm{<mi>\; </mi>}O\left(R{o}^3\right)$ (15)

where $D/Dt$ must be expanded up to $O\left(R{o}^2\right)$ (see subsection 3.4), and with the explicit form of Q:

$\begin{array}{c}Q\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left(1+Ro\left({\partial }_{x}v-{\partial }_{y}u\right)+\beta y\right)\left(1-\frac{Ro}{Bu}{\partial }_z\left(\frac{\rho }{S}-\frac{Ro}{Bu}\frac{{\rho }^2}{2S}{\partial }_z\frac{1}{S}\right)\right)\\ +\frac{R{o}^2}{Bu}\left({\partial }_{z}v{\partial }_x-{\partial }_{z}u{\partial }_y\right)\left(\frac{\rho }{S}\right)+O\left(R{o}^3\right).\end{array}$ (16)

We remark that Q has the form $A\ast B+C$ , where A is the absolute vorticity, B is the “baroclinic” part and C is related to the baroclinicity vector. Neglecting the horizontal variations of vertical velocity, Q can be expressed as the product of the horizontal vorticity by a scalar gradient.

For very large Bu, the barotropic vorticity is retrieved and we note $\zeta ={\partial }_{x}v-{\partial }_{y}u$ .

From here on, equation (C) is replaced by equation (GV). For consistency, remark that the 5 unknowns are still described by 5 equations. The following step relates the generalized vorticity Q and the Ertel vorticity.

3.2. Conservation of the Ertel Vorticity

In dimensionless variables, the Ertel vorticity is written

$\Pi \mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-S\left[\left(1+Ro\zeta +\beta y\right)\left(1-\frac{Ro}{SBu}{\partial }_z\rho \right)+\frac{R{o}^2}{Bu}\left({\partial }_{z}v{\partial }_x-{\partial }_{z}u{\partial }_y\right)\left(\frac{\rho }{S}\right)\right]+O\left(R{o}^3\right).$ (17)

Note that:

$\left(1-\frac{Ro}{Bu}{\partial }_{z}X\right)\frac{\Pi }{S\left(z\right)}=\mathrm{<mi>\; </mi>}-\left(1-\frac{Ro}{SBu}{\partial }_z\rho \right)Q+O\left(R{o}^3\right).$

or

$Q\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}A\Pi +O\left(R{o}^3\right)$

with

$A\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-\left[\frac{1}{S}-\frac{Ro}{Bu}X{\partial }_z\frac{1}{S}+{\left(\frac{Ro}{Bu}\right)}^2\frac{{\rho }^2}{2S^2}{\partial }_z^2\frac{1}{S}\right].$

We remark that A is a function of density only. The conservation of the generalized vorticity, Equation (15), implies that

$A\frac{D\Pi }{Dt}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-\Pi \frac{DA}{Dt}+O\left(R{o}^3\right);$

some straightforward algebra, with the use of relation (13), leads to:

$\frac{DA}{Dt}=O\left(R{o}^3\right)\mathrm{<mo>\; ;\; </mo>}$

which is mass conservation, and thus immediately to:

$\frac{D\Pi }{Dt\mathrm{<mi>\; </mi>}}=\mathrm{<mi>\; </mi>}O\left(R{o}^3\right).$ (18)

with the form of $\Pi$ given by Equation (17). The Ertel vorticity is thus conserved as for the primitive equations, but here only up to $O\left(R{o}^3\right)$ . We remark that, contrary to quasi-geostrophic motions where total vorticity is advected by a two-dimensional geostrophic velocity, here the Ertel vorticity is advected by the three-dimensional velocity.

With less complex algebra, the equivalent relation between Q and $\Pi$ for quasi-geostrophic dynamics is:

$Q_{qg}=A_{qg}{\Pi }_{qg}$ , with

$Q_{qg}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left(1+Ro\left({\partial }_{x}v-{\partial }_{y}u\right)+\beta y-\frac{Ro}{Bu}{\partial }_z\left(\frac{\rho }{S}\right)\right)+O\left(R{o}^2\right),$

${\Pi }_{qg}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-S\left(1+Ro\left({\partial }_{x}v-{\partial }_{y}u\right)+\beta y-\frac{Ro}{SBu}{\partial }_z\rho \right)+O\left(R{o}^2\right),$

$A_{qg}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{-1}{S}\left(1-\frac{Ro}{Bu}\rho {\partial }_z\frac{1}{S}\right)$

Both in the QG and GV models, $\Pi$ contains both the static and dynamic parts of potential vorticity, while Q expresses the dynamic anomaly, relevant to the description of motions (Y. Morel, priv. comm.). A represents the new distribution of fluid masses induced by these motions.

3.3. Integral Invariant Properties

3.3.1. Potential Enstrophy

Equation (18) immediately leads to the conservation of potential enstrophy

$Z\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\displaystyle \int }{\displaystyle \int }{\displaystyle \int }\frac{{\Pi }^2}{2}\text{d}x\text{d}y\text{d}z$

as for the primitive equations, but here up to $O\left(R{o}^3\right)$ only.

3.3.2. Total Energy

Kinetic energy density is expressed in full form as:

$K\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\rho }_s\frac{\left(u^2+v^2\right)}{2}+{\left(\frac{Ro}{Bu}\right)}^2{\rho }_s\frac{w^2}{2}$

At the model order, the kinetic energy density is multiplied by Ro and the following form is kept:

$RoK\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}Ro{\rho }_s\frac{\left(u^2+v^2\right)}{2}+O\left(R{o}^3\right),$ (19)

where the horizontal velocity expansion in Ro is kept up to second order for consistency.

Potential energy density is expressed as (see Appendix A):

$RoP=Ro{\rho }_s\frac{{\rho }^2}{2BuS}\left(1-\frac{Ro}{Bu}\frac{\rho }{3}{\partial }_z\frac{1}{S}\right)+O\left(R{o}^3\right).$ (20)

The energy conservation equation (E) is obtained via:

$\left(E\right)=\mathrm{<mi>\; </mi>}{\rho }_{s}u\cdot QM$

in three-dimensional form. The main steps of the calculation are:

$Ro\left[\frac{DK}{Dt}\right]+\frac{Ro}{Bu}{\rho }_s\rho w\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-u\cdot \nabla \left({\rho }_{s}p\right)+O\left(R{o}^3\right);$

the product of density by vertical velocity is then expanded into

$\frac{Ro}{Bu}{\rho }_s\rho w\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{Ro}{Bu}\frac{DX}{Dt}{\rho }_s\rho +O\left(R{o}^3\right),$

or, with a few algebraic manipulations,

$\frac{Ro}{Bu}{\rho }_s\rho w=\mathrm{<mi>\; </mi>}Ro\frac{DP}{Dt}+O\left(R{o}^3\right)\mathrm{.}$

The final result is, in Lagrangian form:

$Ro\frac{D\left(K+P\right)}{Dt}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-u\cdot \nabla \left({\rho }_{s}p\right)+O\left(R{o}^3\right)$ (21)

which expresses the conventional result that the variation of energy with time is due to the work of pressure forces, or in Eulerian form

$Ro{\partial }_t\left(K+P\right)+\nabla \cdot \left(uB\right)\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}O\left(R{o}^3\right)$ (22)

where B is the total Bernoulli function $B={\rho }_{s}p+Ro\left(K+P\right)$ . This property is analogous to that for primitive equations. Energy is conserved here up to $O\left(R{o}^3\right)$ in any volume for which the integrated flux of the Bernoulli function across the boundary vanishes. In the QG approximation, the right-hand side of Equation (22) vanishes as the velocity is horizontal and geostrophic at first order. The total energy is then exactly conserved.

3.4. Model Equations in Generalized Streamfunction

Here the three-dimensional Lagrangian derivative is expanded in terms of pressure, while conserving the model consistency at order $O\left(R{o}^3\right)$ , and thus preserving integral invariants at this order. The derivative is split into two parts: a horizontal Lagrangian derivative and a vertical advection term.

$\frac{D}{Dt}=\frac{\text{d}}{\text{d}t}+\frac{Ro}{Bu}w{\partial }_z.$

The horizontal derivative is applied to terms of order $O\left(Ro\right)$ at least, in equations (GV), (H) and (QM) and thus needs to be expanded to $O\left(R{o}^2\right)$ to ensure global consistency at $O\left(R{o}^3\right)$ . The vertical advection term is multiplied by a $Ro/Bu$ factor and is thus expanded to $O\left(R{o}^2\right)$ to ensure the same consistency at $O\left(R{o}^3\right)$ . For both the horizontal and vertical components of the 3D derivative, specific operators, $L_h$ and $L_z$ , are defined from the original Equation (8)-(13). The product L of these operators will permute with the time derivative and the Jacobian. This specific property will help transform Equation (16) into a QG-like equation with a generalized streamfunction, acting on LQ. This form of the conservation of generalized vorticity is called the internal form. With the Bernoulli function, this operator does not completely permute, and the equation obtained is called the external form. For legibility, part of the algebra is given in Appendix B.

3.4.1. Vertical Derivative

The expression of the vertical velocity is given by equation (H), that is, the quantity is substracted from the total derivative:

$\frac{D}{Dt}-\left(H\right){\partial }_z\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{\text{d}}{\text{d}t}+\frac{Ro}{Bu}\frac{DX}{Dt}{\partial }_z+O\left(R{o}^2\right)$

Part integration, permutation and summation (see Appendix B1) yield:

$L_z\frac{D}{Dt}=\frac{\text{d}}{\text{d}t}{L}_z+\frac{Ro}{Bu}XJ\left({\partial }_{z}p,\cdot \right)+O\left(R{o}^2\right)$ (23)

where we define the vertical operator $L_z$ as

$L_z\mathrm{<mi>\; </mi>}=\left(1+\frac{Ro}{Bu}X{\partial }_z\right)$

and the Jacobian $J\left(a,b\right)={\partial }_{x}a{\partial }_{y}b-{\partial }_{x}b{\partial }_{y}a=\left(\nabla a\times \nabla b\right)\cdot n_z$ . Using Equation (9) and (13) for $\rho$ and X, we obtain an expression for $L_z$ in terms of pressure only:

$L_z\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>\; 1}-\frac{Ro}{Bu}\frac{{\partial }_{z}p}{S}{\partial }_z+O\left(R{o}^2\right).$

3.4.2. Horizontal Lagrangian Operator

The horizontal derivative is explicited as:

$\frac{\text{d}}{\text{d}t}=\mathrm{<mi>\; </mi>}\frac{R{o}_T}{Ro}{\partial }_t+u{\partial }_x+v{\partial }_y$

and the following combination is performed:

$\left(1+\beta y\right)\frac{\text{d}}{\text{d}t}-n_z\cdot \left(QM\times \nabla \right)\mathrm{.}$

We define the horizontal operator $L_h$ as

$L_h=1+Ro\left(n_z\cdot \left(u\times \nabla \right)\right)\mathrm{,}$

or, in an expansion in pressure,

$L_h\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>\; 1}-Ro{\nabla }_{h}p\cdot {\nabla }_h+O\left(R{o}^2\right),$

where the subscript h refers to the sole horizontal components. This definition and some algebra (see Appendix B2) lead to

$\left(1+\beta y\right)L_h\frac{\text{d}}{\text{d}t}{L}_h^{-1}=\left(1+\beta y\right)\frac{R{o}_T}{Ro}{\partial }_t+J\left(p-Ro\frac{u^2+v^2}{2},\cdot \right)+O\left(R{o}^2\right).$ (24)

Noting $\Phi =p-Ro\left(u^2+v^2\right)/2$ and using the straightforward relation $I=L_h^{-1}{L}_h$ for the exact form of $L_h$ leads to

$\left(1+\beta y\right)L_h\frac{\text{d}}{\text{d}t}=\mathrm{<mi>\; </mi>}\left(1+\beta y\right)\frac{R{o}_T}{Ro}{\partial }_t{L}_h+J\left(\Phi \mathrm{,}{L}_h\right)+O\left(R{o}^2\right)$ (25)

as a second intermediate step.

3.4.3. Final Result for the Total Lagrangian Derivative

Defining $L=L_h{L}_z$ and applying it to the total derivative yields (see the algebra in Appendix B3):

$\begin{array}{c}\left(1+\beta y\right)L\frac{D}{Dt}=\left(1+\beta y\right)\frac{R{o}_T}{Ro}{\partial }_{t}L+J\left(p-Ro\frac{u^2+v^2}{2},L\right)\\ +\frac{Ro}{Bu}XJ\left({\partial }_{z}p,L\right)+O\left(R{o}^2\right),\end{array}$

or, in terms of energy density,

$\left(1+\beta y\right)L\frac{D}{Dt}=\left(1+\beta y\right)\frac{R{o}_T}{Ro}{\partial }_{t}L+J\left(p-Ro\frac{K+P}{{\rho }_s},L\right)+O\left(R{o}^2\right).$ (26)

The generalized streamfunction can now be defined as

$\psi \mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}p-Ro\frac{K+P}{{\rho }_s}=\mathrm{<mi>\; </mi>}p-Ro\left(\tilde{K}\left(p\right)+\tilde{P}\left(p\right)\right),$ (27)

with

$Ro\tilde{K}\left(p\right)=\frac{Ro}{2}{\left[\nabla p\right]}^2+O\left(R{o}^2\right),$

$Ro\tilde{P}\left(p\right)=\frac{Ro}{2BuS}{\left[{\partial }_{z}p\right]}^2+O\left(R{o}^2\right),$

which is the symmetric of the Bernoulli function with respect to p. Operator L is written in terms of $\psi$

$L\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>\; 1}-Ro{\nabla }_h\psi \cdot {\nabla }_h-\frac{Ro}{\left(SBu\right)}{\partial }_z\psi {\partial }_z+O\left(R{o}^2\right).$ (28)

The final form for the Lagrangian derivative, along with Equations (28)-(29), is:

$\left(1+\beta y\right)L\frac{D}{Dt}=\frac{R{o}_T}{Ro}\left(1+\beta y\right){\partial }_{t}L+J\left(\psi ,L\right)+O\left(R{o}^2\right),$ (29)

which applies to any conserved Lagrangian quantity (such as Q or $\Pi$ ), expanded up to $O\left(R{o}^3\right)$ , having a zeroth-order term depending only on z. Note that, if we choose $\beta <O\left(Ro\right)$ as in subsection 2.2, Equation (29) can also be expressed as:

$L\frac{D}{Dt}=\frac{R{o}_T}{Ro}{\partial }_{t}L+\left(1-\beta y\right)J\left(\psi ,L\right)+O\left(R{o}^2\right),$

an equation of interest for the dynamics of intense vortices.

3.4.4. Explicit Form of the Generalized Vorticity in Generalized Streamfunction

With the 3D Lagrangian derivative written in terms of $\psi$ (or p) hereabove, the generalized vorticity Q can also be expressed in terms of $\psi$ . This expression is obtained in two steps: first, the generalized vorticity $Q\left(\zeta \mathrm{,}\rho \right)$ (Equation (16)) is given as a function of pressure p, using the expressions $\zeta \left(p\right)\mathrm{,}\rho \left(p\right)$ . Then the pressure is given as a function of the generalized streamfunction $\psi$ by Equation (27) to yield $Q\left(\psi \right)$ . The same calculations could also be performed on the Ertel vorticity $\Pi$ , but these are not given here since Q directly derives from the original equations. Also note that from here on the $\nabla$ operator is considered horizontal and the subscript h is dropped.

Using

$\zeta \mathrm{<mi>\; </mi>}=\left(1-\beta y\right){\nabla }^{2}p-2RoJ\left({\partial }_{x}p\mathrm{,}{\partial }_{y}p\right)-\beta {\partial }_{y}p+O\left(R{o}^2\right)\mathrm{,}$

and

$\rho \mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-{\partial }_{z}p,$

we get:

$\begin{array}{c}Q\left(p\right)=\left[1+Ro\left(1-\beta y\right){\nabla }^{2}p-2R{o}^{2}J\left({\partial }_{x}p,{\partial }_{y}p\right)-Ro\beta {\partial }_{y}p+\beta y\right]\\ \ast \left[1+\frac{Ro}{Bu}{\partial }_z\left(\frac{1}{S}{\partial }_{z}p+\frac{Ro}{Bu}{\left({\partial }_{z}p\right)}^2\frac{1}{2S}{\partial }_z\frac{1}{S}\right)\right]\\ -\frac{R{o}^2}{Bu}\frac{1}{S}{\left[\nabla {\partial }_{z}p\right]}^2+O\left(R{o}^3\right).\end{array}$

The quasi-geostrophic form of the total vorticity Q is:

$Q_{qg}\left(p\right)=1+Ro{\nabla }^{2}p+\beta y+\frac{Ro}{Bu}{\partial }_z\left(\frac{1}{S}{\partial }_{z}p\right)+O\left(R{o}^2\right).$

As all terms containing $\psi$ in Q will be factor of $Ro$ , the expression of $\psi \left(p\right)$ given in Equation (27) is sufficient to maintain consistency at $O\left(R{o}^3\right)$ . It can easily be inverted to $p\left(\psi \right)$ and substituted into the previous expression for Q to reach

$\begin{array}{c}Q\left(\psi \right)=[1+Ro\left(1-\beta y\right){\nabla }^2\psi -2R{o}^{2}J\left({\partial }_x\psi ,{\partial }_y\psi \right)-Ro\beta {\partial }_y\psi +\beta y\\ +R{o}^2{\nabla }^2\left(\tilde{K}\left(\psi \right)+\tilde{P}\left(\psi \right)\right)]\ast [1+\frac{Ro}{Bu}{\partial }_z(\frac{1}{S}{\partial }_z\psi +\frac{Ro}{Bu}{\left({\partial }_z\psi \right)}^2\frac{1}{2S}{\partial }_z\frac{1}{S}\\ +Ro\frac{1}{S}{\partial }_z\left(\tilde{K}\left(\psi \right)+\tilde{P}\left(\psi \right)\right))]-\frac{R{o}^2}{Bu}\frac{1}{S}{\left[\nabla {\partial }_z\psi \right]}^2+O\left(R{o}^3\right).\end{array}$

We can now apply the operator L to $Q\left(\psi \right)$ to obtain the expression which will be conserved in the evolution equation hereafter:

$\begin{array}{c}LQ\left(\psi \right)=\left[1+Ro{\nabla }^2\psi +\beta y\right]\ast \left[1+\frac{Ro}{Bu}{\partial }_z\left(\frac{1}{S}{\partial }_z\psi \right)\right]-R{o}^2{\left({\nabla }^2\psi \right)}^2\\ -2R{o}^2\nabla \psi \cdot \nabla \left({\nabla }^2\psi \right)+2R{o}^2{\nabla }^2\tilde{K}-R{o}^2\beta y{\nabla }^2\psi -2R{o}^2\beta {\partial }_y\psi \\ +\frac{R{o}^2}{Bu}\frac{1}{S}{\left[\nabla \left({\partial }_z\psi \right)\right]}^2+\frac{R{o}^2}{B{u}^2}{\left[{\partial }_z\left(\frac{1}{S}{\partial }_z\psi \right)\right]}^2+O\left(R{o}^3\right).\end{array}$ (30)

We recognize in this expression the quasi-geostrophic terms at $O\left(Ro\right)$ and the nonlinear coupling terms at $O\left(R{o}^2\right)$ .

In generalized streamfunction, the conservation of generalized vorticity is thus:

$\frac{R{o}_T}{Ro}{\partial }_t\left[LQ\left(\psi \right)\right]+\left(1-\beta y\right)J\left(\psi ,LQ\left(\psi \right)\right)=O\left(R{o}^3\right),$ (31)

which, with expressions (27) and (30), is the final result for Section 3.4. Obviously energy conservation could be derived directly from Equation (33) by multiplying it by the generalized streamfunction, as it is done in QG. It is less tedious and more straightforward to derive it as in Subsection 3.3.2, though the result is the same at $O\left(R{o}^3\right)$ .

This problem is similar in its mathematical form to the QG equations, i.e. it can be initialized either in streamfunction or in vorticity, though the relation between both is a nonlinear elliptic equation. If the streamfunction is given initially, the generalized vorticity is computed via the relation above Equation (30); conversely, if Q is given, $\psi$ is found by inverting the same formula iteratively. For the boundary conditions, $\psi$ and $\zeta ={\nabla }^2\psi$ must be given at the lateral boundaries of the domain, and the density anomaly $\rho ~{\partial }_z\psi$ must be specified at the surface and at the bottom of the domain. These latter conditions are usually a rigid lid at the surface and a given bottom topography. In the layer model, the asymptotic form of an infinitely deep lower layer at rest (reduced-gravity case) can also be used (see part 2).

3.5. An Alternative Form of the Dynamics in Bernoulli Function

We recall here that all tilded quantities are massless (i.e. are divided by ${\rho }_s$ ), as it was already defined in Equation (27). We also use the subscript h for the horizontal quantities, which are relevant for the layer model. Horizontal velocities can also be written with respect to the horizontal Bernoulli function ${\tilde{B}}_h=p+Ro\tilde{K}\left(p\right)$ as:

$u\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{-1}{1+Ro{\nabla }^2{\tilde{B}}_h+\beta y}{\partial }_y{\tilde{B}}_h-Ro{\partial }_t{\partial }_x{\tilde{B}}_h+O\left(R{o}^2\right),$

$v\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{1}{1+Ro{\nabla }^2{\tilde{B}}_h+\beta y}{\partial }_x{\tilde{B}}_h-Ro{\partial }_t{\partial }_y{\tilde{B}}_h+O\left(R{o}^2\right).$

The derivation of the model equations is similar to that in Sections 3.4.1 and 3.4.2. The Lagrangian derivative is written with ${\tilde{B}}_h$ replacing $\Phi =p-Ro\tilde{K}\left(p\right)$ in the Jacobian, and using the $L_h$ operator with ${\tilde{B}}_h$ substituted to p:

$\frac{\text{d}}{\text{d}t}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{L}_h^{-1}\left({\tilde{B}}_h\right)\frac{R{o}_T}{Ro}{\partial }_t{L}_h\left({\tilde{B}}_h\right)+\frac{1}{1+Ro{\nabla }^2{\tilde{B}}_h+\beta y}J\left({\tilde{B}}_h,\right)+O\left(R{o}^2\right).$

The two forms in ${\tilde{B}}_h$ and $\Phi$ are equivalent (see Appendix B4).

As for the horizontal Lagrangian operator in Section 3.4.2, the vertical operator $L_z$ must now be applied to obtain the total Lagrangian derivative:

$\begin{array}{l}L\frac{D}{Dt}=L_h\left(L_z\frac{D}{Dt}\right)\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{L}_h\frac{\text{d}}{\text{d}t}{L}_z+\frac{Ro}{Bu}\left(XJ\left({\partial }_{z}p,\right)\right)+O\left(R{o}^2\right)\\ =\left[\frac{R{o}_T}{Ro}{\partial }_t{L}_h+\frac{L_h}{1+Ro{\nabla }^2{B}_h+\beta y}J\left(B_h\mathrm{,}\right)\right]L_z+\frac{Ro}{Bu}\left(L_{h}XJ\left({\partial }_{z}p\mathrm{,}\right)\right)+O\left(R{o}^2\right)\\ =\frac{R{o}_T}{Ro}{\partial }_t{L}_h{L}_z+\frac{L_h}{1+Ro{\nabla }^2{B}_h+\beta y}J\left(B_h\mathrm{,}{L}_z\right)+\frac{Ro}{Bu}\left(XJ\left({\partial }_{z}p\mathrm{,}\right)\right)+O\left(R{o}^2\right)\mathrm{,}\end{array}$

where $\left(Ro/Bu\right)L_{h}X$ has finally been replaced by $\left(Ro/Bu\right)X$ at $O\left(R{o}^2\right)$ .

Thus the Bernoulli formulation of the evolution equation for the generalized vorticity is:

$L\frac{DQ\left(\tilde{B}\right)}{Dt}=\frac{R{o}_T}{Ro}{\partial }_{t}LQ\left(\tilde{B}\right)+\frac{L_h}{1+Ro{\nabla }^2{\tilde{B}}_h+\beta y}J\left(\tilde{B},L_{z}Q\left(\tilde{B}\right)\right)=O\left(R{o}^3\right).$ (32)

Here $Q\left(\tilde{B}\right)$ has a similar form to $Q\left(\psi \right)$ , and

$\tilde{B}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\tilde{B}}_h-Ro\tilde{P}\left(p\right).$

Note that $\tilde{B}$ is not the usual total Bernoulli function, because of the minus sign in front of the potential energy component. In the layer model, $\tilde{B}$ reduces to ${\tilde{B}}_h$ , the usual horizontal Bernoulli function.

4. Layer Model

The form of the relative vorticity and of the Lagrangian derivative in terms of pressure given in part 1 is projected onto homogeneous layers. The basic property of the layer model is the independence of the horizontal velocities (and of the horizontal vorticity $\zeta$ ) on the vertical coordinate in each homogeneous layer.

Obviously, energy and enstrophy are conserved to the same order in Ro as for the continuously stratified case.

4.1. Layerwise Conservation of Ertel Vorticity

To obtain the conservation of the Ertel vorticity in each layer, successive steps are followed:

1) the hydrostatic equation is replaced by the constraint on the pressure jump at the density interface between layers $i-1$ and $i$ :

$p_i-p_{i-1}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\rho }_0{g^{\prime }}_{i-1/2}{\eta }_{i-1},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{g^{\prime }}_{i-1/2}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}g\left({\rho }_i-{\rho }_{i-1}\right)/{\rho }_0,$ (33)

where ${\eta }_{i-1}$ is the elevation of this density interface from its depth at rest $z=z_{i-1}\le 0$ .

2) Equation (H) is replaced by the equation of the vertical velocity component at the density interfaces

$w_i\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{\text{d}{\eta }_i}{\text{d}t},$

3) equation (GV) in each layer is obtained by vertical integration on the layer thickness:

${\left(GV\right)}_i\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\displaystyle {\int }_{z_i+{\eta }_i}^{z_{i-1}+{\eta }_{i-1}}}\left[\left(\nabla \times Q{M}_i\right)\cdot n_z+\left(1+Ro{Y}_i\right)\left(H\right)\right]\text{d}z.$

The discretization of the model equations of the previous section on layers leads to the same results as the truncation of the multilayer shallow-water equation model to $O\left(R{o}^2\right)$ . Both simplifications conserve potential vorticity and the integral invariant properties at the model order.

Rather than using potential vorticity, its inverse, potential thickness

$Z_i\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{h_i}{1+Ro{\zeta }_i+\beta y}$

is used with

$h_i\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{z}_{i-1}-z_i+{\eta }_{i-1}-{\eta }_i,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\zeta }_i={\partial }_x{v}_i-{\partial }_y{u}_i$

Indeed, in case of strong interface deformation, $h_i$ could locally vanish or become considerably small, whereas the relative vorticity remains positive since $Ro\le 0.5$ (and $\beta$ is usually much smaller). Therefore, $Z_i$ cannot become infinite (inertial instability cannot arise).

In each layer, potential thickness conservation is written:

$\frac{\text{d}{Z}_i}{\text{d}t}=O\left(R{o}^3\right)$ (34)

Potential thickness $Z_i$ is then expanded in $Ro$ :

$Z_i\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{h}_i\left(1-Ro{\zeta }_i-\beta y+R{o}^2{\zeta }_i^2+2Ro\beta y{\zeta }_i+{\beta }^2{y}^2\right)+O\left(R{o}^3\right).$

In terms of pressure, potential thickness conservation in each layer is written:

$\frac{R{o}_T}{Ro}{\partial }_t\left(L_i\left(p_i\right)Z_i\right)+\left(1-\beta y\right)J\left(p_i-\frac{Ro}{2}\nabla p_i\cdot \nabla p_i,L_i\left(p_i\right)Z_i\right)=\mathrm{<mi>\; </mi>}O\left(R{o}^3\right),$ (35)

and the operator $L_i$ is, in terms of pressure:

$L_i\left(p_i\right)\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>\; 1}-Ro\nabla p_i\cdot \nabla +O\left(R{o}^2\right).$ (36)

The term $L_i{Z}_i$ has the following form in layers $i=2$ to $i=N-1$ (intermediate layers)

$L_i{Z}_i\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left(1-Ro\nabla p_i\cdot \nabla \right)\left[\frac{1+Ro\left[\left(F_i^{+}+F_i^{-}\right)p_i-F_i^{-}{p}_{i-1}-F_i^{+}{p}_{i+1}\right]}{1+Ro{\zeta }_i+\beta y}\right],$

${\zeta }_i\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\nabla }^2{p}_i-2RoJ\left({\partial }_x{p}_i,{\partial }_y{p}_i\right)-\nabla \cdot \left[\beta y\nabla p_i\right],$

with

$F_i^{+}=\frac{f_0^2{L}^2}{{g^{\prime }}_{i+1/2}{H}_i},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{F}_i^{-}=\frac{f_0^2{L}^2}{{g^{\prime }}_{i-1/2}{H}_i}.$

For the upper layer, with the rigid-lid assumption

$L_1{Z}_1\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left(1-Ro\nabla p_1\cdot \nabla \right)\left[\frac{1-Ro\left[F_1^{+}\left(p_2-p_1\right)\right]}{1+Ro{\zeta }_1+\beta y}\right],$

and for the bottom layer:

$L_N{Z}_N\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\left(1-Ro\nabla p_N\cdot \nabla \right)\left[\frac{1+Ro{F}_N^{-}\left(p_N-p_{N-1}\right)-F_B^{+}{\eta }_B}{1+Ro{\zeta }_N+\beta y}\right],$

with

$F_B^{+}=\frac{\Delta H_B}{H_N},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}{\eta }_B=\frac{h_B}{\Delta H_B},$

for the coupling with the bottom topography.

Numerical implementation of these equations in the N-layer form, with iterative solution of pressure from potential thickness is described and applied in (Filatoff et al., 1997).

4.2. Comparison with the Generalized Geostrophy Model

We show here that a simple restriction of the two-layer model yields the well-known generalized geostrophy equations ( [11] ). We recall that these equations were derived in the context of a rigid lid and flat bottom ocean. For an easier comparison, we use similar notations to these authors’ (beware that some might differ from ours in the previous sections). We define the interface elevation as:

$\eta =\gamma p_2-p_1,\gamma =R{o}_2/R{o}_1,R{o}_1=\frac{g^{\prime }\Delta H}{f^2{L}^2},R{o}_2=\frac{P}{f^2{L}^2},$

where $\Delta H$ is the magnitude of the interface displacement and P the amplitude of the pressure in the lower layer. We also define a local Burger number and a thickness ratio of the two layers:

$s=\frac{g^{\prime }{H}_1}{f^2{L}^2},\delta =\frac{H_1}{H_1+H_2}.$

We reorder our equations and keep as variables the interface elevation h and the geostrophic streamfunction in the lower layer $\varphi$ . The upper layer dynamical equation is written:

$\begin{array}{l}\frac{R{o}_T}{Ro}{\partial }_t\left(-s{\nabla }^2\eta +s\gamma {\nabla }^2\left(\varphi \eta \right)+\eta \right)+\gamma J\left(\varphi ,\eta \right)-sJ\left(\eta -\gamma \varphi ,{\nabla }^2\left(\eta -\gamma \varphi \right)\right)\\ -RoJ\left(\eta ,\eta {\nabla }^2\eta +\frac{\nabla \eta \cdot \nabla \eta }{2}\right)-sJ\left(\eta -\gamma \varphi ,\beta y\right)-RoJ\left(\eta ,\eta \beta y\right)\\ =R{o}_T{\partial }_t\nabla \cdot \left(\eta \nabla \eta \right)+R{o}_T\gamma {\partial }_t\nabla \cdot \left(\eta \nabla \varphi \right)+R{o}_T\beta y{\partial }_t\eta -Ro\gamma [J\left(\varphi ,\nabla \cdot \left(\eta \nabla \eta \right)\right)\\ +J\left(\eta ,\eta {\nabla }^2\varphi \right)]+Ro{\gamma }^{2}J\left(\varphi ,\nabla \cdot \left(\eta \nabla \varphi \right)\right)+Ro\gamma \beta \left[yJ\left(\varphi ,\eta \right)+J\left(\varphi ,\eta y\right)\right],\end{array}$

for which the r.h.s. terms are on order of $R{o}_T\mathrm{,}R{o}_T\gamma \mathrm{,}R{o}_T\beta \mathrm{,}Ro\gamma \mathrm{,}Ro{\gamma }^2\mathrm{,}Ro\gamma \beta$ and are thus neglected by the aforementioned authors. In the present framework, this amounts to assuming that both $R{o}_T~R{o}^2$ and $\gamma ~Ro$ , i.e. that the time-derivative is slow and that the motion is concentrated in the upper layer (a restriction of the present model). Our model thus agrees with the generalized geostrophic model under its assumptions, and provides the form of the complementary terms one order beyond these limits.

If the lower layer is at rest, the reduced gravity equation is restored

$\begin{array}{l}\frac{R{o}_T}{Ro}{\partial }_t\left(-s{\nabla }^2\eta +\eta \right)-sJ\left(\eta ,{\nabla }^2\eta \right)-RoJ\left(\eta ,\eta {\nabla }^2\eta +\frac{\nabla \eta \cdot \nabla \eta }{2}\right)\\ -sJ\left(\eta ,\beta y\right)-RoJ\left(\eta ,\eta \beta y\right)=\mathrm{<mi>\; </mi>\; 0.}\end{array}$ (37)

In [10] , the supplementary conditions $R{o}_T=R{o}^2$ and $s=Ro$ are used which lead to the well-known frontal geostrophic model in a one-and-a-half layer configuration (with $h=1+\eta$ )

${\partial }_{t}h\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}J\left(h,h{\nabla }^{2}h+\frac{\nabla h\cdot \nabla h}{2}\right)+J\left(h,h\beta y\right).$ (38)

A second equation is obtained by summing the two layer equations and by substracting Ro times the dynamical equation and the continuity equation of the upper layer. It yields:

$\begin{array}{l}\frac{s}{\delta }\frac{R{o}_T}{Ro}\gamma {\partial }_t{\nabla }^2\varphi -s\frac{R{o}_T}{Ro}{\partial }_t{\nabla }^2\eta +R{o}_T\nabla \cdot \left(\eta \nabla {\partial }_t\eta \right)-sJ\left(\eta ,{\nabla }^2\eta \right)\\ -RoJ\left(\eta ,\eta {\nabla }^2\eta +\frac{\nabla \eta \cdot \nabla \eta }{2}\right)+Ro\gamma \eta J\left(\varphi ,{\nabla }^2\eta \right)+Ro\gamma \eta J\left(\eta ,{\nabla }^2\varphi \right)\\ +Ro\gamma J\left(\varphi ,\nabla \eta \right)\cdot \nabla \eta +Ro\gamma J\left(\eta ,\nabla \varphi \right)\cdot \nabla \eta -s\gamma J\left(\varphi ,{\nabla }^2\eta \right)-s\gamma J\left(\eta ,{\nabla }^2\varphi \right)\\ +s\frac{{\gamma }^2}{\delta }J\left(\varphi ,{\nabla }^2\varphi \right)-sJ\left(\eta ,\beta y\right)-RoJ\left(\eta ,\eta \beta y\right)+s\frac{\gamma }{\delta }J\left(\varphi ,\beta y\right)\\ =-R{o}_T\nabla \cdot \left({\partial }_t\eta \nabla \eta \right)-Ro\gamma {\nabla }^2\eta J\left(\varphi ,\eta \right)-Ro\gamma J\left(\varphi ,\nabla \eta \right)\cdot \nabla \eta +Ro\gamma J\left(\eta ,\nabla \varphi \right)\cdot \nabla \eta .\end{array}$

Using algebraic combinations and similar scaling arguments as previously, the r.h.s. terms are negligible according to these authors. When the upper layer is extremely thin (vanishing $\delta$ ), the lower layer motion is decoupled from the upper layer influence and obeys a barotropic vorticity equation:

$\frac{R{o}_T}{Ro}{\partial }_t{\nabla }^2\varphi +J\left(\varphi \mathrm{,}{\nabla }^2\varphi +\beta y\right)=0.$

5. Conclusions

We have presented here a set of equations satisfying the constraints that we had initially imposed. The model retains all terms at $O\left(Ro\right)$ in momentum and at $O\left(R{o}^2\right)$ in vorticity-divergence or integral invariants; it appears under two forms, one in generalized streamfunction, and one in Bernoulli function. The model conserves total energy (kinetic and potential) at the consistent order ( $O\left(R{o}^3\right)$ ), if the total work of pressure forces is zero, or if the total flux of Bernoulli function across the boundaries vanishes; a new general expression of potential energy (up to $O\left(R{o}^2\right)$ ) has been provided. Potential enstrophy (as derived from Ertel’s vorticity) is conserved at $O\left(R{o}^3\right)$ . This conservation of potential enstrophy is essential for the study of vortex flows.

It should be remarked that in the framework initially defined, the present model is formally equivalent to other intermediate models, but for terms at the order in Ro neglected in the present calculation. The model with continuous stratification is of use for comparison with the balance equations, and for general applications of Ertel vorticity conservation for oceanic processes. As for balance equations, a nonlinear elliptic equation has to be solved, here to recover streamfunction from potential vorticity. This nonlinearity is inherent to all intermediate models written in a single variable, if multiple time scale expansion is not used.

The equations derived for continuous stratification have then been projected onto a discrete representation of vertical stratification. In the layer representation, the dynamical model equations conserve potential vorticity for each fluid particle and integral invariants for the whole fluid. For practical purposes, potential thickness was used rather than potential vorticity; this avoids singularity when layers become vanishingly thin. The model equations were compared successfully to the frontal and to the generalized geostrophic equations derived by [11] ; this comparison also identified the terms which appear in the equations at higher order than that considered by these authors.

This model has been implemented in generalized streamfunction, both with pseudospectral and finite-difference horizontal schemes. It has been applied to various studies of vortex dynamics:

・ barotropic vortex merger, showing a cyclone/anticyclone asymmetry when Ro increases [18] ; this asymmetry is similar in the SWE and GV models but does not exist in the QG model.

・ intrathermocline vortex merger (the merger of two vortices in the middle layer of a three-layer ocean); the merger of two distant anticyclones is possible both in the SWE and the GV models, but not in the QG model.

・ the effects of strong bathymetry gradients on vortex dynamics and the effects of large density interface deviation on baroclinic vortex merger, in comparison with a SWE model [19] ; the comparison of asymptotic and complete dynamical model results is then very satisfactory.

・ In summary, the GV model provides similar physical results as the complete shallow-water modem for vortex dynamics, in cases where the QG model misses the parity bias.

Acknowledgements

This work was made possible by a contract from SHOM/CMO (French Navy Oceanographic Research Center) under the program “Processus associes aux modeles intermediaires et regionaux”. Thanks are due to Pr Bach Lien Hua for an insightful suggestion leading to the result in Appendix A3. Discussions with Prs James C. Mc Williams and Vladimir Zeitlin are gratefully acknowledged.

Cite this paper

Filatoff, N. and Carton, X. (2018) A Conservative Model for Nonlinear Dynamics in a Stratified, Rotating Fluid. Open Journal of Marine Science, 8, 253-275. https://doi.org/10.4236/ojms.2018.82013

References

Appendix A: Available Potential Energy

A1. Dynamical Expression

The restoring force acting on a fluid particle displaced from elevation $z_0$ to $z$ in a stably stratified fluid is buoyancy:

$F_z\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}-g\left({\rho }_T\left(z\right)-{\rho }_T\left(z_0\right)\right).$

The work produced by this force from $z_0$ to $z$ is

$\begin{array}{c}W\left(F_z\right)={\displaystyle {\int }_{z_0}^z}{F}_z\text{d}z=-g{\left({\partial }_z{\rho }_s\right)}_{z=z_0}\frac{{\left(z-z_0\right)}^2}{2}-g{\left({\partial }_z^2{\rho }_s\right)}_{z=z_0}\frac{{\left(z-z_0\right)}^3}{6}\\ +g{\rho }_{z=z_0}\left(z-z_0\right),\end{array}$

in a Taylor expansion in ${\rho }_s$ . This expansion is performed to the second order, like that in B2, for two reasons:

1) it will later be expressed in terms of anomaly densities up to ${\rho }^2$ . This is consistent with the model order, with the density decomposition given in Subsection 2.2;

2) the constant and first-order terms in the expansion are the static and quasi-geostrophic contributions to potential energy, that we already know.

We note that the above expression has two contributions, one hydrostatic due to the ambient stratification, one dynamical due to the motion of fluid masses.

A2. Thermodynamical Arguments

The thermodynamical variables of the system are entropy, density and tracer concentration (salinity) for the extensive variables, and temperature, pressure and chemical potential for the corresponding intensive variables. We assume that salinity is constant and that any adiabatic transformation is isentropic as the system is non dissipative. To maintain the fluid particle in an invariant thermodynamical state, density must be conserved for the particle (but obviously not for its environment or else potential energy would vanish). In other words, we look for the elevation z corresponding to a static thermodynamic equilibrium for the fluid particle, identical to its initial (dynamical) thermodynamic equilibrium at elevation $z_0$ . This is written:

${\rho }_T\left(z_0\right)-{\rho }_s\left(z_0\right)\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}{\rho }_s\left(z\right)-{\rho }_s\left(z_0\right),$

or

${\rho }_{z=z_0}={\left({\partial }_z{\rho }_s\right)}_{z=z_0}\left(z-z_0\right)+{\left({\partial }_z^2{\rho }_s\right)}_{z=z_0}\frac{{\left(z-z_0\right)}^2}{2},$

which is inverted to yield

$z-z_0\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\rho {\left({\partial }_z{\rho }_s\right)}^{-1}-\frac{{\rho }^2}{2}\left({\partial }_z^2{\rho }_s\right){\left({\partial }_z{\rho }_s\right)}^{-3},$

truncating terms at order ${\rho }^3$ and higher; this expression is consistent with the model order. We have also omitted the index $z=z_0$ on densities for simplicity.

A3. Final Expression of APE

The $z-z_0$ expression is now introduced into the previous formula for the work of buoyancy forces

$W\left(F_z\right)=\mathrm{<mi>\; </mi>}\frac{g{\rho }^2}{2}{\left({\partial }_z{\rho }_s\right)}^{-1}\left[1-\frac{1}{3}\rho \left({\partial }_z^2{\rho }_s\right){\left({\partial }_z{\rho }_s\right)}^{-2}\right],$

which can be put in dimensionless form and yield the expression

$RoP\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}Ro\frac{{\rho }_s{\rho }^2}{2SBu}\left(1-\frac{Ro}{Bu}\frac{\rho }{3}{\partial }_z\left(1/S\right)\right)+O\left(R{o}^3\right),$

used in Section 3.3.2

Appendix B: Derivation of the Lagrangian Derivative in Terms of Pressure

These appendices give more detail on the algebra of Section 3.4. The total derivative $D/Dt$ is given by Equation (12), the quantities X and Y by Equations (13) and (14), operator $L_z$ is given in Subsection 3.4.1, and $L_h$ in 3.4.2, and the momentum conservation QM is written in Equations (7)-(9).

B1. Vertical Derivative

We expand here the combination given in Subsection 3.4.1, $\frac{D}{Dt}-\left(H\right){\partial }_z$ , which

gives by part integration:

$\frac{\text{d}}{\text{d}t}+\frac{Ro}{Bu}\frac{DX}{Dt}{\partial }_z=\frac{\text{d}}{\text{d}t}+\frac{Ro}{Bu}\frac{D}{Dt}\left(X{\partial }_z\right)-\frac{Ro}{Bu}X\frac{D}{Dt}\left({\partial }_z\right).$

By permutation of Lagrangian and vertical derivatives, and using the definition of $L_z$ , we have:

$L_z\frac{D}{Dt}=\frac{\text{d}}{\text{d}t}{L}_z+\frac{Ro}{Bu}X\left({\partial }_{z}u{\partial }_x+{\partial }_{z}v{\partial }_y\right)+O\left(R{o}^2\right).$

The operator is transformed by the following difference:

$L_z\frac{D}{Dt}-\frac{Ro}{Bu}X{n}_z\cdot {\partial }_z\left(QM\times \nabla \right)$

to lead to

$L_z\frac{D}{Dt}=\frac{\text{d}}{\text{d}t}{L}_z+\frac{Ro}{Bu}XJ\left({\partial }_{z}p,\cdot \right)+O\left(R{o}^2\right),$

which is the form written in Equation (23).

B2. Horizontal Derivative

The combination

$\left(1+\beta y\right)\frac{\text{d}}{\text{d}t}-n_z\cdot \left(QM\times \nabla \right).$

used in subsection 3.4.2 is expressed as:

$\left(1+\beta y\right)\frac{R{o}_T}{Ro}{\partial }_t+J\left(p,\cdot \right)+Ro{n}_z\cdot \left(\frac{\text{d}u}{\text{d}t}\times \nabla \right)+O\left(R{o}^2\right),$

or by part integration,

$\left(1+\beta y\right)\frac{R{o}_T}{Ro}{\partial }_t+J\left(p,\cdot \right)+Ro\frac{\text{d}}{\text{d}t}\left[n_z\cdot \left(u\times \nabla \right)\right]-Ro{n}_z\cdot \left(u\times \frac{\text{d}}{\text{d}t}\nabla \right)+O\left(R{o}^2\right).$

Permutation of horizontal derivatives and horizontal Lagrangian operator yields:

$\begin{array}{l}\left(1+\beta y\right)\left[1+Ro{n}_z\cdot \left(u\times \nabla \right)\right]\frac{\text{d}}{\text{d}t}\left[1-Ro{n}_z\cdot \left(u\times \nabla \right)\right]\\ =\left(1+\beta y\right)\frac{R{o}_T}{Ro}{\partial }_t+J\left(p-Ro\frac{u^2+v^2}{2},\cdot \right)+O\left(R{o}^2\right),\end{array}$

which, with the expression for $L_h$ , is Equation (24).

B3. Final Result for the Total Lagrangian Derivative

With the definition given in 3.4.1 for the $L_z$ operator, we have

$L_z\frac{D}{Dt}\mathrm{<mi>\; </mi>}=\mathrm{<mi>\; </mi>}\frac{\text{d}}{\text{d}t}{L}_z+\frac{Ro}{Bu}X{n}_z\cdot \left[\nabla {\partial }_{z}p\times \nabla \right]+O\left(R{o}^2\right),$

so that by applying the product of operators $L_h{L}_z$ to the total Lagrangian derivative, we obtain

$\begin{array}{l}\left(1+\beta y\right)L_h{L}_z\frac{D}{Dt}=\left(1+\beta y\right)\frac{Ro}{R{o}_T}{\partial }_t{L}_h{L}_z+J\left(p-Ro\frac{u^2+v^2}{2},L_h{L}_z\right)\\ +\left(1+\beta y\right)L_h\frac{Ro}{Bu}X{n}_z\cdot \left[\nabla {\partial }_{z}p\times \nabla \right]+O\left(R{o}^2\right),\end{array}$

which is the equation above Equation (26).

B4. Equivalence between the Expressions of Generalized Vorticity in Bernoulli Function and in Streamfunction

To prove the identity between the two expressions of the generalized vorticity equation, one just needs to check that

$L_h\left({\tilde{B}}_h\right)\frac{\text{d}}{\text{d}t}=\frac{R{o}_T}{Ro}{\partial }_t{L}_h\left({\tilde{B}}_h\right)+\frac{1}{1+Ro{\nabla }^2{\tilde{B}}_h+\beta y}J\left({\tilde{B}}_h,L_h\left({\tilde{B}}_h\right)\right)$

and

$L_h\left(\Phi \right)\frac{\text{d}}{\text{d}t}=\frac{R{o}_T}{Ro}{\partial }_t{L}_h\left(\Phi \right)+\left(1-\beta y\right)J\left(\Phi ,L_h(\Phi )\right)$

are equivalent. For the partial time derivative, this is obvious to order $R{o}^2$ . We compare the Jacobian terms, the difference of which is:

$J\left({\tilde{B}}_h-\Phi ,\right)-Ro\nabla {\tilde{B}}_h\cdot \nabla J\left({\tilde{B}}_h,\right)-Ro{\nabla }^2{\tilde{B}}_{h}J\left({\tilde{B}}_h,\right)+RoJ\left(\Phi ,\nabla {\tilde{B}}_h\cdot \nabla \right),$

or also with ${\tilde{B}}_h-\Phi =Ro\left(\nabla {\tilde{B}}_h\cdot \nabla {\tilde{B}}_h\right)$ ,

$Ro\left[\nabla {\tilde{B}}_h\cdot J\left(\nabla {\tilde{B}}_h,\right)-{\nabla }^2{\tilde{B}}_{h}J\left({\tilde{B}}_h,\right)+J\left({\tilde{B}}_h,\nabla {\tilde{B}}_h\right)\cdot \nabla \right],$

which vanishes for any test function. Thus the two systems are equivalent.