Scattering of internal tides by barotropic quasigeostrophic flows

Abstract

Oceanic internal tides and other inertia–gravity waves propagate in an energetic turbulent flow whose length scales are similar to the wavelengths. Advection and refraction by this flow cause the scattering of the waves, redistributing their energy in wavevector space. As a result, initially plane waves radiated from a source such as a topographic ridge become spatially incoherent away from the source. To examine this process, we derive a kinetic equation which describes the statistics of the scattering under the assumptions that the flow is quasigeostrophic, barotropic and well represented by a stationary homogeneous random field. Energy transfers are quantified by computing a scattering cross-section and shown to be restricted to waves with the same frequency and identical vertical structure, hence the same horizontal wavelength. For isotropic flows, scattering leads to an isotropic wave field. We estimate the characteristic time and length scales of this isotropisation, and study their dependence on parameters including the energy spectrum of the flow. Simulations of internal tides generated by a planar wavemaker carried out for the linearised shallow-water model confirm the pertinence of these scales. A comparison with the numerical solution of the kinetic equation demonstrates the validity of the latter and illustrates how the interplay between wave scattering and transport shapes the wave statistics.

1 Introduction

The propagation of inertia–gravity waves in the ocean has received a great deal of attention, mainly motivated by the role they play in the large- and mesoscale circulation, through wave–mean-flow interaction, mixing and dissipation. The inertia–gravity wave spectrum is dominated by two types of waves: near-inertial oscillations, with frequencies close to the inertial frequency $f$ , which are mainly generated by winds, and internal tides (ITs), primarily at the semi-diurnal lunar frequency, which are generated by the interaction of the barotropic tide with topography (e.g. Ferrari & Wunsch 2009). Near-inertial oscillations have distinctive dynamics, including weak dispersion and weak vertical motion, that stems from their unique place at the low-frequency end of the inertia–gravity wave spectrum (Alford et al. 2016); ITs, in contrast, are generic mid-frequency inertia–gravity waves with their externally imposed frequency as their sole defining property.

The ocean’s highly energetic quasigeostrophic turbulence has a strong impact on the structure of both inertial oscillations and ITs and hence, in the case of ITs, on their signature on the sea-surface height (Rainville & Pinkel 2006; Ray & Zaron 2016). There is by now an extensive literature devoted to this impact, with a recent impetus provided by upcoming high-resolution satellite-altimetry instruments and the need to disentangle ITs from mesoscale (balanced) motion in the observed sea-surface height. We refer the reader to the recent papers by Dunphy et al. (2017) and Wagner, Ferrando & Young (2017) for further background.

A key aspect of the interactions between quasigeostrophic turbulence and both near-inertial oscillations and low-mode ITs is that turbulence and waves share similar horizontal scales, of the order of 100 km. A consequence is that, in such a regime, the Wentzel–Kramers–Brillouin (WKB) approximation on which much of our understanding of inertia–gravity wave propagation is built is not valid. This has prompted the development of simplified, wave-averaged models that rely only on time scale separation to represent the interactions between waves and flow in a simplified manner. Models of this kind include the Young–Ben Jelloul model of near-inertial oscillations in a quasigeostrophic flow (Young & Ben Jelloul 1997) and its extensions accounting for the feedback of the waves on the flow (Xie & Vanneste 2015; Wagner & Young 2016; Thomas, Smith & Bühler 2017). Wagner et al. (2017) recently derived an analogue of the Young–Ben Jelloul model equation for ITs. This model is formulated in physical space and retains a stiff term which enforces the constraint that the ITs’ fixed frequency imposes on their spatial structure and which cannot be eliminated without resorting to a Fourier-space formulation.

The present paper focuses similarly on the impact of a quasigeostrophic flow on ITs, making no asymptotic assumption about the relative horizontal scales of ITs and flow. Our focus is on quantifying the scattering induced by a barotropic (i.e. $z$ -independent) geostrophic turbulent flow which we model as a spatially homogeneous random field. We take advantage of the assumption of barotropic flow to use a vertical-mode expansion and thus reduce the problem to the study of an (equivalent) shallow-water model. By applying the theory of waves in random media developed by Ryzhik, Papanicolaou & Keller (1996), we derive a kinetic equation describing, in a statistically averaged sense, the energy exchanges between ITs with different wavevectors. The theory is formulated in terms of a wavevector-resolving energy density, $a(\boldsymbol{x},\boldsymbol{k},t)$ , which makes it possible to capture spatial variations of the wave energy. The form of the scattering term in the kinetic equation for $a(\boldsymbol{x},\boldsymbol{k},t)$ shows that energy transfers are restricted to waves with the same frequency or, equivalently, the same wavenumber $|\boldsymbol{k}|$ . These transfers result from interactions within resonant triads consisting of two ITs of equal frequencies with a zero-frequency flow (vortical) mode – the so-called catalytic interactions of Lelong & Riley (1991) and Bartello (1995). The rate of these transfers is proportional to the energy spectrum of the geostrophic flow. In the case of an isotropic flow, the scattering leads to the relaxation of the energy density towards a locally isotropic density $a(\boldsymbol{x},|\boldsymbol{k}|)$ . Our theory complements that developed by Ward & Dewar (2010), shifting from a deterministic to a statistical treatment that can be regarded as a version of wave turbulence (e.g. Nazarenko 2011) in which the statistics of the flow are prescribed. It generalises the theory developed by Danioux & Vanneste (2016) for near-inertial oscillations to inertia–gravity waves of arbitrary frequencies. Note that our statistical approach focuses on a homogeneous field of scatterers (resulting from the turbulent flow) and that different analysis techniques apply to waves incident on isolated scatterers (e.g. Olbers 1981).

We analyse the predictions of the kinetic equation, focusing our attention on parameters representative of the first baroclinic mode of the semidiurnal lunar tide $M_{2}$ . These predictions include a time scale for wave isotropisation applicable to statistically homogeneous wave fields (i.e. such that $\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a=0$ ) and in particular to the isotropisation of an initially plane wave examined numerically by Ward & Dewar (2010) in a shallow-water set-up. The kinetic equation applies to more general, non-homogeneous situations in which the energy density is modulated spatially ( $\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a\neq 0$ ). This makes it possible to study the scattering of ITs generated by a localised source such as a topographic ridge. Ponte & Klein (2015) and Dunphy et al. (2017) recently used three-dimensional Boussinesq simulations to study this problem and quantify the temporal incoherence of the ITs that results from the presence of a time-dependent turbulent flow. (See also Kelly et al. (2016) and Kelly & Lermusiaux (2016) for simulations of ITs in realistic configurations.) We carry out shallow-water simulations in a set-up analogous to theirs, and compare the results with direct solutions of the kinetic equation. This provides an estimate for the length scale over which the wave field becomes isotropic and, more broadly, sheds light on the interplay between transport of the wave energy by the group velocity and scattering. We emphasise that our theory concentrates on the statistical properties of the IT energy and makes no predictions for their phase. In the regime we consider, with a flow assumed to vary on a time scale much larger than the tidal period, the stationarity of the turbulent-energy spectrum ensures that the tidal energy remains concentrated at the single wavenumber dictated by the fixed tide frequency.

