Inhomogeneity-induced wavenumber diffusion

Abstract

Inertia–gravity waves are scattered by background flows as a result of Doppler shift by a non-uniform velocity. In the Wentzel–Kramers–Brillouin regime, the scattering process reduces to a diffusion in spectral space. Other inhomogeneities that the waves encounter, such as density variations, also cause scattering and spectral diffusion. We generalise the spectral diffusion equation to account for these inhomogeneities. We apply the result to a rotating shallow-water system, for which height inhomogeneities arise from velocity inhomogeneities through geostrophy, and to the Boussinesq system for which buoyancy inhomogeneities arise similarly. We compare the contributions that height and buoyancy variations make to the spectral diffusion with the contribution of the Doppler shift. In both systems, we find regimes where all contributions are significant. We support our findings with exact solutions of the diffusion equation and with ray tracing simulations in the shallow-water case.

1. Introduction

Inertia–gravity waves (IGWs) propagate in the atmosphere and ocean under the restoring forces of buoyancy and Coriolis effect. As they propagate, they encounter and interact with a variety of inhomogeneities, including background flows, topography and other waves. These inhomogeneities refract, reflect and, in the case of background flows, advect the waves. For example, internal tides – IGWs generated at the semi-diurnal and diurnal frequencies by astronomical forcing – are advected and refracted by background flows (e.g. Park & Watts 2006; Rainville & Pinkel 2006; Chavanne et al. 2010; Pan, Haley & Lermusiaux 2021) and reflected by bottom topography (e.g. Müller & Xu 1992; Buhler & Holmes-Cerfon 2011; Kelly et al. 2013; Lahaye, Gula & Roullet 2020; Pan et al. 2021). The IGWs are also affected by each other (e.g. Eden & Olbers 2014). In the limit of weak, repeated interactions, wave energy is redistributed in spectral space in a scattering process which can be described statistically.

Kinetic equations derived from the governing fluid equations provide this statistical description of scattering by weak, slowly evolving, random inhomogeneities (see Hasselmann (1966) and Ryzhik, Papanicolaou & Keller (1996) for general formulations). For example, Hasselmann (1966) and Ardhuin & Herbers (2002) derive a kinetic equation for surface waves scattered by gently sloping bottom topography; Eden, Pollmann & Olbers (2019) derive a scattering equation for IGWs interacting with a weak, random, slowly evolving wave field (see also Eden, Pollmann & Olbers 2020); and Danioux & Vanneste (2016), Savva & Vanneste (2018) and Savva, Kafiabad & Vanneste (2021) derive kinetic equations for near-inertial waves, internal tides and IGWs, respectively, scattered by weak, slowly evolving turbulence. A key feature of scattering is that if the scattering inhomogeneities evolve sufficiently slowly compared with the waves, wave frequency is preserved.

This paper concerns the Wentzel–Kramers–Brillouin (WKB) limit of large-scale inhomogeneities scattering small-scale waves, in which scattering reduces to diffusion in spectral space. The corresponding induced-diffusion equation can be derived by invoking conservation of wave action density, which holds in the WKB limit. This was first done by McComas & Bretherton (1977) in the context of wave–wave interactions using the ray equations – the characteristic equations for action conservation. (McComas & Bretherton (1977) note that their derivation also applies to diffusion induced by low-frequency currents. A flow-induced diffusivity is first discussed in Müller & Olbers (1975) and expanded on in Müller (1976, 1977). Recently, alternative derivations for the diffusion equation start with the conservation of wave action and use multi-scale asymptotics. They have been carried out for weak geostrophic flows scattering (i) IGWs in a three-dimensional (3-D) Boussinesq system (Kafiabad, Savva & Vanneste 2019, hereafter KSV 2019); (ii) Poincaré waves in a rotating shallow-water system (Dong, Buhler & Smith 2020); and (iii) deep-water surface waves (Villas Bôas & Young 2020). In all these derivations, the geostrophic flow impacts wave propagation solely through the Doppler shift term of the wave dispersion relation. McComas & Bretherton (1977) and Savva (2020) show that induced diffusion is the WKB limit of a scattering integral for wave–wave and wave–flow interactions, respectively. Yang et al. (2023) investigate the relevance of diffusion theories in a realistic ocean simulation.

The aims of this paper are twofold: (i) we argue at the level of the dispersion relation that inhomogeneities other than Doppler shift can be significant in scattering waves in the WKB regime (§ 2) and (ii) we derive the spectral diffusivity induced by any weak inhomogeneity (§ 3), thus generalising the original result of McComas & Bretherton (1977). For, say, waves scattered by bottom topography (e.g. Hasselmann 1966; Müller & Xu 1992; Ardhuin & Herbers 2002), it is widely appreciated that scattering mechanisms other than Doppler shift are significant. However, this point has at times been overlooked in the study of wave scattering by mean flows as we detail below.

In § 2 we use scaling arguments to compare the Doppler shift term of the dispersion relation with two ofttimes neglected inhomogeneities: height fluctuations in a rotating shallow-water system and buoyancy gradients in a Boussinesq system. The height fluctuation effects are neglected in Dong et al. (2020) and the buoyancy fluctuation term is neglected in previous work of the authors, KSV (2019) and Cox, Kafiabad & Vanneste (2023, hereafter CKV 2023). For both systems, we find regimes where these inhomogeneities can be significant. Doppler shift and vertical flow buoyancy gradients are accounted for in tidal ray tracing implemented by Chavanne et al. (2010). They find both to be important and, in their simulations, refraction by these gradients is more significant in the transfer of wave energy than advection by Doppler shift. (The ray tracing formulation used by Chavanne et al. (2010) is outlined in Olbers (1981) for IGWs propagating in geostrophic flows with vertical and horizontal shears.) Doppler shift, refraction through a background flow, buoyancy gradients and topography are all found to be significant in a coupled set of tidal equations derived by Pan et al. (2021).

Our analysis throughout this paper applies to weak scattering of waves by inhomogeneities. In practice, this means that the scattering inhomogeneities induce small perturbations to the wave frequencies. The wave amplitudes are assumed small enough that linear wave theory applies and the wave feedback on the background flow is negligible.

In § 3, we impose further statistical assumptions on the inhomogeneities, assume they are slowly time dependent and follow the perturbation expansion of KSV (2019) to reach a general diffusion equation for any inhomogeneity. We evaluate the general diffusivity for our two example systems (the method is similar for both systems and so the Boussinesq evaluation is relegated to Appendix B.1). We then revise the scaling arguments of the previous section. For the shallow-water system, we support our analysis with ray tracing simulations (§ 3.2), finding good agreement with the exact solution for two-dimensional (2-D) wave action given in Villas Bôas & Young (2020). For the Boussinesq system, we find the forced equilibrium spectrum of wave energy. We also evaluate the Boussinesq diffusivity for a typical quasi-geostrophic flow. We find the revised importance estimate of buoyancy fluctuation to Doppler shift effects obeys negative power laws in horizontal wavenumber and frequency as explained in Appendix B.2.

In § 4, we evaluate the general diffusivity for scattering by topography of (i) rotating shallow-water waves (§ 4.1) and (ii) surface waves propagating on a background current (§ 4.2). This is proof of concept that our general diffusivity has applications outside of waves scattered by mean flows.

Section 5 is a discussion of our results including possible applications and limitations. For completeness, we derive wave action conservation for a rotating shallow-water system with significant height fluctuation effects in Appendix A. Lists of key symbols are included in Appendix C.

2. Scaling arguments

In this section, we introduce a general inhomogeneity term into the wave dispersion relation for waves in the WKB regime. This term encompasses Doppler shift by a background flow but can include additional inhomogeneities. We motivate the inclusion of two such inhomogeneities: flow-induced height fluctuations in a rotating shallow-water system and vertical buoyancy gradients in a 3-D Boussinesq system. We find regimes where the additional inhomogeneities are significant, at points dominant. This has implications to ray tracing simulations which do not always take into account inhomogeneities other than Doppler shift and motivates the general diffusion equation derivation of § 3, because previous diffusion equations of the form first introduced by McComas & Bretherton (1977) have only considered Doppler shift (by a background flow or long waves) as the scattering mechanism.

Figure 1. Sketch of flat-bottom shallow-water set-up. In (a,b), the free surface is the solid blue line, the flat bottom is the solid black line and the dashed black line is the constant mean height $\bar {H}$ . The wave perturbations are given by $h$ . The height of the layer in the absence of waves is $H$ which is equal to $\bar {H}$ in (a), whilst in (b) it includes geostrophic height corrections $\Delta H$ , given by the dashed-dotted line.

Our starting point is the conservation of wave action $a$ in $(\boldsymbol {x}, \boldsymbol {k})$ space in the WKB limit of small-scale waves scattered by large-scale inhomogeneities:

(2.1) \begin{equation} \partial _{t} a + \boldsymbol {\nabla }_{\boldsymbol {k}} \omega \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}} a - \boldsymbol {\nabla }_{\boldsymbol {x}} \omega \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {k}} a = 0. \end{equation}

(This is derived for the rotating shallow-water system in Appendix A but holds more generally.) The absolute frequency of the waves $\omega$ is the sum of a homogeneous part, $\omega _{0}$ , and a weak, inhomogeneity-induced part, $\omega _{1}$ :

(2.2) \begin{equation} \omega (\boldsymbol {x}, \boldsymbol {k}, t) = \omega _{0}(\boldsymbol {k}) + \omega _{1}(\boldsymbol {x}, \boldsymbol {k}, t), \end{equation}

with $\omega _{1} \ll \omega _{0}$ . We call $\omega _{0}$ the bare frequency and assume that it varies sufficiently slowly in time and space as to be considered constant. The frequency correction $\omega _{1}$ includes the Doppler shift $\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}$ induced by a background velocity $\boldsymbol {U}$ , and any other inhomogeneities. We consider two systems where other inhomogeneities inevitably arise in the presence of Doppler shift.

2.1. Shallow water

We consider Poincaré waves propagating in rotating shallow water with flat bottom and a background geostrophic flow. The background velocity $\boldsymbol {U} = (U, V, 0)$ and height $H$ are related through the geostrophic balance

(2.3) \begin{equation} f \boldsymbol {e}_z \times \boldsymbol {U} = - g \boldsymbol {\nabla }_{\boldsymbol {x}} \Delta H, \end{equation}

