Diffusion of inertia-gravity waves by geostrophic turbulence

TL;DR

This paper derives a diffusion equation describing how inertia-gravity waves scatter and spread in wavenumber space due to geostrophic turbulence, supported by numerical simulations and relevant to atmospheric and oceanic spectra.

Contribution

It introduces a new diffusion model for inertia-gravity wave scattering by geostrophic turbulence and validates it with numerical simulations.

Findings

Wave energy diffuses on the constant-frequency cone in wavenumber space.

The model predicts a $k^{-2}$ energy spectrum in forced turbulence scenarios.

Results align with observed atmospheric and oceanic spectral features.

Abstract

The scattering of inertia-gravity waves by large-scale geostrophic turbulence in a rapidly rotating, strongly stratified fluid leads to the diffusion of wave energy on the constant-frequency 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 three-dimensional Boussinesq equations in initial-value and forced scenarios with horizontally isotropic wave and flow fields. In the forced case, wavenumber diffusion results in a $k^{-2}$ wave energy spectrum consistent with as-yet-unexplained features of observed atmospheric and oceanic spectra.

Equations

$${\partial }_{t}a+{\nabla }_{\mathbf{k}}\Omega \cdot {\nabla }_{\mathbf{x}}a-{\nabla }_{\mathbf{x}}\Omega \cdot {\nabla }_{\mathbf{k}}a=0$$

$$\omega =\sqrt{f^2{\cos }^2\theta +N^2{\sin }^2\theta },$$

$${\partial }_{t}a+\mathbf{c}\cdot {\nabla }_{\mathbf{x}}a={\nabla }_{\mathbf{k}}\cdot \left(\text{mathsfbi}D\cdot {\nabla }_{\mathbf{k}}a\right),$$

$$D_{ij}(\mathbf{k})=-\frac{1}{2}{k}_m{k}_n{\int }_{-\infty }^{\infty }\frac{{\partial }^2{\Pi }_{mn}}{\partial x_i\partial x_j}(\mathbf{c}(\mathbf{k})s)\mathrm{d}s,$$

$$D_{ij}(\mathbf{k})\cdot c_j(\mathbf{k})=-\frac{1}{2}{k}_m{k}_n{\int }_{-\infty }^{\infty }\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{\partial {\Pi }_{mn}}{\partial x_i}(\mathbf{c}(\mathbf{k})s)\right)\mathrm{d}s=0.$$

$${\partial }_{t}a=\frac{1}{k^2}{\partial }_k(k^2{D}_{kk}{\partial }_{k}a)+\frac{D_{\varphi \varphi }}{k^2{\sin }^2\theta }{\partial }_{\varphi \varphi }a,$$

$$D_{kk}$$

$$\frac{D_{\varphi \varphi }}{{\sin }^2\theta }$$

$$B=\frac{\omega {\sin }^2\theta }{4{\pi }^3(N^2-f^2)|{\cos }^5\theta |}$$

$${\partial }_{t}e={\partial }_k\left(Q{k}^5{\partial }_k\left(k^{-2}e\right)\right),$$

$$e(k,t)=\frac{1}{2}{k}_{\ast }^{-2}{\int }_0^{\infty }{J}_4(k^{-1/2}\lambda )J_4(k_{\ast }^{-1/2}\lambda ){\mathrm{e}}^{-Q{\lambda }^{2}t/4}\lambda \mathrm{d}\lambda ,$$

$$\frac{Q{k}_{\ast }}{\omega }=\gamma \frac{N^2}{N^2-f^2}\frac{k_{\ast }}{K_{\mathrm{h}\ast }}\mathrm{R}{\mathrm{o}}^2,$$

