SciELO - Scientific Electronic Library Online

 
vol.67 número3Modeling the Electronic structure and stability of three aluminum nitride phasesTemperature profiles due to continuous hot water injection into homogeneous fluid-saturated porous media through a line source índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

Links relacionados

  • No hay artículos similaresSimilares en SciELO

Compartir


Revista mexicana de física

versión impresa ISSN 0035-001X

Rev. mex. fis. vol.67 no.3 México may./jun. 2021  Epub 21-Feb-2022

https://doi.org/10.31349/revmexfis.67.351 

Research

Fluids Dinamics

Multilayer shallow-water model with stratification and shear

F. J. Beron-Vera1 

1Department of Atmospheric Sciences Rosenstiel School of Marine & Atmospheric Science University of Miami, Miami, FL 33145 USA. e-mail: fberon@miami.edu


Abstract

The purpose of this paper is to present a shallow-water-type model with multiple inhomogeneous layers featuring variable linear velocity vertical shear and stratification in horizontal space and time. This is achieved by writing the layer velocity and buoyancy fields as linear functions of depth, with coefficients that depend arbitrarily on horizontal position and time. The model is a generalization of Ripa’s (1995) single-layer model to an arbitrary number of layers. Unlike models with homogeneous layers, the present model can represent thermodynamics processes driven by heat and freshwater fluxes through the surface or mixing processes resulting from fluid exchanges across contiguous layers. By contrast with inhomogeneous-layer models with depth-independent velocity and buoyancy, the model derived here can sustain explicitly at a low frequency a current in thermal wind balance (between the vertical vertical shear and the horizontal density gradient) within each layer. In the absence of external forcing and dissipation, energy, volume, mass, and buoyancy variance constrain the dynamics; conservation of total zonal momentum requires also the usual zonal symmetry of the topography and horizontal domain. The inviscid, unforced model admits a formulation suggestive of a generalized Hamiltonian structure, which enables the classical connection between symmetries and conservation laws via Noether’s theorem. A steady solution to a system involving one Ripa-like layer and otherwise homogeneous layers can be proved formally (or Arnold) stable using the above invariants. A model configuration with only one layer has been previously shown to provide: a very good representation of the exact vertical normal modes up to the first internal mode; an exact representation of long-perturbation (free boundary) baroclinic instability; and a very reasonable representation of short-perturbation (classical Eady) baroclinic instability. Here it is shown that substantially more accurate overall results with respect to single-layer calculations can be achieved by considering a stack of only a few layers. Similar behavior is found in ageostrophic (classical Stone) baroclinic instability by describing accurately the dependence of the solutions on the Richardson number with only two layers.

Keywords: Shallow water equations; inhomogeneous layers; stratification; shear; mixed layer dynamics and thermodynamics

1.Introduction

1.1.Motivation

There is renewed interest to construct models for the study of the dynamics in the upper ocean (i.e., above the main thermocline, including the mixed layer) such that:

  1. are capable of incorporating thermodynamic processes while maintaining the two-dimensional structure of the rotating shallow-water equations, a paradigm of ocean dynamics on scales longer than a few hours 45; and

  2. preserve the geometric (generalized Hamiltonian) structure of the exact three-dimensional models from which they derive 31.

Property 1) promises to deliver a fundamental understanding of ocean processes that is difficult-if not impossible-to be attained using ocean general circulation models. Property 2) enables applying a recent flow-topology-preserving framework 28 to build parametrizations 17 of unresolvable submesoscale motions and this way investigating the contribution of these to transport at resolvable scales, a topic of active research 39.

1.2.Background

Back in the late 1960s and early 1970s and independently by various authors 20 33 42, the rotating shallow-water model was extended by allowing for horizontal and temporal variations of the density field while keeping it as well as the velocity field independent of depth. In the simplest setting, e.g., with one active layer floating atop an abyssal layer of inert fluid, the resulting inhomogeneous-layer model enables the investigation of thermodynamic processes in the upper ocean driven by heat and freshwater fluxes across the surface. Due to the two-dimensional nature of the model, the computational cost involved in such an investigation is considerably much lower than that produced by an ocean general circulation model 2 37.

Following the nomenclature introduced in 51, we will refer to the model above as IL0, which represents an inhomogenous-layer model wherein fields are not allowed to vary in the vertical. The homogeneous-layer shallow-water model will be called HL. Additional, more recent terminology for the IL0 is “thermal rotating shallow-water model” 66 71, which emphasizes the ability of the to include (horizontal) gradients of temperature. The is also being called the “Ripa model” in the literature 15 18 19 41 47 59, in recognition of Pedro Ripa’s contribution to its understanding 49-51 53 55. We will reserve that to refer to the model generalized here, which was introduced in Ref. 51.

The assessment on the computational cost efficiency of the IL0 holds even when more than one active layer is considered 34-36,60,70,72 or when the abyssal layer is activated and rests over irregular topography 6 7 44. Furthermore, due to the simplicity of the IL0 compared to the primitive equations for arbitrarily stratified fluid, referred to herein as IL, it has facilitated conceptual understanding of basic aspects of the upper-ocean dynamics and thermodynamics 13 54 56 57. Due in part to this very important reason, namely, the possibility to gain insight that is difficult to attain with an ocean general circulation model, the IL0 has been recently revisited 13,24,31,70.

A multilayer version of the IL0 was derived in 49, and a low-frequency approximation was developed in 53; cf. recent derivations in 64,65. The no-vertical-variation ansatz cannot be maintained under the exact dynamics produced by the IL when horizontal density gradients are present. The recipe used to keep the dynamical field’s depth-independent is to vertically average the horizontal pressure gradient. (Some authors 22 postulate a turbulent momentum flux that exactly cancels the vertical variation of the horizontal pressure gradient, but this simply is an ad-hoc hypothesis which does not contribute to the understanding of the problem.) While this is an approximation, 49 showed that it does not spoil the integrals of motion and generalized Hamiltonian structure of the problem.

Furthermore, the IL0 possesses a Lie-Poisson Hamiltonian structure (18) and associated with it an Euler-Poincare variational formulation (16) wherein the Hamilton principle’s Lagrangian follows by vertically averaging that of the IL29. When the equations of motion are derived in this formulation, there is a natural way to express three fundamental relations 31. These are 1) the Kelvin circulation theorem, 2) the advection equation for potential vorticity, and 3) an infinite family of conserved Casimir invariants (arising from Noether’s theorem for the symmetry of Eulerian fluid quantities under Lagrangian particle relabelling). The Euler-Poincare formulation provides a means to consistently introduce data-driven parameterizations of stochastic transport using the SALT (stochastic advection by Lie transport) algorithm 29 29, enabling data assimilation in a geometry-preserving context.

The IL0 provides an attractive framework for applying the SALT algorithm to derive parameterizations for unresolved submesoscale motions in the upper ocean. Indeed, numerical simulations of the IL024 43 46 tend to reveal small-scale circulations that resemble quite well 40 submesoscale filament rollups often observed in satellite-derived ocean color images. Such submesoscale motions may be unresolvable in many computational simulations. The extent to which they contribute to fluid transport at resolvable scales is a subject of active investigation 39 that the SALT stochastic version of the IL0 may cast light on.

1.3.Limitations of the IL0

Despite the above geometric properties of the IL0, it has several less attractive aspects, which can be consequential for the production of small-scale circulations in the model. Discussed in detail by (55), these include:

  1. In addition to the classical Poincare and Rossby waves, the IL0 represents variations of the thickness and density that do not change the vertical average of the pressure gradient 51 52. This mode is not present in the IL.

  2. A uniform flow may be unstable 22 52 69. A priori, this phenomenon seems to be something different than baroclinic instability. For instance, unlike Eady’s problem, it experiences an “ultraviolet divergence” in the sense that a short-wave cutoff is lacking.

  3. Since the dynamical fields are kept depth-independent within each layer, there is no explicit representation of the thermal wind balance between the velocity vertical shear and the horizontal density gradient, which dominates at low frequency.