The plan of the paper is as follows. We describe the equations satisfied by linear internal waves propagating on a barotropic quasigeostrophic flow in § 2, expanding them in vertical modes to obtain an equivalent shallow-water system for each mode. We then sketch the derivation of the kinetic equation using the method of Ryzhik et al. (1996), relegating the technical computations to appendix B. We focus on the application of the kinetic equation to the case of an isotropic flow in § 3 where we derive explicit estimates for the time- and length scales over which the IT field becomes isotropic. In § 4 we compare our theoretical predictions with direct simulations of the linearised shallow-water equations and with numerical solutions of the kinetic equation itself in a configuration where a wavemaker forces a plane IT in a turbulent flow. We conclude in § 5 with a discussion.

2 Scattering theory for internal tides

2.1 Model

We model the propagation of ITs through a turbulent quasigeostrophic eddy field using the hydrostatic Boussinesq equations linearised about a slowly evolving barotropic flow. The background flow is time dependent and geostrophically and hydrostatically balanced, given by $\boldsymbol{U}=(U,V,0)=(-\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D713},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713},0)$ in terms of a $z$ -independent streamfunction $\unicode[STIX]{x1D713}$ . With these assumptions, the linearised hydrostatic Boussinesq equations read

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}+f\hat{\boldsymbol{z}}\times \boldsymbol{u}=-\unicode[STIX]{x1D735}p,\\ \unicode[STIX]{x2202}_{z}\,p=b,\\ \unicode[STIX]{x2202}_{t}b+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}b+N^{2}w=0,\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}+\unicode[STIX]{x2202}_{z}w=0,\end{array}\right\}\end{eqnarray}$$

where $(\boldsymbol{u},w)$ denotes the IT velocity, $\hat{\boldsymbol{z}}$ is the vertical unit vector, $p$ is the pressure normalised by a reference density, $b$ the buoyancy, $f$ the Coriolis parameter and $N(z)$ the buoyancy frequency. We use the notation $\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y},0)$ for the horizontal gradient throughout.

Assuming a flat bottom boundary and rigid lid, we project (2.1) onto baroclinic modes to obtain a set of rotating shallow-water equations governing their amplitudes:

(2.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{m}+\boldsymbol{u}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{m}+f\hat{\boldsymbol{z}}\times \boldsymbol{u}_{m}=-g\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}_{m},\\ \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D702}_{m}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}_{m}+h_{m}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{m}=0,\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D702}_{m}$ is the equivalent surface height and $h_{m}$ is the equivalent depth (see appendix A for details). Note that this system differs from the one obtained by linearising the shallow-water equations about a background flow in geostrophic balance since the latter system includes a contribution from the (sloping) background free surface.

For physical applications in later sections, we take parameters corresponding to the first baroclinic mode only, since this contains the majority of the IT energy. We drop the subscript $m$ in (2.2) from this point on. In the ocean, energy is transferred between vertical modes as a result of vertical shear. However, as discussed by Dunphy & Lamb (2014) and Ponte & Klein (2015) the effect is small. Simulations in Dunphy et al. (2017) put the transfer of energy from the first mode to higher modes at 3 % in their most extreme cases, with a highly energetic background flow, and less for typical ocean conditions.

2.2 Derivation of the kinetic equation

We study IT scattering in the distinguished regime where the spatial scale of the flow, $L_{\ast }$ say, is of the same order as the wavelength, that is, $|\boldsymbol{k}|L_{\ast }=O(1)$ , where $\boldsymbol{k}=(k,l)$ is the IT horizontal wavevector. The assumption of a geostrophic flow requires a small Rossby number $Ro=U_{\ast }/(fL_{\ast })\ll 1$ , where $U_{\ast }$ is a typical flow velocity; in turn, this implies that the background flow velocities are small compared with the wave phase speed $\unicode[STIX]{x1D714}/|\boldsymbol{k}|$ , where

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D714}=\sqrt{f^{2}+gh|\boldsymbol{k}|^{2}}\end{eqnarray}$$

is the IT frequency, since $U_{\ast }/(\unicode[STIX]{x1D714}/|\boldsymbol{k}|)=O(U_{\ast }/(\unicode[STIX]{x1D714}L_{\ast }))=O(Ro)$ , given that $\unicode[STIX]{x1D714}=O(f)$ away from the equator. With the flow time scale $T_{\ast }$ taken as the natural advective time scale $L_{\ast }/U_{\ast }$ , this also implies that the flow evolves slowly compared with the IT time scale since $\unicode[STIX]{x1D714}T_{\ast }=O(Ro)\ll 1$ . We further assume that, while the IT phases vary over the length scale $|\boldsymbol{k}|^{-1}$ , their amplitudes vary over a much larger scale $\unicode[STIX]{x1D700}^{-1}|\boldsymbol{k}|$ , where $\unicode[STIX]{x1D700}\ll 1$ . We adopt the scaling $\unicode[STIX]{x1D700}=O(Ro^{2})$ . As emerges below, this is the distinguished scaling that ensures that transport and scattering affect the wave field at the same order and are both captured at leading order by our asymptotic model.

Since our focus is on generic, statistical properties of the IT field, we model the turbulent background flow by a random streamfunction with homogeneous and stationary statistics. With our scaling assumptions, it is then possible to derive a single equation that describes the scattering and transport of IT energy following the theory of Ryzhik et al. (1996). This theory is formulated in terms of the Wigner transform

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D652}(\boldsymbol{x},\boldsymbol{k},t)=\frac{1}{(2\unicode[STIX]{x03C0})^{2}}\int _{\mathbb{R}^{2}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{y}}\unicode[STIX]{x1D753}(\boldsymbol{x}-\boldsymbol{y}/2,t)\unicode[STIX]{x1D753}^{\ast }(\boldsymbol{x}+\boldsymbol{y}/2,t)\,\text{d}\boldsymbol{y},\end{eqnarray}$$

where $\unicode[STIX]{x1D753}=(u,v,\unicode[STIX]{x1D702})^{\text{T}}$ groups the dynamical fields and the asterisk denotes conjugate transpose. The Wigner transform is a Hermitian matrix; it is not necessarily positive definite, although its integral over the wavevector $\boldsymbol{k}$ , simply given by $\unicode[STIX]{x1D753}(\boldsymbol{x},t)\unicode[STIX]{x1D753}^{\ast }(\boldsymbol{x},t)$ , is.

The equation derived by Ryzhik et al. (1996) governs the evolution of a scalar amplitude $a(\boldsymbol{x},\boldsymbol{k},t)$ which appears naturally in an eigenvector decomposition of the matrix $\unicode[STIX]{x1D652}(\boldsymbol{x},\boldsymbol{k},t)$ (see appendix B for details). Physically, this amplitude is interpreted as a wavevector-resolving energy density, related to the (leading-order) energy density of the system by

