3.1. Mass conservation

The Lagrangian average used in the formulation above, with a non-divergent $\boldsymbol{u}^{\scriptscriptstyle{L}}$ in (2.9), is not the generalised Lagrangian mean (GLM) of Andrews & McIntyre (Reference Andrews and McIntyre1978): in GLM, $\boldsymbol{u}^{\scriptscriptstyle{L}}$ has non-zero divergence (McIntyre Reference McIntyre1988; Ardhuin et al. Reference Ardhuin, Rascle and Belibassakis2008). Instead, we use the Lagrangian mean of Soward & Roberts (Reference Soward and Roberts2010) (written as glm in Soward & Roberts (Reference Soward and Roberts2010)) which has a non-divergent Lagrangian mean velocity (Vanneste & Young Reference Vanneste and Young2022).

(Mass is conserved by GLM: with non-zero ${\boldsymbol{\nabla}} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}}$ there is a Lagrangian density variable – denoted $\tilde \rho$ in Andrews & McIntyre (Reference Andrews and McIntyre1978) – ensuring mass conservation. Our decision to use the Lagrangian mean of Soward & Roberts (Reference Soward and Roberts2010), with exact ${\boldsymbol{\nabla}} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}}=0$ , avoids the additional $\rho$ -variable.)

3.2. Momentum conservation

We start with the vertically integrated pseudomomentum equation obtained by forming $\iint (3.1) \boldsymbol{k} \, \mathrm{d} \boldsymbol{k}$ :

(3.3) \begin{align} \partial _t \iint \mathcal{N}(\boldsymbol{x},\boldsymbol{k}) k_{\alpha }\, \mathrm{d} \boldsymbol{k} +\partial _{\beta } \iint \mathcal{C}_{\beta } k_{\alpha }\mathcal{N} \mathrm{d} \boldsymbol{k} + \iint \mathcal{N} \partial _{\alpha }\varOmega \, \mathrm{d} \boldsymbol{k} = \iint k_{\alpha }\mathcal{N}^\circ \, \mathrm{d} \boldsymbol{k} . \end{align}

In (3.3) the Greek indices $\alpha$ and $\beta$ run from $1$ to $2$ and $\partial _{\alpha }$ and $\partial _{\beta }$ denote derivatives with respect to $x_\alpha$ , e.g. $\partial _{\alpha } =\partial /\partial x_\alpha$ . The group velocity, $\mathcal{C}_\beta$ above, is defined in (2.14). With pseudomomentum $\boldsymbol{\mathsf{p}}(\boldsymbol{x},z)$ defined in (2.10), the first term in (3.3) is $\partial _t$ applied to the total vertical integral of $\boldsymbol{\mathsf{p}}(\boldsymbol{x},z)$ . The third term in (3.3), involving $\partial _{\alpha } \varOmega$ , is transformed by switching the order of $\boldsymbol{k}$ - and $z$ -integrals; see Appendix A for details. These manipulations bring (3.3) to the form

(3.4) \begin{align} \partial _t \int \! {\mathsf{p}}_\alpha \, \mathrm{d} z + \partial _{\beta } \iint \mathcal{C}_{\beta } k_{\alpha } \mathcal{N} \, \mathrm{d} \boldsymbol{k} + \int {\mathsf{p}}_{\beta }\partial _{\alpha } u^{\scriptscriptstyle{L}}_{\beta } \, \mathrm{d} z + \iint \frac {\partial \varOmega }{\partial d} \mathcal{N} \mathrm{d} \boldsymbol{k}\, \partial _{\alpha } d = \iint k_{\alpha }\mathcal{N}^\circ \, \mathrm{d} \boldsymbol{k} . \end{align}

Above and hereafter, we lighten notation by suppressing the limits $-d(\boldsymbol{x})$ and $0$ on vertical integrals. In writing the derivative $\partial \varOmega /\partial d$ we abuse notation slightly by treating $\varOmega$ as a function of $d$ rather than $\boldsymbol{x}$ as in (2.2). The interpretation is clear from the explicit dependence of $\sigma$ and $\boldsymbol{\mathcal{U}}$ (via $Q$ ) on $d$ in (2.3) and (2.4)–(2.5). We give an explicit expression for $\partial \varOmega /\partial d$ in (A3) and (A8).

Turning to the Craik–Leibovich equation, the vertical integral of the horizontal components of (2.18) is

(3.5) \begin{align} \partial _t \! \int \!u^{\scriptscriptstyle{L}}_{\alpha } \, \mathrm{d} z -& \partial _t \!\int \!{\mathsf{p}}_\alpha \mathrm{d} z - \!\int {\mathsf{p}}_{\beta } \partial _{\alpha } u^{\scriptscriptstyle{L}}_{\beta } \, \mathrm{d} z + \partial _{\beta }\!\int \!u^{\scriptscriptstyle{L}}_\alpha u^{\scriptscriptstyle{L}}_\beta - {\mathsf{p}}_\alpha u^{\scriptscriptstyle{L}}_\beta \, \mathrm{d} z \nonumber \\ &\quad + \Big (\int \!\boldsymbol{f} {\times } \boldsymbol{u}^{\scriptscriptstyle{L}} \, \mathrm{d} z \Big )_{\alpha }+ \partial _{\alpha } \!\int \tilde \pi \, \mathrm{d} z - [ \tilde \pi ]_{-d} \, \partial _{\alpha } d = \int \!u^\circ _\alpha \, \mathrm{d} z , \end{align}

where $[\,\boldsymbol{\cdot }\,]_{-d}$ denotes evaluation at $z=-d$ , e.g. $ [ \tilde \pi ]_{-d} = \tilde \pi (\boldsymbol{x}, -d(\boldsymbol{x}))$ . Terms stemming from the commutation of $\partial _\beta$ with the vertical integral vanish on using the boundary conditions (2.12) and (2.13). The sole remaining boundary contribution is the topographic form stress $[ \tilde \pi ]_{-d} \, \partial _{\alpha } d$ in (3.5).

