Wave propagation in random two-dimensional turbulence: a multiscale approach

Abstract

To study two-dimensional dispersive waves propagating through turbulent flows, a new and less restrictive fast waves approximation is proposed using a multiscale setting. In this ansatz, large and small scales of the turbulence are treated differently. Correlation lengths of the random small-scale turbulence components can be considered negligible in the wave packet propagating frame. Nevertheless, the large-scale flow can be relatively strong, to significantly impact wavenumbers along the propagating rays. New theoretical results, numerical tools and proxies are derived to describe ray and wave action distributions. All model parameters can be calibrated robustly from the large-scale flow component only. We illustrate our purpose with ocean surface gravity waves propagating in different types of surface currents. The multiscale solution is demonstrated to efficiently document wave trapping effects by intense jets.

1. Introduction

This paper aims to revisit the ray path concept for fast waves propagating over heterogeneous turbulent flows. Considering ocean surface wave propagation, many authors have already discussed the random changes of rays subject to a random current (Voronovich 1991; White & Fornberg 1998; Smit & Janssen 2019), and consequences on wave action distributions. Closures have been derived in the Eulerian setting (Bal & Chou 2002; Klyatskin & Koshel 2015; Borcea, Garnier & Solna 2019; Kafiabad, Savva & Vanneste 2019; Bôas & Young 2020; Garnier, Gay & Savin 2020). Some of these approaches can be traced back to wave–wave interaction models (e.g. McComas & Bretherton (1977); see also Kafiabad et al. (2019), and references therein). In most cases, the central assumption is either time-delta-correlated turbulent velocity (Voronovich 1991; Klyatskin 2005; Klyatskin & Koshel 2015) and/or fast waves in comparison to fluid flow velocities (White & Fornberg 1998; Dysthe 2001; Bal & Chou 2002; Borcea et al. 2019; Kafiabad et al. 2019; Smit & Janssen 2019; Bôas & Young 2020; Garnier et al. 2020; Boury, Bühler & Shatah 2023; Wang et al. 2023). Medium variations may be slow, and delta-correlations are hardly justifiable in a fixed frame. However, attached to a fast-propagating wave group, the medium may seem to vary rapidly, and the delta-correlation assumption makes more sense. Another common assumption is frozen turbulence. In such a case, weak currents also imply conservation along rays of intrinsic frequency, wavenumber and group velocity magnitude in two dimensions (Boury et al. 2023). Subsequently, most wave dynamics models neglect variations and diffusion of frequency or wavenumber.

The diffusion of the wave action at large distance with a multiscale decomposition of the current has already been reported (Bal & Chou 2002). However, an explicit formulation for the diffusivity has been derived solely for a zero large-scale current. More generally, fast wave models rely mostly on either zero or constant current components at larger scales. West (1978), for instance, discussed acoustic waves in two-component random media, but no velocity was involved.

Hereafter, the proposed two-scale velocity decomposition falls into the family of stochastic transport models (Kunita 1997; Mikulevicius & Rozovskii 2004; Resseguier et al. 2020a; Zhen, Resseguier & Chapron 2023), including dynamics under location uncertainty (LU) (Mémin 2014; Resseguier, Mémin & Chapron 2017a) and stochastic advection by Lie transport (SALT) (Holm 2015). Under this framework, the small-scale velocity component is delta-correlated in time (Cotter, Gottwald & Holm 2017). Up to usual source terms, fluid dynamics quantities (temperature, momentum, etc.) are transported by both the large-scale revolved component and that random unresolved turbulence component. The stochastic closures obtained are conservative. Nonlinear wave Hamiltonian dynamics and wave influence on currents (e.g. Stokes drift) have then been derived (e.g. Crisan & Holm 2018; Bauer et al. 2020; Holm 2021; Holm & Luesink 2021; Dinvay & Mémin 2022; Holm, Hu & Street 2023). Considering a single-wavevector current, solutions for a monochromatic shallow-water wave were developed by Mémin et al. (2022). In the present study, our objective is restricted to the influence of turbulent flows on linear waves.

After first recalling the principles of the ray tracing method, we present the multiscale framework for fast wave dynamics, its physical grounds, and a calibration method for the closure. Simplified stochastic equations are then derived for the ray dynamics and the wave action spectrum, in both Lagrangian and Eulerian settings. For illustrative examples, numerical tools, analytic models and proxies are applied to ocean surface gravity waves propagating through two types of two-dimensional (2-D) turbulent flows: a typical slow homogeneous turbulence, and a jet case.

2. Characteristics of wave packet rays

Isolating a single progressive group of quasi-regular wave trains, it follows a form $h(\boldsymbol {x},t) \exp ({\mathrm {i} \phi (\boldsymbol {x},t)})+{\rm c.c.}$, for most properties. Typically, $h$ would be the local wave height, in metres. If a packet is to be followed, then the phase $\phi (\boldsymbol {x},t)$ must vary smoothly along the propagation, i.e. $\phi (\boldsymbol {x},t)$ is differentiable. The relative frequency is then ${\omega = -\partial _t \phi (\boldsymbol {x},t)}$, and the wavenumber vector is $\boldsymbol {k} = \boldsymbol {\nabla } \phi (\boldsymbol {x},t)$, with wavenumber $k=\| \boldsymbol {k}\|$, and direction given by the normalized wavevector, $\tilde {\boldsymbol {k}} = \boldsymbol {k} / k = \left (\begin{smallmatrix} \cos \theta _k \\ \sin \theta _k \end{smallmatrix}\right )$. To first order, such a train of waves is dispersive, and the intrinsic frequency reads

(2.1) \begin{equation} \omega - \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{v}=\omega_0 = \begin{cases} \mathrm{constant}\times\dfrac{1}{\alpha}\,k^\alpha , & \alpha \neq0, \\ \mathrm{constant}\times\log(k), & \alpha = 0, \end{cases}\end{equation}

and propagates with its group velocity $\boldsymbol v_g = \boldsymbol {\boldsymbol {\nabla }}_{\boldsymbol {k}} \omega$, constantly modified by the local velocity of the currents $\boldsymbol v$,

(2.2)\begin{equation} \frac{{\mathrm{d}} \boldsymbol{x}_r}{{\mathrm{d}} t} =\boldsymbol v_g=\boldsymbol v_g^0+ \boldsymbol v, \end{equation}

where $\boldsymbol {x}_r$ is the centroid of a wave group, $\boldsymbol v_g^0 = ({\partial \omega _0 (k) }/{\partial k})\,\tilde {\boldsymbol {k}}$ is the group velocity without currents, i.e. depending solely on the wavevector. For $\alpha =1$, the medium is non-dispersive (e.g. acoustic waves). Parameter $\alpha =1/2$ corresponds to gravity waves over deep ocean ($\omega _0=\sqrt {gk}$). The dominant wavevector $\boldsymbol {k}$ within the group evolves according to

(2.3)\begin{equation} \frac{{\mathrm{d}} \boldsymbol{k}}{{\mathrm{d}} t}={-}\boldsymbol{\boldsymbol{\nabla}} {\boldsymbol v}^{{\rm T}}\boldsymbol{k}. \end{equation}

Equations (2.2)–(2.3) are the Hamilton eikonal equations. Along the propagating ray, velocity gradients induce linear variations. Decelerating currents will, for instance, shorten waves, and reduce the group velocity. Travelling over fields of random velocities $\boldsymbol v$, the wavevector $\boldsymbol {k}$ will also become randomly distributed. Scattering of ocean surface wave packets by random currents can generally be assumed to be weak, with $\| \boldsymbol v \|$ of order $0.5$ m s$^{-1}$, much smaller than $v_g^0 = \| \boldsymbol {v}_g^0 \|$ of order $10$ m s$^{-1}$. Yet cumulative effects of these random surface currents can lead to strong convergence or divergence between initially nearby ray trajectories.