An important additional limitation imposed by the depth independence of the dynamical fields, and particularly the buoyancy, is:

  • 4.The IL0 cannot represent the restratification of the oceanic surface mixed layer resulting from ageostrophic baroclinic instability of lateral density gradients, which tend to slump from the horizontal to the vertical 14 25 65.

1.4.The IL1

To cure the unwanted features of the IL0, Ref. 51 proposed the following improved closure to incorporate thermodynamic processes in a one-layer ocean model not restricted to low frequencies:

in addition to allowing arbitrary velocity and buoyancy variations in horizontal position and time, the velocity and buoyancy fields are also allowed to vary linearly with depth.

Ripa’s single-layer model, denoted IL1, enjoys many properties which make it very promising. For instance:

  1. The IL1 explicitly represents the thermal wind balance which dominates at low frequency.

  2. The free waves supported by the IL1 (Poincaré, Rossby, midlatitude coastal Kelvin, equatorial, etc.) are a very good approximation to the first and second vertical modes in the exact model with an unlimited vertical variation.

  3. The IL1 provides an exact representation of long-perturbation baroclinic instability and a very reasonable representation of short-perturbation baroclinic instability.

1.5.This paper

In this paper; I present a generalization of the IL1 to an arbitrary number of layers, including two possible (mathematically equivalent) vertical configurations (Sec.2). The model obtained incorporates additional flexibility to treat more complicated problems than those that can be tackled with only one layer. With a single layer in a reduced-gravity setting, mixed-layer processes can be minimally modeled. Including additional layers can lead to a more accurate representation of such processes. On the other hand, considering a stack of several layers atop an irregular bottom will enable investigating the influence of the ocean’s interior and even topographic effects. Several aspects of the generalized IL1 are discussed in Sec.3. These include: remarks on submodels derived from the generalized model as special cases (Subsec. 3.1); the nature of the layer boundaries (Subsec. 3.2); the model conservation laws (Subsec. 3.3); a discussion on circulation theorems (Subsec. 3.4); a formulation of the model suggestive of a generalized Hamiltonian structure (Subsec. 3.5); a formal stability theorem (Subsec. 3.6); results on vertical normal modes (Subsec. 3.7) and baroclinic instability (Subsec. 3.8), both quasigeostrophic and ageostrophic, which demonstrate that improved performance with respect to the single-layer results can be attained by considering only a few more layers; and the incorporation of forcing in the model equations (Subsec. 3.9). Section 4 closes the paper with some concluding remarks.

2.The multilayer IL1

Consider a stack of n active fluid layers with thickness hi(x,t), i = 1,…,n, where x is the horizontal position and t stands for time (Fig. 1). The geometry can be either planar or spherical; in the former case, the vertical coordinate, z, is perpendicular to the plane, whereas, in the latter, it is radial. The total thickness is h(x,t) = ∑ j h j (x,t). The stack of inhomogeneous-density layers can be either limited from below by a rigid bottom, z = h0(x), or from above by a rigid lid, z = -h0(x). The usual choice in the rigid lid case is h00; however, laboratory experiments are often designed to have a nonhorizontal top lid. The remaining boundary in the rigid-bottom (resp., rigid-lid) configuration is a soft interface with a passive, infinitely thick layer of lighter (resp., denser) homogeneous fluid of density ρn+1. Although vacuum (ρn+10) is the typical setting in the rigid-bottom configuration, the choice ρn+10 can be useful to study of deep flows over topography.

Figure 1 The two possible vertical configurations of the n-IL1 are rigid bottom (a) and rigid lid (b). Within each layer, the velocity and buoyancy fields not only vary arbitrarily with the horizontalposition and time but also linearly with depth. 

A key element to generalize Ripa’s model is to define a scaled vertical coordinate σ i that varies linearly from ±1 at the base to 1 at the top of the ith layer (Fig.2):

±z=:h~i-1(x,t)+1-σi2hi(x,t)=νi(x,σi,t), (1)

where

h~i(x,t):=h0(x)+j=1ihj(x,t) (2)

[henceforth, an upper (resp., lower) sign will correspond to the rigid-bottom (resp., rigid-lid) configuration]. The scaled vertical coordinate σ defined in Ref. 51 according to

z=:h0(x)+12(σ-1)h(x,t)=ν(x,σ,t) (3)

relates to the ith-layer scaled vertical coordinate σi defined here through

σ=1-2j=1i-1hjh+(1-σi)hih. (4)

Figure 2 Vertical coordinate choice. Within each layer, the rescaled vertical coordinate σ varies linearly from ±1, at the base, to 1, at the top. The upper (resp., lower) sign corresponds to the rigid-bottom (resp., rigid-lid) configuration of Fig. 1

Let an overbar denote vertical average within the ith layer:

a¯i(x,t):=12-1+1a(x,σ,t)dσi=12-1+1ai(x,σi,t)dσi. (5)

Following Ref. 51 closely, the ith-layer horizontal velocity and buoyancy fields are written, respectively, as

ui(x,σi,t)=u¯i(x,t)+σiuiσ(x,t), (6a)

ϑi(x,σi,t)=ϑ¯i(x,t)+σiϑiσ(x,t), (6b)

which can be regarded as a truncation of an expansion in orthogonal polynomials of σi of the form

ai(x,σi,t)=a¯i(x,t)+σiaiσ(x,t)+12(σi2-13)aiσσ(x,t)+16(σi3-35σi)aiσσσ(x,t)+, (7)

where aiσ:=σiai¯, aiσσ:=σiσiai¯, etc. 55. The ith-layer buoyancy is defined as

ϑi(x,σi,t):=±gρi(x,σi,t)-ρn+1ρr, (8)

where the upper (resp., lower) sign corresponds to the rigid-bottom (resp., rigid-lid). Here, g is gravity, ρi(x,σi,t)=ρ¯i(x,t)+σiρiσ(x,t) is the (variable) density in the ith layer, and ρr denotes the (constant) reference density used in the Boussinesq approximation. Physically admissible buoyancy values, i.e., everywhere positive and monotonically increasing (resp., decreasing) with depth in the rigid-bottom (resp., rig-lid) case, are such that

ϑ¯i>ϑiσ>0,ϑ¯i-ϑ¯i+1ϑiσ+ϑi+1σ. (9)

If ni2(x,t)>0 is the square of the instantaneous Brunt-Väisälä frequency within the ith layer, then note that

ϑiσ=1/2ni2hi. (10)

In order to obtain the equations for the n-layer version of Ripa’s, model one must proceed as follows:

  • 1.Substitute ansatz (6) in the inviscid, unforced, primitive equations (namely, rotating, incompressible, hydrostatic, Euler-Boussinesq equations) for arbitrarily stratified fluid (IL), which can be written as

Dϑ=0, (11a)

t|σh+|σhu+hσμ=0, (11b)

Du+fz^×u+|σp+ϑ|σν=0, (11c)

σp-12hϑ=0, (11d)

where

μ:=Dσ=2Dh0+1-σDh2wh. (11e)

In (11),

D:=t|σ+u|σ+μσ (11f)

is the material derivative, where tσ and σ indicate, respectively, that the partial time derivative and the horizontal gradient operate at constant σ [note that tσata and σaa, and thus Data+ua, for any a(x,t) ]; f is the Coriolis parameter (twice the local angular rotation frequency) and z^ is the vertical unit vector. Also, in Eq. (11), (u,w) is the three-dimensional velocity, μ denotes the σ-vertical velocity, ϑ stands for buoyancy, and p is a kinematic pressure; the vertical variation in all these fields is unrestricted. Equations (11a-d) are defined in -1<σ<+1 (i.e., h0<±z<h0+h) and are subject to the boundary conditions

μ=0  at σ={-1,+1}, (11g)

p=0  at  σ=-1. (11h)

