1 Introduction

Weather forecasting, climate change prediction and global ocean circulation all face the same fundamental challenge to create models which incorporate the effects of measurement error and uncertainty due to unresolved scales, unknown physical phenomena and incompleteness of observed data. We tackle this issue of modelling effects of unknown causes in observational science by applying new methods in stochastic data-driven modelling which are designed to predict both future measurements and their uncertainty, based on analysing the available data for the problem at hand.

For example, a common approach for modelling and simulating climate and weather is based on stochastic parametrisation. For recent reviews of stochastic parametrisation in geophysical fluid dynamics (GFD), see e.g. Berner et al. (2012, 2017), Gottwald et al. (2016). The fundamental conclusions of Berner et al. (2012) are twofold:

A posteriori addition of stochasticity to an already tuned model is simply not viable.

Stochasticity must be incorporated at a very basic level within the design of physical process parametrisations and improvements to the dynamical core.

A new approach (Holm 2015) which meets the challenge of incorporating stochastic parametrisation at the fundamental level enunciated in Berner et al. (2012) introduces stochastic transport directly into the loop velocity in Kelvin’s circulation theorem. The dynamical quantities of physical interest are then modelled together with their statistical uncertainty, and data assimilation is used to reduce that uncertainty. This is the SALT approach.

The SALT (stochastic advection by Lie transport) approach combines stochasticity in the velocity of the fluid material loop in Kelvin’s circulation theorem with ensemble forecasting. The ensemble forecasting in SALT has been coordinated with the results of the particle filtering method of data assimilation. A protocol for applying the SALT approach in combination with data assimilation based on comparing fine-scale and coarse-scale computational simulations has recently been established in Cotter et al. (2018, 2019a). These results demonstrate the capability of the SALT approach to successfully reduce forecast uncertainty in a variety of test problems for fluid dynamics in two spatial dimensions. The three-dimensional SALT theory has been developed and analysed to determine their existence, uniqueness and blow-up criterion in Crisan et al. (2019), but it still awaits computational implementation at the present time.

The present paper aims to use the SALT approach for fluid dynamics described above to provide a barotropic (vertically averaged) description of wave–current interaction (WCI) in a stratified incompressible fluid flow, by incorporating stochastic fluid transport and circulation with nonlinear dispersive wave propagation internally and on the free surface. In doing so, this paper combines a variational principle approach with asymptotic analysis to derive simplified models. Historically in ocean modelling, the rapid propagation of the barotropic (or external) mode representing disturbances on the free surface, for example, has required special handling; because otherwise incorporating the simulation of its rapid time scale and multicomponent physical processes would tend to occupy an inordinate amount of computer power (Dukowicz and Smith 1994; Fox-Kemper et al. 2019).

Fig. 1
figure 1

The flow diagram of approximations via vertical averages and asymptotic expansions in Camassa and Holm (1992)

In addressing this challenge, the Camassa–Holm 1992 model (referred to as CH92 hereafter) derived in Camassa and Holm (1992) used vertical averaging to transform the 3D Euler–Boussinesq fluid equations into a family of 2D stratified ‘rotating shallow water equations’ which incorporate effects of weak deviations from hydrostatic balance, weak stratification and strong topography. Through a series of approximations and asymptotic limits, the CH92 model was found to contain the Kadomtsev–Petviashvili (KP) and Korteweg–de Vries (KdV) equations in a rotating frame, as shown in Fig. 1.

The present paper will develop two families of stochastic models of barotropic wave–current interaction for mesoscale and submesoscale ocean dynamics based on the deterministic CH92 model and its further development in Camassa et al. (1996, 1997). Our approach combines dimensional analysis, asymptotic expansions and vertical averaging to obtain the barotropic component of the fluid motion, as done in Camassa and Holm (1992), extended first to the Euler–Poincaré variational approach of Holm et al. (1998) and then to the SALT approach (Holm 2015) for introducing stochasticity. In the Euler–Poincaré version of the SALT approach, the approximation of the Lagrangian is separate from the introduction of stochasticity. The asymptotic expansions are applied to the Lagrangian first, and then, the stochasticity is introduced in the variations of the Eulerian variables, which depend on spatially smooth maps with stochastic time dependence. In the variational step to include stochasticity, one also introduces the Strouhal number. In the process, we handle the barotropic effects by vertically averaging, applied either to the equations of motion as in Wu (1981), or to the variational principle for SALT (Holm 2015). Of course, the vertical averaging procedure eliminates vertical buoyancy gradients. However, horizontal gradients of the vertically averaged buoyancy remain. Here, the equations obtained after vertical averaging which retain horizontal gradients of buoyancy will be called thermal equations. This name applies because the buoyancy plays the role of entropy per unit mass in the equation of state for adiabatic compressible fluid flows. Likewise, the variation of the energy with respect to the buoyancy plays the role of temperature in the adiabatic compressible fluid case. Thus, the present paper aims to incorporate stochasticity into the theory of nonlinear dispersive water waves interacting with horizontal buoyancy gradients, as governed by vertically averaged fluid equations. This stochastic theory of wave–current interaction in thermal shallow water dynamics is expected to be useful for quantifying uncertainty and perhaps even reducing it by using data assimilation in the SALT approach (Holm et al. 2020).

Background. A framework for combining data with existing models in a probabilistic manner was presented in Holm (2015), where a stochastic variational principle for continuum mechanics was introduced. This stochastic variational principle enables one to derive stochastic models of inviscid fluid dynamics which satisfy a Kelvin circulation theorem, starting from the Lagrangian of the corresponding deterministic fluid model and using a Clebsch constraint to introduce the stochastic advection by lie transport (SALT). This approach decomposes the fluid velocity vector field into the sum of a drift velocity and a Stratonovich stochastic velocity. The former is obtained from the constrained variational principle, and the latter is determined by analysing available data according to the protocol established in Cotter et al. (2019a, 2018). The constraints may be introduced either by imposing the advection equations for the relevant physical quantities of the model, or equivalently by imposing the advection equation for the fluid labels.

Recently, in Bethencourt de Leon et al. (2020), the known Euler–Poincaré and Hamilton–Pontryagin stochastic variational principles were reformulated and shown to be equivalent to the Clebsch variant, by proving existence and uniqueness of the solution of the SALT advection constraint. The noise used in the Holm (2015) approach also appears in Cotter et al. (2017), where the decomposition for SALT of the fluid velocity vector field into the sum of a drift velocity and a Stratonovich stochastic velocity was derived by using multi-time homogenisation theory. Many subsequent investigations of the properties of the equations of fluid dynamics with the SALT modification have appeared in the literature over the last four years. In particular, the SALT approach preserves most physical conservation laws by construction, while it also possesses much of the analytical structure of the underlying deterministic model. For example, in Crisan et al. (2019), the three-dimensional SALT Euler equations are shown to have the same local-in-time existence and uniqueness analytical properties as the deterministic version, as well as the same Beale–Kato–Majda (Beale et al. 1984) criterion for blow-up of solutions. In Geurts et al. (2017), the Lorenz 63 equations are derived from Rayleigh–Bénard convection with this type of stochasticity and the rate of convergence towards the attractor is shown to be preserved by this type of noise. From a more operational point of view, in Cotter et al. (2019a), SALT was introduced into the two-dimensional Euler equations and it was shown that the stochastic equations, which are solved on a coarse grid, mimic the deterministic equations, which are solved on a fine grid, for a significant period of time. In Cotter et al. (2018), a similar result was established for the flow in a channel of a two layer quasi-geostrophic system.

In this paper, we are concerned with consistency of SALT under asymptotic expansion and analysis for the simulation of the barotropic mode in ocean dynamics. As mentioned earlier, the barotropic mode in ocean dynamics is the fastest excitation in the free-surface dynamics. It is treated separately (for example, by subcycling) in most 3D simulations of large scale ocean circulation. The issue of the free-surface treatment which motivated the original investigation of the various types of nonlinear wave behaviour in Camassa and Holm (1992) is still of current concern.

A motivating question for introducing SALT into nonlinear dispersive water wave theory to be addressed in the present paper is: How can one use available data to quantify the uncertainty due to the barotropic mode in the free-surface treatment for computational simulations? This work is done in preparation for using the data assimilation methods of Cotter et al. (2018, 2019a) to reduce that uncertainty e.g. by using satellite data.

As in Camassa and Holm (1992), we will combine asymptotic analysis with the vertical averaging principle of Wu (1981) to derive a sequence of two-dimensional barotropic models. This averaging principle will be applied both on the equations and also on the variational principle. The latter turns out to be advantageous in situations where the Strouhal number (the ratio of the chosen time scale over the natural time scale induced by the length and fluid velocity scales) is not equal to unity. The starting point of these derivations is the three-dimensional rotating stratified Euler model, a three-dimensional fluid model that includes the effects of rotation and buoyancy stratification. By making assumptions about the buoyancy stratification, we transition into the Euler–Boussinesq model. Here, we apply the averaging principle to derive two-dimensional models with nonhydrostatic effects, rotation and stratification. The two-dimensional models will be derived with respect to two different time scales: the first time scale is the natural one, and the second is the time scale that corresponds to gravity waves. When the time scale is the natural one, the Strouhal number is equal to unity, which means that the asymptotic analysis applied to the equations and the asymptotic analysis applied inside the variational principle lead to the same result at each order in the asymptotic expansion. The assumption that the free surface amplitude is very small leads to the Great Lake, Lake and Benney long wave equations, first derived in Camassa et al. (1996, 1997); Benney (1973), respectively, although in this paper we also include the effects of rotation, stratification and stochasticity. The second scaling regime is where the Strouhal number is equal to the inverse of the Froude number. This scaling regime leads to equations in the Green–Naghdi (Green and Naghdi 1976) class, if the free surface amplitude is assumed to be small, rather than very small. This derivation was first accomplished in Camassa and Holm (1992), where also a Kadomtsev–Petviashvili equation is derived, augmented by the effects of rotation and bathymetry. In the presence of stochasticity, however, this derivation cannot be done directly. As we shall see, in the situation where the Strouhal number is not equal to unity, asymptotic analysis applied to the equations fails to respect the geometric structure of the problem, but the asymptotic analysis of the variational principle does preserve the geometric structure.

In Holm (2015), the SALT vector field whose characteristic curves generate the stochastic Lagrangian fluid trajectories is defined as

$$\begin{aligned} {\mathsf {d}}{\varvec{\chi }}_t := {\mathbf {u}}({\mathbf {x}},t) dt + \sum _{i=1}^M {\varvec{\xi }}_i({\mathbf {x}})\circ dW_t^i. \end{aligned}$$

Here, \({\mathbf {u}}({\mathbf {x}},t)\) is the fluid velocity field, \({\varvec{\xi }}_i({\mathbf {x}})\) are the vector fields that represent spatial velocity–velocity correlations, \(W_t^i\) denotes independent, identically distributed Wiener processes for each \(i=1,\dots ,M\), and the symbol \(\circ \) means Stratonovich integration. The number M of eigenvectors \({\varvec{\xi }}_i({\mathbf {x}})\) required for a given level of accuracy can be determined via the amount of variance required from a principal component analysis, or via empirical orthogonal function analysis. Via data assimilation procedures, in particular via novel high-dimensional particle filtering methods, the uncertainty may be controlled and reduced dramatically when even a small amount of new data is observed, as shown in Cotter et al. (2018, 2019a). As we shall see, the variational approach of SALT used here has the additional advantage of preserving the Kelvin circulation theorem and the Hamiltonian framework, both of which have been fundamental in the history of studying wave–current interaction and now can be made stochastic.

1.1 Overview of the Paper

The starting point, described in Sect. 1, will be the introduction of a number of tools which are invaluable for this work. First we will introduce the stochastic Euler–Poincaré variational principle, Kelvin’s circulation theorem and an averaging principle. Then, starting with the rotating, stratified Euler equations, we will assume that the buoyancy stratification is weak enough to allow us to work with the Euler–Boussinesq equations. This is a justified assumption when the goal is to model the ocean. The flow in wave–current interaction is primarily incompressible, so the models used here will reflect this property. The ocean is shallow, compared to the horizontal distances of interest. In particular, the characteristic height scale is much smaller than the characteristic horizontal scales. This situation allows a reduction in spatial dimension by vertically integrating the Euler–Boussinesq equations to find the vertical average of the nonlinearity and an unknown vertically averaged pressure. Not surprisingly, these are the two terms which we cannot determine from the averaged equations alone. In order to derive a set of closed equations, we will turn to asymptotic analysis, which we will execute in two different regimes. Within each of those two regimes, we will apply asymptotic analysis in two different ways. In the first regime, called “long time-very small wave scaling”, the time scale is determined by the ratio of the characteristic velocity scale and horizontal length scale, and with very small wave amplitude. The second regime, called the “short time-small wave scaling”, will employ the time scale based on the gravity wave speed and a characteristic horizontal length. The vertical averaging principle of Wu (1981) will be applied, both on the 3D equations and on the corresponding Euler–Poincaré Lagrangian. We will show that in the first regime, the approaches coincide and produce the same equations. In the second regime the asymptotic analysis requires special treatment, as the Strouhal number is not equal to unity in that situation. This difference in Strouhal number means that the material derivative contributes at two different orders in the asymptotic expansion. We shall focus on deriving two-dimensional stochastic fluid models in these two different time-scale regimes, starting from a model for a three-dimensional stochastic fluid with rotation and stratification in a shallow box with bathymetry and a free surface.

In Sect. 2, “the long time-very small wave scaling regime” of the Euler–Boussinesq equations with negligible buoyancy stratification will be derived from asymptotics applied to the equations and to the corresponding Lagrangian. At leading order, this will give rise to the Benney long wave equations, before making the columnar motion assumption. It will produce the stochastic and rotating version of the Lake equations after assuming that the motion is columnar. The Benney equations have an interesting mathematical structure, such as an infinity of conservation laws, as presented in Kupershmidt (2006). From a different perspective, the rotating Lake equations are also obtained after assuming that the rotating shallow water equations have a rigid lid. At the next order, we find the stochastic and rotating version of the Great Lake equations (Camassa et al. 1996, 1997), which can be interpreted as the rigid lid version of the Green–Naghdi equations. The deterministic versions of the Lake and Great Lake equations are both globally well-posed in time, as shown in Levermore et al. (1996a, 1996b).

In Sect. 3, “the short time-small wave scaling” of the Euler–Boussinesq equations with non-negligible buoyancy will be considered. The results of the asymptotics obtained in this regime are quite different from those obtained in the previous section. In this regime, the Strouhal number is not unity and the asymptotics on the equations provides us with a set of equations which does not satisfy the Kelvin circulation theorem. The reader may refer to Camassa and Holm (1992) for the deterministic derivation of these equations and their relation to the Kadomtsev–Petviashvili equation. The corresponding asymptotic analysis on the Lagrangian does give a set of equations that satisfy Kelvin’s circulation theorem, though, as it results in a buoyant version of the Green–Naghdi equations. As it turns out, a variational derivation of equations for the free surface alone is not available. Hence, the corresponding Boussinesq-type water wave equations are not available, unless model assumptions in the variational principle were to be changed. It will be shown that a hierarchy of stochastic Camassa–Holm equations can be derived from this point of view leading to the stochastic Korteweg–De Vries equation, as well.

In Sect. 4, we summarise by diagramming the pathways which relate the sequences of approximations leading to the results obtained in this paper for each of the two families of nonlinear fluid wave equations.

2 Stochastic Variational Principle and Averaging Principle

Central to this work is the stochastic Euler–Poincaré variational principle, presented in Bethencourt de Leon et al. (2020), which is equivalent to the variational principle in Holm (2015). However, the Euler–Poincaré variational principle uses prescribed variations, rather than variations induced by constraints used in Holm (2015). The most general version of the Euler–Poincaré theorem is formulated on the Lie algebra of a semidirect product Lie group and uses the language of differential geometry and representation theory, which first appeared deterministically in Holm et al. (1998). For fluids, the group of interest is the diffeomorphism group, which is the space of differentiable maps whose inverse maps are equally differentiable. The group action is composition of functions. The group of diffeomorphisms is a suitable group for geometric mechanics in the sense of Ebin and Marsden (1970). In order to state the Euler–Poincaré theorem, we first need to introduce some notation.

Notation. The domain of interest for the paper is a three-dimensional box with bathymetry specified by h(xy) and a free surface \(\zeta (x,y,t)\), as illustrated in Fig. 2. The domain, which we will call \(\Omega \), is a subset of \(\mathbb {R}^3\) and equipped with Cartesian coordinates. As we proceed, we will present the Euler–Poincaré theorem and its sequence of implications in \(\mathbb {R}^3\) vector calculus, rather than using the more abstract differential geometric notation.

Two-dimensional and three-dimensional objects will be distinguished by putting a subscript on the three-dimensional objects, as follows:

$$\begin{aligned} {\mathbf {x}}_3 = ({\mathbf {x}},z), \qquad {\mathbf {u}}_3 = ({\mathbf {u}},w), \qquad \nabla _3 = \left( \nabla , \frac{\partial }{\partial z}\right) . \end{aligned}$$
(1.1)

Here, \({\mathbf {x}}_3\) denotes the coordinate system, \({\mathbf {u}}_3\) is the fluid vector field, and \(\nabla _3\) is the gradient. In oceanographic terms, we will work primarily with mesoscale dynamics, where the typical horizontal length scale L is in the order of one hundred kilometres, or more, and the typical depth H is four kilometres. Hence, the domain will be shallow. The rotation of the planet is included by introducing the vector potential \({\mathbf {R}}_3({\mathbf {x}}_3)=({\mathbf {R}}({\mathbf {x}}),0)\) for the Coriolis parameter, \(f({\mathbf {x}}){\hat{\mathbf {z}}}\), so that

$$\begin{aligned} \nabla _3\times {\mathbf {R}}_3({\mathbf {x}}) = f({\mathbf {x}}){\hat{\mathbf {z}}} \end{aligned}$$
(1.2)

and \({\hat{\mathbf {z}}}\) is the unit vector in the vertical direction.

Fig. 2
figure 2

The 3D flow domain, \(\Omega \). The wavy green surface is the free surface \(\zeta (x,y,t)\), and the wavy blue surface is the bathymetry h(xy). This figure is not to scale, as the horizontal length scale of an ocean domain is typically much larger than its height scale. In the paper, we will assume that \(L_x = L_y = L\)

By \(\mathfrak {X}(\Omega )\), we denote the space of vector fields over \(\Omega \) and by \(V^*\) we mean the abstract vector space of advected quantities, which are usually tensor fields of different degrees. In this paper, the elements of \(V^*\) that we will consider are buoyancy b, which is a scalar function, the density D, which is a volume form and later in the two-dimensional setting, we will consider the depth \(\eta ({\mathbf {x}},t):= \zeta ({\mathbf {x}},t) + h({\mathbf {x}})\), which is the volume form in that scenario. The stochastic vector fields which generate the Lagrangian transport is given below. In this paper, we will always work with the same stochastic basis, consisting of a set, a \(\mathbb {P}\)-complete \(\sigma \)-algebra, a probability measure \(\mathbb {P}\) and a right-continuous filtration. All stochastic processes that we will consider are adapted to the filtration. The stochastic vector field for three-dimensional transport is given by

$$\begin{aligned} {\mathsf {d}}{\varvec{\chi }}_{3t} := {\mathbf {u}}_3({\mathbf {x}}_3,t)dt + \sum _{i=1}^M {\varvec{\xi }}_{3i}({\mathbf {x}}_3)\circ dW_t^i\,. \end{aligned}$$
(1.3)

The stochastic vector field \({\varvec{\chi }}_{3t}\) in (1.3) is an example of a semimartingale.

Definition 1.1

(Semimartingale). A cádlág process Y is a semimartingale if it can be written as

$$ Y_t = Y_0 + M_t + V_t \,,$$

where M is a cádlág local martingale, V is a cádlág process of finite variation and \(M_0 = V_0 = 0\).

The adjective cádlág stands for “right continuous with a limit on the left”. The processes we will be considering will have continuous paths almost surely. Continuous semimartingales can be decomposed into a martingale part and a finite variation part uniquely. For more background on semimartingales and stochastic integration, see e.g. Protter (2005). Semimartingales have several nice properties. In particular:

  • For a suitably bounded predictable process X and a semimartingale Y, the stochastic integral \(\int X dY\) is again a semimartingale.

  • For a twice differentiable function f, the quantity f(Y) is again a semimartingale.

Note that the stochastic integral is of the Stratonovich type. This stochastic integral has the benefit that the stochastic calculus associated with it uses the usual chain rule and product rule, which are both required to derive the coordinate free stochastic Euler–Poincaré theorem. The Stratonovich integral relates to the Itô integral via the following transformation. Let X and Y both be continuous semimartingales, then

$$\begin{aligned} \int _0^t X_s\circ dY_s = \int _0^t X_s dY_s + \frac{1}{2}\langle X,Y\rangle _t, \end{aligned}$$
(1.4)