To complete the wave field description, $E(\boldsymbol {x},t)= \frac 12 \rho g h^2(\boldsymbol {x},t)$ and $A(\boldsymbol {x},t)= E(\boldsymbol {x},t)/ \omega _0(k(\boldsymbol {x},t))$ denote energy and action by unit of surface. Here, $E$ is expressed in J m$^{-2}$, and $A$ in J s m$^{-2}$. To avoid spurious notations, we set the multiplicative constant $\frac 12 \rho g$ to unity. The wave action is considered to be an adiabatic invariant in the absence of source terms. Wave action is then crucial to anticipate wave transformations by currents (White 1999). Unlike wave energy, wave action is conserved, in the absence of wave generation or dissipation. This action is the integral over wavevectors of the action spectrum $N$, also related to the wave energy spectrum $E$:

(2.4)\begin{equation} A(\boldsymbol{x},t)=\int {\mathrm{d}} \boldsymbol{k} \, N(\boldsymbol{x},\boldsymbol{k},t)= \int {\mathrm{d}} \boldsymbol{k}\,\frac{E(\boldsymbol{x},\boldsymbol{k},t)}{\omega_0(\boldsymbol{k},t)}.\end{equation}

Action and energy spectrum quantify action and energy by unit of surface (unit of $\boldsymbol {x}$) and by unit of wavevector surface (unit of $\boldsymbol {k}$). Consider the $(\boldsymbol {x},\boldsymbol {k})$ variable change between different times $t_i$ and $t_f$ integrating the characteristic eikonal equations (2.2)–(2.3):

(2.5)\begin{equation} \begin{pmatrix} \boldsymbol{x}_r(t_i) \\ \boldsymbol{k}(t_i) \end{pmatrix}\mapsto \begin{pmatrix} \boldsymbol{x}_r(t_f) \\ \boldsymbol{k}(t_f) \end{pmatrix}. \end{equation}

According to the Liouville theorem for Hamiltonian mechanics (Landau & Lifshitz 1960, § 46), the state space of the ‘packet-by-packet’ approach (the $(\boldsymbol {x},\boldsymbol {k})$ space) does not contract or dilate along time. Readers not familiar with Hamiltonian dynamics may see the divergence free of the four-dimensional flow (2.5) – i.e. $\boldsymbol {\nabla }_{\boldsymbol {x}} \boldsymbol {\cdot } {{\mathrm {d}} \boldsymbol {x}_r}/{{\mathrm {d}} t} + \boldsymbol {\nabla }_{\boldsymbol {k}} \boldsymbol {\cdot } {{\mathrm {d}} \boldsymbol {k}}/{{\mathrm {d}} t} =0$ – as the divergence free of incompressible flow velocities, leading naturally to volume-preserving dynamics. Therefore, if wave dissipation is neglected, then the wave action spectrum $N$ is conserved (Lavrenov 2013), i.e.

(2.6)\begin{equation} N (\boldsymbol{x}_r(t_i),\boldsymbol{k}(t_i),t_i)=N (\boldsymbol{x}_r(t_f),\boldsymbol{k}(t_f),t_f). \end{equation}

This result is extremely useful because it involves only quantities of the characteristics, i.e. each Fourier mode can be modified independently of the others. The wave energy spectrum can be computed from the characteristics

(2.7)\begin{equation} E (\boldsymbol{x}_r(t_f),\boldsymbol{k}(t_f),t_f)=\frac{\omega_0 (\boldsymbol{k}(t_f))}{\omega_0 (\boldsymbol{k}(t_i))}\, E (\boldsymbol{x}_r(t_i),\boldsymbol{k}(t_i),t_i). \end{equation}

starting with an initial incoming wave spectrum $E (\boldsymbol {x}_r(t_i),\boldsymbol {k}(t_i),t_i)$ for every wavevector $\boldsymbol {k}(t_i)$, starting from a small set of spatial points $\boldsymbol {x}_r(t_i)$.

3. A new fast wave assumption

Eikonal equations (2.2)–(2.3) are driven by currents and their gradients. Commonly, the Eulerian current $\boldsymbol {v}$ is decomposed into a low-frequency large-scale component $\bar {\boldsymbol {v}}$ and a transient small-scale unresolved component $\boldsymbol {v}'$:

(3.1)\begin{equation} {\boldsymbol{v}} = \bar{\boldsymbol{v}} + \boldsymbol{v}'. \end{equation}

Current gradients naturally follow the same scale separation. From now on, we will consider divergence-free 2-D currents only.

3.1. The ray Lagrangian correlation time

To better characterize the wave dynamics in such a random environment, the covariance of the fluid velocity can be evaluated in the wave group frame. To take into account the small-scale unresolved component $\boldsymbol {v}'$, its Eulerian spatio-temporal covariance is considered, assuming statistical homogeneity and stationarity for the Eulerian velocity $\boldsymbol {v}'_E(t,\boldsymbol {x})=\boldsymbol {v}'(t,\boldsymbol {x})$:

(3.2)\begin{align} C_{ij}^{v'_{E}}(\delta t,\delta \boldsymbol{x}) =\mathbb{E} (v'_i (t,\boldsymbol{x})\,v'_j (t+\delta t,\boldsymbol{x}+\delta\boldsymbol{x})) =\mathbb{E} (v'_i (t,\boldsymbol{x}_r(t))\,v'_j (t+\delta t,\boldsymbol{x}_r(t)+\delta\boldsymbol{x})), \end{align}

where $\boldsymbol {x}_r$ is a solution of (2.2) with an arbitrary initial position $\boldsymbol {x}_r^0$. Then we define ${v'_R(t) =v' (t,\boldsymbol {x}_r(t))}$, the Lagrangian velocity along the ray $\boldsymbol {x}_r(t)$. The temporal covariance of the small-scale component $\boldsymbol {v}'$ – in the wave group frame – is the covariance of that Lagrangian velocity:

(3.3)\begin{equation} C_{ij}^{v'_R} (\delta t)= \mathbb{E} (v'_i (t,\boldsymbol{x}_r(t))\,v'_j (t+\delta t,\boldsymbol{x}_r(t+\delta t))) =C_{ij}^{v'_E} (\delta t, \boldsymbol{x}_r(t+\delta t)-\boldsymbol{x}_r(t)). \end{equation}

Assuming, for example, a typical isotropic form for the Eulerian covariance,

(3.4)\begin{equation} C^{v'_E} ( \delta t, \delta \boldsymbol{x} ) = C \left(\frac{|\delta t|}{\tau_{v'}} + \frac{\| \delta \boldsymbol{x} \|}{l_{v'}}\right), \end{equation}

the covariance can be evaluated in the wave group frame for small time increment $\delta t$:

(3.5)\begin{align} C^{v'_R} ( \delta t )= C \left(\frac{|\delta t|}{\tau_{v'}} + \frac{\| \boldsymbol{x}_r(t'+t)-\boldsymbol{x}_r(t') \|}{l_{v'}}\right) =C \left(\left( \frac{1}{\tau_{v'}} + \frac{\| \boldsymbol{v}_g \|}{l_{v'}}\right) |\delta t| + O(\delta t^2) \right), \end{align}

since $\boldsymbol {x}_r(t'+t)-\boldsymbol {x}_r(t') = \boldsymbol {v}_g\,\delta t + O(\delta t^2)$. Therefore, $({1}/{\tau _{v'}} + {\| \boldsymbol {v}_g \|}/{l_{v'}})^{-1}$ is the correlation time of $\boldsymbol {v}' (t,\boldsymbol {x}_r(t))$. For fast waves, the along-ray correlation time of the small-scale velocity can be approximated by ${l_{v'}}/{v_g^0 }$. Note that eikonal equations (2.2)–(2.3) involve both velocity and velocity gradients. The above derivation is also valid for the small-scale velocity gradients $(\boldsymbol {\nabla } {\boldsymbol {v}}^{{\rm T}})'(t,\boldsymbol {x}_r(t))$. The ratio $\epsilon$ between that along-ray correlation time and the characteristic time of the wave group properties evolution will then control the time decorrelation assumption of $\boldsymbol {v}'$:

(3.6)\begin{equation} \epsilon =\frac{l_{v'}}{v_g^0 }\,\|\boldsymbol{\nabla} {\boldsymbol{v}}^{{\rm T}}\| \sim\frac{l_{v'}}{l_{v}}\,\frac{\|\boldsymbol{v}\|}{v_g^0}.\end{equation}

This time scale estimation can be obtained from spatio-temporal covariances more general than (3.4) (not shown) even though the derivation is more technical. Note that the Eulerian small-scale velocity $\boldsymbol {v}'$ is not necessarily time-uncorrelated, as assumed in Voronovich (1991) and Klyatskin & Koshel (2015). Yet for small enough $\epsilon$, the Lagrangian small-scale velocity along the ray can be considered time-uncorrelated. From the expression for $\epsilon$, such a condition depends upon:

1. (i) $v_g^0$, the fast wave group velocity;

2. (ii) $\|\boldsymbol {v}\|$, often slow but not always negligible compared to the intrinsic wave group $v_g^0$;

3. (iii) ${l_{v'}}/{l_{v}}$, related to the separation between large scales $\bar {\boldsymbol {v}}$ and small scales $\boldsymbol {v}'$, e.g. the spatial filtering cutoff of the large-scale velocity $\bar {\boldsymbol {v}}$, but also related to its kinetic energy distribution over spatial scales, typically the spectrum slope.

This along-ray partial time-decorrelation assumption is less restrictive than the usual fast wave approximation (White & Fornberg 1998; Dysthe 2001; Bal & Chou 2002; Borcea et al. 2019; Kafiabad et al. 2019; Smit & Janssen 2019; Bôas & Young 2020; Garnier et al. 2020; Boury et al. 2023; Wang et al. 2023) – say ${\|\boldsymbol {v}\|}/{v_g^0 } \ll 1$ – and than the SALT-LU time-decorrelation used for turbulence dynamics (Mémin 2014; Holm 2015; Cotter et al. 2017; Resseguier et al. 2020a) – say ${l_{v'}}/{l_{v}} \ll 1$. Similarly, this last validity criterion can be obtained by replacing $\boldsymbol {x}_r$ in (3.2)–(3.6) by the fluid particle Lagrangian path $\boldsymbol {x}$ (solution of ${{\mathrm {d}} \boldsymbol {x}}/{{\mathrm {d}} t} =\boldsymbol {v}$) and thus $\boldsymbol {v}_g^0$ by $\boldsymbol {v}$. These asymptotic models often rely on averaging or homogenization techniques (Papanicolaou & Kohler 1974; White & Fornberg 1998) to derive Markovian dynamics involving various types of diffusivity.

3.2. Ray absolute diffusivity and turbulence statistics: calibration

Diffusivity is a natural tool to specify statistics of uncorrelated random media. For waves in random media, we will specify multi-point statistics, and the Fourier space is convenient for this purpose. We will first present scalar diffusivity and then distribute it over spatial scales to fully calibrate the random velocity $\boldsymbol {v}'$, i.e. choose some parameter values to set the statistics of that velocity field. As such, we will obtain a closed model to derive analytic results and generate samples for simulations.

The absolute diffusivity (or Kubo-type formula) usually corresponds, in the so-called diffusive regime, to the variance per unit of time of a fluid particle Lagrangian path ${{\mathrm {d}} \boldsymbol {x}(t)}/{{\mathrm {d}} t} =\boldsymbol {v}_L(t)= \boldsymbol {v} (t,\boldsymbol {x}(t))$. It is approximately equal to the velocity variance times its correlation time. The Eulerian velocity covariance (3.4) will thus induce an absolute diffusivity (Piterbarg & Ostrovskii 1997; Klyatskin 2005)

(3.7)\begin{equation} {\tfrac 12}\,a^L=\int_0^{+\infty} {\mathrm{d}} \delta t\, C^{v'_L} (\delta t)=\int_0^{+\infty} {\mathrm{d}}\delta t\, C^{v'_E} (\delta t, \boldsymbol{x}(t+\delta t)-\boldsymbol{x}(t))\approx \tfrac 12\,\tau_{v'}C(0). \end{equation}

This diffusivity well describes effects of fast-varying eddies, but is not appropriate in our case. Indeed, along a propagating wave group ${{\mathrm {d}} \boldsymbol {x}_r(t)}/{{\mathrm {d}} t} = \boldsymbol v_g^0(t) + \boldsymbol v_R(t)$, a ray absolute diffusivity occurs and slightly differs from the usual absolute diffusivity to become

(3.8)\begin{equation} {\frac 12}\,a^{R}=\int_0^{+\infty} {\mathrm{d}} \delta t\, C^{v'_R} (\delta t)\approx\frac 12\left(\frac{1}{\tau_{v'}} + \frac{\| \boldsymbol{v}_g \|}{l_{v'}}\right)^{{-}1} C(0)\approx\frac 12\, \frac{l_{v'}}{v_g^0}\,C(0). \end{equation}

The absolute diffusivity sets the amplitude of the small-scale velocity $\boldsymbol {v}'$. Indeed, since the kinetic energy of a time-continuous white noise is infinite, it has no physical meaning. It is more relevant to deal with absolute diffusivity rather than kinetic energy in order to describe the statistics of the time-uncorrelated velocity. To calibrate its spatial correlations, we may focus on its Fourier transform $\widehat {\boldsymbol {v}'}(\boldsymbol {\kappa },t)$, denoting by ${\boldsymbol {\kappa } = \kappa \left (\!\begin{smallmatrix} \cos \theta _\kappa \\ \sin \theta _\kappa \end{smallmatrix}\!\right )}$, the surface current wavevector. By analogy with the current kinetic energy spectra $E_\kappa = \frac 12 \oint _0^{2{\rm \pi} } {\mathrm {d}} \theta _\kappa \,\kappa ({\|\hat {\boldsymbol {v} }(\boldsymbol {\kappa },t) \|^2}/{(2{\rm \pi} )^2})$, Resseguier, Mémin & Chapron (2017b) and Resseguier, Pan & Fox-Kemper (2020b) decompose the absolute diffusivity scale by scale:

(3.9)\begin{equation} a^R=\int_0^{+\infty} A^{R} _{v'} (\kappa) \,{\mathrm{d}} \kappa. \end{equation}

Referring it to absolute diffusivity spectral density (ADSD), it is defined as the kinetic energy spectra multiplied by the correlation time at each scale, $\tau (\kappa )$. Unlike Resseguier et al. (2017b, 2020b), that correlation time is here imposed by the wave dynamics. Therefore, by analogy with (3.8), we choose a correlation time $\tau ^R(\kappa )={1/ \kappa }/{v_g^0 (k)}$, and then

(3.10)\begin{equation} {\frac 12}\,A^{R} (\kappa)=\frac 12\,\tau^R(\kappa)\,E_\kappa (\kappa) =\frac 12\,\frac {1/ \kappa} { v_g^0 (k) }\,E_\kappa (\kappa), \end{equation}

where $k$ denotes the wave wavenumber and $\kappa$ the current wavenumber.

To calibrate an equivalent noise, we model $\boldsymbol {v}'$ by $\boldsymbol {\sigma } \,{{\mathrm {d}}} B_t/{\mathrm {d}} t$, where ${{\mathrm {d}}} B_t/{\mathrm {d}} t$ is a spatio-temporal white noise, and $\boldsymbol {\sigma }$ denotes a spatial filtering operator that encodes spatial correlations through its ADSD, $A^{R} _{v'}$ and the horizontal incompressibility condition ($\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\sigma }=0$). For incompressibility, we work with the curl of a streamfunction. To generate a homogeneous and isotropic streamfunction, we can filter a one-dimensional white noise $\dot {B}$ with a filter $\breve {\psi }_\sigma$ (Resseguier et al. 2017b), i.e. $\breve {\psi }_\sigma \star \dot {B}$, where $\star$ denotes a spatial convolution. The velocity field is hence

(3.11)\begin{equation} \boldsymbol{v}' = \boldsymbol{\sigma} \,{{\mathrm{d}}} B_t/{\mathrm{d}} t= \boldsymbol{\nabla}^{\bot} \breve{\psi}_{\sigma} \star{{\mathrm{d}}} B_t/{\mathrm{d}} t, \end{equation}

with $\boldsymbol {\nabla }^{\bot }$ the 2-D curl. That formula is easily written and implementable in Fourier space (see (A2)). To define the streamfunction filter, we note that $({{\rm \pi} \kappa ^3}/{(2{\rm \pi} )^2}) | \widehat { \breve {\psi }_{\sigma } }( \kappa ) |^2 = \frac 12 \oint _0^{2{\rm \pi} } {\mathrm {d}} \theta _\kappa \,\kappa ({\| \widehat { \boldsymbol {\sigma } \,{{\mathrm {d}}} B_t }(\boldsymbol {\kappa }) \|^2}/{(2{\rm \pi} )^2\,{\mathrm {d}} t}) = A^{R} _{v'} (\kappa )$, i.e. the filter can be fully defined by the small-scale ADSD $A^{R} _{v'}$. To close our model, we assume an ADSD power law:

(3.12)\begin{equation} A^{R} (\kappa)\approx A^R_0 \kappa^{-\mu}. \end{equation}

It enables automatic closure calibration $A^{R} _{v'} ( \kappa ) = A^R_0 \kappa ^{-\mu } - A^{R}_{\bar {v}} ( \kappa )$, from instantaneous large-scale current statistics $A^{R} _{\bar {v}}$ only (Resseguier et al. 2020b), as illustrated in figure 1.

Figure 1. (a) Kinetic energy (KE) spectrum (m$^2$ s$^{-2}$/ (rad m$^{-1}$)) and (b) ADSD (m$^2$ s$^{-1}$/ (rad m$^{-1}$)) of the resolved high-resolution velocity $A^{R}$ in red, low-resolution velocity $A^{R}_{\bar {v}}$ in blue, and modelled stochastic velocity $A^{R} _{v'} ( \kappa ) = A^R_0 \kappa ^{-\mu } - A^{R}_{\bar {v}} (\kappa )$ in green. For the ADSD power law $A^{R} (\kappa )\approx A^R_0 \kappa ^{-\mu }$, we impose the theoretical kinetic energy spectrum slope $- \frac 5 3$ (black solid line), coherently with homogeneous surface quasi-geostrophic dynamics (see § 5). The residual ADSD (green line) is set to extrapolate that power law at small scales.

4. Statistical wave dynamics

In a stochastic framework, the Stratonovich or Itō notations can both be used (Kunita 1997; Oksendal 1998). Under Stratonovich calculus rules, expressions become similar to deterministic ones. Specifically, stochastic versions of linearized dynamical equations are obtained by replacing $\boldsymbol {v}$ by $\bar {\boldsymbol {v}} + \boldsymbol {\sigma } \circ {{\mathrm {d}}} B_t/{\mathrm {d}} t$. Then the stochastic transport of phase $({{\mathrm {d}}}\phi /{{\mathrm {d}} t}) =\omega _0(\| \boldsymbol {\nabla } \phi \|)$, i.e. – up to that velocity replacement – the Stratonovich dispersion relation, is exactly (2.1). The method of characteristics also applies. Note that one can switch from Stratonovich to Itō notations, where $\boldsymbol {v}'$ corresponds to $\boldsymbol {\sigma } \,{{\mathrm {d}}} B_t/{\mathrm {d}} t$. The characteristics equations (2.2)–(2.3) also remain unchanged for homogeneous and isotropic $\boldsymbol {v}'$:

(4.1)\begin{equation} \begin{cases} {\mathrm{d}} \boldsymbol{x}_r =(\boldsymbol v_g^0+ \bar{\boldsymbol v})\,{\mathrm{d}} t+\boldsymbol{\sigma} \,{{\mathrm{d}}} B_t,\\ {\mathrm{d}} \boldsymbol{k} ={-} \boldsymbol{\nabla} (\bar{\boldsymbol{v}}\,{\mathrm{d}} t+\boldsymbol{\sigma} \,{{\mathrm{d}}} B_t)^{{\rm T}}\boldsymbol{k}. \end{cases} \end{equation}

4.1. Single-ray stochastic differential equations

When studying a single ray in a homogeneous and isotropic turbulence (3.11), the wavevector dynamics simplifies. In the local crest-oriented frame, the influence of small-scale currents can be represented solely by four one-dimensional white noise forcings.

Notably, dynamics of wavevectors (2.3) are similar to tracer gradient dynamics (Bühler 2009; Plougonven & Zhang 2014). Only the coupled ray path dynamics (2.2) differs. Accordingly, we follow the notations and derivations of the mixing analysis from Lapeyre, Klein & Hua (1999), and references therein. Without loss of generality, the large-scale velocity can be parametrized as

(4.2a,b)\begin{equation} \bar{\boldsymbol{v}}=\bar{v} \begin{pmatrix} \cos\bar{\theta} \\ \sin\bar{\theta} \end{pmatrix}\quad \text{and} \quad \boldsymbol{\nabla}{\bar{\boldsymbol{v}}}^{{\rm T}}= \frac{1}{2} \begin{bmatrix} \bar{\sigma} \sin{2 \bar{\phi}} & \bar{\omega} + \bar{\sigma} \cos{2 \bar{\phi}} \\ -\bar{\omega} + \bar{\sigma}\cos{2 \bar{\phi}} & - \bar{\sigma} \sin{2 \bar{\phi}} \end{bmatrix}. \end{equation}

Figure 2 provides a synthetic view of angles involved. The dynamics wave group centroid $\boldsymbol {x}_r$ is driven directly by the large current wave group velocity $\boldsymbol v_g^0 + \bar {\boldsymbol v}$. The influence of the large-scalecurrent gradients on the wavevector dynamics (4.1), expressed in the local crest-oriented frame $(\tilde {\boldsymbol {k}}, \tilde {\boldsymbol {k}}^\bot )$, is straightforward (Lapeyre et al. 1999). The small-scale currents force the ray dynamics through a stochastic noise. For a single ray $(\boldsymbol {x}_r,\boldsymbol {k})=(x_r, y_r,k \cos \theta _k, k \sin \theta _k)$, this noise can be described rigorously by four independent one-dimensional white noises only (see Appendix A), $\dot {B}_t^{(1)}$, $\dot {B}_t^{(2)}$, $\dot {B}_t^{(3)}$ and $\dot {B}_t^{(4)}$, and

(4.3)$$\begin{gather} \frac{{\mathrm{d}} }{{\mathrm{d}} t} x_{r}= v_g^0\cos {\theta}_k + \bar{v} \cos \bar{\theta} + \sqrt{a_0}\,{\dot{B}_t^{(1)}}, \end{gather}$$

(4.4)$$\begin{gather}\frac{{\mathrm{d}}}{{\mathrm{d}} t} y_{r} = v_g^0\sin {\theta}_k + \bar{v} \sin \bar{\theta} + \sqrt{a_0}\,{\dot{B}_t^{(2)}}, \end{gather}$$

(4.5)$$\begin{gather}\frac{{\mathrm{d}} }{{\mathrm{d}} t} \log k ={-}\bar{\sigma} \sin(\zeta) + \gamma_0 + \sqrt{\gamma_0}\,{\dot{B}_t^{(3)}}, \end{gather}$$

(4.6)$$\begin{gather}\frac{{\mathrm{d}} }{{\mathrm{d}} t} {\theta}_{k} = \frac{1}{2}\,(\bar{\omega} - \bar{\sigma} \cos(\zeta) ) + \sqrt{3 \gamma_0}\,{ \dot{B}_t^{(4)}}, \end{gather}$$

where $\zeta = {2({\theta }_{k}+\bar {\phi })}$ and

(4.7)$$\begin{gather} a_0=\frac 1{2\,{\mathrm{d}} t}\,\mathbb{E}\| \boldsymbol{\sigma}\,{{\mathrm{d}}} B_t \|^2= \int_0^{+\infty} A^{R} _{v'} (\kappa) \,{\mathrm{d}} \kappa, \end{gather}$$

(4.8)$$\begin{gather}\gamma_0 = \frac 1{8\,{\mathrm{d}} t}\,\mathbb{E} \| \boldsymbol{\nabla}_{\boldsymbol{x}} (\boldsymbol{\sigma}\,{{\mathrm{d}}} B_t)^{{\rm T}} \|^2 = \frac 14 \int_0^{+\infty} k^2\,A^{R} _{v'} (\kappa) \,{\mathrm{d}} \kappa. \end{gather}$$

Diffusivity constants depend through (3.10) on both the correlation length and the spectrum slope of the small-scale velocity. In contrast to the classical fast wave approximation, the wavenumber does vary. This is due to (i) the finite large-scale strain rate $\bar {\sigma }$, and (ii) the small-scale isotropic velocity model (3.11). This isotropy assumption and its implication are discussed in Appendix C. Note that neither the large-scale component nor the small-scale component is assumed to be steady, even though that Eulerian velocity unsteadiness is only a secondary process in the wave dynamics. The fast temporal variations seen by the wave are driven mainly by the large wave speed and not by the Eulerian velocity unsteadiness. The current unsteadiness can also lead to wavenumber variations (Dong, Bühler & Smith 2020; Boury et al. 2023; Cox, Kafiabad & Vanneste 2023). Given a known wavevector angle, it leads to a wavenumber evolution

(4.9)\begin{equation} k(t) = k(0)\exp \left(- \int_0^t \bar{\sigma} \sin({2({\theta}_{k}+\bar{\phi})}) \,{\rm d} t'\right) \exp (\gamma_0 t + \sqrt{\gamma_0}\,B_t^{(3)}), \end{equation}

and hence to the complete wavevector distribution, i.e. the wave spectrum. The second exponential factor in (4.9) is a geometric Brownian motion. Its mean diverges in time exponentially rapidly. Physically, shear and strain of $\boldsymbol {v}'$ tends to shorten the wavelength (Voronovich 1991; Boury et al. 2023) leading to this exponential divergence. This factor has a log-normal distribution, suggesting possible extreme transient wavenumber events. This generalizes previous results (Voronovich 1991; Klyatskin & Koshel 2015), obtained with neglecting the time-correlated current component $\bar {\boldsymbol {v}}$.

Figure 2. Schematic view of vectors and angles involved in single-ray dynamics, where $\bar {\boldsymbol {S}}_-$ and $\bar {\boldsymbol {S}}_+$ are respectively compression and dilatation axes associated with the large-scale velocity gradient $\boldsymbol {\nabla }{\bar {\boldsymbol {v}}}^{{\rm T}}$.

For completeness, the action distribution over space and wavevector can be derived. Some approaches consider finite-size wave trains either through additional equations (Jonsson 1990; White & Fornberg 1998) or re-meshing (Hell, Fox-Kemper & Chapron submitted). Otherwise, each ray transports its action spectrum (2.6), and we need to numerically combine many rays (Lavrenov 2013), or rely on analytic approximations. Typically, we solve (4.3)–(4.5), exhibiting $p(\boldsymbol {x},\boldsymbol {k} \,|\, \boldsymbol {x}_r^0,\boldsymbol {k}_r^0,t)$, the distribution of the ray $(\boldsymbol {x},\boldsymbol {k})$ at time $t$ given initial conditions $(\boldsymbol {x}_r^0,\boldsymbol {k}^0)$. Then by analogy with tracers in incompressible turbulence (Piterbarg & Ostrovskii 1997, (1.31); see also Appendix D), we can evaluate the wave action spectrum mean – or any pointwise statistics – as

(4.10)\begin{equation} \mathbb{E} N(\boldsymbol{x},\boldsymbol{k},t)= \iint {\mathrm{d}} \boldsymbol{x}_r^0 \,{\mathrm{d}} \boldsymbol{k}^0\,N^0(\boldsymbol{x}_r^0,\boldsymbol{k}_r^0)\,p(\boldsymbol{x},\boldsymbol{k} \,|\, \boldsymbol{x}_r^0,\boldsymbol{k}^0,t), \end{equation}

where $N^0$ is the initial wave action spectrum. Integrating this expression over wavevectors, we note that the distribution inside the integrals changes:

(4.11)\begin{equation} \mathbb{E} A(\boldsymbol{x},t)= \iint {\mathrm{d}} \boldsymbol{x}_r^0 \,{\mathrm{d}} \boldsymbol{k}^0\, N^0(\boldsymbol{x}_r^0,\boldsymbol{k}_r^0)\,p(\boldsymbol{x} \,|\, \boldsymbol{x}_r^0,\boldsymbol{k}^0,t). \end{equation}

The wave action mean depends solely on group positions distribution. Multi-point action statistics – e.g. focusing $\mathbb {E}\|\boldsymbol{\nabla} _x A\|^2$ – rely on multi-ray correlations, encoded in the stochastic characteristic equations (4.1), but not the simplified model (4.3)–(4.6). Alternatively, Eulerian descriptions of wave action dynamics directly provide action distribution over space and wavevector.

4.2. Eulerian dynamics and action diffusion

Wave action spectrum is transported along a four-dimensional volume-preserving stochastic flow (4.1). Again by analogy with incompressible turbulence (Resseguier et al. 2017a), the stochastic transport of wave action spectrum in Itō notations reads

(4.12)\begin{align} & \partial_t N + \left(\boldsymbol{v}_g^0 + \bar{\boldsymbol{v}} + \boldsymbol{\sigma}\,\frac{{{\mathrm{d}}} B_t}{{\mathrm{d}} t}\right) \boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}} N + \left(-\boldsymbol{\nabla}_{\boldsymbol{x}} \left(\bar{\boldsymbol{v}} + \boldsymbol{\sigma}\,\frac{{{\mathrm{d}}} B_t}{{\mathrm{d}} t}\right)^{{\rm T}} \boldsymbol{k}\right) \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{k}} N \nonumber\\ &\quad =\begin{bmatrix} \boldsymbol{\nabla}_{\boldsymbol{x}} \\ \boldsymbol{\nabla}_{\boldsymbol{k}} \end{bmatrix}\boldsymbol{\cdot} \left(\boldsymbol{\mathsf{D}} \begin{bmatrix} \boldsymbol{\nabla}_{\boldsymbol{x}} \\ \boldsymbol{\nabla}_{\boldsymbol{k}} \end{bmatrix}N\right) =\frac 1 2\,a_0\,\Delta_{\boldsymbol{x}}N+\frac 1 2\,\gamma_0\, \frac 1 {k}\,\partial_{k}\left(k^3\,\partial_{k}N\right) +\frac 3 2\,\gamma_0\,\partial^2_{\theta_k}N. \end{align}

The right-hand side is reminiscent of (3.16) in Bôas & Young (2020), and (36) in Smit & Janssen (2019), and more generally of rapid wave models. Nevertheless, (4.12) is not averaged and explicitly involves large-scale currents and noise terms (terms with factor ${{{\mathrm {d}}} B_t}/{{\mathrm {d}} t}$). Differences with Smit & Janssen (2019) and Bôas & Young (2020) for the diffusivity estimates and the detailed computation of the $4\times 4$ diffusion matrix $\boldsymbol{\mathsf{D}}$ can be found in Appendix A. Itō notations for (4.12) explicitly separate mean terms (e.g. diffusion terms) and zero-mean noise terms. Here, the Eulerian Itō notations reveal that coefficients $\frac 1 2 a_0$, $\frac 1 2 \gamma _0$ and $\frac 32 \gamma _0$ act to diffuse wave action in space, wavenumber and wavevector angle, respectively.

5. Numerical experiments

To illustrate these developments, we consider ocean surface gravity waves propagating over a dynamical flow region. Ray tracing through synthetic surface currents will provide a benchmark. It will be shown that a broad range of the current scales can be replaced by the stochastic parametrization (3.11) without affecting ray scattering and action distribution. Theoretical results (4.3)–(4.12) will suggest approximate analytic solutions.

5.1. Surface current dynamics

Simplified upper ocean dynamics are considered to follow

(5.1)\begin{equation} (\partial_t + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) \varTheta= 0, \quad \text{with } \boldsymbol{v}={-} \boldsymbol{\nabla}^{\bot} (-\varDelta)^{-\xi} \varTheta, \end{equation}

where $\varTheta$ stands for the buoyancy, $\boldsymbol {\nabla }^{\bot }$ for the curl, and $\varDelta$ for the Laplacian. Two extreme cases are the surface quasi-geostrophic (SQG) dynamics ($\xi = \frac 12$) (Held et al. 1995; Lapeyre 2017), and the 2-D Euler dynamics ($\xi = 1$). The SQG dynamics has an extreme locality (kinetic energy spectrum slope $-5/3$), whereas 2-D Euler dynamics has an extreme non-locality (kinetic energy spectrum slope $-3$). The objective is to test how the proposed closures apply to both dynamics to be equally useful for any more realistic upper ocean dynamics. Additionally, test cases are developed to assess the multiscale stochastic closure in both homogeneous and heterogeneous propagating media. Moreover, we would like to challenge our closure beyond the validity of rapid wave models. In our first test case, surface fast waves travel in a homogeneous and isotropic SQG turbulence. Then we simulate waves propagating in a spatially heterogeneous 2-D Euler turbulence, mimicking an oceanic jet. For both SQG and 2-D Euler dynamics, a reference simulation is obtained at resolution $512\times 512$ for a $1000$ km$^2$ domain, with the help of a pseudo-spectral code (Resseguier et al. 2017b, 2020b). Once initialized, the current velocity $\boldsymbol {v}$ is approximately 0.1 m s$^{-1}$ for the homogeneous turbulence, and 1 m s$^{-1}$ for the jet (see figure 3).

Figure 3. Current velocity norms of (a) the SQG homogeneous turbulence and (b) the 2-D Euler jet current at high resolution ($512\times 512$).

5.2. Rays scattering in homogeneous SQG turbulence

A wave system enters the bottom boundary, propagating to the top. The carrier incident wave has intrinsic wave group velocity 10 m s$^{-1}$, i.e. wavelength $\lambda = 250$ m. Its envelope is Gaussian with isotropic spatial extension $30 \lambda$. Figures 4(a) and 5(a) illustrate the resulting dynamics, spreading the wavevectors (figure 5) of the incoming waves. From bottom to top, spectral diffusion occurs (figure 5) in the direction orthogonal (here $k_x$) to the propagation (here $k_y$), in line with the additive noise appearing in (4.6). This scattering accelerates – along the propagation – the wave position spread (figure 4). This acceleration is explained by the ray equation (4.3) dominated by the intrinsic wave group velocity.

Figure 4. Swell (wavelength $\lambda = 250$ m) interacting with (a) a high-resolution ($512\times 512$) deterministic SQG current, (b) a low-resolution ($32\times 32$) deterministic SQG current, and (c) a low-resolution ($32\times 32$) deterministic SQG current plus (one realization of) the time-uncorrelated stochastic model – coloured by the corresponding wave amplitude, $h(t)=\sqrt {\omega _0(k(t))\,N(t=0)}$ (right-hand side colour bars) – computed by forward advection and superimposed on the current vorticity $\omega = \boldsymbol {\nabla }^{\bot } \boldsymbol {\cdot } \boldsymbol {v}$. The red crosses indicate where the bidirectional wave spectra of figure 5 are computed.

Figure 5. Bidirectional wave spectra, computed by backward advection, at eight locations along a vertical axis (the mean wave propagation direction) resulting from a swell interacting with (a) a high-resolution ($512\times 512$) deterministic SQG current, (b) a low-resolution ($32\times 32$) deterministic SQG current, and (c) a low-resolution ($32\times 32$) deterministic SQG current plus (one realization of) the stochastic model (3.11). The spatial locations where the spectra are calculated are highlighted in figure 4 by the red crosses.

To mimic a badly resolved $\bar {\boldsymbol v}$ field, ${\boldsymbol v}$ is smoothed at resolution $32\times 32$. Using this coarse-scale current in figures 4(b) and 5(b), the scattering – described in the previous paragraph – is strongly depleting in comparison to ray tracing in fully-resolved turbulence. The spectral diffusion induced by small-scale turbulence is missing. Thus the spatial spreading also is narrower compared to high-resolution simulations. A stochastic current $\boldsymbol {v}'$ is then added for ray tracing (4.1). This stochastic component is divergence-free and has a self-similar distribution of energy across spatial scales (3.11) (see figure 1). The resulting spatial and spectral spreads are now comparable to simulations with high-resolution currents. For this setting, the stochastic closure provides satisfying results for a sufficiently well-resolved large-scale current. The key decorrelation ratio $\epsilon = ({l_{v'}}/{l_{v}}) ({\|\boldsymbol {v}\|}/{v_g^0})$ indeed depends on the resolution through $l_{v'}$. The large-scale current $\bar {\boldsymbol {v}}$ is resolved on a $32\times 32$ grid, i.e. with resolution $l_{v'} = ({\|\boldsymbol {\nabla } {\boldsymbol {v}}^{{\rm T}}\|}/{\|\boldsymbol {\nabla } {\boldsymbol {v}'}^{{\rm T}}\|}) l_{v} = 0.33 l_{v}$. As such, $\epsilon = 4.1 \times 10^{-3}$, computed with $v_g^0 \approx 10$ m s$^{-1}$ and $\|\boldsymbol {v}\|\approx 0.12$ m s$^{-1}$, so ${\|\boldsymbol {v}\|}/{v_g^0} \approx 1.2\times 10^{-2}$, which is sufficiently small to make the proposed model applicable.

From the ADSD estimate (3.10) (illustrated by figure 1) and (4.7)–(4.8), evaluations of the diffusivity coefficients $a_0$ and $\gamma _0$ are straightforward. As discussed previously (Smit & Janssen 2019), the spatial diffusivity is extremely weak: $a_0 = 6.4 \times 10^{-1}$ m$^2$ s$^{-1}$ (spatial variations in ray equations (4.3)–(4.4) of approximately $\sqrt {a_0 t} = 230$ m during $1$ day). In contrast, the spectral angle diffusivity is large: $3 \gamma _0 = 3.0 \times 10^{-8}$ rad$^2$ s$^{-1}$. Along our $1$-day simulation, neglecting large-scale velocity influence, (4.6) leads to Brownian wavevector angle variations $\delta \theta _k =\theta _k - \theta _k(0) = \sqrt {3 \gamma _0}\,B_t^{(4)}$ with standard deviation $\sigma _{\delta \theta _k} = \sqrt {3 \gamma _0 t} = 5.2 \times 10^{-2}$ rad $\approx 3.0^{\circ }$, eventually increasing the wave group spectral maximal extension from $\pm 2\sigma _{ k_x} = \pm 2({2{\rm \pi} }/{30 \lambda }) = \pm 1.7 \times 10^{-3}$ rad m$^{-1}$ to $\pm 2\sigma _{ k_x} \approx \pm 2 \sqrt {({2{\rm \pi} }/{30 \lambda })^2 + (k \sigma _{\delta \theta _k})^2} = \pm 3.1 \times 10^{-3}$ rad m$^{-1}$, confirmed by figure 5. This figure also illustrates the wave action diffusion induced by diffusivity $\gamma _0$, well predicted by the Eulerian wave action model (4.12). In this scattering regime, the increased angle variability leads, by advection, to a spatial spread. The simplified ray equation (4.3) gives $\delta x \approx \int _0^t v_g^0 \cos {\theta }_k \,{\rm d} t' \approx v_g^0 \int _0^t \delta {\theta }_k \,{\rm d} t' \approx v_g^0 \sqrt {3 \gamma _0} \int _0^t B_{t'}^{(4)} \,{\rm d} t'$ with maximal extension $\pm 2\sigma _{x} \approx \pm 2 v_g^0 \sqrt { \gamma _0 t^3} \approx \pm 52$ km, in agreement with figure 4. Finally, we estimate a first-order delay along the propagation

(5.1) \begin{align} \delta t = t - (y - y(0))/v_g^0 \approx \int _0^t ( 1 - \sin {\theta }_k ) {\rm d}t' \approx \frac {1}{2} \int _0^t \delta {\theta }_k^2 \,{\rm d} t' \approx \frac {3 }{2} \gamma _0 \int _0^t ( B_{t'}^{(4)} )^2 \,{\rm d} t', \end{align}

with mean value $\mathbb {E} \delta t = \frac {3}{4} \gamma _0 t^2$.

5.3. Wave groups trapped in a 2-D Euler turbulent jet

Tests are now performed for rays travelling in fast and strongly heterogeneous 2-D Euler flows. Classical fast wave models – assuming flows of weak amplitude and often uniform statistics – are expected to fail here. Jets exhibit strong current gradients (e.g. Kudryavtsev et al. 2017), creating strong ray focusing and possibly rogue events. Passing through localized spatial structures, caustics can appear, but solely from unrealistically collimated wave trains (White & Fornberg 1998; Heller, Kaplan & Dahlen 2008; Wang et al. 2023).

Occurrences strongly reduce for finite directional spread (Slunyaev & Shrira 2023). Here, wave groups are trapped in a jet, but nonlinear wave interactions are neglected. The high-resolution numerical simulations (see figure 6) reveal that even linear wave trains are well trapped in adversarial currents. Freund & Fleischman (2002) observed a similar behaviour for acoustic waves in a three-dimensional turbulent jet. Note that during our simulation, rays cross the domain several times (because of the doubly periodic boundary conditions; see Appendix E for technical details). At the top (resp. bottom) of the jet, the vorticity and thus – at first order – rays curvatures (Dysthe 2001) are negative (resp. positive). Therefore, rays oscillate around the jet. A toy model can explain this behaviour. Following the multiscale stochastic approach (4.3)–(4.6), wave scattering is also taken into account.

Figure 6. Rays facing a high-resolution ($512\times 512$) deterministic 2-D Euler jet current – coloured by the corresponding wave amplitude $h(t)=\sqrt {\omega _0(k(t))\,N(t=0)}$ (right-hand side colour bars) – computed by forward advection and superimposed on the current vorticity $\omega = \boldsymbol {\nabla }^{\bot } \boldsymbol {\cdot } \boldsymbol {v}$ (top colour bars).

For very-coarse-grained ($4\times 4$) current $\bar {\boldsymbol {v}}$, oscillation remains, but most of the scattering vanishes, as illustrated by figure 7. Moreover, the curvature of the jet creates artificial wave focusing at $t=8$ and $10$ days. Introducing a time-uncorrelated model (3.11) corrects the resolution issue in figure 8. Figure 9 plots the current ADSD. The current is strong ($\|\boldsymbol {v}\|\approx 1.4$ m s$^{-1}$), and the usual fast wave approximation cannot be applied (${\|\boldsymbol {v}\|}/{v_g^0 } \approx 1.2\times 10^{-1}$). However, the proposed modified fast wave model is valid, even at the very coarse $4\times 4$ resolution.

Figure 7. Rays facing a low-resolution ($4\times 4$) deterministic 2-D Euler jet current – coloured by the corresponding wave amplitude $h(t)=\sqrt {\omega _0(k(t))\,N(t=0)}$ (right-hand side colour bars) – computed by forward advection and superimposed on the current vorticity $\omega = \boldsymbol {\nabla }^{\bot } \boldsymbol {\cdot } \boldsymbol {v}$ (top colour bars).

Figure 8. Rays facing a low-resolution ($4\times 4$) deterministic 2-D Euler jet current plus (one realization of) the time-uncorrelated stochastic model – coloured by the corresponding wave amplitude $h(t)=\sqrt {\omega _0(k(t))\,N(t=0)}$ (right-hand side colour bars) – computed by forward advection and superimposed on the low-resolution current vorticity $\bar {\omega } = \boldsymbol {\nabla }^{\bot } \boldsymbol {\cdot } \bar {\boldsymbol {v}}$ (top colour bars).

Figure 9. The ADSD (m$^2$ s$^{-1}$/(rad m$^{-1}$)) of the resolved high-resolution jet velocity in red, low-resolution jet velocity in blue, and modelled stochastic velocity in green. The theoretical spectrum slope $- 3$ (black solid line) is imposed, consistent with homogeneous 2-D Euler dynamics. The residual ADSD (green line) is set to extrapolate that power law at small scales.

Indeed, 2-D Euler spectra are steeper than for SQG dynamics, and the length scale ratio is already significant at this resolution, ${l_{v'}}/{l_{v}} = 0.14$, and the derived time-decorrelation ratio is small: $\epsilon = ({l_{v'}}/{l_{v}}) ({\|\boldsymbol {v}\|}/{v_g^0}) = 1.6 \times 10^{-2}$.

Furthermore, by approximating the under-resolved current $\bar {v}$, an analytic stochastic solution can be obtained for a ray travelling against the current. The large-scale pattern of the jet takes a quadratic form

(5.3a,b)\begin{equation} \bar{u} \approx \bar{U}_0 - \frac 12\,\bar{\beta} \left(y-\frac{L_y}{2}\right)^2 \quad \textrm{and} \quad \bar{v} \approx 0,\quad \textrm{with}\ \bar{U}_0,\bar{\beta} <0. \end{equation}

Note that the toy model (5.3a,b) simply considers a straight jet, neglecting its curvature. For weak subgrid currents and a ray $(x_r,y_r'+{L_y}/{2},k,\theta _k)$, propagating mainly to the right, $\theta _k$ is small and the simplified ray equation (4.4) determines the group position with respect to the jet $y_r'$:

(5.4)\begin{equation} \frac{{\mathrm{d}} }{{\mathrm{d}} t} y_r'\approx v_g^0 \sin ( \theta_k) = v_g^0 \theta_k+ O(\theta_k^2). \end{equation}

For frozen turbulence, the wavenumber and hence $v_g^0$ will not vary significantly. The other ray equation (4.3) localizes the group along the jet, $x_r \approx x_r(0) + (v_g^0-\bar {u})t$, dropping the $O(\theta _k^2)$ from now on. Moreover, $\tilde {\boldsymbol {k}}^\bot \boldsymbol {\cdot } \boldsymbol {\nabla } \bar {\boldsymbol {v}}^{{\rm T}} \tilde {\boldsymbol {k}} \approx -\partial _y \bar {u}$, and the dynamics of wavevector angle (4.6) simplifies to a stochastic oscillator equation:

(5.5)\begin{equation} \frac{{\mathrm{d}}^2 }{{\mathrm{d}} t^2} y_r'= v_g^0\, \frac{{\mathrm{d}} }{{\mathrm{d}} t} {\theta}_{k}=s-\partial_y (v_g^0 \bar{u}) + v_g^0 \sqrt{3 \gamma_0}\,\dot{B}_t^{(4)}={-} \bar{\omega}_r^2 y_r'+ v_g^0 \sqrt{3 \gamma_0}\,\dot{B}_t^{(4)}, \end{equation}

with $\bar {\omega }_r = \sqrt {|v_g^0 \bar {\beta }|}$. Here, $v_g^0 \bar {u}$ plays the role of a potential, trapping the rays in the jet vicinity, whereas the noise accounts for wave scattering. The solution of this linear equation is known (e.g. Resseguier et al. (2017a), (51)–(55)):

(5.6)\begin{align} y_r (t) = \underbrace{\frac{L_y}{2} + y_r'(0) \cos(\bar{\omega}_rt) + \frac{v_g^0}{\bar{\omega}_r}\,\theta_k(0)\sin(\bar{\omega}_rt)}_{=\,\mathbb{E} (y_r(t) )} +\underbrace{Y_{\gamma_0} \sqrt{ \bar{\omega}_r} \int_{0}^t \sin(\bar{\omega}_r(t-r))\, {\mathrm{d}} B_r^{(4)}}_{=\,y_r''(t)}, \end{align}

with $Y_{\gamma _0} = v_g^0\sqrt {{ 3 \gamma _0}/{\bar {\omega }_r^3} }$. The wavevector angle solution is similar. The solution ensemble mean $\mathbb {E} y_r$ is a simple coherent deterministic oscillator. This mean solution describes well the interaction between the group and the under-resolved current from figure 7. From the coarse-scale vorticity shear plotted in figure 10 in the vicinity of the jet, we can estimate $\bar {\beta } = - 2.7 \times 10^{-11}$ m$^{-1}$ s$^{-1}$. It yields an oscillation frequency $\bar {\omega }_r = 1.3 \times 10^{-5}$ rad s$^{-1}$, i.e. period $2{\rm \pi} /\bar {\omega }_r = 5.7$ days, in agreement with the ray tracing simulations. Note that the high-resolution vorticity shear in figure 10(a) does not suggest any relevant values to explain the ray oscillations. Only the proposed multiscale current decomposition provides a quantitative explanation for these oscillations, and by extension for trapping rays inside the jet. Added to the mean solution, the random parts $y_r''(t)$ are continuous summations of zero-mean incoherent wave fluctuations. At each time $r$, the additive random forcing introduces an oscillation. But the influence of the past excitations is weighed by a sine wave due to the phase change. The group position and wavevector angle are Gaussian random variables (as linear combinations of independent Gaussian variables). Therefore, their finite-dimensional law (i.e. the multi-time probability density function) is entirely defined by their mean and covariance functions. Specifically,

(5.7)\begin{align} \mathbb{E} (y_r''(t)\,y_r''(t+\tau))=\tfrac{1}{4} Y_{\gamma_0}^2 (\cos(\bar{\omega}_r \tau )\,(2 \bar{\omega}_r t - \sin( 2 \bar{\omega}_r t)) +\sin(\bar{\omega}_r \tau )\,(1 - \cos( 2 \bar{\omega}_r t))). \end{align}

In particular, the variance of the vertical positions reads $\sigma _y^2 (t) =\frac {1}{4} Y_{\gamma _0}^2 (2 \bar {\omega }_r t - \sin ( 2 \bar {\omega }_r t))$. At $t= 2{\rm \pi} / \bar {\omega }_r$, the group has oscillated once around the jet, and the maximal position extension reaches $\pm 2\sigma _y = \pm 2\sqrt {{\rm \pi} }\,Y_{\gamma _0} = \pm 42$ km, well confirmed by ray simulations. In contrast, usual fast wave models (e.g. Smit & Janssen 2019) do not consider the interplay between smooth and rough currents, and hence solely predict a classical scattering with a much faster vertical location spreading: $\pm 2 \sigma _y = \pm 2 \sqrt {(2 {\rm \pi})^3/3}\,Y_{\gamma _0} = \pm 217$ km. For large time, our multiscale approach predicts a scaling in $t$, much slower than the usual scattering $t^3$ scaling.

Figure 10. Vorticity shear $\partial ^2_y u$ of the deterministic 2-D Euler jet current at (a) high-resolution ($512\times 512$) and (b) low-resolution ($4\times 4$), and (c) the corresponding swell system period $2{\rm \pi} /\bar {\omega }_r$. Far from the jet (${\pm }200$ km away), the vorticity shear becomes zero or even positive, so periods larger than $10$ days are cropped.

From the group vertical location and wavevector angle, we can also solve (4.5) analytically to estimate the group wavenumber variations. For small wavevector angles, $-\int _0^t \bar {\sigma } \sin (\zeta ) \,{\rm d} t' \approx 2\int _0^t \bar {\omega } \theta _k \,{\rm d} t' = 2 \bar {\beta }\int _0^t y_r'\theta _k \,{\rm d} t'$, and (4.9) together with the analytic solutions for $y_r'$ and $\theta _k$ give a closed stochastic expression for the group wavenumber. Thus the wavenumber factor $\exp (2 \bar {\beta }\int _0^t y_r'\theta _k \,{\rm d} t')$ oscillates at frequency $2\bar {\omega }_r$, and the oscillations modulate the wave amplitude: $h=\sqrt {E}=\sqrt {\omega _0 N} = \mathrm {constant}\times k^{1/4}$. The modulations are associated with wave–current energy exchanges (Boury et al. 2023), visible in the coloured rays of figures 6, 7 and 8 when the groups enter and exit the jet.

Finally, the conditional ray distribution $p(\boldsymbol {x},\boldsymbol {k} \,|\, \boldsymbol {x}_r^0,\boldsymbol {k}^0,t)$, the action spectrum mean from (4.10) and the action mean from (4.11) can all be derived. For a system initially localized in $(0,{L_y}/{2})$ with action $A^0$, wavenumber $k^0$ and a $\sigma _{\delta \theta _k}^0$-width Gaussian angular spreading, propagating to the right, the action mean reads

(5.8)\begin{equation} \mathbb{E} A(x,y,t)= A^0\delta(x - (v_g^0(k^0)-\bar{u}(y))t)\, \mathcal{N}\left(y-\left.\frac{L_y}{2}\right|\tilde{\sigma}_y^2(t)\right), \end{equation}

with $\mathcal {N} (\cdot \,|\, \tilde {\sigma }_y^2(t))$ a Gaussian function with variance $\tilde {\sigma }_y^2(t) = \sigma _y^2(t) + (({v_g^0}/{\bar {\omega }_r}) \sin (\bar {\omega }_r t) \sigma _{\delta \theta _k}^0)^2$. The action is advected in the horizontal direction, and slowly diffuses along the vertical direction.

6. Conclusion

Developed to generalize the ray path concept for waves propagating over a heterogeneous turbulence, a practical stochastic framework is derived. For fast waves, the smallest scales of a turbulent flow decorrelate along the wave propagation. Flows with steeper spectra decorrelate faster, leading to a broader validity range of fast wave approximations. The proposed framework encodes both large-scale refraction and random scattering effects on wave statistical properties. The mean wave action statistics are directly linked to resolved strain rate and vorticity, but also to unresolved kinetic energy spectral properties. Both Eulerian and Lagrangian views are presented. A convenient calibration method is also proposed for the subgrid parametrization.

As anticipated, random horizontal currents delay wave arrival and augment the initial radiative transport equation with a directional diffusive term. These phenomena are illustrated with numerical simulations, analytical solutions and quantitative proxies describing weak homogeneous turbulence. Using these proxies, measured delays in ray arrivals, estimated wave energy spectral characteristics and decays, and/or varying directional spread, will then be more quantitatively interpreted. It will lead to valuable information about underlying flow properties.

The generalized fast wave approximation does takes into account wavenumber variation and handles strong heterogeneous flows, like localized jets with strong current gradients. As compared to numerical simulations, numerical and theoretical results explain and quantify ray trapping effects by jets, unlike usual fast wave approaches.

Among the fast wave literature, isotropic diffusion and hence wavenumber diffusion may (e.g. Voronovich 1991) or may not (e.g. Bôas & Young 2020) come into play (see Appendix C for details). Future works could adapt our convenient stochastic calculus framework to the second models family. Besides, further analytical developments could consider finite-size wave groups, their dynamics (Jonsson 1990; White & Fornberg 1998) and statistical distributions, or alternatively the Eulerian action dynamics (4.12) with all its multi-point stochastic structure. When achieved, this next theoretical development could provide new means to analyse wave dynamics with subsequent fast simulations of ensembles. Beside comprehension and analysis, our stochastic simulation tools aim to eventually facilitate future ensemble-based data assimilation algorithms (Smit et al. 2021).