Note that boundary conditions (11g) can be expressed as (t|σ+u|σ)(h0+(1/2)[11]hζ)=0 at the base of the layer and (t|σ+u|σ)(h0+(1/2)[1±1])hζ)=0 at the top of the layer. Here, ζ(x,σ,t) is the vertical displacement of a constant-density surface or isopycnal, which, by virtue (11a), relates to the vertical velocity through w=Dζ. These conditions thus indicate that a fluid particle initially on a given boundary remains there at all times, conserving its density. A particular case is one in which all particles on the boundary have the same density, i.e., h0+12[11]hζ=const at the base of the layer and/or h0+12[1±1]hζ=const at the top of the layer.

  • 2.Replace all occurrences of σ2 by its vertical average (i.e., σ2(1/3)) to preserve the linear vertical structure within each layer.

  • 3.Collect terms in powers of σ and equate them to zero afterward.

The equations that result from the above three-step procedure constitute the n-IL2 and are given by:

Diϑi¯=0, (12a)

(Diϑi)σ=0, (12b)

thi+hiu¯i =0,(12c)

Diui¯+fz^×u¯i+p¯i =0, (12d)

(Diui)σ+fz^×uiσ+(pi)σ=0. (12e)

Here,

Diai¯=ta¯i+u¯ia¯i+13hi-1hiaiσuiσ, (12f)

(Diai)σ=taiσ+u¯iaiσ+uiσa¯i, (12g)

are the mean and σ components of the material derivative of any field ai(x,σi,t)=a¯i(x,t)+σiaiσ(x,t) in the ith layer; and

p¯i=(ϑ¯i-13ϑiσ)hi+12hi(ϑ¯i-13ϑiσ)+ϑ¯ih~i-1+j=i+1nhjϑ¯j, (12h)

(pi)σ=12ϑiσhi+12hiϑ¯i+ϑiσh~i-1, (12i)

which are the mean and σ components of the ith-layer pressure gradient force.

System (12) consists of 7n evolution equations in the 7n independent fields (ϑ¯i,ϑiσ,hi,u¯i,uiσ), i = 1,…n. The coupling among different layer quantities is provided by the last terms on the right-hand side of the pressure forces (12h,i). It is important to note that the dynamics in both the rigid-bottom and rigid-lid configurations are described by the system (12); no double signs are needed. The latter must be taken into account, however, in the computation of the total pressure in the ith layer, which, up to the addition of an irrelevant constant, is given by ρrpi±ρn+1gz, where

pi=1/21+σihiϑ¯i-1/4(1-σi2)hiϑiσ+j=i+1nhjϑ¯j. (13)

Finally, (12) are satisfied in some closed but multiply-connected horizontal domain, say D. On D,i.e., the union of each disconnected part of the solid boundary of D, the zero normal flow condition holds:

u¯in^=0=uiσn^  on D (14)

where n^ is normal to D.

3.Discussion of several aspects of the n-IL1

3.1.Submodels

Any initial state with uniform buoyancy inside each layer (ϑ¯i=const and ϑiσ0) and vanishing vertical shear (uiσ0) is readily seen to be preserved by (12); consequently, the n-HL (a model with n homogeneous layers) follows from (12) as a particular case, just as it does it from the (exact, three-dimensional) IL model (11). In other words, the n-HL evolves on an invariant submanifold of both the n-IL1 and IL. Noteworthy, the n-HL is exact for a stepwise density stratification; however, as mentioned above, it is not able to accommodate thermodynamic processes, e.g., due to heat and buoyancy fluxes across the ocean surface. The n-IL0 developed in 49 follows from (12) upon neglecting uiσ and ϑiσ; note that an initial condition with uiσ0 and ϑiσ0 is preserved neither by (12) nor by (11), so the n- IL0 is not a particular solution of neither the n-IL1 nor the IL. Ignoring uiσ in (12) results in a model with uiσ0 but ϑiσ0, which provides a generalization for Schopf and Cane’s [1983] intermediate layer model. Alternatively, the omission of ϑiσ in the system (12) gives a model with ϑiσ0 but uiσ0. This model differs from earlier related models 8 61 68 in that it is not restricted to low-frequency motions and that it explicitly represents vertical shear within each of an arbitrary number of layers.

3.2.Layer boundaries

Consistent with ansatz (6) and the assumption of zero mass transport across layer boundaries, the σ-vertical velocity (11e) in the ith layer reads

μi=1-σi22hihiuiσ, (15)

which vanishes at the base and the top of the layer. Consequently, (t+[u¯iuiσ])(h~i-1+(1/2)[11]hi[ζ¯iζiσ])=0 at the base of the ith layer and (t+[u¯i±uiσ])(h~i-1+(1/2)[1±1]hi[ζ¯i±ζiσ])=0 at the top of the ith layer. Namely, the layer boundaries (interfaces and rigid bottom or lid) of the n- are material surfaces on which each fluid particle retains its density. This includes the particular situation in which all fluid particles on these boundaries have the same density, i.e., h̃i-1+1211hi[ζ¯iζiσ]=const at the base of the ith layer and h̃i-1+121±1hi[ζ¯i±ζiσ]=const at the top of the ith layer. The latter situation, which is most likely to happen far away from the ocean surface, cannot be described by the with only one layer.

3.3.Conservation laws

In a closed horizontal domain, on whose boundary (14) are satisfied, conservation of the ith-layer volume, mass, and buoyancy variance follows, respectively, because of (12c),

t(hiϑ¯i)+hi(ϑ¯iu¯i+13ϑiσuiσ)=0, (16)

and

t(hiϑi2¯)+hi(ϑi2¯u¯i+23ϑ¯iϑiσuiσ)=0. (17)

The total energy (sum of the energies in each layer) is also preserved in a closed horizontal domain since

tjEj+jhj(b¯ju¯j+13bjσujσ)=0, (18a)

where

Ei:=12hiu¯i2+16hi(uiσ)2\na+12hi2(ϑ¯i-13ϑiσ)+hih~i-1ϑ¯i, (18b)

And

b¯i:=12u¯i2+16(uiσ)2+hi(ϑ¯i-13ϑiσ)+h~i-1ϑ¯i+j=i+1nhjϑ¯j, (18c)

biσ:=u¯iuiσ+(h~i-1+12hi)ϑiσ, (18d)

which are the mean and σ components of the ith-layer Bernoulli head. The above result follows upon realizing that j=1nhjϑ¯jth~j-1-j=1nthjk=j+1nhkϑ¯k0, and is largely facilitated by rewriting (12d,e) in the form

tu¯i+μ¯iuiσ+hiz^×(q¯iu¯i+13qiσuiσ)+b¯i=R¯i, (19a)

tuiσ+hiz^×(qiσu¯i+q¯iuiσ)+biσ=Riσ. (19b)

Here,

μ¯i=13hi-1hiuiσ (20)

is the vertical average of the ith-layer σ-vertical velocity (12);

q¯i:=hi-1(f+u¯i×z^),qiσ:=hi-1uiσ×z^ (21)

are the mean and σ components of the ith-layer σ-potential vorticity;a and

R¯i:=h̃i-1ϑ¯i+1/2hiϑ¯i-1/3ϑiσ,  (22a)

Riσ:=(h̃i-1+1/2hi)ϑiσ-1/2hiϑ¯i, (22b)

which are rotational forces that arise as a consequence of the buoyancy inhomogeneities within each layer (ϑ¯i0ϑiσ).

In turn, the local conservation law for the sum of the zonal momenta within each layer is given by

tjMj+jFjM+xh0jhjϑ¯j=0, (23a)

where

FiM:=Miu¯i+13hiuiσuiσ+12γhi2(ϑ¯i-13ϑiσ)x^+γhi+1ϑ¯i+1j=1i-1hjx^ (23b)

with ui denoting the zonal component of ui and x^ the unit vector in the same direction iii . The above result follows upon multiplying by γhi the zonal component of (12d),