$$e(k)=\frac{1}{4Q{k}_{\ast }^2}\{\begin{array}{ll}(k/k_{\ast }{)}^2& \text{for}0<k<k_{\ast }\\ (k_{\ast }/k{)}^2& \text{for}k>k_{\ast }\end{array}.$$

$$e_{\mathrm{h}}(k_{\mathrm{h}})=\iiint \delta (k\sin \theta -k_{\mathrm{h}})\frac{e(k)\delta (\theta -{\theta }_{\ast })}{2\pi k^2\sin \theta }\mathrm{d}\mathbf{k}=\csc {\theta }_{\ast }e(k_{\mathrm{h}}\csc {\theta }_{\ast }),$$

$$k_{\mathrm{h}}\sim \frac{\omega }{⟨|\mathbf{U}{|}^2{⟩}^{1/2}}\sim \frac{K_{\mathrm{h}\ast }}{\mathrm{Ro}}\text{and}{k}_{\mathrm{h}}\sim \frac{K_{\mathrm{h}\ast }}{\mathrm{Ro}}\frac{⟨|\mathbf{U}{|}^2{⟩}^{1/2}}{⟨|\mathbf{u}{|}^2{⟩}^{1/2}},$$

$$a=a^{(0)}(\mathbf{X},\mathbf{k},T)+\epsilon a^{(1)}(\mathbf{x},\mathbf{X},\mathbf{k},t,T)+\cdots$$

$${\partial }_t{a}^{(1)}+c_i{\partial }_{x_i}{a}^{(1)}=k_m{\partial }_{x_i}{U}_m{\partial }_{k_i}{a}^{(0)},$$

$$a^{(1)}(\mathbf{x},\mathbf{X},\mathbf{k},t,T)=k_m{\int }_0^t{\partial }_{x_j}{U}_m(\mathbf{x}-\mathbf{c}s,T)\mathrm{d}s{\partial }_{k_j}{a}^{(0)}.$$

$${\partial }_T{a}^{(0)}+c_i{\partial }_{X_i}{a}^{(0)}=k_n⟨{\partial }_{x_i}{U}_n{\partial }_{k_i}{a}^{(1)}⟩,$$

$${\partial }_T{a}^{(0)}+c_i{\partial }_{X_i}{a}^{(0)}={\partial }_{k_i}\left(D_{ij}{\partial }_{k_j}{a}^{(0)}\right)$$

$$D_{ij}=-k_m{k}_n{\int }_0^{\infty }⟨{\partial }_{x_i}{U}_n(\mathbf{x}){\partial }_{x_j}{U}_m(\mathbf{x}-\mathbf{c}s)⟩\mathrm{d}s.$$

$$D_{ij}=\frac{k_m{k}_n}{8{\pi }^2}{\int }_{{\mathbb{R}}^3}{K}_i{K}_j{\hat{\Pi }}_{mn}(\mathbf{K})\delta (\mathbf{K}\cdot \mathbf{c})\mathrm{d}\mathbf{K},$$

$$\mathbf{c}(\mathbf{k})=\frac{(N^2-f^2)\cos \theta \sin \theta }{k\omega }{\mathbf{e}}_{\theta }.$$

$$\text{mathsfbi}D=D_{kk}{\mathbf{e}}_k\otimes {\mathbf{e}}_k+D_{k\varphi }({\mathbf{e}}_k\otimes {\mathbf{e}}_{\varphi }+{\mathbf{e}}_k\otimes {\mathbf{e}}_{\varphi })+D_{\varphi \varphi }{\mathbf{e}}_{\varphi }\otimes {\mathbf{e}}_{\varphi },$$

$$\mathbf{K}$$

$$+K\sin \Theta \sin \gamma {\mathbf{e}}_{\varphi }+K(\sin \Theta \cos \theta \cos \gamma -\cos \Theta \sin \theta ){\mathbf{e}}_{\theta },$$

$$\delta (\mathbf{K}\cdot \mathbf{c})=\frac{k\omega \left(\delta (\gamma -{\gamma }_{\ast })+\delta (\gamma +{\gamma }_{\ast })\right)}{K(N^2-f^2)\sin \Theta \sin \theta {\cos }^2\theta \sin {\gamma }_{\ast }},$$

$$k_m{k}_n{\hat{\Pi }}_{mn}=(k_1{K}_2-k_2{K}_1{)}^2⟨\hat{\psi }(\mathbf{K})\hat{\psi }(-\mathbf{K})⟩=2k^2{\sin }^2\theta {\sin }^2\gamma E(\mathbf{K}),$$

$$D_{kk}$$

$$D_{\varphi \varphi }$$

$$k^3{f}^{\prime \prime }+k^2{f}^{\prime }+4\left({\lambda }^2/16-k\right)f=0,$$

$$e(k,t)={\int }_0^{\infty }A(\lambda )J_4(\lambda /\sqrt{k}){\mathrm{e}}^{-Q{\lambda }^{2}t/4}\mathrm{d}\lambda ,$$

Full text

1 Introduction

The dynamics of rotating stratified fluids, most notably the atmosphere and ocean, is characterised by the coexistence of vortical flow and inertia-gravity 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 vortical 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 low-Rossby-number, 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 fixed-frequency IGWs. The diffusivity components associated with radial and angular diffusion on the cone are obtained in closed forms involving the IGW 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 initial-value problem (§3) and a steady forced problem (§4), assuming horizontally isotropic IGW fields, and we test the results against numerical simulations of the three-dimensional Boussinesq equations, finding good agreement in both cases. With forcing, the diffusion equation predicts a constant-flux, steady energy spectrum scaling with wavenumber as $k^{-2}$ 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 $k^{-2}$ 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 initial-value 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 vortical flow of much larger spatial scale so that the WKB approximation applies. The distribution of wave energy in the $(\boldsymbol{x},\boldsymbol{k})$ phase space is then governed by the conservation

[TABLE]

of the action density $a(\boldsymbol{x},\boldsymbol{k},t)$. Here $\Omega=\omega+\boldsymbol{U}\cdot\boldsymbol{k}$ is the frequency, which sums the intrinsic frequency

[TABLE]

where $f<N$ are the Coriolis and buoyancy frequencies and $\theta$ is the angle between the wavevector $\boldsymbol{k}$ and the vertical, and the Doppler shift $\boldsymbol{U}\cdot\boldsymbol{k}$, where $\boldsymbol{U}=\boldsymbol{U}(\boldsymbol{x},t)$ is the vortical flow velocity. Assuming that the flow is (i) weak enough that $\omega\gg\boldsymbol{U}\cdot\boldsymbol{k}$, (ii) evolving on a time scale much longer than $\omega^{-1}$, and (iii) well modelled by a homogeneous and stationary random field, we can approximate (1) by

[TABLE]

where $\boldsymbol{c}=\nabla_{\boldsymbol{k}}\omega$ is the intrinsic group velocity and $\mathsfbi{D}$ a $\boldsymbol{k}$-dependent diffusivity tensor (see Appendix A for a derivation). The right-hand side of (3) captures the scattering of wave action that results from small-but-sustained random Doppler shifting by the flow; in the regime considered, this naturally leads to diffusion in $\boldsymbol{k}$-space. In Cartesian coordinates, the diffusivity tensor takes the form

[TABLE]

where $\Pi_{mn}(\boldsymbol{x})=\langle U_{m}(\boldsymbol{y}+\boldsymbol{x})U_{n}(\boldsymbol{y})\rangle$ is the velocity correlation tensor, with $\langle\cdot\rangle$ 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 induced-diffusion regime (see Müller et al., 1986, §5, for a review). Müller & Olbers (1975) and Müller (1976, 1977) discussed a flow-induced diffusivity that differs from (4) to account heuristically for wave–wave interactions and dissipation.

A key property of (4) is that $\mathsfbi{D}(\boldsymbol{k})\cdot\boldsymbol{c}(\boldsymbol{k})=0$ since

[TABLE]

Thus there is no diffusion in the direction of the group velocity $\boldsymbol{c}$. Since $\boldsymbol{c}$ is perpendicular to constant-frequency surfaces, for the IGW dispersion relation (2) diffusion is restricted to the cones $\theta=\mathrm{const}$, see Fig. 1. This is because diffusion in $\boldsymbol{k}$-space stems from resonant-triad interactions between two IGWs and one vortical mode (also termed balanced mode) associated with the flow. The flow is treated as a zero-frequency mode because it evolves slowly compared with $\omega^{-1}$, 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 (4) to IGWs and geostrophic flows using the dispersion relation (2) and the geostrophic balance satisfied by the velocity in $\Pi_{mn}$. It is natural to use spherical polar coordinates $(k,\phi,\theta)$ in $\boldsymbol{k}$-space and a Fourier counterpart to $\Pi_{mn}$ in the form of the vortical flow kinetic energy spectrum $E(K_{\mathrm{h}},K_{\mathrm{v}})$, which we assume to be horizontally isotropic so that it only depends on the horizontal and vertical wavenumbers $K_{\mathrm{h}}$ and $K_{\mathrm{v}}$ (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 (3) to

[TABLE]

under the further assumption of spatial homogeneity $\nabla_{\boldsymbol{x}}a=0$. This makes it plain that there is no diffusion in the direction of $\theta$. Hence, $\theta$, or equivalently $\omega$, can be treated as a fixed parameter. The only non-zero components of the diffusivity tensor are given by

[TABLE]

where

[TABLE]

depends solely on $\theta$, $N$ and $f$.

Eqs. (6)–(7) provide a full description of the diffusion of IGW on the constant-frequency cone in $\boldsymbol{k}$-space for a turbulent flow of given energy spectrum. In the angular $\phi$-direction, this diffusion leads to an isotropisation of the wave field with rate $D_{\phi\phi}/(k^{2}\sin^{2}\theta)$. In the radial $k$-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 downward-propagating IGWs. This is only an approximation; exchanges between upward- and downward-propagating 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 $\phi$, $\partial_{\phi}a=0$, leaving the study of horizontal isotropisation for future work.

3 Initial-value problem

For $a=a(k,t)$, we rewrite (6) as

[TABLE]

where we have introduced $e(k,t)=2\pi k^{2}\,\sin\theta\,\omega a(k,t)$ and the $k$-independent parameter $Q=D_{kk}/k^{3}$. The function $e(k,t)$ is the IGW energy density in $k$, with $e(k,t)\,\mathrm{d}k$ the energy contained within the interval $[k,k+\mathrm{d}k]$. We solve (9) with initial condition $e(k,0)=\delta(k-k_{*})$ corresponding to the excitation of IGWs with a single wavenumber $k_{*}$. (The solution associated with arbitrary initial condition can be deduced by integration over $k_{*}$.) We show in Appendix A.3 that

[TABLE]

where $J_{4}$ is a Bessel function of the first kind (DLMF, 2018). The large-time behaviour of $e(k,t)$ is readily deduced as $e(k,t)\propto k^{-2}t^{-5}$ away from an asymptotically small neighbourhood of $k=0$ (see Appendix A.3). An inverse diffusion time scale can be read off from (10) as $Qk_{*}$. Using (7) this can be written in the dimensionless form

[TABLE]

where $\gamma$ is a dimensionless ‘geometric’ factor that depends only on $\theta$ and the shape (but not the magnitude) of the flow kinetic-energy spectrum and $\mathrm{Ro}=K_{\mathrm{h}*}\langle|\boldsymbol{U}|^{2}\rangle^{1/2}/f$ is a flow Rossby number. The typical horizontal and vertical inverse flow scales $K_{\mathrm{h}*}$ and $K_{\mathrm{v*}}$ are assumed to be related by $K_{\mathrm{v}*}=NK_{\mathrm{h}*}/f$. Eq. (11) 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 $Qk_{*}/\omega\ll 1$ in addition to the WKB conditions $k_{*}\sin\theta\gg K_{\mathrm{h}}$ and $k_{*}\cos\theta\gg K_{\mathrm{v}}$.

We verify the solution of (9) against simulations of the three-dimensional non-hydrostatic Boussinesq equations. These are solved using a code adapted from that in Waite & Bartello (2006b) which relies on a de-aliased pseudospectral method and a third-order Adams–Bashforth scheme with timestep $0.015/f$. The triply-periodic domain, $(2\pi)^{3}$ in the scale coordinates $(x,y,z^{\prime}=Nz/f)$, is discretised uniformly with $768^{3}$ grid points. A hyperdissipation of the form $-\nu(\partial_{x}^{2}+\partial_{y}^{2}+\partial_{z^{\prime}}^{2})^{4}$, with $\nu=2\times 10^{-17}$, is employed in the momentum and density equations. We take $N/f=32$, a representative value of mid-depth 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 vortical flow peaks at $K_{\mathrm{h*}}\simeq 4$ and has an inertial subrange scaling approximately as $K_{\mathrm{h}}^{-3}$ and $K_{\mathrm{v}}^{-3}$. This spectrum evolves slowly over the IGW-diffusion timescale, and an average is used to calculate $D_{kk}$ in (7a), and hence $Q$ in (9). We report experiments with the two Rossby numbers $\mathrm{Ro}=K_{\mathrm{h}*}\langle|\boldsymbol{U}|^{2}\rangle^{1/2}/f=0.057,\,0.117$ (or $\langle\zeta^{2}\rangle^{1/2}/f=0.1,\,0.2$ for the alternative Rossby numbers based on the vertical vorticity $\zeta$), and the two IGW frequencies $\omega=2f,\,3f$. Upward-propagating IGWs are initialised as a ring in $\boldsymbol{k}$-space with $k_{\mathrm{h*}}=16$, $k_{\mathrm{v}}=\mathrm{\cot}\,\theta\,k_{\mathrm{h}}$, random phases, and an initial kinetic energy $\langle|\boldsymbol{u}|^{2}\rangle/2=0.1\langle|\boldsymbol{U}|^{2}\rangle/2$. The IGW spectrum $e(k,t)$ is computed following the normal-mode decomposition of Bartello (1995).

Fig. 1, obtained for the lower $\mathrm{Ro}$ and $\omega=3f$, illustrates the confinement of wave energy on the constant-frequency 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 $2\omega$ and $3\omega$. Fig. 2 shows the evolution of $e(k,t)$ for the four sets of values of $(\mathrm{Ro},\omega)$. The numerical results are compared with the predictions of the diffusion equation obtained by solving (9) initialised with the form of $e(k,t_{\mathrm{a}})$ extracted from the simulation after an adjustement time $t_{\mathrm{a}}>0$. This procedure accounts for the fact that the diffusion equation (3) is only valid after an adjustment period, requiring $t_{\mathrm{a}}\gg(K_{*}|\boldsymbol{c}|)^{-1}$, 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 $t$ is scaled suitably. The decay is slower for $\omega=3f$ than $\omega=2f$, consistent with a decrease in $Q$ obtained when evaluating (7a). Scattering from upward-propagating to downward-propagating IGWs, neglected in the diffusion approximation, occurs; it is more substantial for the larger $\omega$ because the two nappes of the constant-frequency cones are closer together, facilitating energy transfers. Fig. 3 displays the wave energy spectrum $e(k,t)$ obtained for $\omega=2f$ (top row of Fig. 2) in log–log coordinates. It shows that the good agreement between numerical and predicted spectra extends to large wavenumbers for $\mathrm{Ro}=0.057$ but not for the larger Rossby, $\mathrm{Ro}=0.117$, at the later times. We note that the wave energy is then very small and may be affected by a contribution associated with spontaneous generation (cf. Kafiabad & Bartello, 2018)

4 Forced response and observed ocean and atmosphere spectra

We now turn to the steady solution of (9) in the presence of a forcing of the form $\delta(k-k_{*})$. Eq. (10) admits two steady solutions: the no-flux solution $e(k)\propto k^{2}$ and the constant-flux solution $e(k)\propto k^{-2}$. Matching these at $k_{*}$ yields the steady spectrum

[TABLE]

Note that for IGWs with a single frequency and correspondingly a single angle $\theta_{*}$, the horizontal energy spectrum $e_{\mathrm{h}}(k_{\mathrm{h}})$ satisfies the same power laws as $e(k)$ since

[TABLE]

using that the energy density in $\boldsymbol{k}$-space is $e(k)\delta(\theta-\theta_{*})/(2\pi k^{2}\sin\theta)$. Thus, (12) implies a $k_{\mathrm{h}}^{-2}$ horizontal spectrum at large $k_{\mathrm{h}}$. This remains true for a superposition of IGWs with different frequencies, corresponding to an integration over $\theta_{*}$.

We confirm the prediction (12) 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 $\omega=2f$ (see Waite, 2017). The forcing amplitude is adjusted so that the wave energy is about $0.01$ of the vortical flow energy after reaching a stationary state. Fig. 4 shows the stationary spectra for $\mathrm{Ro}=0.057$ and $0.117$. For the low Rossby number, the prediction (12) is well borne out by the simulation results with a clean $k_{\mathrm{h}}^{-2}$ spectrum spanning nearly a decade from the forcing scale down to dissipation. For the high Rossby number, the spectrum shallows a little from wavenumber $40$ or so to take a shape more consistent with $k_{\mathrm{h}}^{-5/3}$. Two mechanisms can be invoked to explain this shallowing: the Doppler term is not weak compared with the intrinsic IGW frequency, invalidating the diffusion approximation, or nonlinear wave–wave interactions become significant. We can roughly estimate the wavenumbers at which each of these mechanisms is important as

[TABLE]

corresponding to order-one Rossby numbers based on the wave lengthscale $k_{\mathrm{h}}^{-1}$ and on the root-mean-square velocity of, respectively, the vortical flow and the IGWs. For the simulation with $\mathrm{Ro}=0.117$, these wavenumbers are about 40 and 400, suggesting that the shallowing of the spectrum is associated with the breakdown of the assumption of weak Doppler shift.

The prediction of a $k_{\mathrm{h}}^{-2}$ spectrum is significant in view of the ubiquity of this scaling in ocean and atmosphere observations. In the ocean, kinetic energy spectra show a $k_{\mathrm{h}}^{-2}$ 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 $k_{\mathrm{h}}^{-2}$ 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 $k_{\mathrm{h}}^{-2}$ dependence. While this spectrum is generally attributed to wave–wave interactions (e.g. Müller et al., 1986; Lvov et al., 2012), interactions with the geostrophic flow may play a significant role, dominating for wavenumbers much smaller than (14b). We emphasise that theories based on linear IGWs, be it the diffusion approximation of this paper or a more general theory accounting for strong Doppler shift, cannot predict the level of IGW spectrum nor its frequency content since both are determined by the forcing.

In the atmosphere, similarly, there is a broad range of scales, from 500 km to 10 km, where the energy spectrum scales approximately as $k_{\mathrm{h}}^{-2}$. This is the shallow, mesoscale part of the celebrated Nastrom & Gage (1985) spectrum, which is traditionally interpreted as a $k_{\mathrm{h}}^{-5/3}$ spectrum but is also consistent with $k_{\mathrm{h}}^{-2}$. 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 power-law spectra …  but an explanation for the spectral shape is so far missing’. Our results provide a possible explanation.

The total spectrum in the high-$\mathrm{Ro}$ simulation, shown in Fig. 5, is reminiscent of atmospheric observations, with a $k_{\mathrm{h}}^{-3}$ range at large scales associated with the vortical flow, a $k_{\mathrm{h}}^{-2}$ range associated with nearly linear IGWs at intermediate scales, and a further shallowing at small scales (best seen in Fig. 4, right panel). While the diffusion approximation explains the $k_{\mathrm{h}}^{-2}$ range in our simulations, a degree of caution is required to draw a similar conclusion for the atmospheric spectrum since some of the underlying assumptions – weak flow with homogeneous statistics and relevance of the equilibrium spectrum in particular – are questionable.

5 Discussion

This paper examines the impact of a turbulent geostrophic flow on the statistics of small-amplitude IGWs. This impact has received less attention than that paid to wave–wave interactions. Yet the timescale found for a substantial effect of the geostrophic flow, of the order of $0.1\,\mathrm{Ro}^{-2}f^{-1}$ (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 $29^{\circ}$, 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).

The present paper focuses on the diffusive regime of IGW scattering that arises for weak flows and small-scale, linear IGWs. A remarkable feature of this regime is the prediction of a $k^{-2}_{\mathrm{h}}$ 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 downward-propagating 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, as required for wavenumbers not small compared with (14a), 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. We conclude by noting that the consistency between predicted and observed spectral slopes is only indicative: further investigations are needed to establish the importance of IGW scattering in determining oceanic and atmospheric spectra.

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, 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. This work used the ARCHER UK National Supercomputing Service.

Appendix A Derivation of the diffusion equation and of its solution

A.1 General wave systems

We introduce a small parameter $\varepsilon\ll 1$ in the action conservation (1) by writing the frequency as $\Omega=\omega+\varepsilon\boldsymbol{U}\cdot\boldsymbol{k}$, 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 $T=\varepsilon^{2}t$ and $\boldsymbol{X}=\varepsilon^{2}\boldsymbol{x}$, we substitute the expansion

[TABLE]

into (1). The first non-trivial equation appears at $O(\varepsilon)$ and is given by

[TABLE]

using Cartesian components and implied summation. Assuming that the velocity field varies on the slow time scale only, the solution is given by

[TABLE]

Averaging the next-order equation to eliminate the terms containing $a^{(2)}$, we find

[TABLE]

since $\langle U_{i}\partial_{x_{i}}a^{(1)}\rangle=\langle\partial_{x_{i}}(U_{i}a^{(1)})\rangle=0$ using incompressibility and spatial homogeneity. Substituting the limit of (17) as $t\to\infty$ as appropriate for the slow dynamics, we obtain the diffusion equation

[TABLE]

with the diffusivity

[TABLE]

This can be written as (4) in terms of the correlation tensor $\Pi_{mn}$ or, alternatively, as

[TABLE]

in terms of the Fourier transform $\hat{\Pi}_{mn}$ of $\Pi_{mn}$. The diffusive approximation (19) is standard for Hamiltonian systems with weak random perturbation and has been obtained in a variety of contexts (e.g., McComas & Bretherton (1977) for wave–wave interactions). The formal derivation above follows Bal et al. (2010, §4.2) who also discuss its rigorous justification.

A.2 IGWs in quasigeostrophic flow

We particularise (21) to the IGW dispersion relation (2) and a velocity field of the form $\boldsymbol{U}=(-\partial_{x_{2}}\psi,\partial_{x_{1}}\psi,0)$ with $\psi$ the geostrophic streamfunction. We use the spherical polar coordinates $(k,\theta,\phi)$ for $\boldsymbol{k}$, with $\boldsymbol{e}_{k}$, $\boldsymbol{e}_{\theta}$ and $\boldsymbol{e}_{\phi}$ the corresponding unit vectors, and express the group velocity as

[TABLE]

The diffusivity can be written in the basis $(\boldsymbol{e}_{k},\boldsymbol{e}_{\theta},\boldsymbol{e}_{\phi})$ as

[TABLE]

where $\mathsfit{D}_{kk}=\boldsymbol{e}_{k}\cdot\mathsfbi{D}\cdot\boldsymbol{e}_{k}$, $\mathsfit{D}_{k\phi}=\boldsymbol{e}_{k}\cdot\mathsfbi{D}\cdot\boldsymbol{e}_{\phi}$, $\mathsfit{D}_{\phi\phi}=\boldsymbol{e}_{\phi}\cdot\mathsfbi{D}\cdot\boldsymbol{e}_{\phi}$, and we have made use of the fact that $\mathsfbi{D}\cdot\boldsymbol{e}_{\theta}\propto\mathsfbi{D}\cdot\boldsymbol{c}=0$ to eliminate all components along $\boldsymbol{e}_{\theta}$.

With $\Theta$ and $\Phi$ the polar and azimuthal angles of the flow wavevector $\boldsymbol{K}$, we have

[TABLE]

where $\gamma=\Phi-\phi$. Hence the delta function in (21) can be written as

[TABLE]

where $0\leq\gamma_{*}=\cos^{-1}(\tan\theta/\tan\Theta)\leq\pi$. We also note that

[TABLE]

where $E(\boldsymbol{K})=K^{2}\sin^{2}\Theta\langle\hat{\psi}(\boldsymbol{K})\hat{\psi}(-\boldsymbol{K})\rangle/2$ is the flow kinetic energy spectrum. We now introduce (25)–(26) into (21) projected onto $\boldsymbol{e}_{k}$ and $\boldsymbol{e}_{\phi}$ to compute the components of $\mathsfbi{D}$ in (23). Assuming that the flow is isotropic in the horizontal so that $E(\boldsymbol{K})$ is independent of $\gamma$, we obtain after some simplifications

[TABLE]

and $D_{k\phi}=0$. The form (7) follows by replacing the kinetic-energy spectrum $E(\boldsymbol{K})$ by its two-dimensional counterpart $E(K_{\mathrm{h}},K_{\mathrm{v}})=2\pi K_{\mathrm{h}}E(\boldsymbol{K})$ and changing the integration variables from $(K,\Theta)$ to $(K_{\mathrm{h}},K_{\mathrm{v}})$, with $K\,\mathrm{d}K\mathrm{d}\Theta=\mathrm{d}K_{\mathrm{h}}\mathrm{d}K_{\mathrm{v}}$.

A.3 Solution of (9) and its long-time approximation

Introducing a solution of the separable form $e(k,t)=\mathrm{e}^{-Q\lambda^{2}t/4}f(k,\lambda)$, with $\lambda\geq 0$ a spectral parameter, into (9) leads to

[TABLE]

where the prime denotes derivative with respect to $k$. Solutions bounded as $k\to 0$ are proportional to the Bessel function $J_{4}(\lambda/\sqrt{k})$. The general solution of (9) follows as

[TABLE]

for an arbitrary function $A(\lambda)$. Imposing the initial condition $e(k,0)=\delta(k-k_{*})$ yields (10) on using the Bessel-function expansion of $\delta(k-k_{*})$ (DLMF, 2018, Eq. 1.17.13).

For large $t$, the integral in (10) is dominated by a neighbourhood of $\lambda=0$. The Bessel functions $J_{4}$ can therefore by replaced by their small-argument approximation, $J_{4}(z)\sim z^{4}/16$ as $z\to 0$ (DLMF, 2018, Eq. 10.2.1), leading to

[TABLE]