(2.5) $$\begin{eqnarray}{\mathcal{E}}_{0}(\boldsymbol{x},t)=\frac{1}{2}\int _{\mathbb{R}^{2}}a(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}.\end{eqnarray}$$

To avoid any confusion, we emphasise that $a(\boldsymbol{x},\boldsymbol{k},t)$ itself represents a wave-energy density and not a wave amplitude; as its definition in (B 29) makes clear, it is a quadratic function of the wave fields, like the Wigner transform $\unicode[STIX]{x1D652}(\boldsymbol{x},\boldsymbol{k},t)$ .

In appendix B, we show that $a(\boldsymbol{x},\boldsymbol{k},t)$ satisfies the kinetic equation

(2.6) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}a+\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a={\mathcal{L}}a-\unicode[STIX]{x1D6F4}a,\end{eqnarray}$$

where the notation $\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}=\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y})$ emphasises that the spatial gradient applies to functions of both $\boldsymbol{x}$ and $\boldsymbol{k}$ , and where $\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}=(\unicode[STIX]{x2202}_{k},\unicode[STIX]{x2202}_{l})$ . Here $\unicode[STIX]{x1D714}$ is determined by the IT dispersion relation (2.3), so that $\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}$ is the group velocity and the left-hand side of (2.6) represents the familiar wave transport. (The term $-\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}a$ would be added if $\unicode[STIX]{x1D714}$ depended explicitly on $\boldsymbol{x}$ .) The right-hand side represents wave scattering by the background flow. The first term, given by

(2.7) $$\begin{eqnarray}{\mathcal{L}}a(\boldsymbol{x},\boldsymbol{k},t)=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })a(\boldsymbol{x},\boldsymbol{k}^{\prime },t)\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$

quantifies the transfers of energy from all wavevectors $\boldsymbol{k}^{\prime }$ into wavevector $\boldsymbol{k}$ that result from interactions with the background flow; the second term, where

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D6F4}=\unicode[STIX]{x1D6F4}(\boldsymbol{k})=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$

is the total scattering cross-section, quantifies the energy lost by wavevector $\boldsymbol{k}$ to all other wavevectors.

The function $\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })$ that appears in (2.7)–(2.8) is the main object of interest. It is known as the differential scattering cross-section and measures the rate at which energy is scattered from $\boldsymbol{k}^{\prime }$ to $\boldsymbol{k}$ at a position $\boldsymbol{x}$ in space. We obtain it in the form

(2.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime }) & = & \displaystyle \frac{2\unicode[STIX]{x03C0}}{gh\unicode[STIX]{x1D714}^{3}|\boldsymbol{k}|^{5}}\{\!|\boldsymbol{k}\times \boldsymbol{k}^{\prime }|^{2}[(\unicode[STIX]{x1D714}^{2}+f^{2})\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime }-f^{2}|\boldsymbol{k}|^{2}]^{2}\nonumber\\ \displaystyle & & \displaystyle +\,f^{2}\unicode[STIX]{x1D714}^{2}[|\boldsymbol{k}\times \boldsymbol{k}^{\prime }|^{2}+\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime }(|\boldsymbol{k}|^{2}-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime })]^{2}\!\}\frac{\widehat{E}(\boldsymbol{k}-\boldsymbol{k}^{\prime })}{|\boldsymbol{k}-\boldsymbol{k}^{\prime }|^{2}}\unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-|\boldsymbol{k}^{\prime }|),\end{eqnarray}$$

where $\boldsymbol{k}\times \boldsymbol{k}^{\prime }=kl^{\prime }-lk^{\prime }$ , and $\widehat{E}$ is the energy spectrum of the flow. We note that $\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })$ is real, positive, and symmetric with respect to the exchange between $\boldsymbol{k}$ and $\boldsymbol{k}^{\prime }$ . In appendix B, we also show that these properties ensure conservation of the leading-order energy density (2.5):

(2.10) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}{\mathcal{E}}_{0}+\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}\boldsymbol{\cdot }{\mathcal{F}}_{0}=0,\end{eqnarray}$$

where

(2.11) $$\begin{eqnarray}{\mathcal{F}}_{0}(\boldsymbol{x},t)=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}(\boldsymbol{k})\,a(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}\end{eqnarray}$$

is the leading-order energy flux (see (B 56)).

The presence of the factor $\unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-|\boldsymbol{k}^{\prime }|)$ in (2.9) implies that energy is only exchanged between wavevectors of the same magnitude, that is, between waves with the same frequency, as a result of the assumed slow time dependence of the background flow. Thus, in the regime considered, the IT energy is confined to the constant-frequency circle $|\boldsymbol{k}|=((\unicode[STIX]{x1D714}^{2}-f^{2})/(gh))^{1/2}$ in the wavevector plane. This can be related to the observation that the background flow only enters $\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })$ through its energy spectrum $\widehat{E}$ , which, for the statistically stationary flows considered, is time independent. The scattering described by (2.6) results from the resonant interactions of two ITs, with wavevectors $\boldsymbol{k}$ and $\boldsymbol{k}^{\prime }$ and identical frequencies, with a vortical flow mode of wavevector $\boldsymbol{k}-\boldsymbol{k}^{\prime }$ and zero frequency. Because of potential-vorticity conservation, these interactions would leave the vortical mode unaffected even if it were allowed to evolve freely; hence they have been termed catalytic interactions (Lelong & Riley 1991; Bartello 1995; Ward & Dewar 2010). We emphasise that (2.6) captures the net effect of multiple triadic interactions acting over long time scales. This is why the time scale of evolution is not linear in the flow amplitude but quadratic, dictated by the energy spectrum of the flow, in a manner familiar from wave turbulence (e.g. Nazarenko 2011).

3 Scattering in isotropic turbulence

3.1 Isotropisation

In this section we use the kinetic equation (2.6) to make predictions about the scattering process and quantify the time and length scales over which ITs lose their spatial coherence. For simplicity, we assume that the flow is isotropic, $\widehat{E}(\boldsymbol{k})=\widehat{E}(|\boldsymbol{k}|)$ . This motivates the use of polar coordinates for the wavevector, such that

(3.1a,b ) $$\begin{eqnarray}\boldsymbol{k}=|\boldsymbol{k}|\left(\begin{array}{@{}c@{}}\cos \unicode[STIX]{x1D703}\\ \sin \unicode[STIX]{x1D703}\end{array}\right)\quad \text{and}\quad \boldsymbol{k}^{\prime }=|\boldsymbol{k}^{\prime }|\left(\begin{array}{@{}c@{}}\cos (\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}^{\prime })\\ \sin (\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}^{\prime })\end{array}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D703}^{\prime }$ is the angle between $\boldsymbol{k}$ and $\boldsymbol{k}^{\prime }$ . The change of coordinates reduces the scattering operator (2.7) to