Adding (3.4) and (3.5) cancels the terms

(3.6) \begin{align} \partial _t \int \!\! {\mathsf{p}}_\alpha \, \mathrm{d} z + \int {\mathsf{p}}_{\beta }\partial _{\alpha }u^{\scriptscriptstyle{L}}_{\beta } \, \mathrm{d} z , \end{align}

leaving

(3.7) \begin{align} \partial _t &\!\int \!u^{\scriptscriptstyle{L}}_{\alpha } \, \mathrm{d} z + \partial _{\beta }\left (\int \!u^{\scriptscriptstyle{L}}_\alpha u^{\scriptscriptstyle{L}}_\beta - {\mathsf{p}}_\alpha u^{\scriptscriptstyle{L}}_\beta \,\mathrm{d} z + \iint k_\alpha \mathcal{C}_\beta \mathcal{N} \mathrm{d} \boldsymbol{k} \right ) + \Big (\int \!\boldsymbol{f} {\times } \boldsymbol{u}^{\scriptscriptstyle{L}}\mathrm{d} z \Big )_{\alpha } \nonumber \\ & + \partial _{\alpha } \!\int \tilde \pi \, \mathrm{d} z + \bigg ( \iint \frac {\partial \varOmega }{\partial d} \mathcal{N} \mathrm{d} \boldsymbol{k} \, - [\tilde \pi ]_{-d} \bigg ) \partial _{\alpha } d = \int \! u^\circ _\alpha \mathrm{d} z + \iint k_\alpha \mathcal{N}^\circ \, \mathrm{d} \boldsymbol{k} . \end{align}

The result

(3.8) \begin{align} \int {\mathsf{p}}_\alpha u^{\scriptscriptstyle{L}}_\beta \, \mathrm{d} z = \iint k_\alpha \, \boldsymbol{\mathcal{U}}_\beta \, \mathcal{N} \mathrm{d} \boldsymbol{k} , \end{align}

deduced from the definition of the Doppler velocity $\boldsymbol{\mathcal{U}}$ in (2.4) and the definition of pseudomomentum $\boldsymbol{\mathsf{p}}$ in (2.10), produces a further partial cancellation in (3.7):

(3.9) \begin{align} \iint k_\alpha \mathcal{C}_\beta \mathcal{N} \mathrm{d} \boldsymbol{k} - \int {\mathsf{p}}_\alpha u^{\scriptscriptstyle{L}}_\beta \, \mathrm{d} z &=\iint \mathcal{G} k_\alpha k_\beta \, \mathcal{N} \, \mathrm{d} \boldsymbol{k} , \end{align}

where $\mathcal{G}(\boldsymbol{x},\boldsymbol{k})$ is defined in (2.15). Thus, (3.7) becomes

(3.10) \begin{align} \partial _t & \!\int \!u^{\scriptscriptstyle{L}}_{\alpha } \, \mathrm{d} z + \partial _{\beta }\left (\int \!u^{\scriptscriptstyle{L}}_\alpha u^{\scriptscriptstyle{L}}_\beta \, \mathrm{d} z + \iint \mathcal{G} \, k_\alpha k_\beta \mathcal{N} \mathrm{d} \boldsymbol{k} \right ) + \Big (\int \!\boldsymbol{f} {\times } \boldsymbol{u}^{\scriptscriptstyle{L}}\mathrm{d} z \Big )_{\alpha }\nonumber \\ & + \partial _{\alpha } \!\int \tilde \pi \, \mathrm{d} z +\bigg ( \iint \frac {\partial \varOmega }{\partial d} \mathcal{N} \mathrm{d} \boldsymbol{k} \, - [ \tilde \pi ]_{-d} \bigg ) \partial _{\alpha } d = \int \! u^\circ _\alpha \mathrm{d} z + \iint k_\alpha \mathcal{N}^\circ \, \mathrm{d} \boldsymbol{k} . \end{align}

This is momentum conservation in the absence of the Coriolis term, depth variation $\partial _{\alpha } d$ , and of forcing and dissipation on the right-hand side.

Although the pseudomomentum $\boldsymbol{\mathsf{p}}$ appears throughout the manipulations above, remarkable cancellations ensure that $\boldsymbol{\mathsf{p}}$ does not feature in the final form (3.10). This is a manifestation of McIntyre’s (Reference McIntyre1981) ‘wave momentum myth’: the total momentum density, namely

(3.11) \begin{align} \!\int \boldsymbol{u}^{\scriptscriptstyle{L}} \, \mathrm{d} z , \end{align}

contains only the Lagrangian mean velocity and the direct effect of the waves in (3.10) is flux of momentum via the radiation stress tensor

(3.12) \begin{align} \iint \mathcal{G} \, k_\alpha k_\beta \mathcal{N} \mathrm{d} \boldsymbol{k} , \end{align}

where $\mathcal{G} (\boldsymbol{x},\boldsymbol{k})$ is defined in (2.15).

As noted by Dewar (Reference Dewar1970) we can regard (3.11) as the sum of the Lagrangian mean momentum of the current

(3.13) \begin{align} \int \boldsymbol{u}^{\scriptscriptstyle{L}} - \boldsymbol{\mathsf{p}} \, \mathrm{d} z \end{align}

and the wave pseudomomentum $\int \boldsymbol{\mathsf{p}} \, \mathrm{d} z$ (often termed wave momentum). This is consistent with a geometric interpretation of the Lagrangian mean momentum (Gilbert & Vanneste Reference Gilbert and Vanneste2018, Reference Gilbert and Vanneste2025). The total momentum (3.11) is the conserved quantity associated with translational symmetry of the system.

In this section, for simplicity, we have focused on the depth-integrated momentum. By forming $\iint (2.1) Q \boldsymbol{k} \, \mathrm{d} \boldsymbol{k}$ and combining with the Craik–Leibovich equation (2.18), it is possible to derive an equation for the total momentum $\boldsymbol{u}^{\scriptscriptstyle{L}}$ at each $z$ . This equation reduces to (3.10) when integrated with respect to $z$ . It provides an alternative to (2.18) in which wave forcing appears through radiation stresses with vertical as well as horizontal components.