where \(\langle X,Y\rangle _t\) denotes the quadratic covariation process of X and Y. In particular, the Itô integral is an adapted process and also the quadratic covariation process is an adapted process. This means that also the Stratonovich integral is an adapted process. The two-dimensional version of (1.3) is given by

$$\begin{aligned} {\mathsf {d}}{\varvec{\chi }}_{t} := {\mathbf {u}}({\mathbf {x}},t)dt + \sum _{i=1}^M {\varvec{\xi }}_i({\mathbf {x}})\circ dW_t^i. \end{aligned}$$
(1.5)

Hereafter, we will use Einstein’s summation convention (summing repeated indices over their range) to keep the notation compact.

Boundary conditions. Having defined these vector fields for fluid transport, we can now specify the boundary conditions for the domain illustrated in Fig. 2 above. One assumes that the free surface at the top is a Lagrangian surface, and that no fluid penetrates the bottom and vertical walls. Consequently, the following stochastic kinematic boundary conditions hold for the vertical velocity

$$\begin{aligned} wdt+{\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\circ dW_t^i= & {} {\mathsf {d}}\zeta + ({\mathsf {d}}{\varvec{\chi }}_t\cdot \nabla )\zeta \quad \text { at } z=\zeta ({\mathbf {x}},t), \quad \text { and }\nonumber \\ wdt+{\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\circ dW_t^i= & {} -({\mathsf {d}}{\varvec{\chi }}_t \cdot \nabla )h \quad \text { at } z = -h({\mathbf {x}})\,. \end{aligned}$$
(1.6)

Here, \({\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\) selects the vertical component of the data vector fields. Since the stochastic flow does not penetrate the lateral boundaries, the horizontal velocity is taken to be tangential to the lateral boundaries

$$\begin{aligned} {\mathsf {d}}{\varvec{\chi }}_{t}\cdot {\hat{\mathbf {n}}} = 0, \quad \text { on any vertical lateral boundary,} \end{aligned}$$
(1.7)

where \({\hat{\mathbf {n}}}\) is the unit vector normal to the lateral boundaries. Finally, we assume the dynamic boundary condition for the pressure, namely,

$$\begin{aligned} p = 0 \qquad \text{ at } \quad z = \zeta ({\mathbf {x}},t), \end{aligned}$$
(1.8)

or, alternatively, one can take \(p=\zeta \) at \(z=0\). This condition means that at the free surface the pressure is purely hydrostatic. In this formulation, surface tension has been neglected and the ambient pressure has been set to be zero at the surface. The lateral boundary condition is consistent with the incompressibility condition

$$\begin{aligned} \nabla _3\cdot {\mathsf {d}}\varvec{\chi }_{3t} = 0. \end{aligned}$$
(1.9)

We want to be able to recover the deterministic fluid equations upon removing the stochastic terms in (1.3) and (1.5), each of which is the sum of a deterministic vector and a stochastic vector. That is, the stochastic fluid equations must return to the deterministic fluid equations when the noise term on the transport velocity is switched off. This type of consideration will be repeated as a ‘sanity check’ throughout the paper. For example, this consideration requires that both terms in the transport vector field in (1.3) must be divergence-free,

$$\begin{aligned} \nabla _3\cdot {\mathbf {u}}_3 = 0, \text { and } \nabla _3\cdot {\varvec{\xi }}_i = 0, \hbox { for all } i=1,\dots ,M. \end{aligned}$$
(1.10)

We will assume that the free surface and the pressure are both semimartingales.

2.1 Stochastic Euler–Poincaré Theorem and Averaging

Variational derivatives of functionals.

Definition 1.2

(Functionals and functional derivatives).

A functional \(F[\rho ]\) is defined as a map \(F: \mathcal {B}\rightarrow \mathbb {R}\), where \(\mathcal {B}\) is a Banach space. The variational derivative of a functional \(F(\rho )\), denoted \(\delta F/ \delta \rho \), is defined by

$$\begin{aligned} \delta F[\rho ] := \lim _{\varepsilon \rightarrow 0}\frac{F[\rho +\varepsilon \phi ]-F[\rho ]}{\varepsilon } =: \frac{d}{d\varepsilon }F[\rho +\varepsilon \phi ]\bigg |_{\epsilon =0} = \int _\Omega \frac{\delta F}{\delta \rho }(x) \phi (x) \; dx =: \left\langle \frac{\delta F}{\delta \rho }\,,\,\phi \right\rangle \end{aligned}$$
(1.11)

where \(\varepsilon \in \mathbb {R}\) is a real parameter, \(\phi \in \mathcal {B}\) is arbitrary and the angle brackets \(\langle \,\cdot \,,\,\cdot \,\rangle \) indicate \(L^2\) pairing on \(\mathcal {B}\). The derivative itself can be interpreted as a Fréchet derivative. The function \( \phi (x)\) above is called the ‘variation of \(\rho \)’, and it will be denoted as \(\delta \rho := \phi (x) \). For notational convenience, we denote the functional derivative \(\delta \) operationally as

$$ \delta F[\rho ] = \left\langle \frac{\delta F}{\delta \rho }\,,\,\delta \rho \right\rangle . $$

Euler–Poincaré theorem. Given the boundary conditions and definitions above, the following form of the Euler–Poincaré theorem with stochastic variations provides the corresponding stochastic equations of motion derived from Hamilton’s principle with a deterministic Lagrangian functional \(\ell :\mathfrak {X}\times V^*\rightarrow \mathbb {R}\) defined on the domain of flow, \(\Omega \). Here, \(\mathfrak {X}\) denotes the Lie algebra of smooth vector fields whose action in three-dimensional space by the Jacobi–Lie bracket is denoted as \([\,\cdot \,,\,\cdot ]:\mathfrak {X}\times \mathfrak {X}\rightarrow \mathfrak {X}\), and is defined for \(u,v\in \mathfrak {X}\) by the commutator relation, which in turn defines the minus adjoint operator, \(\mathrm{ad}\), given by

$$\begin{aligned} \big [u,v\big ] := \big ( ({\mathbf {u}}_3\cdot \nabla _3){\mathbf {v}}_3 - ({\mathbf {v}}_3\cdot \nabla _3){\mathbf {u}}_3 \big )\cdot \nabla _3 =: - \,\mathrm{ad}_uv \,. \end{aligned}$$
(1.12)

Theorem 1.1

(Stochastic Euler–Poincaré equations (Holm 2015; Bethencourt de Leon et al. 2020)). The following two statements are equivalent:

  1. i)

    Hamilton’s variational principle in Eulerian coordinates, with \({\mathbf {u}}_3\in \mathfrak {X}(\Omega )\) and \(b, D\in V^*(\Omega )\),

    $$\begin{aligned} \delta S := \delta \int _{t_1}^{t_2}\ell ({\mathbf {u}}_3,b,D) \,dt= 0, \end{aligned}$$
    (1.13)

    holds on \(\mathfrak {X}(\Omega )\times V^*\), upon using variations of the form

    $$\begin{aligned}&\delta {\mathbf {u}}_3 \,dt = {\mathsf {d}}{\mathbf {v}}_3 + [{\mathsf {d}}{\varvec{\chi }}_{3t}, {\mathbf {v}}_3],\qquad \delta b \,dt = -({\mathbf {v}}_3\cdot \nabla _3) b \,dt,\nonumber \\&\qquad \delta D\, dt= -\nabla _3\cdot (D{\mathbf {v}}_3)dt\,, \end{aligned}$$
    (1.14)

    where the arbitrary vector field \({\mathbf {v}}_3\) is a semimartingale.

  2. ii)

    The stochastic Euler–Poincaré equations hold. These equations are

    $$\begin{aligned}&{\mathsf {d}}\frac{\delta \ell }{\delta {\mathbf {u}}_3} + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)\frac{\delta \ell }{\delta {\mathbf {u}}_3} + (\nabla _3{\mathsf {d}}{\varvec{\chi }}_{3t})\cdot \frac{\delta \ell }{\delta {\mathbf {u}}_3} + \frac{\delta \ell }{\delta {\mathbf {u}}_3} (\nabla _3\cdot {\mathsf {d}}{\varvec{\chi }}_{3t}) = -\frac{\delta \ell }{\delta b}\nabla _3 b \,dt\nonumber \\&+ D\nabla _3\frac{\delta \ell }{\delta D}dt \end{aligned}$$
    (1.15)

    or, equivalently,

    $$\begin{aligned}&{\mathsf {d}}\frac{\delta \ell }{\delta {\mathbf {u}}_3} - {\mathsf {d}}{\varvec{\chi }}_{3t}\times \left( \nabla _3 \times \frac{\delta \ell }{\delta {\mathbf {u}}_3} \right) + \nabla _3 \left( {\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \frac{\delta \ell }{\delta {\mathbf {u}}_3} \right) + \frac{\delta \ell }{\delta {\mathbf {u}}_3} (\nabla _3\cdot {\mathsf {d}}{\varvec{\chi }}_{3t}) \nonumber \\&\quad = -\frac{\delta \ell }{\delta b}\nabla _3 b\, dt + D\nabla _3\frac{\delta \ell }{\delta D}dt \,, \end{aligned}$$
    (1.16)

    with advection equations

    $$\begin{aligned} {\mathsf {d}}b = - \,{\mathsf {d}}{\varvec{\chi }}_{3t} \cdot \nabla b \quad \hbox {and}\quad {\mathsf {d}}D = -\, \nabla _3\cdot (D{\mathsf {d}}{\varvec{\chi }}_{3t}) \,. \end{aligned}$$
    (1.17)

Remark 1.2

The abstract statement of the Euler–Poincaré Theorem 1.1, formulated on general semidirect product Lie groups, is presented in Holm et al. (1998) deterministically and in Holm (2015); Bethencourt de Leon et al. (2020) stochastically.

Remark 1.3

In Theorem 1.1, the operator \(\delta \) in (1.13) is the functional derivative defined in (1.11), the brackets \([\,\cdot \,,\,\cdot \,]\) denote the commutator of vector fields defined in (1.12), and \({\mathbf {v}}_3 \in \mathfrak {X}(\Omega )\) is an arbitrary semimartingale vector field in three dimensions which vanishes at the endpoints in time, \(t_1\) and \(t_2\). Note that the stochasticity is introduced in the variation of the Eulerian velocity in (1.14). This stochasticity in the variation is inherited from the stochasticity in the Lagrangian particle paths, as in Holm (2015).

Remark 1.4

(Newton’s Law interpretation of Euler–Poincaré equation (1.15)). One may interpret the stochastic Euler–Poincaré equation (1.15) as the Newton’s law of motion for a stochastic process. That is, the stochastic rate of change of the covector momentum \({\mathbf {P}}:=\delta \ell /\delta {\mathbf {u}}_3\) equals the sum of forces on the right hand side of Eq. (1.15). Of course, when the stochasticity is removed from the vector field in (1.5), Eq. (1.15) recovers its deterministic version.

Proof

Hamilton’s variational principle implies

$$\begin{aligned} 0&= \int _{t_1}^{t_2}\left[ \left\langle \frac{\delta \ell }{\delta {\mathbf {u}}_3},\delta {\mathbf {u}}_3 dt\right\rangle _\mathfrak {X} + \left\langle \frac{\delta \ell }{\delta b},\delta b dt\right\rangle _{V^*} + \left\langle \frac{\delta \ell }{\delta D},\delta D dt\right\rangle _{V^*}\right] \\&= \int _{t_1}^{t_2}\left[ \left\langle \frac{\delta \ell }{\delta {\mathbf {u}}_3},{\mathsf {d}}{\mathbf {v}}_3 + [{\mathsf {d}}{\varvec{\chi }}_{3t}, {\mathbf {v}}_3]\right\rangle _\mathfrak {X} + \left\langle \frac{\delta \ell }{\delta b},-({\mathbf {v}}_3\cdot \nabla _3)b dt\right\rangle _{V^*} + \left\langle \frac{\delta \ell }{\delta D},-\nabla _3\cdot (D {\mathbf {v}}_3)dt\right\rangle _{V^*}\right] \\&= \int _{t_1}^{t_2}\left[ \left\langle -{\mathsf {d}}\frac{\delta \ell }{\delta {\mathbf {u}}_3}-({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)\frac{\delta \ell }{\delta {\mathbf {u}}_3} - (\nabla _3{\mathsf {d}}\varvec{\chi }_{3t})\cdot \frac{\delta \ell }{\delta {\mathbf {u}}_3} + \frac{\delta \ell }{\delta {\mathbf {u}}_3}(\nabla _3\cdot {\mathsf {d}}{\varvec{\chi }}_{3t}),{\mathbf {v}}_3\right\rangle _\mathfrak {X} \right. \\&\quad \left. + \left\langle -\frac{\delta \ell }{\delta b}\nabla _3 b dt, {\mathbf {v}}_3\right\rangle _\mathfrak {X}+ \left\langle D\nabla _3\frac{\delta \ell }{\delta D} dt, {\mathbf {v}}_3\right\rangle _\mathfrak {X} \right] . \end{aligned}$$

The subscripts \(\mathfrak {X}\) and \(V^*\) on the \(L^2\) pairings indicate over which space that the pairing is defined. Since the semimartingale \({\mathbf {v}}_3\) is arbitrary, except for vanishing at the endpoints \(t_1\) and \(t_2\) in time, the following equation holds,

$$\begin{aligned} {\mathsf {d}}\frac{\delta \ell }{\delta {\mathbf {u}}_3} + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)\frac{\delta \ell }{\delta {\mathbf {u}}_3} + (\nabla _3{\mathsf {d}}{\varvec{\chi }}_{3t})\cdot \frac{\delta \ell }{\delta {\mathbf {u}}_3} + \frac{\delta \ell }{\delta {\mathbf {u}}_3}(\nabla _3\cdot {\mathsf {d}}{\varvec{\chi }}_{3t})= -\frac{\delta \ell }{\delta b}\nabla _3 b \,dt + D\nabla _3\frac{\delta \ell }{\delta D}\,dt. \end{aligned}$$

This finishes the proof of the stochastic Euler–Poincaré equation in (1.15). The equivalent form in Eq. (1.16) follows by means of a standard vector identity. \(\square \)

2.2 Stochastic Kelvin–Noether Circulation Theorem

A straight forward calculation combining equation (1.15) and the second advection equation in (1.17) proves the following.

Lemma 1.5

(Circulation form of the stochastic Euler–Poincaré equation (Holm 2015; Bethencourt de Leon et al. 2020)). The stochastic Euler–Poincaré equation in (1.15) is equivalent to the following,

$$\begin{aligned}&{\mathsf {d}}\left( \frac{1}{D}\frac{\delta \ell }{\delta {\mathbf {u}}_3}\right) + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)\left( \frac{1}{D}\frac{\delta \ell }{\delta {\mathbf {u}}_3}\right) + (\nabla _3{\mathsf {d}}{\varvec{\chi }}_{3t})\cdot \left( \frac{1}{D}\frac{\delta \ell }{\delta {\mathbf {u}}_3}\right) \nonumber \\&= -\,\frac{1}{D}\frac{\delta \ell }{\delta b}\nabla _3 b \,dt + \nabla _3\frac{\delta \ell }{\delta D}dt \,. \end{aligned}$$
(1.18)

One of the main benefits of Theorem 1.1 is that its stochastic Euler–Poincaré equations satisfy the following Kelvin circulation theorem.

Theorem 1.6

(Stochastic Kelvin–Noether circulation theorem (Holm 2015; Bethencourt de Leon et al. 2020)). For an arbitrary loop \(c({\mathsf {d}}{\varvec{\chi }}_{3t})\) which is advected by the stochastic velocity field \({\mathsf {d}}{\varvec{\chi }}_{3t}\), the following circulation dynamics holds

$$\begin{aligned} {\mathsf {I}} := \oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})}\frac{1}{D}\frac{\delta \ell }{\delta {\mathbf {u}}_3}\cdot d{\mathbf {x}}_3 \,,\qquad {\mathsf {d}} {\mathsf {I}} = -\oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})} \left( \frac{1}{D}\frac{\delta \ell }{\delta b}\right) \nabla _3 b\cdot d{\mathbf {x}}_3 \,dt\,. \end{aligned}$$
(1.19)

Proof

The Kelvin circulation law (1.19) follows from Newton’s law of motion obtained from the stochastic Euler–Poincaré equation (1.18) for the evolution of momentum/mass \(D^{-1}\delta \ell /\delta {\mathbf {u}}_3\) concentrated on an advecting material loop, \(c({\mathsf {d}}{\varvec{\chi }}_{3t}) = \phi _t c(0)\), where \(\phi _t\) is the stochastic flow map generated by the stochastic vector field \({\mathsf {d}}{\varvec{\chi }}_{3t}\) defined in Eq. (1.3). Upon changing variables to pull back the integrand to its initial position, the stochastic differential can be moved inside and the Kunita–Itô-Wentzell formula may be applied (Bethencourt de Leon et al. 2020). Then, by inverting the pull-back we have the following

$$\begin{aligned} {\mathsf {d}}\oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})}\frac{1}{D}\frac{\delta \ell }{\delta {\mathbf {u}}_3}\cdot d{\mathbf {x}}_3&= \oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})} ({\mathsf {d}} + {\mathsf {d}}\varvec{\chi }_{3t}\cdot \nabla _3 + (\nabla _3{\mathsf {d}}\varvec{\chi }_{3t})\cdot )\left( \frac{1}{D}\frac{\delta \ell }{\delta {\mathbf {u}}_3}\right) \cdot d{\mathbf {x}}_3 \\&= -\oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})} \frac{1}{D}\frac{\delta \ell }{\delta b}\nabla _3 b \,dt\cdot d{\mathbf {x}}_3 + \oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})}\nabla _3 \frac{\delta \ell }{\delta D} \cdot d{\mathbf {x}}_3\,dt \\&= -\oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})} \left( \frac{1}{D}\frac{\delta \ell }{\delta b}\right) \nabla _3 b \cdot d{\mathbf {x}}_3 \,dt\,. \end{aligned}$$

In the second line, we have used the Euler–Poincaré equation (1.15) and the advection equation for the density. The last step applies the fundamental theorem of calculus to show vanishing of the last loop integral in the second line. For the corresponding proof in the deterministic case, see Holm et al. (1998). For detailed discussion of pull-back by stochastic flow maps, see Bethencourt de Leon et al. (2020). \(\square \)

Corollary 1.7

(Generation of circulation, \({\mathsf {I}}\)). By Stokes Law, Eq. (1.19) in the stochastic Kelvin–Noether circulation theorem 1.6 implies

$$\begin{aligned} {\mathsf {d}} {\mathsf {I}} = -\int \!\!\int _{\partial S = c({\mathsf {d}}{\varvec{\chi }}_{3t})} \nabla _3\left( \frac{1}{D}\frac{\delta \ell }{\delta b}\right) \times \nabla _3 b\cdot d{\mathbf {S}}_3 \,dt\,. \end{aligned}$$
(1.20)

Therefore, circulation is created by misalignment of the gradients of buoyancy b and its dual quantity \(D^{-1}\delta \ell /\delta b\).

Remark 1.8

(A mechanism for cyclogenesis). Formula (1.20) expresses the mechanism for generation of circulation (i.e. convection) driven by misalignment of certain potential gradients with gradients of scalar advected fluid quantities such as the buoyancy, b. In particular, formula (1.20) is the fundamental mechanism for generation of circulation or convection by wave–current interaction in stratified fluids. For the vertically averaged stratified fluid models treated later in the present paper, this formula will express a barotropic mechanism for generating horizontal circulation by misalignment of horizontal gradients of certain barotropic fluid quantities (such as wave elevation or bottom topography) with the horizontal gradient of vertically averaged buoyancy.

In three-dimensional stochastic fluid dynamics, the Lagrangian in the Euler–Poincaré theorem is a functional defined over the volume of flow which, as we will discuss below, involves the kinetic energy density of the fluid relative to the rotating frame and the potential energy density. Our aims in the remainder of the paper are to combine asymptotic expansions and vertical averaging with the stochastic Euler–Poincaré variational theorem to formulate a new approach for developing stochastic parametrisation methods. To achieve these aims, we will apply asymptotic expansions in a vertically averaged (barotropic) stochastic Euler–Poincaré variational principle. For this purpose, we will apply asymptotic expansions to the nondimensionalised Lagrangian for 3D incompressible flows of a stratified and rotating Euler fluid and then evaluate the vertical integral at an appropriate order in the expansion and finally use the Euler–Poincaré theorem to derive the equations of motion and advection we seek. We will then analyse and discuss their solution properties from the viewpoints of Newton’s laws of motion and the Kelvin–Noether circulation theorem. We will also discuss the conservation laws for these equations.

2.3 Nondimensionalising the Lagrangian

The dimensional form of the Lagrangian in Hamilton’s principle for the rotating, stratified Euler equations (rsE) is given by