tu¯i+u¯iu¯i+13hi-1hiuiσuiσ-(f+τu¯i)v¯i+γ-1xp¯i=0, (24)

and realizing that

jhjϑ¯jx(h~j-1-h0)+jhjxk=j+1nhkϑ¯kxjhj+1ϑ¯j+1k=1j-1hk.

At this point, it is crucial to specify whether the geometry is flat or spherical. On the sphere, a=γ-1xa,ya, for any scalar a(x), and a=γ-1[xa+yγb], for any vector a = (a,b), where x=(λ-λ0)cosθR and y=(θ-θ0)R are, respectively, rescaled geographic longitude and latitude on the surface of the Earth whose mean radius is R; and γ(y):=cosθ0cosθ and τ(y):=R-1tanθ-γ-1dγ/dy are coefficients that characterize the geometry of the space (the arclength element square and area element are dx2=γ2dx2+dy2 and d2x=γdxdy, respectively). The ith zonal momentum (angular momentum around the Earth’s axis) is then given by

Mi:=hi[γu¯i-ΩR(cosϑ0-γcosϑ)], (25)

where Ω is the Earth’s angular rotation rate. In the classical β plane, γ = 1 and τ = 0 so that all the operators are Cartesian and Mi=hi(u¯i-f0y-(1/2)βy2). However, the geometry in a consistent β plane cannot be Cartesian; instead γ=1-τ0y, τ=τ0/γ, and Mi=hi[γu¯i-f0y-12β1-R2τ02y2] 54. Finally, conservation of the total zonal momentum (sum over all layers) in a horizontal domain in addition requires, in all cases, that both the topography and coasts be zonally symmetric.

3.4.Circulation theorems

In the IL, the circulation of u+uf, where z^×uf:=f, around a material loop, is constant in time if the latter is chosen to lie on an isopycnic surface ii . This is known as the Kelvin circulation theorem, which via Stokes’ theorem implies conservation of ϑ-potential vorticity. From the Hamiltonian mechanic’s side, the Kelvin theorem is the geometrical statement of invariance of the fluid action integral on level surfaces of ϑ 27. The existence of a Kelvin circulation property is thus closely related to the existence of a (constrained) Hamilton’s principle for the IL. The n-IL1 does not hold such a circulation property. As a consequence, the evolution of the ith-layer ϑ -potential vorticity is not correctly represented. In Ref. 51 it is shown that this is the result of the lack of information on the vertical curvature of the horizontal velocity field. It is easy to show, however, that the evolution of the three components of the vorticity field is correctly represented, and, consistent with the IL, neither q¯i nor qiσ are conserved. The evolution equations of the latter fields and the horizontal vorticity are given by (4.21) in Ref. 51 evaluated in the ith layer (note that evaluation of v in the ith layer does not simply mean replacing h by h i). The nonexistence of a Kelvin circulation property for the n-IL1 suggests that finding Hamilton’s principle for it is, at least, nontrivial. The n-IL1 is nonetheless shown in Sec. 3.5 to admit a formulation suggestive of a generalized Hamiltonian structure. The n-IL0, surprisingly, possesses a Kelvin circulation property since

\d\dtlt(u¯i)(u¯i+uf)dx=lt(u¯i)(h~i-1ϑ¯i+12hiϑ¯i)dx

holds in that model, and the material loop lt(u¯i) can be chosen to lie on an isopycnic surface. Consistent with the presence of this property, in 18 it is shown that the IL0 has a Lie-Poisson Hamiltonian structure which implies an analogous Euler-Poincare variational formulation 31 and, hence, the existence of a Lagrangian functional.

In the n-IL0, the following circulations theorems hold:

ddtDu¯idx=D(R¯i-μ¯iuiσ)dx,

and

ddtDuiσdx=DRiσdx.

This contrasts with the IL for which the circulation of u around D is time-independent. Note that the circulation of uiσ around D would be invariant if both ϑ¯i and ϑiσ were chosen such that n^×ϑ¯i=0=n^×ϑiσ on D iv . However, the latter boundary is not preserved by the n-dynamics. In opposition, the condition n^×ϑ¯i=0 on D is preserved by the n-IL0 dynamics, thereby guaranteeing invariance of the circulation of u¯i around D. This has been shown 49 to have important consequences for the generalized Hamiltonian structure of the IL0.

3.5.A formulation suggestive of a generalized Hamiltonian structure

The Euler equations of fluid mechanics possess what is called a generalized Hamiltonian structure 40. The IL (11), which derives from the Euler equations, is also Hamiltonian in a generalized sense 1. A good sign of the validity of any approximate model derived from IL the is the preservation of the generalized Hamiltonian structure. This section is devoted to showing that the n-IL1 -admits a formulation suggestive of a generalized Hamiltonian structure. A stronger statement was made in Ref. 51 for 1- IL1.

Let φ(x,t)=(φ1(x,t),,φ7n(x,t)) be a point on the infinite-dimensional phase space with coordinates (ϑ¯i,ϑiσ,hi,u¯i,uiσ), i = 1,…,n. Consider the relevant class, say A, of sufficiently smooth real-valued functionals of φ. For any phase functional F[φ]A it is further assumed that its density does not depend explicitly on t, namely, F[φ]=DF(φ,φ,,x)d2x, and that it satisfies the boundary conditions v

δFδu¯in^=0=δFδuiσn^  onD. (25)

A phase functional F[φ]A will be said to be admissible. Introduce then the functional

H[φ]:=jEj, (27)

where

j:=Dd2xj (28)

and E i is the energy in the ith layer (18b); its functional derivatives are given by

δHδϑ¯i=hi(h~i-1+12hi),δHδϑiσ=-hi26,δHδhi=hib¯i,δHδu¯i=hiu¯i,δHδuiσ=hiuiσ3. (29)

The latter and the zero normal flow conditions across D (14) show that H is admissible. Let now

J=jIj+Kj (30)

be a skew-adjoint 7 x 7 block-diagonal matrix operator where

Ii=-(0000000000000()000()q¯iz^×()qiσz^×()000qiσz^×()3q¯iz^×()), (31a)

Ki=-(000hi-1()ϑ¯ihi-1ϑiσ()000hi-1()ϑiσ3hi-1()ϑ¯i00000-hi-1()ϑ¯i-hi-1()ϑiσ00hi-1\uiσ()ϑiσ(hi-1)-3hi-1()ϑ¯i0(hi-1uiσ)0). (31b)

Here, the circle (resp., bullet) in parenthesis indicates operation on a scalar (resp., two-component vector). Define further a bracket operation {,}:A×AA as

{F,G}:=DδFδφJδGδφd2x (32)

F,G[φ]A. Then the layer model (12) can be written in the form

tφ={φ,H}, (33)

which is equivalent to F˙={F,H}F[φ]A.

The bracket operator (32) satisfies {F,G}=-{G,F} (anticommutativity), {F,aG+bK}=a{F,G}+b{F,K} (bilinearity), and {FG,K}=F{G,K}+G{F,K} (Leibniz’ rule), where a, b are arbitray numbers and F,G,K[φ] are any admissible functionals. The anticomutativity property follows from the skew-adjointness of the matrix operator J [boundary terms cancel out by (26)]. The bilinearity property and Leibniz’ rule are direct consequences of the bracket’s definition.

That system (12) can be cast in the form (33) appears to suggest that the n-IL1 is Hamiltonian in a generalized sense, with the functional H and the bracket operator {,} being the Hamiltonian and Poisson bracket, respectively. However, the bracket (32) does not seem to qualify as Poisson since {{F,G},K}+{{G,K},F}+{{K,F},G}=0 (Jacobi’s identity) does not seem to hold.