3.3. Energy conservation

Following Bretherton & Garrett (Reference Bretherton and Garrett1968), we expect that the wave-energy density is

(3.14) \begin{align} \mathcal{E}_{w}(\boldsymbol{x},\boldsymbol{k})=\sigma (\boldsymbol{x},\kappa ) \, \mathcal{N}(\boldsymbol{x},\boldsymbol{k}) , \end{align}

where $\sigma (\boldsymbol{x},\kappa )$ is the intrinsic frequency in (2.3). We also expect that the total energy density,

(3.15) \begin{align} E(\boldsymbol{x}) = \int \tfrac {1}{2} |\boldsymbol{u}^{\scriptscriptstyle{L}} |^2 \, \mathrm{d} z + \iint \mathcal{E}_{w}(\boldsymbol{x},\boldsymbol{k}) \, \mathrm{d} \boldsymbol{k} , \end{align}

satisfies a physical space conservation law

(3.16) \begin{align}\partial _t E + {\boldsymbol{\nabla}} _{\!\boldsymbol{x}} \!\,{\boldsymbol{\cdot}}\, \!\boldsymbol{J} = \, {E}^\circ . \end{align}

The horizontal vector $\boldsymbol{J}(\boldsymbol{x})$ above is the energy flux.

To determine $\boldsymbol{J}$ and ${E}^\circ$ in (3.16), start with the depth integral of $\boldsymbol{u}^{\scriptscriptstyle{L}} \,{\boldsymbol{\cdot}}\, (2.8)$ and use the vertical boundary conditions (2.12) and (2.13) to find

(3.17) \begin{align} \partial _t \int \tfrac {1}{2} |\boldsymbol{u}^{\scriptscriptstyle{L}}|^2 \, \mathrm{d} z +{\boldsymbol{\nabla}} _{\!\boldsymbol{x}} \,{\boldsymbol{\cdot}}\, \int \boldsymbol{u}^{\scriptscriptstyle{L}} \pi \, \mathrm{d} z= \int \boldsymbol{u}^{\scriptscriptstyle{L}}\,{\boldsymbol{\cdot}}\, \partial _t \boldsymbol{\mathsf{p}}\, \mathrm{d} z + \int \boldsymbol{u}^{\scriptscriptstyle{L}}\,{\boldsymbol{\cdot}}\, \boldsymbol{u}^\circ \mathrm{d} z . \end{align}

Equation (3.17) is the energy equation for the currents.

Although the wave energy $\mathcal{E}_{w}$ in (3.14) is related to action density $\mathcal{N}(\boldsymbol{x},\boldsymbol{k})$ by multiplication with the intrinsic frequency $\sigma (\boldsymbol{x},\boldsymbol{k})$ , the wave energy conservation law is obtained with $\iint (3.1) \varOmega \, \mathrm{d} \boldsymbol{k}$ , i.e. multiplication of the action conservation equation by the total frequency $\varOmega(\boldsymbol{x},\boldsymbol{k})=\sigma (\boldsymbol{x},\kappa ) + \boldsymbol{k}\,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}}(\boldsymbol{x},\kappa )$ . This produces

(3.18) \begin{align} \partial _t \iint \mathcal{E}_{w} \, \mathrm{d} \boldsymbol{k} + \iint (\boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}}) \partial _t \mathcal{N} \, \mathrm{d} \boldsymbol{k} + {\boldsymbol{\nabla}} _{\!\boldsymbol{x}} \,{\boldsymbol{\cdot}}\, \iint \boldsymbol{\mathcal{C}} \varOmega \mathcal{N} \mathrm{d} \boldsymbol{k} = \iint \varOmega \mathcal{N}^\circ \mathrm{d} \boldsymbol{k} , \end{align}

where $\boldsymbol{\mathcal{C}}$ is the group velocity in (2.14). The key identity

(3.19) \begin{align} \iint ( \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}} ) \, \partial _t \mathcal{N} \, \mathrm{d} \boldsymbol{k} = \int \!\boldsymbol{u}^{\scriptscriptstyle{L}} \! \,{\boldsymbol{\cdot}}\, \! \iint Q \boldsymbol{k} \, \partial _t \mathcal{N} \, \mathrm{d} \boldsymbol{k} \, \mathrm{d} z = \int \!\boldsymbol{u}^{\scriptscriptstyle{L}} \!\,{\boldsymbol{\cdot}}\, \partial _t \boldsymbol{\mathsf{p}} \, \mathrm{d} z \end{align}

then puts (3.18) in the form

(3.20) \begin{align} \partial _t \iint \mathcal{E}_{w}\, \mathrm{d} \boldsymbol{k} + \int \boldsymbol{u}^{\scriptscriptstyle{L}} \,{\boldsymbol{\cdot}}\, \partial _t \boldsymbol{\mathsf{p}} \, \mathrm{d} z + {\boldsymbol{\nabla}} _{\!\boldsymbol{x}} \,{\boldsymbol{\cdot}}\, \iint \boldsymbol{\mathcal{C}} \varOmega \mathcal{N} \mathrm{d} \boldsymbol{k} = \iint \varOmega \mathcal{N}^\circ \mathrm{d} \boldsymbol{k} . \end{align}

Summing (3.17) and (3.20) eliminates the energy-exchange term

(3.21) \begin{align} \int \boldsymbol{u}^{\scriptscriptstyle{L}}\,{\boldsymbol{\cdot}}\, \partial _t \boldsymbol{\mathsf{p}}\, \mathrm{d} z \end{align}

and produces energy conservation in (3.16) with

(3.22) \begin{align} \boldsymbol{J} = \int \boldsymbol{u}^{\scriptscriptstyle{L}} \pi \, \mathrm{d} z + \iint \boldsymbol{\mathcal{C}} \varOmega \mathcal{N} \, \mathrm{d} \boldsymbol{k} , \end{align}

and