$$\begin{aligned} \ell _{rsE}({\mathbf {u}}_3,b,\zeta ,D) := \int _\Omega \rho _0 (1+b)\left( \frac{1}{2}|{\mathbf {u}}|^2 + \frac{1}{2}w^2 + {\mathbf {u}}\cdot {\mathbf {R}} - gz\right) D\,dx\,dy\,dz. \nonumber \\ \end{aligned}$$
(1.21)

Here, \(\rho _0\) represents the reference density and g represents gravity. The ocean has quite a few small dimensionless numbers which can be used to simplify the rsE Lagrangian and will allow one to access a hierarchy of simplified models. In particular, we want to derive the Lagrangian for the Euler–Boussinesq equations, which requires assumptions on the smallness of buoyancy, in terms of the Rossby number. To derive the equations of motion associated with the Lagrangian, we introduce the following action

$$\begin{aligned} S_{rsE} = \int _{t_1}^{t_2} \ell _{rsE} \,dt - \langle {\mathsf {d}}p, D-1\rangle =: \int _{t_1}^{t_2} c\ell _{rsE}, \end{aligned}$$
(1.22)

where \({\mathsf {d}}p\) is the Lagrange multiplier that enforces the density ratio D to be equal to one, the times \(0\le t_1 \le t_2\) are arbitrary, and the angle brackets refer to the \(L^2\) pairing over the domain \(\Omega \). The notation \(c\ell _{rsE}\) refers to constrained Lagrangian and is introduced to keep the notation similar to the stochastic Euler–Poincaré theorem 1.1. This constraint implies incompressibility and is required because it affects the measure \(D\,dx\,dy\,dz\) in the Lagrangian. The treatment of the stochastic pressure is explained in the following remark.

Remark 1.9

(Semimartingale pressure). At this point one recognises a departure from the stochastic Euler–Poincaré equations without constraints derived in the Euler–Poincaré theorem 1.1. Namely, we have written the Lagrange multiplier \({\mathsf {d}}p\) which imposes the constraint \(D-1=0\). The notation stresses that \({\mathsf {d}}p\) is imposing a constraint that is stochastic. Now, setting \(D=1\) in the advection equation for D by the stochastic vector field \({\mathsf {d}}{\varvec{\chi }}_{3t}\) implies that \(\nabla _3\cdot ( {\mathsf {d}}{\varvec{\chi }}_{3t})=0\). Following the discussion leading to (1.10), this in turn must also imply \(\nabla _3\cdot {\mathbf {u}}_3=0\). By its definition in (1.5), the quantity \({\varvec{\chi }}_{3t}\) is a semimartingale. Therefore, accounting for both the deterministic and stochastic parts of the motion equation in (1.33) will require that the pressure \({\mathsf {d}}p\) must also be a semimartingale, hence the notation. The point is that the semimartingale D cannot be enforced to be a constant by a deterministic Lagrange multiplier. The Lagrange multiplier must also be obtained from a semimartingale equation. In the present case, this can be accomplished by acknowledging that the pressure is a semimartingale and writing its contribution in the motion equation as \({\mathsf {d}}p\), in a notation which implies a sum of both Lebesque and stochastic time integrations. Then, upon imposing the consequence of \(D=1\) in the form \(\nabla _3\cdot {\mathbf {u}}_3=0\) we find a semimartingale Poisson equation for \({\mathsf {d}}p\) which encompasses both the deterministic and stochastic parts of the constrained motion equation. Finally, the time integration of the solution of the Poisson equation for \({\mathsf {d}}p\) determines the semimartingale p. For a treatment of general semimartingale driven variational principles, see Street and Crisan (2020).

The nondimensional versions of all the relevant variables and parameters are given below,

figure a

Here L denotes the horizontal scale, H is the vertical scale, U is the typical horizontal velocity, \(f_0\) is the rotation frequency, \(\zeta _0\) is the typical free surface amplitude and T is the time scale. The dimensionless numbers in the bottom row are, respectively, the aspect ratio \(\sigma \), the wave amplitude \(\alpha \), the Froude number \(\mathrm{Fr}\), the Rossby number \(\mathrm{Ro}\) and the Strouhal number \(\mathrm{Sr}\). Note that we have also scaled the Brownian motion so that in the nondimensional setting, the noise is again a standard Brownian motion. The dimensional factor that arises can be absorbed into the \({\varvec{\xi }}_{3i}\) for each i. The vertical component of the data vector fields \({\varvec{\xi }}_{3i}\) is scaled with the aspect ratio as well, that is \({\varvec{\xi }}_{3i} = ({\varvec{\xi }}_i',\sigma {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}')\). In particular, this means that we can write

$$\begin{aligned} {\mathsf {d}}{\varvec{\chi }}_{3t} = U({\mathsf {d}}{\varvec{\chi }}_t', \sigma {\hat{\mathbf {z}}}\cdot {\mathsf {d}}{\varvec{\chi }}_{3t}'). \end{aligned}$$
(1.24)

We do not make any assumptions on the size of the data vector fields relative to the velocity field itself. Last, but not least, the stratification parameter \(\mathfrak {s}\) is introduced. Since the buoyancy is already dimensionless, it does not appear in the table, but it works as follows

$$\begin{aligned} b = \mathfrak {s}b'. \end{aligned}$$
(1.25)

The purpose of the stratification parameter is to make sure that the buoyancy variable b is an order \(\mathcal {O}(1)\) variable. By controlling the size of the stratification parameter \(\mathfrak {s}\), the Boussinesq approximation can be introduced. The nondimensional rsE Lagrangian is obtained by substituting (1.23) into (1.21) and dropping the primes, which yields

$$\begin{aligned} \ell _{rsE}({\mathbf {u}}_3,b,D) = \int _{\Omega } (1+\mathfrak {s}b)\left( \frac{1}{2}|{\mathbf {u}}|^2 + \frac{\sigma ^2}{2}w^2 + \frac{1}{\mathrm{Ro}}{\mathbf {u}}\cdot {\mathbf {R}}- \frac{1}{\mathrm{Fr}^2}z\right) D\,dx\,dy\,dz, \nonumber \\ \end{aligned}$$
(1.26)

and the dimensionless action is given by

$$\begin{aligned} S_{rsE} = \int _{t_1}^{t_2}\ell _{rsE}\,dt - \left\langle \frac{1}{\mathrm{Fr}^2}{\mathsf {d}}p,D-1\right\rangle =: \int _{t_1}^{t_2}c\ell _{rsE}. \end{aligned}$$
(1.27)

In the ocean, the horizontal scale L is of the order of hundreds of kilometres, whereas the vertical scale H is typically about four kilometres. The free surface amplitude is five metres and the horizontal velocity is about a tenth of a metre per second. Hence, the aspect ratio \(\sigma \ll 1\), the wave amplitude \(\alpha \ll 1\) and the Froude number \(\mathrm{Fr}\ll 1\). The Rossby number at these scales is also small, \(\mathrm{Ro}\ll 1\). Also, the buoyancy stratification is weak, which allows us to apply the Boussinesq approximation. This approximation corresponds to \(\mathfrak {s}\ll 1\), that is, requiring the stratification parameter to be small. When the dimensionless numbers satisfy \(\mathcal {O}(\alpha ) = \mathcal {O}(\mathfrak {s}) = \mathcal {O}(\mathrm{Ro}) = \mathcal {O}(\mathrm{Fr}) = \mathcal {O}(\sigma ^2)\), the Lagrangian can be approximated. Consequently, the rsE Lagrangian simplifies, as the remaining effect of buoyancy is restricted to the potential energy term. This yields the Euler–Boussinesq (EB) Lagrangian, given by

$$\begin{aligned} \ell _{EB}({\mathbf {u}}_3,b,D) = \int _{\Omega } \left( \frac{1}{2}|{\mathbf {u}}|^2 + \frac{\sigma ^2}{2}w^2 + \frac{1}{\mathrm{Ro}}{\mathbf {u}}\cdot {\mathbf {R}} - \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}b)z\right) D\,dx\,dy\,dz. \nonumber \\ \end{aligned}$$
(1.28)

The Euler–Boussinesq equations are obtained by applying the Euler–Poincaré theorem to the action obtained by taking the Lagrangian in (1.28) with the pressure constraint, as in (1.22). The action for the EB equations is then given by

$$\begin{aligned} S_{EB} = \int _{t_1}^{t_2}\ell _{EB} \, dt - \left\langle \frac{1}{\mathrm{Fr}^2}{\mathsf {d}}p, D-1\right\rangle =: \int _{t_1}^{t_2} c\ell _{EB}. \end{aligned}$$
(1.29)

Besides assuming the buoyancy is small, we will assume that the variations of the Coriolis parameter and of the bathymetry profile are also small, of order \(O(\mathrm{Ro})\),

$$\begin{aligned} f({\mathbf {x}}) = 1 + \mathrm{Ro}\, f_1({\mathbf {x}}), \qquad h({\mathbf {x}}) = 1 + \mathrm{Ro}\, h_1({\mathbf {x}}). \end{aligned}$$
(1.30)

These assumptions are made because they are consistent with the assumptions for quasi-geostrophy. The Lagrangian of interest in (1.28) is in dimensionless form, but the constraints in theorem 1.1 are still dimensional. Since \({\mathbf {v}}_3\) is arbitrary, multiplying it by some constant does not change its arbitrary nature. Hence, besides the \(\delta {\mathbf {u}}_3\) constraint, nothing changes upon nondimensionalisation. As said earlier, the \(\delta {\mathbf {u}}_3\) variational constraint does change, as follows,

$$\begin{aligned} \delta {\mathbf {u}}_3 dt = \mathrm{Sr}\,{\mathsf {d}}{\mathbf {v}}_3 - [{\mathsf {d}}{\varvec{\chi }}_{3t},{\mathbf {v}}_3]. \end{aligned}$$
(1.31)

Time does not appear explicitly anywhere in the rsE and EB Lagrangians. Thus, the Strouhal number has not appeared before; but time rescaling has a significant impact on the behaviour of the model. In (1.31), one can see that if the Strouhal number is not unity, advection will no longer be balanced. This observation will be crucial later, when we look at the short time limit. So far, we have obtained a theorem which, for a certain deterministic Lagrangian for three-dimensional fluids, provides us with the corresponding stochastic equations. By explicitly evaluating the vertical integral, when possible, in that theorem, we have a systematic way to obtain the vertically averaged version of the three-dimensional fluid equations of interest. We also have introduced a general nondimensionalisation and identified the scales in the problem which determine the small dimensionless numbers in the ocean. Now, an application of theorem 1.1 to the EB Lagrangian (1.29), with variations given by

$$\begin{aligned} \begin{aligned} \frac{\delta c\ell _{EB}}{\delta {\mathbf {u}}}&= D\left( {\mathbf {u}} + \frac{1}{\mathrm{Ro}}{\mathbf {R}}\right) ,\\ \frac{\delta c\ell _{EB}}{\delta w}&= \sigma ^2 D w,\\ \frac{\delta c\ell _{EB}}{\delta D}&= \frac{1}{2}|{\mathbf {u}}|^2 + \frac{\sigma ^2}{2}w^2 + \frac{1}{\mathrm{Ro}}{\mathbf {u}}\cdot {\mathbf {R}} - \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}b)z - \frac{1}{\mathrm{Fr}^2}{\mathsf {d}}p,\\ \frac{\delta c\ell _{EB}}{\delta b}&= -\frac{\mathfrak {s}}{\mathrm{Fr}^2}Dz,\\ \frac{\delta c\ell _{EB}}{\delta {\mathsf {d}}p}&= \frac{1}{\mathrm{Fr}^2} (D-1). \end{aligned} \end{aligned}$$
(1.32)

implies the following stochastic Euler–Poincaré equations in circulation form (see Lemma 1.5)

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\mathsf {d}}{\mathbf {u}} + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3){\mathbf {u}} + (\nabla \varvec{\xi }_{3i})\cdot {\mathbf {u}}_3\circ dW_t^i&= -\frac{1}{\mathrm{Fr}^2} \nabla {\mathsf {d}}p - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times {\mathsf {d}}{\varvec{\chi }}_t \\&\quad - \frac{1}{\mathrm{Ro}}\nabla ({\varvec{\xi }}_i\cdot {\mathbf {R}})\circ dW_t^i,\\ \sigma ^2\left( \mathrm{Sr}\,{\mathsf {d}}w + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)w+ \Big (\frac{\partial }{\partial z}{\varvec{\xi }}_{3i}\Big )\cdot {\mathbf {u}}_3\circ dW_t^i\right)&= -\frac{1}{\mathrm{Fr}^2}\frac{\partial }{\partial z}{\mathsf {d}}p + \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}b)\,dt,\\ \mathrm{Sr}\,{\mathsf {d}}b + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)b&= 0,\\ \nabla _3\cdot ( {\mathsf {d}}{\varvec{\chi }}_{3t})&= 0. \end{aligned} \end{aligned}$$
(1.33)

The Euler–Boussinesq equations satisfy the following Kelvin circulation theorem, for any closed loop \(c({\mathsf {d}}{\varvec{\chi }}_{3t})\) which is advected with the stochastic velocity \({\mathsf {d}}{\varvec{\chi }}_{3t}\) in Eq. (1.3),

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})}\left( ({\mathbf {u}},\sigma ^2 w) + \frac{1}{\mathrm{Ro}}({\mathbf {R}},0)\right) \cdot d{\mathbf {x}}_3&= -\frac{\mathfrak {s}}{\mathrm{Fr}^2}\oint _{c({\mathsf {d}}{\varvec{\chi }}_{3t})} z\nabla _3 b\cdot d{\mathbf {x}}_3 dt\\&= -\frac{\mathfrak {s}}{\mathrm{Fr}^2}\int \!\!\int _{\partial S=c({\mathsf {d}}{\varvec{\chi }}_{3t})} {\hat{\mathbf {z}}}\times \nabla _3 b\cdot d{\mathbf {S}} dt, \end{aligned} \end{aligned}$$
(1.34)

where the notation \(({\mathbf {u}}, \sigma ^2 w)\) denotes a three-dimensional vector field, two horizontal components from \({\mathbf {u}}\) and the vertical component \(\sigma ^2 w\). As \({\mathbf {R}}\) is strictly horizontal, the vertical component is zero. Hence, the misalignment of the unit vector in the vertical direction and the gradient of buoyancy creates vertical circulation, or convection.

Additionally, the Euler–Boussinesq equations satisfy the Silberstein–Ertel theorem for potential vorticity. This theorem states that the potential vorticity, defined by

$$\begin{aligned} q := \mathfrak {s}\nabla _3 b \cdot \nabla _3\times \left( ({\mathbf {u}},\sigma ^2 w) + \frac{1}{\mathrm{Ro}}({\mathbf {R}},0)\right) , \end{aligned}$$
(1.35)

is conserved along particle trajectories and thus satisfies the following equation

$$\begin{aligned} \mathrm{Sr}\,{\mathsf {d}}q + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)q = 0. \end{aligned}$$
(1.36)

Since the buoyancy and the potential vorticity are constant along particle trajectories, the spatially integrated quantity,

$$\begin{aligned} C_\Phi = \int _{\Omega } \Phi (b,q) \,dx\,dy\,dz, \end{aligned}$$
(1.37)

is also preserved in time for any differentiable function, \(\Phi \), for which the integral exists. The proof is analogous to the deterministic case, which is shown in Holm et al. (1998, 1999). A special case of this statement is the preservation of the enstrophy, which is defined as the \(L^2\) norm of the potential vorticity. Since the flow is divergence free, one can also define the enstrophy in terms of the gradients of the velocity. This shows that the Euler–Boussinesq equations, even in the presence of SALT, have an infinite number of conservation laws. This structure must also be preserved by the vertical averaging. The spatially integrated quantities \(C_\Phi \) are also referred to as Casimirs, as they are the functions whose Lie–Poisson bracket corresponding to the Euler–Boussinesq equations vanishes for any Hamiltonian expressed in the Eulerian fluid variables.

2.4 Averaging of Newton’s Second Law

Besides evaluating the vertical integral in the variational principle, one can also choose to use Newton’s second law to derive the equations of fluid motion in this domain, rather than using the Euler–Poincaré theorem. By means of the method of control volumes, it is possible to derive the equations and also come up with an averaging principle. This is what is shown in Wu (1981) for the deterministic case. The stochastic case is not that different, but there is one issue that requires careful treatment: there is an additional advection term. Let us denote the vertical average by putting a bar over the relevant quantity

$$\begin{aligned} \overline{f} := \frac{1}{\eta }\int _{-h}^{\alpha \zeta } f dz. \end{aligned}$$
(1.38)

The stochastic vector field in the averaged situation is denoted

$$\begin{aligned} \overline{{\mathsf {d}}{\varvec{\chi }}}_t = \overline{{\mathbf {u}}}\,dt + \overline{{\varvec{\xi }}}_i\circ dW_t^i. \end{aligned}$$
(1.39)

For incompressible flows, the advection equation for a scalar and the continuity equation for a density can be written in the same form. That is, the average of a scalar function \(f({\mathbf {x}}_3,t)\) and that of a volume form \(f({\mathbf {x}}_3,t) d^3x\), for incompressible flows, are of the same form,

$$\begin{aligned} \mathrm{Sr}\, {\mathsf {d}}\int _{-h}^\zeta f({\mathbf {x}}_3,t) dz + \nabla \cdot \int _{-h}^{\alpha \zeta } f({\mathbf {x}}_3,t){\mathsf {d}}{\varvec{\chi }}_t dz = 0. \end{aligned}$$
(1.40)

In the deterministic case, it is possible to substitute in the fluid velocity for f in (1.40) and obtain the vertically averaged momentum equation after applying (1.38). The formula above holds for scalars and densities, but fluid velocity is neither. However, the fluid velocity equation obtained in this way is correct, but only in the deterministic case. The explanation for this coincidence is the following. In the deterministic setting, the advective terms in the equation for the fluid velocity for incompressible fluids are \(({\mathbf {u}}\cdot \nabla ){\mathbf {u}} + (\nabla {\mathbf {u}})\cdot {\mathbf {u}}\). The latter term is equal to the gradient of the kinetic energy, so a cancellation occurs in Newton’s second law. When SALT is introduced in this problem, the kinetic energy is the same as in the deterministic situation, but the advective terms are now stochastic; hence, this cancellation no longer occurs.

Applying (1.38) and (1.40) to the Euler–Boussinesq equations (1.33) yields the following vertically averaged nonlinear equations,

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\mathsf {d}} (\eta \overline{{\mathbf {u}}}) + \nabla \cdot (\eta \overline{{\mathsf {d}}{\varvec{\chi }}_{t}\otimes {\mathbf {u}}}) + \eta (\nabla \overline{\varvec{\xi }}_i)\cdot \overline{{\mathbf {u}}} \circ dW_t^i&= -\frac{1}{\mathrm{Fr}^2}\eta \overline{\nabla {\mathsf {d}} p} - \frac{1}{\mathrm{Ro}}\eta f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}_t}\\&\quad - \frac{1}{\mathrm{Ro}}\eta \nabla (\overline{{\varvec{\xi }}}_i\cdot {\mathbf {R}})\circ dW_t^i,\\ \mathrm{Sr}\,{\mathsf {d}}(\eta \overline{b}) + \nabla \cdot (\eta \overline{b {\mathsf {d}}\varvec{\chi }_t})&= 0,\\ \mathrm{Sr}\,{\mathsf {d}}\eta + \nabla \cdot (\eta \overline{{\mathsf {d}}{\varvec{\chi }}_t})&= 0. \end{aligned} \end{aligned}$$
(1.41)

The last equation is obtained by substituting unity into (1.40). It corresponds to conservation of volume in the two-dimensional setting. As the problem is incompressible, the vertical velocity can be expressed in terms of the horizontal velocity field and the vertical component of the data vector fields \({\varvec{\xi }}_{3i}\) can be expressed in terms of the horizontal components as