(3.2) $$\begin{eqnarray}{\mathcal{L}}a(\boldsymbol{x},|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)=\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703}^{\prime })a(\boldsymbol{x},|\boldsymbol{k}|,\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime },t)\,\text{d}\unicode[STIX]{x1D703}^{\prime },\end{eqnarray}$$

where

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703}^{\prime }):=\int _{0}^{\infty }\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })|\boldsymbol{k}^{\prime }|\,\text{d}|\boldsymbol{k}^{\prime }|.\end{eqnarray}$$

Note that we have used the evenness of $\unicode[STIX]{x1D70E}^{\prime }$ in $\unicode[STIX]{x1D703}^{\prime }$ to rewrite (3.2) as a convolution. Note also that $\unicode[STIX]{x1D70E}^{\prime }$ is independent of the direction $\unicode[STIX]{x1D703}$ of $\boldsymbol{k}$ because the scattering process is rotationally invariant. The scattering cross-section (2.9) can be written explicitly as

(3.4) $$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})=\frac{\unicode[STIX]{x03C0}|\boldsymbol{k}|^{2}}{gh\unicode[STIX]{x1D714}^{3}}\{(\!1+\cos \unicode[STIX]{x1D703})[(\unicode[STIX]{x1D714}^{2}+f^{2})\cos \unicode[STIX]{x1D703}-f^{2}]^{2}+\unicode[STIX]{x1D714}^{2}f^{2}[\!1-\cos (3\unicode[STIX]{x1D703})]\}\widehat{E}(2|\boldsymbol{k}\sin (\unicode[STIX]{x1D703}/2)|),\end{eqnarray}$$

where we have removed the prime from $\unicode[STIX]{x1D703}^{\prime }$ for convenience. Similarly, the total scattering cross-section (2.8) reduces to

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D6F4}(|\boldsymbol{k}|)=\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

In the limit $\unicode[STIX]{x1D714}\rightarrow f$ corresponding to near-inertial waves the cross-section reduces to

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})=\frac{2\unicode[STIX]{x03C0}f|\boldsymbol{k}|^{2}}{gh}\widehat{E}(2|\boldsymbol{k}\sin (\unicode[STIX]{x1D703}/2)|),\end{eqnarray}$$

which recovers the result obtained by Danioux & Vanneste (2016) starting from the Young–Ben Jelloul model.

With the scattering cross-section (3.4), the scattering operator (3.2) can be diagonalised using a Fourier series, or more precisely a cosine series since $a$ is even in $\unicode[STIX]{x1D703}$ . Denoting the cosine transform by a hat, with

(3.7) $$\begin{eqnarray}\widehat{a}_{n}(\boldsymbol{x},|\boldsymbol{k}|,t)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\cos (n\unicode[STIX]{x1D703})a(\boldsymbol{x},|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

we find that

(3.8) $$\begin{eqnarray}(\widehat{{\mathcal{L}}a})_{n}=\unicode[STIX]{x1D706}_{n}\widehat{a}_{n},\quad n=0,1,\ldots ,\end{eqnarray}$$

with the eigenvalues

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{n}=\unicode[STIX]{x1D706}_{n}(|\boldsymbol{k}|):=2\unicode[STIX]{x03C0}\widehat{\unicode[STIX]{x1D70E}^{\prime }}=\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})\cos (n\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

Fourier transforming the kinetic equation (2.6) then gives

(3.10)

It follows from (3.9) and the non-negativity of $\unicode[STIX]{x1D70E}^{\prime }$ in (3.4) that

(3.11a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{0}=\unicode[STIX]{x1D6F4}(|\boldsymbol{k}|)\quad \text{and}\quad |\unicode[STIX]{x1D706}_{n\geqslant 1}|<\unicode[STIX]{x1D706}_{0}.\end{eqnarray}$$

Thus the scattering term on the right-hand side of (3.10) vanishes for $n=0$ and represents a damping for $n\geqslant 1$ .

The implications are clearly seen for a wave field that is spatially homogeneous, that is, with $\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a=0$ : the solution of (2.6), with initial condition

(3.12) $$\begin{eqnarray}a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t=0)=A(|\boldsymbol{k}|,\unicode[STIX]{x1D703}),\end{eqnarray}$$

is then simply

(3.13) $$\begin{eqnarray}a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)=\mathop{\sum }_{n=0}^{\infty }\widehat{A}_{n}(|\boldsymbol{k}|)\text{e}^{(\unicode[STIX]{x1D706}_{n}-\unicode[STIX]{x1D6F4})t}\cos (n\unicode[STIX]{x1D703}).\end{eqnarray}$$

This describes the relaxation of the solution towards a stationary, isotropic ( $\unicode[STIX]{x1D703}$ -independent) solution, since

(3.14) $$\begin{eqnarray}\lim _{t\rightarrow \infty }a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)=\widehat{A}_{0}(|\boldsymbol{k}|)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t=0)\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

This is a key feature of the scattering: the main impact of the random isotropic flow is to lead to the isotropisation of the IT field regardless of the initial condition. Note that, with $a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)$ the wave-energy density, $\widehat{A}_{0}$ represents the total, $\unicode[STIX]{x1D703}$ -integrated energy at wavevector $|\boldsymbol{k}|$ , while the $\widehat{A}_{n}$ for $n\neq 0$ capture the energy’s dependence on  $\unicode[STIX]{x1D703}$ .

We can identify two time scales for the scattering process. First, the scattering time

(3.15) $$\begin{eqnarray}T_{\mathit{scat}}=\unicode[STIX]{x1D6F4}^{-1}\end{eqnarray}$$

estimates the time over which energy concentrated at $\boldsymbol{k}$ in spectral space is reduced by a factor of $\text{e}^{-1}$ while converted to waves with other wavevectors. In other words, it is the time scale over which scattering effects become significant. Second, the time scale for convergence to an isotropic wave field is given by

(3.16) $$\begin{eqnarray}T_{\mathit{iso}}=(\unicode[STIX]{x1D6F4}-\unicode[STIX]{x1D706}^{\prime })^{-1},\quad \text{where}\;\unicode[STIX]{x1D706}^{\prime }:=\max _{n\geqslant 1}\unicode[STIX]{x1D706}_{n}.\end{eqnarray}$$

This is the time for the last surviving anisotropic (i.e. $n\neq 0$ ) Fourier mode to decay by a factor of $\text{e}^{-1}$ . Scattering length scales associated with the time scales (3.15) and (3.16) can be defined as

(3.17a,b ) $$\begin{eqnarray}L_{\mathit{scat}}=c_{g}T_{\mathit{scat}},\quad L_{\mathit{iso}}=c_{g}T_{\mathit{iso}},\end{eqnarray}$$

where $c_{g}=|\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}|=gh|\boldsymbol{k}|/\unicode[STIX]{x1D714}(|\boldsymbol{k}|)$ is the group speed.

3.2 Predicted behaviour