where $f$ is the Coriolis frequency, $g$ the gravitational constant and $\boldsymbol {e}_z$ the unit vertical vector. Here we have introduced $\Delta H$ , the geostrophic perturbation to the mean height $\bar H$ such that $H = \bar H + \Delta H$ (see figure 1). The perturbation $\Delta H$ , and hence $\boldsymbol {U}$ , is assumed to vary slowly in time and space compared with the period and wavelength of the Poincaré waves. For completeness we verify that the action conservation (2.1) holds in this case in Appendix A.

Assuming that $\Delta H \ll \bar {H}$ , the frequency of waves with wavevector $\boldsymbol {k} = (k_1, k_2)$ can be approximated as

(2.4) \begin{equation} \omega = \big(f^2 + g (\bar {H} + \Delta H) k_h^2\big)^{1/2} + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k} \approx \omega _{0} + \frac {g \Delta H k_h^2}{2 \omega _{0}} + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}, \end{equation}

where $k_h = (k_1^2 + k_2^2)^{1/2}$ is the wavenumber and $\omega _{0} = (f^2 + g \bar {H} k_h^2)^{1/2}$ is the intrinsic frequency for constant height $\bar {H}$ . Thus, there are two contributions to the frequency inhomogeneity,

(2.5) \begin{equation} \omega _{1} = \underbrace {\frac {g \Delta H k_h^2}{2 \omega _{0}}}_{\text{height fluctuation}} + \underbrace {\vphantom {\frac {g \Delta H k_h^2}{2 \omega _{0}}} \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}}_{\text{Doppler shift}}, \end{equation}

and hence two contributions to the scattering.

Figure 2. Ratio $R_{{sw}}$ (solid lines) defined by (2.8), the estimated relative importance of height fluctuation and Doppler shift terms in the shallow-water system, plotted for different values of the Burger number $Bu$ against non-dimensionalised horizontal wavenumber $k_h/K_* \gt 1$ . Quasi-geostrophic flow is $Bu = 1$ . Burger numbers $Bu = 0.1,\ 0.25$ are associated with the planetary geostrophic regime. Dashed lines are the square of this ratio, $R^{\prime}_{{sw}}$ (3.29), which gives a more accurate relative importance ratio for the diffusion regime, as shown in § 3.1.

We now compare the relative size of the terms in (2.5) to argue that height fluctuation effects can be just as important as Doppler shift in ray tracing and induced diffusion. We introduce the characteristic background flow speed $U_*$ , characteristic horizontal wavenumber $K_*$ and the flow Burger number

(2.6) \begin{equation} Bu = \frac {g\bar {H}}{f^2} K_*^2. \end{equation}

Using the geostrophic balance (2.3) and expressing $\omega _{0}$ in terms of $Bu$ , the height fluctuation term scales like

(2.7) \begin{equation} \frac {g \Delta H k_h^2}{2 \omega _{0}} \sim \frac { U_* k_h^2 }{2 K_* (1 + Bu (k_h/K_*)^2)^{1/2}}. \end{equation}

Then, the ratio between the height fluctuation and Doppler shift terms is given by

(2.8) \begin{equation} R_{{sw}} = \frac {\text {height fluctuation}}{\text {Doppler shift}} \sim \frac {k_h/K_*}{2(1 + Bu \, (k_h/K_*)^2)^{1/2}}. \end{equation}

In figure 2, we plot the ratio $R_{{sw}}$ given by (2.8) against non-dimensionalised wavenumber $k_h/K_*$ for different realistic values of the Burger number: $Bu = O(1)$ corresponds to the quasi-geostrophic regime and $Bu \ll 1$ to the planetary geostrophic regime. We take the limit of (2.8) for large $k_h/K_*$ , a necessary condition for (2.1) to hold:

(2.9) \begin{equation} R_{{sw}} \rightarrow Bu^{-1/2}/2, \quad k_h/K_* \gg 1. \end{equation}

We see that for $Bu = 0.25$ , $R_{{sw}} \rightarrow 1$ and the height fluctuation and Doppler shift terms in (2.5) have the same magnitude. For $Bu = 1$ corresponding to quasi-geostrophic flow, $R_{{sw}} \rightarrow 0.5$ . For smaller $Bu$ associated with planetary geostrophy, height fluctuations dominate over Doppler shift. For larger $Bu$ , Doppler shift dominates over height fluctuations.

In § 3.1, we show that a better estimate for the relative importance of height fluctuation to Doppler shift effects in the diffusion regime is $R^{\prime}_{{sw}} = R_{{sw}}^2$ , which is also plotted in figure 2.

2.2. Three-dimensional Boussinesq

Chavanne et al. (2010) compare the effect of refraction through flow buoyancy gradients with that of Doppler shift on internal tide propagation. Our set-up is similar to theirs – we consider 3-D IGWs of wavevector $\boldsymbol {k} = (k_1, k_2, k_3)$ propagating in a geostrophic flow $\boldsymbol {U} = (U, V, 0)$ with buoyancy $B$ . Unlike the tides, our waves are not in hydrostatic balance.

We first consider waves propagating with no background flow buoyancy gradients. The dispersion relation is

(2.10) \begin{equation} \omega = \left (f^2 \cos ^2 \theta + N^2 \sin ^2 \theta \right )^{1/2} + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}. \end{equation}

The intrinsic frequency is the first term, dependent on $\theta$ , the angle between $\boldsymbol {k}$ and the vertical, and $N$ is the buoyancy frequency. With uniform vertical buoyancy gradients, $N^2 = \bar {N}^2 = \text{const.}$ , the bare frequency coincides with the intrinsic frequency. To obtain (2.10), a WKB ansatz is substituted into the 3-D Boussinesq equations. We omit the full derivation, but it follows the same method as Appendix A for the shallow-water system. (See also derivations of action conservation with vertical buoyancy gradients by (i) Bretherton (1966) for gravity waves in a shear flow without the effect of rotation and (ii) Pan et al. (2021) for internal tides in hydrostatic balance.)

We make a rough argument for the inclusion of inhomogeneous vertical buoyancy gradients associated with the geostrophic flow in (2.10). Horizontal geostrophic balance and vertical hydrostatic balance lead to the thermal wind balance

(2.11) \begin{equation} \boldsymbol {f} \times \partial _{z}\boldsymbol {U} = \boldsymbol {\nabla }_{\boldsymbol {x},h} B. \end{equation}

Subscript $h$ indicates a purely horizontal gradient. This means that horizontal flow buoyancy gradients are induced by vertical shear. If the flow’s vertical shear is nonlinear, vertical buoyancy gradients are also induced.

In deriving (2.10), the buoyancy frequency $N^2$ only appears in the equation for wave buoyancy $b$ which is given, for negligible background flow buoyancy gradients, by

(2.12) \begin{equation} \partial _{t} b + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}} b + N^2 w = 0. \end{equation}

Here, $w$ is the vertical wave velocity and $N^2 = \bar {N}^2$ . Gradients in $\boldsymbol {U}$ have been neglected under the WKB ansatz (see e.g. Olbers 1981).

The full wave velocity is $\boldsymbol {u} = (u, v, w)$ . Introducing flow buoyancy gradients such that $B = B(\boldsymbol {x})$ , (2.12) becomes

(2.13) \begin{equation} \partial _{t} b + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}} b + \boldsymbol {u}_h \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x},h} B + N^2 w = 0, \end{equation}

where

(2.14) \begin{equation} N^2 = \bar {N}^2 + \partial _{z} B \end{equation}

and $\boldsymbol {u}_h = (u,v,0)$ is the horizontal wave velocity. This is (A5) of Appendix A of Chavanne et al. (2010). The vertical buoyancy gradient $\partial _{z}B$ acts as a perturbation to $\bar {N}^2$ . If horizontal gradients are neglected, then (2.10) remains the same with variable $N^2$ defined by (2.14).

Horizontal buoyancy gradients introduce $u$ and $v$ terms into the buoyancy equation. This changes the nature of the eigenvalue problem for $\omega$ . It can be shown, as discussed in Chavanne et al. (2010), that the perturbations to the total frequency (2.10) caused by these $u$ and $v$ terms are imaginary. This corresponds to growing and decaying wave modes which exchange energy with the background flow. This shear instability is because buoyancy is a non-conservative force in a baroclinic fluid (Jones 2001). Following Bretherton (1966), Olbers (1981) and Chavanne et al. (2010), we neglect these horizontal gradients which amounts to an assumption of large Richardson number and leave exploration of this instability and its interaction with wavevector diffusion to future work.

Accounting for weak flow-induced vertical buoyancy gradients, we expand the dispersion relation (2.10) with $N^2$ given by (2.14) to obtain

(2.15) \begin{equation} \omega = \underbrace {\vphantom {\frac {\partial _{z}B \sin ^2 \theta }{2\omega _0}}\left (f^2 \cos ^2 \theta + \bar {N}^2 \sin ^2 \theta \right )^{1/2}}_{\omega _{0}} + \underbrace {\frac {\partial _{z}B \sin ^2 \theta }{2\omega _0} + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}}_{\omega _{1}}, \end{equation}

valid for small $\partial _{z}B/\bar {N}^2$ . We compare the size of the Doppler shift and buoyancy gradient terms in $\omega _{1}$ . We introduce the background flow aspect ratio $\alpha$ defined by

(2.16) \begin{equation} \alpha = \frac {L_*}{H_*} = \frac {K_{v*}}{K_*}, \end{equation}

for characteristic length scales $L_* = 1/K_*$ and $H_* = 1/K_{v*}$ in the horizontal and vertical, respectively. Using the thermal wind relation (2.11) and aspect ratio, the buoyancy fluctuation term scales like

(2.17) \begin{equation} \frac {\partial _{z} B \sin ^2 \theta }{2 \omega _{0}} \sim \frac {\alpha ^2 U_* K_* \sin ^2 \theta }{2 (\cos ^2 \theta + (\bar {N}^2/f^2)\sin ^2 \theta )^{1/2}}. \end{equation}

Thus, the ratio $R_{{B}}$ between the buoyancy fluctuation and Doppler shift terms is roughly

(2.18) \begin{equation} {R_{{B}}} = \frac {\text{buoyancy fluctuation}}{\text{Doppler shift}} \sim \frac {\alpha ^2 \sin \theta }{2(\cos ^2 \theta +(\bar {N}^2/f^2) \sin ^2 \theta )^{1/2}(k/K_*)}, \end{equation}