$$\begin{aligned} \begin{aligned} w({\mathbf {x}},z)&= - \nabla \cdot \int _{-h}^z {\mathbf {u}}({\mathbf {x}},z') dz' = \nabla \cdot \int _z^{\alpha \zeta } {\mathbf {u}}({\mathbf {x}},z')dz',\\ {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}({\mathbf {x}},z)&= - \nabla \cdot \int _{-h}^z {\varvec{\xi }}_i({\mathbf {x}},z')dz' = \nabla \cdot \int _{z}^{\alpha \zeta }{\varvec{\xi }}_i({\mathbf {x}},z')dz'. \end{aligned} \end{aligned}$$
(1.42)

This expression has been derived by vertically integrating the three-dimensional incompressibility condition (1.9), using the uniqueness of the semimartingale decomposition and using the boundary conditions on the vertical velocity to pull the divergence outside of the integral. Importantly, the boundary conditions introduce a dependence between the horizontal components of the vector fields and the vertical component. A horizontal two-dimensional model with bathymetry and a free surface will therefore give some information about what is happening in the vertical direction. This holds also for the \({\varvec{\xi }}_i\). Even though the Newtonian averaging approach is very insightful, there is a drawback. Namely, the averaged equations (1.41) are not closed. Indeed, they contain three terms which are unknown. In the momentum equation, the average of the nonlinear term and the average of the pressure are unknown. In the buoyancy equation, the advection term is unknown. In order to close this set of equation, we will use asymptotic analysis, which we shall employ in two different scaling regimes. These are the long time - very small wave (LT-VSW) scaling regime in Sect. 2 and the short time-small wave (ST-SW) scaling regime in Sect. 3. Here, ‘long time scale’ is \(T = L/U\), the time it takes for a fluid parcel to cross the horizontal length scale; and ‘short time scale’ is \(T = L/\sqrt{gH}\), the time it takes for a gravity wave to cross the horizontal length scale. Likewise, ‘small wave’ means that the wave amplitude is small, but not small enough to consider taking the rigid lid limit; while ‘very small wave’ means that the wave amplitude is the small parameter of interest.

3 Long Time—Very Small Wave Scaling Regime

Long time corresponds to choosing the time scale to be \(T = L/U\) and very small wave means that the amplitude of the wave \(\alpha \) is the small parameter of interest. In this setting, we therefore have the following dimension-free parameters,

figure b

In particular, the velocity field and the data vector fields are scaled in the same way; hence, we have

$$\begin{aligned} {\mathsf {d}}{\varvec{\chi }}_{3t} = U({\mathsf {d}}{\varvec{\chi }}_t', \sigma {\hat{\mathbf {z}}}\cdot {\mathsf {d}}{\varvec{\chi }}_{3t}'). \end{aligned}$$
(2.2)

With these scaling relations and the stratification parameter \(\mathfrak {s}\), the constrained EB Lagrangian in Eq. (1.29) takes the following form

$$\begin{aligned} \begin{aligned} S_{EB}({\mathbf {u}}_3,b,D)&= \int _{t_1}^{t_2}\int _{\Omega } \left( \frac{1}{2}|{\mathbf {u}}|^2 + \frac{\sigma ^2}{2}w^2 + \frac{1}{\mathrm{Ro}}({\mathbf {u}}\cdot {\mathbf {R}}) - \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}b)z\right) D\,dx\,dy\,dz\,dt \\&\quad - \langle {\mathsf {d}}p, D-1\rangle ,\\&=: \int _{t_1}^{t_2} c\ell _{EB} \end{aligned} \end{aligned}$$
(2.3)

Note that no information about the very small free surface amplitude appears in the Lagrangian; it only contains the aspect ratio, which controls the size of the vertical kinetic energy. However, information about the size of the free surface amplitude does appear in the boundary conditions, which are

$$\begin{aligned} \begin{matrix} p = \zeta \qquad &{} \text { at } z = 0,\\ wdt + {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\circ dW_t^i= \alpha \big ({\mathsf {d}}\zeta + ({\mathsf {d}}{\varvec{\chi }}_t\cdot \nabla )\zeta \big ) \qquad &{} \text { at } z = \alpha \zeta ({\mathbf {x}},t),\\ wdt + {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\circ dW_t^i= -({\mathsf {d}}{\varvec{\chi }}_t\cdot \nabla )h \qquad &{} \text { at } z = -h({\mathbf {x}}),\\ {\mathsf {d}}{\varvec{\chi }}_t\cdot {\mathbf {n}} = 0 \qquad &{} \text { on lateral boundaries}. \end{matrix} \end{aligned}$$
(2.4)

An application of the stochastic Euler–Poincaré Theorem 1.1 on the long time scale Lagrangian in (2.3) now yields the following equations

$$\begin{aligned} \begin{aligned} {\mathsf {d}}{\mathbf {u}} + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3){\mathbf {u}} + (\nabla \varvec{\xi }_{3i})\cdot {\mathbf {u}}_3\circ dW_t^i&= -\frac{1}{\mathrm{Fr}^2}\nabla {\mathsf {d}}p - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times {\mathsf {d}}{\varvec{\chi }}_t \\&\quad - \frac{1}{\mathrm{Ro}}\nabla ({\varvec{\xi }}_i\cdot {\mathbf {R}})\circ dW_t^i,\\ \sigma ^2\left( {\mathsf {d}}w + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)w + \Big (\frac{\partial }{\partial z}{\varvec{\xi }}_{3i}\Big )\cdot {\mathbf {u}}_3 \circ dW_t^i\right)&= - \frac{1}{\mathrm{Fr}^2}\frac{\partial }{\partial z}{\mathsf {d}}p - \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}b) dt,\\ \nabla \cdot {\mathsf {d}}{\varvec{\chi }}_{3t}&= 0. \end{aligned} \end{aligned}$$
(2.5)

The equations in (2.5) satisfy the Kelvin circulation theorem as in (1.34) and have conservation of potential vorticity along particle trajectories as in (1.36). These equations also conserve an infinity of integral quantities as in (1.37). In the long-time scaling in (2.1) the Strouhal number is equal to one. In this scaling regime, the equations take a particularly nice form. The dimensionless numbers of interest are the aspect ratio \(\sigma \) and the wave amplitude \(\alpha \), the Rossby number \(\mathrm{Ro}\) shall be left untouched. In particular, we consider \(\alpha \ll \sigma \ll 1\), where we let the wave amplitude tend to zero while holding the aspect ratio fixed.

Rigid lid approximation. The effect of sending the wave amplitude \(\alpha \) to zero is the rigid lid approximation, where the free surface is no longer allowed to vary and becomes a rigid boundary, instead. This removes gravity waves from the problem. However, the leading order dynamics can still be recovered from the dynamic boundary condition on the pressure. The effect of sending \(\alpha \rightarrow 0\) before touching the aspect ratio is that one can derive equations that include the nonhydrostatic effect due to the vertical velocity. The corresponding equations are the so-called Great Lake equations, first derived in Camassa et al. (1996, 1997). Taking \(\sigma \rightarrow 0\) after the rigid lid limit leads to the Lake equations. If one takes \(\sigma \ll \alpha \ll 1\), the result is the same, but the route is slightly different. Upon sending \(\sigma \rightarrow 0\), the vertical component in the Lagrangian (2.3) vanishes and upon assuming columnar motion, one can integrate the Lagrangian vertically. This leads to the Lagrangian for rotating shallow water. Sending \(\alpha \rightarrow 0\) corresponds to putting a rigid lid on top of the rotating shallow water equations and this leads to the Lake equations. Upon taking \(\alpha \rightarrow 0\) while keeping \(\sigma \) fixed, the equations (2.5) do not change, but the boundary conditions in (2.4) do:

$$\begin{aligned} \begin{matrix} p = \zeta \qquad &{} \text { at } z = 0,\\ w = 0 \qquad &{} \text { at } z = 0,\\ wdt + {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\circ dW_t^i = -({\mathsf {d}}{\varvec{\chi }}_t\cdot \nabla )h \qquad &{} \text { at } z = -h({\mathbf {x}}),\\ {\mathsf {d}}{\varvec{\chi }}_t\cdot {\mathbf {n}} = 0 \qquad &{} \text { on lateral boundaries}. \end{matrix} \end{aligned}$$
(2.6)

In the limit \(\alpha \rightarrow 0\), the depth \(\eta = h\), as the contribution of the free surface vanishes. Also, the expressions for the vertical velocity and the vertical component of the data vector field simplify, as the free surface contribution vanishes, and take the form

$$\begin{aligned} \begin{aligned} w&= \nabla \cdot \int _z^0 {\mathbf {u}}\, dz',\\ {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}&= \nabla \cdot \int _z^0 {\varvec{\xi }}_i\, dz'. \end{aligned} \end{aligned}$$
(2.7)

Averaging with the Newtonian approach leads to the following vertically averaged versions of the equations (2.5),

$$\begin{aligned} \begin{aligned} {\mathsf {d}} \overline{{\mathbf {u}}} + \frac{1}{h}\nabla \cdot (h\overline{{\mathsf {d}}{\varvec{\chi }}_{t}\otimes {\mathbf {u}}}) + (\nabla \overline{\varvec{\xi }}_{i})\cdot \overline{{\mathbf {u}}}\circ dW_t^i&= -\frac{1}{\mathrm{Fr}^2}\overline{\nabla {\mathsf {d}}p} - \frac{1}{\mathrm{Ro}} f {\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}}_t\\&\quad - \frac{1}{\mathrm{Ro}}\nabla (\overline{{\varvec{\xi }}}_i\cdot {\mathbf {R}})\circ dW_t^i,\\ {\mathsf {d}}\overline{b} + \nabla \cdot (\overline{b {\mathsf {d}}{\varvec{\chi }}}_t)&= 0, \end{aligned} \end{aligned}$$
(2.8)

and

$$\begin{aligned} \nabla \cdot (h\overline{{\mathsf {d}}{\varvec{\chi }}}_t) = 0. \end{aligned}$$
(2.9)

The continuity equation has become a weighted incompressibility condition (2.9), where the weight is determined by the bathymetry profile. As in the discussion above about the incompressibility condition (1.10), the weighted incompressibility must hold for the velocity field and the \({\varvec{\xi }}_i\) independently. If the bathymetry is flat, one finds the two-dimensional incompressibility condition. However, the momentum equation and the buoyancy equation above still suffer from the problem that terms are present which we, as yet, have not determined.

3.1 Leading Order Expansion in the Long Time: Very Small Wave Scaling Regime

As an initial approach, let us assume a leading order expansion in \(\sigma ^2\). Even though the Rossby number is small as well, we will consider a single scale expansion in \(\sigma ^2\) for the variables:

$$\begin{aligned} \begin{matrix} {\mathbf {u}} = {\mathbf {u}}_0 + o(1), &{} w = w_0 + o(1), &{} {\varvec{\xi }}_{3i} = {\varvec{\xi }}_{0,3i} + o(1),\\ {\mathsf {d}}{\varvec{\chi }}_{3t} = {\mathsf {d}}{\varvec{\chi }}_{0,3t} + o(1), &{} {\mathsf {d}}p = {\mathsf {d}}p_0 + o(1), &{} \zeta = \zeta _0 + o(1), \\ b = 0 + o(1). &{} &{} \end{matrix} \end{aligned}$$
(2.10)

The buoyancy does not contribute in the leading order expansion, since the stratification parameter is required to satisfy \(\mathfrak {s} \ll 1\) for the Boussinesq approximation. Also note that the data vector fields \({\varvec{\xi }}_{3i}\) are expanded in the same way as the velocity itself. No assumptions are made about the size of the data vector fields. Upon substituting (2.10) into (2.5), the vertical velocity equation at leading order implies hydrostatic balance

$$\begin{aligned} \frac{\partial }{\partial z}{\mathsf {d}}p_0 + 1\,dt= 0, \end{aligned}$$
(2.11)

and the dynamic boundary condition (1.8) implies that the leading order pressure is equal to the leading order free surface elevation.

Remark 2.1

Note that there is no stochasticity entering (2.11) explicitly. Due to the assumption of the pressure being a semimartingale, the pressure has the standard semimartingale decomposition. When there is no stochasticity in the equation, the martingale part of the pressure must vanish and we have the expression \({\mathsf {d}}p_0 = p_0 dt\) with a slight abuse of notation.

Interestingly, the substitution of the leading order expansion leads to a closed model even before averaging, when one uses the expression (1.42) for the vertical velocity as an additional equation. Given the boundary conditions in (2.6), the leading order expansion leads to a set of equations reminiscent of the Benney long wave model (Benney 1973). There are a few twists, though, since stochasticity and rotation are also involved. Moreover, the wave amplitude \(\alpha \) is very small, which enforces the rigid lid approximation in the vertical integral. At leading order, there cannot be any confusion as to which order of the expansion we are considering. Consequently, we may drop the subscript o in writing the following set of equations,

$$\begin{aligned} \begin{aligned} {\mathsf {d}}{\mathbf {u}} + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3){\mathbf {u}} + (\nabla \varvec{\xi }_{i})\cdot {\mathbf {u}}\circ dW_t^i&= -\frac{1}{\mathrm{Fr}^2}\nabla {\mathsf {d}}p - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times {\mathsf {d}}{\varvec{\chi }}_t - \frac{1}{\mathrm{Ro}}\nabla ({\varvec{\xi }}_i \cdot {\mathbf {R}})\circ dW_t^i,\\ w&= \nabla \cdot \int _z^0 {\mathbf {u}} \,dz',\\ {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}&= \nabla \cdot \int _z^0 {\varvec{\xi }}_i\, dz'. \end{aligned} \end{aligned}$$
(2.12)

Together with the weighted incompressibility condition in (2.9), the dynamic boundary condition on the pressure (1.8) and the lateral boundary condition (1.7), the Benney-like equations (2.12) form a closed set. The Benney long wave equations are interesting because they have a very rich mathematical structure, including an infinite hierarchy of conservation laws, as shown in Kupershmidt (2006). If we now make the additional assumption that the leading order component of the horizontal velocity field and the leading order component of the horizontal data vector field are independent of the vertical coordinate; that is, if we assume that the leading order component is columnar, then a considerable simplification of (2.12) occurs. Namely, the derivative in the vertical direction drops out. Consequently, it is no longer necessary to determine the vertical velocity and now every term in the equation is horizontal. This set of equations we will refer to as the stochastic, rotating, Lake equations, given by

$$\begin{aligned}&{\mathsf {d}}{\mathbf {u}} + ({\mathsf {d}}{\varvec{\chi }}_t\cdot \nabla ){\mathbf {u}}+ (\nabla \varvec{\xi }_i)\cdot {\mathbf {u}}\circ dW_t^i = -\frac{1}{\mathrm{Fr}^2}\nabla {\mathsf {d}}\zeta \nonumber \\&- \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times {\mathsf {d}}{\varvec{\chi }}_t - \frac{1}{\mathrm{Ro}}\nabla ({\varvec{\xi }}_i \cdot {\mathbf {R}})\circ dW_t^i, \end{aligned}$$
(2.13)

accompanied by the weighted incompressibility condition in (2.9) and the lateral boundary condition (1.7). The dynamic boundary condition can now be used to determine the pressure at the free surface. The deterministic, irrotational version of these equations has been shown by Levermore et al. (1996a, 1996b); Levermore and Oliver (1997) to be globally wellposed. These equations satisfy a Kelvin circulation theorem, namely

$$\begin{aligned} {\mathsf {d}}\oint _{c({{\mathsf {d}}{\varvec{\chi }}_t})}\left( {\mathbf {u}}+\frac{1}{\mathrm{Ro}}{\mathbf {R}}\right) \cdot d{\mathbf {x}} = 0, \end{aligned}$$
(2.14)

where \(c({{\mathsf {d}}{\varvec{\chi }}_t})\) is any fluid loop that is advected by the stochastic vector field \({{\mathsf {d}}{\varvec{\chi }}_t}\). This means that circulation is conserved, as there are no terms on the right hand side to generate circulation. Hence, the enstrophy in this model is conserved as well. The proof of the Kelvin circulation theorem is either a direct computation, or a corollary of the Euler–Poincaré theorem. We will derive these equations from a variational point of view as well, which will prove the Kelvin circulation theorem above.

3.2 Higher-Order Expansion in the Long Time: Very Small Wave Scaling Regime

Let us now consider a higher-order perturbation expansion:

$$\begin{aligned} \begin{matrix} {\mathbf {u}} = {\mathbf {u}}_0 + \sigma ^2{\mathbf {u}}_1 + o(\sigma ^2), &{} w = w_0 +\sigma ^2 w_1 + o(\sigma ^2),&{} {\varvec{\xi }}_{3i} = {\varvec{\xi }}_{0,3i} + \sigma ^2 {\varvec{\xi }}_{1,3i}+ o(\sigma ^2),\\ {\mathsf {d}}{\varvec{\chi }}_{3t} = {\mathsf {d}}{\varvec{\chi }}_{0,3t} + \sigma ^2 {\mathsf {d}}{\varvec{\chi }}_{1,3t}+ o(\sigma ^2), &{} {\mathsf {d}}p = {\mathsf {d}}p_0 + \sigma ^2 {\mathsf {d}}p_1 + o(\sigma ^2), &{} \zeta = \zeta _0 + \sigma ^2 \zeta _1+ o(\sigma ^2),\\ \qquad \mathfrak {s}b = \sigma ^2 b + o(\sigma ^2). &{} &{} \end{matrix} \nonumber \\ \end{aligned}$$
(2.15)

It is natural to assume that the leading order terms satisfy the Lake equations (2.13). Note that the stratification parameter is assumed to satisfy \(\mathcal {O}(\mathfrak {s})=\mathcal {O}(\sigma ^2)\). This will allow us to consider the buoyancy independently from the higher order terms that will appear in the equations to come. Hence, at leading order in the vertical velocity equation, we have hydrostatic balance (2.11) and in the horizontal component we have columnar motion. At the next order, we substitute (1.42) for the vertical velocity and obtain

$$\begin{aligned} \frac{\partial }{\partial z}{\mathsf {d}}p_1 + b_1\,dt = z\big ({\mathsf {d}}\nabla \cdot {\mathbf {u}}_0 + ({\mathsf {d}}{\varvec{\chi }}_{0,t}\cdot \nabla )(\nabla \cdot {\mathbf {u}}_0) - (\nabla \cdot {\mathsf {d}}{\varvec{\chi }}_{0,t})(\nabla \cdot {\mathbf {u}}_0)\big ).\qquad \end{aligned}$$
(2.16)

On the right hand side, everything in the brackets is independent of the vertical coordinate, so integration is particularly simple and leads to

$$\begin{aligned} {\mathsf {d}}p_1 = {\mathsf {d}}\zeta _1 - b_1 z dt + \frac{1}{2}z^2 \big ({\mathsf {d}}\nabla \cdot {\mathbf {u}}_0 + ({\mathsf {d}}{\varvec{\chi }}_{0,t}\cdot \nabla )(\nabla \cdot {\mathbf {u}}_0) - (\nabla \cdot {\mathsf {d}}{\varvec{\chi }}_{0,t})(\nabla \cdot {\mathbf {u}}_0)\big ).\nonumber \\ \end{aligned}$$
(2.17)

This shows that the pressure deviates from hydrostatic balance at order \(\sigma ^2\), as the pressure is a function of free surface elevation, buoyancy and horizontal velocity. The vertical average of the horizontal gradient of the pressure above is

$$\begin{aligned} \begin{aligned} \overline{{\mathsf {d}}\nabla p_1}&= \nabla {\mathsf {d}}\zeta _1 + \frac{1}{2}h\nabla b_1 dt + \frac{1}{6}h^2\big ({\mathsf {d}}\nabla \nabla \cdot {\mathbf {u}}_0 + ({\mathsf {d}}\varvec{\chi }_{0,t} \cdot \nabla )(\nabla \nabla \cdot {\mathbf {u}}_0) + (\nabla {\mathsf {d}}\varvec{\chi }_{0,t})\cdot (\nabla \nabla \cdot {\mathbf {u}}_0)\\&\quad - (\nabla \nabla \cdot {\mathsf {d}}{\varvec{\chi }}_{0,t})(\nabla \cdot {\mathbf {u}}_0)- (\nabla \cdot {\mathsf {d}}{\varvec{\chi }}_{0,t})(\nabla \nabla \cdot {\mathbf {u}}_0)\big ). \end{aligned}\nonumber \\ \end{aligned}$$
(2.18)

By using the weighted incompressibility condition (2.9), the expression above can be simplified and combined into

$$\begin{aligned} \overline{{\mathsf {d}}\nabla p_1} = \nabla {\mathsf {d}}\zeta _1 + \frac{1}{2}h\overline{\nabla b_1} dt + \bigg ({\mathsf {d}} + ({\mathsf {d}}{\varvec{\chi }}_{0,t}\cdot \nabla ) + (\nabla {\mathsf {d}}{\varvec{\chi }}_{0,t})\cdot \bigg )\left( \frac{1}{6}h^2(\nabla \nabla \cdot {\mathbf {u}}_0)\right) . \nonumber \\ \end{aligned}$$
(2.19)

The following observation allows us to deal with the average of the nonlinear term. Namely, if the leading order terms satisfy the stochastic, rotating Lake equations (2.13), then the leading order component of the stochastic velocity field is independent of the vertical coordinate. The higher order component of the stochastic vector field is not independent of the vertical coordinate, though, so its average is not trivial. Hence, the average of the full stochastic velocity field is

$$\begin{aligned} \overline{{\mathsf {d}}{\varvec{\chi }}}_t = {\mathsf {d}}{\varvec{\chi }}_{0,t} + \sigma ^2 \overline{{\mathsf {d}}{\varvec{\chi }}}_{1,t} + o(\sigma ^2). \end{aligned}$$
(2.20)

From this expression, it is clear that the average of the product minus the product of the average is a higher-order term:

$$\begin{aligned} \overline{{\mathsf {d}}{\varvec{\chi }}_t\otimes {\mathbf {u}}} - \overline{{\mathsf {d}}{\varvec{\chi }}}_t\otimes \overline{{\mathbf {u}}} = \mathcal O(\sigma ^4). \end{aligned}$$
(2.21)

Therefore, by adding and subtracting the product of the average in (2.8), we can write a closed system of equations. For notational convenience, we define