(3.23) \begin{align} {E}^\circ = \iint \varOmega \mathcal{N}^\circ \, \mathrm{d} \boldsymbol{k} + \int \boldsymbol{u}^{\scriptscriptstyle{L}} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^\circ \, \mathrm{d} z . \end{align}

The right-hand side of (3.20),

(3.24) \begin{align} \iint \varOmega \mathcal{N}^\circ \mathrm{d} \boldsymbol{k} = \iint \sigma \mathcal{N}^\circ \mathrm{d} \boldsymbol{k} + \iint \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}} \, \mathcal{N}^\circ \mathrm{d} \boldsymbol{k} , \end{align}

involves a non-obvious contribution of the Doppler shift to ${E}^\circ$ . We return to this term in § 5.

Using the dispersion relation (2.2), the final term in (3.22) can be split into physically distinct contributions

(3.25) \begin{align} \iint \boldsymbol{\mathcal{C}} \varOmega \mathcal{N} \, \mathrm{d} \boldsymbol{k} = \iint \boldsymbol{\mathcal{C}} \mathcal{E}_{w} \, \mathrm{d} \boldsymbol{k} + \iint \mathcal{G} \boldsymbol{k} \boldsymbol{k}\,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}} \mathcal{N} \mathrm{d} \boldsymbol{k} + \iint \boldsymbol{\mathcal{U}} \boldsymbol{\mathcal{U}} \,{\boldsymbol{\cdot}}\, \boldsymbol{k} \mathcal{N} \, \mathrm{d} \boldsymbol{k} . \end{align}

The first term on the right-hand side of (3.25) is group-velocity transport of the wave energy $\mathcal{E}_{w}$ . The second term is the wave radiation stress (3.12) acting on $\boldsymbol{\mathcal{U}}$ . The third term is the momentum flux associated with $\boldsymbol{\mathcal{U}}$ acting on the wave pseudomomentum.

The mean kinetic energy in the conserved energy (3.15) is the LKE $|\boldsymbol{u}^{\scriptscriptstyle{L}}|^{2}/2$ rather than the Eulerian kinetic energy (EKE) $|\boldsymbol{u}^{\scriptscriptstyle{E}}|^{2}/2$ . For the CWCM, this is natural: $\boldsymbol{u}^{\scriptscriptstyle{L}}$ is the primary dynamical variable and $\boldsymbol{u}^{\scriptscriptstyle{E}}$ does not appear explicitly. Moreover, the form of the conserved energy is linked to the variational formulation of the CWCM in § 4 via Noether’s theorem.

Wagner et al. (Reference Wagner, Chini, Ramadhan, Gallet and Ferrari2021) and Czeschel & Eden (Reference Czeschel and Eden2023) contrast the simplicity of the LKE budget for the (uncoupled) Craik–Leibovich system with the complexity of the EKE budget. The CWCM explains this simplicity by the privileged role of the LKE as one of the two terms in the conserved total energy (3.15). In § 5 we return to the primacy of $|\boldsymbol{u}^{\scriptscriptstyle{L}}|^2/2$ over $|\boldsymbol{u}^{\scriptscriptstyle{E}}|^2/2$ in a discussion of wave-forced inertial oscillations. In this basic example the CWCM is forced by prescribing $\mathcal{N}^\circ$ in the action transport equation (2.1) and all terms in the energy balance can be calculated explicitly.

3.4. Circulation and potential vorticity

Material conservation of circulation and potential vorticity is a property of the Craik–Leibovich equation (2.8), i.e. these laws are independent of the wave-action equation (2.1). We state the results and defer to Andrews & McIntyre (Reference Andrews and McIntyre1978) or Holm (Reference Holm1996) for proof.

Introduce the ‘absolute Lagrangian mean momentum’

(3.26) \begin{align} \boldsymbol{\upsilon }^{{a}} = \boldsymbol{h} + \boldsymbol{u}^{\scriptscriptstyle{L}} - \boldsymbol{\mathsf{p}} , \end{align}

where

(3.27) \begin{align} \boldsymbol{h}(\boldsymbol{x}_3) = \tfrac {1}{2} \boldsymbol{f} \times \boldsymbol{x}_3 \end{align}

with $\boldsymbol{x}_3 = (\boldsymbol{x},z)$ . A key property of $\boldsymbol{h}$ is that ${{\boldsymbol{\nabla}} \times } \boldsymbol{h}= \boldsymbol{f}$ . The corresponding absolute vorticity is

(3.28) \begin{align} \boldsymbol{\zeta }^{{a}} = {{\boldsymbol{\nabla}} \times } \boldsymbol{\upsilon }^{{a}}=\boldsymbol{f} + {\boldsymbol{\nabla}} \times (\boldsymbol{u}^{\scriptscriptstyle{L}} - \boldsymbol{\mathsf{p}} ). \end{align}

Then conservation of circulation is

(3.29) \begin{align} \frac {\mathrm{d}}{\mathrm{d} t}\oint _{{\mathcal{C}}} \boldsymbol{\upsilon }^{{a}} \,{\boldsymbol{\cdot}}\, \mathrm{d} {\boldsymbol{\ell }}=0 , \end{align}

where $\mathcal{C}$ is a material circuit and we assume $\boldsymbol{u}^\circ = 0$ .

As noted by Andrews & McIntyre (Reference Andrews and McIntyre1978), for $\boldsymbol{f} = 0$ , conservation of circulation implies that ${\boldsymbol{\nabla}} \times (\boldsymbol{u}^{\scriptscriptstyle{L}} - \boldsymbol{\mathsf{p}}) = 0$ for all times if this holds initially. In this irrotational flow regime, $\boldsymbol{u}^{\scriptscriptstyle{L}}$ is diagnosed from $\boldsymbol{\mathsf{p}}$ as $\boldsymbol{u}^{\scriptscriptstyle{L}} = \boldsymbol{\mathsf{p}} - {\boldsymbol{\nabla}} \phi$ for a potential $\phi$ that is uniquely determined by the Poisson equation $\boldsymbol{\nabla} ^2 \phi = {\boldsymbol{\nabla}}\,{\boldsymbol{\cdot}}\, \boldsymbol{\mathsf{p}}$ and the boundary conditions satisfied by $\boldsymbol{u}^{\scriptscriptstyle{L}}$ . The CWCM thus reduces to the action equation (2.1), with $\boldsymbol{\mathcal{U}}(\boldsymbol{x},\kappa )$ related linearly to $\mathcal{N}(\boldsymbol{x},\boldsymbol{k})$ via solution of the Poisson equation.