In addition to independence of the choice of phase space coordinates, the Hamiltonian structure conveys other important properties like the direct linkage of conservation laws with symmetries via Noether’s theorem 62. While the n-IL1 cannot be formally proved to be Hamiltonian, its energy, H, and -M, where M[φ]:=jMj is the zonal momentum of the system, do appear to be generators of t- and x-translations because of (33) and xφ={M,φ}, respectively. The latter assumes that M is an admissible functional, which requires the horizontal domain to be x-symmetric since (δM/δuiσ)0 and (δM/δv¯i)0, but (δM/δu¯i)=γhi0. Then δHH=ε{H,H}=εH0 for the infinitesimal variation δHφ:=ε{φ,H}=εtφ induced by H, and δMH=ε{H,M}=-εM˙=-εjhjϑ¯jxh00 iff xh00 for the infinitesimal variation δMφ:=ε{φ,M}=-εxφ induced by M. Consequently, conservation of H and M are linked, respectively, to t- and x-symmetries of H (horizontal domain and topography in this case included).

A distinguishing feature of generalized Hamiltonian systems is the existence of Casimirs C[φ]A, satisfying {C,F}0 F[φ]A. The Casimirs are thus integrals of motion, yet not related to (explicit) symmetries because {φ,C}0 ( C does not generate any transformation). The ith-layer integrals of volume, mass, and buoyancy variance are all admissible functionals that commute with any admissible function in the bracket in Eq. (32). The n-IL1-does not seem to support additional “Casimir” invariants.

The possibility of deriving a stochastic n-using the SALT approach 28 is constrained to the existence of a Kelvin circulation theorem, which is lacking for the n-IL1. The lack of a Kelvin circulation theorem is tied to the nonexistence of a generalized Hamiltonian structure and associated Euler-Poincare variational formulation for the n-IL1. While building parameterizations of unresolved submesoscale motions does not seem plausible using this flow-topolgy-preserving framework, investigating the contribution of the submesoscale motions to transport at mesoscales is still possible via direct numerical simulation. For this, the apparent generalized Hamiltonian formulation of the n-IL1 can be helpful, as finite-difference schemes that preserve the conservation laws of the system might be sought using the bracket approach developed in 58.

3.6.Arnold stability

In Ref. 51 it was shown that a state of rest (or a steady-state with at most a uniform zonal current) in the 1-IL1 can be shown to be formally stable using Arnold’s [1965; 1966] method if and only if Eq. (9) is satisfied, i.e., if and only if the buoyancy is everywhere positive and increases (resp., decreases) with depth within a layer with the rigid bottom (resp., rigid lid). Arnold’s method for proving the stability of a steady solution of a system consists of searching for conditions that guarantee the sign-definiteness of a general invariant which is quadratic to the lowest-order in the deviation from that state; the resulting conditions are only sufficient 26 38. In the n-IL1 with n>1, however, Arnold’s method fails to provide stable conditions even for a state of rest and with no topography (h00). The lowest-order (quadratic) contribution to that invariant, which can be called a “free energy” because it is defined with respect to a state of rest,

E:=12jHj(δu¯j)2+13Hj(δujσ)2+(gj-12Nj2Hj)(δhj)2+Nj-2Hj(δϑ¯j+12Nj2δhj)2+13Nj-2Hj(δϑjσ-12Nj2δhj)2+(gjδhj+Hjδϑ¯j)δh~j-1, (34)

cannot be proved sign-definite when n>1. Here, Hi, gi and Ni are the ith-layer unperturbed depth, vertically averaged buoyancy, and Brunt-Väisälä frequency, respectively. Similarly, a state of rest in the n-for any n cannot be proved formally stable using Arnold’s method. Surprisingly, it is possible to prove the stability of a steady-state with a uniform zonal current in that model. But the condition of stability is not one of “static” stability like Eq. (9) as in the 1-IL1. Contrarily, it is one of “baroclinic” stability since a uniform current in the n-IL0 has an implicit vertical shear through the thermal-wind balance. These results can all be inferred from 49 and 52.

Nevertheless, there is at least a system, which has one IL1-like layer and n -1 IL1-like layers, for which a state of rest can be proved formally stable. For instance, choosing the uppermost layer to be IL0 -like, the corresponding free energy takes the form

E:=12jHj(δu¯j)2+13Hα(δuασ)2+12Nα-2Hα(δϑ¯α+12Nα2δhα)2\na+13Nα-2Hα(δϑασ-12Nα2δhα)2+(gj-gj+1)(δh~j)2-12Nα2Hα(δhα)2, (35)

where α:=n (resp., α:=1) for the rigid-bottom (resp., rigid-lid) configuration, and Hi, gi, and N i are all constants. The above free energy is positive-definite if and only if (9) if fulfilled. [The n-HL has an infinite set of invariants which are given by jhjF(q¯j) where F() is arbitrary; these include the volume integral, which is the only one needed to obtain the above result.] When all layers are homogeneous, the same result is obtained. When one IL0-like layer is included, however, the free energy cannot be shown of one sign.

That a steady-state (with or without a current) of the n IL0-cannot be proved formally stable does not mean that such a state is unstable; it means that Arnold’s method is not useful to provide sufficient conditions for the stability of that state.

3.7.Waves

The n-IL1 Eqs. (12), linearized with respect to a reference state with no currents, can be shown to sustain the usual midlatitude and equatorial gravity and vortical waves (Poincaré, Kelvin, Rossby, Yanai, etc.) in 2n vertical normal modes. Here I shall concentrate on how well these modes are represented by considering the phase speed of (internal) long gravity waves assuming a rigid-lid setting.

The reference state is characterized by the parameter

S:=Nr2Hr2gr, (36)

which must be such that 0<S<1 11 51. Here, Nr is the reference Brunt-Väisälä frequency within an active layer floating on top of an inert layer; Hr is the total thickness of the active fluid layer, and gr denotes the vertically averaged reference buoyancy within the active layer. All three reference quantities are held constant. The reference buoyancy then varies linearly from gr(1+S) at the top of the active layer to gr(1-S) at the base of the active layer. In Ref. 51, it was shown that the 1-IL1 gives the exact result for the “equivalent” barotropic or external mode phase speed of (internal) long gravity waves for all and a very good approximation to the first internal mode phase speed for all S.

Figure 3 compares, as a function of S, the phase speed as determined by the IL, n-HL, n-IL0, and n-IL1 for various n. The figure shows the results for the external mode (c0), and the first (c1) and second (c2) internal modes. The analytical expression for the ’s phase speed for an arbitrary mode number can be found in Ref. 51; the phase speeds for the layer models are computed numerically. The solutions of the n-HL and n-IL0 coincide because ϑ¯i is constant for a normal mode in the n-IL1. These models can only support n vertical normal modes. In contrast, the n-sustains vertical normal modes up to the (n + 1)th internal mode.

Figure 3 Phase speed of (internal) long gravity waves as a function of the stratification strength in a reduced-gravity reference state with no currents. 

As noted above, the 1-IL1 result coincides with that of for the barotropic mode. To approximate well the exact solution, two HL- or IL0-like layers are needed. The first internal mode solution is very well approximated using two IL1-like layers. Four HL-like layers do not provide a similar degree of approximation. The second internal mode solution is reasonably approximated with two IL1-like layers. The distance between the exact solution and that produced using four HL-like layers is of the same order. However, in every case, the n-HL (or the n- IL0) overestimates the exact phase speeds.

3.8.Baroclinic instability

As one further test of the validity of the n-IL1, the problem of baroclinic instability, particularly upper-ocean baroclinic instability, is considered here. (A subset of the results presented here appeared in 10.) The behavior in both quasigeostrophic and ageostrophic regimes is explored. The n- IL1solutions are compared in all cases with the IL solutions. In some cases, comparisons are also made with n-HL and n- IL0 solutions. In the quasigeostrophic regime, analytical expressions exist for the IL solutions. Analytical or semianalytical formulas for the dispersion relations also exist in this regime for the 1-IL1 and models with one IL0-like or two HL-like layers. The rest of the solutions shown are computed numerically upon finite differencing the corresponding eigenvalue problems.