where $k = k_h/\sin \theta$ is the wavenumber magnitude. It is instructive to consider $R_{{B}}$ as a function of non-dimensionalised frequency and horizontal wavenumber, $\omega _0/f$ and $k_h/K_*$ :

(2.19) \begin{equation} {R_{{B}}} \sim \frac {1}{2} \, \frac {\alpha ^2}{\bar {N}^2/f^2 - 1} \, \frac {(\omega _{0}/f)^2 - 1}{\omega _{0}/f} \, \frac {1}{k_h/K_*} \approx \frac {1}{2} \, \left (\frac {\alpha }{\bar {N}/f}\right )^2 \, \frac {(\omega _{0}/f)^2 - 1}{\omega _{0}/f} \, \frac {1}{k_h/K_*}. \end{equation}

The second expression holds for $(\bar {N}/f)^2 \gg 1$ , which is true in both the atmosphere and ocean. At $\omega _0 = f$ , ${R_{{B}}} = 0$ . This is to be expected: buoyancy effects are absent for vertically propagating inertial waves. The ratio attains a maximum value of ${R_{{B}}} \sim \alpha ^2(2(\bar {N}/f)(k_h/K_*))^{-1}$ when $\omega _{0} = \bar {N}$ . The ratio decays as $(k_h/K_*)^{-1}$ . In the WKB regime we consider, $k_h/K_* \gg 1$ and so it appears justified to assume the Doppler shift term dominates. However, for large values of $\alpha$ or when considering the lower limit of the WKB regime, this may not be the case.

Figure 3 shows the ratio (2.19) for $(\bar {N}/f)^2 \gg 1$ and three values of $\alpha /(\bar {N}/f)$ , including $\alpha \sim \bar {N}/f$ , a realistic regime for geostrophic turbulence. Contours of ${R_{{B}}} = 0.1,\ 1$ and $10$ , corresponding to a negligible, balanced and dominant buoyancy term, respectively, are given for each value of $\alpha$ . The buoyancy term can be equal to or greater than the Doppler shift term, namely for high aspect ratio $\alpha$ , higher frequencies and lower wavenumbers.

Figure 3. Ratio $R_{{B}}$ as given in (2.19), the relative importance of buoyancy fluctuation and Doppler shift terms in the full Boussinesq system, against non-dimensionalised frequency $\omega _0/f$ and horizontal wavenumber $k_h/K_*$ : (a) $\alpha = \bar {N}/(2f)$ , (b) $\alpha=\bar{N}/f$ and (c) $ \alpha = 3\bar {N}/f$ , with $(\bar {N}/f)^2 \gg 1$ . Contours are shown for ${R_{{B}}} = 0.1$ (dotted black), $1$ (solid black) and $10$ (dotted white).

3. Diffusion regime

In this section, we impose statistical assumptions on slowly evolving, weak inhomogeneities to derive a spectral diffusion equation from the action conservation equation (2.1). We expand the derivation of KSV (2019) to include any total frequency of the form (2.2) assuming (i) the inhomogeneity is weak enough, that is, $\omega _{1} \ll \omega _{0}$ (this is the weak interaction limit of statistical scattering theories such as Hasselmann (1966)); (ii) the bare frequency varies slowly over $\boldsymbol {x}$ ; and (iii) the inhomogeneity can be well modelled by a statistically homogeneous and stationary field. As in McComas & Bretherton (1977), Dong et al. (2020) and CKV (2023), we retain time dependence in $\omega _{1}$ but later will simplify to the time-independent case.

We start with the action conservation equation (2.1) with total frequency (2.2). Following § A.1 of KSV (2019), we introduce a small bookkeeping parameter $\epsilon$ , making the substitution $\omega _{1} \rightarrow \epsilon \omega _{1}$ to enforce the assumption that the perturbation terms are weak enough to be dominated by the bare frequency. We define slow time and space scales $T = \epsilon ^2 t$ , $\boldsymbol {X} = \epsilon ^2 \boldsymbol {x}$ and expand the action $a(\boldsymbol {x}, \boldsymbol {X}, \boldsymbol {k}, t, T)$ in powers of $\epsilon$ :

(3.1) \begin{equation} a = a^{(0)}(\boldsymbol {X}, \boldsymbol {k}, T) + \epsilon a^{(1)}(\boldsymbol {x}, \boldsymbol {X}, \boldsymbol {k}, t, T) + \cdots . \end{equation}

Allowing $a^{(0)}$ to vary only on slow time and space scales immediately satisfies the leading-order equation. At $O(\epsilon )$ ,

(3.2) \begin{equation} \partial _{t} a^{(1)} + c_i \partial _{x_i} a^{(1)} = \partial _{x_i} \omega _{1} \partial _{k_i} a^{(0)}, \end{equation}

where $c_i$ is the $i{\text {th}}$ component of

(3.3) \begin{equation} \boldsymbol {c} = \boldsymbol {\nabla }_{\boldsymbol {k}} \omega _{0}, \end{equation}

the (leading-order contribution to the) group velocity of the waves. This has the solution

(3.4) \begin{equation} a^{(1)}(\boldsymbol {x}, \boldsymbol {X}, \boldsymbol {k}, t, T) = \int _0^t \partial _{x_j} \omega _{1} (\boldsymbol {x} - \boldsymbol {c} s,\boldsymbol {k}, t - s) \, {\text {d}}s \, \partial _{k_j} a^{(0)}. \end{equation}

At $O(\epsilon ^2)$ , we average to eliminate $a^{(2)}$ terms and find

(3.5) \begin{equation} \partial _{T} a^{(0)} + c_i \partial _{X_i} a ^{(0)} = \langle \partial _{k_i} ( a^{(1)} \partial _{x_i} \omega _{1} ) \rangle , \end{equation}

where we have used spatial homogeneity. Unlike KSV (2019) and following Villas Bôas & Young (2020), we do not require $\partial _{x_i} \partial _{k_i} \omega _1 = 0$ , equivalent to incompressibility for $\omega _1 = \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}$ . We substitute (3.4) into (3.5), taking the upper bound of the integral to be $t \rightarrow \infty$ , appropriate for slowly evolving inhomogeneities. Then

(3.6) \begin{equation} \partial _{T} a^{(0)} + c_i \partial _{X_i} a ^{(0)} = \partial _{k_i} \left ( {{\unicode{x1D60B}}}_{ij} \partial _{k_j} a^{(0)} \right ), \end{equation}

where

(3.7) \begin{equation} {{\unicode{x1D60B}}}_{ij} = \int _0^\infty \langle \partial _{x_i} \omega _{1}(\boldsymbol {x}, \boldsymbol {k}, t) \partial _{x_j} \omega _{1}(\boldsymbol {x} - \boldsymbol {c} s, \boldsymbol {k}, t - s) \rangle \, {\text {d}}s, \end{equation}

with $\langle \cdot \rangle$ the ensemble average, are the components of the diffusivity tensor $\unicode{x1D63F}$ . Setting the bookkeeping parameter to $1$ gives the diffusion equation:

(3.8) \begin{equation} \partial _{t} a + \boldsymbol {c} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}} a = \boldsymbol {\nabla }_{\boldsymbol {k}} \boldsymbol {\cdot } \left ( {\unicode{x1D63F}} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {k}} a\right ). \end{equation}

Considering the special case of $\omega _{1} = \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}$ , the diffusivity (3.7) reduces to the McComas & Bretherton (1977) diffusivity:

(3.9) \begin{equation} {{\unicode{x1D60B}}}_{ij} = k_m k_n \int _0^\infty \langle \partial _{x_i} U_m(\boldsymbol {x}, t) \partial _{x_j} U_n(\boldsymbol {x} - \boldsymbol {c} s, t - s) \rangle \, {\text {d}}s. \end{equation}

To evaluate the diffusivity, we introduce the correlation function

(3.10) \begin{equation} \Lambda (\boldsymbol {y},\boldsymbol {k}, r)= \langle \omega _{1}(\boldsymbol {x}, \boldsymbol {k}, t) \omega _{1}(\boldsymbol {x} + \boldsymbol {y}, \boldsymbol {k}, t + r)\rangle , \end{equation}

and rewrite (3.7) as

(3.11) \begin{equation} {{\unicode{x1D60B}}}_{ij} = -\frac {1}{2}\int _{-\infty }^\infty \frac {\partial ^2 \Lambda }{\partial y_i \partial y_j}(\boldsymbol {c} s, \boldsymbol {k}, s) \, {\text {d}}s, \end{equation}

where we use $\Lambda (\boldsymbol {y}, \boldsymbol {k}, r) = \Lambda (-\boldsymbol {y}, \boldsymbol {k}, -r)$ to extend the limits of integration. Defining the Fourier-transformed correlation function $\hat {\Lambda }$ through its inverse,

(3.12) \begin{equation} \Lambda (\boldsymbol {x},\boldsymbol {k}, t) = \int _{\mathbb {R}^{n+1}} \hat {\Lambda }(\boldsymbol {K},\boldsymbol {k}, \mathit {\unicode{x1D6FA} }) {\mathrm {e}}^{\mathrm {i}(\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {x} - \mathit {\unicode{x1D6FA} } t)} \,{\text {d}}\boldsymbol {K} {\text {d}}\mathit {\unicode{x1D6FA} }, \end{equation}

the diffusivity becomes

(3.13) \begin{equation} {{\unicode{x1D60B}}}_{ij} = \pi \int _{\mathbb {R}^{n+1}} K_i K_j \hat {\Lambda }(\boldsymbol {K}, \boldsymbol {k}, \mathit {\unicode{x1D6FA} })\delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c} - \mathit {\unicode{x1D6FA} } ) \, {\text {d}}\boldsymbol {K} {\text {d}}\mathit {\unicode{x1D6FA} }, \end{equation}

where we have used that $\int _{\mathbb {R}} \exp {(\mathrm {i}(\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c} - \mathit {\unicode{x1D6FA} } )s)} {\text {d}}s = 2 \pi \delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c} - \mathit {\unicode{x1D6FA} } )$ . For time-independent $\omega _{1}$ , this reduces to

(3.14) \begin{equation} {{\unicode{x1D60B}}}_{ij} = \pi \int _{\mathbb {R}^{n}} K_i K_j \hat {\Lambda }(\boldsymbol {K}, \boldsymbol {k})\delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}) \, {\text {d}}\boldsymbol {K}. \end{equation}

Next, we evaluate this expression for our two example systems and, in the shallow-water case, support our findings with ray tracing simulations. In the Boussinesq case, we find the forced equilibrium wave energy spectrum.