In this section, we use the time and length scales (3.15)–(3.17) to examine how the scattering depends on the Coriolis parameter $f$ , and on the strength and horizontal scales of the eddies as encoded in $\widehat{E}$ . Since we focus on ITs, we regard the frequency $\unicode[STIX]{x1D714}$ as fixed and deduce $|\boldsymbol{k}|$ from the dispersion relation (2.3). We test some of our predictions against numerical simulations in § 4.

We assume an isotropic energy spectrum of the form

(3.18) $$\begin{eqnarray}\widehat{{\mathcal{E}}}(|\boldsymbol{k}|):=2\unicode[STIX]{x03C0}|\boldsymbol{k}|\,\widehat{E}(|\boldsymbol{k}|)=\left\{\begin{array}{@{}ll@{}}c_{1}|\boldsymbol{k}|,\quad & |\boldsymbol{k}|\leqslant \unicode[STIX]{x1D705},\\ c_{2}|\boldsymbol{k}|^{-3.5},\quad & |\boldsymbol{k}|\geqslant \unicode[STIX]{x1D705}.\end{array}\right.\end{eqnarray}$$

This depends on two parameters: $\unicode[STIX]{x1D705}$ , a peak wavenumber which sets the dominant length scale of the flow, and the root-mean-square velocity defined by $v_{\mathit{rms}}^{2}=\int _{0}^{\infty }\widehat{{\mathcal{E}}}\,\text{d}|\boldsymbol{k}|$ . The constants $c_{1}$ and $c_{2}$ are determined by $\unicode[STIX]{x1D705}$ and $v_{\mathit{rms}}$ and the requirement of continuity at $|\boldsymbol{k}|=\unicode[STIX]{x1D705}$ . In practice, we choose $\unicode[STIX]{x1D705}$ so that the correlation length $l_{c}:=\unicode[STIX]{x03C0}/k_{c}$ , where $k_{c}:=\iint |\boldsymbol{k}|\widehat{E}\,\text{d}\boldsymbol{k}/\iint \widehat{E}\,\text{d}\boldsymbol{k}$ , is similar to the wavelength of the IT; calculation of the integrals using (3.18) gives $\unicode[STIX]{x1D705}=9\unicode[STIX]{x03C0}/(10l_{c})$ . Although quasigeostrophic theory predicts a kinetic-energy spectrum that decays as $|\boldsymbol{k}|^{-3}$ for balanced geostrophic turbulence, the slope is often observed to be slightly steeper, a result which is typically attributed to the presence of large-scale coherent structures that emerge in the turbulent flow (McWilliams, Weiss & Yavneh 1994; Bartello 1995; Kafiabad & Bartello 2016). This motivates the form of (3.18) as representative of balanced geostrophic eddy fields in the ocean. Note that we have chosen an energy spectrum with non-zero energy for all $|\boldsymbol{k}|$ . This matters because only the range $[0,2|\boldsymbol{k}|]$ of the energy spectrum contributes to the scattering of ITs with wavenumber $|\boldsymbol{k}|$ , as the factor $\widehat{E}(2|\boldsymbol{k}\sin (\unicode[STIX]{x1D703}^{\prime }/2)|)$ in the scattering cross-section (3.4) indicates. A lower cutoff of the spectrum, say at some wavenumber $k_{cut}$ , would then imply that waves with $|\boldsymbol{k}|\leqslant k_{cut}/2$ are unaffected by scattering. This is related to the resonant triad view that two ITs with the same wavenumber $|\boldsymbol{k}|$ and hence the same frequency can only form a resonant triad with a flow mode if this has a wavenumber in $[0,2|\boldsymbol{k}|]$ .

Figure 1. Scattering cross-section $\unicode[STIX]{x1D70E}^{\prime }$ in (3.4) for the energy spectrum (3.18) as a function of the peak wavenumber $\unicode[STIX]{x1D705}$ and angle $\unicode[STIX]{x1D703}$ treated as polar coordinates. The IT wavenumber is fixed as $|\boldsymbol{k}|=3\times 10^{-5}~\text{m}^{-1}$ corresponding to the mode-1 $M_{2}$ tide at $45^{\circ }$ latitude, and the flow’s root-mean-square velocity as $v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$ . (a) Peak wavenumber in the range $0.2|\boldsymbol{k}|\leqslant \unicode[STIX]{x1D705}\leqslant |\boldsymbol{k}|$ , (b) in the range $0.5|\boldsymbol{k}|\leqslant \unicode[STIX]{x1D705}\leqslant 2.5|\boldsymbol{k}|$ . Three (circular) contours of $\unicode[STIX]{x1D705}$ are labelled in each panel in units of $10^{-5}~\text{m}^{-1}$ .

Figure 1 shows the scattering cross-section $\unicode[STIX]{x1D70E}^{\prime }$ in (3.4) as a function of $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D705}$ for $v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$ and for $f=1.028\times 10^{-4}~\text{s}^{-1}$ corresponding to $45^{\circ }$ latitude. The equivalent depth is set to $h=1.2~\text{m}$ , as appropriate for the first baroclinic mode (Olbers, Willebrand & Eden 2012), and the frequency to $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}/12.42~\text{h}^{-1}$ , corresponding to the $M_{2}$ tide. The horizontal wavenumber is then $|\boldsymbol{k}|=3\times 10^{-5}~\text{m}^{-1}$ corresponding to a wavelength of approximately $200~\text{km}$ . The figure indicates that scattering is local in the angular coordinate, that is, ITs are preferentially scattered into waves with nearby directions. This is especially the case for small values of $\unicode[STIX]{x1D705}$ , corresponding to flows with typical scales much larger than the IT wavelength (figure 1 a), when the values of $\unicode[STIX]{x1D70E}^{\prime }$ are also the largest. For larger values of $\unicode[STIX]{x1D705}$ , that is, for flows with smaller scales, the energy transfers are slower, but less localised in the angular direction (figure 1 b).

The net effect of the scattering depends on both the value of $\unicode[STIX]{x1D70E}^{\prime }$ at fixed $\unicode[STIX]{x1D703}$ and the range of $\unicode[STIX]{x1D703}$ where $\unicode[STIX]{x1D70E}^{\prime }$ is substantial; it is best measured by the scattering and isotropisation time and length scales introduced in (3.15)–(3.17). These scales are deduced from the cosine transform of $\unicode[STIX]{x1D70E}^{\prime }$ which give the eigenvalues of the scattering operator. The eigenvalues are shown in figure 2(a) for $\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$ , corresponding to a flow correlation length $l_{c}\approx 180~\text{km}$ , with all other parameters as in figure 1. The most important eigenvalues are the two largest, $n=0$ and here $n=1$ , since they control the scattering and isotropisation time and length scales. These scales are displayed as functions of $\unicode[STIX]{x1D705}$ in figure 2(b). The figure shows that large-scale flows lead to rapid scattering but slow isotropisation. This can be easily understood: large-scale flows cause rapid energy transfers but, because of the localised nature of the scattering, these transfers are limited to waves of similar directions and a long time is needed for energy to be distributed near uniformly in the angular direction. (The large-scale flow regime can be tackled using ray tracing as has frequently been applied for deterministic flows (e.g. Rainville & Pinkel 2006; Chavanne et al. 2010). For weak random flows as assumed here, the ray equations can be analysed asymptotically using methods developed for noisy Hamiltonian systems (e.g. Bal, Komorowski & Ryzhik 2010) to show that the IT wavevector diffuses along the constant-frequency circle $|\boldsymbol{k}|=\text{const.}$ , consistent with the kinetic-equation description; see Müller (1976, 1977) for early treatments in this spirit.) Isotropisation is most effective when $\unicode[STIX]{x1D705}$ has an order of magnitude similar to $|\boldsymbol{k}|$ : for the chosen energy spectrum, isotropisation is fastest for $\unicode[STIX]{x1D705}\approx 6\times 10^{-5}~\text{m}^{-1}$ corresponding to a flow correlation length of approximately 50 km. Isotropisation slows down for larger values of $\unicode[STIX]{x1D705}$ simply because the total flow energy in the useful range $[0,2|\boldsymbol{k}|]$ decreases with $\unicode[STIX]{x1D705}$ .