$$\begin{aligned} \overline{{\mathbf {V}}}({\mathbf {x}},t) := \overline{{\mathbf {u}}}({\mathbf {x}},t) + \frac{\sigma ^2}{6}h^2\nabla (\nabla \cdot \overline{{\mathbf {u}}}) + o(\sigma ^2), \end{aligned}$$
(2.22)

and use our expression for the average of the pressure (2.19) into (2.8) to write

$$\begin{aligned} \begin{aligned} {\mathsf {d}}\overline{{\mathbf {V}}} + (\overline{{\mathsf {d}}{\varvec{\chi }}_t}\cdot \nabla )\overline{{\mathbf {V}}} + (\nabla \overline{{\mathsf {d}}{\varvec{\chi }}_t})\cdot \overline{{\mathbf {V}}}&= -\frac{1}{\mathrm{Fr}^2}\nabla {\mathsf {d}}\zeta + \frac{1}{2}|\overline{{\mathbf {u}}}|^2dt - \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}h\nabla \overline{b} dt - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}\varvec{\chi }}_t \\&\quad - \frac{1}{\mathrm{Ro}}\nabla (\overline{{\varvec{\xi }}}_i \cdot {\mathbf {R}})\circ dW_t^i,\\ {\mathsf {d}}\overline{b} + (\overline{{\mathsf {d}}{\varvec{\chi }}}_t \cdot \nabla )\overline{b}&= 0. \end{aligned} \end{aligned}$$
(2.23)

The stratification parameter \(\mathfrak {s}\) is assumed to be of the same order as \(\sigma ^2\), the aspect ratio squared. When the buoyancy becomes negligible, the stratification parameter tends to zero. This removes the buoyancy from Eq. (2.23), but the nonhydrostatic pressure terms stay. Taking the shallow water limit by letting the aspect ratio tend to zero also removes the buoyancy contribution, since the buoyancy is linked to the aspect ratio in the expansions introduced in (2.15). Together with the weighted incompressibility condition (2.9) and lateral boundary condition (1.7), the set of equations (2.23) comprises the stochastic, rotating, thermal Great Lake equations. The deterministic, non-rotating version of these equations is presented in Camassa et al. (1996, 1997), together with the elliptic operator that relates \(\overline{{\mathbf {V}}}\) and \(\overline{{\mathbf {u}}}\). To solve for the pressure \({\mathsf {d}}\zeta \), one uses the elliptic operator just mentioned, which is defined by

$$\begin{aligned} \begin{aligned} h\overline{{\mathbf {V}}}&= h\overline{{\mathbf {u}}} + \left[ -\frac{\sigma ^2}{3}\nabla (h^3\nabla \cdot \overline{{\mathbf {u}}}) - \frac{\sigma ^2}{2}\nabla (h^2\overline{{\mathbf {u}}}\cdot \nabla h) + \frac{\sigma ^2}{2}h^2(\nabla \cdot \overline{{\mathbf {u}}})\nabla h + \sigma ^2 h(\overline{{\mathbf {u}}}\cdot \nabla h)\nabla h\right] ,\\&=: \mathfrak L(h)\overline{{\mathbf {u}}}. \end{aligned}\nonumber \\ \end{aligned}$$
(2.24)

This operator is positive-definite and self-adjoint since \(h>0\). An application of the Lax–Milgram theorem guarantees the continuous dependence of \(\overline{{\mathbf {u}}}\) on \(\overline{{\mathbf {V}}}\) (Levermore et al. 1996a, b). By operating with \(\nabla \cdot h\mathfrak L(h)^{-1} h\) on the velocity equation in (2.23) and using the weighted incompressibility condition (2.9), one finds an elliptic problem for \({\mathsf {d}}\zeta \). The Kelvin circulation theorem for the stochastic, rotating, thermal Great Lake equations is given by

$$\begin{aligned} {\mathsf {d}}\oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})}\left( \overline{{\mathbf {V}}}+\frac{1}{\mathrm{Ro}}{\mathbf {R}}\right) \cdot d{\mathbf {x}} = -\frac{\mathfrak {s}}{2}\oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})}h\nabla \overline{b} \cdot d{\mathbf {x}}dt. \end{aligned}$$
(2.25)

Here \(c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})\) is any fluid loop that is being advected by the vertically averaged stochastic vector field \(\overline{{\mathsf {d}}{\varvec{\chi }}_t}\). The right hand side of the circulation theorem reveals that circulation will be generated when the gradients of the buoyancy and the bathymetry are not aligned. This term can be seen as a baroclinic torque. The proof that the rotating, thermal, Great Lake equations satisfy this Kelvin theorem is postponed to end of the next subsection, where we will derive the same set of equations from a variational principle.

Remark 2.2

Note that the small aspect ratio limit \(\sigma \rightarrow 0\) reduces the Great Lake equations in (2.23) to the Lake equations in (2.13). If the bathymetry is flat, then the weighted incompressibility condition in (2.9) reduces to the usual two-dimensional incompressibility condition. In this case, the nonhydrostatic pressure term that is part of \({\mathbf {V}}\) vanishes and one obtains the two-dimensional version of the stochastic, rotating, Euler equations.

3.3 Averaged Euler–Poincaré Lagrangian for Long Time: Very Small Wave Scaling

To apply vertical averaging in the Euler–Poincaré setting, we return to the dimensionless Lagrangian (2.3) with boundary conditions given in (2.4). In line with the derivation of the Great Lake equations from the Newtonian point of view above, we assume that the horizontal velocity is independent of the vertical coordinate. This can be guaranteed upon replacing the horizontal velocity by its vertical average. In that situation, the expression for the vertical velocity in terms of the horizontal velocity in (1.42) can be integrated explicitly and we obtain as before

$$\begin{aligned} w = -\nabla \cdot (z+h)\overline{{\mathbf {u}}} = \nabla \cdot (\alpha \zeta - z)\overline{{\mathbf {u}}}. \end{aligned}$$
(2.26)

The same reasoning applies to the data vector fields, for which we obtain

$$\begin{aligned} {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i} = -\nabla \cdot (z+h)\overline{{\varvec{\xi }}}_i = \nabla \cdot (\alpha \zeta - z)\overline{{\varvec{\xi }}}_i. \end{aligned}$$
(2.27)

Note that in the limit \(\alpha \rightarrow 0\), the expression on the right hand side in (2.26) and (2.27) implies the free surface boundary condition when w and \({\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\) are evaluated on the free surface. However, evaluation on the bottom boundary does not imply the boundary condition (1.6) unless the weighted incompressibility condition (2.9) holds. Substituting (2.26) into the Euler–Boussinesq Lagrangian (2.3) and replacing \({\mathbf {u}}\) by \(\overline{{\mathbf {u}}}\) means that we can evaluate the vertical integral. Hence, we have the approximate Euler–Boussinesq Lagrangian

$$\begin{aligned} \begin{aligned} \ell _{EB}&\approx \int _{CS}\int _{-h}^{\alpha \zeta } \left( \frac{1}{2}|\overline{{\mathbf {u}}}|^2 + \frac{\sigma ^2}{2}\big ((z+h)(\nabla \cdot \overline{{\mathbf {u}}}) + (\overline{{\mathbf {u}}}\cdot \nabla )h\big )^2 + \frac{1}{\mathrm{Ro}}\overline{{\mathbf {u}}}\cdot {\mathbf {R}} -\frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}\overline{b})z\right) dz\,dx\,dy. \end{aligned}\nonumber \\ \end{aligned}$$
(2.28)

This Lagrangian is an integral over the horizontal cross section of the domain \(\Omega \), which we call CS. Evaluating the vertical integral leads to

$$\begin{aligned} \begin{aligned} \ell _{TRGN}&= \int _{CS}\left( \frac{1}{2}|\overline{{\mathbf {u}}}|^2 + \frac{\sigma ^2}{6}\eta ^2(\nabla \cdot \overline{{\mathbf {u}}})^2 + \frac{\sigma ^2}{2}\eta (\nabla \cdot \overline{{\mathbf {u}}})(\overline{{\mathbf {u}}}\cdot \nabla h)+ \frac{\sigma ^2}{2}(\overline{{\mathbf {u}}}\cdot \nabla h)^2 + \frac{1}{\mathrm{Ro}}\overline{{\mathbf {u}}}\cdot {\mathbf {R}}\right. \\&\quad \left. - \frac{1}{2\,\mathrm{Fr}^2}(1+\mathfrak {s}\overline{b})(\eta -2h)\right) \eta \,dx\,dy. \end{aligned} \end{aligned}$$
(2.29)

The subscript on the Lagrangian in (2.29) stands for thermal rotating Green–Naghdi, because the equations that this Lagrangian gives rise to are a thermal and rotating extension to the usual Green–Naghdi equations (Green and Naghdi 1976). The incompressibility constraint has been used to ensure that the expression for the vertical velocity is valid and are thus no longer required. However, the weighted incompressibility condition (2.9) must still hold; so, we introduce a new constraint to make the total depth equal to the bathymetry. Weighted incompressibility has to be enforced via a constraint because it affects the measure \(\eta \,dx\,dy\) in the Lagrangian above. The constraint is equivalent to saying that the free surface elevation is zero, that is, \(\eta -h=0\). Thus, the action for the thermal rotating Great Lake equations is given by

$$\begin{aligned} S_{TRGL} = \int _{t_1}^{t_2}\ell _{TRGN}\,dt + \left\langle \frac{1}{\mathrm{Fr}^2}{\mathsf {d}}\pi ,\eta -h\right\rangle =: \int _{t_1}^{t_2}c\ell _{TRGL}. \end{aligned}$$
(2.30)

The action in (2.30) has been suggestively called the thermal rotating Great Lake action and defines the constrained thermal rotating Great Lake Lagrangian. Note that this Lagrangian features the \(H_{div}\) Sobolev norm in the situation where the bathymetry is flat, which has interesting relations with integrable systems and geometric statistics, as shown in Khesin et al. (2013). When the bathymetry is nontrivial, the norm is more complicated. Here \({\mathsf {d}}\pi \) is a semimartingale Lagrange multiplier, whose purpose is to ensure that the weighted incompressibility condition holds. In order to apply the Euler–Poincaré theorem 1.1 to this Lagrangian, we need to define the variations. By substituting the higher order perturbation expansion (2.15) into the formulas for the variations in the theorem, we obtain

$$\begin{aligned} \begin{aligned} \delta \overline{{\mathbf {u}}} \,dt&= {\mathsf {d}}{\mathbf {v}} - \big [{\mathsf {d}}\overline{{\varvec{\chi }}_t}, {\mathbf {v}}\big ], \end{aligned} \end{aligned}$$
(2.31)

where the arbitrary vector field \({\mathbf {v}}\) is a vector field semimartingale. The variations of the advected quantities are obtained by directly integrating the formulae for the variations in the three-dimensional case. First we notice that the only advected quantities in this problem are scalar functions and volume forms, which due to incompressibility, satisfy the same form of advection equation, as we saw above in the Newtonian averaging principle. The functional derivative and spatial derivatives commute. Hence, if \({\mathbf {u}}_3\) is incompressible, then \(\delta {\mathbf {u}}_3\) must be incompressible, as well. This argument implies that the arbitrary vector field is also incompressible, which means that the constraints for the variations of the buoyancy and the density can be shown to satisfy

$$\begin{aligned} \begin{aligned} \delta \overline{b}\, dt&= - ({\mathbf {v}}\cdot \nabla )\overline{b}\,dt,\\ \delta \int _{-h}^{\alpha \zeta } D dz\, dt&= -\nabla \cdot \left( \int _{-h}^{\alpha \zeta } D dz\,{\mathbf {v}}\right) dt. \end{aligned} \end{aligned}$$
(2.32)

In this paper, \(D = 1\), so the vertical integral of D is the depth \(\eta = \alpha \zeta + h\), showing that the depth \(\eta \) functions as a two-dimensional density; hence, its variation satisfies

$$\begin{aligned} \delta \eta \,dt = -\nabla \cdot \left( \eta {\mathbf {v}}\right) dt. \end{aligned}$$
(2.33)

In the \(\alpha \rightarrow 0\) limit, the depth is given by the bathymetry \(\eta = h\), which is the constraint introduced to imply weighted divergence. The variations of the thermal rotating Great Lake Lagrangian in (2.30) are

$$\begin{aligned} \begin{aligned} \frac{\delta c\ell _{TRGL}}{\delta \overline{{\mathbf {u}}}}&= \eta \overline{{\mathbf {u}}} - \frac{\sigma ^2}{3}\nabla (\eta ^3\nabla \cdot \overline{{\mathbf {u}}}) - \frac{\sigma ^2}{2}\nabla \big (\eta ^2(\overline{{\mathbf {u}}}\cdot \nabla h)\big ) + \frac{\sigma ^2}{2}\eta ^2(\nabla \cdot \overline{{\mathbf {u}}})\nabla h + \sigma ^2\eta (\overline{{\mathbf {u}}}\cdot \nabla h)\nabla h + \frac{1}{\mathrm{Ro}}\eta {\mathbf {R}}\,,\\ \frac{\delta c\ell _{TRGL}}{\delta \eta }&= \frac{1}{2}|\overline{{\mathbf {u}}}|^2 + \frac{\sigma ^2}{2}(\eta \nabla \cdot \overline{{\mathbf {u}}}+\overline{{\mathbf {u}}}\cdot \nabla h)^2 + \frac{1}{\mathrm{Ro}}(\overline{{\mathbf {u}}}\cdot {\mathbf {R}}) - \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}\overline{b})(\eta -h) - \frac{1}{\mathrm{Fr}^2}{\mathsf {d}}\pi , \\ \frac{\delta c\ell _{TRGL}}{\delta \overline{b}}&= -\frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\eta ^2-2\eta h) \,, \\ \frac{\delta c\ell _{TRGL}}{\delta {\mathsf {d}}\pi }&= \frac{1}{\mathrm{Fr}^2}(\eta - h)\,. \end{aligned} \end{aligned}$$
(2.34)

The variational derivative with respect to \(\overline{{\mathbf {u}}}\) of the thermal rotating Great Lake Lagrangian shows that the elliptic operator that relates \(\overline{{\mathbf {V}}}\) and \(\overline{{\mathbf {u}}}\) that we encountered in (2.24) arises naturally in the variational context upon evaluating \(\eta =h\). Note that the variational derivative with respect to \(\eta \) simplifies considerably upon evaluating \(\eta =h\). The second term in the variational derivative with respect to \(\eta \) is the square of the weighted incompressibility condition, which is equal to zero. Hence, this term disappears. The fourth term in the variational derivative with respect to \(\eta \) vanishes since \(\eta -h=0\). The variational derivative with respect to the buoyancy simplifies. An application of the stochastic Euler–Poincaré theorem 1.1 to the Great Lake Lagrangian in (2.30) with these variational derivatives and the variations in (2.31) leads to the stochastic Great Lake equations (2.23), with rotation and buoyancy. For notational convenience, let us use (2.24) to define \(\overline{{\mathbf {V}}}\). Then, the thermal rotating Great Lake equations are given by

$$\begin{aligned} \begin{aligned} {\mathsf {d}}\overline{{\mathbf {V}}} + (\overline{{\mathsf {d}}{\varvec{\chi }}_t}\cdot \nabla )\overline{{\mathbf {V}}} + (\nabla \overline{{\mathsf {d}}{\varvec{\chi }}_t})\cdot \overline{{\mathbf {V}}}&= -\frac{1}{\mathrm{Fr}^2}\nabla {\mathsf {d}}\pi + \frac{1}{2}|\overline{{\mathbf {u}}}|^2 dt - \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}h\nabla \overline{b} dt - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}_t} - \frac{1}{\mathrm{Ro}}\nabla (\overline{{\varvec{\xi }}}_i\cdot {\mathbf {R}})\circ dW_t^i,\\ {\mathsf {d}}\overline{b} + (\overline{{\mathsf {d}}{\varvec{\chi }}_{t}} \cdot \nabla )\overline{b}&= 0,\\ \nabla \cdot (h\overline{{\mathsf {d}}{\varvec{\chi }}_t})&= 0, \end{aligned} \end{aligned}$$
(2.35)

and with the boundary condition

$$\begin{aligned} \overline{{\mathsf {d}}{\varvec{\chi }}}_t \cdot {\mathbf {n}} = 0. \end{aligned}$$
(2.36)

The pressure \({\mathsf {d}}\pi \) is solved for using the elliptic operator defined in (2.24). Hence, one can make the identification \(\pi =\zeta \). This calculation shows that the Great Lake equations with rotation, stratification and stochasticity can be obtained by averaging the equations and using a perturbation series approach, or by taking a variational approach. The results are identically equal. Since the Lagrangian framework implies the Kelvin circulation theorem (2.25), the proof is now immediate that the circulation theorem has the form

$$\begin{aligned} \begin{aligned} {\mathsf {d}}\oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})}\left( \overline{{\mathbf {V}}}+\frac{1}{\mathrm{Ro}}{\mathbf {R}}\right) \cdot d{\mathbf {x}}&= -\frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}\oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})}h\nabla \overline{b}\cdot d{\mathbf {x}} \,dt,\\&= -\frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}\int \!\!\int _{\partial S = c(\overline{{\mathsf {d}}\varvec{\chi }_t})}\nabla h \times \nabla \overline{b} \,d{\mathbf {S}}\,dt. \end{aligned} \end{aligned}$$
(2.37)

Thus, in this scaling regime, applying asymptotics to the equations implies the same result as applying the asymptotics in the variational principle.

Remark 2.3

(Kelvin theorem result for generation of horizontal circulation). The Kelvin circulation theorem in (2.37) shows that any misalignment of the horizontal gradients of the bathymetry and of the vertically averaged buoyancy will generate horizontal circulation in the material loop \(c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})\) which follows the stochastic Lagrangian flow velocity \(\overline{{\mathsf {d}}{\varvec{\chi }}_t}\) in the horizontal plane given in Eq. (1.5). The Kelvin circulation theorem (2.37) implies an evolution equation for potential vorticity, as well.

In the next section, we will extend the comparative asymptotic expansion approach to consider the short time-small wave limit. This extension will be accomplished by first deriving equations using asymptotics in the Euler–Boussinesq equations and later doing asymptotics in the Lagrangian and applying the Euler–Poincaré theorem.

4 Short Time—Small Wave Scaling Regime

Short time corresponds to choosing the time scale to be \(T = L/\sqrt{gH}\), the time it takes for a gravity wave to traverse the horizontal length scale. ‘Small wave’ means that the amplitude of the wave is small, but not small enough to consider taking the rigid lid limit. In this setting, the scales are given by

figure c

In this scaling regime, the EB Lagrangian takes the form

$$\begin{aligned} \ell _{EB}({\mathbf {u}}_3,b,D) = \int _{\Omega }D\left( \frac{1}{2}|{\mathbf {u}}|^2 + \frac{\sigma ^2}{2} w^2 + \frac{1}{\mathrm{Ro}}{\mathbf {u}}\cdot {\mathbf {R}} - \frac{1}{\mathrm{Fr}^2}(1+b)z\right) \,dx\,dy\,dz, \nonumber \\ \end{aligned}$$
(3.2)

so the corresponding action is given by

$$\begin{aligned} S_{EB} = \int _{t_1}^{t_2}\ell _{EB} \,dt - \left\langle \frac{1}{\mathrm{Fr}^2}{\mathsf {d}}p, D-1 \right\rangle =: \int _{t_1}^{t_2} c\ell _{EB}, \end{aligned}$$
(3.3)

with boundary conditions given by

$$\begin{aligned} \begin{matrix} p = \alpha \zeta \qquad &{} \text { at } z = \alpha \zeta ({\mathbf {x}},t),\\ wdt + {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\circ dW_t^i = \alpha \left( \dfrac{1}{\mathrm{Fr}}{\mathsf {d}}\zeta + ({\mathsf {d}}{\varvec{\chi }}_t\cdot \nabla )\zeta \right) \qquad &{} \text { at } z = \alpha \zeta ({\mathbf {x}},t),\\ wdt + {\hat{\mathbf {z}}}\cdot {\varvec{\xi }}_{3i}\circ dW_t^i = -({\mathsf {d}}{\varvec{\chi }}_t\cdot \nabla )h \qquad &{} \text { at } z = -h({\mathbf {x}}),\\ {\mathsf {d}}{\varvec{\chi }}_t\cdot {\mathbf {n}} = 0 \qquad &{} \text { on lateral boundaries}. \end{matrix} \end{aligned}$$
(3.4)

An application of the stochastic Euler–Poincaré theorem 1.1 on the short-time scaled Lagrangian in (3.2) yields the following equations

$$\begin{aligned} \frac{1}{\mathrm{Fr}}{\mathsf {d}}{\mathbf {u}} + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3){\mathbf {u}} + (\nabla {\varvec{\xi }}_{3i})\cdot {\mathbf {u}}_3\circ dW_t^i&= -\frac{1}{\mathrm{Fr}^2}\nabla {\mathsf {d}}p - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times {\mathsf {d}}{\varvec{\chi }}_t \nonumber \\&\quad - \frac{1}{\mathrm{Ro}}\nabla ({\varvec{\xi }}_i\cdot R)\circ dW_t^i,\nonumber \\ \sigma ^2\left( \frac{1}{\mathrm{Fr}} {\mathsf {d}}w + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)w + \Big (\frac{\partial }{\partial z}{\varvec{\xi }}_{3i}\Big )\cdot {\mathbf {u}}_3\circ dW_t^i \right)&= -\frac{1}{\mathrm{Fr}^2}\frac{\partial }{\partial z}{\mathsf {d}}p -\frac{1}{\mathrm{Fr}^2}(1 + b)dt,\nonumber \\ \frac{1}{\mathrm{Fr}}{\mathsf {d}}b + ({\mathsf {d}}{\varvec{\chi }}_{3t}\cdot \nabla _3)b&= 0,\nonumber \\ \nabla _3\cdot {\mathsf {d}}{\varvec{\chi }}_{3t}&= 0. \end{aligned}$$
(3.5)