3.1. Shallow water

To evaluate the shallow-water diffusivity, we first evaluate the correlation function (3.10), then find its Fourier transform and substitute into (3.13). We introduce the background flow stream function $\psi$ such that

(3.15a,b) \begin{equation} \boldsymbol {U} = \boldsymbol {\nabla }_{\boldsymbol {x}}^{\perp } \psi = (-\partial _{x_2} \psi , \partial _{x_1} \psi , 0) \quad \textrm {and} \quad g \Delta H = f \psi , \end{equation}

where $\boldsymbol {\nabla }_{\boldsymbol {x}}^{\perp } = (-\partial _{x_2},\partial _{x_1})$ is the skew gradient and we use geostrophic balance (2.3). Substituting this and the frequency (2.5) into (3.10) yields

(3.16) \begin{equation} \Lambda (\boldsymbol {y},\boldsymbol {k}, r)= \underbrace {\vphantom {\frac {f^2 k_h^4}{4 \omega _{0}^2}}-\big(\boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {y}}^{\perp }\big)^2\langle \psi (\boldsymbol {x}, t) \psi (\boldsymbol {x} + \boldsymbol {y}, t + r)\rangle }_{\text{Doppler shift}} \, + \, \underbrace {\frac {f^2 k_h^4}{4 \omega _{0}^2} \langle \psi (\boldsymbol {x}, t)\psi (\boldsymbol {x} + \boldsymbol {y}, t + r)\rangle }_{\text{height fluctuation}}. \end{equation}

For the Doppler shift term, the skew gradients have been moved through the ensemble average using integration by parts and the symmetry of $\boldsymbol {x}$ and $\boldsymbol {y}$ arguments. Notably, there are no cross terms in (3.16) and the height fluctuation and Doppler shift effects are uncoupled. This is because, using integration by parts under the ensemble average,

(3.17) \begin{align} \text {cross terms} &\propto \langle \psi (\boldsymbol {x} + \boldsymbol {y}, t + r) \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}}^{\perp } \psi (\boldsymbol {x}, t) \rangle + \langle \psi (\boldsymbol {x}, t) \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}}^{\perp } \psi (\boldsymbol {x} + \boldsymbol {y}, t + r) \rangle \notag \\ &\propto - \langle \psi (\boldsymbol {x}, t) \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}}^{\perp } \psi (\boldsymbol {x} + \boldsymbol {y}, t + r) \rangle + \langle \psi (\boldsymbol {x}, t) \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}}^{\perp } \psi (\boldsymbol {x} + \boldsymbol {y}, t + r) \rangle \notag \\ &= 0. \end{align}

The Fourier transform of the correlation function (3.16) is

(3.18) \begin{equation} \hat {\Lambda }(\boldsymbol {K}, \boldsymbol {k}, \mathit {\unicode{x1D6FA} }) = |\boldsymbol {k} \times \boldsymbol {K}|^2 E_\psi (\boldsymbol {K}, \mathit {\unicode{x1D6FA} }) \, + \, \frac {f^2 k_h^4}{4 \omega _{0}^2} E_\psi (\boldsymbol {K}, \mathit {\unicode{x1D6FA} }), \end{equation}

where

(3.19) \begin{equation} E_\psi (\boldsymbol {K}, \mathit {\unicode{x1D6FA} }) = \langle \hat {\psi }(-K, -\mathit {\unicode{x1D6FA} }) \hat {\psi }(K, \mathit {\unicode{x1D6FA} })\rangle \end{equation}

is the spectrum of the stream function. Defining the horizontal energy spectrum of the background flow as

(3.20) \begin{equation} E(\boldsymbol {K}, \mathit {\unicode{x1D6FA} }) = K_h^2 E_\psi /2, \end{equation}

we rewrite (3.18) as

(3.21) \begin{equation} \hat {\Lambda }(\boldsymbol {K}, \boldsymbol {k}, \mathit {\unicode{x1D6FA} }) = 2 k_h^2 \sin ^2 \gamma E(\boldsymbol {K}, \mathit {\unicode{x1D6FA} }) + \frac {f^2 k_h^4}{2 K_h^2 \omega _{0}^2} E(\boldsymbol {K}, \mathit {\unicode{x1D6FA} }), \end{equation}

where $K_h$ is the horizontal wavenumber of the background flow and $\gamma$ is the angle between $\boldsymbol {K}$ and $\boldsymbol {k}$ . Combined with (3.13), we have a diffusivity that accounts for height fluctuation effects for a time-dependent flow.

As proof of concept, we simplify to a time-independent flow. This (i) allows for inexpensive ray tracing simulations in § 3.2 and (ii) results in a simpler form for the diffusivity. Adding time dependence is possible, as done solely for the Doppler shift effect by Dong et al. (2020) and – for the full 3-D set-up – by CKV (2023). Combining the time-independent-flow limit of (3.21) with (3.14), we have

(3.22) \begin{align} {{\unicode{x1D60B}}}_{ij} = \, &2 \pi k_h^2 \int _{0}^{\infty } {\text {d}}K_h\int _{0}^{2\pi } {\text {d}}\gamma \, K_h K_i K_j \sin ^2 \gamma E(\boldsymbol {K}) \delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c})\notag \\ &+ \frac {\pi f^2 k_h^4}{2 \omega _{0}^2} \int _{0}^{\infty } {\text {d}}K_h\int _{0}^{2\pi } {\text {d}}\gamma \, \frac {K_i K_j}{K_h} E(\boldsymbol {K})\delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}), \end{align}

in polar coordinates $\boldsymbol {K} = (K_h, \gamma )$ . Here $E(\boldsymbol {K})$ is the flow kinetic energy spectrum marginalised over frequencies. In the local polar basis $(\boldsymbol {e}_{k_h}, \boldsymbol {e}_{\phi })$ associated with $\boldsymbol {k}$ , this diffusivity has one non-zero component:

(3.23) \begin{equation} {{\unicode{x1D60B}}}_{\phi \phi } = \boldsymbol {e}_{\phi } \boldsymbol {\cdot } {\unicode{x1D63F}} \boldsymbol {\cdot } \boldsymbol {e}_{\phi }. \end{equation}

This is because the constant frequency surface in spectral space is the circle $\omega = \omega _{0}(k_h) + O(\epsilon ) = \text{const}$ . There is no diffusion of action across this circle because for a time-independent linear system, wave frequency does not change. This means there are no radial components of the diffusivity, only (3.23). We see this explicitly by expanding $\boldsymbol {K}$ as

(3.24) \begin{equation} \boldsymbol {K} = K_h \cos \gamma \boldsymbol {e}_{k_h} + K_h \sin \gamma \boldsymbol {e}_{\phi }. \end{equation}

Then, noting $\delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}) \propto \delta (\cos \gamma )$ as $\boldsymbol {c} = \boldsymbol {\nabla }_{\boldsymbol {k}} \omega _0 (k_h) = c \boldsymbol {e}_{k_h}$ shows that any integrand containing $\cos \gamma$ integrates to zero (providing the energy spectrum $E(\boldsymbol {K})$ is well behaved). This is the case for all components of (3.22) bar (3.23), which, using (3.24), is

(3.25) \begin{align} {{\unicode{x1D60B}}}_{\phi \phi } = \, &2 \pi k_h^2 \int _{0}^{\infty } {\text {d}}K_h\int _{0}^{2\pi } {\text {d}}\gamma \, K_h^3 \sin ^4 \gamma E(\boldsymbol {K}) \delta (K_hc\cos \gamma )\notag \\ &+ \frac {\pi f^2 k_h^4}{2 \omega _{0}^2} \int _{0}^{\infty } {\text {d}}K_h\int _{0}^{2\pi } {\text {d}}\gamma \, K_h \sin ^2 \gamma E(\boldsymbol {K})\delta (K_hc \cos \gamma ). \end{align}

Assuming flow isotropy such that $E(\boldsymbol {K}) = E(K_h)$ , we integrate in $\gamma$ to get

(3.26) \begin{equation} {{\unicode{x1D60B}}}_{\phi \phi } = \underbrace {\frac {4 \pi k_h^2}{c \vphantom {\omega ^2_0}} \int _{0}^{\infty } K_h^2 E(K_h) \, {\text {d}}K_h}_{\text{Doppler shift}} \, + \, \underbrace {\frac {\pi f^2 k_h^4}{c \omega _{0}^2} \int _{0}^{\infty } E(K_h) \, {\text {d}}K_h}_{\text{height fluctuation}}. \end{equation}

With this single component, the diffusion equation (3.8) reduces to

(3.27) \begin{equation} \partial _{t} a + \boldsymbol {c} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}} a = \mu \partial _{\phi \phi } a, \end{equation}

where we define the directional diffusivity

(3.28) \begin{equation} \mu = {{\unicode{x1D60B}}}_{\phi \phi } / k_h^2. \end{equation}

To examine the relative importance of each term in (3.26), we consider the ratio between the two:

(3.29) \begin{equation} R^{\prime}_{{sw}} = \frac {\left [{{\unicode{x1D60B}}}_{\phi \phi }\right ]_{\text{height fluctuation}}}{\left [{{\unicode{x1D60B}}}_{\phi \phi }\right ]_{\text{Doppler shift}}} = \frac {f^2 k_h^2}{4 \omega _{0}^2} \frac {\int _{0}^{\infty } E(K_h) \, {\text{d}}K_h}{\int _{0}^{\infty } K_h^2 E(K_h) \, {\text {d}}K_h} \sim \frac {(k_h/K_*)^2}{4 (1 + Bu \, (k_h/K_*)^2)}. \end{equation}

Here, we approximate the integrals by

(3.30a,b) \begin{equation} 2 \pi \int _{0}^{\infty } K_h^2 E(K_h) \, {\text {d}}K_h \sim K_* U_*^2/2 \quad \text {and} \quad 2 \pi \int _{0}^{\infty } E(K_h) \, {\text {d}}K_h \sim U_*^2/(2 K_*) \end{equation}

using the characteristic speed and horizontal wavenumber of the flow. The final expression in (3.29) is the square of the scaling argument estimate $R_{{sw}}$ (2.8) which comes from the wave dispersion relation. This is sensible – the diffusivity (3.7) consists of a square frequency term. This means that in the diffusion regime, the scaling argument of § 2.1 is an underestimate for $R_{{sw}} \gt 1$ and an overestimate for $R_{{sw}} \lt 1$ . In figure 2, the adjustment from $R_{{sw}}$ to $R^{\prime}_{{sw}}$ is shown.

