Stochastic mesoscale circulation dynamics in the thermal ocean

Holm, Darryl D.; Luesink, Erwin; Pan, Wei

doi:10.1063/5.0040026

In analogy with similar effects in adiabatic compressible fluid dynamics, the effects of buoyancy gradients on incompressible stratified flows are said to be “thermal.” The thermal rotating shallow water (TRSW) model equations contain three small nondimensional parameters. These are the Rossby number, the Froude number, and the buoyancy parameter. Asymptotic expansion of the TRSW model equations in these three small parameters leads to the deterministic thermal versions of the Salmon's L1 (TL1) model and the thermal quasi-geostrophic (TQG) model, upon expanding in the neighborhood of thermal quasi-geostrophic balance among the flow velocity and the gradients of free surface elevation and buoyancy. The linear instability of TQG at high wavenumber tends to create circulation at small scales. Such a high-wavenumber instability could be unresolvable in many computational simulations, but its presence at small scales may contribute significantly to fluid transport at resolvable scales. Sometimes, such effects are modeled via “stochastic backscatter of kinetic energy.” Here, we try another approach. Namely, we model “stochastic transport” in the hierarchy of models TRSW/TL1/TQG. The models are derived via the approach of stochastic advection by Lie transport (SALT) as obtained from a recently introduced stochastic version of the Euler–Poincaré variational principle. We also indicate the potential next steps for applying these models in uncertainty quantification and data assimilation of the rapid, high-wavenumber effects of buoyancy fronts at these three levels of description by using the data-driven stochastic parametrization algorithms derived previously using the SALT approach.

I. INTRODUCTION

In this paper, we are dealing with the thermal rotating shallow water (TRSW) equations, which can be regarded as the vertically averaged version of the primitive equations with a buoyancy variable Zeitlin (2018).

In the balanced 2D model hierarchy of TRSW, Salmon's L1 model (TL1), and thermal quasi-geostrophic model (TQG), we are investigating a certain stochastic model of potential vorticity dynamics as a basis for stochastic parametrization of the dynamical creation of unresolved degrees of freedom in computational simulations of upper ocean dynamics. Specifically, we have chosen the SALT (stochastic advection by Lie transport) algorithm introduced in Holm and Luesink (2019) and applied in Cotter et al. (2018; 2019) as our modeling approach. The SALT approach preserves the Kelvin circulation theorem and an infinite family of integral conservation laws. The goal of the SALT algorithm is to quantify the uncertainty in the process of upscaling, or coarse graining of either observed or synthetic data at fine scales, for use in computational simulations at coarser scales. The present work prepares us to take the next step from (ii) to (iii) in the well-known path of discovery in oceanography, weather prediction, and climate science, which is

driven by large datasets and new methods for its analysis;
informed by rigorous mathematical derivations and analyses of stochastic geophysical fluid equations;
quantified using computer simulations, evaluated for uncertainty, variability and model error;
optimized by cutting edge data assimilation techniques, then
compared with new observation datasets to determine what further analysis and improvements will be needed.

The objective in applying the SALT algorithm to coarse grained simulations is to answer the following question, enunciated in Cotter et al. (2018, 2019): “How can one use computationally simulated surrogate data at highly resolved scales, in combination with the mathematics of stochastic processes in nonlinear dynamical systems, to estimate and model the effects on the simulated variability at much coarser scales of the computationally unresolvable, small, rapid, scales of motion at the finer scales?” The present paper will lay the theoretical foundations for addressing this question in the 2D context of the thermal rotating shallow water (TRSW) model and its balanced thermal quasi-geostrophic (TQG) model. Our eventual goal is to apply the SALT algorithm to calibrate our stochastic models for assimilating data, e.g., from satellite observations of the cascade in the upper ocean dynamics of horizontally circulating structures to smaller scales, as shown in Fig. 1 below.

A. The TRSW equations

As shown in Fig. 2, the TRSW equations arise in a series of nested approximations leading from the 3D Euler equations for inhomogeneous, stratified, rotating incompressible fluids, first to the 3D Euler–Boussinesq equations for small stratification, then to the rotating thermal Green–Naghdi equations, which were derived and investigated in Holm and Luesink (2019). Upon further neglecting the nonhydrostatic pressure effects which are present in the Green–Nahgdi model, the TRSW model is obtained. The TRSW equations [which were first called the IL⁰ model in Ripa (1995)] and their quasi-geostrophic approximation, the TQG equations, comprise standard models of thermal effects in geophysical fluid dynamics (GFD), as reviewed, e.g., in Bröcker et al. (2018); Zeitlin (2018). The TQG equations are discussed in the GFD literature by Warneford and Dellar (2013), for example, following earlier work by Ripa (1993; 1995) and Ripa (1999). In fact, the TRSW and TQG equations and their high wavenumber instabilities have been rederived several times, as recounted in Zeitlin (2018). A multilayer extension of the shallow water model with stratification and shear can be found in Beron–Vera (2020), which includes a historical background on the developments of the thermal rotating shallow water model. In this paper, we will derive the SALT stochastic versions of TRSW and TQG, as well as TL1, which is a thermal version of an intermediate theory known in the GFD literature as Salmon's L1 model Salmon (1983). In deriving the SALT versions of these deterministic equations, we will follow the geometric approach of Holm (2015), which is based on Hamilton's variational principle for Eulerian fluid flows Holm et al. (1998).

The TRSW, TL1, and TQG models separate the wave and current aspects of their flows into gravity waves and Rossby waves on the free surface, and fluid circulation in the region between the free surface and the bottom topography. The Kelvin circulation theorems for TRSW, TL1, and TQG show that horizontal gradients of the buoyancy in the fluid region (e.g., at thermal fronts) can couple to either the elevation gradient at the upper interface, or to the bathymetry gradient at the lower interface. Namely, horizontal circulation of the fluid is produced whenever either of the gradients at the upper and lower interfaces are misaligned with the horizontal gradient of the buoyancy. Thus, both waves on the surface and variations of the bottom topography can create horizontal fluid circulation when the buoyancy is spatially inhomogeneous.

The generation of submesoscale circulations involves a wide range of time scales, as well as many couplings among the various degrees of freedom and the boundaries. The “irreducible imprecision” of numerical simulations McWilliams (2007) and the sparsity of observed data in both space and time produce uncertainty in forecasts of rapid, high-wavenumber GFD processes and thereby present a “grand challenge” for data assimilation.

In preparation for meeting this challenge, the present paper develops the stochastic variational principles for the TRSW, TL1, and TQG models. This mathematical foundation is needed in applying the SALT (Stochastic Advection by Lie Transport) approach to the derivation of stochastic fluid equations which preserve the geometric structure of fluid dynamics Holm (2015). The physical effects of stochasticity in the SALT approach for deriving the stochastic TRSW, TL1, and TQG models are revealed in their Kelvin circulation theorems. Namely, the corresponding material loops defining the circulation integrals for these models are shown to move along stochastic Lagrangian paths. The motivations and recent results of applications of the SALT approach for uncertainty quantification and for data assimilation are laid out in Cotter et al. (2018; 2019).

Remark 1.1 (Separation of scales behavior in TQG vs QG). High Reynolds-number two-dimensional Navier–Stokes turbulence relaxes through a combination of vortex merger and filamentation, in which energy tends to flow to large scales and enstrophy tends to flow and be dissipated on fine scales. The combined conservation of two integral conservation laws of energy and enstrophy produces the famous inverse cascade Kraichnan (1967) , which also occurs in QG turbulence. However, TQG turbulence is different. First, TQG does not preserve an enstrophy because there are two active degrees of freedom in TQG, both the momentum and the buoyancy. Second, TQG possesses high-wavenumber instabilities. That is, linear analysis reveals that high wavenumbers in a TQG flow can suddenly become unstable. Turbulence tends to mimic the properties and locations of its energy source. Consequently, numerical simulations of the nonlinear processes in TQG reveal the nonlinear sudden creation of coherent structures at the scale of highest growth rate at the onset of the high-wavenumber, small-lengthscale instabilities. Thus, unlike QG turbulence, a scale-separation exists in TQG turbulence. Therefore, one can expect TQG dynamics to mimic QG dynamics at low wavenumbers and to transform into another type of motion when the flow enters the onset of the high-wavenumber instability for TQG. When this happens, the turbulence will be limited to the scale of the wavelength at which the highest growth-rate occurs, which might be significantly smaller than the scale of the entire domain, see Fig. 3. This means TQG possesses a scale separation in its solution behavior which can be radically different from the solution behavior of other models in the QG family. In particular, one can expect that TQG will require different approaches from QG in the methods of its analysis and parameterization. Indeed, the misalignment of gradients of bathymetry with gradients of buoyancy can create circulation even when there is no initial flow. Because of the scale gap in its solution behavior, one may expect that TQG may require different methods of analysis, parameterization and numerical implementation in the same computational simulation of a given TQG flow.

This paper is organized as follows: In Sec. II, we review the deterministic TRSW model by re-deriving its equations in the Euler–Poincaré variational framework of Holm et al. (1998). In the Euler–Poincaré framework, we prove the Kelvin–Noether circulation theorem and discuss steady solution properties of the deterministic TRSW equations. This derivation of the TRSW equations with stochastic advection by Lie transport (SALT) is intended to be the mathematical foundation for a systematic means of introducing data-driven parametrizations of stochastic transport for uncertainty quantification and data assimilation for upper ocean dynamics. This type of uncertainty quantification and data assimilation has already been accomplished using this approach for the 2D Euler equations in a square domain with fixed boundaries and the 2-layer QG equations in a periodic channel, in Cotter et al. (2018); Cotter et al. (2019), respectively. In Sec. III, we discuss the deterministic thermal Eliassen approximation, or TL1 model, of TRSW, as derived from a combination of the Euler–Poincaré variational approach and asymptotic expansions in the vorticity–divergence representation of the fluid velocity. In the derivation of TL1, we use a modified version of the Euler–Poincaré framework introduced in Allen and Holm (1996) which expresses the approximate momentum in terms of gradients of advected quantities. In Sec. IV, we derive the TL1 equations with stochastic advection by Lie transport (SALT) in the modified Euler–Poincaré variational framework. These stochastic TL1 equations would be useful for uncertainty quantification and data assimilation at this intermediate level of approximation. In Sec. V, we take the next step in the asymptotic expansion to derive and discuss the deterministic version of the thermal quasi-geostrophic (TQG) equations in Sec. V A. Section V B specifies the numerical details of an example implementation of the TQG solution shown in Fig. 4. We then investigate the Hamiltonian framework of the TQG equations. The Hamiltonian formulation of the TQG equations can be used to derive the SALT stochastic TQG equations. This is described in Sec. V C. In Sec. VI, we conclude by outlining a few next steps and open problems to which the present work has led us, but which we feel are beyond the scope of the present paper.

II. THE THERMAL ROTATING SHALLOW WATER (TRSW) MODEL

The thermal rotating shallow water (TRSW) model describes the motion of a single two dimensional layer of fluid with horizontally varying buoyancy and bottom topography (or, bathymetry). The TRSW model is an extension of the rotating shallow water (RSW) model and a simplification of the various three dimensional models such as the primitive equations and the Euler–Boussinesq model, which are commonly used for computationally simulating large-scale ocean and atmosphere circulation dynamics. The thermal rotating shallow water equations may also be interpreted as a model for an upper active layer of fluid on top of a lower inert layer. For that reason, the TRSW model is sometimes called a 1.5 layer model Warneford and Dellar (2013). A stochastic version of this model has already been derived from a variational point of view in Appendix B of Holm and Luesink (2019). For a related deterministic discussion of a fully multilayer variational model with nonhydrostatic pressure, see Cotter et al. (2010).

A. Deterministic TRSW equations

The deterministic TRSW equations in a rotating planar domain $D \in ℝ^{2}$ with boundary $\partial D$ are expressed using the following notation. The depth is denoted $η = η (x, t)$ ⁠, where $x = (x, y)$ is the horizontal vector position, and t is time. The (non-negative) horizontal buoyancy is written as $b (x, t) = ρ (x, t) / \bar{ρ}$ ⁠, where $ρ (x, t)$ is the mass density, $\bar{ρ}$ is the uniform reference mass density. The nondimensional deterministic TRSW equations for the Eulerian horizontal vector velocity $u (x, t)$ ⁠, thickness $η (x, t)$ ⁠, and buoyancy $b (x, t)$ of the active fluid layer are given by

\begin{matrix} \frac{D}{D t} u + \frac{1}{Ro} f \hat{z} \times u = - \frac{α}{{Fr}^{2}} \nabla ((1 + s b) ζ) + \frac{s}{2 {Fr}^{2}} (α ζ - h) \nabla b, \\ \frac{\partial η}{\partial t} + \nabla \cdot (η u) = 0, \frac{D}{D t} b = 0 . \end{matrix}

(2.1)

The other notation is f for the Coriolis parameter, $α ζ = η - h$ for the free surface elevation, where $h (x)$ is the time-independent mean depth, $s$ for the stratification parameter, α for the wave amplitude, $Ro = U / (f_{0} L)$ for the Rossby number, and ${Fr}^{2} = U^{2} / (g H)$ for the Froude number. The material time derivative for scalar advected quantities is denoted by $\frac{D}{D t} : = \partial_{t} + u \cdot \nabla$ ⁠. In defining the dimensionless numbers, U denotes the horizontal velocity scale, L is the horizontal length scale, f₀ is the typical rotation frequency, g is the gravitational acceleration, and H is the typical depth. The stratification parameter $s$ is introduced so that the buoyancy variable has size $O (1)$ and the importance of buoyancy is governed by the size of the stratification parameter. The stratification parameter is particularly important in introducing the Boussinesq approximation in Holm and Luesink [2019]. The Boussinesq approximation is necessary to derive the thermal rotating shallow water equations. The wave amplitude α is the typical free surface elevation ζ₀ divided by the typical depth H and is introduced so that ζ has size $O (1)$ and the importance of free surface waves is governed by the size of the wave amplitude. The notation $\hat{z}$ is used to denote the unit vector perpendicular to the flow domain $D$ ⁠. The boundary conditions are

\hat{n} \cdot u = 0 and \hat{n} \times \nabla b = 0 on the boundary \partial D,

(2.2)

meaning that fluid velocity u is tangential and buoyancy b is constant on the boundary $\partial D$ ⁠, fixed in the frame rotating at time-independent angular frequency $\hat{z} f (x) / 2$ ⁠. The boundary conditions in (2.2) are required to guarantee that circulations along the boundary are constant. This is explained in Ripa (1993), who also explains the Hamiltonian formulation of the thermal rotating shallow water equations and results on the stability of equilibrium flows. In the boundary conditions, $\hat{n}$ denotes the outward unit normal. Periodic boundary conditions may also be considered.

1. Variational formulation

The TRSW equations (2.1) can be derived by means of the Euler–Poincaré variational principle, as is shown in Bröcker et al. (2018). When the equations of motion are derived in this framework, there is a natural way to express three fundamental relations. The first fundamental relation is the Kelvin circulation theorem, the second one is the advection equation for potential vorticity, and the third one is an infinity of conserved integral quantities arising from Noether's theorem for the symmetry of Eulerian fluid quantities under Lagrangian particle relabeling. For example, in rotating shallow water, without a buoyant scalar, the enstrophy is among these integral quantities. This framework turns out to be ideal for introducing stochasticity, as shown in Holm (2015) and de Leon et al. (2020), and one can similarly introduce rough paths Crisan et al. (2020) into continuum mechanics. In applications, this framework provides a means to consistently introduce data-driven parametrizations of stochastic transport into a large class of fluid models Cotter et al. (2018; 2019).

B. The Euler–Poincaré theorem

1. Variational derivatives of functionals

The Euler–Poincaré theorem relies on variational derivatives of functionals. This type of derivative is given by the following definition:

Definition 2.1. A functional $F [ρ]$ is defined as a map $F : ρ \in B \to ℝ$ , where $B$ is a Banach space. The variational derivative of $F (ρ)$ , denoted $δ F / δ ρ$ , is defined by the linear functional

δ F [ρ] : = \lim_{ε \to 0} \frac{F [ρ + ε ϕ] - F [ρ]}{ε} = : \frac{d}{d ε} F [ρ + ε ϕ] |_{ϵ = 0} = \int_{Ω} \frac{δ F}{δ ρ} (x) ϕ (x) d x = : ⟨ \frac{δ F}{δ ρ}, ϕ ⟩ .

(2.3)

In this definition, $ε \in ℝ$ is a real parameter, $ϕ$ is an arbitrary smooth function, and the angle brackets $⟨ \cdot, \cdot ⟩$ indicate L² pairing. The function $ϕ (x)$ above is called the “variation of ρ” and will be denoted as $δ ρ : = ϕ (x)$ . Since the variation is a linear operator on functionals, we can define the functional derivative δ in (2.3) operationally as

δ F [ρ] : = ⟨ \frac{δ F}{δ ρ}, δ ρ ⟩ .

(2.4)

2. Euler–Poincaré theorem

Given the boundary conditions and definitions above, the following form of the Euler–Poincaré theorem will provide the deterministic equations of motion derived from Hamilton's principle. Suppose a deterministic Lagrangian functional $ℓ : X \times V^{*} \to ℝ$ is defined on the domain of flow, $D$ ⁠. Here, $X$ denotes the space of smooth vector fields on $D$ and $V^{*}$ is the vector space of advected quantities. Advected quantities are tensor fields of various types which are preserved along the flow. The space of smooth vector fields is a Lie algebra under the action of the Jacobi–Lie bracket, which is denoted as $[\cdot, \cdot] : X \times X \to X$ ⁠, and is defined for two vector fields $u, v$ by the commutator relation

[u, v] : = ((u \cdot \nabla) v - (v \cdot \nabla) u) .

(2.5)

Theorem 2.1 [Euler–Poincaré equations Holm et al. (1998)].

The following two statements are equivalent:

Hamilton's variational principle in Eulerian coordinates, with $u \in X (D)$ and $b, η \in V^{*} (D)$ ⁠,
$δ S : = δ \int_{t_{1}}^{t_{2}} ℓ (u, b, η) d t = 0,$
(2.6)
holds on $X (D) \times V^{*}$ , using variations of the form
$\begin{matrix} δ u = \frac{\partial}{\partial t} v - [u, v], δ b = - (v \cdot \nabla) b, \\ δ η = - \nabla \cdot (η v), \end{matrix}$
(2.7)
where the vector field $v \in X (D)$ is arbitrary and vanishes on the endpoints t₁ and t₂.
The Euler–Poincaré equations hold. These equations are
$\frac{\partial}{\partial t} \frac{δ ℓ}{δ u} + (u \cdot \nabla) \frac{δ ℓ}{δ u} + (\nabla u) \cdot \frac{δ ℓ}{δ u} + \frac{δ ℓ}{δ u} (\nabla \cdot u) = - \frac{δ ℓ}{δ b} \nabla b + η \nabla \frac{δ ℓ}{δ η},$
(2.8)
or, equivalently, in two-dimensional vector calculus notation,
$\frac{\partial}{\partial t} \frac{δ ℓ}{δ u} + (\nabla^{⊥} \cdot \frac{δ ℓ}{δ u}) u^{⊥} + \nabla (u \cdot \frac{δ ℓ}{δ u}) + \frac{δ ℓ}{δ u} (\nabla \cdot u) = - \frac{δ ℓ}{δ b} \nabla b + η \nabla \frac{δ ℓ}{δ η},$
(2.9)
or, finally, as an embedding in three dimensional space,
$\frac{\partial}{\partial t} \frac{δ ℓ}{δ u} - u \times (\nabla \times \frac{δ ℓ}{δ u}) + \nabla (u \cdot \frac{δ ℓ}{δ u}) + \frac{δ ℓ}{δ u} (\nabla \cdot u) = - \frac{δ ℓ}{δ b} \nabla b + η \nabla \frac{δ ℓ}{δ η},$
(2.10)
with advection equations
$\frac{\partial}{\partial t} b = - u \cdot \nabla b and \frac{\partial}{\partial t} η = - \nabla \cdot (η u) .$
(2.11)

Remark 2.2. The abstract statement of the Euler–Poincaré Theorem 2.1, formulated on general semidirect product Lie groups and its proof, is presented in Holm et al. (1998) deterministically, in Holm (2015) and de Leon et al. (2020 ) stochastically and finally in Crisan et al. (2020) on rough paths. We shall state theorem 2.1 again in stochastic form, where we shall also provide a proof.

Remark 2.3. In Theorem 2.1, the operator δ in (2.6) is the functional derivative defined in (2.3), the brackets $[\cdot, \cdot]$ denote the commutator of vector fields defined in (2.5), and $v \in X (D)$ is an arbitrary vector field in two dimensions which vanishes at the endpoints in time, t₁ and t₂. Equation (2.8) is the advective form of the Euler–Poincaré equation, where the operator $(\nabla u) \cdot a = {(\nabla u)}^{T} a = a_{j} \nabla u^{j}$ . Equations (2.9) and (2.10) are the curl form of the Euler–Poincaré equation. They are equivalent and can be transformed into each other by the conventions $u^{⊥} = \hat{z} \times u$ and $\nabla^{⊥} = \hat{z} \times \nabla$ , where $\hat{z}$ is the outward unit vector perpendicular to the planar domain $D$ . Equations (2.9) and (2.10) are used to derive the vorticity equation and are useful in the energy-Casimir stability method, see Holm et al. (1985).

Remark 2.4. One may interpret the Euler–Poincaré equation (2.8) in terms of Newton's law for fluid motion. Namely, the rate of change of the covector momentum density $P : = δ ℓ / δ u$ along the flow of the fluid velocity vector field u equals the sum of covector force densities on the right-hand side of equation (2.8).

C. Kelvin–Noether circulation theorem

A straightforward calculation using the second advection equation in (2.11) shows that (2.8) may be written equivalently as follows.

Lemma 2.5. The Euler–Poincaré equation in (2.8) is equivalent to the following:

\frac{\partial}{\partial t} (\frac{1}{η} \frac{δ ℓ}{δ u}) + (u \cdot \nabla) (\frac{1}{η} \frac{δ ℓ}{δ u}) + (\nabla u) \cdot (\frac{1}{η} \frac{δ ℓ}{δ u}) = - \frac{1}{η} \frac{δ ℓ}{δ b} \nabla b + \nabla \frac{δ ℓ}{δ η} .

(2.12)

One of the main features of Theorem 2.1 for fluid dynamics is that its Euler–Poincaré equations satisfy the following Kelvin circulation theorem.

Theorem 2.6 (Kelvin–Noether circulation). For an arbitrary loop c(t) which is advected by the velocity field u, the following dynamics holds for the circulation integral $I$ , given by

I : = \oint_{c (u)} \frac{1}{η} \frac{δ ℓ}{δ u} \cdot d x, \frac{d I}{d t} = - \oint_{c (u)} (\frac{1}{η} \frac{δ ℓ}{δ b}) \nabla b \cdot d x .

(2.13)

Remark 2.7. The notation $c (u)$ indicates that the material loop c is transported by the flow $ϕ_{t}$ which is generated by the vector field $u : = u \cdot \nabla$ . To be precise, $c (u) = ϕ_{t *} c (0)$ , where $ϕ_{t *}$ is the pull-back by the inverse of the flow $ϕ_{t}$ , also known as the push-forward by $ϕ_{t}$ ⁠.

Proof. The Kelvin circulation law (2.13) follows from Newton's law of motion obtained from the Euler–Poincaré equation (2.12) for the evolution of momentum per unit mass $η^{- 1} δ ℓ / δ u$ concentrated on an advecting material loop, $c (t) = ϕ_{t} c (0)$ ⁠, where $ϕ_{t}$ is the flow map which is generated by the vector field u. Upon changing variables by pulling back the integrand to its initial position, the time derivative can be moved inside the integral and the product rule may be applied. Then, by inverting the pull-back we obtain the following:

\begin{matrix} \frac{d}{d t} \oint_{c (u)} \frac{1}{η} \frac{δ ℓ}{δ u} \cdot d x & = \oint_{c (u)} (\frac{\partial}{\partial t} + u \cdot \nabla + (\nabla u) \cdot) (\frac{1}{η} \frac{δ ℓ}{δ u}) \cdot d x, \\ = - \oint_{c (u)} \frac{1}{η} \frac{δ ℓ}{δ b} \nabla b \cdot d x + \oint_{c (u)} \nabla \frac{δ ℓ}{δ η} \cdot d x, \\ = - \oint_{c (u)} (\frac{1}{η} \frac{δ ℓ}{δ b}) \nabla b \cdot d x . \end{matrix}

(2.14)

In the second line, we have used the Euler–Poincaré equation (2.8) and the advection equation for the density. The last step applies the fundamental theorem of calculus to prove vanishing of the last loop integral in the second line. For the corresponding proof in the general case, see Holm et al. (1998).

Corollary 2.7.1. By Stokes Law, Eq. (2.14) in the Kelvin–Noether circulation theorem 2.6 implies

\frac{d I}{d t} = - \int \int_{\partial S = c (u)} \hat{z} \cdot \nabla (\frac{1}{η} \frac{δ ℓ}{δ b}) \times \nabla b d x d y .

(2.15)

Therefore, circulation is created by misalignment of the gradients of buoyancy b and its dual quantity $η^{- 1} δ ℓ / δ b$ ⁠.

The thermal rotating shallow water equations (2.1) in a planar domain $D$ are obtained by applying the Euler–Poincaré theorem 2.1 to the Lagrangian

\begin{matrix} ℓ_{TRSW} (u, η, b) = \int_{D} \frac{1}{2} η | u |^{2} + \frac{1}{Ro} η u \cdot R (x) - \frac{1}{2 {Fr}^{2}} (1 + s b) (η^{2} - 2 η h) d x d y . \end{matrix}

(2.16)

This is easily shown once we have computed the variational derivatives of the Lagrangian, since these derivatives can simply be substituted into (2.8), or into one of the equivalent formulations (2.9) or (2.10). These variational derivatives are obtained by using definition 2.1 to find

\begin{matrix} \frac{δ ℓ_{TRSW}}{δ u} = η u + η R, \\ \frac{δ ℓ_{TRSW}}{δ η} = \frac{1}{2} | u |^{2} + \frac{1}{Ro} u \cdot R - \frac{α}{{Fr}^{2}} (1 + s b) ζ, \\ \frac{δ ℓ_{TRSW}}{δ b} = - \frac{s}{2 {Fr}^{2}} (η^{2} - 2 η h) . \end{matrix}

(2.17)

With these variations, we obtain (2.1) and by means of theorem 2.6, we see that the TRSW equations satisfy the following Kelvin circulation theorem.

Theorem 2.8 (Kelvin theorem for deterministic TRSW). The deterministic TRSW equations (2.1) imply the following Kelvin circulation law:

\frac{d}{d t} \oint_{c (u)} (u + \frac{1}{Ro} R) \cdot d x = \frac{s}{2 {Fr}^{2}} \oint_{c (u)} (α ζ - h) \nabla b \cdot d x = \frac{s}{2 {Fr}^{2}} \int \int_{\partial S = c (u)} \hat{z} \cdot \nabla (α ζ - h) \times \nabla b d x d y,

(2.18)

where c(u) is any CLOSED loop moving with horizontal fluid velocity $u (x, t)$ in two dimensions.

Proof. This result follows from the Kelvin–Noether theorem 2.6 for Euler–Poincaré fluid equations.

Remark 2.9. Equation (2.18) implies that misalignment between the horizontal gradients of either the free surface elevation $α ζ$ , or the bathymetry h with the buoyancy b will generate circulation in a horizontal plane. When the free surface elevation is negligible, that is when the wave amplitude $α ≪ 1$ , circulation is still being generated due to gradients of bathymetry misaligning with gradients of buoyancy. When the contribution of buoyancy is negligible, that is when the stratification parameter $s ≪ 1$ , circulation is conserved. In the sections to come, we will use theorem 2.18 to interpret the properties of the approximate equations we will derive.

Corollary 2.9.1 (Circulation on the boundary). The circulation $\oint_{Γ_{k}} (u + \frac{1}{Ro} R) \cdot d x$ on each connected component of the boundary $Γ_{k} \in \partial D$ is conserved by the deterministic TRSW equations.

Proof. Preservation of circulation on each connected component of the boundary follows from the boundary conditions in (2.2) and Kelvin's theorem for TRSW in (2.18). The first boundary condition in (2.2) implies that the velocity is tangent to the boundary. Hence, a circuit $c (u)$ on the boundary remains on the boundary. Consequently, Kelvin's theorem for TRSW in (2.18) applies to a boundary circuit. The second boundary condition in (2.2) implies that $\nabla b \cdot d x = 0$ on the boundary $\partial D$ ⁠. Hence, the right-hand side of (2.18) vanishes for a circuit $c (u)$ on the boundary and the circulation $\oint_{Γ_{k}} (u + \frac{1}{Ro} R) \cdot d x$ is conserved on the boundary.

The potential vorticity for the thermal rotating shallow water equations (2.1) is defined as

q : = \frac{1}{η} \hat{z} \cdot \nabla \times (u + \frac{1}{Ro} R) .

(2.19)

Even though this is the same definition of potential vorticity as for the rotating shallow water equations, in thermal rotating shallow water, the potential vorticity is not conserved along Lagrangian paths. Rather, the potential vorticity satisfies the following equation:

\frac{\partial}{\partial t} q + u \cdot \nabla q = \frac{s}{2 {Fr}^{2} η} \hat{z} \cdot \nabla (α ζ - h) \times \nabla b .

(2.20)

Not unexpectedly, the mechanism responsible for the generation of circulation in the Kelvin circulation theorem 2.8 is also rate of creation of potential vorticity, q, along fluid particle trajectories in Eq. (2.20). When the horizontal buoyancy gradients are negligible, that is $s ≪ 1$ ⁠, potential vorticity is conserved along Lagrangian paths.

1. Conservation laws for deterministic TRSW

The deterministic TRSW equations (2.1) conserve the energy

E_{TRSW} (u, η, b) = \frac{1}{2} \int_{D} η | u |^{2} + \frac{1}{{Fr}^{2}} (1 + s b) (η^{2} - 2 η h) d x d y .

(2.21)

The conservation of energy (2.21) can be proved directly by using the TRSW equations (2.1) and the boundary conditions (2.2). The TRSW equations (2.1) also conserve an infinity of integral conservation laws, determined by two arbitrary differentiable functions of buoyancy $Φ (b)$ and $Ψ (b)$ as

C_{Φ, Ψ} = \int_{D} η Φ (b) + ϖ Ψ (b) d x d y = \int_{D} (Φ (b) + q Ψ (b)) η d x d y,

(2.22)

where $ϖ = η q = \hat{z} \cdot \nabla \times (u + \frac{1}{Ro} R)$ is the total vorticity. That is, for any choice of differentiable $Φ$ and $Ψ$ ⁠, the quantity $C_{Φ, Ψ}$ is conserved in time. The conservation of (2.37) can also be proved as a direct calculation using Eqs. (2.1) and the boundary conditions (2.2).

2. Noether's theorem

Conservation of the integral quantities in Eqs. (2.21) and (2.22) is associated by Noether's theorem with smooth transformations which leave invariant the Eulerian fluid quantities in the Lagrangian Abraham and Marsden (1978). For example, conservation of energy (2.21) arises from invariance of the Lagrangian in (2.16) under translations in time; since this Lagrangian does not depend explicitly on time. Likewise, the conserved quantities in (2.22) are associated by Noether's theorem with the smooth flows which translate the fluid parcels along steady solutions of the equations of motion, because these transformations preserve the Eulerian fluid variables in the Lagrangian Holm et al. (1985). Upon introducing stochasticity via the Euler–Poincaré theorem in Sec. II D, the latter transformations and their Noether conservation laws will persist. However, energy conservation will not persist because the stochastic Lagrangian will depend explicitly on time through the Brownian noise. The geometrical significance of the conservation laws in Eq. (2.22) which persist for stochastic TRSW will be discussed further in remark 2.15.