For potential vorticity, suppose that we have a passive scalar $b(\boldsymbol{x},z,t)$ satisfying wave-averaged material conservation,

(3.30) \begin{align} \partial _t b + \boldsymbol{u}^{\scriptscriptstyle{L}}\,{\boldsymbol{\cdot}}\, {{\boldsymbol{\nabla}} } b = {b}^\circ . \end{align}

Potential vorticity is

(3.31) \begin{align} \varPi = \boldsymbol{\zeta }^{{a}} \,{\boldsymbol{\cdot}}\, {{\boldsymbol{\nabla}} } b , \end{align}

and satisfies

(3.32) \begin{align} \partial _t \varPi + \boldsymbol{u}^{\scriptscriptstyle{L}}\,{\boldsymbol{\cdot}}\, {{\boldsymbol{\nabla}} } \varPi ={\boldsymbol{\nabla}} \,{\boldsymbol{\cdot}}\, \left ( \boldsymbol{u}^\circ {\times } {{\boldsymbol{\nabla}} } b + \boldsymbol{\zeta }^{{a}} {b}^\circ \right ) \!. \end{align}

4.1. Rotating Euler equations

Variational formulations of fluid models typically involve the flow map, $\boldsymbol{\varphi }$ , which sends labels $\boldsymbol{a}$ identifying fluid particles to the time- $t$ positions of the particles, $\boldsymbol{x} = \boldsymbol{\varphi }(\boldsymbol{a},t)$ (e.g. Salmon Reference Salmon1988). Incompressibility implies that the flow map preserves volume, that is, $| \partial \boldsymbol{\varphi } / \partial \boldsymbol{a}| = \mathrm{const}$ . We can take the constant to be $1$ without loss of generality. Flow map and velocity field are related via

(4.1) \begin{align} \partial _t \boldsymbol{\varphi }(\boldsymbol{a},t)=\boldsymbol{u}(\boldsymbol{\varphi }(\boldsymbol{a},t),t). \end{align}

Following Dewar (Reference Dewar1970), Bretherton (Reference Bretherton1971) and Gjaja & Holm (Reference Gjaja and Holm1996) (see also Holm et al. (Reference Holm, Hu and Street2023)), we make minimal use of labels and flow map by adopting an Eulerian variational formulation. In this formulation – an instance of the Euler–Poincaré framework (Holm et al. Reference Holm, Marsden and Ratiu1998) – the flow map only appears implicitly through a constrained form of the variations of $\boldsymbol{u}$ first obtained by Newcomb (Reference Newcomb1962) and Bretherton (Reference Bretherton1970).

In Appendix B, we show that the rotating, incompressible Euler equations, with a free surface at $z = s(\boldsymbol{x},t)$ , follow from Hamilton’s principle $\delta \int \mathcal{L} \, \mathrm{d} t = 0$ applied to the Lagrangian

(4.2) \begin{align} \mathcal{L}[\boldsymbol{u},s] = \iint \mathrm{d} \boldsymbol{x} \left ( \int _{-d}^s \left ( \frac {1}{2}|\boldsymbol{u}|^2 + \boldsymbol{h} \,{\boldsymbol{\cdot}}\, \boldsymbol{u} \right ) \, \mathrm{d} z - \frac {1}{2} g s^2 \right )\!. \end{align}

The Eulerian variations $\delta \boldsymbol{u}$ and $\delta s$ are induced by variations $\delta \boldsymbol{\varphi }$ of the flow map. Introducing an Eulerian vector field $\boldsymbol{q}(\boldsymbol{x},t)=(q_1,q_2,q_3)$ by

(4.3) \begin{align} \delta \boldsymbol{\varphi }(\boldsymbol{a},t) = \boldsymbol{q}(\boldsymbol{\varphi }(\boldsymbol{a},t),t), \end{align}

leads to

(4.4a) \begin{align} \delta \boldsymbol{u} &= \partial _t \boldsymbol{q} + \boldsymbol{u} \,{\boldsymbol{\cdot}}\, \boldsymbol{\nabla }\boldsymbol{q} - \boldsymbol{q} \,{\boldsymbol{\cdot}}\, \boldsymbol{\nabla }\boldsymbol{u}, \\[-12pt] \nonumber \end{align}

(4.4b) \begin{align} \delta s &= q_3 - q_1 \partial _x s - q_2 \partial _y s, \\[10pt] \nonumber \end{align}

where $\boldsymbol{q}$ satisfies ${{\boldsymbol{\nabla}} } \,{\boldsymbol{\cdot}}\, \boldsymbol{q} = 0$ and is otherwise arbitrary (see Appendix B). Incompressibility is enforced by considering only volume-preserving flow maps and hence divergence-free $\boldsymbol{u}$ and $\boldsymbol{q}$ . No explicit constraint enforcing ${\boldsymbol{\nabla}} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}=0$ is therefore required in (4.2).

4.2. Lagrangian mean framework

Following the GLM framework, we introduce a wave–mean flow decomposition at the level of the flow map by writing

(4.5) \begin{align} \boldsymbol{\varphi }(\boldsymbol{a},t) = \bar {\boldsymbol{\varphi }}(\boldsymbol{a},t) + \boldsymbol{\xi }(\bar {\boldsymbol{\varphi }}(\boldsymbol{a},t),t), \end{align}

where $\bar {\boldsymbol{\varphi }}$ and $\boldsymbol{\xi }$ denote the mean flow and perturbation maps. Below, we introduce a linear-wave ansatz for $\boldsymbol{\xi }$ that captures the wavy displacements resulting from surface gravity waves; $\bar {\boldsymbol{\varphi }}(\boldsymbol{a},t)$ captures the mean currents.