3.2. Shallow-water ray tracing

The 2-D diffusion equation (3.27) has an exact solution as outlined in § 4 of Villas Bôas & Young (2020). We use this solution and ray tracing – a numerical method to find constant-wave-action trajectories – to assess the validity of the diffusion approximation and form of the diffusivity (3.26).

We start by multiplying the diffusion equation by $\cos \phi$ and integrating over the entire $(x, y)$ plane and $\phi \in (-\pi , \pi )$ to get

(3.31) \begin{equation} \frac {\text {d}}{{\text {d}}t} \iiint \cos \phi \, a \, {\text {d}}x\, {\text {d}}y\, {\text {d}}\phi = -\mu \iiint \cos \phi \, a \, {\text {d}}x\, {\text {d}}y\, {\text {d}}\phi . \end{equation}

Note that without pre-multiplying by $\cos \phi$ , the left-hand side is zero by action conservation. Here, the $\boldsymbol {c} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {x}} a$ term disappears under the $(x,y)$ integrals because the $a$ we consider is ensemble-averaged and we assume statistical spatial homogeneity in deriving the diffusion equation. Integrating with respect to time yields

(3.32) \begin{equation} \ln {(\langle \cos \phi \rangle _a)} = - \mu t + \ln {(\langle \cos \phi _0 \rangle _a)}, \end{equation}

where $\langle \cdot \rangle _a$ denotes the action-weighted ensemble average and $\phi _0$ is the initial angle. This is an alternative form of (4.3) in Villas Bôas & Young (2020).

We confirm that (3.32) holds by ray tracing: solving the characteristic equations of the action conservation equation (2.1) numerically to give many constant-action trajectories and taking an ensemble average. These ray equations are

(3.33a,b) \begin{equation} \partial _t \boldsymbol {x} = \boldsymbol {\nabla }_{\boldsymbol {k}} \omega \quad \text {and} \quad \partial _t \boldsymbol {k} = -\boldsymbol {\nabla }_{\boldsymbol {x}} \omega . \end{equation}

We define non-dimensional quantities

(3.34a–d) \begin{equation} t^{\prime} = t f, \quad \boldsymbol {k}^{\prime} = \boldsymbol {k}/K_*, \quad \boldsymbol {x}^{\prime} = \boldsymbol {x} K_* \quad \text {and} \quad \boldsymbol {U}^{\prime} = \boldsymbol {U}/U_* \end{equation}

(and $\psi^{\prime} = \psi K_*/U_*$ ). Substituting the approximate frequency for shallow-water waves (2.4) into (3.33a,b ) gives

(3.35a) \begin{align} & \partial _{t^{\prime}} \boldsymbol {x}^{\prime} =\ \ \frac {Bu}{{(1 + Bu k_h^{\prime 2})}^{1/2}}\boldsymbol {k}^{\prime} + \quad Ro \boldsymbol {U}^{\prime} \quad\quad \quad \quad +\quad Ro \frac {2 + Bu k_h^{\prime 2}}{2 {(1 + Bu k_h^{\prime 2})}^{3/2}} \psi^{\prime} \boldsymbol {k}^{\prime} , \end{align}

(3.35b) \begin{align} & \partial _{t^{\prime}} \boldsymbol {k}^{\prime} = \ \ 0 \qquad\qquad \quad\quad\ \ - \quad Ro \boldsymbol {\nabla }_{\boldsymbol {x}^{\prime}}\left (\boldsymbol {U}^{\prime} \boldsymbol {\cdot } \boldsymbol {k}^{\prime}\right ) \ - \quad Ro \frac {k_h^{\prime 2}}{2 (1 + Bu k_h^{\prime 2}) }\boldsymbol {U}^{\prime \perp },\\ & \qquad \quad \ \underbrace {\hphantom {\frac {Bu}{(1 + Bu k_h^{\prime 2})^{1/2}}\boldsymbol {k}^{\prime}}}_{\text{no flow}} \quad \quad \underbrace {\hphantom {Ro \boldsymbol {\nabla }_{\boldsymbol {x}^{\prime}}\left (\boldsymbol {U}^{\prime} \boldsymbol {\cdot } \boldsymbol {k}^{\prime}\right )}}_{\text{Doppler shift}} \quad \quad \underbrace {\hphantom {Ro \frac {2 + Bu k_h^{\prime 2}}{2 {(1 + Bu k_h^{\prime 2})}^{3/2}} \psi^{\prime} \boldsymbol {k}^{\prime}}}_{\text{height fluctuation}}\nonumber \end{align}

where

(3.36) \begin{equation} Ro = \frac {U_*K_*}{f} \end{equation}

is the flow Rossby number and we introduce $\boldsymbol {U}^\perp = (V, -U, 0)$ .

The rays propagate in a flow which is a snapshot of a fully evolved 2-D quasi-geostrophic Navier–Stokes simulation, solved using a pseudo-spectral method and Crank–Nicholson time-stepping with time step ${\text {d}}t = 0.01$ . The doubly periodic domain $\boldsymbol {x} \in [0, L]$ is discretised by $1024^2$ gridpoints. Viscosity is $\nu = 10^{-5}$ . The initial wavenumber distribution of the flow is Gaussian. We prescribe nine values of the Rossby radius of deformation $L_D = Bu^{1/2}/K_*$ and three values of the initial wavenumber $K_*$ to give $3^3$ flows. This gives final snapshots with varying $Ro$ and $Bu$ . A typical flow is shown in figure 4(a).

We perform ray simulations with the inhomogeneity $\omega _{1}$ given by (i) Doppler shift only; (ii) height fluctuation only; and (iii) both Doppler shift and height fluctuation.

We initialise $50^2$ rays equally distributed across the periodic $\boldsymbol {x}$ domain with a constant horizontal wavenumber $k_h \approx 32.2 K_*$ (within the WKB regime) and random initial angle. By giving the rays different initial angles, our rays sample more of the flow and the ensemble average approaches the exact solution more rapidly as the number of rays increases than it would with a constant initial angle. For simplicity, we consider the angle $\phi$ to be the angle of deviation from the ray’s initial angle meaning $\ln {(\langle \cos \phi _0 \rangle _a)} = 0$ . This is possible because of the flow’s statistical isotropy. The angle of deviation for each ray is calculated at each time step, and $\cos \phi$ is averaged over. We choose action $a=1$ for each ray so that $\langle \cdot \rangle _a$ and the ensemble average are identical. Some of the $3^3$ flows have very similar Burger or Rossby numbers; we average over these Burger and Rossby numbers to leave $3^2$ results with distinct Burger and Rossby numbers. We estimate $\mu$ , (3.28), the negative gradient of $\ln {(\langle \cos \phi \rangle _a)}$ against time $t$ with ${{\unicode{x1D60B}}}_{\phi \phi }$ computed from (3.26). We compare this with the simulation gradient.

An example ray tracing simulation is shown in figure 4 in physical and wavevector space. The Rossby number of the flow is the small parameter in § 3, $\epsilon \sim Ro \sim O(0.01)$ , and so the flow is weak and the rays are only slightly deflected from their initial angle of propagation. In $\boldsymbol {k}$ space, the small Rossby number characterises the thickness of the constant-frequency ring. The exact solution and simulation approximation of (3.32) are shown in figure 4(c). This confirms the validity of the diffusion approximation and of the formulas (3.26). As expected from these formulas, the height fluctuations and Doppler shift make comparable contributions to the spectral diffusivity.

Figure 4. Ray trajectories satisfying the characteristic equations (3.35a , b ), propagating in the flow described in § 3.2 with $Ro = 0.03$ , $Bu = 0.32$ and $K_* = 3.10$ . (a) A sample of 10 rays in physical space superimposed on the flow vorticity field (darker sections indicate higher-magnitude vorticity). (b) A sample of 50 rays in spectral space initialised on a constant-frequency circle given by the dotted black line. (c) The exact solution (3.32) approximated by the rays’ ensemble average (solid lines) and the diffusivity (3.26) (dashed lines). Red lines are calculated solely with the height fluctuation terms in (3.26) and (3.35a , b ); green lines are solely Doppler shift and; blue lines are the entire system with both Doppler and height fluctuation effects. The vertical dashed line indicates $\tau$ , the point at which the gradient calculation for figure 5 begins.

The $Ro$ and $Bu$ dependence of $\mu$ can be seen from (3.26). We define the non-dimensional directional diffusivity $\tilde {\mu }$ as

(3.37) \begin{equation} \tilde {\mu } = \frac {Bu^{1/2 }}{f Ro^2}\mu \sim \underbrace {\frac {(\omega _{0}/f)}{Bu^{1/2} (k_h/K_*)}}_{\text{Doppler shift}} \, + \, \underbrace {\frac { (k_h/K_*) }{4 Bu^{1/2} (\omega _{0}/f)}}_{\text{height fluctuation}} \rightarrow \underbrace {\vphantom {\frac {(\omega _{0}/f)}{Bu^{1/2} (k_h/K_*)}}1}_{\text{D.s.}} \,+\, \underbrace {\vphantom {\frac {(\omega _{0}/f)}{Bu^{1/2} (k_h/K_*)}}(4Bu)^{-1}}_{\text{h.f.}}, \end{equation}

where the final limit holds for $(k_h/K_*)^2 \gg 1$ and we have used (3.30a,b ). The scaling of $\mu$ is chosen so that the Doppler shift component of $\tilde {\mu }$ is independent of $Bu$ and $Ro$ . Note that the ratio between height fluctuation and Doppler shift effects is in agreement with $R^{\prime}_{{sw}}$ , (3.29). In figure 5, $\tilde {\mu }$ is displayed for a range of Burger numbers. The Burger dependence of (3.37) is confirmed for each contribution and, because three separate Rossby numbers are used in the plot, so too is the Rossby dependence.

The characteristic flow velocity $U_*$ is chosen such that (3.30a,b ) holds exactly and $\tilde {\mu } = 1$ is the exact limit of the non-dimensionalised, Doppler shift directional diffusivity. Then the height fluctuation $\tilde {\mu }$ shown in figure 5 coincides with the ratio between the height fluctuation and Doppler shift effects $R^{\prime}_{{sw}}$ . We see good agreement between the height fluctuation $\tilde {\mu }$ computed exactly and from simulations and thus confirm our estimate of $R^{\prime}_{{sw}}$ given by (3.29).