These equations satisfy the Kelvin circulation theorem, which for the Euler–Boussinesq equations takes the form of (1.34), and also have conservation of potential vorticity along fluid trajectories, as in (1.36), as well as conservation of an infinity of integral quantities (1.37), but now the Strouhal number is explicitly given in terms of the Froude number. In this scaling, the free surface is small rather than very small. Hence, we will not take the limit of the Froude number going to zero explicitly. Instead, we will introduce a regular perturbation expansion with small parameters \(\epsilon \) and \(\gamma \) whose magnitudes need to be determined with respect to \(\alpha \), \(\mathrm{Fr}\) and \(\sigma \).

$$\begin{aligned} \begin{matrix} {\mathbf {u}} = {\mathbf {u}}_0 + \epsilon {\mathbf {u}}_1 + o(\epsilon ), &{} w = w_0 +\epsilon w_1 + o(\epsilon ), &{} {\varvec{\xi }}_i = {\varvec{\xi }}_{0,i} + \epsilon {\varvec{\xi }}_{1,i}+ o(\epsilon ),\\ {\mathsf {d}}{\varvec{\chi }}_t = {\mathsf {d}}{\varvec{\chi }}_{0,t} + \epsilon {\mathsf {d}}{\varvec{\chi }}_{1,t}+ o(\epsilon ), &{} p = p_0 + \gamma p_1 + \gamma ^2 p_2 + o(\gamma ^2), &{} b = \mathfrak {s} b_1 + \mathfrak {s}^2 b_2 + o(\gamma ^2). \end{matrix} \nonumber \\ \end{aligned}$$
(3.6)

Substitution of (3.6) into (3.5) provides equations of unknown order. By requiring certain balances to hold, the order of the dimensionless numbers can be related to each other. The boundary condition related to the vertical velocity at the free surface in (3.4) implies that \(\alpha = \mathcal O(\mathrm{Fr})\). In the horizontal velocity equation, the leading order velocity \(\mathrm{Fr}\,{\mathsf {d}}{\mathbf {u}}_0\) needs to be of the same order as \(\gamma \nabla {\mathsf {d}} p_1\), which means that \(\gamma = \mathcal O(\mathrm{Fr})\). At the next order, \(\mathrm{Fr}\, \epsilon \,{\mathsf {d}}{\mathbf {u}}_1\) is required to be of the same order as \(\gamma ^2\nabla {\mathsf {d}} p_2\), which implies that \(\epsilon = \mathcal O(\mathrm{Fr})\). In the vertical velocity equation, we want hydrostatic balance to be broken at \(\mathcal O(\gamma ^2)\), which means that \(\mathrm{Fr}\,\sigma ^2 {\mathsf {d}}w_0\) has to be of the same order as \(\gamma ^2\frac{\partial }{\partial z}{\mathsf {d}}p_2 \). It also implies for our ordering scheme that \(\sigma ^2 = \mathcal O(\mathrm{Fr})\). In the Boussinesq approximation, we assumed that \(\mathcal O(\mathfrak {s})=\mathcal O(\mathrm{Fr})\). To summarise, our ordering scheme is now fixed to be

$$\begin{aligned} \mathcal O(\alpha ) = \mathcal O(\mathfrak {s}) = \mathcal O(\gamma ) = \mathcal O(\epsilon ) = \mathcal O(\mathrm{Fr}) = \mathcal O(\sigma ^2). \end{aligned}$$
(3.7)

4.1 Averaging of Newton’s second Law in the Short Time: Small Wave Scaling

Averaging in the Newtonian equations leads to the following vertically averaged version of (3.5),

$$\begin{aligned} \begin{aligned} \frac{1}{\mathrm{Fr}}\,{\mathsf {d}} \overline{{\mathbf {u}}} + \frac{1}{\eta }\nabla \cdot (\eta \overline{{\mathsf {d}}{\varvec{\chi }}_t\otimes {\mathbf {u}}}) + (\nabla \varvec{\xi }_i)\cdot \overline{{\mathbf {u}}}\circ dW_t^i&= -\frac{1}{\mathrm{Fr}^2}\nabla \overline{{\mathsf {d}}p} - \frac{1}{\mathrm{Ro}} f {\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}_t} \\&\quad - \frac{1}{\mathrm{Ro}}\nabla ( \overline{\varvec{\xi }}_i\cdot {\mathbf {R}})\circ dW_t^i,\\ \frac{1}{\mathrm{Fr}}{\mathsf {d}}\overline{b} + \nabla \cdot (\overline{b {\mathsf {d}}{\varvec{\chi }}_t})&= 0,\\ \frac{1}{\mathrm{Fr}}{\mathsf {d}}\eta + \nabla \cdot (\eta \overline{{\mathsf {d}}\varvec{\chi }_t})&= 0, \end{aligned} \end{aligned}$$
(3.8)

where \(\overline{{\mathsf {d}}{\varvec{\chi }}_t} \) is the vertical average of \({\mathsf {d}}{\varvec{\chi }}_t \) in equation (3.6); namely,

$$\begin{aligned} \overline{{\mathsf {d}}{\varvec{\chi }}_t} := {\mathsf {d}}{\varvec{\chi }}_{0,t} + \epsilon \overline{{\mathsf {d}}{\varvec{\chi }}_{1,t}}+ o(\epsilon ). \end{aligned}$$
(3.9)

In this part of our discussion, we will not consider a leading order expansion before doing a higher order expansion. Instead, we work with directly with the expansion introduced in (3.6) and use the ordering scheme (3.7) to apply single scale asymptotics.

Remark 3.1

It is possible to study the system (3.8) on its own. One can simplify the system by dropping the Coriolis terms and assume that the flow is irrotational. The equations (3.8) can then be written in the so-called Zakharov–Craig–Sulem formulation. Alternatively, one can reformulate the system in terms of the free surface elevation and the horizontal discharge. Both of these approaches are explained in great detail in lecture notes by Lannes (2019). See also Lannes (2013) for a comprehensive and complete treatment of the general water wave problem and Lannes (2005) for the wellposedness results on the water wave problem in two and three dimensions.

At leading order in the vertical velocity equation one finds

$$\begin{aligned} \frac{\partial }{\partial z}{\mathsf {d}}p_0 + 1\,dt = 0, \end{aligned}$$
(3.10)

and from the horizontal velocity equation at the same order,

$$\begin{aligned} \nabla {\mathsf {d}}p_0 = 0, \end{aligned}$$
(3.11)

which implies hydrostatic balance. This information determines the leading order pressure, upon integrating in the vertical direction, to find

$$\begin{aligned} {\mathsf {d}}p_0 = (const. - z)dt, \end{aligned}$$
(3.12)

for the leading order pressure. In Remark 2.1 we discussed how to deal with the semimartingale equations when the stochasticity is absent. This allows us to compute the expression for \(p_0\) above. The arbitrary constant is due to integration and will be eliminated later using the boundary condition for the pressure. At the next order in the vertical velocity equation, one finds

$$\begin{aligned} \frac{\partial }{\partial z}{\mathsf {d}}p_1 + b_1\,dt = 0. \end{aligned}$$
(3.13)

Vertical integration of the expression above leads to

$$\begin{aligned} {\mathsf {d}}p_1 = \left( -\int ^z b_1 dz' + \psi ({\mathbf {x}},t)\right) dt, \end{aligned}$$
(3.14)

where \(\psi ({\mathbf {x}},t)\) is an arbitrary function of horizontal coordinates and time, introduced by the integration. From the horizontal velocity equation at the same order, we have

$$\begin{aligned} {\mathsf {d}}{\mathbf {u}}_0 = -\nabla {\mathsf {d}}p_1. \end{aligned}$$
(3.15)

By applying the gradient to (3.14) and taking the vertical derivative of (3.15), we can derive a relation between the horizontal velocity field and the buoyancy,

$$\begin{aligned} \frac{\partial }{\partial z}{\mathsf {d}}{\mathbf {u}}_0 = \nabla b_1 \,dt. \end{aligned}$$
(3.16)

From the buoyancy equation at order \(\mathcal O(\mathfrak {s})\), it is clear that \(b_1\) is independent of time. Upon integrating (3.16) both vertically and in time, one finds

$$\begin{aligned} {\mathbf {u}}_0({\mathbf {x}},z,t) = t\int ^z\nabla b_1({\mathbf {x}},z) dz' + {\mathbf {u}}_0'({\mathbf {x}},t) + \widetilde{{\mathbf {u}}_0}({\mathbf {x}},z). \end{aligned}$$
(3.17)

Unless \(\nabla b_1 = 0\), the first term in (3.17) grows linearly in time. Consequently, we choose the buoyancy \(b_1\) to have the following profile

$$\begin{aligned} b_1(z) = \widetilde{b} - Sz, \end{aligned}$$
(3.18)

where \(\widetilde{b}\) is some constant background buoyancy and S is some \(\mathcal O(1)\) positive constant. Of course, one can choose a more complicated and more realistic dependence on the vertical coordinate, at the cost of making some computations slightly more involved. The first term in (3.17) now vanishes. The third term in (3.17) arose due to integration with respect to time, hence \(\widetilde{{\mathbf {u}}_0}\) plays the role of the initial condition. It is also the only term that has z-dependence. So, let us choose an initial condition which is independent of the vertical coordinate. This choice leaves us with

$$\begin{aligned} {\mathbf {u}}_0({\mathbf {x}},t) = {\mathbf {u}}_0'({\mathbf {x}},t) + \widetilde{{\mathbf {u}}_0}({\mathbf {x}}). \end{aligned}$$
(3.19)

Hence, \({\mathbf {u}}_0\) has no vertical dependence. We can then use the incompressibility condition (1.9) to obtain an expression for the vertical velocity as in (1.42), but now only looking at the leading order component of this relation. This leads to

$$\begin{aligned} w_0 = -(z+h)\nabla \cdot {\mathbf {u}}_0, \end{aligned}$$
(3.20)

provided the variations of the bathymetry are small enough. Substituting the expression for the leading order vertical velocity into the vertical velocity equation at order \(\mathcal O(\gamma ^2)\) yields

$$\begin{aligned} -(z+h){\mathsf {d}}(\nabla \cdot {\mathbf {u}}_0) + \frac{\partial }{\partial z} {\mathsf {d}}p_2 + b_2\,dt = 0. \end{aligned}$$
(3.21)

From the equation above, we can determine an expression for \(p_2\). Rearranging and taking a vertical integral yield

$$\begin{aligned} {\mathsf {d}}p_2 = \left( \frac{1}{2}z^2+zh\right) {\mathsf {d}}(\nabla \cdot {\mathbf {u}}_0) - \int ^z b_2 dz'dt + \psi '({\mathbf {x}},t). \end{aligned}$$
(3.22)

Since the expressions for \({\mathsf {d}}p_1\) and \({\mathsf {d}}p_2\) in (3.14) and (3.22), respectively, involve the unknown functions \(\psi ({\mathbf {x}},t)\) and \(\psi '({\mathbf {x}},t)\), we are not yet in the position to write down the average of the pressure. By means of the dynamic boundary condition (1.8) and the expansion for the pressure in (3.6), though, we can write

$$\begin{aligned} 0&= [{\mathsf {d}}p_0 + \gamma {\mathsf {d}}p_1 + \gamma ^2 {\mathsf {d}}p_2 + \mathcal O(\gamma ^3)]|_{z=\alpha \zeta } dt\nonumber \\&= (const.\,dt - \alpha {\mathsf {d}}\zeta + \gamma \left( -\int ^{\alpha \zeta }b_1 dz' + \psi ({\mathbf {x}},t)\right) dt + \gamma ^2\Bigg [ \left( \frac{1}{2}\alpha ^2 \zeta ^2 + \alpha \zeta h\right) {\mathsf {d}}(\nabla \cdot {\mathbf {u}}_0)\nonumber \\&\quad - \left( \int ^{\alpha \zeta } b_2 dz' + \psi '({\mathbf {x}},t)\right) \Bigg ]dt + \mathcal O(\gamma ^3). \end{aligned}$$
(3.23)

The difference between the pressure at the free surface and elsewhere in the domain can now be evaluated. In particular, functions that are independent of z will be eliminated in this procedure and we are left with

$$\begin{aligned} \begin{aligned} {\mathsf {d}}p&= -z\,dt-\alpha {\mathsf {d}}\zeta + \gamma \int _z^{\alpha \zeta }b_1dz'dt\\&\quad + \gamma ^2\left[ \left( \frac{1}{2}(z^2 - \alpha ^2\gamma ^2)+(z-\alpha \gamma )h\right) {\mathsf {d}}(\nabla \cdot {\mathbf {u}}_0) + \int _z^{\alpha \zeta }b_2 dz'dt\right] + \mathcal O(\gamma ^3). \end{aligned} \end{aligned}$$
(3.24)

We can now determine the gradient of the pressure and collect terms that are of order \(\mathcal O(\gamma ^3)\) or equivalent in the remainder. Since \(b_1\) does not depend on the horizontal coordinates, the gradient of \(b_1\) vanishes and we have

$$\begin{aligned}&\nabla {\mathsf {d}}p = \alpha (1+\widetilde{b})\nabla {\mathsf {d}}\zeta + \gamma ^2\left[ \left( \frac{1}{2}z^2 + zh\right) {\mathsf {d}}\nabla (\nabla \cdot {\mathbf {u}}_0) + \int _z^0 \nabla b_2 dz'dt\right] \nonumber \\&+ \mathcal O(\gamma ^3,\alpha ^2\gamma ,\alpha \gamma ^2), \end{aligned}$$
(3.25)

where the contribution of \(\widetilde{b}\) is due to the evaluation of \(b_1\) at the free surface boundary. By taking the vertical average of the pressure gradient and switching the order of integration on the \(b_2\) term, we obtain

$$\begin{aligned} \overline{\nabla {\mathsf {d}}p} = \alpha (1+\widetilde{b})\nabla {\mathsf {d}}\zeta + \gamma ^2\left( \frac{1}{3}h^2 {\mathsf {d}}\nabla (\nabla \cdot {\mathbf {u}}_0) + \overline{(z+h)\nabla b_2}dt\right) + \mathcal O(\gamma ^3,\alpha ^2\gamma ,\alpha \gamma ^2).\nonumber \\ \end{aligned}$$
(3.26)

At this stage, we can make a choice. We can use the averaged equation for the advection of buoyancy (1.41), or we can use the expanded buoyancy equation and find an equation for the evolution \(\overline{(z+h)\nabla b_2}\). The latter choice dictates that we look at the expanded buoyancy equation at order \(\mathcal O(\gamma ^2)\), where we have

$$\begin{aligned} {\mathsf {d}}b_2 - S(z+h)(\nabla \cdot {\mathbf {u}}_0)dt = 0. \end{aligned}$$
(3.27)

Here we have used (3.18) and (3.20). By taking the gradient, then multiplying by \((z+h)\) and taking the average, we obtain after some algebra

$$\begin{aligned} {\mathsf {d}}\overline{(z+h)\nabla b_2} = S\left( \frac{1}{3}h^2\nabla (\nabla \cdot {\mathbf {u}}_0)\right) dt. \end{aligned}$$
(3.28)

Similar to the derivation of the Great Lake equations, the difference between the average of the nonlinearity and the product of the average is of higher order, since \({\mathbf {u}}_0\) is independent of the vertical coordinate. Therefore, we can also express \(\overline{{\mathbf {u}}} = {\mathbf {u}}_0 + \epsilon \overline{{\mathbf {u}}_1} + \mathcal O(\epsilon ^2)\). At this stage, one follows (Camassa and Holm 1992) to introduce the variables

$$\begin{aligned} \begin{aligned} {\mathbf {A}}:&= \overline{(z+h)\nabla b_2},\\ {\mathbf {D}}:&= \frac{1}{3}h^2\nabla (\nabla \cdot \overline{{\mathbf {u}}}), \end{aligned} \end{aligned}$$
(3.29)

and writes the following set of stochastic partial differential equations (SPDEs),

$$\begin{aligned} \begin{aligned} \frac{1}{\mathrm{Fr}}\,{\mathsf {d}}\overline{{\mathbf {u}}} + (\overline{{\mathsf {d}}{\varvec{\chi }}}_t\cdot \nabla )\overline{{\mathbf {u}}} + (\nabla \varvec{\xi }_i)\cdot \overline{{\mathbf {u}}}\circ dW_t^i&= -\frac{\alpha }{\mathrm{Fr}^2}(1+\widetilde{b})\nabla {\mathsf {d}}\zeta - {\mathbf {A}}dt + {\mathsf {d}}{\mathbf {D}}\\&\qquad - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}}_t - \frac{1}{\mathrm{Ro}}\nabla (\overline{{\varvec{\xi }}}_i\cdot {\mathbf {R}})\circ dW_t^i,\\ \frac{\alpha }{\mathrm{Fr}}{\mathsf {d}}\zeta + \nabla \cdot (\alpha \zeta + h)\overline{{\mathsf {d}}{\varvec{\chi }}}_t&= 0,\\ {\mathsf {d}}{\mathbf {A}}&= S{\mathbf {D}}dt. \end{aligned} \nonumber \\ \end{aligned}$$
(3.30)

where \(\overline{{\mathsf {d}}{\varvec{\chi }}_t} \) is defined in Eq. (3.9).

Equation (3.30) comprises the stochastic version of those obtained in Camassa and Holm (1992), provided one sets the dynamic boundary condition to \(p=\widetilde{p}\), rather than zero.

In the special case of deterministic, irrotational motion around the quiescent state \(\overline{{\mathbf {u}}}=0\), the covector quantities \({\mathbf {A}}\) and \({\mathbf {D}}\) form an oscillator pair which oscillates with the Brunt-Väisälä frequency S. Also, in the deterministic case, an elimination procedure allows one to derive the Kadomtsev–Petviashvili equation and subsequently the Korteweg-De Vries equation for shallow water waves, as is done in Camassa and Holm (1992). The direct approach for the derivations for water wave equations requires the substitution of the velocity field into the free surface equation, which requires time derivatives. In the stochastic case, however, one cannot take these time derivatives; so, the corresponding stochastic shallow water wave equations cannot be derived by using SALT. If instead, one takes a pathwise approach so that at least one time derivative can be taken, then the corresponding water-wave equations can be derived in this framework. In the next subsection, a hierarchy of stochastic water-wave equations is derived from the variational point of view.

The set of equations (3.30) can be solved by observing that the operator F, defined by \(F\overline{{\mathbf {u}}}:= \overline{{\mathbf {u}}} - \frac{\gamma ^2}{3}h^2\nabla (\nabla \cdot \overline{{\mathbf {u}}})\), is a positive definite, self-adjoint and invertible operator. The Kelvin circulation theorem takes the following form for the equations in (3.30),

$$\begin{aligned}&\frac{1}{\mathrm{Fr}}{\mathsf {d}}\oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})}\left( \overline{{\mathbf {u}}} - {\mathbf {D}} + \frac{1}{\mathrm{Ro}}{\mathbf {R}}\right) \cdot d{\mathbf {x}}\nonumber \\&= - \frac{1}{\mathrm{Fr}^2}\oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}}_t)} \left( (\overline{{\mathsf {d}}{\varvec{\chi }}_t}\cdot \nabla ){\mathbf {D}} + (\nabla \overline{{\varvec{\xi }}}_i)\cdot {\mathbf {D}}\circ dW_t^i- {\mathbf {A}}dt\right) \cdot d{\mathbf {x}}. \end{aligned}$$
(3.31)

Note that besides the buoyancy term \({\mathbf {A}}\), also transport terms show up on the right hand side. These transport terms indicate that these fluid equations are not geometric, in the sense that geometric fluid equations will only feature the relevant forces on the right hand side. The reason that these transport terms appear is that strict asymptotics sees the advection constraint (1.31) as two individual terms, rather than as two objects that should always go together. Possibly, a multiscale analysis approach would be able to resolve this problem. This issue is discussed extensively in Gjaja and Holm (1996). We will resolve this issue by linking these two objects in a variational principle for a system closely related to (3.30). First we will investigate the one-dimensional equations related to (3.30).