Differentiating (4.5) with respect to $t$ and using the mean position ${\boldsymbol{x}_3}=(\boldsymbol{x},z)=\bar {\boldsymbol{\varphi }}(\boldsymbol{a},t)$ as independent variable gives

(4.6) \begin{align} \boldsymbol{u}^{\boldsymbol{\xi }} = \boldsymbol{u}^{\scriptscriptstyle{L}} + D_t \boldsymbol{\xi }. \end{align}

Here, $\boldsymbol{u}^{\scriptscriptstyle{L}}(\boldsymbol{x},t)$ defined by

(4.7) \begin{align} \boldsymbol{u}^{\scriptscriptstyle{L}}(\bar {\boldsymbol{\varphi }}(\boldsymbol{a}),t) = \partial _t \bar {\boldsymbol{\varphi }}(\boldsymbol{a},t) \end{align}

is the Lagrangian mean velocity,

(4.8) \begin{align} D_t = \partial _t + \boldsymbol{u}^{\scriptscriptstyle{L}} \,{\boldsymbol{\cdot}}\, {{\boldsymbol{\nabla}} } , \end{align}

and

(4.9) \begin{align} \boldsymbol{u}^{\boldsymbol{\xi }} (\boldsymbol{x}_3,t)= \boldsymbol{u}(\boldsymbol{x}_3 + \boldsymbol{\xi }(\boldsymbol{x}_3,t),t). \end{align}

The decomposition (4.5) involves a degree of arbitrariness. We take advantage of this to impose that the mean flow map preserves volume and leaves the free surface flat, that is,

(4.10) \begin{align} \left | \partial \bar {\boldsymbol{\varphi }} / \partial \boldsymbol{a} \right | = 1 \quad \textrm {and} \quad \bar \varphi _3(a,b,c=0)=0. \end{align}

In contrast with standard GLM, the assumptions in (4.10) imply a non-zero average for the perturbation map, $\overline {\boldsymbol{\xi }} \not =0$ (see Vanneste & Young (Reference Vanneste and Young2022)). Below in (4.12), however, we use a linear approximation for $\boldsymbol{\xi }$ so that $\overline {\boldsymbol{\xi }}$ vanishes up to neglected terms.

Volume preservation implies that $\boldsymbol{u}^{\scriptscriptstyle{L}}$ is divergence free, i.e. (2.9). With this, we can rewrite the Lagrangian (4.2) using the mean position $(\boldsymbol{x},z)= \bar {\boldsymbol{\varphi }}(\boldsymbol{a})$ as an independent variable. This results in

(4.11) \begin{align} \mathcal{L} [\boldsymbol{u}^{\scriptscriptstyle{L}},\boldsymbol{\xi } ] = \iint \mathrm{d} \boldsymbol{x} \left ( \int _{-d}^{0} \frac {1}{2} |\boldsymbol{u}^{\scriptscriptstyle{L}} + D_t \boldsymbol{\xi } |^2 + \boldsymbol{h}^{\boldsymbol{\xi }} \,{\boldsymbol{\cdot}}\, (\boldsymbol{u}^{\scriptscriptstyle{L}} + D_t \boldsymbol{\xi } ) \, \mathrm{d} z - \frac {1}{2} g \zeta ^2 \right ) \!, \end{align}

where $\zeta = \xi _3(\boldsymbol{x},z=0)$ and $\boldsymbol{h}^{\boldsymbol{\xi }}(\boldsymbol{x}_3)=\boldsymbol{h}(\boldsymbol{x}_3 + \boldsymbol{\xi }(\boldsymbol{x}_3))= \boldsymbol{h}(\boldsymbol{x}_3) + ({1}/{2}) \boldsymbol{f} {\times } \boldsymbol{\xi }(\boldsymbol{x}_3)$ . (Without the assumption of volume preservation (4.10), (4.11) would contain a Jacobian factor $\tilde \rho$ , as in the GLM formulation of Andrews & McIntyre (Reference Andrews and McIntyre1978).)

4.3. Averaged Lagrangian

Up to (4.11) no approximations have been made. We now make a Wentzel–Kramers–Brillouin (WKB) ansatz for $\boldsymbol{\xi }$ , assuming a scale separation between rapid phase variations and slow envelope modulations, and Whitham average. Following the prescription of Whitham (Reference Whitham1965, Reference Whitham1974), the displacement field $\boldsymbol{\xi }$ is approximated by that of a linear wavetrain,

(4.12) \begin{align} \boldsymbol{\xi } \approx \textit{Re}\, a(\boldsymbol{x},t) \hat {\boldsymbol{\xi }}(\boldsymbol{x},z,\boldsymbol{k}){\mathrm{e}}^{\mathrm{i} \theta (\boldsymbol{x},t)} , \end{align}

where $a(\boldsymbol{x},t)$ is a slowly varying amplitude, $\theta (\boldsymbol{x},t)$ the phase and $\boldsymbol{k}=(k,l) = {\boldsymbol{\nabla}} _{\boldsymbol{x}} \theta$ . The vertical structure is dictated by