D. TRSW with stochastic advection by Lie transport (SALT)

By modifying the fluid transport vector field in the Euler–Poincaré theorem 2.1, one can derive the stochastic equations of motion which preserve the geometric properties of their deterministic counterparts. Following Holm (2015), we introduce the stochastic vector field for fluid transport in semimartingale form

d χ_{t} : = u (x, t) d t + \sum_{k = 1}^{M} ξ_{k} (x) ○ d W_{t}^{k},

(2.23)

where $x \in D \subset ℝ^{2}$ ⁠. The time-independent vector fields $ξ_{k}$ ⁠, with $k = 1, 2, \dots, M$ ⁠, represent spatially correlated sources of temporal uncertainty. The vector fields $ξ_{k}$ must satisfy the same boundary condition as u. The circle notation $°$ means that the stochastic integral is to be understood in the Stratonovich sense. Stratonovich calculus possesses the ordinary chain rule and product rule, which are crucial for defining the variational derivative. These properties are written in integral form, though, because stochastic equations are not differentiable with respect to time. The sources of the stochasticity are the independent, identically distributed Brownian motions $W_{t}^{k}$ associated with each $ξ_{k}$ ⁠. The Brownian motions are defined with respect to the standard probability space, see Ito (1984). One may regard the $ξ_{k} (x)$ as eigenvectors of the velocity–velocity correlation tensor. In practice, the eigenvectors $ξ_{k} (x)$ are expected to be obtained using the SALT algorithm developed in Cotter et al. (2018; 2019). This algorithm is based on empirical orthogonal function analysis, and the number M of $ξ_{k}$ needed in (2.23) would be decided by how much of the variance is required to be represented.

The SALT version of Theorem 2.1 may be stated, as follows:

Theorem 2.10 [Stochastic Euler–Poincaré equations Holm (2015); de Leon et al. (2020)].

The following two statements are equivalent:

The stochastic Hamilton's variational principle in Eulerian coordinates, with $u \in X (D)$ and $b, η \in V^{*} (D)$ ⁠,
$δ S : = δ \int_{t_{1}}^{t_{2}} ℓ (u, b, η) d t = 0,$
(2.24)
holds on $X (D) \times V^{*}$ , using variations of the form
$\begin{matrix} δ u d t = d v - [d χ_{t}, v], δ b d t = - (v \cdot \nabla) b d t, \\ δ η d t = - \nabla \cdot (η v) d t, \end{matrix}$
(2.25)
where the vector field $v \in X (D)$ is arbitrary and vanishes on the endpoints t₁ and t₂ and the semimartingale vector field $d χ_{t}$ is defined in (2.23).
The stochastic Euler–Poincaré equations hold. These equations are
$d \frac{δ ℓ}{δ u} + (d χ_{t} \cdot \nabla) \frac{δ ℓ}{δ u} + (\nabla d χ_{t}) \cdot \frac{δ ℓ}{δ u} + \frac{δ ℓ}{δ u} (\nabla \cdot d χ_{t}) = - \frac{δ ℓ}{δ b} \nabla b d t + η \nabla \frac{δ ℓ}{δ η} d t,$
(2.26)
or, equivalently, either in two dimensional vector calculus notation,
$d \frac{δ ℓ}{δ u} + (\nabla^{⊥} \cdot \frac{δ ℓ}{δ u}) d χ_{t}^{⊥} + \nabla (d χ_{t} \cdot \frac{δ ℓ}{δ u}) + \frac{δ ℓ}{δ u} (\nabla \cdot d χ_{t}) = - \frac{δ ℓ}{δ b} \nabla b d t + η \nabla \frac{δ ℓ}{δ η} d t,$
(2.27)
or as an embedding in three dimensional space,
$d \frac{δ ℓ}{δ u} - d χ_{t} \times (\nabla \times \frac{δ ℓ}{δ u}) + \nabla (d χ_{t} \cdot \frac{δ ℓ}{δ u}) + \frac{δ ℓ}{δ u} (\nabla \cdot d χ_{t}) = - \frac{δ ℓ}{δ b} \nabla b d t + η \nabla \frac{δ ℓ}{δ η} d t,$
(2.28)
with advection equations
$d b = - d χ_{t} \cdot \nabla b and d η = - \nabla \cdot (η d χ_{t}) .$
(2.29)

For the proof of the general version of the theorem and the technical details, we refer to Holm (2015) and de Leon et al. (2020). The proof for the theorem above is given below.

Proof. Hamilton's variational principle implies

\begin{matrix} 0 = \int_{t_{1}}^{t_{2}} [{⟨ \frac{δ ℓ}{δ u}, δ u ⟩}_{X} + {⟨ \frac{δ ℓ}{δ b}, δ b ⟩}_{V^{*}} + {⟨ \frac{δ ℓ}{δ η}, δ η ⟩}_{V^{*}}] d t, \\ = \int_{t_{1}}^{t_{2}} [{⟨ \frac{δ ℓ}{δ u}, d v - [d χ_{t}, v] ⟩}_{X} + {⟨ \frac{δ ℓ}{δ b}, - (v \cdot \nabla) b ⟩}_{V^{*}} + {⟨ \frac{δ ℓ}{δ η}, - \nabla \cdot (η v) ⟩}_{V^{*}}] d t, \\ = \int_{t_{1}}^{t_{2}} [{⟨ - \frac{\partial}{\partial t} \frac{δ ℓ}{δ u} - (d χ_{t} \cdot \nabla) \frac{δ ℓ}{δ u} - (\nabla d χ_{t}) \cdot \frac{δ ℓ}{δ u} + \frac{δ ℓ}{δ u} (\nabla \cdot d χ_{t}), v ⟩}_{X} + {⟨ - \frac{δ ℓ}{δ b} \nabla b, v ⟩}_{X} + {⟨ η \nabla \frac{δ ℓ}{δ η}, v ⟩}_{X}] d t . \end{matrix}

The subscripts $X$ and $V^{*}$ on the L² pairings indicate over which space that the pairing is defined. Since v is arbitrary and vanishes at the endpoints t₁ and t₂ in time, the following equation holds:

\frac{\partial}{\partial t} \frac{δ ℓ}{δ u} + (d χ_{t} \cdot \nabla) \frac{δ ℓ}{δ u} + (\nabla d χ_{t}) \cdot \frac{δ ℓ}{δ u} + \frac{δ ℓ}{δ u} (\nabla \cdot d χ_{t}) = - \frac{δ ℓ}{δ b} \nabla b + η \nabla \frac{δ ℓ}{δ η} .

This finishes the proof of the stochastic Euler–Poincaré equation in (2.26). The equivalent forms in Eqs. (2.27) and (2.28) follow by means of a standard vector identity.

By taking the variational derivatives of the Lagrangian for thermal rotating shallow water as in (2.17), we obtain the stochastic TRSW equations

\begin{matrix} d u + (d χ_{t} \cdot \nabla) u + \sum_{k} (\nabla ξ_{k}) \cdot u ○ d W_{t}^{k} = - \frac{α}{{Fr}^{2}} \nabla ((1 + s b) ζ) d t + \frac{s}{2 {Fr}^{2}} (α ζ - h) \nabla b d t \\ - \frac{1}{Ro} f \hat{z} \times d χ_{t} - \frac{1}{Ro} \sum_{k} \nabla (ξ_{k} \cdot R) ○ d W_{t}^{k}, \end{matrix}

\begin{matrix} - \frac{1}{Ro} f \hat{z} \times d χ_{t} - \frac{1}{Ro} \sum_{k} \nabla (ξ_{k} \cdot R) ○ d W_{t}^{k}, \\ d η + \nabla \cdot (η d χ_{t}) = 0, \\ d b + (d χ_{t} \cdot \nabla) b = 0 . \end{matrix}

(2.30)

The boundary conditions are given by

\hat{n} \cdot u = 0 and \hat{n} \cdot ξ_{k} = 0 and \hat{n} \times \nabla b = 0 on the boundary \partial D .

(2.31)

The boundary condition on $ξ_{k}$ is required to be satisfied for each k. The Kelvin circulation theorem has now become stochastic, because the circulation loop is transported by the stochastic vector field $d χ_{t}$ ⁠, rather than by the deterministic vector field u. Specifically, we have

Theorem 2.11. The stochastic Kelvin circulation law associated with the stochastic Euler–Poincaré theorem is

d \oint_{c (d χ_{t})} \frac{1}{η} \frac{δ ℓ}{δ u} \cdot d x = - \oint_{c (d χ_{t})} \frac{1}{η} \frac{δ ℓ}{δ b} \nabla b \cdot d x d t,

(2.32)

where $c (d χ_{t})$ is a closed loop that is transported by the flow generated by the stochastic fluid velocity $d χ_{t}$ in two dimensions.

Proof. By following the proof (2.14) of the deterministic Kelvin circulation theorem 2.6 and using the product rule and chain rule for the stochastic time differential d, we have

\begin{matrix} d \oint_{c (d χ_{t})} \frac{1}{η} \frac{δ ℓ}{δ u} \cdot d x = \oint_{c (d χ_{t})} (d + d χ_{t} \cdot \nabla + (\nabla d χ_{t}) \cdot) (\frac{1}{η} \frac{δ ℓ}{δ u}) \cdot d x d t, \\ = - \oint_{c (d χ_{t})} \frac{1}{η} \frac{δ ℓ}{δ b} \nabla b \cdot d x d t + \oint_{c (d χ_{t})} \nabla \frac{δ ℓ}{δ η} \cdot d x d t, \\ = - \oint_{c (d χ_{t})} (\frac{1}{η} \frac{δ ℓ}{δ b}) \nabla b \cdot d x d t . \end{matrix}

(2.33)

Remark 2.12. For the stochastic TRSW equations (2.30), we have

d \oint_{c (d χ_{t})} (u + \frac{1}{Ro} R) \cdot d x = \frac{s}{2 {Fr}^{2}} \oint_{c (d χ_{t})} (α ζ - h) \nabla b \cdot d x d t = \frac{s}{2 {Fr}^{2}} \int \int_{\partial S = c (d χ_{t})} \hat{z} \cdot \nabla (α ζ - h) \times \nabla b d x d y d t .

(2.34)

One sees in Eq. (2.34) that misalignment of the horizontal gradient of either the free surface elevation ζ, or the bathymetry h with the horizontal gradient of the buoyancy b will generate circulation, cf. the corresponding deterministic TRSW Kelvin circulation theorem in Eq. (2.18).

Remark 2.13. The evolution of potential vorticity on fluid parcels for the TRSW equations in (2.30) is given by

d q + (d χ_{t} \cdot \nabla) q = \frac{s}{2 {Fr}^{2} η} \hat{z} \cdot \nabla (α ζ - h) \times \nabla b d t,

(2.35)

where the potential vorticity is defined by

q : = \frac{ϖ}{η}, and ϖ : = \hat{z} \cdot \nabla \times (u + \frac{1}{Ro} R) .

(2.36)

Remark 2.14. The stochastic TRSW equations (2.30) have an infinite number of conserved integral quantities

C_{Φ, Ψ} = \int_{D} (Φ (b) + q Ψ (b)) η dxdy,

(2.37)

for the boundary conditions given in (2.2) and any differentiable functions $Φ$ and $Ψ$ . This can be proved directly by computing the stochastic differential of (2.37), integrating by parts and using the boundary conditions (2.31).

Remark 2.15 (Stochastic Hamiltonian formulation). The Legendre transform which determines the Hamiltonian $d h_{TRSW}$ for the stochastic TRSW equations is defined as

d ℏ_{TRSW} (μ, η, b) : = ⟨ μ, d χ_{t} ⟩ - ℓ_{TRSW} (u, η, b) d t,

(2.38)

in which the angle brackets in the definition of the Legendre transform denote the L² pairing over the domain $D$ . Notice that the Hamiltonian $d h_{TRSW}$ in (2.38) is a semimartingale. SeeStreet and Crisan (2020 ) for variational principles driven by semimartingales. The Hamiltonian form of the stochastic TRSW equations is given for a functional $F (μ_{i}, η, b)$ by

d F = {F, d ℏ_{TRSW}} = - \int_{D} {[\begin{matrix} δ F / δ μ_{i} \\ δ F / δ η \\ δ F / δ b \end{matrix}]}^{T} [\begin{matrix} μ_{j} \partial_{i} + \partial_{j} μ_{i} & η \partial_{i} & - b_{, i} \\ \partial_{j} η & 0 & 0 \\ b_{, j} & 0 & 0 \end{matrix}] [\begin{matrix} δ (d ℏ_{TRSW}) / δ μ_{j} \\ δ (d ℏ_{TRSW}) / δ η \\ δ (d ℏ_{TRSW}) / δ b \end{matrix}] d x d y .

(2.39)

In this notation, repeated indices are summed over. The conserved integral quantities $C_{Φ, Ψ}$ defined in (2.37) are Casimirs of the Lie–Poisson bracket in (2.39). That is, the vector of variational derivatives of $C_{Φ, Ψ}$ comprises a null eigenvector of the Lie–Poisson bracket in (2.39). Consequently, their conservation persists when the Hamiltonian is made stochastic. This means that the solutions of these equations describe stochastic coadjoint motion in function space on level sets of the Casimir functionals $C_{Φ, Ψ}$ . Thus, the introduction of SALT into the TRSW equations preserves the Lie–Poisson bracket in their Hamiltonian formulation and thereby preserves their geometric interpretation as coadjoint motion Holm et al. (2009).

III. BALANCED INTERPRETATIONS OF TRSW

There exist several approximations of the rotating shallow water (RSW) equations, the most famous one being the quasi-geostrophic (QG) approximation. By assuming the motion to take place in a particular scaling regime, it can be shown that the largest component of the velocity field, called the geostrophic velocity field, is determined by a diagnostic equation, rather than a prognostic equation. The QG approximation is a small perturbation around this geostrophic velocity field. There exists an intermediate model which is more accurate than QG, but is still an approximation of RSW. In this section, we will derive the thermal geostrophic balance by identifying the correct scaling regime and use asymptotic expansions to simplify the TRSW equations. Next, we will show that the thermal rotating shallow water equations can be approximated geometrically to derive a class of equations which was first proposed by Eliassen (1949) and made into a variational theory by Salmon (1983), where it is called L1. The Lagrangian corresponding to the equations proposed by Eliassen can be obtained via two approaches. The first approach involves the Helmholtz decomposition and the second approach follows Allen and Holm (1996). The methods of Allen and Holm (1996) will be applied in the Euler–Poincaré framework to derive the corresponding equations of motion. Finally, the stochastic thermal L1 (TL1) equations will be derived via the stochastic Euler–Poincaré theorem.

A. Thermal geostrophic balance

To obtain the thermal geostrophic balance relation, we return to the nondimensional deterministic TRSW equations for the Eulerian horizontal vector velocity $u (x, t)$ ⁠, thickness $η (x, t)$ ⁠, buoyancy $b (x, t)$ ⁠, and free surface elevation $α ζ = η - h$ ⁠, with mean depth $h (x)$ ⁠, given in (2.1) by

\begin{matrix} \frac{D}{D t} u + \frac{1}{Ro} f \hat{z} \times u = - \frac{α}{{Fr}^{2}} \nabla ((1 + s b) ζ) + \frac{s}{2 {Fr}^{2}} (ζ - h) \nabla b, \\ \frac{\partial η}{\partial t} + \nabla \cdot (η u) = 0, \frac{D}{D t} b = 0, \end{matrix}

(3.1)