Upper-ocean baroclinic instability, e.g., above the ocean thermocline, is studied in 11 using the IL and the 1-IL1 in a reduced-gravity setting. A basic state with a parallel current U=U(z)x^ is considered in that work to lie in an infinite channel on the f plane, to have a uniform vertical shear, and to be in thermal-wind balance with the across-channel buoyancy gradient. The basic velocity is further set to vary (linearly) from U¯+Uσ at the top of the active layer to U¯-Uσ at the base of the active layer. Accordingly, the basic buoyancy field Θ(y,z) varies from gr(1-2fUσy/Hr+S) at the top of the active layer to gr(1-2fUσy/Hr-S) at the base of the active layer (y is the across-channel coordinate). A nonvanishing velocity at the base of the active layer implies that the latter has a linear y-slope given by gr-1fUσ-U¯/(1-S). This basic state is a steady solution of the IL to the lowest order in the Rossby number, Ro:=U¯/L|f|Uσ/L|f| where L is the relevant length scale, which is assumed to be an infinitesimal parameter. In the limit of weak stratification (S0), the horizontal scales

RE:=grHrf,RI:=NrHrf (37)

are well separated (RERI), and thus long and short normal-mode perturbations to this state can be identified. Under long small-Rossby-number normal-mode perturbations, the base of the active layer behaves as a free boundary. For short small-Rossby-number normal-mode perturbations, this interface is effectively rigid. When the vertical shear is assumed strong, U¯/UσO(S-1), the short-perturbation limit corresponds to the classical Eady problem of baroclinic instability, in whose case solutions are insensitive to U¯/Uσ.

The left panel of Fig. 4 shows the minimum along-channel wavenumber, k, for instability as a function of U¯/Uσ in the long-perturbation and strong-shear limits (free-boundary baroclinic instability). The 1-gives the exact result for all U¯/Uσ . To provide a close approximation to this result for all U¯/Uσ with the n-HL, a fairly large n (cir. 25) is needed. Note that the 1-predicts, incorrectly, stability for U¯/Uσ<0 (the vertical shear in this model is implicit through the thermal-wind relation).

Figure 4 (left panel) Minimum wavenumber for long-perturbation and strong-shear (i.e., free-boundary) baroclinic instability as a function of the slope of the lower interface in the basic state. (right panel) Growth rate of the most unstable perturbation as a function of the wavenumber in short-perturbation, strong-shear (i.e., classical Eady) baroclinic instability. 

The right panel of Fig. 4 depicts, as a function of the along-channel wavenumber k, the growth rate of the most unstable perturbation in the short-perturbation and strong-shear limits (classical baroclinic instability). The comparison of the maximum growth rate predicted by the 1-IL1 with the IL’s maximum growth rate is less satisfactory in this limit. However, and very importantly, a high wavenumber cutoff of baroclinic instability is present. The 1- IL0 model only gives the k = 0 value of the growth rates of this figure, and thus it cannot be used to describe this regime ( Nr0 in this model). Three -like layers are enough to approximate well the exact maximum growth rate for all k. To obtain a similar result using HL-like layers, at least six must be considered.

For the basic state considered above, the Richardson number

Ri:=NrzU2SRE22Ro2L2. (38)

In classical baroclinic instability, for which Ro,S0, L=RI2SRE and U¯/UσO(S-1), the well-known result Ri holds. In free-boundary baroclinic instability, for which Ro,S0, L=RE and U¯/UσO(S-1), Rican acquire any value because a proper way to achieve the S0 limit is to set S=O(Roν) for any v. Unlike quasigeostrophic baroclinic instability, ageostrophic baroclinic instability is characterized by a dependence of the solutions on Ri 63,64. This dependence is checked in the layer model by considering infinitesimal non-geostrophic normal-mode perturbations to the above basic state but with U¯Uσ=U, and assuming Ro=10-1 and L=RE.

The left panel of Fig. 5 shows, as a function of Ri, the maximum growth rate maxk{kIc} of the perturbation. The right panel of the figure shows, also as a function of Ri, the wavenumber, kmax, at which the latter value is attained. Shown for reference is an IL asymptotic solution, valid up to O(Ro3). The asymptotic formulas for maxk{Ikc} and kmax are those given in (4.27) and (4.28) of 63. The n-IL1 fares very well even with n = 1. A model with a single IL0-like layer, however, cannot describe this regime because of the dependence on Ri (for the 1-IL0S0). With two IL1-like layers, the maximum growth rates and corresponding wavenumbers at which they are achieved are in very close agreement with the IL predictions in the range of Ri values explored, which was much wider than that shown in Fig. 5. Note, however, that observations indicate that typical values of Ri in the upper ocean are close to unity 65.

Figure 5 (left panel) Maximum normal-mode perturbation growth rate for ageostrophic (classical Stone) baroclinic instability as a function of the Richardson number. (right panel) Wavenumber for maximum growth rate. 

3.9 Forcing

In Ref. 51, forcing (wind stress, interfacial drag, and buoyancy/heat input) was introduced in the 1-IL1 model equations in a way that was compatible with the conservation laws of energy, momentum, and mass/heat content. The same approach is adopted here to include, also, freshwater fluxes through the surface following the conservation law of salt content. The possibility for the exchange of fluid across the other interfaces is also considered.

Let τ(x,t) be wind stress acting at the surface of the ocean ( ρn+10 must be the setting in the rigid-bottom configuration and typically h00 in the rigid-lid one). Assume further that there is a friction force acting at the interface between contiguous layers. Introduction of these forces in Newton’s equations (12d,e) in the form

tu¯i+=δiατ/hα-ri(u¯i±uiσ), (39a)

tuiσ+=3δiατ/hα+3ri(u¯i±uiσ), (39b)

implies that the work done by the wind stress is proportional to the velocity at the top of the uppermost layer, u¯αuασ, and that one done by the friction force in the ith layer is proportional to the velocity at the base of that layer, u¯i±uiσ. Namely,

tjEj+=τ(u¯αuασ)-jrjhj(u¯j±ujσ)2, (40a)

tjMj+=τx^-jrjhj(u¯j±ujσ)x^. (40b)

In the above equations, δij is the Kronecker delta, and r i is a friction coefficient that can be taken as a constant or as some function of h i and |u¯i±uiσ|. [Recall that α:=n (resp., α:=1) for the rigid-bottom (resp., rigid-lid) configuration.]

Let now Γ(x,t) be a buoyancy input through the surface and write the buoyancy (12a,b) in the form

tϑ¯i+=δiαΓ/hα, (41a)

tϑiσ+=ηδiαΓ/hα, (41b)

where η is any constant. Consider, also, the possibility of fluid crossing the interface between consecutive layers; then the volume conservation Eq.(12c) can be rewritten as

thi+=wib-wit. (41c)

Here, the quantities wit(x,t) and wib(x,t) are volume fluxes per unit area through the top and base of the ith layer, respectively. The set (41), for any value of η, is compatible with the mass conservation equation

t(hiϑ¯i)+=δiαΓ+ϑ¯i(wib-wit). (42)

At the surface wαt(x,t)=E(x,t)-P(x,t), which represents the imbalance of evaporation minus precipitation vi . Away from the surface, some parametrization must be adopted. In models with IL0-like layers, it is commonly set 34

wib-wit=(-1)i+1(hi-1-Hi-1e)2Hi-1etieθ(Hi-1e-hi-i). (43)

Here, Hie and tie are constants with units of length and time, respectively, that characterize the “entrainment” process, and θ() is the Heaviside step function. In the present case, an algorithm may be designed such that condition (9) is fulfilled at all times. This would allow for a more natural representation of mixing processes, including the possibility of representing localized mixing events, e.g., characterized by ϑ¯i+1+ϑi+1σ<ϑ¯i-ϑiσ instantaneously at a certain position. This subject deserves to be studied in detail.