(4.13) \begin{align} \hat {\boldsymbol{\xi }}(\boldsymbol{x},z,\boldsymbol{k}) = \frac {1}{\kappa \sinh (\kappa d(\boldsymbol{x}))} \left ( \begin{array}{c} \mathrm{i} k \cosh (\kappa (z+d(\boldsymbol{x})) \\[3pt] \mathrm{i} l \cosh (\kappa (z+d(\boldsymbol{x})) \\[3pt] \kappa \sinh (\kappa (z+d(\boldsymbol{x}))) \end{array} \right )\!, \end{align}

with $\kappa = |\boldsymbol{k}| = |{{\boldsymbol{\nabla}} }_{\boldsymbol{x}} \theta |$ . This structure is that obtained for surface waves propagating in a locally uniform current and with locally uniform depth as is usual with the WKB approximation. We also neglect the effect of rotation on the waves, assuming $f \ll \sigma$ .

With $\hat {\boldsymbol{\xi }}$ in (4.13), the displacement $\boldsymbol{\xi }(\boldsymbol{x},t)$ in (4.12) is approximately incompressible and irrotational:

(4.14) \begin{align} {{\boldsymbol{\nabla}} } \,{\boldsymbol{\cdot}}\, \boldsymbol{\xi } \approx 0 , \qquad \text{and} \qquad {{\boldsymbol{\nabla}} \times } \boldsymbol{\xi } \approx 0 . \end{align}

This approximation is associated with Whitham’s variational derivation of the WKB approximation: applying ${{\boldsymbol{\nabla}} }=({\boldsymbol{\nabla}} _{\!\boldsymbol{x}},\partial _z)$ to (4.12) the derivatives act only on the phase $\theta (\boldsymbol{x},t)$ and on the $z$ -dependence of $\hat {\boldsymbol{\xi }}$ . For instance,

(4.15) \begin{align} D_t \boldsymbol{\xi } \approx \textit{Re}\, a(\boldsymbol{x},t) \left ( \hat {\boldsymbol{\xi }}(\boldsymbol{x},z,\boldsymbol{k}) (-\mathrm{i} \omega + \mathrm{i} \boldsymbol{u}^{\scriptscriptstyle{L}}\,{\boldsymbol{\cdot}}\, \boldsymbol{k} ) + w^{\scriptscriptstyle{L}}(\boldsymbol{x},z)\hat {\boldsymbol{\xi }}_z(\boldsymbol{x},z,\boldsymbol{k}) \right ) \mathrm{e}^{\mathrm{i} \theta (\boldsymbol{x},t)} . \end{align}

We further neglect the term involving $w^{\scriptscriptstyle{L}} \hat {\boldsymbol{\xi }}_z$ in (4.15) on the grounds that $w^{\scriptscriptstyle{L}} \ll u^{\scriptscriptstyle{L}}, v^{\scriptscriptstyle{L}}$ .

Introducing (4.12) into the Lagrangian (4.11) and averaging over the phase $\theta$ leads to the averaged Lagrangian

(4.16) \begin{align} \overline {\mathcal{L}} [\boldsymbol{u}^{\scriptscriptstyle{L}},\theta ,a ] &= \iint \mathrm{d} \boldsymbol{x} \left ( \int _{-d}^0 \left ( \frac {1}{2} |\boldsymbol{u}^{\scriptscriptstyle{L}} |^2 + \boldsymbol{h}(\boldsymbol{x}) \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} \right ) \, \mathrm{d} z \right . \nonumber \\& \qquad \qquad + \left . \frac {1}{4} \left ( \int ^0_{-d} \frac {\cosh (2 \kappa (z+d)) (\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} )^2}{\sinh ^2(\kappa d)} \, \mathrm{d} z - g \right )|a|^2 \right ), \end{align}

where $\omega = -\partial _t \theta$ and we neglect the wave contributions to the Coriolis term.

Using the definition of $Q(\boldsymbol{x},z,\kappa )$ in (2.5), we rewrite the second $z$ -integral in (4.16) as

(4.17) \begin{align} \int ^0_{-d} \frac {\cosh (2 \kappa (z+d))(\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} )^2}{\sinh ^2(\kappa d)} \, \mathrm{d} z = \frac {1}{\kappa \tanh (\kappa d)} \int ^0_{-d} (\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} )^2 Q(\boldsymbol{x},z,\kappa ) \, \mathrm{d} z . \end{align}

Justified by the smallness of the Doppler shift against the wave frequency, we make the approximation

(4.18a) \begin{align} \int _{-d}^0 (\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} )^2 Q \, \mathrm{d} z &= \omega ^2 - 2 \omega \boldsymbol{k}\,{\boldsymbol{\cdot}}\, \int ^0_{-d}\!\! \boldsymbol{u}^{\scriptscriptstyle{L}} Q \, \mathrm{d} z + \int ^0_{-d}\! (\boldsymbol{k}\,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} )^2 Q \, \mathrm{d} z\, \\[-12pt] \nonumber \end{align}

(4.18b) \begin{align} &\approx (\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}})^2 ,\\[10pt] \nonumber \end{align}

where the Doppler velocity is

(4.19) \begin{align} \boldsymbol{\mathcal{U}} =\int ^0_{-d} \boldsymbol{u}^{\scriptscriptstyle{L}} Q \, \mathrm{d} z , \end{align}

as in (2.4). This manoeuvre simplifies (4.16) to

(4.20) \begin{align} \overline {\mathcal{L}} [\boldsymbol{u}^{\scriptscriptstyle{L}},\theta ,a ] = \iint \mathrm{d} \boldsymbol{x} \left ( \int _{-d}^0 \left ( \frac {1}{2} |\boldsymbol{u}^{\scriptscriptstyle{L}} |^2 + \boldsymbol{h}(\boldsymbol{x}) \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} \right ) \, \mathrm{d} z + \frac {1}{4} \left ( \frac {\left (\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}} \right )^2}{\kappa \tanh (\kappa d)} - g \right )|a|^2 \right ). \end{align}

The approximation in (4.18) is not necessary: one can proceed from the Lagrangian in (4.16). The resulting model has a more complicated dispersion relation than (2.2), and a more complicated expression for the pseudomomentum than (2.10). Alternatively, one obtains a superficially simpler model by boldly dropping the final term in (4.18a ). But this is one step too many: we prefer the dispersion relation (2.2) with a physically intuitive Doppler shift that preserves Galilean invariance.

The Lagrangian (4.20) is the surface-wave analogue of the Lagrangian obtained by Bretherton (Reference Bretherton1971), Gjaja & Holm (Reference Gjaja and Holm1996) and Salmon (Reference Salmon2016) for acoustic and gravity waves. Hamilton’s principle can be applied using variations $\delta \boldsymbol{u}^{\scriptscriptstyle{L}}$ similar to (4.4a ), $\delta a$ and $\delta \theta$ . This leads to the Craik–Leibovich system and to equations governing the phase $\theta (\boldsymbol{x},t)$ and amplitude $a(\boldsymbol{x},t)$ of a single wavetrain. We do not detail these equations here. Instead, we follow Dewar (Reference Dewar1970) to relax the assumption of a single wavetrain and consider a spectrum of linear waves.