with boundary conditions in (2.2). In order to find an asymptotic balance among these equations, a number of assumptions are necessary. First, we assume that the free surface elevation $α ζ = η - h$ is small, that is $α ≪ 1$ ⁠. Second, in line with the Boussinesq approximation in three dimensional fluids, we assume that the buoyancy stratification is also small, meaning that the stratification parameter $s ≪ 1$ ⁠. Third, we assume that the all dimensionless numbers have equal magnitude

O (α) = O (s) = O (Ro) = O (Fr) .

(3.2)

Now, in the QG approximation one also assumes the gradients of bathymetry and Coriolis parameter are small, of order $O (Ro)$ ⁠. Taken together, this amounts to the following set of assumptions:

\begin{matrix} f (x) = 1 + Ro f_{1} (x), \\ h (x) = 1 + Ro h_{1} (x) . \end{matrix}

(3.3)

The beta plane approximation is included in the notation used in (3.3), when $β f_{0}^{- 1} = O (Ro)$ and $f_{1} (x) = y$ ⁠. These assumptions were also made in derivation of this model in Holm and Luesink (2019) and they are sufficient to derive the TRSW balance relation, which we will show now. Note that since all dimensionless numbers have the same magnitude as in (3.2), we can continue the analysis with a single small parameter $ϵ ≪ 1$ ⁠. First, we multiply the momentum equation in (3.1) by ϵ, then we substitute the assumptions (3.3) into the momentum equation to find

\begin{matrix} \hat{z} \times u + ϵ \frac{D}{D t} u + ϵ f_{1} \hat{z} \times u & = - \nabla ζ_{1} - \frac{1}{2} \nabla b - ϵ b \nabla ζ - \frac{ϵ}{2} (ζ + h_{1}) \nabla b_{1} . \end{matrix}

(3.4)

Thus, Eq. (3.4) implies the following relation

\hat{z} \times u = - \nabla ζ - \frac{1}{2} \nabla b + O (ϵ) .

(3.5)

By operating with $\hat{z} \times$ on the TRSW balance relation (3.5), we find the defining expression for the divergence-free thermal geostrophic velocity field, denoted as $u_{T G}$ ⁠,

u_{T G} = \hat{z} \times \nabla (ζ + \frac{1}{2} b) = : \hat{z} \times \nabla ψ .

(3.6)

Here, $ψ (x, t)$ plays the role of the stream function for the divergence-free leading order thermal geostrophic velocity vector field $u_{T G}$ ⁠. Relation (3.5) allows the total fluid velocity to be represented as the sum of the leading order thermal geostrophic velocity $u_{T G}$ and higher order terms

u = u_{T G} + ϵ υ .

(3.7)

In particular, (3.7) implies that the difference, referred to as the ageostrophic component of the velocity field, satisfies $u - u_{T G} = ϵ υ = O (ϵ)$ ⁠. This decomposition of the velocity field into a leading order divergence-free part plus higher order parts is similar to the Helmholtz decomposition, except the divergence-free component is allowed to have both a leading order part and a higher order part, while the rotation-free component has only higher order parts.

B. Moving into the balanced frame

One may transform the thermal rotating shallow water equations in (2.1) into a time-dependent local frame moving with the thermal geostrophically balanced velocity, $u_{T G} (x, t)$ ⁠, by inserting the decomposition (3.7) into the Lagrangian for the TRSW equations (2.16)

\begin{matrix} ℓ_{TRSW} (u, b, η) = \int_{D} \frac{1}{2} η | u |^{2} + \frac{1}{ϵ} η u \cdot R - \frac{1}{2 ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) d x d y, \\ = \int_{D} ϵ η (\frac{ϵ}{2} | υ |^{2} + υ \cdot u_{T G} + \frac{1}{ϵ} υ \cdot R) + η (\frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} u_{T G} \cdot R), \\ - \frac{1}{2 ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) d x d y . \end{matrix}

(3.8)

One can apply Hamilton's variational principle to the Lagrangian (3.8) $0 = δ S$ with $S = \int_{t_{1}}^{t_{2}} ℓ_{TRSW} d t$ ⁠. We substitute the velocity decomposition in (3.7) to define the variation $δ u = δ u_{T G} + ϵ δ υ$ with

δ u_{T G} = \hat{z} \times \nabla (δ η + \frac{1}{2} δ b) .

(3.9)

Hamilton's principle then yields the Euler–Poincaré equation (2.12) in the form

\frac{\partial}{\partial t} (\frac{1}{η} \frac{δ ℓ}{δ υ}) + (u \cdot \nabla) (\frac{1}{η} \frac{δ ℓ}{δ υ}) + (\nabla u) \cdot (\frac{1}{η} \frac{δ ℓ}{δ υ}) = - \frac{1}{η} \frac{δ ℓ}{δ b} \nabla b + \nabla \frac{δ ℓ}{δ η} .

(3.10)

Thus, the relative motion equation for TRSW dynamics in the frame moving with the thermal geostrophic balance velocity $u_{T G} (x, t)$ in (3.6) keeps its Euler–Poincaré form (2.12). Upon eliminating $\partial_{t} u_{T G}$ in (3.10) by using the advection equations for $(b, η)$ in (3.1), the system closes and thereby transforms the TRSW equations into the new variables $(υ, b, η)$ in the reference frame moving with velocity $u_{T G} (x, t)$ ⁠.

1. Stationary thermal geostrophic balance as “mean dynamic topography (MDT)”

To a good approximation, much of upper ocean dynamics is well-approximated by a mean dynamic topography (MDT), which is monitored continuously with in situ instruments and satellites, see, e.g., Maximenko et al. (2009). Ocean dynamics is then envisioned as time-dependent variations in the steady moving frame of the MDT. To apply this idea to TRSW dynamics, we envision TRSW dynamics as taking place in the moving reference frame defined by a time-independent mean thermal geostrophic velocity ${\bar{u}}_{T G} (x)$ ⁠. In this situation, the Kelvin theorem 2.18 for deterministic TRSW derived from Eq. (3.10) takes the following form.

Theorem 3.1 (Kelvin theorem for deterministic TRSW in a stationary balanced frame). The deterministic TRSW equations (2.1) imply the following Kelvin circulation law in a stationary thermal geostrophically (TG) balanced frame moving with time-independent velocity ${\bar{u}}_{T G} (x)$ ⁠,

\frac{d}{d t} \oint_{c (υ)} (ϵ υ + u_{T G} (x) + \frac{1}{ϵ} R (x)) \cdot d x = \frac{1}{2 ϵ} \oint_{c (υ)} (ϵ ζ - h) \nabla b \cdot d x = \frac{1}{2 ϵ} \int \int_{\partial S = c (υ)} \hat{z} \cdot \nabla (ϵ ζ - h) \times \nabla b d x d y,

(3.11)

where $c (υ)$ is any closed loop moving with horizontal fluid velocity $υ (x, t)$ relative to the frame of motion whose velocity is ${\bar{u}}_{T G} (x) + ϵ^{- 1} R (x)$ in two horizontal dimensions.0

Remark 3.2. Thus, in the mean thermal geostrophic balance frame, the frame velocity ${\bar{u}}_{T G} (x)$ simply adds another contribution to the momentum per unit mass. In turn, this contributes an additional “Coriolis” force in the dynamics of the relative velocity $ϵ υ = u - u_{T G}$ . The corresponding SALT version in this case would simply replace $u (x, t)$ by $υ (x, t)$ as the drift velocity of the stochastic vector field $d χ_{t}$ defined in (2.23).

In Secs. III C and III D, we will consider the thermal versions of two of the classic GFD approximations of RSW developed previously in the absence of buoyancy. Namely, we will consider thermal versions of the Eliassen approximation and the quasi-geostrophic approximation.

C. The Eliassen approximation

The starting point in deriving the thermal Eliassen approximation is the Lagrangian for the thermal rotating shallow water equations (2.16)

ℓ_{TRSW} = \int \frac{1}{2} η | u |^{2} + \frac{1}{Ro} η u \cdot R - \frac{1}{2 {Fr}^{2}} (1 + s b) (η^{2} - 2 η h) d x d y .

(3.12)

Since the thermal geostrophic velocity field (3.6) is divergence-free, it is useful to transform the velocity variables inside the Lagrangian to vorticity and divergence by using the Helmholtz decomposition. The forward transformation is $(u, η, b) \mapsto (ω, D, η, b)$ ⁠, which amounts to

\begin{matrix} ω = \hat{z} \cdot \nabla \times u, \\ D = \nabla \cdot u, \\ η = η, \\ b = b . \end{matrix}

(3.13)

The inverse transformation is unique if the kernel of the Laplacian is trivial, which is the same as saying that there are no harmonic functions for the domain $D$ for the boundary conditions on $\partial D$ ⁠. In this case, the inverse transformation $(ω, D, η, b) \mapsto (u, η, b)$ is given by

\begin{matrix} u = \hat{z} \times \nabla Δ^{- 1} ω + \nabla Δ^{- 1} D, \\ η = η, \\ b = b . \end{matrix}

(3.14)

The inverse transformation (3.14) uniquely defines the vector u in terms of its divergence and curl. The inverse of the Laplacian can be interpreted in terms of the appropriate Green's function in two dimensions. Boundary conditions need to be dealt with carefully when taking the Green's function approach. Assuming that this is the case, the symbol $Δ^{- 1}$ denotes the correct Green's function. For doubly periodic boundary conditions, the Green's function for the Laplacian takes the form

Δ^{- 1} F = - \frac{1}{2 π} \int ln | | x - x' | | F (x') d x' .

(3.15)

Changing variables using the classical Helmholtz decomposition leads to the following formulation of the Lagrangian for thermal rotating shallow water

\begin{matrix} ℓ_{TRSW} = \int \frac{1}{2} η {| \nabla Δ^{- 1} ω |}^{2} + η J (Δ^{- 1} ω, Δ^{- 1} D) + \frac{1}{2} η | \nabla Δ^{- 1} D |^{2} \\ + \frac{1}{Ro} η (\hat{z} \times \nabla Δ^{- 1} ω + \nabla Δ^{- 1} D) \cdot R - \frac{1}{2 {Fr}^{2}} (1 + s b) (η^{2} - 2 η h) d x d y . \end{matrix}

(3.16)

It should be noted that the vorticity and divergence are not orthogonal in a weighted L² space, so the Jacobian term $\int η J (Δ^{- 1} ω, Δ^{- 1} D) d x d y \neq 0$ ⁠. The Jacobian term does vanish in the standard L² space, in which $\int J (Δ^{- 1} ω, Δ^{- 1} D) d x d y = 0$ ⁠. By means of the ordering $O (α) = O (s) = O (Ro) = O (Fr)$ made in (3.4) and the decomposition of u into a thermal geostrophic part and a higher order part (3.7), one can take the following asymptotic expansions for the vorticity and the divergence

\begin{matrix} ω = ω_{0} + ϵ ω_{1} + o (ϵ), \\ D = ϵ D_{1} + o (ϵ) . \end{matrix}

(3.17)

The thermal geostrophic balance implies a decomposition of the velocity field in terms a leading order divergence-free component and a higher order general component. This means that one can identify ω₀ with the curl of the thermal geostrophic velocity field (3.6), but it also means that one must keep a higher order vorticity term around. Substituting (3.17) into the Lagrangian yields

\begin{matrix} ℓ_{TRSW} = \int \frac{1}{2} η {| \nabla Δ^{- 1} (ω_{0} + ϵ ω_{1}) |}^{2} \\ + η J (Δ^{- 1} (ω_{0} + ϵ ω_{1}), Δ^{- 1} ϵ D_{1}) + \frac{1}{2} η | \nabla Δ^{- 1} ϵ D_{1} |^{2} \\ + \frac{1}{ϵ} η (\hat{z} \times \nabla Δ^{- 1} (ω_{0} + ϵ ω_{1}) + \nabla Δ^{- 1} ϵ D_{1}) \cdot R \\ - \frac{1}{2 ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) d x d y + o (ϵ) . \end{matrix}

(3.18)

By expanding and collecting all terms that are of higher order than $O (1)$ ⁠, the Lagrangian can be written as

\begin{matrix} ℓ_{TRSW} = \int \frac{1}{2} η {| \nabla Δ^{- 1} ω_{0} |}^{2} + \frac{1}{ϵ} η (\hat{z} \times \nabla Δ^{- 1} (ω_{0} + ϵ ω_{1}) + ϵ \nabla Δ^{- 1} D_{1}) \cdot R - \frac{1}{2 ϵ^{2}} b (η^{2} - 2 η h) d x d y + o (1) . \end{matrix}

(3.19)

We now use the fact that the velocity field also decomposes into a thermal geostrophic part and an ageostrophic part and apply the inverse transformation to recover the original fluid variables. This yields

ℓ_{T L 1} = \int \frac{1}{2} η | u_{T G} |^{2} + \frac{1}{ϵ} η u \cdot R - \frac{1}{2 ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) d x d y .

(3.20)

The subscript TL1 refers to the thermal L1 model, which is an extension of Salmon (1983) to include horizontal variations in buoyancy and bathymetry. However, at this stage we do not have enough information to execute the variational principle. To obtain the required information, we will introduce a higher order term which completes the Lagrangian, by enabling us to vary with respect to the full velocity field u, interpreted as a Lagrange multiplier

ℓ_{T L 1} = \int u \cdot (η u_{T G} + \frac{1}{ϵ} η R) - \frac{1}{2} η | u_{T G} |^{2} - \frac{1}{2 ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) d x d y + O (ϵ) .

(3.21)

Equivalently, one can truncate the Lagrangian for TRSW rewritten in the balanced frame (3.8) at $O (1)$ to obtain this Lagrangian. This emphasizes the fact that a component of the velocity field which can be expressed in terms of the other variables in the problem can be used to change the reference frame. At this point, several important questions arise. How does one take variations of this Lagrangian? Can the equation for the Lagrange multiplier u be found? To answer these questions, we use the methods of Allen and Holm (1996).

D. The Allen–Holm approach

The Allen and Holm (1996) approach is based on the following observation. A Lagrangian $ℓ$ leads to the corresponding Hamiltonian $ℏ$ via the Legendre transform (which is assumed to be invertible for the given $ℓ$ ⁠)

ℏ (m, η) = \int (m - \frac{δ ℓ}{δ u}) \cdot u + \bar{E} (u, η, b, \nabla η, \nabla b, etc .),

(3.22)

where $\bar{E}$ can be interpreted as the energy density. The momentum density m is given in terms of the other fluid variables by the condition $δ ℏ / δ u = 0$ ⁠, where $ℏ$ is the Hamiltonian defined by the Legendre transform in (3.22). In the Legendre transform, the fluid velocity u appears as a Lagrange multiplier which enforces the relation of m to the other fluid variables as a dynamically preserved constraint. This definition is usually taken for granted, but in what follows, we shall model the momentum density as a prescribed function of the other fluid variables. This means that we will define

m = \bar{m} (η, b, \nabla η, \nabla b, etc .) = : \bar{m} [η, b] .

(3.23)

In this type of modeling, it is necessary to have the explicit enforcement of the momentum definition (3.23), both as a constraint as well as a means of determining the fluid velocity for the model by using Lagrange multipliers. We rearrange the Lagrangian in (3.12) using $\bar{m}$ as the momentum density, defined by

\bar{m} : = \frac{δ ℓ_{TRSW}}{δ u} = η u + \frac{1}{Ro} η R .

(3.24)

The Lagrangian in (3.12) can then be written as

ℓ_{TRSW} = \int_{D} \bar{m} \cdot u - \frac{1}{2} η | u |^{2} - \frac{1}{2 {Fr}^{2}} (1 + s b) (η^{2} - 2 η h) d x d y .

(3.25)

Substitution of the decomposition of the velocity field into its geostrophic and ageostrophic components in (3.7) with $u_{T G}$ defined in (3.6) into the Lagrangian (3.25) now leads, without approximation, to the following Lagrangian, which is linear in the velocity u

