Stochastic Closures for Wave–Current Interaction Dynamics

Wave–current interaction (WCI) dynamics energizes and mixes the ocean thermocline by producing a combination of Langmuir circulation, internal waves and turbulent shear flows, which interact over a wide range of time scales. Two complementary approaches exist for approximating different aspects of WCI dynamics. These are the Generalized Lagrangian Mean (GLM) approach and the Gent–McWilliams (GM) approach. Their complementarity is evident in their Kelvin circulation theorems. GLM introduces a wave pseudomomentum per unit mass into its Kelvin circulation integrand, while GM introduces an additional ‘bolus velocity’ to transport its Kelvin circulation loop. The GLM approach models Eulerian momentum, while the GM approach models Lagrangian transport. In principle, both GLM and GM are based on the Euler–Boussinesq (EB) equations for an incompressible, stratified, rotating flow. The differences in their Kelvin theorems arise from differences in how they model the flow map in the Lagrangian for the Hamilton variational principle underlying the EB equations. A recently developed approach for uncertainty quantification in fluid dynamics constrains fluid variational principles to require that Lagrangian trajectories undergo Stochastic Advection by Lie Transport (SALT). Here, we introduce stochastic closure strategies for quantifying uncertainty in WCI by adapting the SALT approach to both the GLM and GM approximations of the EB variational principle. In the GLM framework, we introduce a stochastic group velocity for transport of wave properties, relative to the frame of motion of the Lagrangian mean flow velocity and a stochastic pressure contribution from the fluctuating kinetic energy. In the GM framework, we introduce a stochastic bolus velocity in addition to the mean drift velocity by imposing the SALT constraint in the GM variational principle.