In deriving the diffusivity (3.7), we approximate an integral between $0$ and $t$ by one between $0$ and $\infty$ and so there is a delay between the start of the ray simulation to the point at which the exact solution (3.32) is well approximated. Therefore, our $\mu$ estimates in figure 5 begin a short time $\tau /f$ after the simulations have begun, as indicated in figure 4(c).

Figure 5. Non-dimensional directional diffusivity $\tilde {\mu }$ , as defined by (3.28) and (3.37), plotted against flow Burger number $Bu$ . Red lines correspond to the height fluctuation component of the directional diffusivity, green lines to the Doppler shift component and blue lines to the full directional diffusivity. Crosses indicate the estimates of $\tilde {\mu }$ from the $50^2$ -ray simulations and exact solution (3.32), pluses are computed from the diffusivity expression (3.26) and dotted lines are the power laws predicted by (3.37), fitted to the first of the exact data points. The nine different Burger numbers cover three sets of Rossby number, labelled $1$ ( $Ro \approx 0.023$ ), $2$ ( $Ro \approx 0.030$ ) and $3$ ( $Ro \approx 0.057$ ). The Doppler shift directional diffusivity $\tilde {\mu } \approx 1$ and so the height fluctuation directional diffusivity corresponds to the ratio between the two, $R^{\prime}_{{sw}}$ (3.29).

We stress that the good agreement in figures 4(c) and 5 between the $\mu$ of (3.28) computed from the diffusivity (3.26) and the $\mu$ of (3.32) estimated from the ensemble average of rays validates the diffusion approximation (3.27) of the action conservation equation (2.1) with inhomogeneities of the form (2.5).

3.3. Three-dimensional Boussinesq

The method to evaluate the general diffusivity (3.13) for a Boussinesq system with vertical buoyancy gradients is similar to that for the shallow-water system. We give the result here, providing a full derivation in Appendix B.1. For a Doppler-shift-induced diffusivity, CKV (2023) show that it is well justified to assume a slowly evolving flow is time-independent. The same likely applies to a buoyancy-gradient-induced diffusivity, so we focus on the time-independent case.

The two non-zero components of the diffusivities for a Boussinesq system with vertical buoyancy gradients and dispersion relation (2.15) are

(3.38) \begin{align} & \displaystyle\left.\begin{array}{ll}{{\unicode{x1D60B}}}_{kk} & = \dfrac {4\pi k^3 \omega _{0} \sin ^2 \theta }{(\bar {N}^2 - f^2) |\cos ^5 \theta |}\nonumber\\[10pt]& \qquad \displaystyle\times \int _{0}^{\infty } \int _{\theta }^{\pi - \theta } K^3 \cos ^2 \mathit {\unicode{x1D6E9} } (\cot ^2 \theta - \cot ^2 \mathit {\unicode{x1D6E9} })^{1/2} E(K, \mathit {\unicode{x1D6E9} }) \, {\text {d}}K {\text {d}}\mathit {\unicode{x1D6E9} } \end{array} \right\} {\text{ Doppler shift}}\\[2pt] & \quad \displaystyle\left.\begin{array}{ll} & \quad\ \ + \dfrac {\pi k f^2 \sin ^2 \theta }{ \omega _{0} (\bar {N}^2 - f^2)|\cos ^3 \theta |}\\[10pt] & \!\!\qquad\quad \displaystyle\times \int _{0}^{\infty } \int _{\theta }^{\pi - \theta } \dfrac {K^5 \cos ^6 \mathit {\unicode{x1D6E9} } E(K, \mathit {\unicode{x1D6E9} })}{ \sin ^2 \mathit {\unicode{x1D6E9} } (\cot ^2 \theta - \cot ^2 \mathit {\unicode{x1D6E9} })^{1/2}} \, {\text {d}}K {\text {d}}\mathit {\unicode{x1D6E9} } \end{array} \right\} \text{buoyancy fluctuation} \end{align}

and

(3.39) \begin{align} & \!\!\displaystyle\left.\begin{array}{ll}{{\unicode{x1D60B}}}_{\phi \phi } & = \dfrac {4\pi k^3 \omega _{0} \sin ^4 \theta }{(\bar {N}^2 - f^2) |\cos ^5 \theta |}\\& \quad\displaystyle\times \int _{0}^{\infty } \int _{\theta }^{\pi - \theta } K^3 \sin ^2 \unicode{x1D6E9} (\cot ^2 \theta - \cot ^2 \mathit {\unicode{x1D6E9} })^{3/2}E(K, \mathit {\unicode{x1D6E9} })\, {\text {d}}K {\text {d}}\mathit {\unicode{x1D6E9} }\end{array} \right\} {\text{ Doppler shift}} \\[3pt] & \quad \displaystyle\left.\begin{array}{ll}& \quad\ + \dfrac {\pi k f^2 \sin ^4 \theta }{ \omega _{0} (\bar {N}^2 - f^2) |\cos ^3 \theta |}\\[9pt]& \!\!\!\qquad \displaystyle\times\! \int _{0}^{\infty }\!\! \int _{\theta }^{\pi - \theta }\!\!\! K^5 \cos ^4 \mathit {\unicode{x1D6E9} } (\cot ^2 \theta {-} \cot ^2 \mathit {\unicode{x1D6E9} })^{1/2} E(K,\mathit {\unicode{x1D6E9} }) {\text {d}}K {\text {d}}\mathit {\unicode{x1D6E9} } \end{array} \!\right\}\!\text{buoyancy fluctuation}, \end{align}

where we use polar coordinates $\boldsymbol {K} = (K, \gamma , \mathit {\unicode{x1D6E9} })$ for the background flow and $E$ is the flow kinetic energy spectrum, as defined in the shallow-water case (3.19)–(3.20). Here, $\gamma = \mathit {\Phi } - \phi$ is the angle between the horizontal components of $\boldsymbol {K}$ and $\boldsymbol {k}$ , with $\mathit {\Phi }$ the azimuthal angle of $\boldsymbol {K}$ . The Doppler shift and buoyancy fluctuation diffusivities are uncoupled, as with the inhomogeneities in the shallow-water system. This is true of a time-dependent background flow also (see Appendix B.1).

Only the kinetic energy spectrum of the flow appears in (3.38)–(3.39) because we assume thermal wind balance. More generally, the diffusivity will be expressed in terms of both potential and kinetic energy spectra.

The $(\cot ^2\theta -\cot ^2\mathit {\unicode{x1D6E9} })^{-1/2}$ factor in the integrand of the buoyancy fluctuation term of ${{\unicode{x1D60B}}}_{kk}$ behaves like $\delta ^{-1/2}$ a small distance $\delta$ from the singularities at $\mathit {\unicode{x1D6E9} } = \theta , \pi -\theta$ . Therefore, it is integrable and does not cause the diffusivity to diverge.

Only the part of the flow spectrum with wavevectors of polar angle $\mathit {\unicode{x1D6E9} } \in (\theta , \pi - \theta )$ contributes to the diffusivity integrals. Inertia–gravity wave diffusion is a subregime of IGW scattering. This triadic interaction occurs between an incident wave $\boldsymbol {k}$ and the background flow $\boldsymbol {K}$ , resulting in a scattered wave $\boldsymbol {k}^{\prime}$ of the same frequency. In wavevector space, $\boldsymbol {k} + \boldsymbol {K} = \boldsymbol {k}^{\prime}$ . As the incident wave and scattered wave are of the same frequency, $\boldsymbol {K}$ connects two points on the cone of constant frequency. Thus, the range of $\boldsymbol {K}$ ’s polar angle $\mathit {\unicode{x1D6E9} }$ is the range of angles connecting any two points on the cone, i.e. $(\theta , \pi - \theta )$ . This is demonstrated in figure 6. Significantly, waves of frequency $\omega _{0}(\theta )$ will not be scattered by flow fluctuations with wavevectors outside of these angle bounds, regardless of the flow’s energy.

Figure 6. Schematic of the scattering interaction between an incident wave $\boldsymbol {k}$ and the background flow $\boldsymbol {K}$ resulting in a scattered wave $\boldsymbol {k}^{\prime}$ . The two waves have the same frequency $\omega _{0}(\theta )$ and thus the background flow connects two points on the constant-frequency cone. For an arbitrary scattering interaction, angle $\mathit {\unicode{x1D6E9} }$ of the flow’s wavevector is bounded between $\theta$ and $\pi - \theta$ , i.e. the angular range of a vector between any two points on the cone.

These invisible flow regions make it difficult to apply the scaling argument of § 2.2. The polar angle of the characteristic wavevector of the flow $\boldsymbol {K}_*$ does not necessarily lie within $(\theta , \pi -\theta )$ , in which case $R_{{B}}$ (2.18) computed at $\boldsymbol {K}_*$ is not meaningful. Instead, a dominant wavevector – the characteristic wavevector of the portion of the flow which scatters the waves – must be used.

We define ${R^{\prime}_{{B}}}$ as the ratio between the buoyancy fluctuation and Doppler shift diffusivites for radial and azimuthal components:

(3.40a,b) \begin{equation} {R^{\prime}_{B,kk}} = \frac {[{{\unicode{x1D60B}}}_{kk}]_{\text{buoyancy fluctuation}}}{[{{\unicode{x1D60B}}}_{kk}]_{\text{Doppler shift}}} \quad \text {and} \quad R^{\prime}_{B,\phi \phi } = \frac {[{{\unicode{x1D60B}}}_{\phi \phi }]_{\text{buoyancy fluctuation}}}{[{{\unicode{x1D60B}}}_{\phi \phi }]_{\text{Doppler shift}}}. \end{equation}

Both ratios scale like $(k_h/K_*)^{-2}$ through the diffusivity prefactors’ $k$ dependence. The frequency behaviour is more complicated, requiring the diffusivity integrals to be evaluated. In Appendix B.2, we show that ${R^{\prime}_{{B},kk}} \rightarrow 0$ for $\omega _{0}/f \rightarrow 1$ and for larger frequencies and high aspect ratios $\alpha$ (2.16), ${R^{\prime}_{{B},kk}} \sim (\omega _{0}/f)^{-2}$ . This means that for waves in the WKB regime with high aspect ratios and frequencies, vertical buoyancy gradients can be neglected.

