Scattering of internal tides by barotropic quasigeostrophic flows

Oceanic internal tides and other inertia-gravity waves propagate in an energetic turbulent flow whose lengthscales 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 wavefield. 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.


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 stem from their unique place at the low-frequency end of the inertia-gravity wave spectrum (Alford et al. 2016); ITs, in contrast, are generic midfrequency 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 Wagner et al. (2017) and Dunphy et al. (2017) for further background.
A key aspect of the interactions between quasigeostrophic turbulence and both nearinertial oscillations and ITs is that turbulence and waves share similar horizontal scales, of the order of 100 km. A consequence is that the 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 timescale 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). 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 et al. (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(x, 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 (x, k, t) shows that energy transfers are restricted to waves with the same frequency or, equivalently, the same wavenumber |k|. 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(x, |k|).
Our theory complements that developed by Ward & Dewar (2010), extending their analysis of the scattering induced by a single Fourier mode of geostrophic flow to the full spectrum of realistic turbulence, and shifting from a deterministic to a statistical treatment in the process. It generalises the theory developed by Danioux & Vanneste (2016) for near-inertial oscillations to inertia-gravity waves of arbitrary frequencies.
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 wavefields (i.e. such that ∇ x a = 0) and in particular to the isotropisation of an initially plane wave examined numerically by Ward & Dewar (2010) in a shallowwater setup. The kinetic equation applies to more general, non-homogeneous situations in which the energy density is modulated spatially (∇ x a = 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. We carry out shallow-water simulations in a setup 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 timescale 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 A. 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 lengthscales 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.

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 U = (U, V, 0) = (−∂ y ψ, ∂ x ψ, 0) in terms of a z-independent streamfunction ψ. With these assumptions, the linearised hydrostatic-Boussinesq equations read where (u, w) denotes the IT velocity,ẑ 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 ∇ = (∂ x , ∂ y , 0) for the horizontal gradient throughout.
Since the background flow is barotropic, we project (2.1) onto a vertical-mode basis according to where the F n are eigenfunctions of the Sturm-Liouville problem and H is the ocean depth (Olbers et al. 2012). The eigenvalues r n are the Rossby radii of deformation. Orthonormality implies Substituting (2.2) into (2.1) leads to a system for the amplitudes u n , v n , etc. of each baroclinic mode n. Defining the equivalent depth h n = f 2 r 2 n /g and the equivalent surface height η n = b n /g, we rewrite this system in the shallow-water-like form ∂ t u n + u n · ∇U + U · ∇u n + fẑ × u n = −g∇η n , ∂ t η n + U · ∇η n + h n ∇ · u n = 0.
(2.3) 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. We drop the subscript n from this point onwards.
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. 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.

Derivation of the kinetic equation
We study IT scattering in the distinguished regime where the spatial scale of the flow, L * say, is of the same order as the wavelength, that is, |k|L * = O(1), where k is the IT horizontal wavevector. The assumption of a geostrophic flow requires a small Rossby number Ro = U * /(f L * ) 1, where U * is a typical flow velocity; in turn, this implies that the background flow velocities are small compared with the wave phase speed ω/|k|, where is the IT frequency, since U * /(ω/|k|) = O(U * /(ωL * )) = O(Ro). With the flow timescale T * taken as the natural advective timescale L * /U * , this also implies that the flow evolves slowly compared with the IT timescale since ωT * = O(Ro) 1. We further assume that, while the IT phases vary over the lengthscale |k| −1 , their amplitudes vary over a much larger scale ε −1 |k|, where ε 1. We adopt the scaling ε = O(Ro 2 ) which ensures that transport and scattering affect the wave field at the same order.
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 where φ = (u, v, η) 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 k, simply given by φ(x, t)φ * (x, t), is.
The equation derived by Ryzhik et al. (1996) governs the evolution of a scalar amplitude a(x, k, t) which appears naturally in an eigenvector decomposition of the matrix W (x, k, t) (see Appendix A for details). Physically, this amplitude is interpreted as a wavevector-resolving energy density, related to the (leading-order) energy density of the system by In Appendix A, we show that a(x, k, t) satisfies the kinetic equation Here ω is determined by the IT dispersion relation (2.4), so that ∇ k ω is the group velocity and the left-hand side of (2.6) represents the familiar wave transport. (The term −∇ x ω·∇ k a would be added if ω depended explicitly on x.) The right-hand side represents wave scattering by the background flow. The first term, given by quantifies the transfers of energy from all wavevectors k into wavevector k that result from interactions with the background flow; the second term, where is the total scattering cross section, quantifies the energy lost by wavevector k to all other wavevectors. The function σ(k, k ) 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 k to k at a position x in space. We obtain it in the form where k × k = k 1 k 2 − k 2 k 1 , and E is the energy spectrum of the flow. We note that σ(k, k ) is real, positive, and symmetric with respect to the exchange between k and k . In Appendix A, we also show that these properties ensure conservation of the leading-order energy density (2.5): is the leading-order energy flux (see (A 24)). The presence of the factor δ(|k| − |k |) 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 |k| = ((ω 2 − f 2 )/(gh)) 1/2 in the wavevector plane. This can be related to the observation that the background flow only enters σ(k, k ) through its energy spectrum E, which, for the statisically stationary flows considered, is time independent.

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, E(k) = E(|k|). This motivates the use of polar coordinates for the wavevector, such that where θ is the angle between k and k . The change of coordinates reduces the scattering operator (2.7) to Note that we have used the evenness of σ in θ to rewrite (3.1) as a convolution. Note also that σ is independent of the direction θ of k because the scattering process is rotationally invariant. The scattering cross section (2.9) can be written explicitly as 2) where we have removed the prime from θ for convenience. Similarly, the total scattering cross section (2.8) reduces to (3.3) In the limit ω → f corresponding to near-inertial waves the cross section reduces to which recovers the result obtained by Danioux & Vanneste (2016) starting from the Young-Ben Jelloul model. With the scattering cross section (3.2), the scattering operator (3.1) can be diagonalised using a Fourier series, or more precisely a cosine series since a is even in θ. Denoting the cosine transform by a hat, with we find that ( La) n = λ n a n , n = 0, 1, · · · , with the eigenvalues Fourier transforming the kinetic equation (2.6) then gives We note that σ in (3.2) is a non-negative and smooth function of θ; as a result, the eigenvalues (3.4) satisfy λ 0 = Σ(|k|) and |λ n 1 | < λ 0 .
Thus the scattering term on the right-hand side of (3.5) vanishes for n = 0 and represents a damping for n 1. The implications are clearly seen for a wave field that is spatially homogeneous, that is, with ∇ x a = 0: the solution of (2.6), with initial condition This describes the relaxation of the solution towards a stationary, isotropic (θindependent) solution, since 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. We can identify two timescales for the scattering process. First, the scattering time esimates the time over which energy concentrated at k in spectral space is reduced by a factor of e −1 while converted to waves with other wavevectors. In other words, it is the timescale over which scattering effects become significant. Second, the timescale for convergence to an isotropic wavefield is given by This is the time for the last surviving anisotropic (i.e. n = 0) Fourier mode to decay by a factor of e −1 . Scattering lengthscales associated with the timescales (3.6) and (3.7) can be defined as

Predicted behaviour
In this section, we use the time-and lengthscales (3.6)-(3.8) to examine how the scattering depends on the Coriolis parameter f , and on the strength and horizontal scales of the eddies as encoded in E. Since we focus on ITs, we regard the frequency ω as fixed and deduce |k| from the dispersion relation (2.4). We test some of our predictions against numerical simulations in §4.
We assume an isotropic energy spectrum of the form (3.9) This depends on two parameters: κ, a peak wavenumber which sets the dominant lengthscale of the flow, and the root-mean-square velocity defined by v 2 rms = ∞ 0 E d|k|. The constants c 1 and c 2 are determined by κ and v rms and the requirement of continuity at |k| = κ. In practice, we choose κ so that the correlation length l c = π/k c , where k c := |k| E dk/ E dk, is similar to the wavelength of the IT. Although quasigeostrophic theory predicts a kinetic energy spectrum that decays as |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 et al. 1994;Bartello 1995;Kafiabad & Bartello 2016). This motivates the form of (3.9) as representative of balanced geostrophic eddy fields in the ocean. Note that we have chosen an energy spectrum with non-zero energy for |k| < κ. This turns out to be important: because of the factor E(2|k sin(θ /2)|), the scattering cross section (3.2) depends on the energy in the range [0, 2|k|]. In particular, if |k| κ/2, a spectrum with E = 0 for |k| κ would lead to no scattering. Figure 1 shows the scattering cross section σ in (3.2) as a function of θ and κ for v rms = 0.25 m s −1 and for f = 1.028 × 10 −4 s −1 corresponding to 45 • latitude. The equivalent depth is set to h = 1.2 m, as appropriate for the first baroclinic mode (Olbers et al. 2012), and the frequency to ω = 2π/12.42 hours, corresponding to the M 2 tide. The horizontal wavenumber is then |k| = 3 × 10 −5 m −1 corresponding to a wavelength of about 200 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 κ, corresponding to flows with typical scales much larger than the IT wavelength (left panel), when the values of σ are also the largest. For larger values of κ, that is, for flows with smaller scales, the energy transfers are slower, but less localised in the angular direction (right panel).
The net effect of the scattering depends on both the value of σ at fixed θ and the range of θ where σ is substantial; it is best measured by the scattering and isotropisation time-and lengthscales introduced in (3.6)-(3.8). These scales are deduced from the cosine transform of σ which give the eigenvalues of the scattering operator. The eigenvalues are shown in Figure 2a for κ = 1.45 × 10 −5 m −1 , corresponding to a flow correlation length l c ≈ 180 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 lengthscales. These scales are displayed as functions of κ in Figure 2b. 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. Isotropisation is most effective when κ has an order of magnitude similar to |k|: for the chosen energy spectrum, isotropisation is fastest for κ ≈ 6 × 10 −5 m −1 corresponding to a flow correlation length of about 50 km. Isotropisation slows down for larger values of κ simply because the total flow energy in the useful range [0, 2|k|] decreases with κ. Figure 2c shows the scattering and isotropisation times and lengths as functions of v rms and for κ = 1.45 × 10 −5 m −1 . The dependence is simply in v −2 rms . The figure suggests that full isotropisation of ITs generated at localised topographical features is rare in the ocean since the lengthscales 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 kilometers, 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 timescale of a week or so is shorter than the characteristic timescales of nonlinear wave-wave interactions except, perhaps, for the special case of parametric subharmonic instability at the critical latitude of 29 • . 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. Figure 2d explores the dependence of scattering on latitude. Latitude affects the cross section σ through f and also through |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 2d we have taken the same spectrum for each latitude, keeping κ 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|k|] constant as f changes. A general conclusion we can draw from the form of σ and our parameter-dependence study is the fact that the quasigeostrophic energy spectrum is the key factor determining the strength of the scattering.