4.2 Stochastic Benjamin–Bona–Mahony Equations

From the stochastic CH92 equations in (3.30), one cannot derive the stochastic Kadomtsev–Petviashvili equation and further simplify to obtain the stochastic Korteweg–De Vries equation. This is due to the fact that an elimination procedure involving time derivatives was used. However, by restricting to one-dimensional motion, we do obtain the stochastic versions of familiar one-dimensional water wave models. To be able to restrict to one dimension, we ignore the effect of rotation. The variable \({\mathbf {A}}\) is related to the buoyancy at higher order. By replacing \(b_2\) with the vertical average \(\overline{b}_2\) in the definition of \({\mathbf {A}}\) in (3.29), we can explicitly evaluate the integral. In calculating the integral, we keep in mind that the equations are written up to order \(\mathcal O(\gamma ^2)\). This requires us to drop the free surface terms that arise due to the vertical integral. The equations that we obtain from (3.30) are

$$\begin{aligned} \begin{aligned} \frac{1}{\mathrm{Fr}} {\mathsf {d}}\overline{u} - \frac{\gamma ^2}{3\,\mathrm{Fr}^2}h^2 {\mathsf {d}}\overline{u}_{xx} + \overline{{\mathsf {d}}\chi }_t\, \overline{u}_x + \overline{u} \,(\overline{\xi }_i)_x\circ dW_t^i&= -\alpha (1+\widetilde{b}){\mathsf {d}}\eta _x - \frac{\gamma ^2}{2} h^2 (\overline{b}_2)_x dt,\\ \frac{1}{\mathrm{Fr}}{\mathsf {d}}\eta + (\eta \,\overline{{\mathsf {d}}\chi }_t)_x&= 0,\\ \frac{1}{2}h^2 {\mathsf {d}}\overline{b}_2&= \frac{S}{3}h^2 u_{xx}dt. \end{aligned} \end{aligned}$$
(3.32)

The set of equations given by (3.32) can be interpreted as a non–unidirectional, stochastic version of the Benjamin–Bona–Mahony (BBM) equation, first derived in Benjamin et al. (1972), that includes the effects of depth and buoyancy stratification. Since this set (3.32) consists of three equations, we will refer to this set as BBM3. Upon ignoring the effect of buoyancy stratification, we will obtain the two component version of BBM3, which we will call BBM2. This set of equations is given by

$$\begin{aligned} \begin{aligned} \frac{1}{\mathrm{Fr}} {\mathsf {d}}\overline{u} - \frac{\gamma ^2}{3\,\mathrm{Fr}^2}h^2 {\mathsf {d}}\overline{u}_{xx} + \overline{{\mathsf {d}}\chi }_t\, \overline{u}_x + \overline{u} \,(\overline{\xi }_i)_x\circ dW_t^i&= -\alpha (1+\widetilde{b}){\mathsf {d}}\eta _x,\\ \frac{1}{\mathrm{Fr}}{\mathsf {d}}\eta + (\eta \,\overline{{\mathsf {d}}\chi }_t)_x&= 0,\\ \end{aligned} \end{aligned}$$
(3.33)

The two component version (3.33) still is affected by the variations of the free surface. We assume that the bathymetry is flat, which means that we let \(h\mapsto h_0\) and \(h_0\) is constant in space and in time. We also assume that the free surface elevation is zero. These assumptions lead to the stochastic BBM equation, given by

$$\begin{aligned} \frac{1}{\mathrm{Fr}}{\mathsf {d}}\overline{u} - \frac{\gamma ^2}{3\,\mathrm{Fr}^2}h_0^2{\mathsf {d}}\overline{u}_{xx} + \overline{{\mathsf {d}}\chi }_t\, \overline{u}_x + \overline{u}\,(\overline{\xi }_i)_x\circ dW_t^i = 0. \end{aligned}$$
(3.34)

Upon including linear wave speed in formulation of (3.34) and ignoring stochasticity, we arrive at the celebrated BBM equation (Benjamin et al. 1972),

$$\begin{aligned} \frac{1}{\mathrm{Fr}}\overline{u}_t - \frac{\gamma ^2}{3\,\mathrm{Fr}^2}h_0^2 \overline{u}_{xxt} + \overline{u}\,\overline{u}_x + \kappa \overline{u}_x = 0. \end{aligned}$$
(3.35)

Here \(\kappa \) is a positive constant that enforces unidirectionality. The deterministic unidirectional BBM equation (3.35) is similar in shape to the Korteweg–De Vries equation, but is not completely integrable. Next, we consider the averaging procedure in this section from the Euler–Poincaré perspective.

4.3 Averaged Euler–Poincaré Lagrangian for Short Time: Small Wave Scaling

In the previous section, we used direct asymptotics to derive the stochastic version of the equations in Camassa and Holm (1992). These equations failed to satisfy the Kelvin circulation theorem in a reasonable form. This difficulty will be overcome in the Euler–Poincaré approach, because the variational approach is able to cope with arbitrary Strouhal number. The starting point is the thermal rotating Green–Naghdi Lagrangian in (2.29). This time, we are not interested in the rigid lid limit, so our action is given by

$$\begin{aligned} S_{TRGN} = \int _{t_1}^{t_2} \ell _{TRGN}\,dt. \end{aligned}$$
(3.36)

We will now take variations in much the same way as done for the Great Lake equations in the Euler–Poincaré approach. However, there is a crucial difference. In the present scaling regime, the Strouhal number \(\mathrm{Sr}\) is not equal to unity. Instead, we have \(Sr=1/\mathrm{Fr}\), which is the inverse Froude number. Consequently, in the present case, the Euler–Poincaré variations of the velocities are taken as,

$$\begin{aligned} \begin{aligned} \delta \overline{{\mathbf {u}}}\,dt&= \mathrm{Sr}\,{\mathsf {d}}{\mathbf {v}} - [\overline{{\mathsf {d}}{\varvec{\chi }}_t}, {\mathbf {v}}] = \frac{1}{\mathrm{Fr}}{\mathsf {d}}{\mathbf {v}} - [\overline{{\mathsf {d}}{\varvec{\chi }}_t}, {\mathbf {v}}]. \end{aligned} \end{aligned}$$
(3.37)

The averaging has already occured in deriving the Lagrangian (2.29), where the Strouhal number does not explictly appear. In the variational approach, the Strouhal number appears in the variation of the velocity field. We stick with the \(\mathrm{Sr}\) notation to show the flexibility that one has with the variational approach. By selecting the value of the Strouhal number later, the results of the previous section can be recovered by truncating higher order terms. The variational derivatives of the nondimensional Lagrangian \(\ell _{rtGN}\) in equation (2.29) are the following:

$$\begin{aligned} \begin{aligned} \frac{\delta \ell _{TRGN}}{\delta \overline{{\mathbf {u}}}}&= \eta \overline{{\mathbf {u}}} - \frac{\sigma ^2}{3}\nabla (\eta ^3\nabla \cdot \overline{{\mathbf {u}}}) - \frac{\sigma ^2}{2}\nabla \big (\eta ^2(\overline{{\mathbf {u}}}\cdot \nabla h)\big ) + \frac{\sigma ^2}{2}\eta ^2(\nabla \cdot \overline{{\mathbf {u}}})\nabla h + \sigma ^2\eta (\overline{{\mathbf {u}}}\cdot \nabla h)\nabla h + \frac{1}{\mathrm{Ro}}\eta {\mathbf {R}}\,, \\ \frac{\delta \ell _{TRGN}}{\delta \eta }&= \frac{1}{2}|\overline{{\mathbf {u}}}|^2 + \frac{\sigma ^2}{2}(\eta \nabla \cdot \overline{{\mathbf {u}}}+\overline{{\mathbf {u}}}\cdot \nabla h)^2 + \frac{1}{\mathrm{Ro}}(\overline{{\mathbf {u}}}\cdot {\mathbf {R}}) - \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}\overline{b})(\eta -h) - \frac{1}{\mathrm{Fr}^2}{\mathsf {d}}\pi , \\ \frac{\delta \ell _{TRGN}}{\delta \overline{b}}&= -\frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\eta ^2-2\eta h) \,, \\ \end{aligned} \end{aligned}$$
(3.38)

For notational convenience, we define

$$\begin{aligned} h\overline{{\mathbf {V}}} = h\overline{{\mathbf {u}}} + \left[ -\frac{\sigma ^2}{3}\nabla (\eta ^3\nabla \cdot \overline{{\mathbf {u}}}) - \frac{\sigma ^2}{2}\nabla (\eta ^2\overline{{\mathbf {u}}}\cdot \nabla h) + \frac{\sigma ^2}{2}\eta ^2(\nabla \cdot \overline{{\mathbf {u}}})\nabla h + \sigma ^2 \eta (\overline{{\mathbf {u}}}\cdot \nabla h)\nabla h\right] .\nonumber \\ \end{aligned}$$
(3.39)

In the rigid lid case, one recovers (2.24). A careful application of the Lax–Milgram theorem is able to show that \(\overline{{\mathbf {u}}}\) depends continuously on \(\overline{{\mathbf {V}}}\). Using this notation, an application of the stochastic Euler–Poincaré theorem 1.1 with the velocity variations given in (3.37) and the variational derivatives in (3.38) of the Lagrangian \(\ell _{TRGN}\) in (2.29) yields the following SPDEs,

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\overline{{\mathbf {V}}} + (\overline{{\mathsf {d}}{\varvec{\chi }}_t}\cdot \nabla )\overline{{\mathbf {V}}} + (\nabla \overline{{\mathsf {d}}{\varvec{\chi }}_t})\cdot \overline{{\mathbf {V}}}&= - \frac{\alpha }{\mathrm{Fr}^2}\nabla \big ((1+\mathfrak {s}\overline{b})\zeta \big )dt + \frac{1}{2}\nabla |\overline{{\mathbf {u}}}|^2 dt + \frac{\sigma ^2}{2}\nabla (\eta \nabla \cdot \overline{{\mathbf {u}}} + \overline{{\mathbf {u}}}\cdot \nabla h)^2\,dt \\&\quad + \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\alpha \zeta - h)\nabla \overline{b}\,dt - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}_t} \\&\quad - \frac{1}{\mathrm{Ro}}\nabla (\overline{\varvec{\xi }}_i\cdot {\mathbf {R}})\circ dW_t^i ,\\ \mathrm{Sr}\,\alpha {\mathsf {d}}\zeta + \nabla \cdot \big ((\alpha \zeta +h) \overline{{\mathsf {d}}{\varvec{\chi }}_t}\big )&= 0 ,\\ \mathrm{Sr}\,{\mathsf {d}}\overline{b} + \overline{{\mathsf {d}}{\varvec{\chi }}_t}\cdot \nabla \overline{b}&= 0. \end{aligned} \end{aligned}$$
(3.40)

where \(\overline{{\mathsf {d}}{\varvec{\chi }}_t} \) is defined in equation (3.9). It is useful to note that \(\eta ^{-1}\delta \ell _{TRGN}/\delta \overline{b} = (\mathfrak {s}/2)(\eta - 2h) = (\mathfrak {s}/2)(\alpha \zeta - h)\), since \(\eta = \alpha \zeta + h\). These equations do satisfy a Kelvin circulation theorem, as they have been derived from the Euler–Poincaré variational principle. The circulation theorem takes the following form

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\mathsf {d}} \oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})} \left( \overline{{\mathbf {V}}} + \frac{1}{\mathrm{Ro}}{\mathbf {R}} \right) \cdot d{\mathbf {x}}&= \oint _{c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})} \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\alpha \zeta - h) \nabla \overline{b} \cdot d{\mathbf {x}} \\ {}&= \int \!\!\int _{\partial S = {c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})}} \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2} \nabla (\alpha \zeta - h)\times \nabla \overline{b}\cdot d{\mathbf {S}} \,dt\,. \end{aligned} \end{aligned}$$
(3.41)

As expected from equations (1.20) and (1.34) for the Kelvin circulation theorem which follows from the Euler–Poincaré equation (1.15) in three dimensions, circulation is created by misalignment of the gradients of vertically averaged buoyancy \(\overline{b}\) and its dual quantity \(\eta ^{-1}\delta \ell _{TRGN}/\delta \overline{b}\), for the thermal rotating Green–Naghdi Lagrangian in equation (2.29). This is a balanced statement, because gradients of the bathymetry are assumed to be small. Interestingly, the misalignment of the gradient of vertically averaged buoyancy \(\overline{b}\) and the difference \((\alpha \zeta - h)\) generates horizontal circulation (vertical vorticity). This represents a barotropic mechanism for cyclogenesis (emergence of horizontal circulation, or eddies) in the ocean. The dispersion relation that corresponds to the linearised, deterministic version of equations (3.40) is discussed in Appendix A. A Kelvin circulation theorem similar to that in (3.41) holds for the thermal rotating shallow water (TRSW) equations, as discussed in Appendix B.

Remark 3.2

(Comparison with JEBAR for ocean currents). For the deterministic case, one replaces \(c(\overline{{\mathsf {d}}{\varvec{\chi }}_t})\rightarrow c(\overline{{\mathbf {u}}})\) and the circulation theorem in (3.41) recalls an aspect of the JEBAR (Joint Effect of Baroclinicity and Bottom Relief) approach for modelling the dynamics of ocean currents (Sarkisyan and Ivanov 1971; Cane et al. 1998; Mellor 1999; Sarkisyan 2006; Colin de Verdière and Ollitrault 2016). Namely, the creation of circulation in (3.41) occurs when the gradients of certain fluid properties are not aligned with the gradient of the bottom topography, \(\nabla h({\mathbf {x}})\).

There are also may differences of (3.41) from JEBAR. In particular, the circulation dynamics in (3.41) represents Kelvin’s theorem as derived from a vertically averaged and asymptotically expanded Hamilton’s principle for Euler’s fluid equations for the stochastic dynamics of an incompressible, thermal, rotating fluid flow with a free upper surface moving under the influence of gravity. Nonetheless, many of the physical principles underlying the derivation of (3.41) also relate to principles which could be applied in the oceanographic setting for JEBAR. Hence, it may be advisable to investigate the utility of the present stochastic, asymptotic, vertically averaged variational approach for some applications in oceanography.

Potential vorticity. In the circulation theorem for the rotating, thermal, Great Lake equations in equation (2.25), the circulation is generated by the misalignment between the horizontal gradient of the bathymetry and the horizontal gradient of the buoyancy. Here, we have seen that the misalignment of horizontal gradients of the free surface height with the horizontal gradient of the buoyancy also contributes to the generation of circulation. In terms of the potential vorticity given by

$$\begin{aligned} q := \eta ^{-1}\big ({{\hat{\mathbf {z}}}}\cdot \nabla \times (\overline{{\mathbf {V}}}+\mathrm{Ro}^{-1}{\mathbf {R}})\big ) \,, \end{aligned}$$
(3.42)

the generation of circulation is accompanied by the following

$$\begin{aligned} \mathrm{Sr}\,{\mathsf {d}}q + (\overline{{\mathsf {d}}{\varvec{\chi }}}_t\cdot \nabla )q = \frac{\mathfrak {s}}{2\mathrm{Fr}^2\eta }\,{{\hat{\mathbf {z}}}}\cdot \nabla (\alpha \zeta -h)\times \nabla \overline{b}. \end{aligned}$$
(3.43)

This shows that PV will also be generated by this misalignment of horizontal gradients. Equations (3.40) also possess an infinity of conserved integral quantities of the following form

$$\begin{aligned} C_{f,g} = \int _{CS} \big ( f(\overline{b}) + qg(\overline{b})\big )\eta \, dx dy, \end{aligned}$$
(3.44)

for arbitrary differentiable functions f, g and for boundary conditions \({\mathsf {d}}\overline{{\varvec{\chi }}}_t \cdot {\mathbf {n}} = 0\), \(\nabla \overline{b}\times {\mathbf {n}}=0\). Invariance of the vertically averaged buoyancy \(\overline{b}\) as it is advected along the tangential stochastic flow on the boundary is consistent with the latter condition, which requires the boundary to be a level set of \(\overline{b}\). This can be shown by means of a direct computation using the equations of motion and the boundary conditions.

4.4 Stochastic Camassa–Holm Equations

This section considers a sequence of reductions of the Lagrangian \(\ell _{TRGN}(\overline{{\mathbf {u}}}, \eta , \overline{b})\) in equation (2.29) in one spatial dimension which will eventually lead to the stochastic Camassa–Holm (CH) equation, considered in Holm and Tyranowski (2016), Crisan and Holm (2018)

$$\begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\overline{m} + \big (\overline{m}\partial _x + \partial _x\overline{m}\big )\overline{{\mathsf {d}}\chi _t} = 0\,. \end{aligned}$$
(3.45)

In one dimension, we assume a flat bathymetry profile \(h_0\) and ignore the effect of rotation. Applying these approximations to the thermal Green–Naghdi Lagrangian \(\ell _{TRGN}(\overline{{\mathbf {u}}}, \eta , \overline{b})\) in equation (2.29) yields the following Lagrangian at order \(\mathcal O(\sigma ^2)\),

$$\begin{aligned} \ell _{CH3} = \int _{-\infty }^{\infty } \left( \frac{1}{2}\overline{u}^2 + \frac{\sigma ^2}{6} \eta ^2 \overline{u}_x^2 - \frac{1}{2\,\mathrm{Fr}^2}(1+\mathfrak {s}\overline{b})(\eta - 2h_0)\right) \eta \,dx, \end{aligned}$$
(3.46)

where we have completed the square on the potential energy term. The domain of flow is taken to be the entire real line, rather than a compact line between two lateral boundaries as illustrated in Fig. 2. Boundary conditions on the real line require the vertically averaged velocity \(\overline{u}\) and its horizontal spatial derivative \(\overline{u}_x\) to vanish in the limit \({|x|\rightarrow \infty }\). The variational derivatives of the Lagrangian \(\ell _{CH3}\) in (3.46) are given by

$$\begin{aligned} \begin{aligned} \overline{m}:= \frac{\delta \ell _{CH3}}{\delta \overline{u}}&= \eta \overline{u} - \frac{\sigma ^2}{3}(\eta ^3\overline{u}_{x})_x, \\ \frac{\delta \ell _{CH3}}{\delta \eta }&= \frac{1}{2}\overline{u}^2 + \frac{\sigma ^2}{2}\eta ^2\overline{u}_x^2 - \frac{1}{\mathrm{Fr}^2}(1+\mathfrak {s}\overline{b})(\eta -h_0), \\ \frac{\delta \ell _{CH3}}{\delta \overline{b}}&= -\frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\eta ^2-2\eta h_0). \end{aligned} \end{aligned}$$
(3.47)

An application of the stochastic Euler–Poincaré theorem 1.1 then leads to the following set of three stochastic equations

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\overline{m} + \big (\overline{m}\partial _x + \partial _x\overline{m}\big )\overline{{\mathsf {d}}\chi }_t&= -\frac{1}{\mathrm{Fr}^2}\eta \big ((1+\mathfrak {s}\overline{b})(\eta -h_0)\big )_x\,dt + \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\eta ^2-2\eta h_0)\overline{b}_x\,dt\\&\quad + \eta \left( \frac{1}{2}\overline{u}^2 + \frac{\sigma ^2}{2}\eta ^2\overline{u}_x^2\right) _x \,dt, \\ \mathrm{Sr}\,{\mathsf {d}}\eta + (\eta \, \overline{{\mathsf {d}}\chi }_t)_x&= 0, \\ \mathrm{Sr}\,{\mathsf {d}}\overline{b} + \overline{{\mathsf {d}}\chi }_t\,\overline{b}_x&= 0. \end{aligned} \end{aligned}$$
(3.48)

The set of equations (3.48) defines the three-component stochastic Camassa–Holm system (CH3). The stochastic evolution equation for momentum \(\overline{m}\) includes the effects of varying depth and horizontal variations of the buoyancy. There follows a continuity equation for depth, \(\eta \), and a scalar advection equation for buoyancy, \(\overline{b}\).

Remark 3.3

(Is the deterministic CH3 case completely integrable?). An investigation is underway elsewhere to determine whether the Lie–Poisson Hamiltonian system of CH3 equations in (3.48) is completely integrable in the deterministic case, where it simplifies to

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\partial _t}\overline{m} + \big (\overline{m}\partial _x + \partial _x\overline{m}\big )\overline{u}&= -\frac{1}{\mathrm{Fr}^2}\eta \big ((1+\mathfrak {s}\overline{b})(\eta -h_0)\big )_x + \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\eta ^2-2\eta h_0)\overline{b}_x \\&\quad + \eta \left( \frac{1}{2}\overline{u}^2 + \frac{\sigma ^2}{2}\eta ^2\overline{u}_x^2\right) _x\,, \\ \mathrm{Sr}\,{\partial _t}\eta + (\eta \, \overline{u} )_x&= 0, \\ \mathrm{Sr}\,{\partial _t}\overline{b} + \overline{u}\,\overline{b}_x&= 0. \end{aligned} \end{aligned}$$
(3.49)