ℓ_{TRSW} = \int_{D} u \cdot (η u_{T G} + Ro η u_{A} + \frac{1}{Ro} η R) - \frac{1}{2} η | u_{T G} + Ro u_{A} |^{2} - \frac{1}{2 {Fr}^{2}} (1 + s b) (η^{2} - 2 η h) d x d y .

(3.26)

In line with the ordering scheme $O (α) = O (s) = O (Ro) = O (Fr)$ ⁠, we formulate the Lagrangian in terms of a single parameter ϵ. When $ϵ ≪ 1$ ⁠, one could simply drop the u_A terms in the Lagrangian to obtain

ℓ_{T L 1} = \int_{D} u \cdot (η u_{T G} + \frac{1}{ϵ} η R) - \frac{1}{2} η | u_{T G} |^{2} - \frac{1}{2 ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) d x d y .

(3.27)

One recalls that this is the Lagrangian obtained in (3.21) by using the Helmholtz decomposition to decompose u into vorticity and divergence.

Remark 3.3. Since we will use (3.6) as our definition for $u_{T G}$ , we should keep in mind that according to strict asymptotics the potential energy term should also be expanded using the same assumptions that led to (3.6). By keeping the Lagrangian in the form (3.27), we have included higher order terms, but not all of them, since we have truncated the kinetic energy. In terms of strict asymptotics, this means that we do not have a balance among terms in the Lagrangian. A benefit of not expanding the potential energy at this stage is that the variational derivatives can be taken in the usual way and are thus closer to the variational derivatives of the TRSW system.

The approximate Lagrangian (3.27) is also linear in the velocity u, since it has the form

ℓ_{T L 1} = \int \bar{m} [η, b] \cdot u - \bar{E} [η, b] d x d y,

(3.28)

with

\bar{m} [η, b] = \frac{δ ℓ_{T L 1}}{δ u} = \frac{1}{ϵ} η R + η u_{T G} and \bar{E} [η, b] = \frac{1}{2} η | u_{T G} |^{2} + \frac{1}{2 ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) .

(3.29)

Note that $u_{T G}$ can be expressed in terms of η and b. At this stage, in Allen and Holm (1996), the next step after having obtained the Lagrangian in the form above would have been to take the Legendre transformation and obtain the Hamiltonian. Then, by requiring the first variation of the Hamiltonian to vanish, one would obtain the equations of motion. In Holm et al. (1998) it was shown, however, that one can obtain the same equations of motion by applying the variational principle on the Lagrangian side, by means of the Euler–Poincaré theorem. The first variation of the Lagrangian is given by

\begin{matrix} δ ℓ_{T L 1} = \int_{D} (η u_{T G} + \frac{1}{ϵ} η R) \cdot δ u + (η (u - u_{T G})) \cdot δ u_{T G} \\ + (u_{T G} \cdot u + \frac{1}{ϵ} u \cdot R - \frac{1}{2} | u_{T G} |^{2} - \frac{1}{ϵ} (1 + ϵ b) ζ) δ η \\ + (- \frac{1}{2 ϵ} (η^{2} - 2 η h)) δ b d x d y . \end{matrix}

(3.30)

From the definition of $u_{T G}$ in (3.6), we now substitute

δ u_{T G} = \frac{1}{ϵ} \hat{z} \times \nabla δ η + \frac{1}{2} \hat{z} \times \nabla δ b .

(3.31)

Integration by parts in (3.30) then yields

\begin{matrix} δ ℓ_{T L 1} = \int_{D} (η u_{T G} + \frac{1}{ϵ} η R) \cdot δ u + (u_{T G} \cdot u + \frac{1}{ϵ} R \cdot u - \frac{1}{2} | u_{T G} |^{2} - \frac{1}{ϵ} (1 + ϵ b) ζ - \frac{1}{ϵ} \hat{z} \cdot \nabla \times (η (u - u_{T G}))) δ η \\ + (- \frac{1}{2 ϵ} (η^{2} - 2 η h) - \frac{1}{2} \hat{z} \cdot \nabla \times (η (u - u_{T G}))) δ b d x d y \\ - \oint_{\partial D} \hat{z} \times (η (u - u_{T G})) (\frac{1}{ϵ} δ η + \frac{1}{2} δ b) \cdot d x . \end{matrix}

(3.32)

The integral over the boundary vanishes provided that the ageostrophic velocity field has no tangential component on the boundary. This boundary condition is satisfied since the ageostrophic velocity field can be represented by a velocity potential, as in (3.14). By means of the Euler–Poincaré theorem, we find the equations of motion given in a form first proposed by Eliassen Eliassen (1949) as