Simulations
In this section we analyse numerical simulations of the linearised equivalent shallowwater system (2.3) and compare them with the theoretical predictions of the previous section and with direct simulations of the kinetic equation (2.6).

Shallow-water simulations
We solve (2.3) 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 km × 1024 km channel on an f -plane centred at 45°N. We use a spatial resolution of 1792 × 256 with periodic boundary conditions in the y-direction, and absorbing layers 30-gridpoints wide at each end of the domain in the xdirection, and take timesteps of ∆t = 4000 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  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. For the background flow we take a homogeneous isotropic Gaussian random field, tapered in the region 0 < x < 1000 km so as not to interfere with the wavemaker. The energy spectrum is that in (3.9), with v rms = 0.25 ms −1 and κ = 1.45 × 10 −5 m −1 leading to flows with a correlation length of about 150 km. The other parameters are those of the mode-1 M 2 tide at 45 • , 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 the Appendix that IT scattering is essentially unaffected by the time dependence of stationary random flows with timescales 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.
The top panel of Figure 3 shows one realisation of the IT height field η 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. As predicted by our scattering theory, the wave field retains a single lengthscale -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), Wagner et al. (2017), Ponte & Klein (2015) and Dunphy et al. (2017), the latter two in a three-dimensional setup. 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 λ n for the IT and flow parameters chosen are those shown in Figure 2(a) and correspond to L scat = 420 km and L iso = 5, 600 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(x, k, t) from the simulations. To this end, we take Fourier transforms of the wave fields in five 1024 km × 1024 km square boxes spanning the length of the domain. For each realisation of the flow, we compute the Fourier transform of u, v and η in each box at the end of the simulation, project onto the IT eigenmode then average over the ensemble to obtain an approximationã(k, t), say, of a(x, k, t) in the box. The details of this are given in Appendix B. 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 k = (k 0 , 0), where k 0 = (ω 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 |k| = k 0 . There is also a weak spread transverse to the circle, a result of off-resonant interactions between waves and flow.
We obtain a clearer view of the distribution of energy as a function of θ by integrating a(x, k, t) in the 90 angular sectors 2(n − 1)π/90 θ 2nπ/90, n = 1, · · · , 90. The results are shown in Figure 5a. They enable a better assessment of the validity of the lengthscale estimates L scat = 420 km and L iso = 5600 km. Note that these lengths should be measured from the point where the background flow starts, which is at approximately x = 1000 km, so that we would expect to see the fields isotropise at x > 6000 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 iso is only an order-of-magnitude estimate which, by converting timescale into lengthscale 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.

Kinetic equation simulations
We simulate the kinetic equation (2.6) under the assumption of homogeneity in the ydirection, consistent with the periodic boundary conditions, and using the angle θ, with k = k 0 (cos θ, sin θ), as independent variable. This reduces the number of independent variables to 3, with a(x, θ, t), so that the kinetic equation becomes where L and Σ are given by (3.1) and (3.3), and F(x, θ) is a forcing term mimicking the wavemaker of the shallow-water simulations. We take F to be a Gaussian centred about x = 400 km, θ = 0, with width parameters ∆ x = 40 km and ∆ θ = 0.1, and an amplitude that is scaled to match the initial energy peak from the shallow-water simulation data. We simulate (4.1) 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 km × 2π, with a resolution of 1792 × 256, and time step ∆t = 900 s up to a final time of 80 days. We apply periodic boundary conditions in the θ-direction, and place absorbing layers 30-gridpoints wide at each end of the domain in the x-direction, as in the shallow-water simulation. The evolution of a(x, θ, t) is illustrated in Figure 6. The wave energy, initially concentrated at (x, θ) = (400, 0), gets advected by the group velocity in the x-direction and spreads in the angular direction. Once it reaches |θ| > π/2, energy propagates to the left, leading to the weak signal for k < 0 observed in the first panel of Figure 4. At the end of the simulation, we evaluate a(x, θ, t) at different points along the channel to get a set of curves a(x i , θ, t = 80 days), 1 i 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 km to approximately x = 1000 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 , θ) 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.
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ã 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.

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.2) 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 wavefield, this field has a single spatial scale determined by the forcing frequency. This is obvious for a timeindependent 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 wavefield 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(x, 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 wavefield over a time scale that we can estimate from the Fourier transform of the scattering cross section (3.2). At equilibrium, and in the absence of spatial modulations, the wave energy density is a(k) = a(|k|) ∝ δ(|k| − k 0 ) with k 0 = (ω 2 − f 2 )/gh, corresponding to a correlation function proportional to the Bessel function J 0 (k 0 |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 setup 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.

Acknowledgements. This research is funded by the UK Natural Environment Research
Council under the NSFGEO-NERC programme (grant NE/R006652/1). M.A.C.S was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016508/01), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.

Appendix A. Kinetic equation derivation
In this appendix we apply the formalism of Ryzhik et al. (1996) to derive a kinetic equation for the equivalent shallow-water system (2.3) (see also Bal et al. (2010) and Powell & Vanneste (2005)). This formalism exploits the weakness of the flow as measured by the small Rossby number, as well as the scale separation between wavelength and flow scale on the one hand, and the scale of variation of the wave amplitude on the other hand. The associated small parameter, ε, is assumed to satisfy ε = O(Ro 2 ).

A.1. Shallow-water equations and scaling regime
We start by rewriting (2.3) in the more abstract notation where the vector φ(x, t) = (u, v, η) T groups the dynamical variables and we have introduced the matrix operators In (A 1), we have made the relative importance of the various terms explicit by scaling them with the relevant power of the small parameter ε. For convenience, we keep the equations in their dimensional form and treat ε as a bookkeeping parameter that can be set to 1 at the end of the calculation. Note that we could use Ro as an alternative to ε in this bookeeping role and, in this manner, avoid fractional powers; using ε keeps our derivation in closer correspondence to that of Ryzhik et al. (1996). Note also that the dependence of N on √ εt arises through the slow time dependence of ψ. The depth-averaged energy density for the shallow-water system without background flow is given by (A 3) The matrix L in (A 2) is known as the dispersion matrix, since its eigenvalues give the dispersion relation. Taking the definition in (A 2) and replacing ∂ x by ik, we find that the solutions to the eigenvalue equation L( ik)b(k) = iω(k)b(k) are given by which is the usual dispersion relation for the rotating shallow-water model. The three eigenvectors are given by These are orthonormal in the sense that The zero eigenvalue ω 0 corresponds to the vortical mode of the system. Since this mode is accounted for by the quasigeostrophic background flow, it does not appear in any of the following derivation. In addition to the right eigenvectors b i , we also use the left eigenvectors c i which satisfy

A.2. Scaled Wigner transform
It is convenient to rescale the space and time coordinates as (x, t) → (x/ε, t/ε), and define φ ε (x, t) := φ(x/ε, t/ε). Under this rescaling, the equations of motion (A 1) take the form In the new coordinates, the wavevector is O(ε −1 ). To account for this and to obtain a well-defined Wigner transform in the limit ε → 0, we use the rescaled definition We note that this function remains positive semi-definite as ε goes to zero, which is important as it is closely related to the energy of the system. Indeed it can readily be shown that and so the energy density is found from the Wigner transform of the wave fields as A useful alternative to (A 6) expresses the Wigner transform in terms of Fourier transforms as

A.3. Wigner-transform evolution equation
We obtain an evolution equation for the Wigner transform by differentiating (A 6) with respect to t and substituting (A 5). This gives where c.c. denotes the complex conjugate of the term preceding it. Rewriting φ ε and N in terms of their Fourier transforms and making use of (A 8), we find that Scattering effects due to interactions of waves with the background flow are controlled by the third term, P ε W ε . To derive it, we have introduced the fast-space variable ξ := x/ε and the Fourier transform N of N defined by N (ξ, ·) = R 2 e − ip·ξ N (p, ·) dp.
We now derive the asymptotic limit of (A 9) using a multiscale expansion. Defining the intermediate time variable τ := t/ √ ε to cater for the time dependence of the streamfunction, we expand ) where we have anticipated that the leading-order term depends on the slow variables only. The differential operators are then expanded as where x and ξ, t and τ are treated as independent variables, leading to the expansion of the operators in (A 9). It turns out that only the leading order term P 0 is required for the derivation of the kinetic equation. The operators in (A 11) can be written explicitly through their action on an arbitrary function Z(x, ξ, k): e − ip·ξ N (p, i(k + p 2 ) + 1 2 ∂ ξ , τ )Z(x, ξ, k + p 2 ) dp + c.c.
We have decorated the operators with a tilde to highlight the presence of ∂ ξ in their definition; the tildes will be removed whenever this dependence disappears. Substituting the operators into (A 9) gives us the evolution equation for the Wigner function as 1 Introducing the expansion (A 10) then leads to a hierarchy of equations to be solved at each order in ε.
The leading-order equation is The eigenvalues of L are purely imaginary, so this equation is satisfied by taking W (0) in the form of a linear combination of the eigenvectors of the dispersion matrix. Defining the matrices , the leading order Wigner function is thus given by The so-far undetermined amplitudes a j (x, k, t) are real because the Wigner function is Hermitian.
At O(ε −1/2 ), we find where we have used that ∂ τ W (0) = 0. To solve (A 15), we rewrite W (1) in terms of its Fourier transform with respect to ξ, Substituting this into (A 15) yields where we have suppressed dependencies on x, t and τ for conciseness. Following Ryzhik et al. (1996), we have introduced a regularisation parameter θ > 0 which will be taken to zero at a later stage.
Since the Wigner transform is Hermitian, W (1) (p, k) = W (1) * (−p, k). Using this, expanding W (0) according to (A 14), and pre-and post-multiplying the resulting expression by c n (k − p/2) and c * m (k + p/2) (with c n the left eigenvector defined in (A 4)) gives It is convenient to extract the linear dependence of N on the streamfunction by defining a matrix U (p, iq) such that We now decompose W (1) using the vectors b i (k), which form a complete basis, as Using this along with the orthonormality of the eigenvectors we finally write the solution where we have taken into account that ψ(p) = ψ * (−p). We note that this solution shows W (1) is linear in the random field ψ.
The slow evolution of the leading-order Wigner function W (0) is controlled by the O(1) term in the expansion of (A 12), given by We assume that the random streamfunction is a stationary process in τ and homogeneous in ξ, with zero mean, ψ(ξ, τ ) = 0, and covariance where · denotes an ensemble average, or equivalently an average over ξ. In terms of Fourier transforms, this implies that where the streamfunction power spectrum R is the Fourier transform of R. The more familiar kinetic energy spectrum is then We now take the average of (A 17). The slow time derivative term on the right-hand side disappears since W (1) = 0. Since ∂ ξ W (2) = 0, Q 0 W (2) = Q 0 W (2) , where the removal of the tilde corresponds to setting ∂ ξ to 0 in Q 0 . This leads to an inhomogeneous version of (A 13). The matrix Q 0 has a non-trivial null space, spanned by the matrices B j (k); the righthand side of (A 20) must therefore satisfy a solvability condition. Since iQ 0 is self-adjoint with respect to the matrix inner product ⟪X, Y ⟫ := tr(M X * M Y ), this condition is obtained by applying ⟪B j (k), ·⟫ to (A 20). We deal with the resulting terms one by one. First, by orthogonality and (A 14) we have x, k, t).
In order to evaluate the remaining term, we note that using (A 16) and (A 18), we have N αβ (p, iq) N γδ (p , iq ) = U αβ (p, iq) U γδ (p , iq ) R(p)δ(p + p ), (A 21) where Greek indices are used for matrix elements to make the following derivation clearer, and summation over repeated Greek indices is implied. Expanding all terms, and making use of the delta function in (A 21), we have µν (p , k + p 2 ) dp dp + c.c.
We find that α(k, k ) = k × k ω 2 |k| 2 (f 2 + ω 2 )k · k − f 2 |k| 2 , and β(k, k ) = f ω ω 2 |k| 2 |k × k | 2 + k · k (|k| 2 − k · k ) , and we see that α(k, k ) = −α(k , k) and β(k, k ) = β(k , k). Using these symmetries makes it straightforward to show that We may thus define the differential scattering cross section as such that the amplitudes satisfy the kinetic equation Expanding the terms in (A 22), replacing the power spectrum with the energy spectrum using (A 19), and expressing the delta function instead in terms of wave vectors, we find that the differential scattering cross section for the shallow-water system (2.3) is given by (2.9). The conservation of the leading-order energy is established by noting that the equivalent of (A 7) for the scaled Wigner function is Integrating (A 23) with respect to k and noting that the right-hand side vanishes because of the symmetry σ(k, k ) = σ(k , k), we find leading order-energy conservation with the leading-order energy flux

Appendix B. Projection of simulation data onto modes
We describe how the energy density a(x, k, t) can be estimated from local Fourier transforms of the wave fields. Using the Fourier representation of the Wigner function given in (A 8) (with ε = 1), it is easily verified that the projection property R 2 W (x, −k, t) dx = | φ(k)| 2 holds. In order to discriminate between the energy contributions from the different modes, we expand the fields in the eigenvectors basis according to φ(k, t) = j=± A (j) (k, t)b j (k), so that the energy for each wavenumber is given bỹ with · , · M as defined in (A 3). Orthonomality of the eigenvectors means that we can extract the modal energy contributions by projection to find In terms of the energy density, the leading-order energy is given bỹ Thus, we may track the energy from the '+' mode by projecting the Fourier transform of the wave fields: with · denoting the ensemble average. This relates to the leading-order energy density of a single wave mode to the sea-surface height.