We compute the ratio of diffusivity components for a typical quasi-geostrophic flow and find good agreement with the $-2$ power laws. The geostrophic energy spectrum used is extracted from a snapshot of the full Boussinesq simulation described in CKV (2023) and is pictured in figure 7. For this spectrum, $N/f = 32.0$ and the aspect ratio, (2.16), $\alpha \approx 16.0$ .

Figure 7. The geostrophic flow component $E$ , smoothed and scaled by a maximum value $E_M$ , of the full Boussinesq simulation in CKV (2023). The horizontal and vertical flow wavenumbers ( $K_h, K_v$ ) are scaled by the characteristic wavevector $(K_h, K_v) = (K_*, \alpha K_*)$ – with $\alpha$ the aspect ratio of the flow (2.16) – marked by a white cross, at which the geostrophic energy is maximum.

Both ratios of diffusivities are shown in figure 8, computed with the energy spectrum of figure 7. In figure 9, for the radial diffusivity ratio, we plot cross-sections of constant frequency and horizontal wavenumber and find good agreement with the $-2$ power laws.

Figure 8. Radial and azimuthal ratios $(a)$ $R^{\prime}_{{B},kk}$ and $(b)$ $R^{\prime}_{{B},\phi \phi }$ as given in (3.40a,b ), the ratio of buoyancy fluctuation and Doppler shift diffusivities for a snapshot of the full Boussinesq simulation of CKV (2023), against non-dimensionalised frequency $\omega _0/f$ and horizontal wavenumber $k_h/K_*$ . Contours are shown for ${R_{{B}}}' = 0.01$ (dashed black), $0.1$ (solid black) and $1$ (dashed white).

Figure 9. Cross-sections of radial ratio $R^{\prime}_{{B},kk}$ (3.40a ), as shown in figure 8, for constant $(a)$ $\omega _{0}/f$ and $(b)$ $k_h/K_*$ . The $(\omega _{0}/f/\cos ^2\theta )^{-2}$ prediction of Appendix B.2 is, for small $\theta$ , an $(\omega _{0}/f)^{-2}$ power law.

The values of ${{R_{{B}}}}_{,kk}$ and ${{R_{{B}}}}_{,\phi \phi }$ are, for the WKB regime of $k_h/K_* \gg 1$ , $\lesssim 0.1$ . Thus, at least for the waves in CKV (2023), we predict that a weak vertical buoyancy gradient induces negligible spectral diffusion in the WKB regime.

Figure 8, a comparison of diffusivity magnitudes, is markedly different from the comparison of dispersion relation terms (figure 3). Although both predict a decrease in the buoyancy fluctuation effect as horizontal wavenumber increases, the scaling argument of § 2.2 predicts that the effect dominates for higher frequencies. The $(\omega _{0}/f)^{-2}$ power law predicted from diffusivities and confirmed for an example spectrum in figure 9 contradicts this. We attribute this conflicting prediction to the large part of the flow spectrum not contributing to the diffusivity at higher wave frequencies. This makes the scaling arguments of § 2.2 unreliable and calls for the exact evaluation of the diffusivities (3.38)–(3.39).

3.4. Forced equilibrium spectrum

We solve the diffusion equation (3.8) exactly in the time-independent case corresponding to the equilibrium spectrum obtained under constant forcing. Our aim is to assess how the $k^{\pm 2}$ power-law spectra obtained by KSV (2019) when diffusion is solely caused by Doppler shift are modified when accounting for vertical buoyancy gradients. As in KSV (2019), we focus on radial diffusion, assuming wave statistics independent of $\phi$ such that $\partial _{\phi } a = 0$ . To concisely contrast this work to theirs, we consider only the forced, equilibrium solution to (3.8). We consider the energy density for ease of interpretation, defined as $e(k;\theta ) = 2\pi k^2 \sin \theta \omega a(k;\theta )$ , such that $e \, \text {d}k\, \text {d}\theta$ is the energy contained in the box $[k, k + \text {d}k]$ and $[\theta , \theta + \text {d} \theta ]$ . Thus, ignoring unimportant prefactors on the right-hand side, (3.8) becomes

(3.41) \begin{equation} \partial _{k} \left (k^2 {{\unicode{x1D60B}}}_{k k} \partial _{k} \frac {e}{k^2} \right ) = - \delta (k - k_*). \end{equation}

The forcing in a circle at $k = k_*$ can be generalised via integration. The minus sign ensures a positive energy spectrum. Defining

(3.42) \begin{equation} Q = k^3[{{\unicode{x1D60B}}}_{kk}]_{\text{Doppler shift}} \quad \text {and} \quad P = k[{{\unicode{x1D60B}}}_{kk}]_{\text{buoyancy fluctuation}} \end{equation}

as the $k$ -independent parts of the Doppler shift and buoyancy fluctuation diffusivities, this simplifies to

(3.43) \begin{equation} Q\partial _{k} \left (\left ( k^5 + \beta k^3\right ) \partial _{k} \frac {e}{k^2} \right ) = - \delta (k - k_*), \end{equation}

where $\beta = P/Q$ . This equation has the piece-wise solution found, for example, by partial fractions:

(3.44) \begin{align} e(k) =\left \{ \begin{array}{ll} A\left (1 - \frac {k^2}{\beta } \ln {\left (1 + \frac {\beta }{k^2}\right )}\right ) + B k^2 & \textrm {for}\ 0 \lt k \lt k_* ,\\[3pt] C\left (1 - \frac {k^2}{\beta } \ln {\left (1 + \frac {\beta }{k^2}\right )}\right ) + D k^2 & \textrm {for} \ k \gt k_* .\end{array} \right . \end{align}

We require $e(k)$ is bounded as $k \rightarrow \infty$ which means $D = 0$ . We require zero energy at $k=0$ ; therefore $A = 0$ . Continuity at $k = k_*$ gives $B = C(1/k_*^2 - \ln {(1 + \beta /k_*^2)}/\beta )$ . The jump condition $[Q (k^5 + \beta k^3 ) \partial _{k} \frac {e}{k^2}]_{k_*^-}^{k_*^+} = -1$ gives $C = 1/(2\beta Q)$ . Thus,

(3.45) \begin{align} e(k) =\frac {1}{2 \beta Q}\left \{ \begin{array}{ll} \frac {k^2}{k_*^2} \left (1 - \frac {k_*^2}{\beta } \ln {\left (1 + \frac {\beta }{k_*^2}\right )}\right )& \textrm {for}\ 0 \lt k \lt k_* ,\\ 1 - \frac {k^2}{\beta } \ln {\left (1 + \frac {\beta }{k^2}\right )} & \textrm {for} \ k \gt k_* .\end{array} \right . \end{align}

Note that $Q \rightarrow 0$ is an artificial singularity introduced by the choice of factorisation and forcing in (3.43). The Doppler shift limit of $\beta \rightarrow 0$ is

(3.46) \begin{align} e(k) =\frac {1}{4 Q k_*^2}\left \{ \begin{array}{ll} (k/k_*)^2 & \textrm {for}\ 0 \lt k \lt k_* ,\\ (k_*/k)^2 & \textrm {for} \ k \gt k_* , \end{array} \right . \end{align}

which is exactly (4.1) of KSV (2019). The buoyancy fluctuation limit of $\beta \rightarrow \infty$ is

(3.47) \begin{align} e(k) =\frac {1}{2 P}\left \{ \begin{array}{ll} (k/k_*)^2 & \textrm {for}\ 0 \lt k \lt k_* ,\\ 1 & \textrm {for} \ k \gt k_* .\end{array} \right . \end{align}

In figure 10, (3.45)–(3.47) are displayed for two values of $\beta /k_*^2$ .

Figure 10. Forced equilibrium spectra (3.45)–(3.47) against non-dimensionalised wavenumber for two different values of $\beta /k_*^2 = P/(Qk_*^2)$ (3.42), the $k/k_*$ -independent ratio of buoyancy-induced and Doppler-shift-induced diffusivities. (a) Ratio $\beta /k_*^2 = 1$ corresponds to a diffusivity with equal contributions from buoyancy fluctuations and Doppler shift and (b) ratio $\beta /k_*^2 = 100$ corresponds to a buoyancy-fluctuation-dominated diffusivity. Diffusion by both effects is given in blue (3.45), whilst red and green lines are spectra derived from only the buoyancy and Doppler shift terms, respectively, given by (3.47) and (3.46).

The finite energy at $k \rightarrow \infty$ of (3.47) is unphysical. However, this solution is unstable in the sense that a vanishingly small Doppler shift contribution will result in $e(k) \rightarrow 0$ as $k\rightarrow \infty$ . This is because the limit of (3.45) as $k \rightarrow \infty$ , $\beta = o(k)$ is $e(k) \rightarrow 1/(4Qk^2)$ , i.e. the Doppler shift limit (3.46) for $k\gt k_*$ . This can be seen in figure 10(b). If a buoyancy fluctuation is present, so too is the Doppler shift term by the thermal wind balance (2.11) which means that the $k\gt k_*$ limit of (3.47) never occurs and the energy spectrum will always decay for large $k$ under Doppler shift.

Figure 10 shows how buoyancy fluctuations affect the spectrum amplitude, mainly for $k\lt k_*$ . For a Doppler-shift-dominated diffusivity ( $\beta \rightarrow 0$ ), this amplitude change is negligible. For $k \gt k_*$ , buoyancy fluctuations shallow the Doppler-shift-induced energy spectrum of KSV (2019) for a small range of intermediate wavenumbers.

The scaling argument of § 2.2 and the ratio of radial diffusivity (3.38) predict the buoyancy-induced diffusivity decays to zero for large $k/K_*$ and thus the effect of vertical buoyancy gradients is small. Our findings in this section compound this prediction because (i) a small buoyancy-induced diffusivity has negligible impact on the wave energy spectrum for any wavenumber and (ii) any buoyancy-induced diffusivity has negligible impact on the wave energy spectrum for large $k_h/K_*$ .

4. Diffusion induced by bottom topography

In this section, we evaluate the general diffusivity of § 3 for two systems. We extend the rotating shallow-water diffusivity in § 3.1 to include bottom topography (§ 4.1) and evaluate the general diffusivity for surface gravity waves scattered by bottom topography and a background current (§ 4.2). These diffusivities are evidence that the general diffusion framework can be applied to systems beyond those already discussed and we leave further exploration (ray tracing, comparison to observation etc.) to future work.

4.1. Rotating shallow-water system

We relax the assumption of a flat bottom made in §§ 2.1, 3.1 and 3.2 by writing

(4.1) \begin{equation} \Delta H = \Delta H_1 + \Delta H_2, \end{equation}

where $\Delta H_1$ is induced by geostrophy and satisfies (2.3) and $\Delta H_2$ is the fluctuation to mean height due to bottom topography. This is shown in figure 11.

Figure 11. Shallow-water set-up with bottom topography, a modified version of figure 1(b). The free surface is the solid blue line, the bottom is the solid black line and the horizontal dashed black lines enclose the constant mean height $\bar {H}$ . The wave perturbations are given by $h$ and the fluctuation in height due to geostrophy, the dashed-dotted line, is given by $\Delta H_1$ . The topographic variation in height is $\Delta H_2$ , which is a negative quantity at the labelled location in this figure.

As with the height fluctuation induced by geostrophy, we assume topography varies over longer length scales than the waves so that the WKB assumption holds and action is conserved (the WKB analysis of Appendix A does not rely on the flat-bottom assumption). The amplitude of topographic variation is small relative to the mean height so that the weak interaction assumption of statistical scattering is satisfied. The perturbation (2.5) to the bare frequency has an additional term induced by topography:

(4.2) \begin{equation} \omega _{1} = \underbrace {\frac {g \Delta H_1 k_h^2}{2 \omega _{0}}}_{\text{geostrophic height fluctuation}} + \underbrace {\frac {g \Delta H_2 k_h^2}{2 \omega _{0}}}_{\text{topography}} + \underbrace {\vphantom {\frac {g \Delta H k_h^2}{2 \omega _{0}}} \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}}_{\text{Doppler shift}}. \end{equation}

If the background flow and topography are uncorrelated, the diffusivity (3.26) has an additional term:

(4.3) \begin{equation} {{\unicode{x1D60B}}}_{\phi \phi } = \underbrace {\frac {4 \pi k_h^2}{c \vphantom {\omega ^2_0}} \int _{0}^{\infty }\!\! K_h^2 E(K_h) \, {\text {d}}K_h}_{\text{Doppler shift}} \, + \, \underbrace {\frac {\pi f^2 k_h^4}{c \omega _{0}^2} \int _{0}^{\infty }\!\! E(K_h) \, {\text {d}}K_h}_{\text{geostrophic height fluctuation}} + \, \underbrace {\frac {\pi g^2 k_h^4}{2 c \omega _{0}^2} \int _{0}^{\infty }\!\! K_h^2\mathcal {B}(K_h) \, {\text {d}}K_h}_{\text{topography}}, \end{equation}

where $\mathcal {B}$ is the bottom topography spectrum:

(4.4) \begin{equation} \mathcal {B}(\boldsymbol {K}) = \langle \widehat {\Delta H_2}(-\boldsymbol {K}) \widehat {\Delta H_2}(\boldsymbol {K})\rangle \end{equation}

which we assume to be horizontally isotropic such that $\mathcal {B}(\boldsymbol {K}) = \mathcal {B}(K_h)$ . The topography term in (4.3) is derived by replacing $E(K_h)$ in the height fluctuation term of (3.26) with $ (gK_h/f)^2 \mathcal {B}(K_h)/2$ , which comes from expressing $E(K_h)$ as a function of $\Delta H_1$ using geostrophic balance (2.3), then swapping $\Delta H_1$ for $\Delta H_2$ .

Should the background flow and topography be correlated, the diffusivity will contain geostrophic height fluctuation–topography and Doppler shift–topography cross terms.

4.2. Surface waves

Bretherton & Garrett (1968) prove that wave action is conserved for a wide class of physical systems. One of the examples they consider is surface gravity waves propagating on a mean flow with negligible rotation effects. In this case, the absolute frequency of the waves is

(4.5) \begin{equation} \omega = \left (g k_h \tanh {(H k_h )} \right )^{1/2} + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k} \end{equation}

with $\boldsymbol {U} = (U, V, 0)$ the background surface velocity. As in figure 11, $H$ is the wave-free height of the water column including both topographic variations $\Delta H_2$ and other variations in the mean height $\Delta H_1$ . Unlike in previous sections, we do not necessarily attribute these variations to geostrophy.

Providing the fluctuations to mean height $\bar {H}$ are small,

(4.6) \begin{equation} \omega _{1} \approx \frac {k_h(gk_h - \omega _{0}^2)\Delta H}{2\omega _{0}} + \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {k}. \end{equation}

For small-amplitude, horizontally isotropic height variations which vary over longer length and time scales than the waves, the azimuthal diffusivity is

(4.7) \begin{equation} {{\unicode{x1D60B}}}_{\phi \phi } = \underbrace {\frac {4 \pi k_h^2}{c \vphantom {\omega ^2_0}} \int _{0}^{\infty } K_h^2 E(K_h) \, {\text {d}}K_h}_{\text{Doppler shift}} \, + \, \underbrace {\frac {\pi k_h^2(gk_h - \omega _{0}^2)^2}{2 c \omega _{0}^2} \int _{0}^{\infty } K_h^2\mathcal {H}(K_h) \, {\text {d}}K_h}_{\text{height fluctuations}}, \end{equation}

which we find by comparing the surface wave dispersion relation (4.5) with that of shallow-water waves (4.2) and switching $g^2 k_h^4 \mathcal {B}(K_h)$ for $k_h^2(gk_h - \omega _{0}^2)^2\mathcal {H}(K_h)$ in the topography-induced diffusivity of the rotating shallow-water system (4.3). Here, $\mathcal {H}$ is the spectrum of fluctuations in total water depth. The diffusivity (4.7) is valid for a background velocity $\boldsymbol {U}$ not correlated with topography or other height-fluctuation effects.

The expressions given in this section are for diffusivities induced by an approximately time-independent $\omega _{1}$ which is clearly valid for bottom topography. Tolman (1990) considers surface waves propagating through temporally and spatially varying tidal currents which induce variations in depth. We refer the interested reader to CKV (2023) for work on the time-dependent case.

It should be possible to show that the topography-induced diffusivity in (4.7) is a limiting case of kinetic equations for Bragg scattering of surface gravity waves by bathymetry such as those derived by Hasselmann (1966) and Ardhuin & Herbers (2002). We leave this proof to future work.

5. Discussion

Scattering by turbulent flow leads to the diffusion of IGW energy in spectral space. In previous work on the topic (McComas & Bretherton 1977; Kafiabad et al. 2019; Savva 2020; Dong et al. 2020; Cox, Kafiabad & Vanneste 2023) the mechanism for spectral wave diffusion is a Doppler shift term in the waves’ dispersion relation. In this paper, we argue at the level of the dispersion relation that other mechanisms can cause significant wave diffusion. We provide two main examples: a fluctuation in the mean height of a rotating shallow-water system due to the background geostrophic flow and an effective variation in buoyancy frequency for the 3-D Boussinesq system, caused by vertical background flow buoyancy gradients. There is precedent for this – Chavanne et al. (2010) find refraction by vertical buoyancy gradients to be more significant than the Doppler shift effect in the context of internal tides.

We generalise the derivation of Kafiabad et al. (2019) to give a diffusion equation valid for any slowly evolving, weak inhomogeneities that the waves encounter. For our two examples, we evaluate the corresponding diffusivity.

In the Boussinesq case, we confirm the buoyancy-fluctuation effect can be significant. However, we find this effect is greatly reduced as waves diffuse to higher wavenumbers. We also find a reduction for higher-frequency waves and background flows with large aspect ratios between horizontal and vertical motions, as is the case in the ocean.

We solve the steady-state diffusion equation and find that the resulting energy spectrum agrees with this conclusion – as the waves move further into the WKB regime, the ratio between buoyancy fluctuation and Doppler shift effects decays to zero.

Scaling predictions from the Boussinesq dispersion relation of the relative importance of Doppler shift and vertical buoyancy gradients differ, with respect to wave frequency, from those found from the complete computation of the diffusivities. This is because large parts of the flow energy spectrum which do not meet the relevant resonance condition cannot contribute to scattering the waves.

In the shallow-water system, we find that for small Burger numbers the height-fluctuation diffusivity is comparable to or larger than the Doppler-shift diffusivity. This is supported by ray simulations which (i) validate the diffusion approximation of the action conservation equation and (ii) confirm the relative magnitude of the height fluctuation and Doppler shift effects. The Burger number does not have to be vanishingly small for the height fluctuation effect to be significant – at $Bu = O(1)$ , the two effects differ only by a factor of $1/4$ .

Remarkably, the Doppler shift effect is decoupled from the other diffusion mechanisms in the shallow-water and Boussinesq diffusion equations. This occurs despite the Doppler shift being related to these other effects through the geostrophic and thermal wind balances. Uncorrelated effects will also be uncoupled from the Doppler shift. For example, shallow-water waves and surface gravity waves scattered by topography, for which the corresponding diffusivities are evaluated in § 4.

Our work is relevant to internal tides which, alongside near-inertial waves, dominate the IGW energy spectrum (Ferrari & Wunsch 2009). Previous work uses ray tracing to model internal tide energy distributions and finds good agreement with observation despite marginal scale separation between the tides and background eddies (Park & Watts 2006; Rainville & Pinkel 2006; Chavanne et al. 2010). Furthermore, if internal tides propagate through a barotropic background flow, they form a set of uncoupled shallow-water equations (e.g. Savva & Vanneste 2018). Beyond this, it is clear that our approach applies for any waves propagating through large-scale inhomogeneities.

A possible application of the formula we obtain for the spectral diffusivity is to the representation of the impact of unresolved turbulence on IGWs parameterised in atmosphere and ocean models. In ray tracing modules such as MS-GWaM (Bölöni et al. 2021; Kim et al. 2021), the diffusion could be included by means of additional white-noise terms in the wavevector ray equation, leading to a stochastic parameterisation grounded in the physics of wave–flow interactions.

The diffusion regime is characterised by the weak-interaction assumption of statistical scattering and the WKB approximation. For diffusion induced by Doppler shift, the weak interaction assumption is equivalent to a weak-flow approximation: the group speed of the waves is much greater than the characteristic background flow speed, $c \gg U_*$ . Depending on the wave type, this may not hold consistently throughout the evolution. Conditions which ensure that $c \gg U_*$ for IGWs in the Boussinesq system are discussed in Appendix A of Cox et al. (2023). Dong, Buhler & Smith (2023) demonstrate the breakdown of the weak-flow assumption in this context.