1. Introduction
Our understanding of the large-scale dynamics of the atmosphere and ocean rests on the concept of geostrophic balance: motion with time scales much longer than the inertial time scale 
$f^{-1}$, with 
$f$ the Coriolis parameter, is assumed to dominate, resulting in the familiar balance between Coriolis force, buoyancy and pressure gradients. The idea of dominant slow motion has stimulated decades of research extending balance beyond geostrophy to define dynamics which filter out fast motion to a high degree of accuracy – that is, dynamics restricted to a suitably defined slow manifold (e.g. Machenhauer Reference Machenhauer1977; Leith Reference Leith1980; Allen & Newberger Reference Allen and Newberger1993; Warn et al. Reference Warn, Bokhove, Shepherd and Vallis1995; McIntyre & Norton Reference McIntyre and Norton2000; Mohebalhojeh & Dritschel Reference Mohebalhojeh and Dritschel2001; Vanneste Reference Vanneste2013; Kafiabad & Bartello Reference Kafiabad and Bartello2017, Reference Kafiabad and Bartello2018). However, the amplitude of the fast motion, e.g. in the form of inertia-gravity waves, can often be large, corresponding to states well away from the slow manifold. For instance, the velocities associated with near-inertial waves (NIWs) in the ocean frequently exceed those of the (vortical) balanced motion (e.g. Alford et al. Reference Alford, MacKinnon, Simmons and Nash2016). There is, therefore, a need to reassess the notion of balance to account for the impact of waves. This can be achieved using reduced models which, rather than filtering the waves, average over their rapidly oscillating phases to describe the slow interactions between waves and vortical flow. We refer to these as wave-averaged models. Bretherton (Reference Bretherton1971), Grimshaw (Reference Grimshaw1975) and, more recently, Wagner & Young (Reference Wagner and Young2015) showed that wave averaging naturally leads to a generalised-Lagrangian-mean (GLM) description of the flow, in which averaging is carried out at fixed particle label rather than fixed Eulerian position. Wave-averaged models can therefore be interpreted as approximations to the GLM equations of Andrews & McIntyre (Reference Andrews and McIntyre1978) (see Bühler (Reference Bühler2014) for an introduction to GLM). The wave-averaged models of interactions between NIWs and quasi-geostrophic flow derived by Xie & Vanneste (Reference Xie and Vanneste2015) and Wagner & Young (Reference Wagner and Young2016) fall into this category. Salmon (Reference Salmon2016) provides an interesting variational perspective into this class of models in which weakly nonlinear internal waves interact perturbatively with quasi-geostrophic flow.
One of the most striking predictions of wave-averaged and GLM theories is that a form of hydrostatic and geostrophic balance continues to hold in the presence of strong waves, but that this balance applies to the Lagrangian mean flow, in the sense that
\begin{equation} f \boldsymbol{z} \times \bar{\boldsymbol{u}}^{L} = - \boldsymbol{\nabla} \bar{\rm \pi} + \bar{b}^{L} \boldsymbol{z}, \end{equation}
where 
$\boldsymbol {z}$ is the vertical unit vector, 
$\bar {\boldsymbol {u}}^{L}$ and 
$\bar {b}^{L}$ are the Lagrangian-mean velocity and buoyancy, and 
$\bar {{\rm \pi} }$ is a mean pressure-like scalar that includes a quadratic wave-averaged contribution (Moore Reference Moore1970; Wagner & Young Reference Wagner and Young2015; Xie & Vanneste Reference Xie and Vanneste2015; Gilbert & Vanneste Reference Gilbert and Vanneste2018; Thomas, Bühler & Smith Reference Thomas, Bühler and Smith2018). The balance in (1.1) follows expeditiously by simplification of the exact GLM momentum equation (Andrews & McIntyre (Reference Andrews and McIntyre1978), theorem I) in the small Rossby number limit. In § 3 we provide an alternative derivation of (1.1) proceeding from the Eulerian equations.
The Lagrangian-mean velocity, 
$\bar {\boldsymbol {u}}^{L}$ in (1.1), is the sum of Eulerian-mean and Stokes velocities,
\begin{equation} \bar{\boldsymbol{u}}^{L} = \bar{\boldsymbol{u}} + \bar{\boldsymbol{u}}^{S},\quad \textrm{where}\ \bar{\boldsymbol{u}}^{S} = \overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'}, \end{equation}
with 
$\boldsymbol {\xi }'$ and 
$\boldsymbol {u}'=\boldsymbol {\xi }'_t$ denoting wave displacement and velocity and the overline denoting (Eulerian) time averaging. The appearance of the Stokes velocity 
$\bar {\boldsymbol {u}}^{S}$ in the Coriolis term implies the existence of the Stokes–Coriolis force 
$f \boldsymbol {z} \times \bar {\boldsymbol {u}}^{S}$ whose importance for shallow currents driven by surface gravity waves has long been recognised (e.g. Ursell & Deacon Reference Ursell and Deacon1950; Hasselmann Reference Hasselmann1970; Huang Reference Huang1979; Leibovich Reference Leibovich1980; Bühler Reference Bühler2014). Here we are concerned with the Stokes–Coriolis force associated with internal gravity waves which operates in the ocean's interior. Equation (1.1) is a remarkable generalisation of the familiar geostrophic balance, with practical implications, e.g. for the interpretation of ocean velocities inferred from satellite altimetry data, yet (1.1) is not widely known. Another prediction of wave-averaged and GLM theories is the material conservation of a form of potential vorticity (PV) which combines the PV of the mean flow with a wave-averaged contribution. Together with (1.1) this leads to a wave-averaged form of quasi-geostrophic dynamics which captures the feedback exerted by waves on the mean flow.
The aim of this paper is to illustrate the validity and usefulness of the wave-averaged geostrophic balance and PV conservation by testing these relations against numerical solutions of the three-dimensional Boussinesq equations. Applying wave-averaged relations to numerical simulations poses difficulties related to the definition of a suitable time average and the estimation of particle displacements 
$\boldsymbol {\xi }'$ in (1.2). We sidestep these issues by focussing on a simple configuration: an NIW, initially with no horizontal dependence, and the vertical structure of a plane travelling wave, is superimposed on an axisymmetric (Gaussian) barotropic anticyclonic vortex. This configuration has the advantage that, because the NIW is approximately linear, it retains a simple vertical and temporal structure proportional to 
$\exp ({\mathrm {i}(mz-ft)})$, with 
$m$ the vertical wavenumber. Time averaging can therefore be replaced by a straightforward vertical averaging and the Eulerian-mean component of the flow can be identified with the barotropic component. Moreover, it turns out that for NIWs, the Stokes drift 
$\bar {\boldsymbol {u}}^{S}$ can be deduced from the wave kinetic energy, thus circumventing the need to estimate the displacements 
$\boldsymbol {\xi }'$.
We describe our numerical simulation of the interaction between NIWs and an anticyclonic vortex in § 2. The wave evolution follows a broadly understood scenario: after a rapid adjustment, wave energy concentrates in the vortex core, giving rise to an approximately axisymmetric trapped structure which is modulated periodically in time (Llewellyn Smith Reference Llewellyn Smith1999). Our focus is on the response of the mean flow to this wave pattern. In § 3, we briefly sketch a derivation of the wave-averaged balance (1.1) and of the wave-averaged PV. Particularising this to the case of vertically planar inertial wave, following Rocha, Wagner & Young (Reference Rocha, Wagner and Young2018), we relate (i) mean vertical vorticity to mean pressure and wave energy, and (ii) wave-induced change of mean vorticity to wave energy using PV conservation. We assess the accuracy of these predictions in numerical solutions and find remarkably good agreement in spite of the complexity of the full three-dimensional Boussinesq model and of the fact that the strong scaling assumptions required by the theory are only marginally satisfied for the parameters chosen. Our results illustrate the value of wave-averaged theories for the analysis of geophysical flows in the context of three-dimensional nonlinear dynamics.
2. Numerical solution of the Boussinesq equations
We analyse solutions of the non-hydrostatic Boussinesq equations for an initial condition consisting of an NIW superimposed on a barotropic anticyclonic vortex with the initial vertical vorticity
\begin{equation} \zeta_0(r) = - {\textit{Ro}}\, f \, \mathrm{e}^{-r^2/a^2}, \end{equation}
where 
$r$ is the radial coordinate, 
$a$ the vortex radius and 
${\textit {Ro}}$ the maximum Rossby number,
\begin{equation} {\textit{Ro}} = |\zeta_0|_{{max}} / f. \end{equation}
The NIW is initially horizontally uniform, with vertical structure 
$\mathrm {e}^{\mathrm {i} m z}$ with vertical wavenumber 
$m$ and initial energy 
$E_0$. In addition to the Rossby number, the flow is characterised by the Burger number
\begin{equation} {\textit{Bu}} = {N}^2/(\,{f m a} )^2, \end{equation}
where the inverse wavenumber 
$m^{-1}$ and vortex radius 
$a$ are used to define the characteristic vertical and horizontal scales, respectively. We focus on the regime where both 
${\textit {Ro}}$ and 
${\textit {Bu}}$ are small, as required for nearly geostrophic mean dynamics and NIW dynamics (Young & Ben Jelloul Reference Young and Ben Jelloul1997).
The Boussinesq equations are solved in a triply periodic domain, using a code adapted from that of Waite & Bartello (Reference Waite and Bartello2006), which relies on a dealiased pseudospectral method and a third-order Adams–Bashforth scheme with time step 
$0.06/f$. The 
$(2 {\rm \pi})^2 \times 2 {\rm \pi}/36$ domain is discretised using a 
$1152^2 \times 96$ uniform grid. A hyperdissipation of the form 
$\nu _{h} (\partial _x^2 + \partial _y^2)^4 +\nu _{z} \partial ^8_z$ is used for the momentum and density equations. The simulation parameters are listed in table 1. The domain size is much larger than the vortex radius (
${2{\rm \pi} }/a \approx 14$) to mitigate the effect of waves re-entering from the boundaries. The NIW is initialised with a horizontal velocity
\begin{equation} (u'_0,v'_0) = (\cos m z, \sin mz); \end{equation}
thus the initial NIW kinetic energy density is 
$E_0 = ( {u'_0}^2 + {v'_0}^2)/2=1/2$. The vortex strength is such that the Eulerian mean-flow kinetic energy density is 
$0.1$.
Table 1. Simulation parameters. The horizontal domain size, 
$L = 2{\rm \pi}$, determines the unit of length; the unit of time is defined so that the initial energy density of the NIW is 
$E_0=1/2$. The vortex strength is such that the initial Eulerian domain-averaged mean-flow kinetic energy density is 
$0.1$; the Gaussian vortex in (2.1) has maximum azimuthal velocity 
$0.32{\textit {Ro}}\, f a = 1.45$ at 
$r = 1.13 a$.
The evolution of the NIW in the presence of a vortex is illustrated by the sequence of snapshots of horizontal slices shown in figure 1. The NIW kinetic energy 
$({u'}^2+{v'}^2)/2$ is displayed in the top row, the vertical vorticity in the middle row and the (barotropic) mean vorticity in the bottom row (this is also the Eulerian mean vorticity). The wave-energy concentration is modulated periodically in time, with a period much larger than the inertial period. The dynamics of the vertical vorticity 
$\zeta$ is dominated by the fast oscillation at the inertial frequency induced by NIW advection of 
${\zeta}_{0}$. This oscillation is filtered out by vertical averaging, leaving an Eulerian mean vorticity 
$\bar {\zeta }$ that is almost axisymmetric and oscillates slowly in unison with the wave energy.
Figure 1. (a) Horizontal slices of wave kinetic energy 
$({u'}^2 + {v'}^2)/2$, (b) vertical vorticity 
$\zeta$ and (c) change in Eulerian-mean vertical vorticity 
$\bar {\zeta }-\zeta _0$ at times (from left to right) 
$t=0$, 20.4, 43.2, 76.3, 99.2 and 122.1 inertial periods 
$2{\rm \pi} /f$.
In the top row of figure 1 the NIW energy level in the vortex core increases by over a factor of five above that of the initial condition (2.4). This focussing by the mean flow results in an approximately axisymmetric trapped structure; this is a finite-amplitude, near-inertial normal mode of the vortex (Llewellyn Smith Reference Llewellyn Smith1999). In anticipation of the strong wave expansion in (3.2) below, we emphasise that this concentration of NIW energy happens even if the mean-flow velocity is much less than the NIW velocity: weak mean flows shift the lowest frequency of the internal wave band to 
$f + \bar {\zeta }/2$ (Kunze Reference Kunze1985; Young & Ben Jelloul Reference Young and Ben Jelloul1997). Thus vortices with 
$\bar {\zeta }/2<0$ support trapped modes with frequencies slightly less than 
$f$, and the initial condition (2.4) projects onto this mode. We now turn to the mean-flow effects of this trapped mode.
3. Wave-averaged analysis
We first derive the wave-averaged geostrophic balance (1.1) and wave-averaged PV conservation for general inertia-gravity waves, then consider their application to NIWs and the numerical solution of § 2.
3.1. The strong-wave expansion
We provide only an informal sketch of the derivation and refer the reader to Wagner & Young (Reference Wagner and Young2015) for details (and Thomas et al. (Reference Thomas, Bühler and Smith2018) for related shallow-water results). We start with the Boussinesq equations
\begin{gather} \boldsymbol{u}_t + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + f \boldsymbol{z} \times \boldsymbol{u} + \boldsymbol{\nabla} p= b \boldsymbol{z}, \end{gather}
\begin{gather}b_t + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} b + N^2 w =0, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} =0. \end{gather}
For simplicity we take the Brunt–Väisälä frequency 
$N$ constant (see Wagner & Young (Reference Wagner and Young2015) for the non-constant case).
We assume that the flow consists of fast waves, with small amplitude 
$\varepsilon \ll 1$, interacting with a slow flow, with smaller amplitude 
$\varepsilon ^2$. This is the strong-wave assumption which makes it possible to capture the impact of the waves on the flow without the need to carry out the perturbation expansion to an unwieldy high order. The flow and wave amplitudes are assumed to vary on the slow time scale of quasi-geostrophic dynamics, which is 
$O(\varepsilon ^{-2} f)$ in our notation. We expand all dynamical fields as
\begin{equation} \boldsymbol{u} = \varepsilon \boldsymbol{u}' + \varepsilon^2 \bar{\boldsymbol{u}} + \cdots . \end{equation}
Here the prime identifies the wave component, the overbar the Eulerian-mean component obtained by averaging over the fast wave time scale, and the dots include a rapidly varying 
$O(\varepsilon ^2)$ term as well as 
$O(\varepsilon ^3)$ terms. Introducing (3.2) into (3.1) gives, at order 
$\varepsilon$, the linear inertia-gravity-wave equations
\begin{gather} \boldsymbol{u}'_t + f \boldsymbol{z} \times \boldsymbol{u}' + \boldsymbol{\nabla} p' = b' \boldsymbol{z}, \end{gather}
\begin{gather}b'_t + N^2 w' =0 , \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}' =0. \end{gather}Averaging the next-order equations gives
\begin{gather} \overline{\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'} + f \boldsymbol{z} \times\bar{ \boldsymbol{u} }+ \boldsymbol{\nabla} \bar{p} = \bar{b} \boldsymbol{z} , \end{gather}
\begin{gather}\overline{\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} b'} + N^2 \bar{w} =0, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \bar{\boldsymbol{u}} =0 . \end{gather}
We can rewrite the nonlinear terms in (3.4a) in terms of the wave displacement 
$\boldsymbol {\xi }'$. For instance, we have
\begin{align} \overline{\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'} &= - \overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'_t} \end{align}
\begin{align} &= f\boldsymbol{z} \times \overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'} + \overline{(\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{\nabla} p'} - \overline{(\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla})b} \boldsymbol{z} \end{align}
\begin{align} &= f \boldsymbol{z} \times \bar{\boldsymbol{u}}^{S} + \boldsymbol{\nabla} \tfrac{1}{2} \bar{p}^{S} - \bar{b}^{S} \boldsymbol{z}. \end{align}In passing from (3.5b) to (3.5c) we have used the remarkable result
\begin{equation} (\overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{\nabla} p'} = \tfrac{1}{2} \boldsymbol{\nabla} (\overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla} p'}) = \tfrac{1}{2} \boldsymbol{\nabla} \bar{p}^{S}, \end{equation}
where 
$\bar {p}^{S} = \overline {\boldsymbol {\xi }' \boldsymbol {\cdot } \boldsymbol {\nabla } p'}$ is the Stokes pressure. The identity (3.6) is established in Wagner & Young (Reference Wagner and Young2015) and a simplified proof is given in appendix A. Using (3.5c) the mean momentum equations in (3.4a) are reduced to the wave-averaged balance equations (1.1) with
\begin{equation} \bar{b}^{L} = \bar{b} + \bar{b}^{S} \quad \textrm{and} \quad \bar{\rm \pi} = \bar{p} + \tfrac{1}{2} \bar{p}^{S}, \end{equation}
where 
$\bar {b}^{S}= \overline {\boldsymbol {\xi }' \boldsymbol {\cdot } \boldsymbol {\nabla } b'}$ is the Stokes buoyancy. Finally, with 
$\overline {\boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla } b'}= - \overline {\boldsymbol {\xi }'\boldsymbol {\cdot } \boldsymbol {\nabla } b'_t} = N^2 \overline {\boldsymbol {\xi '} \boldsymbol {\cdot } \boldsymbol {\nabla } w'}$, one finds from (3.4b) that 
$\bar {w}^{L} = \bar {w} + \bar {w}^{S} = 0$.
It is natural to introduce the stream function 
$\bar {\psi }^{L} = \bar {{\rm \pi} }/f$ to write the components of the wave-averaged balance relations (1.1) as
\begin{equation} \bar{\boldsymbol{u}}^{L} = (-\bar{\psi}^{L}_y, \bar{\psi}^{L}_x,0) \quad \textrm{and} \quad \bar{b}^{L} = \bar{\psi}^{L}_z, \end{equation}mirroring the familiar geostrophic and hydrostatic relations of quasi-geostrophy. The vorticity of the Lagrangian-mean flow is then
\begin{equation} \bar{\zeta}^{L} = \bar{v}^{L}_x - \bar{u}^{L}_y = {\rm \Delta} \bar{\psi}^{L}, \end{equation}
where 
${\rm \Delta} = \partial _x^2 + \partial _y^2$ is the horizontal Laplacian. We emphasise that 
$\bar {\psi }^{L}$ and 
$\bar {\zeta }^{L}$ are the stream function and vorticity corresponding the Lagrangian-mean velocity 
$(\bar {u}^{L},\bar {v}^{L})$, but are not themselves Lagrangian means of any specific fields. This is evident from the factor 
$1/2$ in (3.7a,b) and the fact that 
$\bar {\zeta }^{L} \not = \bar {\zeta } + \overline {\boldsymbol {\xi }' \boldsymbol {\cdot } \boldsymbol {\nabla } \zeta '}$ (the difference is related to the curl of the pseudomomentum, see e.g. Bühler (Reference Bühler2014) or Gilbert & Vanneste (Reference Gilbert and Vanneste2018)).
3.2. Wave-averaged balance with NIWs
The horizontal velocity 
$(u',v')$ of the vertically planar NIWs can be written in the complex form
\begin{equation} u'(x,y,z,t) + \mathrm{i} \, v'(x,y,z,t) = \phi(x,y,t) \exp({\mathrm{i} (m z - f t)}), \end{equation}
where 
$\phi$ is a slowly varying complex amplitude, known as the back-rotated velocity. The form in (3.10) is an excellent approximation for the waves in our simulation. It makes clear that vertical averaging is equivalent to time averaging over the fast wave time scale. The wave kinetic energy is 
$({u'}^2 + {v'}^2)/2=|\phi |^2/2$ and thus the top row of figure 1 shows the evolution of 
$|\phi |^2/2$. The back-rotated velocity 
$\phi$ evolves on the slow time scale 
$\varepsilon ^{-2} t$, where 
$\varepsilon$ is the order parameter in the strong-wave expansion (3.2). This long-term evolution is obtained by proceeding beyond the leading-order wave equation in (3.3) so that mean-flow effects, such as advection and 
$\zeta$-refraction, are revealed in the Young and Ben Jelloul (known as YBJ) equation (Young & Ben Jelloul Reference Young and Ben Jelloul1997; Asselin & Young Reference Asselin and Young2019, Reference Asselin and Young2020).
Using the back-rotated velocity 
$\phi$ dramatically simplifies the wave-averaged balance (1.1). First, 
$p' \approx 0$ for NIWs, so 
$\bar {{\rm \pi} } = \bar {p}$ and taking the horizontal divergence of (1.1) gives
\begin{equation} f \bar{\zeta}^{L} = {\rm \Delta} \bar{p}. \end{equation}Second, for vertically planar NIWs the horizontal part of the Stokes velocity can be computed as
\begin{equation} (\bar{u}^{S},\bar{v}^{S})=(\mathcal{A}_y,-\mathcal{A}_x), \quad \textrm{where}\ \mathcal{A} = |\phi|^2/(2f) \end{equation}is the wave action density, whose negative plays the role of a stream function for the Stokes velocity (see Rocha et al. (Reference Rocha, Wagner and Young2018) and appendix A). Combining (3.11) and (3.12) relates the Eulerian-mean vertical vorticity to the pressure according to
\begin{equation} f\left( \bar{\zeta} - {\rm \Delta} \mathcal{A} \right)= {\rm \Delta} \bar{p}. \end{equation}
This is the form of wave-averaged geostrophy that we test in our numerical simulation. Figure 2 shows radial profiles (obtained by azimuthal averaging) of the Eulerian-mean vorticity 
$\bar {\zeta }$, vorticity of the Lagrangian-mean flow 
$\bar {\zeta }^{L}$ and horizontal Laplacian of mean pressure 
${\rm \Delta} \bar {p}$ evaluated from the simulation. Figure 2 shows snapshots at two times corresponding to a maximum in the wave-energy concentration (panel a) and to a nearly uniform wave field (panel b). It is clear that, when the wave energy is non-uniform (panel a), geostrophic balance holds to a much better accuracy in the wave-averaged sense (3.11) than in the conventional sense: it is the Lagrangian-mean velocity, rather than the Eulerian-mean velocity, that is in geostrophic balance.
Figure 2. Radial profiles of 
$f$ times the Lagrangian- and Eulerian-mean vertical vorticity and of 
${\rm \Delta} \bar {p}$. The times (a) 
$t = 43.2$ and (b) 
$t=76.3$ inertial periods correspond to maximally concentrated and nearly uniform wave energy, respectively.
We have verified that the small difference between 
$\bar {\zeta }^{L}$ and 
${\rm \Delta} \bar { p}$ in figure 2 results from the cyclostrophic correction to vortex geostrophy (not shown). One restriction on the utility of wave-averaged balance is that the Stokes–Coriolis force, proportional to 
$f \bar {\boldsymbol {u}}^{S}$, should be larger than non-wave ageostrophic terms such as the aforementioned cyclostrophic correction. The formal justification of this restriction is in (3.2): the Stokes–Coriolis force is of order 
$\varepsilon ^2$ while the non-wave ageostrophic terms are of order 
$\varepsilon ^4$. For the numerical solution summarised in table 1 there is not a large difference between the mean and NIW kinetic energies, e.g. the maximum vortex velocity is only a factor of approximately two smaller than the maximum NIW velocity attained at 
$t = 43.2$. Nonetheless, in this numerical solution the wave-averaged correction to geostrophic balance is much more important than cyclostrophic effects and other higher-order non-wave corrections to balance.
3.3. Potential vorticity
We can go beyond the diagnostic relation (3.8a,b) and make dynamical predictions using the conservation of PV. After removing the constant contribution 
$N^2 f$ and scaling by 
$N^2$, the PV can be written as
\begin{equation} q = v_x - u_y + f b_z/N^2 + \left( \boldsymbol{\nabla} \times \boldsymbol{u} \right) \boldsymbol{\cdot} \boldsymbol{\nabla} b/N^2. \end{equation}
Because the linearised inertia-gravity wave equations (3.3) imply that 
$v'_x - u'_y + f b'_z/N^2=0$, the PV is an 
$O(\varepsilon ^2)$ quantity. Its material conservation 
$\partial _t q + \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } q =0$ shows that its leading-order approximation is evolving on the slow, 
$O(\varepsilon ^{-2} f^{-1})$, time scale and hence that it can be approximated by its time average. Using this we obtain, to leading order,
\begin{equation} q = \bar{\zeta}^{L} + {f^2} \bar{\psi}^{L}_{zz} /{N^2} \underbrace{ - \bar{v}^{S}_x + \bar{u}^{S}_y - f \bar{b}^{S}_z/N^2 + \overline{\left( \boldsymbol{\nabla} \times \boldsymbol{u}' \right) \boldsymbol{\cdot} \boldsymbol{\nabla} b'}/N^2}_{q^{\mathrm{w}}}, \end{equation}
where (3.8a,b) has been used. The first two terms on the right-hand side of (3.15) make up the familiar quasi-geostrophic PV, here involving the stream function of the Lagrangian-mean flow. The remaining terms constitute the wave PV denoted 
$q^{{w}}$. We note that the wave-averaged geostrophic balance and PV conservation under advection by 
$\bar {\boldsymbol {u}}^{L}$ also follow directly from GLM theory (Andrews & McIntyre Reference Andrews and McIntyre1978; Bühler & McIntyre Reference Bühler and McIntyre1998; Holmes-Cerfon, Bühler & Ferrari Reference Holmes-Cerfon, Bühler and Ferrari2011; Xie & Vanneste Reference Xie and Vanneste2015; Gilbert & Vanneste Reference Gilbert and Vanneste2018).
For NIWs, 
$q^{{w}}$ can be expressed in terms of the back-rotated velocity 
$\phi$; retaining only leading-order terms in the small Burger number the result is
\begin{equation} q = \bar{\zeta}^{L} + {f^2} \bar{\psi}^{L}_{zz} /{N^2} + {\rm \Delta} |\phi|^2/(4f) + \mathrm{i} \partial(\phi,\phi^*)/(2f), \end{equation}
where 
$\partial (\cdot ,\cdot )$ denotes the Jacobian operator (see Xie & Vanneste Reference Xie and Vanneste2015; Wagner & Young Reference Wagner and Young2016). For an axisymmetric wave field 
$\phi (r,t)$ and barotropic mean flow (3.16) reduces to
\begin{equation} q = \bar{\zeta}^{L} + {\rm \Delta} \mathcal{A}/2, \end{equation}
where 
${\rm \Delta} (\cdot ) = r^{-1} \partial _r(r \partial _r \cdot )$. For waves that are initially uniform in the horizontal, as in the initial condition (2.4), the Stokes drift and horizontal Laplacian vanish initially, so the conservation of 
$q$ – which holds pointwise since 
$q$ is axisymmetric and 
$\bar {\boldsymbol {u}}^{L}$ has no radial component – implies that
\begin{equation} \bar{\zeta}^{L} = \zeta_0 - {\rm \Delta} \mathcal{A}/2, \end{equation}
where 
$\zeta _0(r)$ is the initial Gaussian vertical vorticity in (2.1). Using the expression for the Stokes velocity in (3.12), the Eulerian-mean vorticity is
\begin{equation} \bar{\zeta} = \zeta_0 + {\rm \Delta} \mathcal{A}/2. \end{equation}
This result explains the correlation between the evolution of the wave kinetic energy 
$({u'}^2 + {v'}^2)/2 = f \mathcal {A}$ and 
$\bar {\zeta }$ observed in figure 1. We test the prediction in (3.19) quantitatively in figure 3 by comparing the change in mean vorticity, 
$\bar {\zeta }-\zeta _0$, with 
${\rm \Delta} \mathcal {A}/2$ at three different times corresponding to different phases in the oscillation of wave energy. The good match between the two quantities confirms the validity of (3.19). Thus conservation of wave-averaged PV predicts the dynamics of the mean flow in the presence of waves with substantial amplitude.
Figure 3. Comparison of mean-vorticity change 
$\bar {\zeta }-\zeta _0$ (
$\circ$, red) with 
${\rm \Delta} \mathcal {A}/2$ (solid line) along a radial line and for (a) 
$t=12.7$, (b) 
$t=30$ and (c) 
$t=43.2$ inertial periods.
4. Discussion
This paper examines the way geostrophic balance – a cornerstone of geophysical fluid dynamics – is altered in the presence of inertia-gravity waves to become a balance between the Coriolis force associated with the Lagrangian-mean velocity and a wave-modified mean pressure. The wave-induced correction to geostrophy has long been known (e.g. Moore Reference Moore1970) but it has received much less attention than the issue of non-wave finite-Rossby-number corrections associated with higher-order balance (e.g. Machenhauer Reference Machenhauer1977; Leith Reference Leith1980; Allen & Newberger Reference Allen and Newberger1993; Warn et al. Reference Warn, Bokhove, Shepherd and Vallis1995; McIntyre & Norton Reference McIntyre and Norton2000; Mohebalhojeh & Dritschel Reference Mohebalhojeh and Dritschel2001; Vanneste Reference Vanneste2013; Kafiabad & Bartello Reference Kafiabad and Bartello2017, Reference Kafiabad and Bartello2018). By focusing on numerical simulations in a very simple set-up, in which vertically planar NIW waves are superimposed on a barotropic vortex, we are able to show unambiguously that the wave-induced correction matters and, for NIWs, can be estimated from the wave energy, consistent with theoretical predictions, at least in the simple case considered here. We further show that the material conservation of a wave-averaged PV can be exploited to predict mean-flow changes from the wave energy. We emphasise that the concepts of wave-averaged balance and wave-averaged PV are not limited to the present case of vertically planar NIWs propagating through a barotropic mean flow. They apply in the presence of vertically sheared mean flows, as in the numerical NIW simulations in Asselin & Young (Reference Asselin and Young2020), Thomas & Arun (Reference Thomas and Arun2020) and Thomas & Daniel (Reference Thomas and Daniel2020), and to generic inertia-gravity waves. In the case of baroclinic mean flows, however, the waves cannot be extracted by removing the vertical average as done here. The standard normal-mode decomposition (e.g. Bartello Reference Bartello1995) also falls short for this task: in the case of strong waves, because of the importance of the Stokes–Coriolis force 
$f \boldsymbol {z} \times \bar {\boldsymbol {u}}^{S}$, balanced flow is poorly approximated by the linearised part of PV, and an appropriate form of time filtering is required to separate waves from balanced flow.
To assess the importance of wave-induced corrections to geostrophic balance in an oceanic context, we need to compare the size of the Stokes velocity induced by inertia-gravity waves with typical mean-flow velocities. We focus here on the case of NIWs in the ocean's mixed layer represented, as is standard, by a slab model (e.g. Alford et al. Reference Alford, MacKinnon, Simmons and Nash2016). We note that the relation 
$\bar {\boldsymbol {u}}^{S} = (-\mathcal {A}_y,\mathcal {A}_x)$ applies to (vertically independent) waves in a slab model as well as to the vertically planar waves considered so far. We estimate 
$\mathcal {A}$ from the values of inertial-wave kinetic energy inferred by Chaigneau, Pizarro & Rojas (Reference Chaigneau, Pizarro and Rojas2008) from drifter data. In parts of the ocean with strong wind forcing, they found this kinetic energy per unit area to be of the order of 
$10^3$ J m
$^{-2}$. This corresponds to 
$({u'}^2+{v'}^2)/2 = 10^3/(\rho h) \approx 10^{-2}$ m
$^2$ s
$^{-2}$, taking the water density 
$\rho =10^3$ kg m
$^{-3}$ and mixed-layer depth 
$h=100$ m, and hence 
$\mathcal {A} \approx 10^2$ m
$^2$ s
$^{-1}$ for 
$f = 10^{-4}$ s
$^{-1}$. Using a spatial scale of 10 km based on recent simulations of NIWs (Asselin & Young Reference Asselin and Young2020), we finally estimate 
$|\bar {\boldsymbol {u}}^{S}| \approx 10^{-2}$ m s
$^{-1}$, a substantial fraction of typical geostrophic velocities at the ocean surface. Larger values are obtained for smaller 
$f$ (lower latitudes) and a shallower mixed layer. The NIW Stokes velocity is also comparable in magnitude to the Stokes drift of surface gravity waves. However, unlike the latter, near-inertial Stokes drift extends across the mixed layer and into the ocean interior, hence it is potentially more important for transport. The contribution from other types of inertia-gravity waves, including internal tides, also deserves attention.
Our order-of-magnitude estimate raises the possibility that wave-averaged corrections to geostrophic balance affect the inference of sea-surface geostrophic velocity from satellite altimetric measurements. Assuming that sea-surface-height measurements provide an approximation to the Eulerian-mean height and hence, for NIWs, also of the Lagrangian-mean height (since NIWs have little effect on the sea-surface height), the surface velocity inferred by standard geostrophy is a Lagrangian mean rather than an Eulerian mean, as normally assumed. It would be useful to examine whether and how this difference can be accounted for.
We conclude by noting that, while the present paper concentrates on the relation between wave energy and wave-induced mean flow, it is also desirable to analyse the mechanism that leads to the slow wave-energy oscillations seen in figure 1. We leave this analysis for another publication (Kafiabad, Vanneste & Young Reference Kafiabad, Vanneste and Young2020).
Acknowledgements
H.A.K. and J.V. are supported by the UK Natural Environment Research Council grant NE/R006652/1. W.R.Y. is supported by the National Science Foundation Award OCE-1657041. This work used the ARCHER UK National Supercomputing Service.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation details
We establish (3.6) by rewriting the wave momentum equation (3.3a) as
\begin{equation} \boldsymbol{\nabla} p' = \mathsf{L} \boldsymbol{\xi}',\quad \textrm{where}\ \mathsf{L} = \begin{pmatrix} - \partial_t^2 & f \partial_t & 0 \\ -f \partial_t & - \partial_t^2 & 0 \\ 0 & 0 & -\partial_t^2 - N^2 \end{pmatrix}. \end{equation}
The linear operator 
$\mathsf {L}$ is self-adjoint in the sense that 
$\overline {\boldsymbol {a} \boldsymbol {\cdot } \mathsf {L} \boldsymbol {b}} = \overline {\boldsymbol {b} \boldsymbol {\cdot } \mathsf {L} \boldsymbol {a}}$ for arbitrary time-periodic vectors 
$\boldsymbol {a}$ and 
$\boldsymbol {b}$. Using this we compute
\begin{equation} \partial_i ( \overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla} p'}) = \partial_i ( \overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \mathsf{L} \boldsymbol{\xi}'}) = 2 \overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \mathsf{L} \partial_i \boldsymbol{\xi}'} = 2 \overline{\boldsymbol{\xi}' \boldsymbol{\cdot} \partial_i \boldsymbol{\nabla} p'} = 2 \overline{(\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla}) \partial_i p'} \end{equation}and (3.6) follows.
We now turn to (3.12). The unaveraged Stokes velocity is
\begin{equation} (\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla}) {{\boldsymbol {u}}} '= {\partial}_t\big[ \tfrac{1}{2} (\boldsymbol{\xi}' \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{\xi}' \big] - \tfrac{1}{2} \boldsymbol{\nabla} \times {\boldsymbol {h}}, \end{equation}
where 
${\boldsymbol {h}} =\boldsymbol {\xi }'\times {{\boldsymbol {u}}}'$ contributes to the absolute angular momentum density, and we used the standard vector identity for the curl of a cross product and that 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\xi }' = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}' = 0$. The time average of (A3) gives the explicit solenoidal representation of the Stokes velocity
\begin{equation} \bar{\boldsymbol{u}}^{S} = - \tfrac{1}{2} \boldsymbol{\nabla} \times \bar{{\boldsymbol {h}}}. \end{equation}
For vertically planar waves 
$\partial _z \bar {{\boldsymbol {h}}} = 0$, hence
\begin{equation} (\bar{u}^{S},\bar{v}^{S})= \tfrac{1}{2} (-\bar{h}_y,\bar{h}_x), \end{equation}
where 
$\bar {h}= \overline {\xi 'v' - \eta 'u'}$ is the vertical component of 
$\bar {{\boldsymbol {h}}}$. For NIWs, 
$u'_t - f v' =v'_t + f u' = 0$, so 
$\bar {h}=(\overline {\xi ' u'_t + \eta ' v'_t})/f = - (\overline {{u'}^2 +{v'}^2})/f = - \mathcal {A}$ and (A5) gives (3.12).
Open access