Figure 2. (a) Eigenvalues $\unicode[STIX]{x1D706}_{n}$ of the scattering operator, given by (3.9), for the energy spectrum in (3.18) with $\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$ and $v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$ , and for an IT with $|\boldsymbol{k}|=3\times 10^{-5}~\text{m}^{-1}$ and $f=1.028\times 10^{-4}~\text{s}^{-1}$ . (b) Scattering and isotropisation length and time scales $L_{scat}$ , $L_{iso}$ , $T_{scat}$ and $T_{iso}$ as functions of the peak wavenumber $\unicode[STIX]{x1D705}$ , with all the other parameters as in (a). (c) As in (b) but as functions of $v_{\mathit{rms}}$ and for $\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$ . (d) As in (b) but as functions of latitude.

Figure 2(c) shows the scattering and isotropisation times and lengths as functions of $v_{\mathit{rms}}$ and for $\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$ . The dependence is simply in $v_{\mathit{rms}}^{-2}$ . The figure suggests that full isotropisation of ITs generated at localised topographical features is rare in the ocean since the length scales required exceed the basin scales even for strong flows. On the other hand, scattering is effective over much shorter spatial scales, of the order of a few hundreds of kilometres, and over time scales of a week or so, comparable to other dynamical time scales in the ocean. The conclusion, then, is that typical ITs are strongly influenced by the quasigeostrophic flow, though not to the extent that they become completely isotropic. As highlighted by Ward & Dewar (2010), the time scale of a week or so is shorter than the characteristic time scales of nonlinear wave–wave interactions except, perhaps, for the special case of parametric subharmonic instability at the critical latitude of $29^{\circ }$ (MacKinnon & Winters 2005). It is likely, then, that scattering by the geostrophic flow plays a more important role than wave–wave interactions in determining the characteristics of oceanic inertia–gravity waves. We should note, however, that the most energetic regions of the ocean, such as western boundary currents, where scattering is most effective, are also strongly inhomogeneous so that out theory does not apply in a strict sense.

Figure 2(d) explores the dependence of scattering on latitude. Latitude affects the cross-section $\unicode[STIX]{x1D70E}^{\prime }$ through $f$ and also through $|\boldsymbol{k}|$ if we consider a fixed frequency as is done here. In the ocean, different latitudes may also lead to different energy spectra. For simplicity, in plotting figure 2(d) we have taken the same spectrum for each latitude, keeping $\unicode[STIX]{x1D705}$ fixed. The figure shows that the scattering time increases with latitude. The isotropisation time, however, decreases with latitude with, as far as we can tell, no obvious interpretation; the scattering is determined as the difference between two eigenvalues and is hence difficult to intuit. The scattering and isotropisation lengths both decrease with latitude, partly as a result of a decrease of the group velocity. We note that Ward & Dewar (2010) conclude from simulations that scattering and isotropisation weaken with latitude, leading to longer propagation distances (see their figure 11). This apparent contradiction is likely resolved by the fact that their non-dimensional formulation implies that their energy spectrum also changes with latitude, keeping the energy in the range $[0,2|\boldsymbol{k}|]$ constant as $f$ changes. A general conclusion we can draw from the form of $\unicode[STIX]{x1D70E}^{\prime }$ and our parameter dependence study is the fact that the quasigeostrophic energy spectrum is the key factor determining the strength of the scattering.

4 Simulations

In this section we analyse numerical simulations of the linearised equivalent shallow-water system (2.2) and compare them with the theoretical predictions of the previous section and with direct simulations of the kinetic equation (2.6).

4.1 Shallow-water simulations

We solve (2.2) numerically, adding a harmonic forcing term to generate a coherent plane wave. The numerical scheme relies on pseudospectral and splitting methods: the terms independent of the background flow are integrated exactly in Fourier space, while the terms that depend on the background flow are integrated using an Euler scheme in physical space. The domain is a $7168~\text{km}\times 1024~\text{km}$ channel on an $f$ -plane centred at $45\,^{\circ }\text{N}$ . We use a spatial resolution of $1792\times 256$ , with that $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=4~\text{km}$ , with periodic boundary conditions in the $y$ -direction, and absorbing layers 30 grid points wide at each end of the domain in the $x$ -direction, and take time steps of $\unicode[STIX]{x0394}t=4000~\text{s}$ . We run an ensemble of 100 simulations with random realisations of the background flow in order to study statistics. Each simulation corresponds to 80 days, which is long enough to study isotropisation in a moderately energetic flow. A wavemaker forces an IT through a term of the form

(4.1) $$\begin{eqnarray}{\mathcal{F}}=A\sin (\unicode[STIX]{x1D6FA}t)\text{e}^{-(x-x_{0})^{2}/\unicode[STIX]{x1D6E5}^{2}},\end{eqnarray}$$

added to the continuity equation (cf. Ponte & Klein 2015). Here, $x_{0}=400~\text{km}$ is the position of the wavemaker in the $x$ -direction and $\unicode[STIX]{x1D6E5}=10~\text{km}$ its width; $\unicode[STIX]{x1D6FA}$ is the tidal frequency. The amplitude $A$ is arbitrary since we solve a linear system. The forcing is ramped up slowly to reach its maximum amplitude over approximately 1 week. The resulting plane waves that are generated have a wavelength of approximately 150 km, as expected for a first baroclinic mode wave at $45^{\circ }$ latitude.

Figure 3. (a) Sea-surface elevation $\unicode[STIX]{x1D702}$ in an equivalent shallow-water simulation of the mode-1 $M_{2}$ tide in a turbulent flow with $v_{rms}=0.25~\text{m}~\text{s}^{-1}$ at $45^{\circ }$ latitude. (b) Vorticity field of the turbulent flow.