Let finally assume a linear state equation, i.e., ϑi=gαT(Ti-Tn+1)-gαS(Si-Sn+1). Here, αT and αS are the thermal expansion and salt contraction coefficients, respectively; Ti(x,σ,t)=T¯i(x,t)+σTiσ(x,t) and Si(x,σ,t)=S¯i(x,t)+σSiσx,t are the ith layer temperature and salinity, respectively; and Tn+1 and Sn+1 are the inactive layer (constant) temperature and salinity, respectively. Let also write the buoyancy input as

Γ=gαT(ρrCp)-1Q+gαSS¯α(P-E), (44)

Cp is the specific heat at constant pressure, and Q(x,t) is the heat input through the surface. Equation (42) can then be split into heat and salt content conservation equations, namely,

t(hiT¯i)+=δiα(ρrCp)-1Q+T¯i(wib-wit), (45a)

t(hiS¯i)+=δiαS¯α(E-P)+S¯i(wib-wit). (45b)

If fluid across the surface is allowed only, the choice (44) enforces, on the one hand 9,

ddtjhjS¯j0, (46a)

and, on the other 12,

ddtT=V-1D(ρrCp)-1Qd2x+(T¯α-T)(P-E), (46b)

where

V:=jhjDd2x

h is the total volume and

T:=V-1jhjT¯j

is the average temperature in V. Note that (46b), unlike the equation satisfied by

jhjT¯j,

is independent-as it should-of the choice of the origin of the temperature scale 67.

4.Concluding remarks

This paper describes a multilayer extension of the single-layer primitive-equation model for ocean dynamics and thermodynamics introduced in 51. Inside each layer, the velocity and buoyancy fields can vary not only arbitrarily in the horizontal position and time but also linearly with depth.

In the absence of external forcing and dissipation, the model conserves volume, mass, buoyancy variance, energy, and zonal momentum for zonally symmetric horizontal domains and topographies. Unlike models with depth-independent velocity and buoyancy fields within each layer, the model generalized here can represent the thermal wind balance explicitly at low frequency inside each layer. In this sense, the model of this paper has “better” physics than a model with depth-independent fields. For a fixed number of layers, the model of this paper can sustain one more vertical normal mode than the homogeneous-layer models, which, on the other hand, are not able to incorporate thermodynamic processes (e.g., due to heat and buoyancy fluxes across the air-sea interface or associated with localized vertical mixing events). In this other sense, the present model has “more” physics than a model with homogeneous layers. Last but not least, overall improved results in both quasigeostrophic (free-boundary and classical Eady) and ageostrophic (classical Stone) baroclinic instability with respect to the single-layer calculations are attained with the addition of a small number of layers.

The present generalization enriches Ripa’s single-layer model by providing it enough flexibility to approach problems for which a single-layer structure is too idealized. Configurations with a small number of layers are particularly useful for the insight they provide into physical processes. Configurations with more layers may provide the basis for an accurate numerical circulation model.

Finally, and returning to the motivation for revisiting the construction of models with reduced thermodynamics, the requirement on the two-dimensional structure of the models is satisfied by the model derived here. A different strategy than that taken here is needed to fulfill the requirement on the geometric structure of the models if the goal is to pursue flow-topology-preserving parameterizations of unresolved scales using the SALT (stochastic advection by Lie transport) framework 28 29. The desired result might follow from plugging Ripa’s ansatz in the Hamilton principle’s Lagrangian of the primitive equations for continuously stratified fluid. This is currently under investigation. A stochastic parameterization framework that can be applied to the model derived here is location uncertainty (LU) 48. Unlike SALT dynamics, which preserve Kelvin circulation, the LU framework conserves energy, so it can be immediately applied to the present model and is a natural fit to considering the parameterizations based on the extraction of available potential energy 5 21 23. Building stochastic parameterizations using the generalized Ripa’s model is left for future work.

Acknowledgements

A stimulating epistolary exchange with Darryl Holm provided the incentive to revisit this work and finish it.

References

1. Abarbanel H., Holm D., Marsden J., and Ratiu T. Philos. Trans. R. Soc. London, A 318 (1986) 349. [ Links ]

2. Anderson D. L. T., and McCreary J. P. J. Atmos. Sci. 42 (1985) 615. [ Links ]

3. Arnold V. I. (1965). Dokl. Akad. Nauk. SSSR 162, 975-978, engl. transl. Sov. Math. 6 (1965) 773. DOI:10.1007/978-3-642-31031-74 [ Links ]

4. Arnold V. I. Izv. Vyssh. Uchebn. Zaved Mat. 54 (1966) 3 engl. transl. Am. Math. Soc. Transl. Ser. 2 79 (1969) 267. [ Links ]

5. Bachman S., Fox-Kemper B., Taylor J., and Thomas L. Ocean Modelling 109 (2017) 72. https://doi.org/10.1016/j.ocemod.2016.12.003 [ Links ]

6. Beier E. J. Phys. Oceanogr. 27 (1997) 615. https://doi.org/10.1175/1520-0485(1997)027h0615:ANIOTAi2.0.CO;2 [ Links ]

7. Beier E., and Ripa P. J. Phys. Oceanogr. 29 (1999) 305. https://doi.org/10.1175/1520-0485(1999)029h0305:SGITNGi2.0.CO;2 [ Links ]

8. Benilov E. J. Fluid Mech. 251 (1993) 501. doi:10.1017/S0022112093003490 [ Links ]

9. Beron-Vera F. J., Ochoa J., and Ripa P. Ocean Modell. 1 (1999)111. https://doi.org/10.1016/s1463-5003(00)00003-2 [ Links ]

10. Beron-Vera F. J., Olascoaga M. J., and Zavala-Garay J. (2004)In ICTAM04 Abstract Book and CD-ROM Proceedings. ISBN 83-89687-01-1, IPPT PAN, Warsaw. [ Links ]

11. Beron-Vera F. J., and Ripa P. J. Fluid Mech. 352 (1997) 245. https://doi.org/10.1017/S0022112097007222 [ Links ]

12. Beron-Vera F. J., and Ripa P. J. Geophys. Res. 105 (2000) 11441. [ Links ]

13. Beron-Vera F. J., and Ripa P. J. Geophys. Res. 107 (2002) (C8), 10.1029/2000JC000769. https://doi.org/10.1029/2000JC900038 [ Links ]

14. Boccaletti G., Ferrari R., and Fox-Kemper B. J. Phys.Oceanogr. 37 (2007) 2228. https://doi.org/10.1175/JPO3101.1 [ Links ]

15. Britton J., and Xing Y. Journal of Scientific Computing 82 (2020) 2. https://doi.org/10.1007/s10915-020-01134-y [ Links ]

16. Bröcker J., et al. Mathematics of Planet Earth: A Primer. World Scientific (2018). [ Links ]

17. Cotter C., Crisan D., Holm D., Pan W., and Shevchenko I. J. Stat. Phys. 179 (2020) 1186. doi:10.3934/fods.2020010 [ Links ]

18. Dellar P. J. Phys. Fluids 15 (2003) 292. https://doi.org/10.1063/1.1530576 [ Links ]

19. Desveaux V., Zenk M., Berthon C., and Klingenberg C. Mathematics of Computation 85 (2015) 1. https://doi.org/10.1002/fld.4177 [ Links ]

20. Dronkers J. J. Hydrau. Div. 95 (1969) 44. [ Links ]

21. Fox-Kemper B., Ferrari R., and Hallberg R. Journal of Physical Oceanography 38 (2008) 1145. [ Links ]

22. Fukamachi Y., McCreary J. P., and Proehl J. A. J. Geophys. Res. 100 (1995) 2559. [ Links ]

23. Gent P. R., and Mcwilliams J. C. Journal of Physical Oceanography 20 (1990) 150. https://doi.org/10.1175/1520-0485(1990)020h0150:IMIOCMi2.0.CO;2 [ Links ]