4.4. Variational form of the CWCM

A spectrum of waves can be regarded as a superposition of independent wavetrains indexed by a two-dimensional (continuous) variable, $\boldsymbol{\alpha }$ say (see Dewar’s (Reference Dewar1970) discussion of the ‘weak turbulence’ assumption). Phase $\theta (\boldsymbol{x},\boldsymbol{\alpha },t)$ and amplitude $a(\boldsymbol{x},\boldsymbol{\alpha },t)$ then depend on $\boldsymbol{\alpha }$ , and the Lagrangian (4.20) acquires an additional integration over $\boldsymbol{\alpha }$ to become

(4.21) \begin{align} \overline {\mathcal{L}} [\boldsymbol{u}^{\scriptscriptstyle{L}},\theta ,a ] &= \iint \mathrm{d} \boldsymbol{x} \left ( \int _{-d}^0 \left ( \frac {1}{2} |\boldsymbol{u}^{\scriptscriptstyle{L}} |^2 + \boldsymbol{h}(\boldsymbol{x}) \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} \right ) \, \mathrm{d} z \right . \nonumber \\& \qquad \qquad + \left . \frac {1}{4} \iint \left (\frac {\left (\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}} \right )^2}{\kappa \tanh (\kappa d)}- g \right ) |a|^2 \, \mathrm{d} \boldsymbol{\alpha } \right ) \! . \end{align}

The representation using $(\boldsymbol{x},\boldsymbol{\alpha })$ as independent variables and $(\theta ,a)$ as dependent variables is inconvenient. It is preferable to use instead $(\boldsymbol{x},\boldsymbol{k})$ as independent variables and $(\omega ,\mathcal{N})$ as dependent variables, where

(4.22) \begin{align} \mathcal{N} = \frac {\sigma |a|^2}{2 \kappa \tanh (\kappa d)} \left |\frac {\partial \boldsymbol{k}}{\partial \boldsymbol{\alpha }} \right |^{-1} \!. \end{align}

This specific choice ensures that $\mathcal{N}$ is dual to the frequency $\omega$ in the sense that $\mathcal{N} = \delta \overline {\mathcal{L}} / \delta \omega$ (see Appendix C). This identifies $\mathcal{N}$ as the wave-action density in $(\boldsymbol{x},\boldsymbol{k})$ -space (Dewar Reference Dewar1970). Making the change of variables in (4.20) puts the averaged Lagrangian in its final form

(4.23) \begin{align} \mathcal{L}_{CWC}[\boldsymbol{u}^{\scriptscriptstyle{L}},\mathcal{N},\omega ] &= \iint \mathrm{d} \boldsymbol{x} \left ( \int _{-d}^0 \left ( \frac {1}{2} |\boldsymbol{u}^{\scriptscriptstyle{L}} |^2 + \boldsymbol{h}(\boldsymbol{x}) \,{\boldsymbol{\cdot}}\, \boldsymbol{u}^{\scriptscriptstyle{L}} \right ) \, \mathrm{d} z \right . \nonumber \\& \qquad \qquad \qquad \left . + \, \frac {1}{2} \iint ( \sigma ^{-1} (\omega - \boldsymbol{k} \,{\boldsymbol{\cdot}}\, \boldsymbol{\mathcal{U}})^2 - \sigma ) \mathcal{N} \, \mathrm{d} \boldsymbol{k} \right ) \!. \end{align}

The Lagrangian $\mathcal{L}_{CWC}$ is the Lagrangian of the CWCM. The term on the second line is a standard Whitham-averaged Lagrangian governing wave dynamics. Because of its dependence on the current variable $\boldsymbol{u}^{\scriptscriptstyle{L}}$ via $\boldsymbol{\mathcal{U}}$ , this term contributes to the dynamics of both waves and currents. This is key to consistency of the CWCM.

Hamilton’s principle can now be applied to (4.23). The variations are

(4.24) \begin{align} \delta \boldsymbol{u}^{\scriptscriptstyle{L}} = \partial _t \boldsymbol{q} + \boldsymbol{u}^{\scriptscriptstyle{L}} \,{\boldsymbol{\cdot}}\, {\boldsymbol{\nabla}} \boldsymbol{q} - \boldsymbol{q} \,{\boldsymbol{\cdot}}\, {\boldsymbol{\nabla}} \boldsymbol{u}^{\scriptscriptstyle{L}}, \end{align}

with $\boldsymbol{q}$ satisfying ${{\boldsymbol{\nabla}} } \,{\boldsymbol{\cdot}}\, \boldsymbol{q} =0$ but otherwise arbitrary (in analogy with (4.4a )), $\delta \mathcal{N}$ and

(4.25) \begin{align} \delta \omega = \partial _t r + {\boldsymbol{\nabla}} _{\boldsymbol{x}} r \,{\boldsymbol{\cdot}}\, {\boldsymbol{\nabla}} _{\boldsymbol{k}} \omega - {\boldsymbol{\nabla}} _{\boldsymbol{k}} r \,{\boldsymbol{\cdot}}\, {\boldsymbol{\nabla}} _{\boldsymbol{x}} \omega \end{align}

with $r$ arbitrary. This leads to the Craik–Leibovich equation in the form (2.16), dispersion relation (2.2) and action transport equation (2.1). Details are given in Appendix C. The form (4.25) of the variation of $\omega$ , reminiscent of (4.24), emerges when changing variables from the mixed Eulerian–Lagrangian $\theta (\boldsymbol{x},\boldsymbol{\alpha })$ to the Eulerian $\omega (\boldsymbol{x},\boldsymbol{k})$ (see Dewar (Reference Dewar1970) and Appendix D).