\begin{matrix} \frac{\partial}{\partial t} u_{T G} - u \times (\nabla \times (u_{T G} + \frac{1}{ϵ} R)) + \nabla (\frac{1}{ϵ} (1 + ϵ b) ζ + \frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} B) - (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \nabla b = 0, \\ \frac{\partial}{\partial t} η + \nabla \cdot (η u) = 0, \frac{\partial}{\partial t} b + u \cdot \nabla b = 0. \end{matrix}

(3.33)

The notation B in the first of these equations is defined by

B : = \hat{z} \cdot \nabla \times (η (u - u_{T G})) .

(3.34)

The function B keeps track of effects that are generated by higher order vorticity terms, since B includes the curl of the velocity difference $u - u_{T G}$ ⁠, which is of order $O (ϵ)$ ⁠. These higher order vorticity terms will contribute in the Kelvin circulation theorem 3.4 arising from Eqs. (3.33). Equations (3.33) extend Salmon's L1 model Salmon (1983) to include horizontal buoyancy variations and bottom topography. The boundary conditions carry over from thermal rotating shallow water and are given by

\hat{n} \cdot u = 0 and \hat{n} \times \nabla b = 0 on the boundary \partial D .

(3.35)

Having been derived from the Euler–Poincaré variational principle Holm et al. (1998), the deterministic TL1 equations in (3.33) satisfy the following Kelvin–Noether circulation theorem.

Theorem 3.4 (Kelvin theorem for the deterministic TL1 model). The deterministic TL1 equations (3.33) imply the following Kelvin circulation law

\begin{matrix} \frac{d}{d t} \oint_{c (u)} (u_{T G} + \frac{1}{ϵ} R) \cdot d x = \oint_{c (u)} (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \nabla b \cdot d x \\ = \int \int_{\partial S = c (u)} \hat{z} \cdot \nabla (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \times \nabla b d x d y . \end{matrix}

(3.36)

Proof. This result follows from the Kelvin–Noether theorem 2.6 for Euler–Poincaré fluid equations.

Remark 3.5. The Kelvin circulation theorem 3.4 for the deterministic TL1 model (3.33) implies that the misalignment of the horizontal gradients of the free surface elevation and the bathymetry with the horizontal gradient of the buoyancy generates circulation. This result is similar to the corresponding Kelvin circulation theorem 2.8 for the deterministic thermal rotating shallow water (TRSW) model. An additional contribution to the generation of circulation relative to theorem 2.8 is made by the misalignment of the gradient of the quantity $(B / η)$ defined in (3.34) with the gradient of the buoyancy. This additional contribution is due to misalignment of the horizontal gradients of the ageostrophic vorticity and the buoyancy.

The potential vorticity q for TL1 is defined as

q = \frac{1}{η} (ω + \frac{1}{ϵ} f) = \frac{1}{η} (Δ (ζ + \frac{1}{2} b) + \frac{1}{ϵ} + f_{1}) .

(3.37)

Note that the potential vorticity q in (3.37) contains a term in the Coriolis parameter which is order $O (ϵ^{- 1})$ ⁠. This feature will become important in the asymptotic expansion of q later, in deriving the thermal QG model at the beginning of Sec. V. In the presence of buoyancy, the evolution of potential vorticity along Lagrangian fluid trajectories is not conserved. Instead, potential vorticity is generated, as indicated in the circulation theorem (3.36) via misalignment of gradients in

\frac{\partial}{\partial t} q + u \cdot \nabla q = \frac{1}{η} \hat{z} \cdot \nabla (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \times \nabla b .

(3.38)

Although the potential vorticity is not conserved along Lagrangian fluid trajectories, the TL1 equations do preserve energy, as well as an infinity of integral conservation laws involving buoyancy and potential vorticity.

1. Conservation laws for deterministic TL1

The deterministic TL1 equations (3.33) conserve the energy

E_{T L 1} (u_{T G}, η, b) = \frac{1}{2} \int_{D} η | u_{T G} |^{2} + \frac{1}{ϵ^{2}} (1 + ϵ b) (η^{2} - 2 η h) d x d y .

(3.39)

Equations (3.33) also conserve an infinity of integral conservation laws, determined by two arbitrary differentiable functions of buoyancy $Φ (b)$ and $Ψ (b)$ as

C_{Φ, Ψ} = \int_{D} η Φ (b) + ϖ Ψ (b) d x d y = \int_{D} (Φ (b) + q Ψ (b)) η d x d y,

(3.40)

where ω and q are defined in Eq. (3.37). Notice that this family of integral conservation laws for the TL1 equations has the same form as the family of integral conserved quantities for the TRSW equations, defined in (2.22). The proof that $C_{Φ, Ψ}$ is conserved in time, for any choice of differentiable $Φ$ and $Ψ$ ⁠, follows from a direct calculation involving the boundary conditions (3.35). The integral conserved quantities in Eqs. (5.15) and (3.40) are associated with smooth transformations which leave invariant the Eulerian fluid quantities in the Lagrangian. As with the TRSW equations, upon introducing stochasticity via the Euler–Poincaré theorem 2.10, the latter conservation laws persist. However, energy conservation does not persist because the stochastic Lagrangian depends explicitly on time through the Brownian motion.

In order to use the TL1 equations (3.33) as a predictive model, one needs to be able to determine the Lagrange multiplier u from the other variables in the model. This will be our next task.

2. Determining the Lagrange multiplier

By operating with $\hat{z} \times$ on the momentum equation in (3.33) and using the definition of $u_{T G}$ in (3.6), we obtain

- \frac{\partial}{\partial t} \nabla (\frac{1}{ϵ} η + \frac{1}{2} b) - η q u = \hat{z} \times \nabla (\frac{1}{ϵ} (1 + ϵ b) ζ + \frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} B) - (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \hat{z} \times \nabla b .

(3.41)

The first term above follows from using the definition of $u_{T G}$ ⁠, by noting that the bathymetry has no time derivative. This allows us to rewrite $u_{T G}$ in terms of η and b. Taking the time derivative through the gradient in (3.41) allows us to use the continuity equation and the advection equation to obtain an elliptic equation. Substituting the definition for B in (3.34) then leads to the following equation which determines the Lagrange multiplier u

\begin{matrix} \frac{1}{ϵ} \nabla (\nabla \cdot η u) + \frac{1}{2} \nabla (u \cdot \nabla b) + \hat{z} \times \nabla (\frac{1}{ϵ} (1 + ϵ b) ζ + \frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} \hat{z} \cdot \nabla \times (η (u - u_{T G}))) - (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} \hat{z} \cdot \nabla \times (η (u - u_{T G}))) \hat{z} \times \nabla b - η q u = 0 . \end{matrix}

(3.42)

Before going to the general case, let us consider the case in which the horizontal gradient of buoyancy vanishes.

3. No horizontal buoyancy gradients

In this case, Eq. (3.42) reduces to

\frac{1}{ϵ} \nabla (\nabla \cdot η u) + \hat{z} \times \nabla (\frac{1}{ϵ} ζ + \frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} \hat{z} \cdot \nabla \times (η (u - u_{T G}))) - η q u = 0 .

(3.43)

This is the diagnostic partial differential equation (PDE) used to determine u in Salmon (1983) when variations in bathymetry are absent and it is identical to Eq. (3.16) in Allen and Holm (1996). After applying the identity

Δ u = \nabla^{⊥} (\nabla^{⊥} \cdot u) + \nabla (\nabla \cdot u)

(3.44)

in Eq. (3.43), we can rewrite the diagnostic equation (3.42) in simpler form. Here we have used the perpendicular “ $⊥$ ” notation for brevity, see remark 2.3. The Laplacian identity (3.44) implies that the equation which determines u is a linear nonautonomous elliptic partial differential equation (PDE), given by

\frac{1}{ϵ} Δ (η u) - η q u = - \nabla^{⊥} (\frac{1}{ϵ} ζ + \frac{1}{2} | u_{T G} |^{2} - \frac{1}{ϵ} \nabla^{⊥} \cdot η u_{T G}) .

(3.45)

Note that the coefficient ηq in (3.45) is the total vorticity, since $η q u = (\frac{1}{ϵ} + f_{1} + Δ ζ) u$ ⁠. The solution behavior of the elliptic equation (3.45) for the quantity $η u$ depends on the sign of the potential vorticity, q, in the following three cases:

q > 0. The equation for $η u$ is a screened Poisson equation.
q < 0. The equation for $η u$ is an inhomogeneous Helmholtz equation.
q = 0. The equation for $η u$ is a Poisson equation.

In the situation being considered at the moment, there are no horizontal buoyancy gradients. Consequently, the potential vorticity q is preserved along Lagrangian particle trajectories and does not change sign during the calculation. This means that such changes in the solution behavior of (3.45) do not occur. Thus, the “equator,” where f changes sign, acts as a boundary between the “northern and southern hemispheres,” in the absence of horizontal buoyancy gradients. In this case, Lagrangian particles which start in the northern hemisphere stay in the northern hemisphere, because their potential vorticity is conserved in the absence of horizontal buoyancy gradients. The case with nonzero horizontal buoyancy gradients is the general case, which we will discuss now.

4. General case

When horizontal gradients of buoyancy are nonzero, Eq. (3.45) for the determination of the Lagrange multiplier u becomes considerably more extensive

\begin{matrix} \frac{1}{ϵ} Δ (η u) + \frac{1}{2} \nabla (u \cdot \nabla b) - \frac{1}{2} (\nabla^{⊥} \cdot u) \nabla^{⊥} b - \frac{1}{2 η} (u \cdot \nabla^{⊥} η) \nabla^{⊥} b - η q u = - \frac{1}{ϵ} \nabla^{⊥} ((1 + ϵ b) ζ) + \frac{1}{2 ϵ} (ϵ ζ - h) \nabla^{⊥} b - \frac{1}{2} \nabla^{⊥} | u_{T G} |^{2} - \frac{1}{2 η} (\nabla^{⊥} \cdot η u_{T G}) \nabla^{⊥} b + \frac{1}{ϵ} \nabla^{⊥} \nabla^{⊥} \cdot (η u_{T G}) . \end{matrix}

(3.46)

The coefficient for the u term now has an additional contribution from the perpendicular gradient of the buoyancy, which changes the conditions for the type of PDE. The zeroth order terms in the presence of horizontal buoyancy gradients indicate that the equator is no longer a stationary boundary between the northern and southern hemispheres. Indeed, the perpendicular gradients of buoyancy in combination with perpendicular gradients of the depth have removed the “equatorial boundary.” Likewise, when the horizontal gradient of buoyancy is included, the potential vorticity is not conserved along Lagrangian particle trajectories. Moreover, the effects of the first order terms at this point remain unexamined. At this point, we shall defer further discussion of these elliptic equations until Sec. V and leave the discussion of the interpretation of the effects of horizontal buoyancy gradients on the solution behavior of the elliptic equation (3.46) for the quantity $η u$ for the TL1 model as an open problem.

Before substituting the asymptotic expansions and truncating the equations to achieve thermal geostrophic balance in terms of strict asymptotics in Sec. V, we will first derive the stochastic thermal L1 equations.

IV. THE ELIASSEN APPROXIMATION OF STOCHASTIC TRSW

The equation sets for the deterministic and stochastic TRSW models in Sec. II and the deterministic TL1 model in Sec. III D have all been derived in the variational framework of the Euler–Poincaré theorem introduced in Sec. II B. The corresponding Kelvin circulation laws for each of these theories follow from their Kelvin–Noether theorem 2.6, proved in Sec. II C. Let us now derive the stochastic version of the TL1 equations and their corresponding Kelvin circulation law by following the SALT formulation in the Euler–Poincaré variational framework. To do so, we first investigate the balance relation in the presence of stochasticity.

A. Stochastic thermal geostrophic balance

To obtain the stochastic thermal geostrophic balance, we start from the TRSW equations with SALT, given in (2.30) by

\begin{matrix} d u + (d χ_{t} \cdot \nabla) u + \sum_{k} (\nabla ξ_{k}) \cdot u ○ d W_{t}^{k} = - \frac{α}{{Fr}^{2}} \nabla ((1 + s b) ζ) d t + \frac{s}{2 {Fr}^{2}} (α ζ - h) \nabla b d t \\ - \frac{1}{Ro} f \hat{z} \times d χ_{t} - \frac{1}{Ro} \sum_{k} \nabla (ξ_{k} \cdot R) ○ d W_{t}^{k}, \\ d η + \nabla \cdot (η d χ_{t}) = 0, \\ d b + (d χ_{t} \cdot \nabla) b = 0, \end{matrix}

(4.1)

with boundary conditions in (2.31). We recall the assumptions that led to deterministic thermal geostrophic balance. As before, we assume that the magnitudes of the dimensionless numbers in the problem is $O (α) = O (s) = O (Ro) = O (Fr)$ ⁠. This asymptotic regime allows us to continue with a single small parameter, ϵ. So, we formulate (3.3) as

\begin{matrix} f (x) = 1 + ϵ f_{1} (x), \\ h (x) = 1 + ϵ h_{1} (x), \\ R (x) = R_{0} (x) + ϵ R_{1} (x), \end{matrix}

(4.2)

where the additional relations $\hat{z} \cdot \nabla \times R_{0} = 1$ and $\hat{z} \cdot \nabla \times R_{1} = f_{1}$ hold. Upon substituting the asymptotic expansions (4.2) into the stochastic TRSW equations (4.1) and collecting all terms of $O (ϵ^{- 1})$ ⁠, we find

(\nabla ζ + \frac{1}{2} \nabla b + \hat{z} \times u) d t + \sum_{k} (\hat{z} \times ξ_{k} + \nabla (ξ_{k} \cdot R_{0})) ○ d W_{t}^{k} = 0 .

(4.3)

The drift part of this stochastic partial differential equation is the deterministic thermal geostrophic balance (3.5) and the diffusion part provides us with a relation between the noise amplitude $ξ_{k}$ and the vector potential for the Coriolis parameter

\begin{matrix} u = \hat{z} \times \nabla (ζ + \frac{1}{2} b) + O (ϵ), \\ ξ_{k} = \hat{z} \times \nabla (ξ_{k} \cdot R_{0}) + O (ϵ) . \end{matrix}

(4.4)

Since the Brownian motions are assumed to be independent, (4.4) needs to be satisfied for each k. We can identify the thermal geostrophic balance velocity field as $u_{T G} = \hat{z} \times \nabla (ζ + \frac{1}{2} b)$ and expand the velocity field u as in the deterministic case

u = u_{T G} + O (ϵ),

(4.5)

and following the same reasoning, the $ξ_{k}$ can be expanded as

ξ_{k} = ξ_{T G_{k}} + O (ϵ) .

(4.6)

We can now investigate the stochastic thermal L1 model.

B. Stochastic TL1

The equations governing the stochastic TL1 model are obtained in the SALT formulation by applying the stochastic Euler–Poincaré theorem 2.10 to the TL1 Lagrangian, given by (3.27). This incorporates the deterministic geostrophic balance. We find the stochastic version of the TL1 equations (3.33), given by

\begin{matrix} d u_{T G} - d χ_{t} \times (\nabla \times (u_{T G} + \frac{1}{ϵ} R)) + \nabla (\frac{1}{ϵ} (1 + ϵ b) ζ + \frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} B) d t \\ + \nabla (\sum_{k} ξ_{k} \cdot (u_{T G} + \frac{1}{ϵ} R)) ○ d W_{t}^{k} - (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \nabla b d t = 0 . \\ d η + \nabla \cdot (η d χ_{t}) = 0, \\ d b + d χ_{t} \cdot \nabla b = 0 . \end{matrix}

(4.7)

Here, the function B is defined by

B : = \hat{z} \cdot \nabla \times (η (u - u_{T G})) .

(4.8)

The boundary conditions are given by

\hat{n} \cdot u = 0 and \hat{n} \cdot ξ_{k} = 0 and \hat{n} \times \nabla b = 0 on the boundary \partial D .

(4.9)

The Kelvin–Noether theorem 2.6 for the stochastic TL1 model is given by the following theorem.

Theorem 4.1 (Kelvin theorem for the stochastic TL1 model). The stochastic TL1 equations (4.7) imply the following Kelvin circulation law

\begin{matrix} d \oint_{c (d χ_{t})} (u_{T G} + \frac{1}{ϵ} R) \cdot d x = \oint_{c (d χ_{t})} (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \nabla b \cdot d x d t \\ = \int \int_{\partial S = c (d χ_{t})} \hat{z} \cdot \nabla (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \times \nabla b d x d y d t . \end{matrix}

(4.10)

Proof. The proof follows the pattern of the standard Kelvin–Noether theorem 2.6, modulo an application of the Kunita–Itô–Wentzell theorem which provides the chain rule for the Lie derivatives of differential forms by stochastic vector fields, as proved in de Leon et al. (2020). The loop $c (d χ_{t})$ does not explicitly require the evaluation of a stochastic integral, as it is the push-forward of a stationary loop by the flow which is generated by the vector field $d χ_{t}$ ⁠, see remark 2.7.

The stochastic TL1 equations do not conserve energy due to their explicit dependence on time via the noise. From the Kelvin circulation theorem associated with the stochastic TL1 equations, an evolution equation for potential vorticity can be derived. This equation shows that potential vorticity is not conserved along Lagrangian particle trajectories, but is generated by the effect also present on the right-hand side in the Kelvin circulation theorem 4.10. Recall that the potential vorticity q is defined by

q = \frac{1}{η} (ω + \frac{1}{ϵ} f) = \frac{1}{η} (Δ (ζ + \frac{1}{2} b) + \frac{1}{ϵ} + f_{1}) .

(4.11)

The evolution equation for q is given by

d q + d χ_{t} \cdot \nabla q = \frac{1}{η} \hat{z} \cdot \nabla (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \times \nabla b d t .

(4.12)

Even though potential vorticity is not a Lagrangian invariant, the stochastic TL1 equations have an infinite family of integral conservation laws, given by

C_{Φ, Ψ} = \int_{D} η Φ (b) + η q Ψ (b) d x d y .

(4.13)

The proof for these conservation laws is a direct calculation that uses the boundary conditions (4.9). In order to use (4.7) as a predictive model, one must be able to determine u from the other variables in the model. We can proceed as in the deterministic case and derive an elliptic equation for u.

1. Determining the Lagrange multiplier

By operating with $\hat{z} \times$ on the momentum equation in (4.7) and using the definition of $u_{T G}$ (3.6), we obtain

\begin{matrix} d \nabla (\frac{1}{ϵ} η + \frac{1}{2} b) - η q d χ_{t} = \hat{z} \times \nabla (\frac{1}{ϵ} (1 + ϵ b) ζ + \frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} B) d t + \nabla (\sum_{k} ξ_{k} \cdot (u_{T G} + \frac{1}{ϵ} R)) ○ d W_{t}^{k} \\ - (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \hat{z} \times \nabla b d t . \end{matrix}

(4.14)

We continue by taking the stochastic differential through the gradient and substitute the continuity equation and the advection equation for the buoyancy from (4.7). By using the definition for B in (4.8), we can then use the vector calculus identity $Δ u = \nabla^{⊥} \nabla^{⊥} \cdot u + \nabla \nabla \cdot u$ ⁠. This leads to two linear nonautonomous elliptic partial differential equations, one for the drift part and one for the diffusion part. The elliptic equation for the drift part is given by

\begin{matrix} \frac{1}{ϵ} Δ (η u) + \frac{1}{2} \nabla (u \cdot \nabla b) - \frac{1}{2} (\nabla^{⊥} \cdot u) \nabla^{⊥} b - \frac{1}{2 η} (u \cdot \nabla^{⊥} η) \nabla^{⊥} b - η q u = - \frac{1}{ϵ} \nabla^{⊥} ((1 + ϵ b) ζ) + \frac{1}{2 ϵ} (ϵ ζ - h) \nabla^{⊥} b - \frac{1}{2} \nabla^{⊥} | u_{T G} |^{2} \\ - \frac{1}{2 η} (\nabla^{⊥} \cdot η u_{T G}) \nabla^{⊥} b + \frac{1}{ϵ} \nabla^{⊥} \nabla^{⊥} \cdot (η u_{T G}), \end{matrix}

(4.15)

and the elliptic equation for the diffusion part is given by

\frac{1}{ϵ} \nabla (\nabla \cdot η ξ_{k}) + \frac{1}{2} \nabla (ξ_{k} \cdot \nabla b) + \nabla^{⊥} (ξ_{k} \cdot (u_{T G} + \frac{1}{ϵ} R)) - η q ξ_{k} = 0.

(4.16)

These elliptic equations at leading order provide the deterministic and the stochastic geostrophic balance conditions.

We will not pursue the properties of the stochastic TL1 equations any further here. Instead, we will apply further asymptotic analysis to derive the thermal quasi-geostrophic (TQG) equations from the TL1 model, investigate the properties of their deterministic solutions and then derive the corresponding stochastic TQG equations by using the SALT approach.

V. THERMAL QG MODEL

In Secs. III and IV, we made an approximation to the kinetic energy by assuming that the velocity field can be decomposed into a thermal geostrophic part and a higher order part. This led to the thermal L1 (TL1) equations, given by (3.33) in the deterministic case, and by (4.7) in the stochastic case. These TL1 equations, in turn, can be approximated further to yield the motion equations for thermal quasi-geostrophy (TQG), which will be the subject of this section. We will start with the deterministic TL1 case, continuing from where we stopped in Sec. III and discuss the solution properties of a numerical example of the TQG equations. We will finish this section by looking into the Hamiltonian formulation of the deterministic TQG equations. We will then use the Hamiltonian framework to derive the stochastic TQG equations.

A. Deterministic TQG model

To obtain the thermal quasi-geostrophic model, one could expand the TL1 equations (3.33) and find that the continuity equation features the divergence of u, which can then be solved for by substitution. This approach has the disadvantage that equations are expanded before substitution, which loses accuracy. Instead we follow the derivation for the elliptic equation which determines the Lagrange multiplier u. By operating with $\hat{z} \times$ on the TL1 momentum equation in (3.33) and using the definition of $u_{T G}$ ⁠, we obtain

\frac{\partial}{\partial t} \nabla (\frac{1}{ϵ} η + \frac{1}{2} b) + η q u = \hat{z} \times (\nabla (\frac{1}{ϵ} (1 + ϵ b) ζ + \frac{1}{2} | u_{T G} |^{2} + \frac{1}{ϵ} B) - (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B) \nabla b) .

(5.1)

We now take the divergence of (5.1) and substitute the continuity equation for the depth, which yields

\begin{matrix} \frac{\partial}{\partial t} Δ (\frac{1}{ϵ} η + \frac{1}{2} b) - q \frac{\partial}{\partial t} η + η u \cdot \nabla q = J (\frac{1}{2 ϵ} (ϵ ζ - h) + \frac{1}{2 η} B, b), \end{matrix}

(5.2)

so the potential energy term $\partial η / \partial t$ also appears in the vorticity equation. Equivalently, one can take the two dimensional curl, or the perpendicular divergence, to arrive directly at (5.2). In a moment, we will expand these equations in the asymptotic regime introduced in (3.3) and truncate at $O (1)$ to obtain the thermal quasi-geostrophic equations (TQG). In the asymptotic expansions to follow, it will be helpful to note that the potential vorticity, defined in (3.37), contains a term that is of order $O (ϵ^{- 1})$ ⁠, since it allows for the substitution of ζ.

1. Asymptotic expansion to the deterministic TQG regime

By means of the asymptotic expansions (3.3) which were used to derive an expression for $u_{T G}$ ⁠, one can expand Eq. (5.2) and collect terms of the same order. Truncating at $O (1)$ then leads to the following motion equation for thermal quasi-geostrophy (TQG) Warneford and Dellar (2013) and Zeitlin (2018):

\frac{\partial}{\partial t} (Δ (ζ + \frac{1}{2} b) - ζ) + u_{T G} \cdot \nabla (Δ (ζ + \frac{1}{2} b) + f_{1}) = \frac{1}{2} \hat{z} \times \nabla (ζ - h_{1}) \cdot \nabla b .

(5.3)

Because the depth equation in (3.33) was substituted into Eq. (5.2), the potential energy term $\partial ζ / \partial t$ remains in the vorticity equation (5.3). In the asymptotic expansion, the deterministic buoyancy equation keeps its form, as

\frac{\partial}{\partial t} b + u_{T G} \cdot \nabla b = 0,

(5.4)

and the boundary conditions at this order become

\hat{n} \cdot u_{T G} = 0 and \hat{n} \times \nabla b = 0 on the boundary \partial D .

(5.5)

Equations (5.3) and (5.4) together with the boundary conditions (5.5) form a closed model which approximates the TL1 model (3.33). This completes the present derivation of the thermal quasi-geostrophic (TQG) model, cf. Warneford and Dellar (2013) and Zeitlin (2018).

The vorticity equation (5.3) for TQG can be rewritten in terms of the stream function $ψ : = ζ + \frac{1}{2} b$ and the Jacobian operator $J (ψ, a) : = \hat{z} \cdot \nabla ψ \times \nabla a = u_{T G} \cdot \nabla a$ for a function $a (x)$ to obtain the more compact form

\frac{\partial}{\partial t} (Δ ψ - ψ + b) + J (ψ, Δ ψ + f_{1}) = - \frac{1}{2} J (h_{1}, b) .

(5.6)

Here, we have added and subtracted $\frac{1}{2} b$ in the time derivative in vorticity equation (5.3). Another $\frac{1}{2} b$ contribution comes from the forcing term on the right-hand side, upon replacing the free surface elevation ζ by the stream function ψ and then using the buoyancy equation (5.4), written now as

\frac{\partial}{\partial t} b + J (ψ, b) = 0.

(5.7)

The boundary conditions are

\hat{n} \times \nabla ψ = 0 and \hat{n} \times \nabla b = 0 on the boundary \partial D .

(5.8)

When the bathymetry is flat and the Coriolis parameter is constant, Eq. (5.6) reduces to the TQG equation found in Warneford and Dellar (2013) and Zeitlin (2018). Since the Jacobian operator is zero when the arguments are functionally related, it is possible to write the deterministic TQG equation in terms of a type of potential vorticity variable, which we will call q. This notation forms a close link between QG without buoyancy and TQG, with

q : = Δ ψ - ψ + f_{1} .

(5.9)

The definition of q in (5.9) allows us to formulate the vorticity equation in (5.6) as

\frac{\partial}{\partial t} (q + b) + J (ψ, q) = - \frac{1}{2} J (h_{1}, b) .

(5.10)

The formulation in terms of q as in (5.10) is particularly useful in showing that these equations conserve energy, but also to note that the TQG equations can be related to Rayleigh–Bénard convection.

Remark 5.1 (Rayleigh–Bénard convection). One can write (5.10) in such a way that all terms that depend on the buoyancy variations appear on the right-hand side,

\begin{matrix} \frac{\partial}{\partial t} q + J (ψ, q) = \frac{1}{2} J (ψ, b) + \frac{1}{2} J (ζ - h_{1}, b), \\ \frac{\partial}{\partial t} b + J (ψ, b) = 0. \end{matrix}

(5.11)

This notation reveals that Eqs. (5.11) are reminiscent of the ideal Rayleigh–Bénard convection equations. The Rayleigh–Bénard convection problem in a vertical xz plane, formulated in terms of vorticity and stream function, is given by

\begin{matrix} \frac{\partial}{\partial t} \tilde{q} + J (ψ, \tilde{q}) = α g T_{z}, \\ \frac{\partial}{\partial t} T + J (ψ, T) = 0. \end{matrix}

(5.12)

Here, $\tilde{q} = Δ ψ$ is the vorticity, T is the temperature, ψ is the stream function, g is gravity, and α is the thermal expansion coefficient. For the typical Rayleigh–Bénard convection problem, in the vertical direction there are two solid boundaries and in the horizontal direction, one either uses periodic boundary conditions or solid boundaries. The bottom boundary is heated and the top boundary is cooled, in such a way that the temperature difference is constant. Similar boundary conditions can be established for the thermal quasi-geostrophic equations. The main differences between the two models is that the forcing terms in TQG involve derivatives in every direction, whereas Rayleigh–Bénard convection only involves derivatives of the temperature in the vertical z-direction. Also, TQG is obtained as a model based on thermal geostrophic balance, whereas the Rayleigh–Bénard model does not impose any balance relation.

2. Elliptic equation for TL1

Upon substituting the asymptotic expansions (3.3) and collecting the leading order terms, the elliptic equation for TL1 (3.46) reads

\begin{matrix} \frac{1}{ϵ} Δ u_{T G} - \frac{1}{ϵ} u_{T G} = - \frac{1}{ϵ} \nabla^{⊥} ζ + \frac{1}{2 ϵ} \nabla^{⊥} b + \frac{1}{ϵ} \nabla^{⊥} \nabla^{⊥} \cdot u_{T G} + O (1) . \end{matrix}

(5.13)

By means of the identity $Δ u_{T G} = \nabla \nabla \cdot u_{T G} + \nabla^{⊥} \nabla^{⊥} \cdot u_{T G}$ and using the fact that $u_{T G}$ is divergence-free, we obtain at leading order the definition for $u_{T G} = \nabla^{⊥} (ζ + \frac{1}{2} b)$ ⁠. At the next order, the elliptic equation will provide an expression for u that is consistent with the asymptotic regime.

B. Numerical TQG example

We implemented the TQG equations (5.11) using finite element methods (FEM) for the spatial variables. The FEM algorithm we used is an adaptation of the algorithm formulated in Bernsen et al. (2006) and was implemented using Firedrake (http://www.firedrakeproject.org/index.html), see Rathgeber et al. (2017). In particular, we approximate the vorticity and buoyancy fields in first order discrete Galerkin finite element space, and approximate the stream function in first order continuous Galerkin finite element space. For the time step, we used an optimal third order strong stability preserving Runge–Kutta method, see Gottlieb (2005); Cotter et al. (2019).

Figure 3 shows a snapshot at a certain time taken from a high resolution numerical run of the TQG equations. In this numerical example, we used the following boundary and initial conditions. The domain is ${[0, 2 π]}^{2}$ discretized at a resolution of 512². The boundary conditions are periodic in the vertical direction and walls in the horizontal direction. The parameters and initial conditions are

\begin{matrix} h_{1} (x, y) = cos x + \frac{1}{2} cos 2 x + \frac{1}{3} cos 3 x, \\ f_{1} (x, y) = 0, (f - plane) \\ b (x, y, 0) = - \frac{1}{1 + exp (- x + π)} + \frac{1}{2}, \\ q (x, y, 0) = - exp (- 5 {(y - π)}^{2}), \\ ζ (x, y, 0) = ψ (x, y, 0) - \frac{1}{2} b_{1} (x, y, 0) . \end{matrix}

(5.14)

The stream function $ψ = ζ + b / 2$ is calculated from the potential vorticity q by means of the elliptic problem given in (5.9).

C. Hamiltonian formulation of TQG

In this section, we will investigate the Hamiltonian structure of TQG. The Hamiltonian formulation of baroclinic QG and QG was analyzed in Holm (1986); Holm and Zeitlin (1998), respectively.

1. Conservation laws for deterministic TQG

The deterministic TQG system (5.6) and (5.7) conserves the energy

H_{TQG} (q, b) = - \frac{1}{2} \int_{D} (q - f_{1}) {(Δ - 1)}^{- 1} (q - f_{1}) + h_{1} b) d x d y .

(5.15)

The proof that the TQG equations conserve energy is a direct calculation that requires the boundary conditions (5.8). The deterministic TQG equations also conserve an infinity of integral conservation laws, determined by two arbitrary differentiable functions of buoyancy $Φ (b)$ and $Ψ (b)$ as

C_{Φ, Ψ} = \int_{D} Φ (b) + q Ψ (b) d x d y .

(5.16)

The proofs that the integral quantities H_TQG in (5.15) and $C_{Φ, Ψ}$ in (5.16) are conserved by the TQG equations in (5.7) and (5.10) are both direct calculations which invoke the boundary conditions (5.8). The integral conservation laws $C_{Φ, Ψ}$ for the TQG equations have the same form as the family of integral conserved quantities for the TRSW equations, defined in (2.22), and likewise the integral conserved quantities (3.40) for the TL1 equations with the corresponding potential vorticity variable for TL1 in (3.37). The persistence of these integral conservation laws for the deterministic TQG equations is best explained in terms of their Hamiltonian formulation.

2. Hamiltonian formulation of deterministic TQG

The Hamiltonian form of the deterministic TQG equations in (5.7) and (5.10) is given by

\frac{d F}{d t} = {F, H} = - \int {[\begin{matrix} δ F / δ q \\ δ F / δ b \end{matrix}]}^{T} [\begin{matrix} J (q - b, \cdot) & J (b, \cdot) \\ J (b, \cdot) & 0 \end{matrix}] [\begin{matrix} δ H_{TQG} / δ q = - ψ \\ δ H_{TQG} / δ b = - h_{1} / 2 \end{matrix}] d^{2} x

(5.17)

for the energy Hamiltonian in Eq. (5.15). The Poisson matrix in (5.17) is a Poisson deformation of the Lie–Poisson Hamiltonian matrix by the invertible linear transformation $(q, b) \to (q - b, b)$ ⁠. This deformation preserves the integral conservation laws $C_{Φ, Ψ}$ in Eq. (5.16) because their variational derivatives lie in the kernel of the Lie–Poisson Hamiltonian matrix operator. That is the Poisson bracket ${C_{Φ, Ψ}, H}$ vanishes for every choice of functional H(q, b), not just for the energy Hamiltonian $H_{TQG} (q, b)$ in (5.15), which generates the deterministic TQG equations in (5.7) and (5.10) under the action of the Poisson bracket in (5.17). The Hamiltonian structure allows one to investigate linear and nonlinear stability of the TQG model. This will be done in future work.

3. Linear stability analysis of TQG

Equilibrium solutions of the TQG equations (5.10) satisfy $J (ψ_{e}, q_{e}) = J (ψ_{e}, b_{e}) = J (h_{1}, b_{e}) = 0$ ⁠. For example, in a channel geometry with streamwise coordinate x and spanwise coordinate y, a TQG equilibrium which satisfies the boundary conditions (5.8) is given by setting

\begin{matrix} \nabla ψ_{e} = - U \hat{y}, \nabla q_{e} = (U - β) \hat{y}, \nabla f = - β \hat{y}, \nabla b_{e} = - B \hat{y}, \nabla h = - H \hat{y}, \end{matrix}

(5.18)

where the equilibrium parameters $U, β, B, H \in ℝ$ in (5.18) are all constants. By linearizing (5.10) around the steady state whose gradients of $ψ_{e}, q_{e}, b_{e}$ and h_e given in (5.18) all point in the y-direction, one finds the following evolution equations for the perturbations from equilibrium $q', ψ', b'$ for TQG given by

\begin{matrix} \partial_{t} q' + U q_{x}^{'} + (U + B - β) ψ_{x}^{'} = (U - H / 2) b_{x}^{'}, \\ \partial_{t} b' + U b_{x}^{'} + B ψ_{x}^{'} = 0, \end{matrix}

(5.19)

where we have used $\hat{z} \times (- \hat{y}) = \hat{x}$ ⁠. For later use, let us also note that for q defined in (5.10), we have $q' = Δ ψ' - ψ'$ ⁠. Upon inserting the 2D plane wave assumption

q' = e^{i (k \cdot x - ν t)} \hat{q}, ψ' = e^{i (k \cdot x - ν t)} \hat{ψ}, b' = e^{i (k \cdot x - ν t)} \hat{b},

(5.20)

with $\hat{q} = - (1 + | k |^{2}) \hat{ψ} = - (1 + k^{2} + l^{2}) \hat{ψ}$ ⁠, we find first from (5.19) that

\begin{matrix} (ν - k U) \hat{q} - k (U + B - β) \hat{ψ} = - k (U - H / 2) \hat{b}, \\ (ν - k U) \hat{b} = - k B \hat{ψ} . \end{matrix}

(5.21)

Consequently, upon eliminating $\hat{q}$ and $\hat{b}$ in favor of $\hat{ψ}$ in (5.21) we find

(ν - k U) (1 + | k |^{2}) \hat{ψ} + k (U + B - β) \hat{ψ} = k (U - H / 2) \hat{b} = - \frac{k^{2} (U - H / 2) B \hat{ψ}}{(ν - k U)} .

(5.22)

Hence, for $\hat{ψ} \neq 0$ we have a relation for the Doppler-shifted phase speed $C = (ν - k U) / k$ from

{(ν - k U)}^{2} (1 + | k |^{2}) + k (U + B - β) (ν - k U) + k^{2} (U - H / 2) B = 0,

(5.23)

via the quadratic formula, as $(1 + k^{2}) C^{2} + X C + Y = 0$ ⁠, with

C : = (ν - k U) / k = \frac{- X \pm \sqrt{X^{2} - 4 (1 + | k |^{2}) Y}}{2 (1 + | k |^{2})},

(5.24)

with X and Y defined as

X : = (U + B - β) and Y : = (U - H / 2) B,

(5.25)

for linearized TQG.

Thus, up to a shift in the definition of the constant parameter X in Eq. (5.24), the TQG system predicts the same qualitative wave properties of their corresponding linearized systems. For B = 0, the Doppler-shifted phase speed C in Eq. (5.24) yields the dispersion relation for Rossby waves. For $B \neq 0$ ⁠, we find a $| k |$ -dependent modification of the Doppler-shifted phase speed which becomes complex at high wavenumber $| k |$ ⁠. When the phase speed becomes complex, the TQG wave may decay on one branch and exponentiate on the other. It is worth investigating whether this branched structure of the phase speed may represent the conversion of one type of wave into another.

According to Eq. (5.24), the solutions of the TQG PDE are linearly well-posed. At high wavenumber $| k | ≫ 1$ ⁠, the Doppler-shifted phase speed becomes complex (so the solution is indeed linearly unstable). However, the growth-rate decays to zero as $1 / \sqrt{1 + | k |^{2}}$ ⁠, so the TQG PDE system is not ill-posed. That is, the solutions of TQG are well-posed in the sense of Hadamard, because this linear analysis has shown that both linearized systems have continuous dependence on initial conditions. The computational simulations of nonlinear TQG solutions show a transition to various kinds of instability as the simulated solutions develop high wavenumbers. The present linear analysis predicts a maximum growth rate at a finite wavenumber k_max. One would expect that the simulated solution at the length-scale corresponding to k_max should be very active. However, the wave activity falls off at wavenumbers greater than k_max as the Doppler-shifted phase speed C tends to zero, indicating loss of wave activity. One might expect that constant potential vorticity damping would control the wave activity at high wavenumber, so that viscosity might be unnecessary.

4. Hamiltonian formulation of stochastic TQG

As with the TRSW equations, introducing stochasticity via the Euler–Poincaré theorem will preserve the conservation laws in (5.16). However, introducing stochasticity would not preserve energy in (5.15), because the stochastic Lagrangian depends explicitly on time through the Brownian motion. Nonetheless, as discussed in remark 2.15 for the stochastic TRSW equations, the stochastic TQG equations may still possess a Hamiltonian formulation. By coupling the potential vorticity to noise, we can construct the noise Hamiltonians as follows:

H_{i} ○ d W_{t}^{i} = \int_{D} q ϑ_{i} (x) d x d y ○ d W_{t}^{i},

(5.26)

where $ϑ_{i} (x)$ are functions of space but not of time. Each $ϑ_{i} (x)$ is associated with an independent Brownian motion $W_{t}^{i}$ ⁠. For each i, we have $\nabla^{⊥} ϑ_{i} (x) = \hat{z} \times \nabla ϑ_{i} (x) = ξ_{i} (x)$ ⁠. Thus, by definition, the divergence of $ξ_{i} (x)$ vanishes for each i. By taking the sum of the Hamiltonian $H_{TQG} d t$ in (5.15) and the noise Hamiltonians $H_{i} ○ d W_{t}^{i}$ in (5.26), we obtain a semimartingale Hamiltonian. Inserting this augmented Hamiltonian into the Poisson bracket (5.17) yields the following stochastic TQG equations:

\begin{matrix} d q + J (ψ, q - b) d t + J (ϑ_{i}, q - b) ○ d W_{t}^{i} = - \frac{1}{2} J (h_{1}, b) d t, \\ d b + J (ψ, b) d t + J (ϑ_{i}, b) ○ d W_{t}^{i} = 0. \end{matrix}

(5.27)

Since we changed only the Hamiltonian to obtain the stochastic TQG equations, the conservation laws that correspond to (5.27) are immediate. We should remark that (5.27) are not the equations one would obtain by applying asymptotic analysis to the stochastic thermal L1 equations. Following the same methods as in the deterministic case leads to a stochastic set of equations with a much smaller family of integral quantities. Since adding SALT, besides the energy, preserves the conservation laws, proceeding via asymptotic analysis does not produce the SALT stochastic version of TQG.

5. Conservation laws for stochastic TQG

The stochastic TQG equations (5.27) do not conserve energy. This is because the Hamiltonian that generates the dynamics depends on time explicitly, due to the presence of the Brownian motions. However, Eqs. (5.27) do conserve the same set of integral quantities $C_{Φ, Ψ}$ ⁠, defined in (5.16), since the Poisson bracket for the deterministic TQG equations and the stochastic TQG equations is the same.

VI. CONCLUSION AND OUTLOOK

Our motivation in this paper has been to prepare the mathematical framework for our ongoing investigations of Stochastic Transport in Upper Ocean Dynamics (STUOD) by using the stochastic data assimilation algorithms developed and applied previously to determine the eigenvectors $ξ_{i} (x)$ in the cases of the stochastic Euler fluid equation and the 2-layer stochastic QG model in Cotter et al. (2018; 2019). This framework has been established by deriving a sequence of realistic 2D models of Upper Ocean Dynamics with buoyancy effects by using nested asymptotic expansions with a shared stochastic variational structure. The process of developing these sequential derivations has also revealed several open mathematical problems at each level of approximation for these new nonlinear stochastic partial differential equations, as listed below.

An extensive computational simulation study will be needed for classifying the solution behavior of these new stochastic TRSW and TQG equations. This computational study has been left as a future step, after having established its efficacy in Sec. V B.
These computational simulations will be required in the calibration of the eigenvectors $\nabla^{⊥} ϑ_{i} (x) = ξ_{i} (x)$ for the noise in Eq. (5.27) and their subsequent use in the new framework for data calibration, uncertainty quantification and data assimilation using particle filters following the SALT algorithm developed in Cotter et al. (2018; 2019). This future simulation study will prepare these models for applications in the analysis of the observed detailed upper ocean dynamics as seen in Fig. 1.
Investigation of the numerous potential effects of the horizontal buoyancy gradients appearing in the elliptic equation for the thermal L1 Lagrange multiplier u (3.46) has been left as an open mathematical problem for further analysis and computational simulation.
The issue of well-posedness of these new nonlinear stochastic partial differential equations (2.34) for TRSW and (5.27) for TQG has also been left for future mathematical investigation.
Finally, we recall that our derivation of the stochastic barotropic TQG balanced model in (5.27) has neglected the potentially important effects of baroclinic instabilities which tend to re-stratify the fluid. In particular, a future study with baroclinic TQG would extend the QG analysis of baroclinic instability of the currents around the Lofoten Basin given in Isachsen (2015) to include thermal effects. For an in-depth discussion of baroclinic effects in comparison with balanced models, see Callies et al. (2016).

ACKNOWLEDGMENTS

We are grateful for constructive suggestions offered in discussions with C. J. Cotter, B. Chapron, D. Crisan, S. R. Ephrati, B. Fox-Kemper, R. Hu, O. Lang, J. M. Leahy, J. C. McWilliams, E. Mémin, O. Street, and S. Takao. We are also grateful to the European Space Agency for access to the Sentinel 3B satellite data shown in Fig. 1. During this work, D.D.H. and W.P. were partially supported by ERC Synergy Grant No. 856408—STUOD (Stochastic Transport in Upper Ocean Dynamics). E.L. was supported by EPSRC Grant (No. EP/L016613/1). E.L. is also grateful for the warm hospitality shown to him at the Imperial College London EPSRC Centre for Doctoral Training in the Mathematics of Planet Earth mpecdt.org.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1.

Abraham

,

R.

, and

Marsden

,

J. E.

,

Foundations of Mechanics

(

Benjamin/Cummings Publishing Company

,

Reading, Massachusetts

,

1978

), Vol.

36

.

Google Scholar

2.

Allen

,

J. S.

, and

Holm

,

D. D.

, “

Extended-geostrophic Hamiltonian models for rotating shallow water motion

,”

Phys. D

98

(

2–4

),

229

–

248

(

1996

).

https://doi.org/10.1016/0167-2789(96)00120-0

Google Scholar

Crossref

3.

Bashmachnikov

,

I.

,

Belonenko

,

T.

,

Kuibin

,

P.

,

Volkov

,

D.

, and

Foux

,

V.

, “

Pattern of vertical velocity in the Lofoten vortex (the Norwegian Sea)

,”

Ocean Dyn.

68

(

12

),

1711

–

1725

(

2018

).

https://doi.org/10.1007/s10236-018-1213-1

Google Scholar

Crossref

4.

Bashmachnikov

,

I. L.

,

Kozlov

,

I. E.

,

Petrenko

,

L. A.

,

Glock

,

N. I.

, and

Wekerle

,

C.

, “

Eddies in the North Greenland Sea and Fram Strait from satellite altimetry, SAR and high-resolution model data

,”

J. Geophys. Res.: Oceans

125

,

e2019JC015832

,

https://doi.org/10.1029/2019JC015832

(

2020

).

Google Scholar

Crossref

5.

Bashmachnikov

,

I. L.

,

Sokolovskiy

,

M. A.

,

Belonenko

,

T. V.

,

Volkov

,

D. L.

,

Isachsen

,

P. E.

, and

Carton

,

X.

, “

On the vertical structure and stability of the Lofoten vortex in the Norwegian Sea

,”

Deep Sea Res., Part I

128

,

1

–

27

(

2017

).

https://doi.org/10.1016/j.dsr.2017.08.001

Google Scholar

Crossref

6.

Bernsen

,

E.

,

Bokhove

,

O.

, and

van der Vegt

,

J. J. W.

, “

A (dis)continuous finite element model for generalized 2D vorticity dynamics

,”

J. Comput. Phys.

211

(

2

),

719

–

747

(

2006

).

https://doi.org/10.1016/j.jcp.2005.06.008

Google Scholar

Crossref

7.

Beron-Vera

,

F. J.

, “

Multilayer shallow-water model with stratification and shear

,”

Rev. Mex. Fis.

(to be published); arXiv:physics/0312083 (

2020

).

Google Scholar

8.

Bröcker

,

J.

,

Calderhead

,

B.

,

Cheraghi

,

D.

,

Cotter

,

C.

,

Crisan

,

D.

,

Holm

,

D.

,

Kuna

,

T.

,

Pelloni

,

B.

,

Shepherd

,

T.

, and

Weller

,

H.

,

Mathematics of Planet Earth: A Primer

(

World Scientific

,

2018

).

Google Scholar

9.

Callies

,

J.

,

Flierl

,

G.

,

Ferrari

,

R.

, and

Fox-Kemper

,

B.

, “

The role of mixed-layer instabilities in submesoscale turbulence

,”

J. Fluid Mech.

788

,

5

–

41

(

2016

).

https://doi.org/10.1017/jfm.2015.700

Google Scholar

Crossref

10.

Cotter

,

C.

,

Dan

,

C.

,

Holm

,

D. D.

,

Wei

,

P.

, and

Shevchenko

,

I.

, “

Modelling uncertainty using stochastic transport noise in a 2-layer quasi-geostrophic model

,”

Found. Data Sci.

2

(

2

),

173

(

2020

).

https://doi.org/10.3934/fods.2020010

Google Scholar

Crossref

11.

Cotter

,

C.

,

Crisan

,

D.

,

Holm

,

D. D.

,

Pan

,

W.

, and

Shevchenko

,

I.

, “

Numerically modeling stochastic Lie transport in fluid dynamics

,”

Multiscale Model. Simul.

17

(

1

),

192

–

232

(

2019

).

https://doi.org/10.1137/18M1167929

Google Scholar

Crossref

12.

Cotter

,

C. J.

,

Holm

,

D. D.

, and

Percival

,

J. R.

, “

The square root depth wave equations

,”

Proc. R. Soc. A

466

(

2124

),

3621

–

3633

(

2010

).

https://doi.org/10.1098/rspa.2010.0124

Google Scholar

Crossref

13.

Crisan

,

D.

,

Holm

,

D. D.

,

Leahy

,

J.-M.

, and

Nilssen

,

T.

, “

Variational principles for fluid dynamics on rough paths

,” arXiv:2004.07829 (

2020

).

Google Scholar

14.

de Leon

,

A. B.

,

Darryl Holm

,

D. D.

,

Erwin

,

L.

, and

So

,

T.

, “

Implications of Kunita–Itô–Wentzell formula for k-forms in stochastic fluid dynamics

,”

J. Nonlinear Sci.

30

,

1421

(

2020

).

https://doi.org/10.1007/s00332-020-09613-0

Google Scholar

Crossref

15.

Eliassen

, A.

,

The Quasi-Static Equations of Motion with Pressure as Independent Variable

(

Grøndahl & sons boktr., I kommisjon hos Cammermeyers boghandel

,

1949

).

Google Scholar

16.

Fedorov

,

A. M.

, and

Belonenko

,

T. V.

, “

Interaction of mesoscale vortices in the Lofoten basin based on the GLORYS database

,”

Russ. J. Earth Sci.

20

,

1

(

2020

).

https://doi.org/10.2205/2020ES000694

Google Scholar

Crossref

17.

Filyushkin

,

B. N.

,

Sokolovskiy

,

M. A.

, and

Lebedev

,

K. V.

, “

Evolution of an intrathermocline lens over the Lofoten basin

,” in

The Ocean in Motion

(

Springer

,

2018

), pp.

333

–

347

.

Google Scholar

Crossref

18.

Gottlieb

,

S.

, “

On high order strong stability preserving Runge–Kutta and multistep time discretizations

,”

J. Sci. Comput.

25

(

1

),

105

–

128

(

2005

).

https://doi.org/10.1007/s10915-004-4635-5

Google Scholar

Crossref

19.

Holm

,

D. D.

, “

Hamiltonian formulation of the baroclinic quasigeostrophic fluid equations

,”

Phys. Fluids

29

(

1

),

7

–

8

(

1986

).

https://doi.org/10.1063/1.865956

Google Scholar

Crossref

20.

Holm

,

D. D.

, “

Variational principles for stochastic fluid dynamics

,”

Proc. R. Soc. A

471

(

2176

),

20140963

(

2015

).

https://doi.org/10.1098/rspa.2014.0963

Google Scholar

Crossref

21.

Holm

,

D. D.

, and

Luesink

,

E.

, “

Stochastic wave-current interaction in thermal shallow water dynamics

,”

J. Nonlinear Sci.

31

,

29

(

2021

).

https://doi.org/10.1007/s00332-021-09682-9

Google Scholar

Crossref

22.

Holm

,

D. D.

, and

Zeitlin

,

V.

, “

Hamilton's principle for quasigeostrophic motion

,”

Phys. Fluids

10

(

4

),

800

–

806

(

1998

).

https://doi.org/10.1063/1.869623

Google Scholar

Crossref

23.

Holm

,

D. D.

,

Marsden

,

J. E.

,

Ratiu

,

T.

, and

Weinstein

,

A.

, “

Nonlinear stability of fluid and plasma equilibria

,”

Phys. Rep.

123

(

1–2

),

1

–

116

(

1985

).

https://doi.org/10.1016/0370-1573(85)90028-6

Google Scholar

Crossref

24.

Holm

,

D. D.

,

Marsden

,

J. E.

, and

Ratiu

,

T. S.

, “

The Euler–Poincaré equations and semidirect products with applications to continuum theories

,”

Adv. Math.

137

(

1

),

1

–

81

(

1998

).

https://doi.org/10.1006/aima.1998.1721

Google Scholar

Crossref

25.

Holm

,

D. D.

,

Schmah

T.

, and

Stoica

,

C.

,

Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions, Volume 12

(

Oxford University Press

,

2009

).

Google Scholar

26.

Isachsen

,

P. E.

, “

Baroclinic instability and the mesoscale eddy field around the Lofoten basin

,”

J. Geophys. Res.: Oceans

120

(

4

),

2884

–

2903

,

https://doi.org/10.1002/2014JC010448

(

2015

).

Google Scholar

Crossref

27.

Ito

,

K.

,

An Introduction to Probability Theory

(

Cambridge University Press

,

1984

).

Google Scholar

Crossref

28.

Kraichnan

,

R. H.

, “

Inertial ranges in two-dimensional turbulence

,”

Phys. Fluids

10

(

7

),

1417

–

1423

(

1967

).

https://doi.org/10.1063/1.1762301

Google Scholar

Crossref

29.

Maximenko

,

N.

,

Niiler

,

P.

,

Centurioni

,

L.

,

Rio

,

M.-H.

,

Melnichenko

,

O.

,

Chambers

,

D.

,

Zlotnicki

,

V.

, and

Galperin

,

B.

, “

Mean dynamic topography of the ocean derived from satellite and drifting buoy data using three different techniques

,”

J. Atmos. Oceanic Technol.

26

(

9

),

1910

–

1919

(

2009

).

https://doi.org/10.1175/2009JTECHO672.1

Google Scholar

Crossref

30.

McWilliams

,

J. C.

, “

Irreducible imprecision in atmospheric and oceanic simulations

,”

Proc. Natl. Acad. Sci.

104

(

21

),

8709

–

8713

(

2007

).

https://doi.org/10.1073/pnas.0702971104

Google Scholar

Crossref

31.

McWilliams

,

J. C.

, “

A survey of submesoscale currents

,”

Geosci. Lett.

6

(

1

),

3

(

2019

).

https://doi.org/10.1186/s40562-019-0133-3

Google Scholar

Crossref

32.

Rathgeber

,

F.

,

Ham

,

D. A.

,

Mitchell

,

L.

,

Lange

,

M.

,

Luporini

,

F.

,

Mcrae

,

A. T. T.

,

Bercea

,

G.-T.

,

Markall

,

G. R.

, and

Kelly

,

P. H. J.

, “

Firedrake: Automating the finite element method by composing abstractions

,”

ACM Trans. Math. Software

43

(

3

),

1

–

27

(

2017

).

https://doi.org/10.1145/2998441

Google Scholar

Crossref

33.

Ripa

,

P.

, “

Conservation laws for primitive equations models with inhomogeneous layers

,”

Geophys. Astrophys. Fluid Dyn.

70

(

1–4

),

85

–

111

(

1993

).

https://doi.org/10.1080/03091929308203588

Google Scholar

Crossref

34.

Ripa

,

P.

, “

Low frequency approximation of a vertically averaged ocean model with thermodynamics

,”

Rev. Mex. Fís.

42

(

1

),

117

–

135

(

1995

), see https://rmf.smf.mx/ojs/rmf/article/view/2559/2527.

Google Scholar

35.

Ripa

,

P.

, “

On the validity of layered models of ocean dynamics and thermodynamics with reduced vertical resolution

,”

Dyn. Atmos. Oceans

29

(

1

),

1

–

40

(

1999

).

https://doi.org/10.1016/S0377-0265(98)00056-6

Google Scholar

Crossref

36.

Salmon

,

R.

, “

Practical use of Hamilton's principle

,”

J. Fluid Mech.

132

,

431

–

444

(

1983

).

https://doi.org/10.1017/S0022112083001706

Google Scholar

Crossref

37.

Street

,

O. D.

, and

Crisan

,

D.

, “

Semi-martingale driven variational principles

,”

Proc. R. Soc. A

477

,

20200957

(

2021

).

https://doi.org/10.1098/rspa.2020.0957

Google Scholar

Crossref

38.

Volkov

,

D. L.

,

Kubryakov

,

A. A.

, and

Lumpkin

,

R.

, “

Formation and variability of the Lofoten basin vortex in a high-resolution ocean model

,”

Deep Sea Res., Part I

105

,

142

–

157

(

2015

).

https://doi.org/10.1016/j.dsr.2015.09.001

Google Scholar

Crossref

39.

Warneford

,

E. S.

, and

Dellar

,

P. J.

, “

The quasi-geostrophic theory of the thermal shallow water equations

,”

J. Fluid Mech.

723

,

374

–

403

(

2013

).

https://doi.org/10.1017/jfm.2013.101

Google Scholar

Crossref

40.

Zeitlin

,

V.

,

Geophysical Fluid Dynamics: Understanding (Almost) Everything with Rotating Shallow Water Models

(

Oxford University Press

,

2018

).

Google Scholar

Crossref

2021

Author(s)

Stochastic mesoscale circulation dynamics in the thermal ocean

I. INTRODUCTION

A. The TRSW equations

II. THE THERMAL ROTATING SHALLOW WATER (TRSW) MODEL

A. Deterministic TRSW equations

1. Variational formulation

B. The Euler–Poincaré theorem

1. Variational derivatives of functionals

2. Euler–Poincaré theorem

C. Kelvin–Noether circulation theorem

1. Conservation laws for deterministic TRSW

2. Noether's theorem

D. TRSW with stochastic advection by Lie transport (SALT)

III. BALANCED INTERPRETATIONS OF TRSW

A. Thermal geostrophic balance

B. Moving into the balanced frame

1. Stationary thermal geostrophic balance as “mean dynamic topography (MDT)”

C. The Eliassen approximation

D. The Allen–Holm approach

1. Conservation laws for deterministic TL1

2. Determining the Lagrange multiplier

3. No horizontal buoyancy gradients

4. General case

IV. THE ELIASSEN APPROXIMATION OF STOCHASTIC TRSW

A. Stochastic thermal geostrophic balance

B. Stochastic TL1

1. Determining the Lagrange multiplier

V. THERMAL QG MODEL

A. Deterministic TQG model

1. Asymptotic expansion to the deterministic TQG regime

2. Elliptic equation for TL1

B. Numerical TQG example

C. Hamiltonian formulation of TQG

1. Conservation laws for deterministic TQG

2. Hamiltonian formulation of deterministic TQG

3. Linear stability analysis of TQG

4. Hamiltonian formulation of stochastic TQG

5. Conservation laws for stochastic TQG

VI. CONCLUSION AND OUTLOOK

ACKNOWLEDGMENTS

DATA AVAILABILITY

References

Contents

Figures & Tables

Supplementary Material

References

Related

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only