24. Gouzien E., Lahaye N., Zeitlin V., and Dubos T. Physics of Fluids 29 (2017) 101702. https://doi.org/10.1063/1.4996981 [ Links ]

25. Haine T.W., and Marshall J. J. Phys. Oceanogr. 28 (1998) 634. https://doi.org/10.1175/1520-0485(1998)028h0634:GSABIOi2.0.CO;2 [ Links ]

26. Holm D., Marsden J., Ratiu T., and Weinstein A. Phys. Lett. A 98 (1983) 15. https://doi.org/10.1016/0375-9601(83)90534-0 [ Links ]

27. Holm D. D. Physica D 98 (1996) 379. https://doi.org/10.1016/0167-2789(96)00121-2 [ Links ]

28. Holm D. D. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences (2015) 471 20140963. [ Links ]

29. Holm D. D., and Luesink E. (2020). arXiv:1910.10627. [ Links ]

30. Holm D. D., Luesink E., and Pan W. (2020). arXiv:2006.05707. [ Links ]

31. Holm D. D., Marsden J. E., and Ratiu T. S. (2002). In Large-Scale Atmosphere-Ocean Dynamics II: Geometric Methods and Models (ed. J. Norbury and I. Roulstone), pp. 251-299. Cambridge University. [ Links ]

32. Lahaye N., Zeitlin V., and Dubos T. Ocean Modelling 153 (2020) 101673. https://doi.org/10.1016/j.ocemod.2020.101673 [ Links ]

33. Lavoie R. J. Atmos. Sci. 29 (1972) 1025. [ Links ]

34. McCreary J. P., Fukamachi Y., and Kundu P. J. Geophys. Res. 96 (1991) 2515. [ Links ]

35. McCreary J. P., et al. J. Geophys. Res. 106 (2001) 7139. [ Links ]

36. McCreary J. P., and Kundu P. J. Mar. Res. 46 (1988) 25. https://doi.org/10.1357/002224088785113711 [ Links ]

37. McCreary J. P., Zhang S. and Shetye S. R. J. Geophys. Res. 102 (1997) 15535. [ Links ]

38. McIntyre M., and Shepherd T. J. Fluid Mech. 181 (1987) 27. https://doi.org/10.1017/S0022112087002209 [ Links ]

39. McWilliams J. C. Proc R Soc A 472 (2016) 20160117. https://doi.org/10.1098/rspa.2016.0117 [ Links ]

40. Morrison P. J. (1982). In Mathematical Methods in Hydrodynamics and Integrability in Dynamical Systems (ed. M. Tabor and Y. Treve), pp. 13-46. Institute of Physics Conference Proceedings 88. [ Links ]

41. Mungkasi S., and Roberts S. G. J. Phys.: Conf. Ser. 693 (2016) 012011. [ Links ]

42. O’Brien J. J., and Reid R. O. J. Atmos. Sci. 24 (1967) 197. [ Links ]

43. Ochoa J. L., Sheinbaum J., and Pavía E. G. J. Geophys. Res. 103 (1998) 24869. [ Links ]

44. Palacios-Hernández E., Beier E., Lavín M. F., and Ripa P. J. Phys. Oceanogr. 32 (2002) 705. https://doi.org/10.1175/1520-0485(2002)032h0705:TEOTSVi2.0.CO;2 [ Links ]

45. Pedlosky J. Geophysical Fluid Dynamics, (2nd edn. Springer 1987). [ Links ]

46. Pinet R., and Pavía E. J. Fluid Mech. 416 (2000) 29. https://doi.org/10.1017/S0022112000008569 [ Links ]

47. Rehman A., Ali I., and Qamar S. Results in Physics 8 (2018) 104. [ Links ]

48. Resseguier V., Pan W., and Fox-Kemper B. Nonlinear Processes in Geophysics 27 (2020) 209. https://doi.org/10.5194/npg-27-209-2020 [ Links ]

49. Ripa P. Geophys. Astrophys. Fluid Dyn. 70 (1993) 85. https://doi.org/10.1080/03091929308203588 [ Links ]

50. Ripa P. In Modelling of Oceanic Vortices (ed. G. V. Heist 1994),pp. 151. [ Links ]

51. Ripa P. J. Fluid Mech. 303 (1995) 169. https://doi.org/10.1017/S0022112095004228 [ Links ]

52. Ripa P. J. Geophys. Res. C 101 (1996) 1233. https://doi.org/10.1029/95JC02899 [ Links ]

53. Ripa P. Rev. Mex. Fís. 42 (1996) 117. [ Links ]

54. Ripa P. J. Phys. Oceanogr. 27 (1997) 597. https://doi.org/10.1175/1520-0485(1997)027h0597:TAPEOTi2.0.CO;2 [ Links ]

55. Ripa P. Dyn. Atmos. Oceans 29 (1999) 1. https://doi.org/10.1016/S0377-0265(98)00056-6 [ Links ]

56. Ripa P. (2001). In Proceedings of the 13th Conference on Atmospheric and Oceanic Fluid Dynamics, pp. 1-4. American Meteorological Society. [ Links ]

57. Ripa P. (2003). In Nonlinear Processes in Geophysical Fluid Dynamics: A Tribute to the Scientific Work of Pedro Ripa (ed. O. U. Velasco-Fuentes, J. Ochoa and J. Sheinbaum), pp. 103-126. Kluwer. Plublished post mortem. [ Links ]

58. Salmon R. J. Atmos. Sci. 61 (2004) 2016. [ Links ]

59. Sánchez-Linares C., de Luna T. M., and Castro Díaz M. J. Applied Mathematics and Computation 272 (2016) 369. https://doi.org/10.1016/j.amc.2015.05.137 [ Links ]

60. Schopf P., and Cane M. J. Phys. Oceanogr. 13 (1983) 917.https://doi.org/10.1175/1520-0485(1983)013h0917:OEDMLPi2.0.CO;2 [ Links ]

61. Scott R. B., and Willmott A. J. Dyn. Atmos. Oce. 35 (2002) 389. [ Links ]

62. Shepherd T. G. Adv. Geophys. 32 (1990) 287. https://doi.org/10.1016/S0065-2687(08)60429-X [ Links ]

63. Stone P. H. J. Atmos. Sci. 23 (1966) 390. https://doi.org/10.1175/1520-0469(1966)023h0390:ONGBSi2.0.CO;2 [ Links ]

64. Stone P. H. J. Atmos. Sci. 27 (1970) 721. https://doi.org/10.1175/1520-0469(1970)027h0721:ONGBSPi2.0.CO;2 [ Links ]

65. Tandon A., and Garrett C. J. Phys. Oceanogr. 24 (1994) 1419. https://doi.org/10.1175/1520-0485(1994)024h1419:MLRDTAi2.0.CO;2 [ Links ]

66. Warnerford E. S., and Dellar P. J. J. Fluid Mech. 723 (2013)374. [ Links ]

67. Warren B. A. J. Geophys. Sci. 104 (1999) 7915. [ Links ]

68. Young W. R. J. Phys. Oceanogr. 24 (1994) 1812. [ Links ]

69. Young W. R., and Chen L.. J. Phys. Oceanogr. 25 (1995) 3172. [ Links ]

70. J. J. Zavala-Hidalgo, A. Pares-Sierra, and J. Ochoa Atmosfera 15 (2002) 81. [ Links ]

71. Zeitlin V. Geophysical fluid dynamics: understanding (almost) everything with rotating shallow water models. (Oxford University Press 2018). [ Links ]

72. E. Moreles, J. Zavala-Hidalgo, B. Martinez-Lopez, and A. Ruiz-Angulo, Influence of Stratification and Transport on the Shedding Process, Journal of Geophysical Research: Oceans 126 (2021) e2020JC016315. [ Links ]

Received: October 12, 2020; Accepted: December 02, 2020

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License