Introduction
The wind drives gravity waves on the ocean surface. Over time, the collective action of these wind-driven gravity waves on the ocean surface generates Langmuir circulations (LC) which transport heat and mix material properties deeper into the ocean. The presence of LC is seen as 'lines on the sea surface' marked by flotsam trapped between roughly parallel, horizontally counter-circulating pairs of Langmuir vortex rolls. Eventually, these wave-current interactions energize and mix the ocean surface boundary layers (OSBLs) which occupy the upper few hundred meters of the ocean. In turn, the well-mixed region at the top of the OSBLs comprises the thermocline. Just below it, the stratified regions propagate internal waves which further transmit and disperse wave activity.
The turbulent wave-current mixing by Langmuir circulation seen in the OSBL is important in climate modeling, because it controls the exchange of heat and trace gases between the atmosphere and ocean through the mix layer. However, a difficulty arises in numerically simulating the regional effects of Langmuir circulation on turbulent mixing in OSBL, because of the huge disparities among length and times scales of the waves, currents, regional flows and their effects on climate. Such huge disparities make direct numerical simulations (DNS) of turbulent mixing by wave, current interaction intractable for any existing or projected computer for decades to come.
For comprehensive reviews of modern approaches for quantifying the dynamics of Lagrangian flows such as Langmuir circulations coupled to surface and internal waves, see, e.g., Sullivan and McWilliams (2010), Phillips (2003), Fujiwara et al. (2018) as well as references therein.
Current parameterizations of turbulent mixing in numerical simulations of the OSBL lead to substantial systematic errors, for example, in predicting the depth of the OSBL for a given wind stress. These errors, in turn, lead to further uncertainty in predictions of sea surface temperature and rate of exchange of gases such as CO 2 between the ocean and the atmosphere, (Belcher et al. 2012).
Because of the computational intractability due to the enormous scale disparity and the space-time distributed nature of wave-current interactions with weather and climate dynamics, simulations of turbulent mixing in OSBL are always carried out in regions of parameter space which are far from the observed values, either with: (a) an unphysical lack of scale separation between the energy-containing, inertial, and dissipative scales while parameterizing the missing physics, or with (b) a study of the processes at much smaller length scales, often with periodic boundaries (unphysical at large scales but used under the hypothesis of spatial homogeneity of the flows). Moreover, because of the nonlinear nature of turbulent flows and the ensuing multiscale, space-time distributed interactions, the physics of the unresolvable, rapid, small scales may differ significantly from the properties (e.g., statistics) of the resolvable large scales. For example, the regime of asymptotic expansions for the large-scale computational models occurs at small Rossby number, which enforces hydrostatic and geostrophic balances. However, for wave-current interaction (WCI) at the submesoscale length scales below the Rossby radius where Langmuir circulations develop, the Rossby number is order O(1) and neither of these large-scale balances is enforced. This imbalance requires another model. Given this situation, there is clearly a need for enhanced methods for parameterizing the effects on the resolvable scales of the unresolvable small scales in space and time. Two main parameterization approaches have been developed over the years to model the effects of the unresolvable small scales in turbulence on the scales resolved in the simulations.
The first parameterization approach is primarily computational, via Large Eddy Simulations (LES). LES is widely used in engineering, in atmospheric sciences, and to a lesser extent in astrophysics. However, in the LES approach, many important physical parameters for the Langmuir circulations are not scale-appropriate. For example, in the LES approach, the Reynolds number Re is not known at the Langmuir scale. Instead, one may attempt modeling the behavior of the Langmuir flow in the limit that Re is very large. LES is an important tool for phenomenological discovery and quantification in wave-current interactions. However, it is known to be vulnerable to significant uncertainty in its subgrid-scale modeling (Sullivan and McWilliams 2010;Pearson et al. 2017). For a comprehensive review of parameterization in computational ocean modeling, see Haidvogel et al. (2017). For considerations of LES design for computational studies of global ocean circulation, see Hecht et al. (2008).
The second parameterization approach is primarily theoretical. The theory is traditionally based on the work of Craik and Leibovich (1976); Craik (1982aCraik ( , b, 1985 with later extensions by Leibovich (1977Leibovich ( , 1980Leibovich ( , 1983, Leibovich and Tandon (1993). In the Craik-Leibovich (CL) model of Langmuir circulation, wave-induced fluid motion affects the OSBL at local scales via the 'Stokes mean drift velocity' through a 'vortex force' as well as material advection. These two effects combine to produce the instability which creates the Langmuir circulation.
In WCI, the waves are propagating through the moving fluid at a speed comparable to the fluid velocity itself. This means that the wave frequency is Doppler shifted by the fluid motion. However, the wave interaction is by no means frozen into the fluid motion. Instead, the wave-current interaction (WCI) is distributed along the path of the wave through the comparably moving fluid. In particular, the Eulerian mean group velocity of the wave is defined relative to the frame moving with fluid Holm (1996), and the Eulerian mean WCI dynamics at a given fixed point in space depends on the history of wave interaction all along the entire Lagrangian path of the fluid parcel currently occupying that point. Mathematically, this implies that the description of WCI must be formulated in terms of the Eulerian mean of the pull-back of the fluid properties by the Lagrange-to-Euler map, which is assumed to be a smooth invertible map. This is a hybrid description in which the wave activity takes place in the frame of motion of the fluid.
The WCI situation is addressed directly by the Generalized Lagrangian Mean (GLM) approach formulated in Andrews and McIntyre (1978a, b). GLM generalizes the CL approach by decomposing the flow into its fast and slow components, then taking various types of time-averages, phase-averages and asymptotic approximations of the wave-current interaction dynamics at which the Rossby number is order O(1). In GLM, another dynamical variable is introduced, called the pseudomomentum, in addition to the Stokes mean drift velocity appearing in the CL approach. Thus, in GLM, the fast-slow split in time is performed at a single spatial scale. No differences in spatial scale of the waves and mean flow need to be assumed. References relevant to our purposes here are Andrews and McIntyre (1978a), Gjaja and Holm (1996), Gilbert and Vanneste (2018).
Aims of the Paper This paper aims to lay down a mathematical foundation which has the potential for both quantifying and reducing the uncertainties in the numerical simulation of ocean-atmosphere mixing layer dynamics, by developing new methods of enhanced modeling of subgrid-scale (SGS) circulation effects in the OSBL produced by wave-current interactions (WCI). Our approach is based on structure-preserving approaches in data-driven stochastic modeling for quantifying these uncertainties, combined with data assimilation methods for reducing uncertainty. Recent applications of this approach for data analysis and simulation for two-dimensional confined fluid flows are reported in Cotter et al. (2018a, b). Specifically, we lay foundations for extending the approach of Holm (2015Holm ( , 2018, Crisan et al. (2018), Cotter et al. (2018a, b) from incompressible flows in fixed domains to incompressible rotating stratified flows driven by subgrid-scale dynamics represented by stochastic processes in three dimensions. Our approach via averaged variational principles is designed to preserve the fundamental nonlinear structure of fluid dynamics. Above all, it introduces stochasticity while preserving the nature of fluid transport, the Kelvin circulation theorem and the geometric structure of fluid dynamics, including its Lie-Poisson Hamiltonian structure. In particular, our approach takes advantage of the recent developments in stochastic fluid dynamics based on geometric mechanics in Holm (2015Holm ( , 2018, Bethencourt de Léon et al. (2019), Drivas and Holm (2019) to introduce a stochastic closure procedure which preserves the geometrical structure of GLM.
The present paper also provides the derivation of a certain stochastic wave-current interaction (SWCI) model. The SWCI model is based on a stochastic closure of the well-known GLM description of the Euler-Boussinesq (EB) equations for a rotating, stratified, incompressible fluid flow. Its derivation is based on GLM averaging of a constrained Hamilton's principle for the EB equations in the Eulerian representation, leading to Euler-Poincaré variational equations for the GLM description, coupled to an Eulerian mean description of the fluctuation dynamics. This formulation is developed via a Legendre transformation into a Lie-Poisson Hamiltonian description (Holm and Kupershmidt 1983;Holm et al. 1998).
In this Hamiltonian setting, two natural stochastic closures of the GLM theory present themselves. The first closure assumes that the unknown GLM group velocity and the GLM kinematic pressure in the Hamiltonian are each temporally stochastic in the Stratonovich sense, with separate stationary spatial correlations. This closure amounts to a stochastic parameterization of the GLM group velocity and the GLM kinematic pressure whose spatial correlations must be calibrated from observed or simulated data. However, these data for the GLM stochastic closure appear to be rather inaccessible.
The elusiveness of data for the two GLM wave closures suggests the formulation of an alternative closure which directly separates the time scales of the fluid transport velocity into its slow fluid and fast wave parts. This approach is reminiscent of the introduction of the bolus velocity in the celebrated Gent-McWilliams (GM) parameterization of subgrid-scale transport (Gent 2011;Gent andMcWilliams 1990, 1996), which is generally used in computational simulations of ocean circulation.
After discovering the elusiveness of the data required in formulating the stochastic WCI closure for GLM, in the second part of the paper, we propose an alternative stochastic closure for WCI. The alternative stochastic closure proposed here is a variant of the existing theory in Gay-Balmaz and Holm (2018) of Stochastic Advection by Lie Transport (SALT) (Holm 2015(Holm , 2018Cotter et al. 2017;Crisan et al. 2018) which introduces Hamiltonian stochastic transport into the material fluid evolution as a constraint in Hamilton's variational principle for fluid dynamics. The SALT approach separates the slow and fast time scales of the fluid transport velocity into drift and stochastic parts. Implementation of this closure has already been tested in Cotter et al. (2018a, b) and found to be quite accessible for calibration by observational data from both high-resolution computational simulations. Because it deals with enhanced transport, the SALT approach is naturally compatible with formulating a data-driven stochastic version of GM parameterization of transport by unresolved time scales.
Stochastic parameterizations have been commonly used in both atmosphere and ocean sciences, ever since the break-through results of Buizza et al. (1999). Indeed, other stochastic versions of the GM already exist, as reviewed in Grooms and Kleiber (2019), and the future comparisons of these approaches with the two stochastic approaches presented here for GLM and GM are bound to be interesting. Plan In Sect. 2, we will review some background information from the GLM theory relevant to the remainder of the paper. We refer to Appendix A for the details in deriving the deterministic GLM equations for the Euler-Boussinesq equations in the Euler-Poincaré variational framework Holm et al. (1998), and the passage to the Lie-Poisson Hamiltonian side as an arena for seeking a natural stochastic closure.
Section 3 introduces stochastic closures for the GLM equations. By way of preparation, Sect. 3.1 provides a summary of the Kunita-Itô-Wentzell (KIW) theorem, which proves the key formula in stochastic transport. Section 3 then uses the KIW formula to investigate stochastic closures of the GLM Euler-Boussinesq equations due to both pressure and displacement fluctuations. Section 3 also advocates an alternative closure based on taking Stochastic Advection by Lie Transport (SALT) as a general strategy, rather than proliferating the possible sources of stochasticity for the various types of subgrid-scale physics for which our knowledge is incomplete.
In Sect. 4, the Gent-McWilliams (GM) transport scheme is reviewed and adapted to the variational SALT strategy in Gay-Balmaz and Holm (2018).
Section 5 summarizes our conclusions and outlook for open problems.

Brief Review of GLM Theory for Euler-Boussinesq Fluids
The Generalized Lagrangian Mean (GLM) theory of nonlinear waves on a Lagrangian mean flow is formulated in two consecutive papers of Andrews and McIntyre (1978a, b). The present section reviews what we shall need later from the rather complete description given in these papers. See also the textbook by Bühler (2014) for an accessible update on the GLM theory. Even now, these fundamental papers still make worthwhile reading and they are taught in many atmospheric science departments, because they represent an exceptional accomplishment in formulating averaged motion equations for fluid dynamics.

Defining Relations for Lagrangian Mean and Stokes Correction in Terms of Eulerian Mean
The GLM equations are based on defining fluid quantities at a displaced fluctuating position x ξ := x + ξ(x, t). In the GLM description, χ denotes the Eulerian mean of a fluid quantity χ = χ + χ while χ L denotes the Lagrangian mean of the same quantity, defined by Here, x ξ ≡ x + ξ(x, t) is the current position of a Lagrangian fluid trajectory whose current mean position is x. Thus, ξ(x, t) with vanishing Eulerian mean ξ = 0 denotes the fluctuating displacement of a Lagrangian particle trajectory about its current mean position x.
Remark 2.1 Fortunately, this GLM notation is also standard in the stability analysis of fluid equilibria in the Lagrangian picture. See, e.g., the classic works of Bernstein et al. (1958); Frieman and Rotenberg (1960) and Newcomb (1962). See Jeffrey and Taniuti (2000) for a collection of reprints showing applications of this approach in controlled thermonuclear fusion research. For insightful reviews, see Bernstein (1983); Chandrasekhar (1987) and, more recently, Hameiri (1998). Rather than causing confusion, this confluence of notation encourages the transfer of ideas between traditional Lagrangian stability analysis for fluids and GLM theory.
In GLM theory, the difference χ ξ − χ L = χ is called the Lagrangian disturbance of the quantity χ . One finds χ = 0, since the Eulerian mean possesses the projection property χ = χ for any quantity χ (and, in particular, it possesses that property for χ ξ ). 1 Andrews and McIntyre (1978a) show that, provided the smooth map x → x + ξ(x, t) is invertible (that is, provided the vector field ξ(x, t) generates a diffeomorphism), then the Lagrangian disturbance velocity u may be expressed in terms of ξ by Consequently, the Lagrangian disturbance velocity u is a genuine fluctuation quantity satisfying u = 0, since u ξ − u L = u ξ − u ξ = 0, by the projection property. Alternatively, u = D L ξ/Dt = 0 also follows, since the Eulerian mean commutes with D L /Dt and ξ(x, t) has mean zero.
To summarize, GLM sets u ξ (x, t) := u(x + ξ (x, t)) where x is evaluated as the current position on a Lagrangian mean path and One then defines the Lagrangian mean velocity as u ξ (x, t) = u L (x, t), where ( · ) is a time, or phase average at fixed Eulerian coordinate x.

The Pull-Back Representation of Fluctuations in Fluid Motion
Here, we briefly explain the GLM approach from the viewpoint of Cotter et al. (2017), whose multi-time homogenization analysis led to a stochastic formulation of the type proposed in the present paper. We will use the slightly expanded notation of Cotter et al. (2017) in this remark and then revert later to GLM notation. For this purpose, we will need to employ the action on functions f , k-forms α and vector fields X of smooth maps φ via pull-back, denoted φ * and defined as the composition of functional dependence from the right. For example, the expression is called the pull-back of the function f by the smooth map φ. This notation will also be applied to k-forms and vector fields. The inverse of the pull-back is called the push-forward. It is the pull-back by the inverse map. The GLM theory can be described (Cotter et al. 2017) as the Eulerian mean with respect to fast time dependence of the pull-back of the fluid properties by an evolutionary fluid flow map g t = g t/ε • g t with two time scales, one slow and one fast. This map is defined as the composition of a mean flow map g t depending on slow time t and a rapidly fluctuating flow map g t/ε associated with the evolution of the fast time scales t/ε, with ε 1. The GLM notation is recovered by defining the flow map associated with the fast scales as the (spatially) smooth invertible map with smooth inverse (i.e., a diffeomorphism, or diffeo for short) given by the sum, The full flow map is taken to be the composition of g t and g t/ε , as The Lagrangian trajectory of a fluid parcel is then given by q(x 0 , t) = g t x 0 , so that where the vector x 0 denotes the fluid label, e.g., the initial condition of a fluid parcel.
Equation (2.6) is equivalent to the displaced fluctuating position denoted as x ξ := x + ξ(x, t), in the GLM notation. That is, the rapidly fluctuating vector displacement field is defined along the slow, large-scale, resolved trajectory, q. At this point, (2.6) may be taken as exact, since it follows directly from the definition of the map ζ t/ε in (2.4). Thus, we have The tangent to the composite flow map g t in (2.5) at q(x 0 , t) along the Lagrangian trajectory (2.6) defines the Eulerian velocity vector field u, written aṡ (2.9) Differentiation of the Lagrangian trajectory (2.8) including the assumed fluctuating displacement field (2.7) yields This is equivalent to the definition of u ξ in Eq. (2.3), in the GLM notation. See Cotter et al. (2017) for more discussion of the pull-back representation of fluctuations in fluid dynamics, including results of multi-time homogenization leading to a stochastic representation of the Lagrangian trajectory in the limit that the ratio of slow and fast time scales diverges. In this case, the decomposition (2.5) becomes a composition of a stochastic map and a deterministic map.

Summary of Natural Operations on Differential k-Forms (3 k )
Differential forms are objects you can integrate. Manifolds are spaces on which the rules of calculus apply. A k-form α ∈ k on a smooth manifold M is defined by the antisymmetric wedge product of k differential basis elements, as in which the function α i 1 ...i k (x) is totally antisymmetric under exchange of any two neighboring indices. Three basic operations are commonly applied to differential forms defined on a smooth manifold, M. The three operations are: exterior derivative (d), contraction ( ) and Lie derivative (£ X ) in the direction of a vector field X . These three operations act as follows.
• Contraction (X α) with a vector field X lowers the degree: X k → k−1 . • Lie derivative (£ X α) by vector field X preserves the degree, £ X k → k .

Remark 2.2
For a k-form α ∈ k , the Lie derivative £ X α is defined geometrically by This geometric definition of the Lie derivative is called Cartan's formula.
Note that the Lie derivative commutes with the exterior derivative. That is, This useful property may be proved via a direct calculation which uses Cartan's formula and the property of the exterior derivative d that d 2 = 0.

How Pull-Back Dynamics Leads to Lie Derivatives
The pull-back φ * t of a spatially smooth flow φ t on a smooth manifold M generated by a smooth vector field X ∈ X(M) commutes with the exterior derivative d, wedge product ∧ and contraction . For an introduction to geometric fluid mechanics based on these standard concepts, see Holm (2008).
A smooth time-dependent invertible map with a smooth inverse (i.e., a diffeomorphism) φ t ∈ Diff(M) acting on a smooth manifold M may be generated by integration along the characteristic curves of a smooth vector field X (x, t) ∈ X(M) via dφ t /dt = X t • φ t . Under the action of such a smooth invertible map φ t on k-forms α, β ∈ k (M), at a point x ∈ M, the pull-back φ * t is natural for the three operations d, ∧ and . That is, in which the last equality in (2.12) is Cartan's geometric formula for the Lie derivative. The equivalence of the dynamic and geometric definitions of the Lie derivative in the last equality may be proved directly. This equivalence can be quite informative. For example, in the case α(x) = u i (x)dx i for x ∈ R 3 , i = 1, 2, 3, this equivalence implies a well-known vector calculus identity, namely The equality of these two expression, of course, yields the fundamental vector calculus identity of fluid dynamics. This calculation turns out to be the basis of the Kelvin circulation theorem.

Definition 2.3 (Pull-back and push-forward Lie derivative formulas)
The mathematical basis for analysis of fluid transport is the following textbook formula (Marsden and Hughes 1983) which relates the pull-back to the Lie derivative: (2.14) In words, the tangent to the pull-back φ * t α of a time-dependent differential k-form α ∈ k (M) by a smooth invertible flow map φ t is the pull-back φ * t of the Lie derivative of the k-form α with respect to the vector field X that generates the flow, φ t .
Likewise, for the push-forward, which is the pull-back by the inverse, (2.15) Equation (2.15) is the push-forward Lie derivative formula. Note the opposite sign from the pull-back formula in (2.14).

Definition 2.4 (Advected quantity)
An advected quantity is invariant along a flow trajectory. Thus, an advected quantity satisfies the pull-back relation which implies the transport formula, where the vector field X =φ t φ −1 t generates the flow map φ t . Equivalently, via the push-forward relation, an advected quantity satisfies (2.17)

Pull-Backs, Push-Forwards and Lie Derivatives for GLM
The GLM theory introduces a composition of maps, in which φ t,t/ε = g t/ε • g t and whose pull-back satisfies the relation,

Advection by the composition of maps
satisfies the pull-back formula for the action of the composite transformation Equivalently, the pull-back of the composition satisfies the relation Expanding out the time derivatives gives the following composite advective transport equation Recall that the pull-back, g * , is the inverse of the push-forward, g * . Hence, the pullback of the previous formula by g * implies the following version of the composite Lie transport formula, cf. Gilbert and Vanneste (2018),

GLM Advective Transport Relations for Euler-Boussinesq
For GLM, the smooth fast time flow map on the manifold M is taken to be g t/ε (M) := where γ t/ε is a smooth invertible map parameterized by the fast time, t/ε. This yields the familiar GLM fluctuation expression, g t/ε x = x + ξ (x, t/ε) = x ξ , when M is taken to be R 3 . Consequently, formula (2.18) expands out in the GLM notation, to become Thus, the expansion of the composite advective Lie transport formula (2.18) implies the following advective transport formula for a k form α, By a final transformation of variables, we will write the advection law (2.19) as This can be done by making the following chain rule calculation for the transformation of the tensor basis of α ξ (x, t) in (2.19), (2.21) Here, de(x) is the basis of the advected differential form or tensor, the quantity a(x, t) is its tensor coefficient in Eulerian coordinates and the centered dot denotes contraction of tensor indices. Equation (2.21) implies that if α ξ (x, t) is advected by u ξ , then α(x, t) will be advected by u L . This is because the fluctuating quantity α(x, t) defined above is merely a change of variables of α ξ (x, t) from x ξ to x via the chain rule. Moreover, the Eulerian mean of the relation (2.21) yields In taking this Eulerian mean, we keep in mind that x is an average quantity, so the right-hand side is already an average quantity. Thus, α satisfies α = α in (2.22) and we note that α = α L , in general, except for the case that α ξ is a scalar. The difference is that the tensor basis must be transformed to fixed Eulerian variables before applying the Eulerian time average, and a scalar function has no tensor basis.
Remark 2.5 (Road map for the remainder of the paper.) In principle, the fast-slow time mean considerations underlying GLM described above could be generalized to the class of stochastic perturbations in Holm (2015Holm ( , 2018) whose analytical properties were examined in Crisan et al. (2018) by using the method of multi-time homogenization (Gottwald and Melbourne 2013a, b) and by invoking the procedure for transition from a fast-slow description to the stochastic description for fluid dynamics developed in Cotter et al. (2017). However, instead of launching into such an investigation by starting over from a stochastic viewpoint, we will build on the deterministic theory described in the Appendix to reach the point of introducing stochastic closures for the deterministic GLM description in Sect. 3. In Sect. 3.1, we will take advantage of the result of Bethencourt de Léon et al. (2019) which proves the stochastic version of the pull-back formula (2.14) for the Lie derivative with respect to a stochastic vector field. This result will allow us to introduce a class of stochastic closures of GLM in Sect. 3.2, each of which preserves its transport structure and fits into earlier work on stochastic fluid dynamics (Holm 2015(Holm , 2018Drivas and Holm 2019). In Sects. 3.3 and 3.4, we will suggest a simplified version of one of the closure models which we expect will be convenient in potential applications for analysis of GLM investigations of WCI elsewhere.
The GM approach is adapted to the SALT framework in Sect. 4. Unlike the GLM model, which requires some development to cast it into the SALT framework, once the Gent-McWilliams (GM) approach is derived from a variational principle in Sect. 4.1, it rather easily adapts to the SALT framework for uncertainty quantification in Sect. 4.2.

GLM Circulation Transport
As an example, we shall apply the composite Lie transport formula in (2.18) to calculate the composite rate of change of the Kelvin circulation integral for the case α = u(x, t) · dx, as follows Now, if g t/ε := I d + γ t/ε , then g t/ε x = x + ξ (x, t/ε) = x ξ , and the previous formula expands out in the GLM notation, to become where J i j is the GLM fluctuating Jacobian matrix Consequently, the 1-form in the integrand of (2.24) becomes, upon assuming that whose Eulerian time average is Thus, we may conclude the following formula for the rate of change of the fast time average of the Kelvin circulation integral, (2.28) As we shall see, formula (2.28) is the basis for the definition of the pseudomomentum in the GLM theory.

GLM Scalar Advection Relations
Now that we have explained how the pull-back formula (2.18) implies the Lie derivative description of advective transport for GLM, we may return to the classic notation of GLM to discuss examples. At fixed position x the GLM velocity decomposition u ξ = u L + u is the sum of the Lagrangian mean velocity u L and the Lagrangian disturbance velocity u . Thus, and for any scalar field χ(x, t) one has, Because u L appearing in the advection operator D L /Dt = ∂ t + u L · ∇ is a mean quantity, one then finds, as expected, that the Lagrangian mean ( · ) L commutes with the original material time derivative D/Dt for a scalar function. That is, where χ = χ ξ − χ L is the Lagrangian disturbance of χ satisfying χ = 0. Hence, one finds several equivalence relations for scalars, cf. formulas (2.20) and (2.21), For example, in the Euler-Boussinesq stratified incompressible flow, consider the In this case, the buoyancy b is advected as a scalar function. That is, it satisfies Db/Dt = 0 and, by the relations (2.29), the average yields Remark 2.6 Of course, this identification is also obvious physically for scalars, since the Lagrangian mean b L and the current value b ξ refer to the same Lagrangian fluid

Mass Conservation: The GLM Continuity Equation
The instantaneous mass conservation relation D ξ (x, t) d 3 x ξ = D(x 0 )d 3 x 0 transforms into current Eulerian coordinates as follows, cf. Eq. (2.21), where one defines the Jacobian, As in the previous section, in taking the Eulerian mean of the relation D ξ J d 3 x = D d 3 x, we keep in mind that x is an average quantity, so the right-hand side is already an average quantity. Thus, D = D ξ J satisfies D = D and we note that D = D L , in general. The mean mass conservation relation for advection, D(x, t)d 3 x = D(x 0 )d 3 x 0 , then implies the continuity equation for D, upon recalling that u L is the velocity tangent to the mean Lagrangian position x. Consequently, the Lagrangian mean D ξ = D L is not the density advected in the GLM theory. Rather, it is the average density, D ξ J = D. As discussed in the previous section, except for scalars such as the buoyancy, b, this observation applies to all advected quantities. That is, the basis of any differential form or tensor field evolves under advection by the flow map, as well as its instantaneous coefficient.

Remark 2.7
For a fluid with constant density, D ξ = 1, the GLM theory gives Hence, for constant instantaneous density, the Lagrangian mean velocity u L has an as shown in Andrews and McIntyre (1978a).

Stochastic Transport via the Kunita-Itô-Wentzell Formula
The remainder of this section will introduce stochastic closure schemes for the GLM and Gent-McWilliams (GM) models of mesoscale and submesoscale transport. The key ingredient for these stochastic closure schemes will be the generalization to stochastic processes of the pull-back formula for the Lie derivative in Eq. (2.14). This stochastic generalization is given by The corresponding pull-back formula for k-forms which are spatially smooth and stochastic in time is proven in Bethencourt de Léon et al. (2019). Namely, in the standard differential notation for stochastic flows, we have where dx t is the stochastic spatially smooth vector field defined by which generates the stochastic flow φ t in Eq. (3.1).
Equation (3.1) is the Kunita-Itô-Wentzell formula (Kunita 1981(Kunita , 1984(Kunita , 1997) which determines the evolution of a k-form-valued stochastic process φ * t K . This result generalizes the classic formula for a stochastic scalar function by allowing K to be any smooth-in-space, stochastic-in-time k-form on R n . Omitting the technical regularity assumptions provided in the more detailed statement of the theorem in Bethencourt de Léon et al. (2019) for our purposes here, we now state a simplified version of the main theorem proved in that paper, as follows.
Theorem 3.1 (Kunita-Itô-Wentzell formula for k-forms, simplified version) Consider a spatially smooth k-form K (t, x) which is a semimartingale in time Formulas (3.4) and (3.6) are compact forms of the equations derived in Bethencourt de Léon et al. (2019), where these equations are written in integral notation to make the stochastic processes more explicit. However, the compact differential stochastic notation used here will suffice to explain the main ideas in the next section. For more details and proofs, see Bethencourt de Léon et al. (2019).

Stochastic Closures for GLM Approximation of the Euler-Boussinesq Equations
So far, the WCI system in the GLM equation sets (A.36) and (A.42) has not been closed. This is because the mean fluctuation quantities comprising the kinematic pressure π in (A.34) and the relative group velocity v G in (A.30) have not yet been parameterized. In this section, following Holm (2015) and Gay-Balmaz and Holm (2018), we consider two different classes of closure options for modeling these unknown quantities stochastically. Simply put, the two different classes of closure are either (1) data-driven, or (2) model-driven. In more detail, the options are: (1) apply prescribed noise which has been calibrated from observations and simulations, or (2) postulate a theoretical model for the dynamics of the noise amplitude depending on advected state variables, such as the buoyancy. In either case, the result would provide an estimate of the uncertainty in the model computations, which in turn would provide opportunities for reduction of uncertainty by using data assimilation.

Stochastic GLM Closure #1a
A very interesting approximation of the kinematic fluctuation pressure is discussed in Andrews and McIntyre (1978b); namely, Both this approximation and the relative group velocity v This observation suggests that one could close the WCI system by introducing a stochastic parameterization of these undetermined time mean correlations among the fluctuating degrees of freedom appropriate to the variable over which one is averaging. For example, the stochastic parameterization could comprise a pair of Stratonovich stochastic process, In turn, this idea suggests a new type of Hamiltonian stochastic closure which has been studied recently for fluid dynamics in Holm (2015Holm ( , 2018, Cotter et al. (2017), Crisan et al. (2018), Cotter et al. (2018a, b). It amounts to changing the GLM Hamiltonian in Eq. (A.32) into the following stochastic process, (3.8) Hamiltonian properties The stochastic GLM Euler-Boussinesq equations may be expressed in Hamiltonian Lie-Poisson matrix operator form as follows, in which the dynamics of the wave variables p and N acquires a stochastic component of its transport velocity, as (3.9) The dynamics of the material variables m j , D and b L acquires a stochastic component of its pressure force, as This stochastic pressure force does not affect the fluid circulation in Kelvin's theorem in Eq. (A.23).
In the stochastic representation of fluctuating wave effects in the GLM picture, the stochastic pressure fluctuations in (3.10) might arguably be dropped because they cause no circulation. In that case, the stochasticity of the GLM group velocity in (3.9) would coincide with the existing theory of Stochastic Advection by Lie Transport (SALT) (Holm 2015(Holm , 2018Cotter et al. 2017;Crisan et al. 2018;Cotter et al. 2018a, b) which introduces the same type of Hamiltonian stochastic transport into the material fluid evolution.

Stochastic GLM Closure #1b
Perhaps the straightest way toward the introduction of stochastic effects in WCI for use in uncertainty quantification and future data assimilation would be to consolidate the stochasticity of the GLM group velocity with the known SALT approach of adding a stochastic vector field to the Lagrangian mean transport drift velocity, u L dt, rather than proliferating the possible sources of uncertainty by making the GLM group velocity independently stochastic. In the SALT procedure, one takes the noise to be a N a=1 ζ a (x) • dW a t in which each stochastic spatial 'mode' ζ a (x) is associated with a different Brownian motion dW a t and must be calibrated, for example, by comparison of high-resolution data from either observation or computational simulation. To simplify the notation in this section, we neglect the option to include modal spatial structure in the noise by ignoring the sum over indices for the individual Brownian motions.
In the closure strategy #1b, both wave and fluid dynamics would acquire the same fluctuating component in the GLM transport velocity, as for the waves, and in which the material loop moves along stochastic Lagrangian trajectories given by the characteristics of the following stochastic vector field (3.14) Adding the stochastic vector field into (3.14) amounts to modifying the final term in the stochastic GLM Hamiltonian in Eq. (3.8), as follows, (3.15) Remark 3.2 In the class of closures #1a and #1b, with prescribed noise, it still remains to determine the set of vectors {ζ a (x t )} in the stochastic part of the Lagrangian trajectory given by dx t in Eq. (3.14). For this, it may be advisable to model the effects of wave fluctuations in the GLM equations (3.13) and (3.14) the same way as for any other high frequency transport effect in the SALT modeling approach of Cotter et al.
(2017), Holm (2015Holm ( , 2018. This approach would also simplify the calibration procedure for the correlation eigenvectors in ζ (x)•dW t , which is required in the application of SALT, because it would consolidate the stochastic effects of the wave transport with those of the material transport. Distinguishing between these two types of stochastic effects in the total transport by using observation data might be problematic, to say the least. For recent developments using the SALT approach to material transport and the description of the use of data assimilation in determining the stochastic amplitudes in two-dimensional flows, see Cotter et al. (2018a, b).

Brief Review of the Deterministic GM Approach
The SALT approach could be regarded as a data-driven stochastic version of the Gent-McWilliams (GM) parameterization of subgrid-scale transport (Gent 2011;Gent andMcWilliams 1990, 1996), which is commonly used in both ocean and atmospheric sciences. In a landmark paper, Gent and McWilliams (1990) modified passive tracer advection by adding a term meant to model eddy transport. The GM term introduced an anisotropic model of fluid transport which depends on the local gradients of the buoyancy. This term is still used today in the large majority of ocean models. Since the wave component of the GLM theory fundamentally depends on buoyancy, one can imagine that the two approaches could interact with each other synergistically. For this purpose, we will first briefly review the GM approach in the present notation. Then, we will discuss how Model 3 in Gay-Balmaz and Holm (2018) enables one to build on the GM approach and construct a stochastic closure of the motion equation in which the spatial correlations of the stochasticity depend on the quantities advected by the flow.

Geometry of the GM approach Let u(x, t) be a fluid velocity variable, and let a(x, t)
be an advected variable. The GM approach begins by introducing a modified transport equation for advection, as where L U a is the Lie derivative of the advected variable a with respect to the vector field U , and the GM model bolus velocity u * (a, a , j , a , jk ) is a prescribed vector function of a and its first two spatial derivatives. In particular, the GM model takes the advected quantity a to be the buoyancy, b, which is a scalar function.
To find the effect on the motion equation of modifying the advection law in (4.1), one may use a Lagrange multiplier to constrain Hamilton's variational principle for ideal fluids δS = 0 with S = (u, a) dt to satisfy the modified auxiliary advection equation (4.1). Before taking variations, one defines the following useful notational constructs. · , · V : V * × V → R, · , · X : X * × X → R.
3. Define the diamond operator ( ) in terms of these two pairings and the Lie derivative, as π , L δu a V = − π a , δu X , for a ∈ V , π ∈ V * and δu ∈ X(M). Thus, L δu a is the Lie derivative of the advected quantity a in the direction of the velocity variation δu.
To determine the effect on the motion equation of modifying the advection law in (4.1), we apply the auxiliary equation (4.1) as a constraint in Hamilton's principle for ideal fluids. Namely, we constrain Hamilton's variational principle δS = 0 with S = (u, a) dt to advect the quantity a by a total U = u +u * (a), by pairing Eq. (4.1) with a Lagrange multiplier, π . Thus, we set We then take variations to find, 3 δu : and manipulate further to obtain the following Euler-Poincaré equations (Holm 2015;Holm and Maddison 2019;Holm et al. 1998), (4.4) The variation δu * /δa of the prescribed bolus velocity u * (a) with respect to the advected variable a results in a differential operator in the γ -term, which arises from integration by parts in the δa-variations, contracted with the variation δ /δu i , for example, as The GM choice for u * (b) in terms of the advected buoyancy b(x, t) is linearly proportional to the local isopycnal slope s = −(∇ H b)/b z , namely, where ∇ H is the horizontal gradient. Note that divu * (b) = 0. Consequently, the pressure is determined by taking the divergence of the motion equation for incompressible flow, as usual. Upon denoting δ /δu = m, one evaluates the operator in Eq. (4.5) for constant scalar κ as The scalar advection U · ∇ δ δu part of the momentum transport L U δ δu in Eq. (4.4) appeared in Eqs. (8) and (9) of Gent and McWilliams (1996), where its magnitude was estimated as order the Rossby number, ε, so that U = u + εu * (b). Thus, according to Gent (2011), this term would make little difference in computational simulations at non-eddy-resolving resolution; so, it has never been implemented in an ocean climate computation. However, it could make a difference in ageostrophic situations, where finer resolution is required. See, e.g.   (Gent andMcWilliams 1990, 1996) in its momentum balance, energetics, Kelvin's circulation theorem and potential vorticity conservation on fluid particles. Then, we will introduce stochastic transport in the VGM setting.

Example: Euler-Boussinesq Equations
For a = (b, D) ∈ V × V * for scalar buoyancy b ∈ 0 and mass density D ∈ 3 in 3D, the diamond operations in these equations may be expressed as follows (4.8) The Lagrangian in Hamilton's principle for the Euler-Boussinesq equations is with rotation vector potential R(x) satisfying curlR(x) = 2 (x). This formula provides the variational derivatives which go into the motion equations in (4.17). For this case, the general equations in (4.4) become, e.g., for the Euler-Boussinesq equations, with a = (b, D), we choose to modify only the advected buoyancy equation, as in Gent andMcWilliams (1990, 1996). Consequently, one finds (4.10) Thus, the particle momentum density, mass density and buoyancy are all transported by the sum U = u + u * (b) of the flow velocity and the bolus velocity. Note that the quantity δ δu = δ δu · dx ⊗ d 3 x is a 1-form density, while γ ∈ V * introduced in Eq. (4.3) lies in the dual space of the advected quantity a ∈ V . The pressure p in Eq. (4.9) is a Lagrange multiplier which enforces D − 1 = 0. This constraint leads via the continuity equation to incompressibility of the augmented velocity U = u + u * (b). Because the GM choice for u * (b) in Eq. (4.6) is already divergence-free, the pressure p can then be determined from the motion equation by preservation of the divergence-free condition divu = 0, in the usual way, for appropriate boundary conditions. The divergence-free condition for a bolus velocity u * (a) depending on any other advected quantities besides the buoyancy would also be required for the theory to close.
Useful formulas for putting the general equations (4.4) into familiar calculus form for this example are, (4.11) These formulas allow one to write the VGM EB motion equation in (4.10) in standard hydrodynamics form as (4.12) with GM bolus velocity in (4.6). On the right-hand side of (4.12) three additional forces appear, all of which are bi-linear in the bolus velocity and the total circulation velocity, v. First, a Lorentz-type force appears, which is reminiscent of the Craik-Leibovich 'vortex force' in the study of Langmuir circulations. Here, the bolus velocity plays the role of the particle velocity in the Lorentz force. Second, a kinetic pressure force appears depending on higher order gradients of the buoyancy. Third, the action of the differential operator in (4.7) on the total circulation velocity, v = m/D contributes a force along the buoyancy gradient. The first two terms on the right-hand side of Eq. (4.12) can be combined as we did in Eq. (A.17) of Remark A.1 to compare the Stokes drift u S in the Craik-Leibovich (CL) theory with the pseudovelocity p/ D in GLM. Namely, we compare (4.13) This relation affords a comparison among the Kelvin circulation theorems for the CL, GLM and VGM theories. In the CL and GLM theories, the Lagrangian mean velocity transports the corresponding Kelvin circulation integrands which contain additional contributions from the fluctuations. However, in the VGM theory the circulation loop is transported by the sum U = u + u * (b) of the flow velocity and the bolus velocity.
Thus, the GM model contribution is in the Kelvin loop velocity, while the model contributions in CL and GLM are in the corresponding Kelvin circulation integrands.
Next, we survey the solution properties of the class of EB VGM equations.

Kelvin Circulation Theorem
The Kelvin circulation theorem for these equations is (4.14) Here, the circulation loop moves with the sum of the fluid velocity and the bolus velocity, U = u + u * (b).
Proof Relation (4.14) appears, upon substituting the right-hand side of the motion equation in (4.10) into the following relation The integration of the pressure gradient(s) in (4.10) around the circulation loop vanishes, and the remainder recovers Eq. (4.14).

PV Conservation
Potential vorticity (PV) is conserved, since ∂ t q + U · ∇q = 0 with q := D −1 ∇b · curlv and v = 1 D δ δu . (4.16) That is, the PV is conserved along characteristic curves (Lagrangian advection paths) of the sum of the fluid velocity and the bolus velocity.
Proof The proof can be accomplished either by using the Stokes theorem in the Kelvin theorem (4.14), or perhaps more explicitly, by first casting the EB-type equations in (4.10) into a convenient form for taking differentials, as where we have used commutation of differential d and Lie derivative L U in taking the differential of the b-equation. Now taking the differential of the (v · dx)-equation and using d(v · dx) = curlv · dS yields (4.18) Then, using the D-equation as (∂ t + L U )(Dd 3 x) = 0 yields the PV conservation equation in (4.16).

Energetics in the Hamiltonian Formulation
The Legendre transform of the constrained Lagrangian produces an extra term in the Hamiltonian with a = (D, b) for the Euler-Boussinesq Hamiltonian The semidirect-product Lie-Poisson bracket for the Euler-Boussinesq equations remains the same. Hence, the following Hamiltonian formulation of the GM transport scheme results, for the choice that the bolus velocity depends only on the advected buoyancy variable b and its derivatives, as follows (4.19) The Poisson bracket for this Hamiltonian formulation of the GM transport scheme may be expressed as For f = h G M , we find energy conservation, dh G M /dt = 0, by antisymmetry of the Lie-Poisson bracket in (4.20).

Stochastic VGM Equations
This section makes a stochastic modification of the variational , by taking the bolus velocity to be stochastic in the Stratonovich sense. Namely, Here, the differential notation dx t refers to stochastic evolution of the Lagrangian This stochastic version of the VGM transport scheme also appears in Model 3 of Gay-Balmaz and Holm (2018), although that paper did not explicitly allow the bolus velocity to depend on gradients of the advected quantities. The difference between the present scheme and the strategy of simply taking the bolus velocity to be stochastic appears in the stochastic term of the motion equation in (4.21). Otherwise, the modeling assumptions agree with those in Gent and McWilliams (1996), although they are implemented stochastically and variationally.

Stochastic Hamiltonian Formulation for the GM Transport Scheme
The Legendre transform of the constrained Lagrangian produces an extra term in the Hamiltonian The semidirect-product Lie-Poisson bracket remains the same. However, now the transport velocity vector field is stochastic, Consequently, we find the following stochastic VGM transport equations for the Euler-Boussinesq equations, when the advected quantity is chosen to be the buoyancy, b, as for Gent andMcWilliams (1990, 1996), As in the deterministic VGM transport scheme in the previous section, the constraint D − 1 = 0 enforces div(u dt + u * (b) • dW t ) = 0 in the stochastic case, as well. This implies divu = 0, since we already have divu * (b) = 0 by (4.6). This result makes the determination of the pressure p systematic and straightforward for stochastic VGM, as well.

Remark 4.2 (Noether's theorem)
The presence of explicit time and space dependence in the stochastic part of the Hamiltonian dh(m, a) in (4.22) precludes conservation of energy and momentum in the VGM transport scheme, respectively. However, the Kelvin circulation theorem in Eq. (4.14) and the PV conservation in Eq. (4.16) both still persist in the presence of the stochastic transport, modulo replacement of the deterministic advective transport velocity by its stochastic counterpart. These two conservation laws result from Noether's theorem for invariance under relabeling of Lagrangian particles and conservation of advected quantities along Lagrangian particle trajectories. To the extent that the initial spatial distributions of the advected quantities reduce the relabeling symmetry to the isotropy subgroup of the diffeomorphisms which preserves the initial distributions of advected quantities, the Kelvin circulation integral is not preserved in time. The Kelvin-Noether theorem for the Euler-Poincaré equation developed in Holm et al. (1998) represents the evolution of the Kelvin circulation resulting from breaking the relabeling symmetry. This is the converse of the Noether theorem for fluid dynamics with advected quantities. The Legendre transform of the constrained Lagrangian (4.9) in Hamilton's principle for the Euler-Boussinesq equations, for example, produces the extra term in the Hamiltonian in (4.22). Thus, the additional transport velocity introduced in the advective constraint on the variations in Hamilton's principle (4.9) is responsible for the choice of the Hamiltonian in Eq. (4.22). This extra transport velocity is also responsible for the additional forcing of the circulation in Eq. (4.14).
(1996) for wave packets, in which GLM may be closed at various asymptotic orders. These basic results were reviewed from the viewpoint of geometric mechanics, particularly via the Euler-Poincaré formulation of Lagrangian reduction by the symmetry of particle relabeling for continuum mechanics in Holm et al. (1998). In the geometric mechanics framework, the Lie-Poisson structure of GLM emerges as a classical Hamiltonian field theory with particle relabeling symmetry. However, the theory is not closed until further assumptions have been made about the group velocity of the waves and the solution for the pressure due to fluctuations.
Several closure procedures have been introduced previously. In the WKB representation of WCI interaction in Euler-Boussinesq fluids (Gjaja and Holm 1996), the closure was supplied at various asymptotic orders via the dispersion relation and phaseaveraged pressure contributions of the waves. By applying slow manifold reduction (MacKay 2004) to dynamics in the space of loops, a broader class of variational nonlinear WKB closures for WCI in plasmas was derived in Burby and Ruiz (2019), and expressed in the standard Eulerian frame, rather than the displaced GLM Eulerian frame. In previous work, similar ideas were applied in both turbulence modeling (Holm and Tronci 2012a, b) and in shape analysis (Bruveris et al. 2011). In earlier work on fluid turbulence modeling, a similar type of closure was based on invoking the Taylor hypothesis, that fluctuating quantities would be carried along in the fluid, e.g., Holm (2002a, b).
In the geometric mechanics setting here, we have added considerations of stochastic modeling of the indeterminate quantities in GLM, based on recent advances in stochastic transport (Bethencourt de Léon et al. 2019), stochastic variational principles and the Hamiltonian formulations of their results (Holm 2015(Holm , 2018Gay-Balmaz and Holm 2018). This variational stochastic formulation seems to promise many future opportunities for the combination of stochastic variational modeling and data assimilation, which in this setting has already had promising results, both in mathematical analysis ) and in uncertainty quantification (Cotter et al. 2018a, b). In particular, the analysis in Crisan et al. (2018) showed that the presence of the stochastic transport in Euler's fluid equation preserves its analytical properties in the deterministic case. Namely, the stochastic transport version of Euler's fluid equation has local-intime existence and uniqueness, while also satisfying the Beale-Kato-Majda criterion for blow-up of the solution.
Section 3 considers data-driven and model-driven classes of stochastic closure options for GLM. The purpose of these stochastic closures would be to provide an estimate of the uncertainty in the model computations, which in turn would provide opportunities for reduction of uncertainty by using data assimilation. The data-driven closure option invokes the SALT method of Cotter et al. (2018a, b), while the modeldriven closure option invokes the familiar Gent-McWilliams approach, as generalized to the stochastic case in Gay-Balmaz and Holm (2018).
Because of the close relation of wave propagation to buoyancy dynamics, we chose the stochastic Gent-McWilliams approach in Sect. 3 to illustrate the example of stochastic transport in the Euler-Boussinesq equations, rather than taking the full GLM equations. One may regard the GM discussion as a first step toward making the buoyancy dynamics in the wave components of GLM fully stochastic, in the sense of making the noise-mean flow interaction dynamical.
The GM step also opens the opportunity to quantify the uncertainty of the GM transport scheme, itself. The GM scheme is widely used in computational ocean science (Gent 2011). Here, we note that applying either the deterministic or stochastic GM advective transport scheme in the buoyancy equation in computational simulations while neglecting both its contributions in the motion equation and in the modified incompressibility condition imposed via the continuity equation could be expected to produce errors in the momentum balance. In turn, these errors will cascade into errors in the circulation and PV transport. It would be interesting to quantify the effects of those types of uncertainties, as well.
Finally, the introduction of stochastic channels into WCI may provide a means of parameterizing wave breaking. For example, in the GLM setting, under wave forcing at the surface, one could introduce a jump process which would stochastically transfer a certain amount of pseudomomentum p to material momentum m while keeping the sum of the two momenta p + m constant at a given point. For example, this sort of bursting event in momentum transfer could be triggered by a threshold in wavenumber steepness ( z · ∇k) 2 /| z × ∇k| 2 > 1, where k = |k| = |p|/N . GLM wave breaking has not been widely considered, and this approach to it has not been tried in applications yet. Likewise, in the stochastic GM setting, since the bolus velocity u * (b) figures dynamically in both the buoyancy equation and the momentum equation in (4.23), one might consider jump processes which induce stochastic impulses into the momentum balance which are triggered by the steepness of the local isopycnal slope s = −(∇ H b)/b z . Thus, the loss of momentum conservation due to the spatial dependence of the noise would be regarded as stochastic forcing. These are only preliminary thoughts which must continue to develop and be investigated elsewhere.
where (x ξ ) is a potential for external or centrifugal forces. If desired, the rotation frequency can be allowed to depend on position along the fluctuating path x ξ as 2 (x ξ ) = (curl R) ξ . The corresponding rotation potential is decomposed in standard GLM fashion as Upon substituting the defining relation into (A.2), the definition of D in (2.30) allows one to write the corresponding Eulerian mean expression of the averaged Lagrangian for the Euler-Boussinesq stratified fluid as Here, the last term introduces the Lagrange multiplier to impose the constraint that the fluctuation velocity u must satisfy its definition via the material derivative of the fluctuation vector displacement field ξ in Eq. (A.2). The relative buoyancy defined by the mass density ratio b ξ = (ρ ref − ρ ξ )/ρ ref is advected as a scalar in the Boussinesq approximation, in the case that D ξ = 1. Most of the important properties of the GLM equations are discussed in Andrews and McIntyre (1978a, b). Many of these properties arise from general mathematical structures that are shared by all exact nonlinear ideal fluid theories; namely, as an Euler-Poincaré (EP) equation (Holm et al. 1998), which is expressed in terms of variational derivatives of an averaged Lagrangian, (u L , D, b L ) and obtained from Hamilton's principle for the Lagrangian mean variables, See Holm et al. (1998) for an exposition of the mathematical structures which arise in the EP theory of ideal fluids which possess advected quantities, such as buoyancy, entropy and magnetic field. In Eq. (A.4), for example, the right-hand side is the usual baroclinic source term.
In particular, the EP equation (A.4) for GLM implies the following Kelvin circulation theorem for the GLM Euler-Boussinesq flow, for any closed loop γ L (t) moving with the Lagrangian mean flow velocity u L . The proof of (A.5) follows immediately by noting that .6) and that the GLM EP equation (A.4) may be written as after using the advection law for D in Eq. (2.32).

Variational Derivatives and the EP Equation for GLM Euler-Boussinesq Stratified Fluid
The mean Lagrangian 3) has been derived via a straight transcription from the standard Lagrangian for Euler-Boussinesq fluids into the GLM formalism, followed by taking the Eulerian mean. Its variational derivatives are given by The last term in the k equation arises from a spatial integration by parts of the variation p ξ δJ in which δJ = K j k (∂ δξ k /∂ x j ) with cofactor Thus, the variations in the fluctuating quantities imply the following quasilinear equations with vanishing mean, The variations with respect to δu L and δu each provides a momentum map. Combining them yields, (A.10) in which the last step defines the pseudomomentum density, p. The average of a combination of the second and third equation in (A.9) will provide the dynamical equation we need for the pseudomomentum density in order to close the equations. We may also refer to the ratio v := p/ D := − (u j + R j )∇ξ j (A.11) as the pseudovelocity, v, see formula (2.28). The Boussinesq potential B arising in (A.8) under the variation of with respect to D is defined by (A.13) and, finally, p L = p ξ is the Lagrangian mean pressure.
Upon substituting these variational derivatives into the Euler-Poincaré (EP) equation (A.4), one finds the following GLM motion equation governing u L for a stratified Boussinesq fluid in Cartesian coordinates, One could also write this equation to mimic a 'vortex force' in Lorentz form E+u L ×B as (A.15) For convenience, the equations for the advected quantities b L and D are recalled from above as Remark A.1 (Comparison of GLM pseudomomentum dynamics with the Craik-Leibovich theory) Without the 'E-field' term on its right side, Eq. (A.15) would seem to have the same form as the Craik-Leibovich theory, except that the Stokes mean drift velocityū S would have been replaced by the pseudovelocity v. Formally, then, the GLM Euler-Boussinesq stratified fluid equations might appear to comprise a dynamical version of the Craik-Leibovich theory. However, the pseudovelocity v is by no means the same as the Stokes mean drift velocity, u S . In fact, their difference has nonzero circulation. This is because the pseudovelocity, v = p/ D, and the Stokes mean drift velocity, u S , are complementary quantities in the Eulerian mean of L ξ (u · dx), which is the Lie derivative of the fluctuating circulation 1-form u · dx with respect to the fluctuation vector field, ξ . Namely, So the two 'velocities' meet here in the Lie derivative. They are so different that their difference means something. The Stokes mean drift velocity, u S , is the rate of distortion of the fluctuating velocity covector by the fluctuating disturbance in the Lagrangian path away from its mean, as if the covector were an array of scalars. The pseudovelocity v is (minus) the corresponding rate of distortion of its covector basis. The place where all this comes together is in the GLM Kelvin's theorem when we bring in the Eulerian mean velocity u E to transform from Lagrangian mean to Eulerian mean quantities in the integrand as For further discussion of the geometric and Hamiltonian properties of the Craik-Leibovich theory, see Holm (1996).

Remark A.2
We still need an equation for the pseudomomentum density p in Eq. (A.10) in order to close the GLM Euler-Boussinesq motion equation in (A.14). However, before deriving that equation, let us make a few remarks about the properties of the (as yet unclosed) GLM equations for the Euler-Boussinesq stratified fluid which have been obtained, so far.

Relation to the EP Kelvin Circulation Theorem for GLM Boussinesq Stratified Fluid
The GLM average of Kelvin's circulation integral is defined as, where the contourγ L (t) moves with velocity u L , since it follows the fluid parcels as the average is taken. Thus, the Lagrangian mean leaves invariant the form of the Kelvin integral, while averaging the velocity of its contour. In addition, the pseudovelocity covector v defined in (A.10) appears in the integrand of the GLM averaged Kelvin integral I (t).
The time derivative of the GLM averaged Kelvin circulation integral is, cf. formula (2.28), where curl R L (x) = 2 (x). The combination of terms in the integrand defines the transport structure of the GLM theory under the Lie derivative L u L along the mean velocity vector, u L . From the GLM motion equation (A.14) one now finds the GLM Kelvin circulation theorem for Boussinesq incompressible flow, Evaluating this for the GLM Boussinesq stratified fluid with given in (A.3) yields,

Local Potential Vorticity Conservation for GLM Boussinesq Stratified Fluid
Invariance of the Lagrangian under diffeomorphisms (interpreted physically as Lagrangian particle relabeling) implies the local conservation law for EP potential vorticity, For the GLM case, the potential vorticity is given explicitly as The EP framework explains the relation of the potential vorticity to the Kelvin circulation theorem. However, there remains the question of the evolution of the pseudovelocity, v.

Fluctuation Equations
Hamilton's principle for the Lagrangian mean variables {u L , D, b L } has already been calculated in Eq. (A.8). We now apply Hamilton's principle for the fluctuation variable ξ k using the original Lagrangian (u ξ , D ξ , b ξ , ξ, ∂ t ξ) in Eq. (A.1).
The result for the momentum density k canonically conjugate to ξ k is Wave action density To introduce the wave action density N and explain how it is related to the GLM pseudomomentum density, p, we take the Eulerian mean of the following pre-canonical transformation, If ξ and π are averaged over a phase parameter, φ, we may write the phase-averaged differential relation as where the wavevector k is defined by dφ = ∇φ · dx = k · dx and the wave action density N is given by Thus, the wave action density N = − k ∂ φ ξ k is related to the GLM pseudomomentum by p = N k.
For the WKB wavepacket ξ(x, t) = 1 2 (a(x, t)e iφ(x,t)/ + a * (x, t)e −iφ(x,t)/ ), one finds the formula for constant Coriolis parameter 2 , (Gjaja and Holm 1996), in which the quantity is the Doppler-shifted wave frequency. As a result of the symmetry under translations in φ induced by phase-averaging the Lagrangian, the corresponding Euler-Lagrange equation implies the conservation law, upon using the variational derivatives in equation (A.8). Andrews and McIntyre (1978b) obtain the same conservation law by directly manipulating the GLM motion equation. This is also Noether's theorem for symmetry of the Lagrangian under phase shifts. For more discussion from a variational viewpoint in the case that the fluctuations are single-frequency wave packets with slowly varying envelopes, see also Gjaja and Holm (1996). Of course, Noether's theorem always applies in averaging Hamilton's principle, since such averaging always produces a continuous symmetry of the Lagrangian. In general, Noether's theorem implies the following about the relation of averaging to local conservation laws, (Hayes 1970;Andrews and McIntyre 1978b;Holm 2002a, b).
Lemma A.5 When Lagrangian averaging introduces an ignorable coordinate in fluid dynamics, the average of the corresponding canonically conjugate momentum is locally conserved; that is, the corresponding quantity is conserved in a shifted frame of motion relative to Lagrangian fluid parcels. In this case, the locally conserved quantity is the wave action density N in (A.28), which is the phase-averaged quantity (momentum map) whose canonical Poisson bracket generates phase shifts. The spatial integral over the domain of. flow D N d 3 x is conserved globally, for appropriate boundary conditions.
We do not vary H with respect to the parameters ω, k and v G . The term (p − N k) · v G vanishes for arbitrary v G , as a consequence of the variation in u L . Moreover, the expected 'wave conservation relation' ∂ t k = −∇ω will follow as a result of the other dynamical equations. We note that the constraints on the averaged Lagrangian will still apply for the Hamiltonian, since they are not Legendre transformed. We may now compute the variations of the Hamiltonian as where tot is given by tot = δ H δ D = p L + gzb L + L (x) − 1 2 |u | 2 + u · R =: L + π . (A.34) Vanishing of the other variations of the averaged Lagrangian in (A.32) still enforces the constraints (A.9) since the corresponding variables were not Legendre transformed.
Wave component We now write the equations of motion for the pseudomomentum density and wave action density in Lie-Poisson form, following the lead of Gjaja and Holm (1996) and we can may choose v j G = ( p ξ K j i ∂ φ ξ i i ) to agree with the definition in (A.30).
Remark A.6 (Wave conservation) Note that Eq. (A.35) and the relation p = N k imply the wave conservation relation ∂ t k = −∇ω.
Lie-Poisson Hamiltonian structure The wave field's semidirect-product Lie-Poisson Hamiltonian structure may be revealed by its Poisson operator, given in matrix form by . (A.36) Expanding out the matrix product yields the Lie-Poisson bracket between two functionals F and H as, The Lie-Poisson bracket in Eq. (A.37) is defined on the dual of the semidirectproduct Lie algebra X 0 of vector fields X ∈ X(M) and functions f ∈ 0 (M) on the domain of flow, M. The corresponding Lie algebra commutator is given by where [X , X ] is the commutator of vector fields and X ( f ) is the Lie derivative of vector fields acting on functions. The dual coordinates are: the pseudomomentum 1form density, p = p · dx ⊗ d 3 x, dual to vector fields; and the wave action density, N d 3 x, dual to functions. Thus, the Lie-Poisson bracket in Eq. (A.37) may be written as (A.39) In other standard notation (Holm et al. 1998), this is (A.40) The corresponding forms of their equations of motion in (A.35) when written in terms of Lie derivatives are The corresponding forms of the fluid equations in (A.42) may then be written in terms of Lie derivatives are Thus, the particle momentum density, mass density and buoyancy are all transported by the same Lagrangian mean velocity. The geometric similarities pervading the equations for the dynamics of the wave and material components of the WCI system argues that it should be treated as a twofluid system, e.g., as for HeII. If so, then one should note that, just as for HeII, the two fluids interpenetrate one another, since the wave and material properties are transported at different velocities. The material component of the GLM fluid is transported at the Lagrangian mean velocity, u L , while the wave component of the GLM fluid is transported at the sum of velocities, u L + v G .
The Lie-Poisson bracket for the WCI system is the sum of two Lie-Poisson brackets. That is, the Lie-Poisson bracket for WWCI is dual to the direct-sum Lie algebra p will also be affected by the wave pseudomomentum p equation, via a type of Lorentz force reminiscent of the 'vortex force' in the Craik-Leibovich theory, except that the Stokes mean drift velocity u S in the CL theory will be replaced by the pseudovelocity v = p/ D in Eq. (A.11). The corresponding Lie-Poisson structure can be obtained by a linear change of variables.