[
Abstract
The scattering of inertiagravity waves by largescale geostrophic turbulence in a rapidly rotating, strongly stratified fluid leads to the diffusion of wave energy on the constantfrequency cone in wavenumber space. We derive the corresponding diffusion equation and relate its diffusivity to the wave characteristics and the energy spectrum of the turbulent flow. We check the predictions of this equation against numerical simulations of the threedimensional Boussinesq equations in initialvalue and forced scenarios with horizontally isotropic wave and flow fields. In the forced case, wavenumber diffusion results in a wave energy spectrum consistent with asyetunexplained features of observed atmospheric and oceanic spectra.
Diffusion of inertiagravity waves]Diffusion of inertiagravity waves by geostrophic turbulence H. Kafiabad, M. A. C. Savva & J. Vanneste]Hossein Kafiabad, Miles A. C. Savva and Jacques Vanneste
1 Introduction
The dynamics of rotating stratified fluids, most notably the atmosphere and ocean, is characterised by the coexistence of vortical flow and inertiagravity waves (IGWs). These evolve independently at a linear level but interact to an increasing degree as flow strength and wave amplitude increase. In the weakly nonlinear regime, corresponding to small Rossby and/or Froude numbers, the flow has a ‘catalytic’ role, enabling the scattering of energy between IGWs through resonant triad interactions while remaining unaffected (Lelong & Riley 1991; Bartello 1995; Ward & Dewar 2010). The qualitative impact of this catalytic interaction has been considered: an isotropic turbulent flow causes the isotropisation of the IGW field (Lelong & Riley 1991; Savva & Vanneste 2018) and a cascade of wave energy to small scales (Bartello 1995; Waite & Bartello 2006a).
Here we provide a quantitative description by deriving a simplified model for the dynamics of IGWs in a lowRossbynumber, homogeneous and horizontally isotropic turbulent flow in geostrophic balance. The derivation (in §2) assumes linear IGWs with small spatial scales relative to the flow. It yields a diffusion equation that captures the spreading of IGWs in wavenumber space or, more precisely, on a cone in this space corresponding to fixedfrequency IGWs. The diffusivity components associated with radial and angular diffusion on the cone are obtained in closed forms involving the IGWs parameters and the energy spectrum of the geostrophic flow. Early versions were proposed by Müller & Olbers (1975) and Müller (1976, 1977).
We solve the diffusion equation for an initialvalue problem (§3) and a steady forced problem (§4), assuming horizontally isotropic IGW fields, and we test the results against numerical simulations of the threedimensional Boussinesq equations, finding good agreement in both cases. With forcing, the diffusion equation predicts a constantflux, steady energy spectrum scaling with wavenumber as which is realised numerically.
Our results are relevant to important open questions about the nature of submesoscale motion in the ocean and mesoscale motion in the atmosphere. Recent data analyses by Bühler et al. (2014) and Callies et al. (2014, 2016) led them to hypothesise these motions are dominated by almost linear IGWs. The prediction of a spectrum lends support to this hypothesis by identifying a robust mechanism – diffusion by turbulence – that produces a spectrum consistent with observations (see §4). As for the initialvalue predictions, they provide estimates for the time scale of the scale cascade of the IGWs that leads ultimately to their dissipation.
2 Diffusion in wavenumber space
We consider the dynamics of IGWs propagating in a turbulent flow of much larger spatial scale so that the WKB approximation applies. The distribution of wave energy in the phase space is then governed by the conservation
(2.0) 
of the action density . Here is the frequency, which sums the intrinsic frequency
(2.0) 
where are the Coriolis and buoyancy frequencies and is the angle between the wavevector and the vertical, and the Doppler shift , where is the flow velocity. Assuming that the flow is (i) weak enough that , (ii) evolving on a time scale much longer than , and (iii) well modelled by a homogeneous and stationary random field, we can approximate (2) by
(2.0) 
where is the intrinsic group velocity and a dependent diffusivity tensor (see Appendix A for a derivation). The righthand side of (2) captures the scattering of wave action that results from smallbutsustained random Doppler shifting by the flow; in the regime considered, this naturally leads to diffusion in space. In Cartesian coordinates, the diffusivity tensor takes the form
(2.0) 
where is the velocity correlation tensor, with denoting ensemble average, and summation over repeated indices is implied. An analogous expression was obtained by McComas & Bretherton (1977) in the context of wave–wave interactions in the induceddiffusion regime (see Müller et al. 1986, §5, for a review). Müller & Olbers (1975) and Müller (1976, 1977) discussed a flowinduced diffusivity that differs from (2) to account heuristically for wave–wave interactions and dissipation.
A key property of (2) is that : there is no diffusion in the direction of the group velocity . Since is perpendicular to constantfrequency surfaces, for the IGW dispersion relation (2) diffusion is restricted to the cones , see Fig. 1. This is because diffusion in space stems from resonanttriad interactions between two IGWs and one flow mode (often termed vortical or balanced mode); the flow is treated as a zerofrequency mode because it evolves slowly compared with , so the resonance condition implies that the interacting IGWs have the same frequency. The restriction to a single frequency means that wave action and wave energy only differ by a constant multiple and can be identified with one another.
We particularise (2) to IGWs and geostrophic flows using the dispersion relation (2) and the geostrophic balanced satisfied by the velocity in . It is natural to use spherical polar coordinates in space and a Fourier counterpart to in the form of the flow kinetic energy spectrum , which we assume to be horizontally isotropic so that it only depends on the horizontal and vertical wavenumbers and (for clarity we systematically use lowercase symbols for coordinates in the IGW wavenumber space and uppercase symbols for coordinates in the flow wavenumber space). Computations detailed in Appendix A then reduce (2) to
(2.0) 
under the further assumption of spatial homogeneity . This makes it plain that there is no diffusion in the direction of . Hence, , or equivalently , can be treated as a fixed parameter. The only nonzero components of the diffusivity tensor are given by \linenomath
(2.0a)  
(2.0b) 
where depends solely on , and .
Eqs. (2)–(2) provide a full description of the diffusion of IGW on the constantfrequency cone in space for a turbulent flow of given energy spectrum. In the angular direction, this diffusion leads to an isotropisation of the wave field with rate . In the radial direction, the diffusion leads to a forward cascade of the wave energy to high wavenumbers where it is efficiently dissipated by viscous processes. Note that wave energy remains confined to one nappe of the cone corresponding to either upward or downwardpropagating IGWs. This is only an approximation; exchanges between upward and downwardpropagating waves do occur, but they are asymptotically small and not captured by the WKB approximation. In what follows, we concentrate on radial diffusion by assuming wave statistics independent of , , leaving the study of horizontal isotropisation for future work.
3 Initialvalue problem
For , we rewrite (2) as
(3.0) 
where we have introduced and the independent parameter . The function is the IGW energy density in , with the energy contained within the interval . We solve (3) with initial condition corresponding to the excitation of IGWs with a single wavenumber . (The solution associated with arbitrary initial condition can be deduced by integration over .) We show in Appendix A.3 that
(3.0) 
where is a Bessel function of the first kind (DLMF 2018). The largetime behaviour of is readily deduced as away from an asymptotically small neighbourhood of . An inverse diffusion time scale can be read off from (3) as . Using (2) this can be written in the dimensionless form
(3.0) 
where is a dimensionless ‘geometric’ factor that depends only on and the shape (but not the magnitude) of the flow kineticenergy spectrum and is a flow Rossby number. The typical horizontal and vertical inverse flow scales and are assumed to be related by . Eq. (3) captures the dependence of the diffusion time scale on the Rossby number and on the scale separation between IGWs and flow. The diffusion approximation requires in addition to the WKB conditions and .
(a)  (b) 
(c)  (d) 
We verify the solution of (3) against simulations of the threedimensional nonhydrostatic Boussinesq equations. These are solved using a code adapted from that in Waite & Bartello (2006b) which relies on a dealiased pseudospectral method and a thirdorder Adams–Bashforth scheme with timestep . The triplyperiodic domain, in the scale coordinates , is discretised uniformly with grid points. A hyperdissipation of the form , with , is employed in the momentum and density equations. We take , a representative value of middepth ocean stratification. The initial condition is the superposition of a turbulent flow, obtained by running a quasigeostrophic model to a statistically stationary state, and IGWs. The initial spectrum of the flow peaks at and has an inertial subrange scaling approximately as and . This spectrum evolves slowly over the IGWdiffusion timescale, and an average is used to calculate in (2), and hence in (3). We report experiments with the two Rossby numbers (or for the alternative Rossby numbers based on the vertical vorticity ), and the two IGW frequencies . Upwardpropagating IGWs are initialised as a ring in space with , , random phases, and an initial kinetic energy . The IGW spectrum is computed following the normalmode decomposition of Bartello (1995).
Fig. 1, obtained for the lower and , illustrates the confinement of wave energy on the constantfrequency cone, one of the keys to the validity of the diffusion approximation. The confinement is of course not perfect and some energy appears around the cones associated with the harmonic frequencies and . Fig. 2 shows the evolution of for the four sets of values of . The numerical results are compared with the predictions of the diffusion equation obtained by solving (3) initialised with the form of extracted from the simulation after an adjustement time . This procedure accounts for the fact that the diffusion equation (2) is only valid after an adjustment period, requiring , the time to traverse typical eddies at the IGW group speed (cf. Müller et al. 1986, §5). The agreement between the numerical simulation and the diffusion approximation is remarkable considering the complexity of the full Boussinesq dynamics and the moderate separation of scales between IGWs and flow. As the diffusion approximation predicts, the simulations with different Rossby numbers behave similarly when is scaled suitably. Scattering from upwardpropagating to downwardpropagating IGWs, neglected in the diffusion approximation, occurs; it is more substantial for the larger because the two nappes of the constantfrequency cones are closer together, facilitating energy transfers.
4 Forced response and observed ocean and atmosphere spectra
We now turn to the steady solution of (3) in the presence of a forcing of the form . Eq. (3) admits two steady solutions: the noflux solution and the constantflux solution . Matching these at yields the steady spectrum
(4.0) 
Note that for IGWs with a single frequency and correspondingly a single angle , the horizontal energy spectrum satisfies the same power laws as since
(4.0) 
using that the energy density in space is . Thus, (4) implies a horizontal spectrum at large . This remains true for a superposition of IGWs with different frequencies, corresponding to an integration over .