For the background flow we take a homogeneous isotropic Gaussian random field, tapered in the region $0<x<1000~\text{km}$ so as not to interfere with the wavemaker. This field is generated numerically as a Fourier series with random coefficients. The energy spectrum is that in (3.18), with $v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$ and $\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$ leading to flows with a correlation length of approximately 200 km. The other parameters are those of the mode-1 $M_{2}$ tide at $45^{\circ }$ , as in § 3.2. The results presented below use a time-independent flow. We carried out additional simulations with slowly time-varying flows to confirm the theoretical prediction of appendix B that IT scattering is essentially unaffected by the time dependence of stationary random flows with time scales $O(Ro^{-1})$ longer than the IT period. Note that the modelling of the background flow by a Gaussian random field is a choice motivated by practicality and the fact that, according to our theory, the only statistical property of the flow that influences the scattering is its energy spectrum. An alternative would be to carry out a large number of quasigeostrophic simulations to generate an ensemble of flows with more realistic statistics. This would however be computationally expensive; it would also require great care to control the flow parameters and to ensure stationary statistics.

Figure 4. Energy density in the $\boldsymbol{k}$ plane in different regions along the channel, from an ensemble of 100 shallow-water simulations. The regions are five $1024~\text{km}\times 1024~\text{km}$ boxes centred about the midpoints $x=\{1000,2250,3500,4750,6000\}~\text{km}$ , as indicated in figure 3.

Figure 3(a) shows one realisation of the IT height field $\unicode[STIX]{x1D702}$ at $t=80$ days, when plane waves generated by the wavemaker (indicated by a dashed line) have propagated across the eddy field shown in the bottom panel. Close to the wavemaker the wave field has a plane-wave structure which becomes scrambled in appearance as the phases randomise due to flow scattering. In agreement with our scattering theory, the wave field retains a single length scale – the wavelength set by the tidal frequency – throughout its evolution: scattering does not lead to a scale cascade. This is consistent with the earlier simulation results of Ward & Dewar (2010), Ponte & Klein (2015), Dunphy et al. (2017) and Wagner et al. (2017), the latter two in a three-dimensional set-up. Note that, since we solve the linearised system, there is no harmonic generation and consequent formation of smaller wave scales as described by Ward & Dewar (2010). The eigenvalues $\unicode[STIX]{x1D706}_{n}$ for the IT and flow parameters chosen are those shown in figure 2(a) and correspond to $L_{scat}=420~\text{km}$ and $L_{iso}=5600~\text{km}$ . This is qualitatively consistent with the wave field in figure 3.

To assess our theoretical results in more detail, we need to estimate the energy density $a(\boldsymbol{x},\boldsymbol{k},t)$ from the simulations. To this end, we take Fourier transforms of the wave fields in five $1024~\text{km}\times 1024~\text{km}$ square boxes spanning the length of the domain. For each realisation of the flow, we compute the Fourier transform of $u,v$ and $\unicode[STIX]{x1D702}$ in each box at the end of the simulation, project onto the IT eigenmode then average over the ensemble to obtain an approximation $\tilde{a}(\boldsymbol{k},t)$ , say, of $a(\boldsymbol{x},\boldsymbol{k},t)$ in the box. The details of this are given in appendix C. The results are shown in figure 4. For the first box, located immediately to the right of the wavemaker, most of the energy is concentrated at the single point $\boldsymbol{k}=(k_{0},0)$ , where $k_{0}=\sqrt{(\unicode[STIX]{x1D714}^{2}-f^{2})/gh}$ , indicating a pure plane wave propagating to the right. (There is also a faint signal of left-propagating waves resulting from scattering at larger $x$ .) In the next boxes, energy spreads around the circle of constant radius $|\boldsymbol{k}|=k_{0}$ . The spectrum is in fact distributed over a finite-width annulus rather than a circle, as a result of off-resonant interactions between waves and flow.

We obtain a clearer view of the distribution of energy as a function of $\unicode[STIX]{x1D703}$ by integrating $a(\boldsymbol{x},\boldsymbol{k},t)$ in the 90 angular sectors $2(n-1)\unicode[STIX]{x03C0}/90\leqslant \unicode[STIX]{x1D703}\leqslant 2n\unicode[STIX]{x03C0}/90$ , $n=1,\ldots ,90$ . The results are shown in figure 5(a). They enable a better assessment of the validity of the length scale estimates $L_{\mathit{scat}}=420~\text{km}$ and $L_{\mathit{iso}}=5600~\text{km}$ . Note that these lengths should be measured from the point where the background flow starts, which is at approximately $x=1000~\text{km}$ , so that we would expect to see the fields isotropise at $x>6000~\text{km}$ along the channel, corresponding to the rightmost box of figure 4. We see however that the field is not fully isotropic in that region. An explanation is that the numerical simulation includes an absorbing layer near the right boundary of the domain, which allows energy to exit the channel but not to re-enter it. As a result, there are no left-propagating waves at the end of the channel. In addition, we emphasise that $L_{\mathit{iso}}$ is only an order-of-magnitude estimate which, by converting time scale into length scale using the group speed, ignores the directional properties of the transport of wave energy with the group velocity. We next go beyond this order-of-magnitude estimate and make direct predictions for the scattering by solving the kinetic equation numerically.

Figure 5. Estimated energy density $\tilde{a}$ as a function of $\unicode[STIX]{x1D703}$ from (a) the ensemble of shallow-water simulations, with the $x_{i}$ taken as the box midpoints, and (b) the kinetic-equation simulation. The parameters are those of figure 3 and the time corresponds to the end of the simulation.

4.2 Kinetic-equation simulations

We simulate the kinetic equation (2.6) under the assumption of homogeneity in the $y$ -direction, consistent with the periodic boundary conditions, and using the angle $\unicode[STIX]{x1D703}$ , with $\boldsymbol{k}=k_{0}(\cos \unicode[STIX]{x1D703},\sin \unicode[STIX]{x1D703})$ , as an independent variable. This reduces the number of independent variables to 3, with $a(x,\unicode[STIX]{x1D703},t)$ , so that the kinetic equation becomes

(4.2) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}a+c_{g}\cos \unicode[STIX]{x1D703}\unicode[STIX]{x2202}_{x}a=({\mathcal{L}}-\unicode[STIX]{x1D6F4})a+{\mathcal{F}}(x,\unicode[STIX]{x1D703}),\end{eqnarray}$$

where ${\mathcal{L}}$ and $\unicode[STIX]{x1D6F4}$ are given by (3.2) and (3.5), and ${\mathcal{F}}(x,\unicode[STIX]{x1D703})$ is a forcing term mimicking the wavemaker of the shallow-water simulations. We take ${\mathcal{F}}$ to be a Gaussian centred about $x=400~\text{km}$ , $\unicode[STIX]{x1D703}=0$ , with width parameters $\unicode[STIX]{x1D6E5}_{x}=40~\text{km}$ and $\unicode[STIX]{x1D6E5}_{\unicode[STIX]{x1D703}}=0.1$ , and an amplitude that is scaled to match the initial energy peak from the shallow-water simulation data.