We proceed farther now in the stochastic case by assuming that the vertically averaged buoyancy \(\overline{b}\) is constant in both space and time, so that we may replace \(\overline{b}(x,t)\mapsto b_0\); a constant, or equivalently, by letting the stratification parameter tend to zero, \(\mathfrak {s}\rightarrow 0\). Under this assumption, the Lagrangian \(\ell _{CH3}\) simplifies, since the buoyancy term no longer contributes to the dynamics, and we arrive at the following Lagrangian \(\ell _{CH2}\) for the stochastic two component Camassa–Holm (CH2) system:

$$\begin{aligned} \ell _{CH2} = \int _{-\infty }^{\infty } \left( \frac{1}{2}\overline{u}^2 + \frac{\sigma ^2}{6}\eta ^2 \overline{u}_x^2 - \frac{1}{2\,\mathrm{Fr}^2}(\eta - 2h_0)\right) \eta \, dx. \end{aligned}$$
(3.50)

The variational derivatives of the Lagrangian \(\ell _{CH2}\) in (3.50) are given by

$$\begin{aligned} \begin{aligned} \overline{m}:=\frac{\delta \ell _{CH2}}{\delta \overline{u}}&= \eta \overline{u} - \frac{\sigma ^2}{3}(\eta ^3\overline{u}_x)_x, \\ \frac{\delta \ell _{CH2}}{\delta \eta }&= \frac{1}{2}\overline{u}^2+\frac{\sigma ^2}{2}\eta ^2\overline{u}_x^2-\frac{1}{\mathrm{Fr}^2}(\eta -h_0). \end{aligned} \end{aligned}$$
(3.51)

An application of the stochastic Euler–Poincaré theorem 1.1 with these variational derivatives yields the following motion equation and advection law,

$$\begin{aligned} \begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\overline{m} + (\overline{m}\partial _x+\partial _x\overline{m}){\mathsf {d}}\overline{\chi }_t&= -\frac{1}{\mathrm{Fr}^2}\eta \eta _x\,dt + \eta \left( \frac{1}{2}\overline{u}^2 + \frac{\sigma ^2}{2}\eta ^2\overline{u}_x^2\right) _x\,dt,\\ {\mathsf {d}}\eta + (\eta \,{\mathsf {d}}\overline{\chi }_t)_x&= 0. \end{aligned} \end{aligned}$$
(3.52)

The set of equations (3.52) is closely related to the stochastic two component Camassa–Holm (CH2) system. The difference is that the usual CH2 system does not include the kinetic energy term, that is, the last term in on the right hand side of the momentum equation in (3.52). Moreover, the definition of the momentum variable \(\overline{m}\) in (3.51) features an \(\eta ^3\)-weighted Helmholtz operator, whereas in the usual CH2 equations, the Helmholtz operator does not include a weight. In the deterministic case, this set of equations is a completely integrable Hamiltonian system, as shown first by Chen et al. (2006).

Finally, we will assume that the free surface elevation in the CH2 Lagrangian \(\ell _{CH2}\) in (3.50) is negligible. This assumption neglects the potential energy term in \(\ell _{CH2}\), which then reduces to

$$\begin{aligned} \ell _{CH} = \int _{-\infty }^{\infty } \left( \frac{1}{2}\overline{u}^2 + \frac{\sigma ^2}{6}h_0^2 \overline{u}_x^2\right) h_0 \, dx. \end{aligned}$$
(3.53)

The variation of the CH Lagrangian (3.53) with respect to the velocity \(\overline{u}\) yields

$$\begin{aligned} \overline{m}:=\frac{\delta \ell _{CH}}{\delta \overline{u}} = h_0\overline{u} - \frac{\sigma ^2}{3}h_0^3\overline{u}_{xx} \end{aligned}$$
(3.54)

An application of the stochastic Euler–Poincaré theorem 1.1 then implies the SPDE,

$$\begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\overline{m} + \big (\overline{m}\partial _x + \partial _x\overline{m}\big )\overline{{\mathsf {d}}\chi }_t = 0. \end{aligned}$$
(3.55)

Equation (3.55) is the dispersionless stochastic Camassa–Holm equation, whose singular ‘peakon’ solutions have been studied in Holm and Tyranowski (2016); Crisan and Holm (2018). Including cubic linear dispersion in the stochastic Camassa–Holm equation yields

$$\begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\overline{m} + \big (\overline{m}\partial _x + \partial _x\overline{m} + \gamma \partial _x^3\big )\overline{{\mathsf {d}}\chi }_t = 0\,. \end{aligned}$$
(3.56)

The solution properties of this equation has been studied in Holm and Tyranowski (2016); Bendall et al. (2019). When terms of order \(O(\sigma ^2)\) are neglected in equation (3.56), it reduces further to the stochastic KdV equation,

$$\begin{aligned} \mathrm{Sr}\,{\mathsf {d}}\overline{u} + \big (\overline{u}\partial _x + \partial _x\overline{u} + \gamma \partial _x^3\big )\overline{{\mathsf {d}}\chi }_t = 0\,, \end{aligned}$$
(3.57)

which has been studied in Woodfield (2019). The deterministic CH equation was first derived in Camassa and Holm (1993); Camassa et al. (1994), by using asymptotics on the Hamiltonian side. Here, the stochastic CH equation has been derived by means of asymptotics in the Lagrangian for the rotating, thermal, Green–Naghdi equations (2.29) followed by applying the stochastic Euler–Poincaré theorem to the approximated Lagrangian at a variety of levels.

4.5 Differences Between the Newtonian and Variational Approaches

There are several striking differences between the equations that one derives from the Newtonian approach and from the Euler–Poincaré approach, as illustrated with underbraces below. The most important difference is that the time derivative of \({\mathbf {D}}\) no longer appears explicitly in the equations above. Instead, the dynamical variable \({\mathbf {V}}\) appears naturally, as it did for the Great Lake equations in (2.25). The pressure and the buoyancy term also take slightly different forms. The averaged equations (3.8) indicate that the usage of the buoyancy equation is natural. In the Newtonian approach, the buoyancy only has dynamics at order \(\sigma ^4\), since \(b_1\) was calculated explicitly and shown only to depend on the vertical coordinate. This explains the sole appearance of \(b_2\) in the buoyancy equation. In the variational approach, we do not calculate the explicit profile of \(b_1\), but instead we introduce a vertically averaged buoyancy in the Lagrangian. This means that the buoyancy is still allowed to vary horizontally, which can be seen in the equation for the buoyancy. The effect of the horizontal dependence of the buoyancy is important for the generation of horizontal circulation, as noticed in (3.41). Below we have expressed the two sets of equations in terms of the same variables so that the differences and similarities are clear.

CH92 equations:

$$\begin{aligned} \begin{aligned} \frac{1}{\mathrm{Fr}}{\mathsf {d}}\overline{{\mathbf {u}}} - \frac{\sigma ^2}{3\,\mathrm{Fr}}h^2{\mathsf {d}}\nabla (\nabla \cdot \overline{{\mathbf {u}}})&+ (\overline{{\mathsf {d}}{\varvec{\chi }}}_t\cdot \nabla )\overline{{\mathbf {u}}} + (\nabla \overline{{\mathsf {d}}\varvec{\chi }}_t)\cdot \overline{{\mathbf {u}}} \\ {}&= -\frac{\alpha }{\mathrm{Fr}^2}\nabla \big ((1+\widetilde{b})\zeta \big )\,dt + \frac{1}{2}\nabla |\overline{{\mathbf {u}}}|^2dt - \overline{(z+h)\nabla b_2}\,dt \\&\quad - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}}_t - \frac{1}{\mathrm{Ro}}\nabla (\overline{{\varvec{\xi }}}_i\cdot {\mathbf {R}})\circ dW_t^i, \\ \frac{\alpha }{\mathrm{Fr}}{\mathsf {d}}\zeta + \nabla \cdot \big ((\alpha \zeta + h)\overline{{\mathsf {d}}{\varvec{\chi }}}_t\big )&= 0, \\ {\mathsf {d}}\overline{(z+h)\nabla b_2}&= \frac{S}{3}h^2\nabla (\nabla \cdot \overline{{\mathbf {u}}})dt. \end{aligned} \end{aligned}$$
(3.58)

Thermal rotating Green–Naghdi equations:

$$\begin{aligned} \begin{aligned} h\overline{{\mathbf {V}}} = h\overline{{\mathbf {u}}} + \left[ -\frac{\sigma ^2}{3}\nabla (\eta ^3\nabla \cdot \overline{{\mathbf {u}}}) \right.&\left. - \frac{\sigma ^2}{2}\nabla (\eta ^2\overline{{\mathbf {u}}}\cdot \nabla h) + \frac{\sigma ^2}{2}\eta ^2(\nabla \cdot \overline{{\mathbf {u}}})\nabla h + \sigma ^2 \eta (\overline{{\mathbf {u}}}\cdot \nabla h)\nabla h\right] ,\\ \mathrm{Sr}\,{\mathsf {d}}\overline{{\mathbf {V}}} + (\overline{{\mathsf {d}}{\varvec{\chi }}}_t\cdot \nabla )\overline{{\mathbf {V}}} + (\nabla \overline{{\mathsf {d}}{\varvec{\chi }}}_t)\cdot \overline{{\mathbf {V}}}&= - \frac{\alpha }{\mathrm{Fr}^2}\nabla \big ((1+\mathfrak {s}\overline{b})\zeta \big )dt + \frac{1}{2}\nabla |\overline{{\mathbf {u}}}|^2 dt + \frac{\sigma ^2}{2}\nabla (\eta \nabla \cdot \overline{{\mathbf {u}}} + \overline{{\mathbf {u}}}\cdot \nabla h)^2\,dt \\&\quad + \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\alpha \zeta - h)\nabla \overline{b}\,dt - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}}_t - \frac{1}{\mathrm{Ro}}\nabla (\overline{\varvec{\xi }}_i\cdot {\mathbf {R}})\circ dW_t^i ,\\ \mathrm{Sr}\,\alpha {\mathsf {d}}\zeta + \nabla \cdot \big ((\alpha \zeta +h) \overline{{\mathsf {d}}{\varvec{\chi }}}_t\big )&= 0 ,\\ \mathrm{Sr}\,{\mathsf {d}}\overline{b} + \overline{{\mathsf {d}}{\varvec{\chi }}}_t\cdot \nabla \overline{b}&= 0. \end{aligned} \nonumber \\ \end{aligned}$$
(3.59)

where \(\overline{{\mathsf {d}}{\varvec{\chi }}_t} \) is defined in equation (3.9). Evaluating the Strouhal number \(\mathrm{Sr} = \frac{1}{\mathrm{Fr}}\) and truncating at order \(\mathcal {O}(1)\) in (3.59) provides the following set of equations

$$\begin{aligned} \begin{aligned} \frac{1}{\mathrm{Fr}}{\mathsf {d}}\overline{{\mathbf {u}}} - \frac{\sigma ^2}{3\,\mathrm{Fr}}h^2{\mathsf {d}}\nabla (\nabla \cdot \overline{{\mathbf {u}}})&+ (\overline{{\mathsf {d}}{\varvec{\chi }}}_t\cdot \nabla )\overline{{\mathbf {u}}} + (\nabla \overline{{\mathsf {d}}{\varvec{\chi }}}_t)\cdot \overline{{\mathbf {u}}} \\&= - \frac{\alpha }{\mathrm{Fr}^2}\nabla \big ((1+\mathfrak {s}\overline{b})\zeta \big )dt + \frac{1}{2}\nabla |\overline{{\mathbf {u}}}|^2 dt + \frac{\mathfrak {s}}{2\,\mathrm{Fr}^2}(\alpha \zeta - h)\nabla \overline{b}\,dt \\ {}&\quad - \frac{1}{\mathrm{Ro}}f{\hat{\mathbf {z}}}\times \overline{{\mathsf {d}}{\varvec{\chi }}}_t - \frac{1}{\mathrm{Ro}}\nabla (\overline{\varvec{\xi }}_i\cdot {\mathbf {R}})\circ dW_t^i ,\\ \frac{\alpha }{\mathrm{Fr}}{\mathsf {d}}\zeta + \nabla \cdot \big ((\alpha \zeta +h) \overline{{\mathsf {d}}{\varvec{\chi }}}_t\big )&= 0 ,\\ \frac{1}{\mathrm{Fr}}{\mathsf {d}}\overline{b} + \overline{{\mathsf {d}}{\varvec{\chi }}}_t\cdot \nabla \overline{b}&= 0. \end{aligned} \end{aligned}$$
(3.60)

There are still some differences between (3.58) and (3.60). In the variational approach, we introduce the vertically averaged buoyancy which gives rise to terms that create horizontal circulation, rather than introducing an explicit profile. The original CH92 equations in (3.58) were derived in Camassa and Holm (1992) by applying vertical averaging and strict asymptotics in the unapproximated equations in the form of Newton’s force law for the fluid. Asymptotics in Strouhal number breaks the Kelvin circulation theorem. The thermal rotating Green–Naghdi equations in (3.59) have the following advantages over the CH92 equations

  1. 1.

    They introduce a dynamical equation for the vertically averaged buoyancy, \(\overline{b}\);

  2. 2.

    The dynamics of the vertically averaged buoyancy, \(\overline{b}\), contributes to the pressure terms;

  3. 3.

    They restore the Kelvin circulation theorem seen in equation (3.40);

  4. 4.

    They reveal a barotropic mechanism for horizontal circulation (cyclogenesis), as seen in equation (3.40); and

  5. 5.

    They allow for a hierarchy of Camassa–Holm equations to be derived, see Secti. 3.4.

5 Conclusion

Summary. This paper has extended the work of Camassa and Holm (1992) and Camassa et al. (1996, 1997) by casting it into the framework of Hamilton’s variational principle and including the multi-time effects of the Strouhal number and the barotropic effects of vertically integrated buoyancy with horizontal gradients. As a result, a variety of new terms representing new effects relative to Camassa and Holm (1992) and Camassa et al. (1996, 1997) have appeared in the resulting equations. For example, in the variational CH92 equations (3.40) written in Kelvin circulation form in (3.41) one sees how horizontal circulation (convection) is generated by an misalignment of horizontal gradient of vertically averaged buoyancy with the horizontal gradients of bathymetry and/or surface elevation.

Having extended the earlier work of Camassa and Holm (1992) and Camassa et al. (1996, 1997) in a variational setting and expressed the results in Kelvin circulation form, the paper has also taken advantage of the variational framework of Holm (2015) to include the effects of stochastic advective Lie transport (SALT). Including the effects of SALT introduces a new capability to quantify the uncertainty and then use data assimilation to reduce the uncertainty of the solutions of these equations due to unmodelled, or unresolved effects. A protocol for doing this has been been developed in Cotter et al. (2018, 2019a, 2019b). This protocol regards SALT as a type of ‘informed randomness’ described by spatially correlated noise obtained from observed or simulated high-resolution data. This protocol may be applied to the present class of fluid equations. In order to reduce the investigation of these equations to their simplest form, the paper has derived the unidirectional version of the equation set in (3.40) in the variational setting. This reduction has yielded stochastic versions of a family of CH equation, including the one derived in Camassa and Holm (1993); Camassa et al. (1994). These stochastic CH equations describe the interaction of solitons with noise. The first developments in this direction for the stochastic CH equation have already been studied in Holm and Tyranowski (2016, 2018), Crisan and Holm (2018) and Bendall et al. (2019).

Two diagrams sketched below provide ‘roadmaps’ of the two routes of simplification we have taken in this paper by using asymptotic expansions in the various small parameters for the ordering scheme in equation (3.7). The Newtonian approach is shown in Fig. 3. The corresponding road map for the variational approach is shown in Fig. 4.

Fig. 3
figure 3

Diagram of derivations from the direct (or Newtonian) point of view. Each blue box refers to the set of equations that corresponds to the model referred to in the box. Clicking on the box will take the reader to the corresponding section. Above each arrow is the approximation that is necessary to transition from one set of equations to the next. Note that the short time-small wave approximation does not lead to rotating thermal Green–Naghdi, but to the CH92 equations. These lead to Benjamin-Bona-Mahony type equations when restricted to one-dimensional motion. The rotating, stratified Euler equations are not linked because these equations have not been written down in this document

Fig. 4
figure 4

Diagram of derivations from the variational point of view. Each blue box refers to the Lagrangian that corresponds to the model referred to in the box. By clicking on the box the reader is taken to the corresponding section. Above each arrow is the approximation that is necessary to transition from one Lagrangian to the next

In Sect. 1, we investigated whether the SALT approach was compatible with the asymptotic expansions. It was shown that an additional assumption on the magnitude of the gradient of the bathymetry was required for the SALT version to be consistent with the deterministic situation. Except for this additional assumption, SALT was verified to be compatible with the methods of asymptotic analysis. From the variational point of view, this was to be expected. Any fluid model which has a corresponding Lagrangian can be made stochastic with the approach of Holm (2015). However, boundary conditions need to be made consistent with the derivation of the equations. A simpler, but also important ‘sanity check’ was passed, by confirming that the stochastic Lake and Great Lake equations successfully recover the deterministic Lake and Great Lake equations when the noise terms are absent.

In Sect. 2, we showed that the Great Lake equations in (2.25) may be derived using a direct approach, by combining vertical averaging of the nondimensional Euler–Boussinesq equations with asymptotic analysis in a long time - very small wave scaling regime. The resulting averaged equations can be closed. One may also derive the same equations by vertically averaging the Lagrangian and applying the Euler–Poincaré theorem. In both situations, an averaging principle is required which respects the boundary conditions for the Euler–Boussinesq equations. The road map of these derivations is sketched on the right-hand branches of Figs. 3 and 4.

In Sect. 3, we worked in a short time - small wave scaling regime, following the left-hand branches of Figs. 3 and 4. In this scaling regime, the Strouhal number does not equal unity. Instead, the Strouhal number is the inverse of the Froude number, which was taken to be small in this scaling regime. Consequently, the material derivative was no longer balanced in the asymptotic expansion. Because of this imbalance, the direct asymptotic expansion approach failed to derive the rotating thermal Green–Naghdi equations in this scaling regime. However, the variational approach was able to take an arbitrary Strouhal number into account. In this scaling regime, the variational approach provided a set of equations reminiscent of the Green–Naghdi equations, and which had the geometric structure required to possess a Kelvin circulation theorem. Thus, the Strouhal number played a crucial role in determining the differences between the direct approach and the variational approach in the short time - small wave scaling regime. In addition, by further approximating the asymptotic expansion of the wave Lagrangian in Hamilton’s principle, in Sect. 3.4 we derived several stochastic variants of the Camassa–Holm equation and the Korteweg–de Vries equation for one-dimensional unidirectional propagation. Finally, in Sect. 3.5 we discussed the differences between the Newtonian and variational approaches in this scaling regime by making a detailed comparison of the equations and explaining the implications of the additional terms in the variational approach which were missing in the direct approach.

5.1 Outlook and Open Problems: What to Do?

This paper has integrated several methodologies into a research framework for investigating the various effects of wave–current interaction in thermal shallow water flows. Several methodologies were required because wave–current interaction involves several elements. Different time scales exist for flow and wave propagation, as indicated by the different regimes of Strouhal number. This means that simultaneous interactions take place among various physical effects with different times scales. For example, we have seen that nonlinear interactions arise among advective transport, dispersive nonlinear wave propagation, stratification and generation of circulation in the interplay of waves, topography and stratification. This is not to even mention the effects of shear on the propagation of waves and the effects of wave perturbations on unstable flow equilibria.

Because of these various interacting elements, modelling the wave–current interaction process involves many uncertainties. These uncertainties arise from the combination of incomplete sparse observations and the ‘irreducible imprecision’ of numerical simulations arising because of under-resolution and the wide variety of choice in numerical simulation algorithms. In the hopes of providing a methodology for systematically quantifying these uncertainties, this paper has introduced stochastic advection by Lie transport (SALT) in the derivation of the various new equations arising in the ramifications of the asymptotic expansions studied here. We believe that the SALT approach could eventually be made useful for stochastic parameterisation and uncertainty quantification of wave–current interaction, for example, in describing the effects of sub-mesoscale unresolved ocean dynamics on the larger, slower, resolvable oceanic flow. Combined with judicious data assimilation approaches based on the earlier work of Cotter et al. (2018, 2019a, 2019b), one can hope that in some cases these uncertainties may even be reduced. The progress made here suggests that further pursuit of the SALT approach for stochastic parameterisation may soon be fruitful in the context of wave–current interaction of dispersive nonlinear waves in shallow water with horizontal buoyancy gradients. In the mean time, the present paper has combined asymptotic expansions and vertical averaging with the stochastic variational framework to formulate the SALT approach for the various thermal shallow water equations which descend from Euler’s three-dimensional fluid equations under approximation by asymptotic expansions and vertical averaging.