We confirm the prediction (4) by the simulation of the Boussinesq equations in the presence of forcing. In the simulation reported, all the specifications are the same as in §3 except for the initial condition, which is devoid of IGWs. Instead, an Ornstein–Uhlenbeck forcing with short correlation time (3 timesteps) is applied to the waves with (see Waite 2017). The forcing amplitude is adjusted so that the wave energy is about of the flow energy after reaching a stationary state. Fig. 3 shows the stationary spectra for and . For the low Rossby number, the prediction (4) is well borne out by the simulation results with a clean spectrum spanning nearly a decade from the forcing scale down to dissipation. For the high Rossby number, the spectrum shallows from wavenumber or so to take a shape consistent with . Several mechanisms may cause this shallowing: the Doppler term is not weak compared with the intrinsic IGW frequency, invalidating the diffusion approximation, or nonlinear wave–wave interactions or nonresonant wave–flow interactions become significant.
The prediction of a spectrum is significant in view of the ubiquity of this scaling in ocean and atmosphere observations. In the ocean, kinetic energy spectra show a dependence in the submesoscale range, say below 20 km, in regions of high mesoscale activity and in a larger range, below 200 km, in less active regions (see Callies & Ferrari (2013) for a comprehensive discussion). Recent analyses by Bühler et al. (2014) and Rocha et al. (2016) which separate the contribution of IGWs from that of geostrophic motion indicate that the IGW part of the spectrum follows a scaling in almost the entirety of its range. Our results above suggest that this may result from IGW energy diffusion by the geostrophic flow. Scales below 10 km or so are the realm of the Garrett & Munk (1972) spectrum, also associated with a dependence. While this spectrum is generally attributed to wave–wave interactions (e.g. Müller et al. 1986; Lvov et al. 2012), we hypothesise that interactions with the geostrophic flow play an important role.
In the atmosphere, similarly, there is a broad range of scales, from 500 km to 10 km, where the energy spectrum scales approximately as . This is the shallow, mesoscale part of the celebrated Nastrom & Gage (1985) spectrum, which is traditionally interpreted as a spectrum but is also consistent with . There is ongoing debate about the nature of this part of the spectrum: Callies et al. (2014, 2016) attribute it to nearly linear IGWs on the basis of their separation between IGWs and geostrophic motion, but this interpretation is controversial (see Li & Linborg (2018) for a recent critique). Callies et al. (2016) note that ‘the wave interpretation is … not inconsistent with the observed powerlaw spectra … but an explanation for the spectral shape is so far missing’. Our results provide a possible explanation. The shallowing of the spectrum at small scales in the high simulation in Fig. 3, which is reminiscent of observations, however suggests that accounting for nonlinear phenomena may be essential to explain atmospheric spectra in their entirety.
5 Discussion
This paper examines the impact of a turbulent geostrophic flow on the statistics of smallamplitude IGWs. This impact has received little attention compared with that paid to wave–wave interactions. Yet the timescale found for a substantial effect of the geostrophic flow, of the order of (see Fig. 2) corresponding to tens of days for ocean parameters, is similar to that of the fastest wave–wave interaction process (parametric subharmonic instability of internal tides at the critical latitude , MacKinnon & Winters (2005)) This confirms the conclusions of Ward & Dewar (2010) and Savva & Vanneste (2018) that scattering by the flow dominates over wave–wave interactions in many ocean circumstances. A similar conclusion has been drawn from numerical simulations (Waite & Bartello 2006a) and likely applies to the atmosphere.
The present paper focuses on the diffusive regime of IGW scattering that arises for weak flows and smallscale, linear IGWs. A remarkable feature of this regime is the prediction of a energy spectrum consistent with observations in both the ocean and atmosphere. When the assumption of small scales is relaxed, the wave energy obeys a kinetic equation generalising the equations obtained by Danioux & Vanneste (2016) and Savva & Vanneste (2018) in the case of inertial waves and IGWs in a barotropic flow. The kinetic equation captures the transfer of energy between upward and downwardpropagating IGWs which is negligible in the diffusive regime. The derivation and analysis of this equation are the subject of ongoing work. When the assumption of weak flow is relaxed, IGWs are in the eikonal regime considered by Henyey & Pomphrey (1983) in the context of wave–wave interactions (see also Müller et al. 1986, §5). It would be desirable to study the scattering by geostrophic flow in this regime.
Acknowledgements. This research is funded by the UK Natural Environment Research Council under the NSFGEONERC programme (grant NE/R006652/1). M.A.C.S was supported by the Maxwell Institute Graduate School in Analysis and its Applications, funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016508/01), the Scottish Funding Council, HeriotWatt University and the University of Edinburgh. This work used the ARCHER UK National Supercomputing Service.
A Derivation of the diffusion equation and of its solution
a.1 General wave systems
We introduce a small parameter in the action conservation (2) by writing the frequency as , indicating that the velocity field is weak enough for the intrinsic frequency to dominate over the Doppler shift. Defining slow time and spatial scales by and , we substitute the expansion
(A 0) 
into (2). The first nontrivial equation appears at and is given by
(A 0) 
using Cartesian components and implied summation. Assuming that the velocity field varies on the slow time scale only, the solution is given by
(A 0) 
Averaging the nextorder equation to eliminate the terms containing , we find
(A 0) 
since using incompressibility and spatial homogeneity. Substituting the limit of (A.1) as as appropriate for the slow dynamics, we obtain the diffusion equation
(A 0) 
with the diffusivity
(A 0) 
This can be written as (2) in terms of the correlation tensor or, alternatively, as
(A 0) 
in terms of the Fourier transform of . The diffusive approximation (A.1) can be justified rigorously (see Bal et al. 2010, §4.2, and references therein).
a.2 IGWs in quasigeostrophic flow
We particularise (A.1) to the IGW dispersion relation (2) and a velocity field of the form with the geostrophic streamfunction. We use the spherical polar coordinates for , with , and the corresponding unit vectors, and express the group velocity as
(A 0) 
The diffusivity can be written in the basis as
(A 0) 
where , , , and we have made use of the fact that to eliminate all components along .
With and the polar and azimuthal angles of the flow wavevector , we have \linenomath
(A 0) 
where . Hence the delta function in (A.1) can be written as
(A 0) 
where . We also note that
(A 0) 
where is the flow kinetic energy spectrum. We now introduce (A 0)–(A 0) into (A.1) projected onto and to compute the components of in (A.2). Assuming that the flow is isotropic in the horizontal so that is independent of , we obtain after some simplifications \linenomath
(A 0a)  
(A 0b) 
and . The form (2) follows by replacing the kineticenergy spectrum by its twodimensional counterpart and changing the integration variables from to , with .
a.3 Solution of (3)
Introducing a solution of the separable form , with a spectral parameter, into (3) leads to
(A 0) 
where the prime denotes derivative with respect to . Solutions bounded as are proportional to the Bessel function . The general solution of (3) follows as
(A 0) 
for an arbitrary function . Imposing the initial condition yields (3) on using the Besselfunction expansion of (DLMF 2018, Eq. 1.17.13).
 Bal et al. (2010) Bal, G., Komorowski, T. & Ryzhik, L. 2010 Kinetic limits for waves in a random medium. Kinet. Related Models 3, 529–644.
 Bartello (1995) Bartello, P. 1995 Geostrophic adjustment and inverse cascades in rotating stratified turbulence. J. Atmos. Sci. 52, 4410–4428.
 Bühler et al. (2014) Bühler, O., Callies, J. & Ferrari, R. 2014 Wave–vortex decomposition of onedimensional shiptrack data. J. Fluid Mech. 756, 1007–1026.
 Callies et al. (2016) Callies, J., Bühler, O. & Ferrari, R. 2016 The dynamics of mesoscale winds in the upper troposphere and lower stratosphere. J. Atmos. Sci. 73, 4853–4872.
 Callies & Ferrari (2013) Callies, J. & Ferrari, R. 2013 Interpreting energy and tracer spectra of upperocean turbulence in the submesoscale range (1Ð200 km). J. Phys. Oceanogr. 43, 2456–2474.
 Callies et al. (2014) Callies, J., Ferrari, R. & Bühler, O. 2014 Transition from geostrophic turbulence to inertia–gravity waves in the atmospheric energy spectrum. Proc. Nat. Acad. Sci. 111, 17033–17038.
 Danioux & Vanneste (2016) Danioux, E. & Vanneste, J. 2016 Propagation of nearinertial waves in random flows. Phys. Rev. Fluids 1, 0033701.
 DLMF (2018) DLMF 2018 NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.0.21 of 20181215, F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.
 Garrett & Munk (1972) Garrett, C. J. & Munk, W. 1972 Space–time scales of internal waves. Geophys. Fluid Dyn. 2, 255–264.
 Henyey & Pomphrey (1983) Henyey, F. S. & Pomphrey, N. 1983 Eikonal description of internal wave interactions: nondiffusive picture of “induced diffusion”. Dyn. Atmos. Oceans 7, 189–219.
 Lelong & Riley (1991) Lelong, M.P. & Riley, J. J. 1991 Internal wavevortical mode interactions in strongly stratified flows. J. Fluid Mech. 232, 1–19.
 Li & Linborg (2018) Li, Q. & Linborg, E. 2018 Weakly or strongly nonlinear mesoscale dynamics close to the tropopause? J. Atmos. Sci. 94, 177–220.
 Lvov et al. (2012) Lvov, Y. V., Polzin, K. L. & Yokoyama, N. 2012 Resonant and nearresonant internal wave interactions. J. Phys. Oceanogr. 42, 669–691.
 MacKinnon & Winters (2005) MacKinnon, J. A. & Winters, K. B. 2005 Subtropical catastrophe: Significant loss of lowmode tidal energy at . Geophys. Res. Lett. 32 (15), L15605.
 McComas & Bretherton (1977) McComas, C. H. & Bretherton, F. P. 1977 Nonlinear interactions of internal gravity waves. J. Geophys. Res. 82, 1397–1412.
 Müller (1976) Müller, P. 1976 On the diffusion of mass and momentum by internal gravity waves. J. Fluid Mech. 77, 789–823.
 Müller (1977) Müller, P. 1977 Spectral features of the energy transfers between a largerscale shear flow. Dynam. Atmos. Oceans 2, 49–72.
 Müller et al. (1986) Müller, P., Holloway, G., Henyey, F. & Pomphrey, N. 1986 Nonlinear interactions among internal gravity waves. Rev. Geophys. 24, 493–536.
 Müller & Olbers (1975) Müller, P. & Olbers, D. J. 1975 On the dynamics of internal waves in the deep ocean. J. Geophys. Res. 75, 3848–3860.
 Nastrom & Gage (1985) Nastrom, G. D. & Gage, K. S. 1985 A climatology of atmospheric wavenumber spectra of wind and temperature observed by commercial aircraft. J. Atmos. Sci. 42, 950–960.
 Rocha et al. (2016) Rocha, C. B., Chereskin, T. K. & Gille, S. T. 2016 Mesoscale to submesoscale wavenumber spectra in Drake passage. J. Phys. Oceanogr. 46, 601–620.
 Savva & Vanneste (2018) Savva, M. A. C. & Vanneste, J. 2018 Scattering of internal tides by barotropic quasigeostrophic flows. J. Fluid Mech. 856, 504–530.
 Waite (2017) Waite, M. L. 2017 Random forcing of geostrophic motion in rotating stratified turbulence. Phys. Fluids 29, 126602.
 Waite & Bartello (2006a) Waite, M. L. & Bartello, P. 2006a Stratified turbulence generated by internal gravity waves. J. Fluid Mech. 546, 313–Ð339.
 Waite & Bartello (2006b) Waite, M. L. & Bartello, P. 2006b The transition from geostrophic to stratified turbulence. J. Fluid Mech. 568, 89–Ð108.
 Ward & Dewar (2010) Ward, M. L. & Dewar, W. K. 2010 Scattering of gravity waves by potential vorticity in a shallowwater fluid. J. Fluid Mech. 663.