We simulate (4.2) using a pseudospectral splitting method. The advection term is integrated using a semi-Lagrangian finite-difference scheme, and the scattering and forcing terms on the right-hand side are integrated exactly in Fourier space. The domain is $7168~\text{km}\times 2\unicode[STIX]{x03C0}$ , with a resolution of $1792\times 256$ and time step $\unicode[STIX]{x0394}t=900~\text{s}$ up to a final time of 80 days. We apply periodic boundary conditions in the $\unicode[STIX]{x1D703}$ -direction, and place absorbing layers 30 grid points wide at each end of the domain in the $x$ -direction, as in the shallow-water simulation. The evolution of $a(x,\unicode[STIX]{x1D703},t)$ is illustrated in figure 6. The wave energy, initially concentrated at $(x,\unicode[STIX]{x1D703})=(400,0)$ , gets advected by the group velocity in the $x$ -direction and spreads in the angular direction. Once energy reaches $|\unicode[STIX]{x1D703}|>\unicode[STIX]{x03C0}/2$ , it propagates to the left, leading to the weak signal for $k<0$ observed in figure 4(a).

Figure 6. Wave-energy density $a(x,\unicode[STIX]{x1D703},t)$ at $t=8$ , 16, 32 and 80 days obtained by solving the kinetic equations numerically with parameters matching those of figure 3. The $x$ -axis is in units of 1000 km.

At the end of the simulation, we evaluate $a(x,\unicode[STIX]{x1D703},t)$ at different points along the channel to get a set of curves $a(x_{i},\unicode[STIX]{x1D703},t=80~\text{days})$ , $1\leqslant i\leqslant 5$ . The points $x_{i}$ are taken to be spaced along the channel in the same way as the windows shown in figure 4. Note that we have to take into account the fact that there is no flow in the shallow-water simulations from the point of generation at $x=400~\text{km}$ to approximately $x=1000~\text{km}$ (see figure 3), whereas the kinetic equation assumes a flow is present throughout the domain. To resolve the discrepancy, the values of $x_{i}$ are taken 600 km less than the midpoints of the boxes used for the shallow-water simulations. The functions $a(x_{i},\unicode[STIX]{x1D703})$ are shown next to the equivalent shallow-water estimates in figure 5. There is a remarkably good agreement between the solution of the kinetic equation and the shallow-water simulation results. This demonstrates the value of the kinetic equation in predicting the generic properties of the scattering and their dependence on the various parameters in the problem.

Figure 7. Energy density from 10 realisations of shallow-water simulations evaluated (a) in the second box of figure 4, centred about $x_{2}=2250~\text{km}$ , and (b) in the fourth box, centred about $x_{4}=4750~\text{km}$ .

We emphasise that the kinetic equation yields only ensemble average predictions and cannot describe the details of the effect of a single flow realisation on the IT. To illustrate how the scattering fluctuates between realisations, we show in figure 7 the energy density $\tilde{a}$ estimated from single shallow-water simulations. While fluctuations can be large, the typical behaviour is a redistribution of energy in the angular direction that is well captured by the ensemble-averaged predictions. Furthermore, ergodicity implies that the ensemble average deductions of the kinetic equation apply accurately to quantities that are spatial averages over many eddy scales.

5 Discussion

This paper examines the scattering of oceanic ITs – or indeed of any inertia–gravity wave – caused by the turbulent mesoscale flow in which they propagate. Assuming that the flow is barotropic, weak (small Rossby number) and random with stationary and homogeneous statistics, we derive the kinetic equation (2.6) governing energy exchanges between waves travelling in different directions. A key outcome is the scattering cross-section (2.9), or (3.4) for an isotropic flow, which measures the energy transfer rate as a function of the wavevectors involved and other parameters. The scattering cross-section depends linearly on the energy spectrum of the flow.

The form of the scattering cross-section shows that the energy exchanges between waves are restricted to waves with the same frequency and hence the same horizontal wavenumber. Therefore, while scattering results in a complex random wave field, this field has a single spatial scale determined by the forcing frequency. This is obvious for a time-independent flow but perhaps less so for the time-dependent flows we consider for which it is a consequence of the statistical stationarity of the flows. Note that this does not imply that the wave field is completely phase locked in time: slow phase variations result from the interaction with a time-dependent flow (Ponte & Klein 2015; Dunphy et al. 2017), but these are not described by our analysis, which focuses on the wave amplitude as measured by the energy density $a(\boldsymbol{x},\boldsymbol{k},t)$ . It would be of interest to study the phase variations from the statistical viewpoint taken here.

For an isotropic flow, scattering leads to an equilibrium isotropic wave field over a time scale that we can estimate from the Fourier transform of the scattering cross-section (3.4). At equilibrium, and in the absence of spatial modulations, the wave energy density is $a(\boldsymbol{k})=a(|\boldsymbol{k}|)\propto \unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-k_{0})$ with $k_{0}=\sqrt{(\unicode[STIX]{x1D714}^{2}-f^{2})/gh}$ , corresponding to a correlation function $\int a(|\boldsymbol{k}|)\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\text{d}\boldsymbol{k}\propto \int \unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-k_{0})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\text{d}\boldsymbol{k}$ that is proportional to the Bessel function $\text{J}_{0}(k_{0}|\boldsymbol{x}|)$ . Thus, while the flow controls the speed of convergence towards the isotropic wave-energy distribution, it has no effect on the form of this distribution.

The time scale necessary for scattering to significantly alter the wave field is deduced from the scattering cross-section and found to be of the order of a few days to a week. This is short compared with time scales associated with nonlinear wave–wave interactions (see Ward & Dewar 2010), which raises the possibility that scattering is as crucial as the more widely considered wave–wave interactions in shaping the inertia–gravity wave spectrum in the ocean. Note however that, as our asymptotic treatment makes clear, the large time scale separation between inertia–gravity waves and the quasigeostrophic flow implies that scattering causes little frequency broadening and so cannot by itself explain the continuum of observed frequencies.

The conclusion that scattering simply relaxes wave energy towards an isotropic equilibrium depends crucially on the two-dimensional set-up implied by our assumption of barotropic flow. It holds because scattering redistributes wave energy in Fourier space over constant-frequency curves which, in this case, are just circles, hence compact. This picture changes radically in the presence of vertical shear since this causes energy exchanges between different vertical wavenumbers. The relevant constant-frequency set is then a cone in wavevector space. Because this cone is unbounded, no finite-energy equilibrium exists, and the energy of an initially plane wave can be expected to cascade to small scales, both horizontally and vertically and with a fixed aspect ratio, as it spreads on the cone. This scenario, already envisioned by Lelong & Riley (1991) and Bartello (1995), is potentially important for the dissipation of oceanic inertia–gravity waves, with implications for their impact on mixing and mesoscale dynamics and for the maintenance of balance. The framework we have adopted, which regards the flow as a prescribed random field and examines the wave statistics on the basis of a kinetic equation, generalises to the case of vertically sheared flow. It is well suited to describe and quantify the scale cascade that results from what is then a fully three-dimensional scattering. This is the subject of future work.